Electrodynamics and Radiation from Rotating Neutron Star Magnetospheres

Abstract

Neutron stars are compact objects rotating at high speed, up to a substantial fraction of the speed of light (up to 20% for millisecond pulsars) and possessing ultra-strong electromagnetic fields (close to and sometimes above the quantum critical field of 4.4 $\times {10}^9$ $\mathrm{T}$). Moreover, due to copious $e^{\pm }$ pair creation within the magnetosphere, the relativistic plasma surrounding the star is forced into corotation up to the light cylinder where the corotation speed reaches the speed of light. The neutron star electromagnetic activity is powered by its rotation which becomes relativistic in the neighborhood of this light cylinder. These objects naturally induce relativistic rotation on macroscopic scales about several thousands of kilometers, a crucial ingredient to trigger the central engine as observed on Earth. In this paper, we elucidate some of the salient features of this corotating plasma subject to efficient particle acceleration and radiation, emphasizing several problems and limitations concerning current theories of neutron star magnetospheres. Relativistic rotation in these systems is indirectly probed by the radiation produced within the magnetosphere. Depending on the underlying assumptions about particle motion and radiation mechanisms, different signatures on their light curves, spectra, pulse profiles and polarization angles are expected in their broadband electromagnetic emission. We show that these measurements put stringent constraints on the way to describe particle electrodynamics in a rotating neutron star magnetosphere.

1. Introduction

Neutron stars are compact objects produced by the explosion of a massive star or by the collapse of an accreting white dwarf reaching the Chandrasekhar limit of about $1.44M_{\odot }$ [1,2] where $M_{\odot }$ is the solar mass. They represent the ultimate fate of the stellar evolution of massive stars before the black hole stage. During their birth, their angular speed and magnetic field are amplified by several orders of magnitude. It is not yet clear what mechanisms are able to produce the expected fields in the range of $B\approx {10}^8-{10}^{11}\mathrm{T}$, but a dynamo effect and magnetic flux freezing during the collapse are certainly key processes. Neutron star rotation periods span four decades from several milliseconds to tenths of seconds. The stellar magnetic field drags charged particles into corotation with the star. Relativistic corotating speeds are reached at the light cylinder defined by

$$r_{\mathrm{L}}=\frac{c}{\omega }=\frac{cP}{2\pi }=48\mathrm{k}\mathrm{m}\left(\frac{P}{1\mathrm{m}\mathrm{s}}\right)$$

(1)

where c is the speed of light, $P=2\pi /\omega$ the pulsar period and $\omega$ its rotation rate. Moreover, their compactness places them closest to the black hole stage because

$$\frac{R_{\mathrm{s}}}{R}=0.345\left(\frac{M}{1.4M_{\odot }}\right){\left(\frac{R}{12\mathrm{k}\mathrm{m}}\right)}^{-1}$$

(2)

where $R_{\mathrm{s}}=2GM/c^2$ is the Schwarzschild radius, M and R the neutron star mass and radius respectively and G the gravitational constant. Neutron stars are therefore places in the universe where general relativity and quantum electrodynamics act together to sustain their electromagnetic activity.

Simple but realistic orders of magnitude for neutron star rotation periods and magnetic field strengths are easily derived from conservation of angular momentum and magnetic flux during the collapse of the progenitor. These conservation laws imply

$$\begin{array}{cc}\hfill M\omega R^2& =M_{*}{\omega }_{*}{R}_{*}^2\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill B{R}^2& =B_{*}{R}_{*}^2\hfill \end{array}$$

(3b)

where the star subscript * refers to quantities relative to the progenitor. Magnetic field and rotation frequency therefore increase by a factor as large as $R_{*}^2/R^2\approx {10}^{10}$ if dissipation processes are neglected and all angular momentum and magnetic flux of the progenitor go to the neutron star. The magnetic flux conservation argument, Equation (3b), was first discussed by [3]. The above estimates are however probably largely overestimated because not all the progenitor is collapsing and because its mass is not conserved during the implosion [4]. Only the iron core of a massive star produces a neutron star star whereas the outer shells are expelled. Electrodynamics in a relativistically rotating frame is a crucial ingredient in our understanding and modelling of neutron star magnetospheres. In this paper, we discuss some issues about rotating magnetospheres and their radiative properties.

First we briefly remind the rotating coordinate system used and its implication for field transformations and especially for Maxwell equations in Section 2. Next we discuss the electromagnetic field expected inside and outside a neutron star in Section 3. The motion of charged particles in this field is exposed in Section 4. Such trajectories can be computed in solutions found from numerical simulations as explained in Section 5. However, corotation is not compulsory and differentially rotating magetospheres have been found as shown in Section 6. Some clues about the electrodynamics of pulsar magnetospheres can be gained from their radiation as explained in Section 7. Possible extensions to general relativity and multipolar magnetic fields are discussed in Section 8. We conclude with a brief summary in Section 9.

2. Rotating vs. Inertial Frame

Neutron stars are mainly observed through their pulsed emission detected in the radio wavelength [5] but also at very high-energy by gamma-ray photons in the MeV/GeV band [6] or through thermal X-ray emission from hot spots on the surface [7]. Such emission is attributed to ultra-relativistic charged particles flowing inside the magnetosphere. The very stable and periodic pulsation is explained by the stellar rotation. It is believed that these particles corotate with the star almost up to the light cylinder and generate a relativistic magnetized outflow outside the light cylinder known as the pulsar wind [8]. Description of the flow taking properly into account this corotation is therefore central to our understanding of neutron star magnetospheric emission. Before going into the dynamics of this plasma and its underlying particle trajectories let us briefly review relativistic rotating frames and related electrodynamics problems. An excellent review about relativistically rotating frames can be found in [9].

2.1. Metric

The space-time geometry of a rotating system is best described in a cylindrical coordinate system labeled by $(t,r,\phi ,z)$. The coordinate transformation from an inertial observer $(t,r,\phi ,z)$ to an observer rotating at the stellar angular speed $(t^{\prime },r^{\prime },{\phi }^{\prime },z^{\prime })$ is usually written as

$$t^{\prime }=t;r^{\prime }=r;{\phi }^{\prime }=\phi -\omega t;z^{\prime }=z.$$

(4)

It is essential to realize that this transformation does not lead to a new orthonormal basis in the rotating frame. Therefore, physical quantities measured locally by an observer in this coordinate system cannot be directly read off from these basis vectors. Indeed, the space-time geometry of a uniformly rotating frame is given in the rotating observer frame by the metric

$$d{s}^2=c^{2}d{t}^2-d{r}^2-\frac{r^{2}d{\phi }^2}{1-r^2{\omega }^2/c^2}-d{z}^2.$$

(5)

The coordinates $(t^{\prime },r^{\prime },{\phi }^{\prime },z^{\prime })$ are not related to any observer because they do not form an orthonormal basis. To relate measurements between the inertial and the rotating observers, it is more convenient to introduce two orthonormal bases attached to each observer. The transformation from the inertial to the rotating frame is simply a Lorentz boost as is always the case when performing reference frame transformations between orthonormal bases. Explicitly, the orthonormal basis vectors of the rotating frame are

$$\begin{array}{cc}\hfill {\mathbf{e}}_0^{\prime }& =\Gamma \left({\mathbf{e}}_t+\beta (-sin\Phi {\mathbf{e}}_x+cos\Phi {\mathbf{e}}_y)\right)=\Gamma (1,-\beta sin\Phi ,\beta cos\Phi ,0)\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill {\mathbf{e}}_1^{\prime }& =cos\Phi {\mathbf{e}}_x+sin\Phi {\mathbf{e}}_y=(0,cos\Phi ,sin\Phi ,0)\hfill \end{array}$$

(6b)

$$\begin{array}{cc}\hfill {\mathbf{e}}_2^{\prime }& =\Gamma \left(\beta {\mathbf{e}}_t-sin\Phi {\mathbf{e}}_x+cos\Phi {\mathbf{e}}_y\right)=\Gamma (\beta ,-sin\Phi ,cos\Phi ,0)\hfill \end{array}$$

(6c)

$$\begin{array}{cc}\hfill {\mathbf{e}}_3^{\prime }& ={\mathbf{e}}_z=(0,0,0,1)\hfill \end{array}$$

(6d)

where the phase is $\Phi =\omega t$, the relative speed is $\beta =\frac{\omega r}{c}$, and the Lorentz factor $\Gamma ={(1-{\beta }^2)}^{-1/2}$. This coordinate transformation is of Lorentz type, going from a Minkowskian metric to another Minkowskian metric. Building on this orthonormal basis, it is easy to deduce the relation between the electromagnetic fields measured by the two observers as we now show.

2.2. Electrodynamics

The electromagnetic field tensor in an inertial frame is derived from the electromagnetic quadri-potential $A_i=(\varphi /c,-\mathbf{A})$ ($\varphi$ being the scalar potential and $\mathbf{A}$ the vector potential adopting the metric with signature (+,-,-,-)) by $F_{ik}={\partial }_i{A}_k-{\partial }_k{A}_i$. The electromagnetic field tensors in the inertial and rotating frame are related by a special relativistic transformation according to the previous section from the basis transformation Equation (6). They are also found from the definition using the observer 4-velocity $\mathbf{u}$ by projection of the electromagnetic field tensor $F_{ik}$ and its dual ${{}^{*}F}_{ik}=\frac{1}{2}{ϵ}_{ikmn}{F}^{mn}$ [10] onto the observer world line such that

$$\begin{array}{cc}\hfill E_i/c& =F_{ik}{u}^k\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill B_i& ={{}^{*}F}_{ik}{u}^k.\hfill \end{array}$$

(7b)

Explicit computations of these transformations show that it is simply the Lorentz transformation of the electromagnetic field between two inertial observers with relative velocity $\mathbf{V}=r\omega {\mathbf{e}}_{\phi }$. Introducing the normalized velocity by $\beta =\mathbf{V}/c$, the transformations of the electric and magnetic field vectors are

$$\begin{array}{cc}\hfill {\mathbf{E}}^{\prime }& =\Gamma \left[\mathbf{E}-\frac{\Gamma }{\Gamma +1}\left(\beta ·\mathbf{E}\right)\beta +c\beta \wedge \mathbf{B}\right]\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill {\mathbf{B}}^{\prime }& =\Gamma \left[\mathbf{B}-\frac{\Gamma }{\Gamma +1}\left(\beta ·\mathbf{B}\right)\beta -\beta \wedge \mathbf{E}/c\right]\hfill \end{array}$$

(8b)

where primed quantities are defined in the rotating frame. There is nothing special about rotation if transformations are made locally between inertial observers. This holds for aberration and Doppler effects when emission emanates from within the light cylinder ($r<r_{\mathrm{L}}$). These results are easily extended to general relativity when gravitation is included. However special relativity applies if transformations are made locally between inertial observers [11].

2.3. Doppler Effect

The Doppler effect is subject to the same transformation as the electromagnetic field. A simple Lorentz transformation between inertial frames holds to relate photon propagation direction $\mathbf{n}=(n_r,n_{\theta },n_{\phi })$ and frequency $\nu$ in the observer frame and in the instantaneous inertial frame coinciding locally with the corotating frame. It can be checked by a direct derivation from the transformation of coordinates between both frames as given by Equation (6). For the sake of completeness, they are given in a spherical coordinate systems $(r,\vartheta ,\phi )$ by

$$\begin{array}{cc}\hfill \nu & ={\gamma }_{\mathrm{obs}}{\nu }^{\prime }(1+{\beta }_{\mathrm{obs}}{n}_{\phi }^{\prime })\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill n_r& =\frac{n_r^{\prime }}{{\gamma }_{\mathrm{obs}}(1+{\beta }_{\mathrm{obs}}{n}_{\phi }^{\prime })}\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill n_{\vartheta }& =\frac{n_{\vartheta }^{\prime }}{{\gamma }_{\mathrm{obs}}(1+{\beta }_{\mathrm{obs}}{n}_{\phi }^{\prime })}\hfill \end{array}$$

(9c)

$$\begin{array}{cc}\hfill n_{\phi }& =\frac{{\beta }_{\mathrm{obs}}+n_{\phi }^{\prime }}{1+{\beta }_{\mathrm{obs}}{n}_{\phi }^{\prime }}.\hfill \end{array}$$

(9d)

This aberration formula also holds in general relativity if the Lorentz factor ${\gamma }_{\mathrm{obs}}$ and velocity ${\beta }_{\mathrm{obs}}$ are properly defined as the true physical quantities measured by a corotating observer. Mathematically, this requires switching from a general curvilinear coordinate system, like the metric of a rotating disk, to an orthonormal coordinate system associated with the Minkowskian metric. Only the latter coordinates have a clear physical interpretation, the former being only appropriate (or not) coordinates to describe the problem.

To summarize, in any local Lorentzian frame, using an orthonormal basis, the Doppler factor is written geometrically as

$$\mathcal{D}=\frac{1}{\Gamma (1-\beta ·\mathbf{n})}.$$

(10)

It enables relating the photon propagation directions in both frames as

$$\mathbf{n}=\frac{1}{\mathcal{D}}\left[{\mathbf{n}}^{\prime }+\Gamma \left(\frac{\Gamma }{\Gamma +1}\left(\beta ·{\mathbf{n}}^{\prime }\right)+1\right)\beta \right]$$

(11)

with the usual redshift phenomenon relating the frequency in the rotating frame ${\nu }^{\prime }$ to the frequency in the observer frame $\nu$ by

$$\nu =\mathcal{D}{\nu }^{\prime }.$$

(12)

The aberration effect can only be computed for a physical realization of a corotating frame. It would fail right at the light cylinder or outside it. Unfortunately, a smooth transition from the magnetosphere $r\le r_{\mathrm{L}}$ to the wind $r\ge r_{\mathrm{L}}$ is required for modelling the pulsar broad band emission. This tells us that the introduction of a corotating frame where the electromagnetic field and plasma motion are both stationary will not help in advancing our understanding of neutron star electrodynamics. It is preferable to keep the physical quantities expressed in an inertial frame even if the coordinate system can be advantageously described in a corotating frame (not related to any comoving observer). We expose such a technique in the next section for the evolution of the electromagnetic field.

2.4. Maxwell Equations in a Rotating Coordinate System

Measuring the electromagnetic field in the corotating system is possible up to the light cylinder. However, outside this surface, the metric has no physical significance anymore. It is impossible to describe the neutron star electrodynamics in whole space with the field locally measured by a corotating observer because such an observer does not exist when $r\ge r_{\mathrm{L}}$.

Nevertheless, it is relevant and useful to keep the definition of the electromagnetic field as measured in the inertial reference frame but using the rotating coordinate system to localize it with $(t^{\prime },r^{\prime },{\phi }^{\prime },z^{\prime })$. In such a case the time derivative of any vector field $\mathbf{A}$ is given by [12]

$$\frac{\partial \mathbf{A}}{\partial t}=\frac{\partial \mathbf{A}}{\partial t^{\prime }}+\mathrm{curl}({\mathbf{V}}_{\mathrm{rot}}\wedge \mathbf{A})-{\mathbf{V}}_{\mathrm{rot}}\mathrm{div}\mathbf{A}.$$

(13)

The solid body corotation velocity, expressed in the inertial frame, is simply

$${\mathbf{V}}_{\mathrm{rot}}=\omega \wedge \mathbf{r}=r\omega {\mathbf{e}}_{\phi }.$$

(14)

Please note that it is not restricted to remain less than the speed of light. There is no singularity in Equation (13) when crossing the light cylinder. With the correspondence established in Equation (13), in the rotating coordinate system, Maxwell equations become

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{B}}{\partial t^{\prime }}& =-\mathrm{curl}(\mathbf{E}+{\mathbf{V}}_{\mathrm{rot}}\wedge \mathbf{B})\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{E}}{\partial t^{\prime }}& =\mathrm{curl}(c^2\mathbf{B}-{\mathbf{V}}_{\mathrm{rot}}\wedge \mathbf{E})-\frac{\mathbf{j}}{{\epsilon }_0}+{\mathbf{V}}_{\mathrm{rot}}\mathrm{div}\mathbf{E}.\hfill \end{array}$$

(15b)

Note however the subtleties that $\mathbf{E}$ and $\mathbf{B}$ are still defined as observed in the inertial frame. No particular problem arises at the light cylinder when Maxwell equations are written in this way. In Section 5 about numerical simulations, we use this mathematical formulation to solve for the force-free (FFE) and radiative magnetosphere for an oblique rotator as a typical example. However, first we have to define the electromagnetic field inside and outside the star in vacuum.

3. Electromagnetic Field Inside and Outside the Star

To first approximation, a neutron star can be assimilated to a very good conductor. We therefore assume that the electric field inside the star as seen by an observer at rest with respect to the star, the so-called comoving observer, vanishes. In the inertial frame of a distant observer, the electric field $\mathbf{E}$ is given by the usual Lorentz transformation

$$\mathbf{E}+{\mathbf{V}}_{\mathrm{rot}}\wedge \mathbf{B}=\mathbf{0}$$

(16)

where $\mathbf{B}$ is the magnetic field in the same observer frame. At this stage, already several implicit assumptions must be done. Should we assume a prescribed magnetic field in the observer frame $\mathbf{B}$ or in the corotating frame ${\mathbf{B}}^{\prime }$? Both situations are obviously not identical. If the magnetic field is fixed in the observer frame, then in the star corotating frame we find

$${\mathbf{B}}^{\prime }=\frac{\mathbf{B}}{{\Gamma }_{\mathrm{rot}}}+\frac{\left({\beta }_{\mathrm{rot}}·\mathbf{B}\right)}{{\Gamma }_{\mathrm{rot}}+1}{\beta }_{\mathrm{rot}}$$

(17)

with ${\Gamma }_{\mathrm{rot}}={(1-{\beta }_{\mathrm{rot}}^2)}^{-1/2}$. If the magnetic field is fixed in the rest frame of the star, which seems more reasonable, then the observer will measure a magnetic field

$$\mathbf{B}={\Gamma }_{\mathrm{rot}}\left[{\mathbf{B}}^{\prime }-\frac{{\Gamma }_{\mathrm{rot}}}{{\Gamma }_{\mathrm{rot}}+1}\left({\beta }_{\mathrm{rot}}·{\mathbf{B}}^{\prime }\right){\beta }_{\mathrm{rot}}\right].$$

(18)

Please note that both assumptions leads to ${\beta }_{\mathrm{rot}}\wedge \mathbf{B}={\Gamma }_{\mathrm{rot}}{\beta }_{\mathrm{rot}}\wedge {\mathbf{B}}^{\prime }$ so Equation (16) remains valid in any case inside the perfectly conducting star. Practically, the corotation speed inside the star is always weakly or mildly relativistic, therefore ${\beta }_{\mathrm{rot}}\ll 1$ reducing both approaches to $\mathbf{B}\approx {\Gamma }_{\mathrm{rot}}{\mathbf{B}}^{\prime }$. An exact analytical solution for the radiating electromagnetic field has been given by [13] assuming a dipolar magnetic field in the inertial frame and neglecting relativistic effects. Close to the light cylinder, relativistic effects have been taken into account as investigated by [14]. Moreover, the current induced inside the star by its rotation influences the magnetic field itself as shown by [15].

Consequently, the description of the magnetic field inside the star is already biased by some assumptions on its geometry in the rotating or inertial frame. What happens outside, in its magnetosphere? Also, at large distances, outside the light cylinder $r>r_{\mathrm{L}}$ there exist no more physical frame in corotation with the star. The Lorentz transformation of the electromagnetic field does not apply anymore. It is preferable to stay in the observer frame without any reference to a corotating frame.

The strong electric field induced by the rotating star pulls out charged particles, filling the magnetosphere with relativistic electron/positron pairs [8]. To good accuracy, this magnetosphere is assumed to be in force-free equilibrium, neglecting particle inertia and fluid temperature, as well as gravity because of the electromagnetic field strength producing forces several orders of magnitude stronger than the gravitational attraction. The magnetosphere is therefore set into corotation with the star up to the light cylinder at a radius $r=r_{\mathrm{L}}$. Outside this cylinder, plasma corotation cannot be maintained. The flow of leptons generates space charges and currents that significantly modify the electromagnetic structure. The charge density ${\rho }_{\mathrm{e}}$ produced by the electric field, following Maxwell-Gauss law and the perfect conductor hypothesis Equation (16) diverges right at the light cylinder because

$${\rho }_{\mathrm{e}}={\epsilon }_0\mathrm{div}\mathbf{E}=-2{\epsilon }_0\frac{\omega ·\mathbf{B}}{1-r^2/r_{\mathrm{L}}^2}$$

(19)

unless $\omega ·\mathbf{B}=0$ on this surface. The magnetic field must therefore adjust itself in order to keep the constrain $\omega ·\mathbf{B}=0$ at $r=r_{\mathrm{L}}$ which implies $B_{\mathrm{z}}\left(r_{\mathrm{L}}\right)=0$. Another possibility would be to break the corotation approximation by introducing some dissipation through an ad hoc resistivity or through a radiative dissipation mechanism as shown in subsequent sections. Once the electromagnetic field settled down, we need to investigate particle motion in these fields.

4. Particle Motion in the Magnetosphere

To predict the radiation emanating from neutron star magnetospheres, a good understanding of particle trajectories in the electromagnetic field produced by these stars is required. Several attempts focused exclusively on motion restricted to within the light cylinder. This limitation avoids the problem of transformations involving larger than the speed of light relative velocities between inertial frames. Radiation from outside the light cylinder therefore requires another description. This artificial transition between the magnetosphere and the wind is far from satisfactory. Claims have been made about a new technique to compute emission everywhere, for instance in the vacuum Deutsch field [16], but we will show that at least in some cases it fails too to smoothly join inside and outside light cylinder electrodynamics. Moreover, their model based on the assumption that Equation (16) is valid within the magnetosphere is inconsistent with the vacuum assumption of a Deutsch field. A more satisfactory solution includes the plasma feedback self-consistently as proposed by [17] in order to avoid superluminal particle speed everywhere. In the following subsections, we summarize the most studied maybe not the most effective models of particle trajectories.

Three different approximations for the particle motion have been tried. Indeed a charged particle can be seen as

• following magnetic field lines ${\mathbf{B}}^{\prime }$ in the corotating frame.
• following magnetic field lines $\mathbf{B}$ in the inertial frame in addition to a corotation imposed by the stellar rotation.
• following the ultra-relativistic radiation reaction limit leading to the so-called Aristotelian dynamics. Actually, this limit can be explained by Newtonian dynamics in a stationary regime balancing electric acceleration and radiation friction.

These different views are not equivalent to each other because they assume different electromagnetic field structures, either fixed in the rotating frame or in the inertial frame. Let us discuss these approaches in depth starting with the corotating frame view.

4.1. Corotating Frame

Viewed from the corotating frame, particles are assumed to follow magnetic field lines along ${\mathbf{B}}^{\prime }$. Therefore in this frame the particle velocity is given by

$${\mathbf{v}}^{\prime }=v^{\prime }\frac{{\mathbf{B}}^{\prime }}{B^{\prime }}=v^{\prime }{\mathbf{n}}_{\mathrm{B}}^{\prime }$$

(20)

where $v^{\prime }$ is the particle speed along the field line ${\mathbf{B}}^{\prime }$ in the normalized direction ${\mathbf{n}}_{\mathrm{B}}^{\prime }={\mathbf{B}}^{\prime }/B^{\prime }$. How then to choose this field ${\mathbf{B}}^{\prime }$? Some authors used in the past the relation ${\mathbf{B}}^{\prime }=\mathbf{B}$ which is only correct to second order in ${\beta }_{\mathrm{rot}}$, a results derived from the coordinate transformation between inertial frame and rotating coordinate systems. However, this equality was used by several authors to compute pulsar high-energy light curves at high altitude, up to a substantial fraction of the light cylinder [18,19,20,21]. Light curves and sky maps derived from this model are sensitive to the upper boundary of the radiating zone. However, this dependence is undesirable. Moreover, the rotating coordinate system is not an orthonormal basis, therefore ${\mathbf{B}}^{\prime }$ should not be interpreted as the local magnetic field measured by a rotating observer. It must be computed according to the Lorentz transformation [16]. Nevertheless, both descriptions agree to good accuracy for non-relativistic corotating speeds.

Going into the rotating frame synchronous with the neutron star rotation can lead to misinterpretation of the physical electromagnetic field measured by a local observer. Moreover, the corotating frame cannot be extended beyond the light cylinder radius $r_{\mathrm{L}}$. Such description therefore faces severe difficulties to deal with the entire neutron star magnetosphere and is inadequate to efficiently model them from the surface to large distances within the striped wind $r\gtrsim r_{\mathrm{L}}$. Much better we think is to perform all calculations in the inertial frame of a distant observer as we now describe in the next two sections.

4.2. Corotating Velocity

When staying in the observer frame, without reference to any rotating frame, the velocity is described by a velocity component along the magnetic field lines $\mathbf{B}$, now expressed in the inertial frame and denoted by $v_{\Vert }$, and a velocity component due to the dragging by the star denoted by ${\mathbf{V}}_{\mathrm{rot}}$. Therefore we write

$$\mathbf{v}=v_{\Vert }\frac{\mathbf{B}}{B}+{\mathbf{V}}_{\mathrm{rot}}.$$

(21)

The velocity along the field line $v_{\Vert }$ must be chosen in order to keep the total speed smaller than the speed of light, $v<c$. It requires a special configuration of the magnetic field $\mathbf{B}$ with an increasing toroidal component to compensate for the linear increase in corotation speed as given by ${\mathbf{V}}_{\mathrm{rot}}$. Therefore not all magnetic field configurations are permitted to fulfil this constrain.

By assumption, in some models [18,19,20,21], particles follow magnetic field lines in the corotating frame. Their distribution function is isotropic in the rest frame of the fluid. Following the previous prescriptions by [16], we assume that their Lorentz factor $\Gamma$ is constant in the observer frame such that the velocity, being a combination between propagation along field lines and corotation at speed ${\mathbf{V}}_{\mathrm{rot}}$, is

$$\mathbf{v}=v_{\Vert }^{\mathrm{c}}\mathbf{t}+{\mathbf{V}}_{\mathrm{rot}}$$

(22)

where $\mathbf{t}=\pm \mathbf{B}/B$ is the outward pointing tangent vector to the field line. Solving for the parallel velocity $v_{\Vert }^{\mathrm{c}}$ the only real and positive solution is

$$v_{\Vert }^{\mathrm{c}}=-\mathbf{t}·{\mathbf{V}}_{\mathrm{rot}}+\sqrt{{\left(\mathbf{t}·{\mathbf{V}}_{\mathrm{rot}}\right)}^2+v^2-V_{\mathrm{rot}}^2}.$$

(23)

$V_{\mathrm{rot}}$ exceeds the speed of light outside the light cylinder by definition. The term $v^2-V_{\mathrm{rot}}^2$ in the square root becomes negative and must be compensated by the term ${\left(\mathbf{t}·{\mathbf{V}}_{\mathrm{rot}}\right)}^2$ meaning that the magnetic field must be strongly bend toward the azimuthal direction ${\mathbf{e}}_{\phi }$. The Deutsch field does not satisfy this requirement and cannot be used to study photon emission within the wind if this view is adopted. Knowing the velocity, we get the Doppler factor for radiation as explained in Equation (10). This velocity field assumes that the electric field vanishes in the corotating frame. However, this requires a large amount of plasma to screen the electric field, in contradiction with the vacuum assumption made in [16]. Therefore, the aberration formula Equation (22) can only be an approximation in this case. Moreover, this approximation also fails at sufficiently large distances because the Deutsch field solution [13] possesses a magnetic field structure for which the poloïdal component does not decay fast enough with respect to the toroidal component. Real solutions to Equation (23) do not exists at several light-cylinder radii because the square root in Equation (23) becomes negative. Indeed, taking an orthogonal rotator, it can be shown that in the equatorial plane the term in the square root of Equation (23) tends to $v^2-4c^2<0$ for $r\to +\infty$ on the spiral given by $\phi +r/r_{\mathrm{L}}-\omega t=\pi /2$. In the most favorable case for which $v=c$, it actually becomes negative already at the light cylinder. Using the corotating frame does not help to go beyond the light cylinder for vacuum fields. Nevertheless, the description exposed in this section is applicable to force-free magnetospheres that exactly cancel the electric field in the frame comoving with the plasma at the electric drift speed. Only in such FFE models can this prescription be correctly applied in whole space, within the magnetosphere $(r\le r_{\mathrm{L}}$) and within the wind $(r\ge r_{\mathrm{L}}$).

Is it possible to find a formulation alleviating the need for special magnetic field configurations? In our opinion, there exist a simple and efficient way to compute particle trajectories in any electromagnetic field when moving at the radiation reaction limit. We detail this last approach in the next section.

4.3. Aristotelian Dynamics

Particles in the neutron star magnetosphere are ultra-relativistic. They copiously radiate photons during their motion. This has to be taken into account. The simplest approximation is given by the radiation reaction limit, where the radiation force, acting as a damping working against the motion, a kind of radiative friction, exactly compensates for the electric acceleration. It is sometime called Aristotelian electrodynamics because the velocity is completely and solely determined by the electromagnetic field felt locally by the particles although we believe that it is more appropriate to speak about radiation reaction motion because it can be derived from the Lorentz force with radiative friction and this according to Newtonian dynamics. The expression for the velocity has been derived in [22], but see also [23]. Assuming that particles move exactly at the speed of light (which is an excellent approximation in neutron star magnetospheres), depending on the sign of their charge, their velocity reads

$${\mathbf{v}}_{\pm }=\frac{\mathbf{E}\wedge \mathbf{B}\pm (E_0\mathbf{E}/c+c{B}_0\mathbf{B})}{E_0^2/c^2+B^2}$$

(24)

where the plus sign corresponds to positive charges and the minus sign to negative charges. Actually, the velocity is independent of the mass m over charge q ratio $q/m$, it only depends on the sign of its charge. Moreover, we introduced the electromagnetic field strengths $E_0$ and $B_0$ according to the two electromagnetic invariants $({\mathcal{I}}_1,{\mathcal{I}}_2)$ such that

$$\begin{array}{cc}\hfill {\mathcal{I}}_1& ={\mathbf{E}}^2-c^2{\mathbf{B}}^2=E_0^2-c^2{B}_0^2\hfill \end{array}$$

(25a)

$$\begin{array}{cc}\hfill {\mathcal{I}}_2& =c\mathbf{E}·\mathbf{B}=c{E}_0{B}_0\hfill \end{array}$$

(25b)

with the subsidiary condition $E_0⩾0$ ensuring that the radiation reaction force is always directed oppositely to the velocity direction. As explained in [24] these invariants are related to the electromagnetic field strength in a frame where $\mathbf{E}$ and $\mathbf{B}$ are parallel. The lepton motion can be decomposed into an electric drift part along the vector $\mathbf{E}\wedge \mathbf{B}$, a motion along magnetic field lines $\mathbf{B}$ and a motion along electric field lines $\mathbf{E}$. This last part of the motion is responsible for dissipation because the power of the Lorentz force is $q(\mathbf{E}+{\mathbf{v}}_{\pm }\wedge \mathbf{B})·{\mathbf{v}}_{\pm }=q{\mathbf{v}}_{\pm }·\mathbf{E}=Zec{E}_0\ge 0$ where $q=\pm Ze$ depending on the charge Z of the particle: positrons and electrons have $Z=1$ whereas ions have Z arbitrary. The velocity field (24) is regular in whole space, nothing singular happens at the light cylinder. It can be implemented to compute realistic pulsar light curves and spectra even in vacuum Deutsch solution as shown by [25]. This latest work serves as a starting point to investigate more deeply pulsar magnetospheric radiation by including for instance the plasma feedback onto the Deutsch field as will be shown in Section 5.

In the near field zone, i.e., close to the neutron star surface, where $E\ll cB$, the particle velocity simplifies into a motion solely along $\mathbf{B}$ such that

$${\mathbf{v}}_{\pm }=\pm c\frac{\left(\mathbf{E}·\mathbf{B}\right)\mathbf{B}}{E_0(E_0^2/c^2+B^2)}.$$

(26)

This expression can be reduced to

$${\mathbf{v}}_{\pm }=\pm c\mathrm{sign}\left(B_0\right)\frac{\mathbf{B}}{B}$$

(27)

by noting that in this weak electric field limit the magnitude of $\mathbf{B}$ is almost equal to the invariant $B_0$, namely $B^2\approx B_0^2$. Particles are accelerated mostly by the electric component parallel to the magnetic field. The surfaces $\mathbf{E}·\mathbf{B}=0$ are of particular interest because the velocity changes sign when the particle crosses this region. It is called a force-free surface and represents trapping regions for those particles [26,27,28].

The concept of magnetic field line in vacuum is misleading and specifying motion along a particular field line is not well defined in the general case. This requires some caution about the interpretation of the corotation speed ${\mathbf{V}}_{\mathrm{rot}}$. The way to follow the particle trajectory replaces this velocity by a special frame in which the electric field is parallel to the magnetic field leading to the Aristotelean dynamics discussed before. Indeed, the velocity ${\beta }_{\Vert }c$ required by the Lorentz transformation to get this condition is [29]

$$\frac{{\beta }_{\Vert }}{1+{\beta }_{\Vert }^2}=\frac{c\mathbf{E}\wedge \mathbf{B}}{E^2+c^2{B}^2}$$

(28)

neglecting all other curvature, gradient and polarization drifts in the limit of vanishing Larmor radius which is correct in a super strong magnetic field. In that frame, where quantities are denoted by a prime, motion is along the common direction of ${\mathbf{E}}^{\prime }$ and ${\mathbf{B}}^{\prime }$. To get the useful solution, we write the frame velocity as

$${\mathbf{V}}_{\Vert }=\frac{\mathbf{E}\wedge \mathbf{B}}{E_0^2/c^2+B^2}.$$

(29)

The electric and magnetic fields in the frame moving at speed ${\mathbf{V}}_{\Vert }$ are found by a special relativistic Lorentz boost of the electromagnetic field and gives

$$\begin{array}{cc}\hfill {\mathbf{E}}^{\prime }& =\Gamma \frac{E_0}{E_0^2/c^2+B^2}\left[\frac{E_0}{c^2}\mathbf{E}+B_0\mathbf{B}\right]\hfill \end{array}$$

(30a)

$$\begin{array}{cc}\hfill {\mathbf{B}}^{\prime }& =\Gamma \frac{B_0}{E_0^2/c^2+B^2}\left[B_0\mathbf{B}+\frac{E_0}{c^2}\mathbf{E}\right].\hfill \end{array}$$

(30b)

Electric and magnetic fields are indeed collinear because $E_0{\mathbf{B}}^{\prime }=B_0{\mathbf{E}}^{\prime }$. In this frame, particles move along the common direction of ${\mathbf{E}}^{\prime }$ and ${\mathbf{B}}^{\prime }$. Thus, the local tangent vector to the trajectory becomes ${\mathbf{t}}_{\Vert }^{\prime }=\pm {\mathbf{E}}^{\prime }/E^{\prime }=\pm {\mathbf{B}}^{\prime }/B^{\prime }$, the sign being chosen such that particles flow outwards. Therefore we replace $\beta$ by ${\mathbf{V}}_{\Vert }/c$ in Equation (11) to get a velocity field that should not be confused or seen as motion along field lines because this concept is usually ill defined for non-ideal plasmas when $\mathbf{E}·\mathbf{B}\ne 0$. Our expression for the particle velocity resembles to the Aristotelian expression given by [30]. Our velocity prescription is however more general because we do not assume that particles travel exactly at the speed of light. The speed along the common $\mathbf{E}$ and $\mathbf{B}$ direction is unconstrained and fixed by the “user” contrary to Aristotelian electrodynamics.

If particles exactly move at the speed of light, in the comoving frame this velocity becomes ${\mathbf{v}}^{\prime }=\pm {\mathbf{E}}^{\prime }/E^{\prime }=\pm {\mathbf{B}}^{\prime }/B^{\prime }$, the sign depending on the charge. Also note that the electromagnetic field strengths are $E^{\prime }=E_0$ and $B^{\prime }=B_0$. Doing the Lorentz transformation to the observer frame, noting that ${\mathbf{V}}_{\Vert }$ and ${\mathbf{v}}^{\prime }$ are orthogonal, this is nothing but Aristotelian electrodynamics. Our treatment is more general because we do not enforce the speed of light in this frame. The prescription for the velocity impacts the high-energy light curves from pulsars. This has been shown in depth by [24].

The parallel velocity ${\mathbf{V}}_{\Vert }$ in Equation (29) generalizes the electric drift approximation to field configurations with an electric field E exceeding the magnetic field $cB$. There is no need to impose the condition $E<cB$ to respect the force-free condition. However, it reduces to force-free if $E_0=0$, meaning no radiation reaction and no dissipation meanwhile requiring $\mathbf{E}·\mathbf{B}=0$. To look for plasma filled magnetospheres, we have to resort to numerical simulations in the force-free regime or in a dissipative regime because of resistivity and/or radiation damping. In the next section, we show some new results for radiative pulsar magnetospheres to be compared with the standard force-free solution for oblique rotators.

5. Numerical Simulation of Rotating Magnetospheres

Neutron stars cannot be surrounded by vacuum because particles are expelled from the surface and accelerated in the surrounding strong electromagnetic field. This is indirectly deduced from their broad band electromagnetic spectrum for which the Crab pulsar is an archetypal example [31]. Although an exact analytical solution for a rotating dipole in vacuum exists, known as Deutsch solution [13], realistic magnetospheres require the presence of plasma producing charges and currents that retroact to the electromagnetic field. Because the problem is highly nonlinear, numerical simulations are compulsory. Two-dimensional neutron star magnetospheres have been computed in the force-free regime two decades ago starting with the aligned FFE case [32] and followed several years later by the general three-dimensional oblique cases by [33]. Since then, these results have been retrieved by several other authors using different numerical approaches like finite difference/finite volume methods [34,35,36] or pseudo-spectral methods [37,38,39]. Even a combined spectral/discontinuous Galerkin method has been tried including general-relativistic effects for a monopole [40] or a dipole [41]. Some extension to dissipative magnetospheres was undertaken by [39,42,43] assuming an ad hoc prescription for the dissipation.

Here we show three models of pulsar magnetosphere for an oblique rotator with obliquity (angle between rotation axis and magnetic axis) $\chi =\{0^{\circ },{30}^{\circ },{60}^{\circ },{90}^{\circ }\}$, namely the vacuum, the force-free and the radiative cases. Simulations are performed in the observer inertial frame but using the corotating coordinate system leading to Maxwell equations written as Equation (15). This particular frame ensures that the solution relaxes to a time independent solution where the current sheet remains at a fixed position in space in order to ease its location for subsequent purposes. In other words, the time derivatives in Equation (15) must vanish when the solution becomes stationary.

The three models correspond to three prescriptions for the electric current density $\mathbf{j}$. In vacuum, for the Deutsch solution it is obviously $\mathbf{j}=\mathbf{0}$. This is our reference solution for checking our algorithm and accuracy of the computed solution. Simulations are performed using our pseudo-spectral Maxwell solver explained in depth in [37]. Before discussing our new results, we remind the essential features of our pseudo-spectral code in the following paragraph.

5.1. Numerical Schemes

Spectral and pseudo-spectral numerical schemes convert a system of partial differential equations (PDE) into a larger system of ordinary differential equations (ODE) much easier to integrate numerically with standard ODE integration techniques like the explicit Runge–Kutta and Adams-Bashforth schemes. See [44] for a detailed review of these techniques. We emphasize that spectral methods do not approximate the equations of the problem but the solution itself. Therefore, the numerical problem exactly reflects the mathematical problem with the same boundary conditions which need to be properly imposed without any under or over-determinacy. Please note that finite volume/finite difference codes are prone to large (with respect to spectral codes) diffusion/dissipation and are therefore able to damp boundary conditions that are not exactly identical to the mathematical problem making it analytically an ill-posed problem (mathematically speaking not from a numerical point of view). Spectral methods are primarily dealing with expansion coefficients of the unknown quantities not their value themselves at the grid points. This expansion possesses the great advantage of removing singularities of differential operators like the gradient, the divergence and the curl in spherical coordinates along the polar axis. We use this flexibility to solve Maxwell equations in polar spherical coordinates $(r,\theta ,\phi )$ with no special care about the polar axis. Boundary conditions on the stellar surface can thus be properly and exactly imposed as required by the original mathematical problem.

Specifically, the components of the electromagnetic field are expanded onto a real Fourier–Legendre–Chebyshev basis. The azimuthal dependence is expanded into a standard Fourier series in $cos\left(m\phi \right)$ and $sin\left(m\phi \right)$ whereas the latitude is expanded into Legendre functions $P_{\ell }^m\left(\theta \right)$ where ℓ and m are integers related to the spherical harmonics $Y_{\ell ,m}(\theta ,\phi )$ [45]. The radial part is expanded into Chebyshev polynomials $T_n\left(x\left(r\right)\right)$ where $r\in [R_1,R_2]$ is mapped into the normalized range $x\in [-1,1]$ by a linear transformation. The straightforward implementation of this mapping accumulates the discrete grid found from the Chebyshev-Gauss-Lobatto points unevenly near the boundary points where the resolution becomes prohibitively high. The constrain on the time step is therefore to severe. To distribute more evenly the grid points, we use the Kozloff/Tal–Ezer mapping [46]. See also [47] for a similar implementation of this technique for axisymmetric neutron star magnetospheres. Derivatives are computed in the Fourier–Legendre–Chebyshev space by simple algebraic operations, instead of pure function derivatives, and then transformed back to real space on grid points. The outer boundary conditions are outgoing waves with a sponge layer absorbing spurious reflections. The inner boundary conditions enforce the tangential part of the electric field and the normal component of the magnetic field at the stellar surface. To keep a mathematically well-posed problem, we employ the characteristic compatibility method described in [44]. Time integration is performed via a standard third order Runge–Kutta scheme.

Spectral methods are known to converge to the exact solution faster than finite difference or finite volume schemes for sufficiently smooth problems without discontinuities. They require less resolution for the same accuracy [48]. Because spectral methods rely on Fourier-like series expansions, they are also sensitive to the Gibbs phenomenon [45], spoiling the solution with overshoot possibly leading to unphysical quantities like negative densities or pressures. In our strong electromagnetic field limit; however, no positivity constrain is required for the unknown field. However, in order to stabilize the algorithm, tending to put more and more energy into small scales because of the Gibbs effects, we need to filter the highest frequencies by applying for instance an exponential filter damping the highest order coefficients in the Fourier–Legendre–Chebyshev expansion. Eventually, we check a posteriori that the simulation has converged to the desired solution to good accuracy by performing a resolution analysis, meaning that increasing by a factor two the grid resolution in each direction, the solution does not significantly change. We found that for the simulations shown below, a resolution $N_r\times N_{\theta }\times N_{\phi }=257\times 32\times 64$ already gave reasonable results. We checked on a few cases that increasing by a factor 2 the resolution in all directions did not change the results (but drastically increased the computational time on a single core). Consequently, we adopted a resolution of $N_r\times N_{\theta }\times N_{\phi }=257\times 32\times 64$ for accurate and converged results. In the special case of an aligned rotator, the Gibbs phenomenon is strongest. We had to resort to higher resolution of $N_r\times N_{\theta }\times N_{\phi }=513\times 64\times 1$ for accurate and converged results.

5.2. Force-Free Magnetospheres

In the force-free regime where $\mathbf{E}·\mathbf{B}=0$ and $E<cB$, the electric current density is uniquely defined by [49]

$$\mathbf{j}={\rho }_{\mathrm{e}}\frac{\mathbf{E}\wedge \mathbf{B}}{B^2}+\frac{\mathbf{B}·\nabla \times \mathbf{B}/{\mu }_0-{\epsilon }_0\mathbf{E}·\nabla \times \mathbf{E}}{B^2}\mathbf{B}.$$

(31)

This current is decomposed into an electric drift part, first term on the right hand side, depending only on the total electric charge density ${\rho }_{\mathrm{e}}$, and a part along the magnetic field that is not constrained but deduced a posteriori from the simulation output. Because all particles drift with the same velocity, contribution to the drift part of the electric current arises solely from the non-neutrality of the plasma, meaning ${\rho }_{\mathrm{e}}\ne 0$.

The electric drift speed must remain strictly less than the speed of light. If the condition $E<cB$ is violates, force-free breaks down and the plasma becomes dissipative. In force-free simulations however, we enforce by hand the condition $E<cB$ everywhere in space in order to stay in the sub-relativistic drift speed limit.

An example of magnetic field lines in the equatorial plane for an orthogonal rotator with $\chi ={90}^{\circ }$ is shown in Figure 1 with $R/r_{\mathrm{L}}=0.2$. The force-free regime tries to put the field lines out of the current sheet which becomes singular, but due to numerical resistivity, field lines also tend to close by crossing this singular surface.

5.3. Radiative Magnetospheres

For the radiative magnetosphere, we introduce an additional free parameter represented by the pair multiplicity factor $\kappa$ such that the electric current derived from the Aristotelian electrodynamics becomes

$$\mathbf{j}={\rho }_{\mathrm{e}}\frac{\mathbf{E}\wedge \mathbf{B}}{E_0^2/c^2+B^2}+\left|{\rho }_{\mathrm{e}}\right|(1+2\kappa )\frac{E_0\mathbf{E}/c^2+B_0\mathbf{B}}{E_0^2/c^2+B^2}.$$

(32)

It is decomposed into a $\mathbf{E}\wedge \mathbf{B}$ drift similar to force-free but without the additional constraint $E<cB$ and a part along $\mathbf{E}$ and $\mathbf{B}$ which reduces in the drift frame to a motion along the common direction of ${\mathbf{E}}^{\prime }$ and ${\mathbf{B}}^{\prime }$.

Figure 1 shows some field lines in the equatorial plane for the orthogonal rotator with $\chi ={90}^{\circ }$ in the radiative regime with pair multiplicity $\kappa =\{0,1,2,5\}$. In the most dissipative case corresponding to $\kappa =0$, field lines cross the current sheet at smaller distances compared to less dissipative cases with $\kappa =2$ or $\kappa =5$.

Figure 2 shows the associated radial dependence of the Poynting flux for force-free and radiative cases with $\kappa =\{0,1,2,5\}$. The radiative magnetosphere dissipates a small fraction of the Poynting flux into particle acceleration and radiation, most efficiently when $\kappa =0$, corresponding to a charge separated plasma. Increasing the pair multiplicity factor $\kappa$ to higher values shifts the radiative model towards the force-free limit. In the aligned case, the decrease in Poynting flux is abrupt right at the light cylinder. It is most prominent for $\kappa =0$. However, due to the intrinsic dissipation of our algorithm, even in the FFE case there some Poynting flux dissipation is observed. This is due to the infinitely thin current sheet with discontinuous toroidal magnetic field that is smeared by our spectral methods (Gibbs phenomenon). The situation improves for oblique cases as the displacement current takes over some fraction of the electric current within the sheet.

In Figure 3 we show the Poynting flux crossing the light cylinder for force-free and radiative cases with $\kappa =\{0,1,2,5\}$ and depending on the inclination angle $\chi$. All cases can be fitted with a single formal expression summarized as

$$L=L_{\perp }(a+b{sin}^2\chi )$$

(33)

with different coefficients depending on the regime considered. The fitted values extracted from the numerical simulations are listed in

Table 1. Fitting coefficients a and b for the spindown luminosity as fitted in Equation (33).

Regime	a	b
FFE	1.42	1.73
$\kappa =0$	1.36	1.75
$\kappa =1$	1.39	1.76
$\kappa =2$	1.40	1.74
$\kappa =5$	1.42	1.73

. The most dissipative case $\kappa =0$ slightly decreases the Poynting flux for the aligned rotator already inside the light cylinder. The decrease is accurately quantified by the fitting parameter a. The FFE normalized Poynting flux is 1.42 whereas for the radiative $\kappa =0$ case it is 1.36. The fitting parameter b seems less dependent to the regime considered. The aligned rotator also shows the most prominent gradual decrease in the Poynting flux with respect to distance. Dissipation starts at the light cylinder but goes on at several light cylinder radii. For oblique rotators, the slope of this radial decrease slowly diminishes, becoming negligible for the orthogonal rotator.

In all regimes, the electromagnetic fluxes are very similar while inside the light cylinder. The discrepancies occur outside the light cylinder, in regions where the electric field is dominant and not fully screened by the plasma because of the too low pair multiplicity. A corotative ideal and dissipationlessness magnetosphere inside the star is therefore a good approximation, whereas outside, efficient dissipation sets in right at the light cylinder, around the current sheet.

Figure 4 shows a summary of the Poynting flux crossing the light cylinder (larger markers) and crossing a sphere of radius $4r_{\mathrm{L}}$ (smaller markers) for oblique rotators in force-free and radiative regimes. The dissipation going on at large distances is most visible for the aligned rotator with green triangles.

Some fraction of the electromagnetic flux goes into particle acceleration and radiation. Quantitatively, this dissipation of the electromagnetic energy is computed as a work done on the plasma such that

$$\mathbf{j}·\mathbf{E}=|{\rho }_e|(1+2\kappa )c{E}_0\ge 0.$$

(34)

This dissipation rate, for $\kappa =\{0,1,2,5\}$, is shown in Figure 5 on a log scale. It shows the location of largest dissipation for an orthogonal rotator according to the dissipation rate controlled by $\kappa$. Poynting flux goes into particle acceleration and radiation mainly outside the light cylinder along the current sheet starting from the Y-point. We expect therefore gamma-rays to be produced along this sheet, emitting pulses at the neutron star rotation frequency. Such models have already been put forward and known as the striped wind. See for instance [50] for the production of high-energy emission and [51] for demonstrating the pulsation.

We conclude that energy conversion occurs mainly around the current sheet. Within the light cylinder, the electric field is always less than the magnetic field $E<cB$. Therefore, the force-free condition can be maintained without resorting to artificial damping of $\mathbf{E}$. However, dissipation sets in right at the light cylinder, where magnetic field lines start to cross the light cylinder. The dissipation region follows a spiral pattern with decreasing amplitude with distance from the star. These new simulations offer for the first time a fully self-consistent description of a dissipative and radiative magnetosphere, where feedback between plasma flow, particle radiation and electromagnetic field is included. Please note that emission occurs only along the current sheet outside the light cylinder. This conclusion supports the idea of the striped wind model introduced by [52] and by [53]. It also explains pulsed high-energy emission from gamma-ray pulsars as demonstrated by [51,54,55].

Contrary to the vacuum case, by construction, particles cannot move faster than the speed of light, even if the corotating velocity Equation (22) is used. This is because the magnetic field is now sufficiently bent to counterbalance the effect of adding the corotation velocity given by Equation (14). Note also that radiative magnetospheres presented in our study do not tend to the vacuum solution when the pair multiplicity vanishes $\kappa =0$ because inside the light cylinder we enforce force-free conditions by construction.

Dissipative losses in the current sheet, also called striped wind, have also been proposed by other authors. For instance, ref. [56] found a new standard solution for the aligned rotator, free of separatrix current layer within the light cylinder. Dissipation occurs only in the equatorial current sheet where acceleration and radiation of particle is allowed. They found an increase of 23% of the spindown with respect to $L_{\perp }$, 40% of which goes into the current sheet dissipation. In our solutions, we found a spindown increase from 36% to 42% depending on the pair multiplicity, see Table 1. The crux of the matter is the microphysical description of this current sheet that conditions the whole magnetospheric solution. To prescribe the electric current in this sheet, ref. [56] assumed a null-like current everywhere, a prescription which is questionable. Moreover, 60% of the magnetic flux crossing the light cylinder opens up to infinity. It is not clear how this percentage is controlled by the solution. A better solution would get all magnetic flux dissipated sooner or later in the equatorial current sheet. Ref. [57] used another approach, performing Particle In Cell (PIC) simulations of pulsar magnetospheres. Here the sensitive parameter is the unconstrained pair injection rate, from the surface or from the whole magnetosphere. The stationary solution crucially depends on this injection mechanisms, going from an electrosphere to an almost force-free magnetosphere. They found that less than 15% of the Poynting flux is dissipated within $2r_{\mathrm{L}}$. It is not clear how much additional decrease is expected if the solution would have been computed to larger distances. A partial answer is given in [58] where the dissipation is as high as 35% at $5r_{\mathrm{L}}$. Comparing both models is difficult because they are not performed with the same set up. The most critical variable being the pair multiplicity which is not fixed by the user and not easily controlled. We showed that $\kappa$ strongly affects the asymptotic large distance dissipation in the axisymmetric case. These different approaches can only be reconciled in light of the pair content within the magnetosphere.

To summarize, all results performed with different numerical codes and different assumptions demonstrated that the magnetosphere relaxes automatically to a state where corotation with the star is enforce by the electric current prescription. However, while this picture is simple and easily understood, nothing forbids solutions with differentially rotating plasmas. Such solutions are discussed in the next section.

6. Differentially Rotating Magnetospheres

The neutron star magnetosphere is often described as perfectly corotating with the star, dragged by the electromagnetic field to enforce strict corotation as in the simulations performed in the previous section. However, it is well known from Ferraro isorotation law [59] that to keep corotation, the plasma must be connected magnetically to the star everywhere in the magnetosphere. If some vacuum gaps exist between the surface and the magnetospheric plasma, corotation becomes impossible. The plasma around the equatorial plane will start to rotate differentially, leading to a much complexer variety of physical processes like non-neutral plasma instabilities [60] and efficient particle diffusion across magnetic field lines [61,62].

A perfectly corotating magnetospheric plasma, as simple as it could be, does not represent a realistic pulsar magnetosphere. Differential rotation or lagging of particle motion is permitted when vacuum gaps are allowed. Indeed, in the electrospheric solution found by [63] and detailed by [64] for an aligned rotator, the domes are corotating because magnetically connected to the star but the equatorial disc over-rotates with respect to the star. This differential rotation is induced by a charge density given for an aligned rotator by

$$\rho =-{\epsilon }_0(2\Omega ·\mathbf{B}+r^2{B}^2{\Omega }^{\prime }\left(a\right))$$

(35)

where a is the magnetic flux function. The quantity ${\rho }_{\mathrm{gj}}=-2{\epsilon }_0\Omega ·\mathbf{B}$ is usually referred as the Goldreich–Julian charge density [8]. The charge density $\rho$ in absolute value is much higher than the one required for corotation [65] because ${\Omega }^{\prime }\left(a\right)>0$. Therefore, the location of the light cylinder shifts nearer towards the surface. Actually it is no more a cylinder but an azimuthally symmetric surface. In the general case, it is more suitable to speak about a light-surface rather than about a light cylinder.

This new structure has profound consequence on the secular evolution of the magnetosphere. Indeed, ref. [66] showed that the equatorial disc is unstable with respect to the diocotron instability. A quasilinear theory has been developed by [67] demonstrating the possibility to transport charges across magnetic field lines. This is of paramount importance for the magnetosphere. Later [68] investigated through electrostatic PIC simulations the fully non-linear evolution of the diocotron instability. He found that particles are transported radially outwards in the equatorial plane when the system is feed with fresh electron/positron pairs from the innermost part of the magnetosphere (produced for instance by magnetic photon absorption or photon-photon collisions). However, Pétri [69] proved that the relativistic rotation tends to stabilize the diocotron instability. Moreover, particle inertia becomes significant close to the light cylinder, but the instability still survives, switching to the magnetron case [70].

The presence of such instabilities destroys the picture of a stationary and corotating magnetosphere. By nature, the pulsar electrodynamics is non-stationary as can be witnessed from the highly erratic emission feature of single radio pulse profiles [71]. This is inherent to the relativistic plasma flow and to the pair creation process occurring close to the stellar surface and/or close to the light cylinder. This remark leads us to the last topic concerning radiative signatures within the magnetosphere.

7. Radiation from the Magnetosphere

Pulsars show a broadband emission from radio through optical up to X-rays and gamma-rays. Photons are produced within the magnetosphere and wind. They must indirectly carry information about their production site, therefore showing an imprint of relativistic rotation if produced close to the light cylinder.

The location of high-energy emission from gamma-ray pulsars is poorly constrained by observations. Several competing models interpret the measurements within the magnetosphere or wind, with either curvature, synchrotron or inverse Compton radiation. So far, there is little hope to get a clear insight about the magnetosphere from these observations. Much more interesting in our view is the wealth of data in radio pulse profiles and their associated polarization feature. Radio emission height in normal radio pulsars with slow periods, larger than 100 ms, is well constrained to lie at altitudes around several hundreds of kilometers [72]. This is deduced from the shift in polarization angle traverse with respect to the middle of the pulse profile. This shift $\Delta \varphi$ is explained in the framework of aberration/retardation effects of photons propagating in the magnetosphere [73] and amounts to

$$\Delta \varphi \approx 4r/r_{\mathrm{L}}.$$

(36)

This shift cannot be explained by emitting a photon along a field line in the corotating frame and then using aberration formulas because the pulse profile would be subject to the same shift in phase as the polarization angle traverse and therefore cancelling the possible time delay between the middle of the pulse profile and the steepest gradient in the polarization angle. The only way to correctly catch this shift, which is a well-defined fact seen in many observations of radio pulsars, is through photons emitted along field lines dragged into corotation as seen in the observer inertial frame. This corotation velocity leads naturally to a time lag between the middle of the pulse profile and the polarization angle inflexion point. Nevertheless, the magnetic field topology has to adjust to compensate for the corotation velocity in order to keep the particle velocity less than the speed of light. This is not always possible outside the light cylinder as already explained in Section 4.2. The approximate estimate given in Equation (36) has been checked for off-center and rotating dipoles by [74]. It is, therefore, a very robust result, sharply constraining the radio emission heights.

Consequently, the view presented in Section 4.1 must be rejected because it cannot reproduce aberration/retardation effects in pulsar radio polarization observations. However, the particle motion described in Section 4.2 seems more appropriate. However, the best choice in our view is represented in Section 4.3 where trajectories are computed according to the full electromagnetic field taking into account radiative effects. The latter option gives the simplest plausible scenario where radiation reaction acts efficiently and self-consistently backwards onto the particle motion and onto the electromagnetic field.

8. Discussion

Pulsar magnetospheres have been extensively computed for stellar centered dipolar fields in special relativity. However, there is increasing evidence for off-center or even multipolar components anchored in the neutron star crust. Indeed, joined modelling of pulsed radio emission and thermal X-rays from the polar cap hot spots requires decentered dipoles [74]. Moreover, detailed investigations of X-ray light curves of millisecond pulsars also favors off-center and quadrupolar components according to recent observations from NICER [75]. Nevertheless, as shown by [76], we do not expect drastic changes in the spindown luminosity and magnetic field structure outside the light cylinder for slowly rotating pulsars with period $P>100$ ms. In the case of radiative magnetospheres, we expect a similar trend because radiation occurs outside the light cylinder where the multipolar components have sufficiently decreased to become negligible. Indeed, a dipole magnetic field decreases like $1/r^3$ whereas a multipole of order ℓ decreases like $1/r^{\ell +2}$. Therefore a multipole of strength $B_{\mathrm{m}}$ at the surface contributes only a ratio $(B_{\mathrm{m}}/B_{\mathrm{d}}){(R/r_{\mathrm{L}})}^{\ell }$ compared to a dipole of strength $B_{\mathrm{d}}$ at the surface. We can draw the same conclusions for general-relativistic radiative magnetospheres. Indeed, for force-free magnetospheres, ref. [41] already showed that qualitatively the picture does not vary and that the spindown luminosity scales like the magnetic field strength at the light cylinder. Extrapolating to the radiative case, general-relativistic effects remain very weak outside the light cylinder for any pulsar, millisecond or second and the overall picture discussed above remains valid.

9. Conclusions

Rotation in neutron stars plays a central role to sustain their electromagnetic activity of particle acceleration, pair creation and the subsequent broadband radiation from radio wavelengths to very high energies. A corotating magnetosphere is often used as a good approximation to describe its electrodynamics. However, due to their fast rotation, relativistic speeds are reached already close to the stellar surface, around the light cylinder. This hypothetical surface separates the inner magnetosphere from the wind. We showed that the transition zone between the quasi-static magnetosphere and the wave zone in the wind is difficult to treat satisfactorily and smoothly because of the absence of a physical frame rotating with the star outside this light cylinder. Several prescriptions where exposed, starting from different assumptions. We also showed that pulsar radio observations give some hint about promising paths to follow particle trajectories and their emission.

Clearly, a deeper and better understanding of the neutron star electrodynamics is required to faithfully explain the wealth of data about their radiation. Neutron stars are one of the only macroscopic objects that produce relativistic rotation on a length scale of the order of Earth radius about several hundreds to several thousands of kilometers. Investigating jointly theoretical models and observational facts will reveal relativistic rotation effects in astrophysics in an unusual way, hoping to solve half a century of mystery about their functioning.