RESISTIVE SOLUTIONS FOR PULSAR MAGNETOSPHERES

Pulsar magnetospheres are filled with plasma, and the presence of plasma affects all magnetospheric properties. Although this fact has been well appreciated since the early days of pulsar research, the ability to quantitatively model plasma effects has emerged mainly in the last decade. Generally speaking, the models of pulsar magnetospheres can be classified according to the amount of plasma supply they assume. At one extreme is the vacuum magnetosphere, which is modeled with the magnetic field of a spinning inclined dipole in vacuum (Deutsch 1955). As this solution has no plasma, it has no possibility of producing any pulsar-like emission. However, the fact that this field is described by an analytic formula has made it the most widely used framework for calculating pulsar properties. For example, the characteristic spin-down energy loss which yielded the fiducial pulsar field strength of 10¹² G or the polarization sweep of the rotating-vector model (Radhakrishnan & Cooke 1969) are guided by the vacuum field shape. Slot-gap and outer-gap models of gamma-ray emission from pulsars are also based on this field (Cheng et al. 1986; Dyks et al. 2004). The next order of approximation is models that allow plasma emission from the surface of the star. These include the original charge-separated model of Goldreich & Julian (1969, hereafter GJ) and the space charge-limited models with pair formation (e.g., Ruderman & Sutherland 1975; Arons & Scharlemann 1979). These models envision both the regions where plasma shorts out the accelerating electric fields and the vacuum-like regions where acceleration is present (so-called gaps). This approach allows for more realism in studying plasma creation and acceleration, but at the price of being decoupled from the global structure of the magnetosphere, as the modification of the vacuum field due to plasma currents is typically not included. Finally, the models that concentrate on the global magnetospheric properties assume that abundant plasma exists throughout the magnetosphere and in the wind. This plasma shorts out the accelerating electric fields and provides the corotation of field lines with the star. Such models include the relativistic magnetohydrodynamic (MHD) description of the magnetosphere and its limit for low-inertia magnetically dominated plasmas, the "force-free" description. Which of these different regimes is applicable to real pulsars may ultimately depend on the physics of plasma supply in the magnetosphere.

Currently, quantitative solutions of the global magnetospheric structure exist only for the vacuum limit (Deutsch 1955) and for the limit of abundant plasma in force-free electrodynamics (see Contopoulos et al. 1999; Gruzinov 2005; Timokhin 2006; McKinney 2006; Komissarov 2006 for aligned rotators, and Spitkovsky 2006, hereafter S06; Contopoulos & Kalapotharakos 2010 for pulsars with arbitrary inclinations).³ The real pulsar magnetosphere is likely operating somewhere in between these limits, with various accelerating gaps, regions of pair production, and strong current sheets likely causing local violations of the ideal MHD constraint, ${\bm E}\,{\bm \cdot}\, {\bm B}=0$. Knowing the structure of the magnetosphere, including such non-ideal effects, would be very useful for calculating the properties of pulsar emission. Indeed, currently the ideal force-free models that include the back-reaction of plasma currents on the field structure lack any accelerating fields by construction, and thus cannot be used to directly predict the spectra of gamma-ray radiation observed by Fermi GST.

One way to reintroduce accelerating electric fields in the magnetosphere is to allow finite resistivity of the plasma. Several formulations of resistive force-free equations have been proposed, most notably by Lyutikov (2003) and Gruzinov (2007, 2008). In resistive MHD, the Ohm's law can be unambiguously defined in the proper frame of the fluid by relating the current in that frame to the electric field through ${\bm j}_{\rm fluid}=\sigma {\bm E}_{\rm fluid}$, where σ is plasma conductivity (see, e.g., Lichnerowicz 1967; Palenzuela et al. 2009). In the force-free system, however, the fluid velocity along the magnetic field is unknown, and only the transverse velocity can be obtained from the electromagnetic fields. This introduces some freedom in prescribing the Ohm's law in the force-free picture. In the prescription proposed by Lyutikov (2003), the parallel velocity along the field was taken to be zero. In "Strong-Field Electrodynamics" (hereafter SFE; Gruzinov 2007, 2008), the Ohm's law was formulated in the frame that moves along the field lines with such speed that the charge density in that frame vanishes. Despite the fact that such a frame formally exists only in space-like regions, SFE prescription appears to give smooth numerical solutions throughout the magnetosphere (Gruzinov 2011), and, most importantly, does not explode in current sheets where magnetic fields can go through zero. Both formulations are arbitrary, however, because the real fluid velocity does not have to follow either frame assumption. As an additional constraint, it is useful, therefore, to construct a resistive formulation that reproduces physical solutions expected at the extremes of very large and very small conductivity of the plasma, namely, the force-free limit (σ → ∞) and the limit of vacuum electromagnetism (σ → 0). Below we describe a resistive prescription that generalizes these schemes and combines the correct limiting behavior of Lyutikov's scheme with the current-sheet stability of Gruzinov's formulation. We then apply this prescription to numerically calculate the structure of resistive magnetospheres in pulsars.

In the resistive force-free picture of the pulsar magnetosphere, the magnetized neutron star is thought to be surrounded by an abundant massless plasma with finite conductivity, so that not all accelerating fields are shorted out. For simplicity and as a proof of principle, we only consider the unrealistic case of constant conductivity throughout space. More complicated prescriptions will be studied elsewhere. Our finding is that using our formulation of the resistive force-free electrodynamics we can construct a family of magnetospheres that smoothly transition from the Deutsch vacuum solution to the ideal force-free magnetosphere as the conductivity of the plasma is increased. Such intermediate magnetospheres possess interesting properties that we discuss in this paper. We study the variation of the spin-down power with magnetic inclination angle as a function of plasma conductivity and relate it to physical conditions such as the effective potential drop on the open field lines. In Section 2, we discuss the derivation of the resistive force-free prescription. In Section 3, we describe our code for solving the equations and present sample magnetospheric solutions. Discussion and potential applications to pulsar physics are in Section 4.

We describe here our prescription for resistive current as used in our numerical code. Ohm's law is defined in the fluid rest frame. In this frame the electric and magnetic fields are parallel; otherwise, there would be a particle drift across magnetic field lines. The laboratory frame can be connected to the fluid rest frame through two boosts. One boost is in the E × B direction in the lab frame and brings the electric field parallel to the magnetic field. The other boost is along the parallel electric and magnetic fields and transforms the four-current to the fluid rest frame while leaving the electric and magnetic fields parallel. We choose a simple Ohm's law to relate the current and electric field in the fluid frame: j_fluid = σE_fluid, where σ is the plasma conductivity. More general formulations of the Ohm's law that take into account time-dependent currents, inertial effects, pressure, and the Hall effect have been derived (see, e.g., Meier 2004). It is likely safe to ignore most of these effects in the bulk of the strongly magnetized cold pair plasma, as envisioned for typical pulsar magnetospheres; however, singular current sheets may require a more elaborate treatment. In this work, we use a constant uniform conductivity throughout the domain and restrict ourselves to the simple Ohm's law to elucidate the basic physics. Boosting back to the laboratory frame, the current vector can be expressed as

$$\begin{eqnarray} \bm{j}&=& \frac{\rho c\bm{E}\times \bm{B}}{B^2+E^2_0} \nonumber \\ && +\,\frac{\left(-\beta _{||}\rho c+\sqrt{\frac{B^2+E^2_0}{B^2_0+E^2_0}\left(1-\beta _{||}^2\right)}\sigma E_0\right)(B_0\bm{B}+E_0\bm{E})}{B^2+E^2_0}.\nonumber\\ \end{eqnarray} \tag{ 1 }$$

See Appendix A for a full derivation of this expression and the subsequent limits that we discuss below. Here, β_|| is the magnitude of the boost along the parallel electric and magnetic fields to the fluid rest frame, ρ is the charge density in the laboratory frame, and E₀ and B₀ are the magnitudes of the electric and magnetic fields in all frames in which they are parallel. E₀ and B₀ are defined by the expressions (Gruzinov 2007, 2008)

$$\begin{eqnarray} B_0^2&=\frac{{\bf B}^2-{\bf E}^2+\sqrt{({\bf B}^2-{\bf E}^2)^2+4(\bm{E}\, {\bm \cdot}\, \bm{B})^2}}{2}, \nonumber \\ E_0&=\sqrt{B_0^2-{\bf B}^2+{\bf E}^2}, \nonumber \\ B_0&={\rm sign}(\bm{E}\, {\bm \cdot}\, \bm{B})\sqrt{B^2_0}, \end{eqnarray} \tag{ 2 }$$

where we allow B₀ to be positive or negative, depending on whether the magnetic field is aligned or antialigned with the electric field. Note that in addition to the advection of charge in the E × B direction, there is conduction current along both the laboratory frame magnetic field and the electric field (the last term in Equation (1)). There is still a principal ambiguity here in defining the fluid rest frame, as we have not yet specified the magnitude of the parallel boost β_||. In fact, this speed cannot be obtained from a purely electrodynamic standpoint without the inclusion of gas dynamics, so a suitable choice has to be made.

One particularly interesting choice of fluid frame is the slowest moving frame that has electric fields parallel to magnetic fields, i.e., the frame with β_|| → 0. If particles start out on the stellar surface with non-relativistic velocity, then the fluid essentially satisfies this condition in the high-conductivity limit (see, e.g., Contopoulos & Kazanas 2002; Tchekhovskoy et al. 2009; Beskin 2010). In this minimal velocity limit the current becomes

$$\begin{equation} \bm{j}=\frac{\rho c\bm{E}\times \bm{B}+\sqrt{\frac{B^2+E^2_0}{B^2_0+E^2_0}}\sigma E_0(B_0\bm{B}+E_0\bm{E})}{B^2+E^2_0}. \end{equation} \tag{ 3 }$$

This form of the current has the especially useful property that in the limit of vanishing conductivity the conduction current along the electric and magnetic fields goes to zero. Such a magnetosphere does not develop space charge due to the flow of charge from the star along the field lines. Thus, the current in the E × B direction also vanishes. This is exactly what we expect in the vacuum limit. Since we are interested in producing a resistive transition from ideal force free to vacuum, we use this form of the current in our numerical investigation. Lyutikov (2003) derived an expression for resistive current in magnetically dominated plasma, choosing the same minimal velocity frame for the formulation of the Ohm's law as in our derivation. However, his boost in the E × B direction to the fluid frame contains an error and does not bring the electric and magnetic fields parallel to one another. As a result, our final expressions differ.

We can alternatively express the minimal velocity current as

$$\begin{equation} \bm{j}=\rho \bm{v}+\sigma \bm{E}_{\rm fluid}, \end{equation} \tag{ 4 }$$

where the fluid velocity v = c(E × B)/(B² + E²₀) is the generalized drift velocity, E_fluid = γ(E + v × B), and γ = (1 − v²/c²)^−1/2. The drift velocity contains the term E²₀ in the denominator to account for the non-zero electric field in the fluid frame (see Equation (A14)). The presence of this term allows the current to remain nonsingular in current sheets, where the magnetic field vanishes.

It is instructive to consider other limits for the current. Gruzinov (2007, 2008) derived an alternate limit of the general current in Equation (1), known as SFE. He postulated that Ohm's law should be formulated in the frame where charge density vanishes, which gives parallel boost

$$\begin{equation} \beta _{||}=\frac{-\rho }{\big(\gamma ^2_x\sigma ^2E^2_0\big/c^2+\rho ^2\big)^{1/2}} \end{equation} \tag{ 5 }$$

and current

$$\begin{equation} \bm{j}=\frac{\rho c\bm{E}\times \bm{B}+\left(\gamma ^2_x\sigma ^2E^2_0+\rho ^2c^2\right)^{1/2}(B_0\bm{B}+E_0\bm{E})}{B^2+E^2_0}, \end{equation} \tag{ 6 }$$

where γ²_x ≡ (B² + E²₀)/(B²₀ + E²₀). We see no special reason as to why the fluid frame charge density must vanish. Further, this current does not reduce to vacuum as the conductivity drops to zero, as it does in Equation (3). Consequently, even for vanishing conductivity the SFE solutions resemble ideal non-resistive solutions.

In the limit of infinitely conductive strongly magnetized plasma, which we refer to as the ideal force-free limit, the plasma satisfies transverse force balance, ρE + j × B = 0, and the parallel electric field is shorted out, E₀ → 0. However, as σ → ∞, the product of σE₀ stays finite, and the minimal velocity current becomes (Osherovich & Gliner 1988; Gruzinov 1999; Blandford 2002)

$$\begin{equation} \bm{j}=\frac{\rho c\bm{E}\times \bm{B}}{B^2}+\frac{c}{4\pi }\frac{(\bm{B}\cdot \nabla \times \bm{B}-\bm{E}\cdot \nabla \times \bm{E})\bm{B}}{B^2}. \end{equation} \tag{ 7 }$$

In order to test the influence of global resistivity on the structure of pulsar magnetospheres, we employ a three-dimensional numerical code (see S06) that implements the finite difference time-domain scheme (FDTD; Taflove & Hagness 2005) to evolve electromagnetic fields from Maxwell's equations,

$$\begin{eqnarray} \partial _t \bm{E} &=c \bm{\nabla }\times \bm{B} - 4\pi \bm{j},\quad \partial _t \bm{B} &= -c\bm{\nabla }\times \bm{E}, \end{eqnarray} \tag{ 8 }$$

where the current is given by Equations (2) and (3). Electric and magnetic fields are decentered on the Cartesian Yee mesh, as in standard FDTD. Instead of the leapfrog time integration typically used in FDTD, we employ a third-order Runge–Kutta integrator as in S06. This is done to remove the interpolation in time that would otherwise be needed to bring the E and B fields to the same time step when computing the source current. We run our simulations on a grid of size 1024³ at Courant number 0.5. The central region of our grid is occupied by a conducting spherical star of radius R_*, rotating at angular velocity Ω, with embedded dipole field of magnetic moment μ inclined relative to the rotation axis by angle α (see S06 for more details). The electric field inside the star is set to ensure rigid corotation, and we smooth fields across the stellar boundary as in S06 in order to minimize stair-stepping of a sphere on a Cartesian grid. The spatial resolution of our simulations is such that the light cylinder (LC) radius, R_LC = c/Ω, is resolved with 80 cells, and the simulation box is 12.8 R_LC on a side. The outer boundary condition is periodic, which limits our run time to the light-travel time across the grid. We set R_* = 30 cells, allowing us to resolve the star while keeping the star sufficiently small compared to the LC. We run our code for a range of different values of the dimensionless conductivity, σ/Ω, to produce a transition from the vacuum solution at low conductivity to the ideal force-free solution at high conductivity. We have verified that our solutions are converged with spatial resolution, as well as run sufficiently long so as to reach a steady state in the frame corotating with the pulsar.

In addition to the resistive solutions obtained with the above method, we also consider two special cases of the current prescription that we call "force free" and "vacuum." By force-free solutions we mean "ideal force-free" magnetospheres; however, instead of using the full Equation (7) for the current, our numerical implementation follows the approach in S06, where at each time step we evolve the perpendicular part of the current from Equation (7) and then correct the electric field by removing any electric field component parallel to the magnetic field. Further field limiters are used to enforce E < B in current sheets. As a result, the force-free scheme of S06 has an effective small resistivity because the cleaning of parallel fields can be interpreted as the effect of very large conductivity along field lines. Also, as all parallel fields are cleaned on every step, such a scheme does not have the usual convergence behavior with time step. The resistive scheme using Equation (3) gives us a more well-defined method for handling magnetospheric dissipation. In addition, this resistive scheme does not require special limiters to handle current sheets.⁴ The "vacuum" solutions described below evolve Maxwell's equations without volume charges or currents outside the star, allowing for significant potential drops along field lines.

In Figure 1, we show magnetic field lines in the μ − Ω plane for 60° inclined force-free dipole (panel (a)), resistive dipoles at (σ/Ω)² = 40 (panel (b)), (σ/Ω)² = 4 (panel (c)), (σ/Ω)² = 0.4 (panel (d)), (σ/Ω)² = 0.04 (panel (e)), and vacuum dipole (panel (f)). The field lines are shown after 1.5 turns, sufficiently long to reach a steady state out to several LC radii in the corotating frame. We show only the central five LCs of our simulation. Color represents the out-of-plane magnetic field into (red) and out of (blue) the page. The maximum magnitude of the out-of-plane magnetic field value is different in each panel, with maximum value increasing with conductivity. Table 1 gives the maximum values normalized to the maximum vacuum value.⁵ The maximum values occur at the stellar surface and fall on the blue end of the color table. To improve contrast, the color table shows values up to 30% of the maximum, and a square root stretching is applied to the data. Gradients in color reflect the sum of in-plane components of conduction and displacement currents.

Zoom In Zoom Out Reset image size

Figure 1(a) shows the force-free magnetosphere. The gross features of this solution are the same as were discussed in S06,⁶ Contopoulos & Kalapotharakos (2010), and Bai & Spitkovsky (2010). The polar cap consists of the footpoints of open field lines that extend out to infinity. Field lines with both ends attached to the star form the closed field line region, which extends out to a cylindrical distance of 1R_LC. A large-scale conduction current flows outward from the polar cap along open field lines and returns through the current sheet, current layer, and a fraction of the neighboring open field lines. The current sheets appear as the sharp transition in the out-of-plane magnetic field, starting at the Y-point at the LC. The current layer flows along the boundary of the closed field line region. The conduction current is closed in a circuit by a surface current flowing across the polar cap, a point we will return to when discussing pulsar spin-down.

At (σ/Ω)² = 40, shown in Figure 1(b), the qualitative picture of the magnetosphere looks quite similar to force free. We find that values (σ/Ω)² > 40 give fairly good estimates of highly conducting magnetospheres for R_*/R_LC = 3/8. The dimensionless quantity σ/Ω = (c/Ω)/(c/σ) is the relativistic analog of the Elsasser number, defined as the ratio of the Alfvén radius in a rotating system, v_A/Ω, to the resistive diffusion scale (Elsasser 1946). The above threshold corresponds to resistive diffusion length scale c/σ small compared to the LC, c/σ < 0.16R_LC ≪ R_LC, while satisfying c/σ ≲ R_pc = R_*(R_*/R_LC)^1/2. We have also run simulations with our current prescription at higher conductivity, (σ/Ω)² = 400, and find spin-down luminosities that match the force-free values quite well. Our explicit scheme requires very small time steps at such high conductivities though, and it is perhaps more appropriate to use an implicit–explicit (Palenzuela et al. 2009) or fully implicit scheme to explore this regime. We note that the kink in the current sheet seen in the force-free solution near 2R_LC disappears for resistive solutions even at high conductivities. It is likely then that this kink is a numerical artifact of the cleaning procedure used in our force-free formulation.

Displacement currents play an important role in the highly conducting magnetospheres of Figures 1(a) and (b). The displacement currents in the current sheets actually dominate over conduction currents for the 60° rotator. Displacement currents are not present in the magnetospheres of aligned rotators, but the strength of the displacement current in the current sheets of highly conducting magnetospheres increases monotonically with increasing inclination angle. This is in contrast to the conduction currents, whose strength is independent of inclination angle for fixed conductivity. Both displacement and conduction currents weaken considerably with decreasing conductivity, reflected by smaller gradients in the out-of-plane magnetic field in panels (b) through (e) of Figure 1. Recall that the maximum in the color table corresponds to weaker out-of-plane magnetic field at lower conductivities. The extent of the closed field line region also expands with decreasing conductivity, evident in the same sequence of panels. We further expect the current sheet to thicken with decreasing conductivity. At (σ/Ω)² = 0.04, shown in Figure 1(e), the current sheet should have characteristic width c/σ of order 5R_LC.

There is an obvious increase in the radius of toroidal field sign reversal beyond the LC as one transitions from the force-free to the vacuum magnetosphere. This can be understood from the increase in the winding radius of the spiral pattern for lower conductivities. The winding of the spiral is determined by the characteristic speed of the radial outflow multiplied by the pulsar period. In the force-free case, the radial outflow is given by the radial component of the E × B velocity, which increases with radius so that the Lorentz factor grows linearly with cylindrical radius: $\gamma _{\bm E\times \bm B}=\sqrt{1+(R/R_{\rm LC})^2}$ (Contopoulos & Kazanas 2002). Thus, the force-free outflow approaches the speed of light only asymptotically. In the vacuum case, the "outflow" is always at the speed of light, so we expect the sign transition of the field to occur near the half-wavelength of the vacuum wave, or πR_LC, as is confirmed by Figure 1(f). The force-free case of Figure 1(a) reverses near 2R_LC, indicative of the smaller outflow speed near the LC.

The surface currents flowing across the polar cap of the pulsar exert a spin-down torque on the star. The rotational energy loss is connected to infinity by an outward-directed Poynting flux, which can be thought of as the ejection of toroidal field. We compute the Poynting flux at the LC for each of our simulation runs and show the results in Figure 2. All spin-down values are calculated for R_* = 3/8R_LC. The spin-down curves have been normalized to L₀, defined as 3/2 times the power of the orthogonal vacuum rotator with finite R_* = 3/8R_LC. Two-dimensional axisymmetric calculations for small star tell us that the spin-down of the aligned force-free rotator should be 3/2 times larger than the power of the orthogonal point vacuum dipole, L₁ = 2μ²Ω⁴/3c³ (Contopoulos et al. 1999; Gruzinov 2005; McKinney 2006; S06). The short-dashed line shows the function L/L₀ = 1 + sin ²α, the expected force-free curve in the limit of small star with inclined dipole (S06). The gray band around the force-free curve indicates the uncertainty in the measurement due to boundary effects and numerical dissipation of Poynting flux in the magnetosphere (see Appendix B for thorough discussion). Similar bands are implied but not shown for resistive spin-down curves. The long dashed line shows the analytic vacuum Deutsch field solution (Michel & Li 1999). One of our principal results is that we see a smooth monotonic transition from the force-free spin-down curve to the vacuum spin-down curve for decreasing conductivity.

Zoom In Zoom Out Reset image size

All spin-down curves in Figure 2 show a strong dependence on both inclination angle and conductivity. The radial Poynting flux carrying the spin-down power is proportional to the product of poloidal and toroidal magnetic fields (E_θB_ϕ ∼ ΩRB_pB_ϕ/c), both of which are affected by the strength of displacement and conduction currents in the magnetosphere. One can get a qualitative feel for the relative significance of two contributions by considering the limiting cases. In vacuum, the increase of spin-down with inclination is solely due to rising displacement current. This current is likely responsible for much of the angular dependence of the resistive and force-free solutions. The increase in spin-down with increasing conductivity is due to the additional sweepback and opening of the poloidal field brought on by the increasing conduction current. We find that spin-down luminosity at any conductivity has an angle-dependent term proportional to sin ²α. We parameterize the spin-down curves in Figure 2 with functions f(σ/Ω) and g(σ/Ω) such that L/L₀ = f(σ/Ω) + g(σ/Ω)sin ²α. We find the piecewise linear fit

$$\begin{eqnarray} \frac{L}{L_0}&=&0.3+0.3\log (\sigma /\Omega)^2+1.2\sin ^2\alpha \,, \,\, (\sigma /\Omega)^2>0.4; \nonumber \\ \frac{L}{L_0}&=&0.2+0.08\log (\sigma /\Omega)^2+(1.3+0.2\log (\sigma /\Omega)^2)\sin ^2\alpha, \nonumber \\ &&\times\,0.004<(\sigma /\Omega)^2<0.4. \end{eqnarray} \tag{ 9 }$$

The amplitude of the angular dependence, g(σ/Ω), is constant for (σ/Ω)² > 0.4 and begins to transition to the vacuum value below (σ/Ω)² = 0.4.

We have thus far shown how spin-down luminosity depends on plasma conductivity, but the physical meaning of the conductivity is not entirely clear. It is instructive to reinterpret the conductivity parameter, σ/Ω, in terms of the potential drop along open field lines in the corotating frame. This gives us a handle on the deviation of the magnetosphere from ideal force free, which has vanishing potential drops along field lines. The electromagnetic fields in the frame corotating with the pulsar are obtained via a coordinate transformation from the laboratory frame (Schiff 1939; Grøn 1984):

$$\begin{equation} {{\bm E}^{\prime }}=\bm{E}+\frac{\bm{\Omega }\times \bm{r}}{c}\times \bm{B} \end{equation} \tag{ 10 }$$

and

$$\begin{equation} {{\bm B}^{\prime }}=\bm{B}. \end{equation} \tag{ 11 }$$

Since the fields are steady in the corotating frame, ${\bm \nabla} \times \bm{E^{\prime }}=0$ and the corotating electric field can be written as the gradient of a scalar potential, i.e., $\bm{E^{\prime }}={\bm \nabla} \chi$. Taking the line integral of the corotating electric field along a magnetic field line, l, we find

$$\begin{equation} \Delta \chi =\int _{l} {{\bm E}^{\prime }}\cdot \bm{d}l=\int _{l} \bm{E}\cdot \bm{d}l \equiv V_{\rm drop}. \end{equation} \tag{ 12 }$$

We see that the potential drop along field lines in the corotating frame can be computed directly from the laboratory frame fields. Although in resistive solutions particles will actually drift across magnetic field lines, in addition to accelerating along them, we choose to study the field-aligned potential drop V_drop as a fiducial measure of particle energy gain. In fact, the exact particle trajectories we choose make little difference when estimating potential drops because the electric field is potential in the frame corotating with the pulsar.

Consider field lines starting on the stellar surface in the μ − Ω plane separated by 15° in latitude from pole to equator. For every such field line we integrate the electric field to find the maximum potential drop along each field line. We then determine the field line with the largest overall potential drop for a given magnetic inclination. Integrating field lines separated by 15° increments on the stellar surface is sufficient to give a good estimate of the maximum potential drop along field lines. Figure 3 shows the maximal potential drop as a function of dipole inclination angle for different conductivities. All results have been normalized to the potential drop from the pole to the equator of an aligned vacuum rotator in the laboratory frame, V₀ = |μ|/(R_LCR_*). As our models do not prescribe a high conductivity to the closed field line region, the available accelerating potential is generally limited by the pole-to-equator potential difference, rather than the smaller polar cap potential, V_pc = V₀R_*/R_LC, with a few notable exceptions. Low-conductivity solutions, (σ/Ω)² ≲ 0.04, at high inclination angle, α > 45°, have V_drop scaling intermediate between V₀ and V_pc. The orthogonal rotator drop at (σ/Ω)² = 0.04 scales roughly with V_pc. Using the vacuum Deutsch fields, one can obtain that the aligned vacuum rotator potential drop scales exactly proportional to V₀, whereas the orthogonal vacuum rotator drop scales closer to V_pc. From Figure 3 we see that for (σ/Ω)² > 0.04 the potential drop is roughly independent of inclination angle for fixed σ, i.e., there is a one-to-one map between conductivity and maximal potential drop along a field line originating in the μ − Ω plane for all inclination angles. The maximal potential drop increases from zero for force free (or σ/Ω → ∞) to V_drop/V₀ ∼ 0.2 for (σ/Ω)² = 0.04.

Zoom In Zoom Out Reset image size

The fact that potential drop is independent of inclination angle for (σ/Ω)² > 0.04 allows us to relate luminosity to potential drop in a very simple manner. We present the result here and explain its origin in more detail below. Figure 4 shows the spin-down luminosity as a function of potential drop for inclination angles α = 0°, 30°, 60°, 90°. Spin-down luminosity increases with increasing inclination angle and with decreasing field line potential drops. We fit the spin-down curves in Figure 4 with the linear relation

$$\begin{equation} \frac{L}{L_0}=0.9 \left(1-\frac{V_{\rm drop}}{0.2 V_0}\right)+1.1\sin ^2\alpha \,, \,\, 0<V_{\rm drop}< 0.2 V_0. \end{equation} \tag{ 13 }$$

This formula applies only to the domain for which we have data shown with solid lines in Figure 4, i.e., 0 < V_drop/V₀ < 0.2. A similar spin-down formula was constructed by Contopoulos & Spitkovsky (2006) for the case of a finite gap at the base of the open field lines in a force-free magnetosphere. The spin-down was also found to be linear in the potential drop in the gap but had a different slope and angular dependence of the form (1 − V_drop/V_pc)cos ²α + sin ²α. This angular dependence was not a rigorous derivation, but a prediction based on aligned force-free and orthogonal vacuum limits. We note that our spin-down formula reduces to the force-free spin-down formula from S06 within error bars when V_drop/V₀ = 0.

Zoom In Zoom Out Reset image size

For V_drop/V₀ > 0.2, the dependence of spin-down on the potential drop enters a different regime (see dotted lines in Figure 4). We do not provide a fit here, as potential drop no longer scales with V₀ at high inclination angle. The potential drop along field lines in this regime is due to the reemergence of the vacuum electric fields as conductivity is reduced. There are three contributions to the vacuum electric field which are provided by the quadrupolar surface charge on a finite-radius star and by central and surface monopolar charges. The central monopolar component provides part of the radial electric field needed for corotation of the magnetized stellar interior. If the star initially is uncharged (as is assumed in our simulations), a surface charge of opposite sign compensates the net charge of the interior of the star (see Michel & Li 1999 for a thorough discussion). This surface charge, together with the induced quadrupolar surface charge from corotation, can leave the star and be redistributed throughout the magnetosphere when conductivity is turned on, or when the work function of the surface is low. Such redistribution lowers the maximum available potential drop.

The available monopolar surface charge varies with inclination as cos α (Michel & Li 1999), so the potential drop associated with it disappears for orthogonal rotators. Hence, the value V_drop ∼ 0.2V₀ for the vacuum 90° rotator in Figure 3 represents the quadrupolar contribution to the potential drop. At conductivity (σ/Ω)² ∼ 0.04, the maximal potential drop is close to the orthogonal vacuum rotator value of V_drop ∼ 0.2V₀ at all inclination angles. Qualitatively, magnetospheres with finite conductivity (σ/Ω)² > 0.04 can effectively redistribute the surface charges and maintain the linear relationship between the maximal potential drop and the spin-down luminosity. Below (σ/Ω)² = 0.04, the surface charge redistribution is incomplete (the monopolar surface charge is concentrated closer to the star in the corotating steady state solution) and the curves in Figure 3 gain a tilt more reminiscent of the vacuum solution. The potential drops are not constant with inclination angle at a given conductivity, and the linear relation (13) no longer holds. We have explicitly verified that by the end of the simulation these solutions reached steady state in the corotating frame.

Redistribution of surface charge in near-vacuum magnetospheres also affects the spin-down power. Looking back at Figure 2 and the piecewise linear relations (9), we see that the final 30% of the spin-down transition between force-free and vacuum rotators occurs for (σ/Ω)² < 0.04. Despite the fact that potential is no longer constant with angle for each value of the conductivity in this regime, the transition in spin-down is still continuous and smooth.

We have presented a continuous one-parameter family of pulsar magnetosphere solutions that span the range between the vacuum and force-free limits. Each solution of the family is characterized by the value of the conductivity parameter, σ/Ω, which is related to the maximum potential drop that a test particle can experience as it moves along field lines, V_drop. The zero conductivity limit, σ/Ω → 0, V_drop → V_{drop, vac}, yields the vacuum magnetosphere. It shows a substantial effective potential drop that makes up a significant fraction of the rotation-induced potential difference between the pole and the equator of the star (20%–60% depending on the inclination). The high-conductivity limit, σ/Ω → ∞, yields essentially an ideal force-free magnetosphere. Unlike in vacuum, abundant magnetospheric charges short out any potential differences along field lines, leading to vanishing effective potential drop, V_{drop, ff} → 0. While in the vacuum case nearly all field lines, even those that extend beyond the LC, eventually return to the star and are formally closed, in the ideal force-free case a fraction of field lines open up and reach infinity. Our simulations show that resistive high-σ pulsars spin down at least three times faster than resistive low-σ and vacuum pulsars (for the same value of inclination), in agreement with earlier ideal force-free simulations (S06). Our resistive solutions bridge the gap between the force-free and vacuum limits by having intermediate values of major magnetospheric parameters, such as effective potential drop, 0 < V_drop < V_{drop, vac}, the fraction of open field lines, and the spin-down rate.

Before discussing possible observational implications of these solutions, we note that our approach to modeling resistive pulsar magnetospheres is quite simplistic: we use a form of Ohm's law in which we neglect several terms (Meier 2004) and we assume constant conductivity σ/Ω throughout all space. The finite bulk conductivity used in this work can be thought of as not due to collisions, but due to the insufficient plasma supply in the magnetosphere that cannot short out all the accelerating fields. It is possible that real magnetospheres have a lower (anomalous) conductivity in the current layer and the current sheet than in the rest of the magnetosphere. We might also expect a lower value of conductivity on the open field lines than on the closed field lines. Since our models do not explicitly model the plasma fluid (Section 1), the plasma velocity component along the direction of magnetic field is unconstrained, which necessitates the choice of a frame in which to write the Ohm's law. We choose the so-called minimum velocity frame, which has a number of attractive properties (see Section 2), yet this is not a covariant choice. With these important caveats in mind, let us now for the sake of argument take our model seriously as a description of pulsar magnetosphere and consider the consequences.

Intermittent pulsars offer a unique test bed of pulsar theory (e.g., PSR B1931+24, Kramer et al. 2006 and PSR J1832+0029, Lyne 2009; see also Wang et al. 2007, Zhang et al. 2007, and Timokhin 2010). Such pulsars switch between an "on," radio-loud, state in which they behave like normal radio pulsars, and an "off," radio-quiet, state in which they produce no detectable radio emission. The two intermittent pulsars, for which published data exist, have quite different duty cycles: PSR B1931+24, with a period P ≈ 0.8 s, cycles through the "on"–"off" sequence of states approximately once a month, whereas PSR J1832+0029, with a period P ≈ 0.5 s, kept quiet for nearly two years. The spin-down rate for each of these pulsars is larger in the "on" state than in the "off" state by a factor f_{on → off} ≃ 1.5. Such a substantial difference in spin-down rates suggests that the pulsar magnetosphere undergoes a dramatic reconfiguration as it transitions between the "on" and "off" states, yet such a transition was reported for PSR B1931+24 to take place in just over 10 pulsar periods.

Kramer et al. (2006) propose that in the "on" state, plasma fills the pulsar magnetosphere and supports (poorly understood) plasma processes that produce radio emission. In this picture, the pulsar transitions to the "off" state when, due to some unknown trigger, magnetospheric pair production shuts off: the remaining plasma leaks off the open field lines and the pulsar goes radio-quiet. In this scenario, we expect the pulsar in the "on" state to spin down due to force-free energy losses and in the "off" state due to vacuum dipole losses. However, this picture is ruled out by observations: the force-free spin-down rate is larger than the vacuum spin-down rate by a factor, f_{ff → vac} ⩾ 3 (S06), that is clearly incompatible with the observed value, f_{on → off} ≃ 1.5 (e.g., Beskin & Nokhrina 2007).

Our resistive magnetospheres provide a possible resolution to this problem since they have intermediate spin-down rates between vacuum and force free. Let us associate the "on" state with the force-free magnetosphere and the "off" state with one of our intermediate resistive magnetospheres. According to Equation (13), for the ratio of spin-down powers in these two states to equal f_{on → off}, the potential drop has to be V_drop = V₀(f_{on → off} − 1)(0.9 + 1.1sin ²α)/4.5f_{on → off}. For f_{on → off} ≈ 1.5, the ratio V_drop/V₀ increases monotonically with increasing inclination angle, α, from 0.07, for an aligned rotator (α = 0°), to 0.15, for an orthogonal rotator (α = 90°). Plugging in for V₀, using the measured spin-down parameters of PSR B1931+24 (we associate the observed spin-down rate in the "on" state with the spin-down rate of the force-free magnetosphere), and assuming a neutron star mass, M_* = 1.4 M_☉, we obtain V_drop = 1.5 × 10¹⁶(0.9 + 1.1sin ²α)^1/2[V]. This value significantly exceeds the characteristic value of the polar cap potential drop expected for pulsars (Goldreich & Julian 1969; Arons 2009), so the pulsar should have no difficulty in creating pairs and shorting out the large potential. The reason for this discrepancy is our simplifying assumption that the conductivity is constant throughout the magnetosphere. This assumption causes the potential drop not just across the polar cap but across the whole stellar surface. In real pulsars, we expect that the closed field line region is filled with plasma and is highly conductive. This plasma shorts out parallel electric field in the closed field line region and screens the potential drop everywhere except across the polar cap. The potential drop across the polar cap, V_pc, is ≈(R_*/R_LC) times the pole-to-equator potential drop, V₀, which gives

$$\begin{eqnarray} V_{\rm drop} &\sim V_{\rm pc}\frac{(f_{{\rm on}\rightarrow {\rm off}}-1)}{4.5 f_{{\rm on}\rightarrow {\rm off}}} =3.8\times 10^{12}[{\rm V}], \end{eqnarray} \tag{ 14 }$$

for PSR B1931+24, where we assumed R_* = 10 km.

Addition of resistivity can, in principle, allow the detailed study of pulsar spin-down and the determination of pulsar braking indices. For a spin-down law Ω = Ω(t), the braking index is defined as $n \equiv {\Omega \ddot{\Omega }}/{\dot{\Omega }^2}$. In order to measure a braking index numerically, we run a series of simulations for different values of angular frequency, Ω, and map out the dependence, L(Ω). We find that L(Ω)∝Ω⁴ for both purely force-free and vacuum solutions, which translates into a value of the braking index, n = 3, as we will see below. This is larger than the range of observed values, n ≃ 2.5–2.8 (Livingstone et al. 2007). One way to lower the braking index is to allow the evolution of the extent of the magnetosphere (the radius of the Y-point) to lag behind the outward expansion of the LC due to spin-down (e.g., Contopoulos & Spitkovsky 2006). The extent of the Y-point is controlled by the rate of reconnection at the edge of the magnetosphere. We calculated the braking index in our simulations by comparing the magnetosphere with constant conductivity spun up to different periods, and we do not find significant deviations from n = 3. This rapid spin-up biases our answer, however, because the Y-point extends to the LC within one rotation period. It is not entirely clear that we would obtain the same braking index if we allowed the pulsar period to increase during a single simulation as the pulsar spins down and the LC slowly recedes. The extent of the Y-point may depend on the past history of the pulsar. To do this calculation, we need to increase the pulsar period several times within a single simulation and measure the resulting spin-down luminosity after each increase. We must also ensure that the reconnection at the Y-point is controlled by the resistivity we impose and not by numerical effects. A more accurate braking index calculation would require much longer simulations, and possibly higher spatial resolution near the Y-region.

Another possibility that can affect the braking index is the evolution of conductivity with time. Such time-evolution translates into evolution in pulsar dimensionless luminosity, ℓ = L/L₀ (see Figure 2), and leads to n ≠ 3:

$$\begin{equation} n = \frac{\Omega \ddot{\Omega }}{\dot{\Omega }^2} = 3 - 2 \frac{d\ell /dt}{\ell /\tau }, \end{equation} \tag{ 15 }$$

where $\tau =P/2\dot{P}$ is the pulsar age, and we used the scaling, L₀∝Ω⁴. If dimensionless plasma conductivity does not change in time, σ/Ω = constant, then ℓ = constant and Equation (15) gives n = 3. Purely force-free and purely vacuum pulsars, discussed above, fall into this category. As a pulsar ages, it is natural to expect that its plasma supply depletes and the pulsar becomes progressively more vacuum-like: in other words, both σ/Ω and ℓ decrease in time (see Figure 2). For such pulsars, dℓ/dt < 0, and Equation (15) gives n > 3. Of course, the converse is also true, and the braking index lower than 3 would result if a pulsar becomes progressively more force-free-like, i.e., if pulsar dimensionless luminosity increases in time, dℓ/dt > 0. Presently, however, it is not clear why such behavior would be physically expected. Thus, the addition of bulk resistivity does not seem to easily solve the braking index puzzle, and a more elaborate explanation is still required.

Our resistive magnetospheres naturally lend themselves to modeling the sub-pulse drift phenomena observed in a number of pulsars (see, e.g., Askegar & Deshpande 2001; Deshpande & Rankin 2001; Rankin et al. 2003; Contopoulos & Spitkovsky 2006). The following discussion takes place in the corotating frame where we define electromagnetic fields by Equations (10) and (11). In this frame there is transverse particle drift with velocity of order ${\bm v}_{\bf drift}\sim \bm{E}_{\perp }\times \bm{B}/B^2$. This drift velocity can be thought of as the minimal velocity of particles, ignoring motion along the magnetic field lines. The transverse electric field E_⊥ reverses sign through the axis of maximum potential drop, and so there is rotation of field lines about this axis. In general, this axis is not in the direction of the magnetic moment, and in fact does not even need to be a straight line. If we consider a simplified picture of the magnetosphere in which there is no parallel electric field in the closed field line region, the potential drop across open field lines is given roughly by the longitudinal potential drop along field lines. The differential rotation can then be roughly related to the longitudinal potential drop by

$$\begin{equation} \Delta \Omega = \Delta V_{\rm drop}/\left(B_p R_\perp ^2\right) = \Delta V_{\rm drop} / \Psi _{\rm cap}. \end{equation} \tag{ 16 }$$

In the laboratory frame, the field lines are undergoing rotation about the pulsar spin axis (with the angular frequency Ω), and in addition they are rotating differentially around the maximum potential axis in a retrograde sense. The differential rotation of plasma can explain sub-pulse emission features that drift with respect to the otherwise periodic light curve of the pulsar. Our result is in contrast to the standard picture in which the plasma precesses about the magnetic axis (Ruderman & Sutherland 1975). The magnetic colatitude parameter in the cartographic transformation of, e.g., Deshpande & Rankin (2001), must then be modified accordingly.

To better visualize the differential rotation in the context of our resistive magnetospheres, consider the full drift velocity

$$\begin{equation} {\bm v}_{\bf drift}=\frac{\bm{E}\times \bm{B}}{B^2+E_0^2}c, \end{equation} \tag{ 17 }$$

again in the corotating frame. Figure 5 shows magnetic field lines in the μ − Ω plane for a 60° inclined dipole at conductivity (σ/Ω)² = 4. Color indicates out-of-plane drift velocity, with red (blue) representing velocity into (out of) the page. The velocity has been rescaled by raising its magnitude to the 1/2 power. The axis of maximum potential extends from the star roughly halfway between the rotation and magnetic axes. The differential rotation of plasma about this axis is denoted by the transition from red to blue across the axis, indicating a reversal in the sign of the out-of-plane drift velocity. Note that the axis bends off in the direction of open field lines beyond the LC. These corotating frame drift velocities can be quite substantial. The maximum illustrated drift velocity occurs near the current sheets and corresponds to a velocity of 0.8c. The differential rotation is much stronger than that implied by sub-pulse drift, but we expect the qualitative features of differential rotation to persist in solutions with more realistic open field line potential drops. At lower conductivities the magnitude of differential rotation about the axis of maximum potential drop increases, and the conductivity is a parameter that can be tuned to attempt to match observed sub-pulse emission features.

Zoom In Zoom Out Reset image size

Incidentally, the bundle of field lines surrounding the axis of the maximum potential corresponds to the field lines that carry the largest distributed current density on the polar cap (see Figure 4 in Bai & Spitkovsky 2010). If the strength of the current is associated with radio emission, then the core emission would have a centroid that is offset from the magnetic pole. The offset could be as large as half of the polar cap radius. This would introduce potentially significant modifications to the polarization sweep of the core radio emission and cause deviations from the expected S curve of Radhakrishnan & Cooke (1969). For a number of pulsars, it would imply differences in the inclination and viewing angles inferred from the shape of the polarization sweep.

We conclude with prospects for future research. We will generalize our assumption that σ/Ω is constant in the magnetosphere. In addition to having higher conductivity in the closed zone, we will experiment with prescriptions for anomalous resistivity in the current sheet. These improvements will yield a more realistic magnetospheric structure, that can be accurate enough for geometrical modeling of gamma-ray light curves from pulsars. The gamma-ray pulse formation is sensitive to the geometry of the magnetosphere, and, in particular, to the field lines near the current sheet (Bai & Spitkovsky 2010). Deviations from the force-free geometry may be important for modeling the light curves of older pulsars which require wider gaps in the outer gap models (Watters et al. 2009).

A.T. acknowledges support by the Princeton Center for Theoretical Science fellowship and by the National Science Foundation through TeraGrid resources provided by NICS Kraken under grant number TG-AST100040. A.S. is supported by NSF grant AST-0807381 and NASA grants NNX09AT95G and NNX10A039G. We thank the referee, Ioannis Contopoulos, for comments that helped improve this paper. The simulations presented in this paper were performed on computational resources supported by the PICSciE-OIT High Performance Computing Center and Visualization Laboratory. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under contract No. DE-AC02-05CH11231.

We derive here in full generality our resistive current prescription. We start in the fluid frame, with charge ρ₀ and current flowing along the electric field with magnitude σE₀. For convenience, we pick the current and electric field to point along the positive z axis. The magnetic field can point along the positive or negative z axis, depending on whether B₀ is positive or negative. If we boost along the z axis with β₁ = (0, 0, β_z), the current four-vector in the boosted frame satisfies

$$\begin{equation} \left[\begin{array}{@{}c@{}} \rho ^{\prime }c \\ \\ \\ j^{\prime }_z \end{array} \right] = \left[\begin{array}{@{}cc@{}} \gamma _z & -\beta _z \gamma _z\\ \\ \\ -\beta _z \gamma _z & \gamma _z \end{array}\right] \times \left[ \begin{array}{@{}c@{}} \rho _0c \\ 0 \\ 0 \\ \sigma E_0 \end{array} \right]. \end{equation} \tag{ A1 }$$

Boosting again in the x-direction transverse to the electric and magnetic fields with β₂ = (β_x, 0, 0),

$$\begin{equation} \left[\begin{array}{@{}c} \rho c \\ j_x \\ \\ j_z \end{array} \right] = {\left[\begin{array}{@{}cc@{}} \gamma _x & \\ -\beta _x \gamma _x & \\ \\ & 1\end{array}\right]} \times \left[ \begin{array}{@{}c@{}}\rho ^{\prime }c \\ 0 \\ 0 \\ j^{\prime }_z\end{array} \right], \end{equation} \tag{ A2 }$$

we obtain for lab frame quantities the system of equations

$$\begin{equation} \rho =\gamma _x \gamma _z \rho _0 -\gamma _x \beta _z \gamma _z \sigma E_0/c, \end{equation} \tag{ A3 }$$

$$\begin{equation} j_x = -\beta _x c\rho, \end{equation} \tag{ A4 }$$

$$\begin{equation} j_z = -\beta _z \gamma _z c\rho _0 + \gamma _z \sigma E_0. \end{equation} \tag{ A5 }$$

Rearranging Equation (A3),

$$\begin{equation} \rho _0=\frac{\rho }{\gamma _x \gamma _z}+\beta _z \sigma E_0/c, \end{equation} \tag{ A6 }$$

and plugging this expression into Equation (A5), yields

$$\begin{eqnarray} j_z &=& {-}\beta _z \gamma _z c\left(\frac{\rho }{\gamma _x \gamma _z}+\beta _z \sigma E_0/c\right) + \gamma _z \sigma E_0 \nonumber \\ &=& {-}\frac{\beta _z}{\gamma _x}c\rho +\frac{1}{\gamma _z}\sigma E_0. \end{eqnarray} \tag{ A7 }$$

We next determine the magnitude of the perpendicular boost in the x direction. The electromagnetic fields are invariant under the boost β₁. After applying the boost β₂, the lab frame electromagnetic fields become

$$\begin{eqnarray} E_x &=&0 \quad\nonumber\\ E_y &=&\gamma _x(E^{\prime }_y-\beta _x B^{\prime }_z)=-\gamma _x\beta _x B_0 \nonumber\\ E_z &=& \gamma _x(E^{\prime }_z+\beta _x B^{\prime }_y)=\gamma _x E_0 \nonumber \\ B_x &=&0 \quad \nonumber\\ B_y &=&\gamma _x(B^{\prime }_y+\beta _x E^{\prime }_z)=\gamma _x\beta _x E_0 \nonumber \\ B_z &=& \gamma _x(B^{\prime }_z+\beta _x E^{\prime }_y)=\gamma _x B_0. \end{eqnarray} \tag{ A8 }$$

Hence,

$$\begin{equation} B^2=\frac{\beta ^2_xE^2_0+B^2_0}{1-\beta ^2_x}, \end{equation} \tag{ A9 }$$

$$\begin{equation} \beta ^2_x=\frac{B^2-B^2_0}{B^2+E^2_0}, \end{equation} \tag{ A10 }$$

and

$$\begin{equation} \gamma ^2_x=\frac{B^2+E^2_0}{B^2_0+E^2_0}. \end{equation} \tag{ A11 }$$

The lab frame x direction corresponds to the −E × B direction. The y and z coordinates are rotated with respect to the B and E_⊥ directions. We let θ denote the angle between the z axis and the laboratory frame magnetic field and assume without loss of generality that β_x > 0. The current in the E × B direction is then

$$\begin{equation} j_{\widehat{\vec{E}\times \vec{B}}}=\beta _x c\rho. \end{equation} \tag{ A12 }$$

Noting that

$$\begin{equation} E^2_{\perp }=\bm{E}^2-(E_0B_0/B)^2=B^2+E^2_0-B^2_0-E^2_0B^2_0/B^2,\quad \end{equation} \tag{ A13 }$$

$$\begin{equation} \beta _x=\sqrt{\frac{B^2-B^2_0}{B^2+E^2_0}}=\frac{E_{\perp }B}{B^2+E^2_0}. \end{equation} \tag{ A14 }$$

The boost −β₁ = (− β_x, 0, 0) is the slowest boost from the lab frame that brings the electric and magnetic fields parallel to one another, and β_x can be thought of as a generalized drift velocity. From β_x, we immediately obtain the transverse current

$$\begin{equation} j_{\widehat{\bm{E}\times \bm{B}}}=\frac{\rho E_{\perp }B}{B^2+E^2_0}c. \end{equation} \tag{ A15 }$$

The current in the direction of the magnetic field is

$$\begin{eqnarray} j_{\widehat{\bm{B}}}&=&j_z {\rm cos} \theta =\left(-\frac{\beta _z}{\gamma _x}\rho c+\frac{1}{\gamma _z}\sigma E_0\right)\frac{\gamma _x B_0}{B} \nonumber \\ &=&\left(-\beta _z\rho c+\frac{\gamma _x}{\gamma _z}\sigma E_0\right)\frac{B_0B+E_0 \bm{E}\cdot \bm{B}/B}{B^2+E^2_0}. \end{eqnarray} \tag{ A16 }$$

The current in the direction E_⊥ is

$$\begin{eqnarray} j_{\widehat{{{\bm E}_{\perp }}}}&=&j_z {\rm sin} \theta =\left(-\frac{\beta _z}{\gamma _x}\rho c+\frac{1}{\gamma _z}\sigma E_0\right)\frac{\gamma _x \beta _x E_0}{B} \nonumber \\ &=&\left(-\beta _z\rho c+\frac{\gamma _x}{\gamma _z}\sigma E_0\right)\frac{E_0E_{\perp }}{B^2+E^2_0}. \end{eqnarray} \tag{ A17 }$$

Hence, the laboratory frame current vector can be expressed in full generality as

$$\begin{eqnarray} &&\bm{j}{=}\frac{\rho c\bm{E}\times \bm{B}+\left(\!{-}\beta _z\rho c+\sqrt{\frac{B^2+E^2_0}{B^2_0+E^2_0}(1{-}\beta _z^2)}\sigma E_0\right)(B_0\bm{B}+E_0\bm{E})}{B^2+E^2_0}.\nonumber\\ \end{eqnarray} \tag{ A18 }$$

In Section 2, β_z has been replaced by β_|| to emphasize that it is the magnitude of the boost along the parallel electric and magnetic fields to the fluid rest frame.

In principle we should be measuring spin-down values at the stellar surface, but there is uncertainty in the stellar spin-down measurement due to stair-stepping at the spherical inner boundary on our Cartesian grid. To avoid this issue, we measure spin-down at the LC. Dissipation of Poynting flux inside the LC artificially suppresses our spin-down estimates measured at the LC, though. We quantify these uncertainties by computing the force-free spin-down values measured both on the star and at the LC. The true force-free spin-down luminosity likely falls within the shaded gray region in Figure 2 bounded by the stellar and LC spin-down values. There is also a deviation from the vacuum Deutsch solution for our zero conductivity solution (compare dashed line to (σ/Ω)² = 0 line). This difference is due to our boundary conditions at the star. The smoothing of the fields across the stair-stepped boundary leads to small charge bleed-off from the inside of the star to the outside. The first term in Equation (3) is then not identically zero just outside the star even in the σ = 0 case. This numerical artifact does not affect the solutions with high conductivity when the physical current exceeds the numerical smoothing current, but below (σ/Ω)² = 4 × 10⁻³ it can influence the field structure and slightly modify the spin-down power.

To better understand these numerical issues it is helpful to look at the run of Poynting flux with radius. In Figure 6, we show the Poynting flux integrated over spherical shells of varying radii for force free, vacuum, and a range of resistive solutions with α = 60°. The results are normalized to the force-free value at the star, located at R_* = 0.375R_LC, and we show the run of Poynting flux with spherical radius out to 2.25R_LC. The vacuum and (σ/Ω)² = 0 curves are flat, indicating negligible dissipation of Poynting flux with increasing radius, but there is an offset between the curves that we attribute to our imperfect inner boundary. Force free is in principle dissipationless and should have the run of Poynting flux flat with radius. As was argued above, our method for cleaning parallel electric field in force-free simulations is resistive and causes Poynting flux to slope downward with increasing radius. Inside the LC, the drop is due to artificial volume $\bm{j}\,{\bm \cdot}\, \bm{E}$ dissipation above the polar caps. This dissipation varies with magnetic inclination angle and is responsible for the varying width of the gray error band with angle for force-free solutions in Figure 2. Outside the LC, the drop in Poynting flux is due to it disappearing into the current sheets.

Zoom In Zoom Out Reset image size

The resistive runs in Figure 6 have Poynting flux dropping with increasing radius, as we expect, indicating dissipation throughout the magnetosphere. Inside the LC the dissipation is primarily physical $\bm{j}\,{\bm \cdot}\, \bm{E}$ dissipation above the polar caps. There is also some dissipation due to the disappearance of Poynting flux into the current layer, though significantly weaker in magnitude as compared to the volume $\bm{j}\,{\bm \cdot}\, \bm{E}$ dissipation. Outside the LC the Poynting flux is dissipated in the current sheets. The dissipation of Poynting flux is strongest at intermediate conductivities, σ/Ω ∼ 1. At higher conductivities, the currents are more ideal and dissipationless, and at lower conductivities, dissipation is reduced because the conduction currents are weaker. Volume dissipation introduced in our resistive simulations implies that the LC measurements of Poynting flux in Figure 2 underestimate the true spin-down power of the resistive solutions, and should be treated as being within the error bands similar to the gray band of the force-free solution.