AB INITIO PULSAR MAGNETOSPHERE: THREE-DIMENSIONAL PARTICLE-IN-CELL SIMULATIONS OF OBLIQUE PULSARS

In the last 15 years the pulsar magnetosphere was intensively studied in the magnetohydrodynamic (MHD) limit, including ideal force-free electrodynamics (Contopoulos et al. 1999; Gruzinov 2005; McKinney 2006; Spitkovsky 2006; Timokhin 2006; Kalapotharakos & Contopoulos 2009; Pétri 2012), resistive force-free (Kalapotharakos et al. 2012; Li et al. 2012), and full relativistic MHD (Komissarov 2006; Tchekhovskoy et al. 2013). While these solutions agree on the general shape of the plasma-filled magnetosphere, the origin of the magnetospheric plasma and the properties of particle acceleration necessary for production of the observed radiation cannot be addressed within the MHD approach, necessitating a kinetic treatment.

Recently, the aligned pulsar magnetosphere was studied with first-principles relativistic particle-in-cell (PIC) simulation (Cerutti et al. 2014; Chen & Beloborodov 2014, CB14 hereafter; Philippov & Spitkovsky 2014, PS14 hereafter). PS14 showed that injecting abundant pairs everywhere in the magnetosphere drives the solution toward the force-free configuration. Cerutti et al. (2014) studied the transition between the disk–dome configuration and the filled solution by varying the rate of particle injection at the stellar surface. Beloborodov (2008) and Timokhin & Arons (2013) studied one-dimensional polar cap cascade and found that in the region where the current density along the magnetic field, ${{j}_{\parallel }}$, is below the Goldreich–Julian value, ${{j}_{{\rm GJ}}}=-{{{\bf \Omega }}_{*}}\cdot {\boldsymbol{B}} /2\pi$ (here, ${{{\bf \Omega }}_{*}}$ is the stellar angular velocity vector, and ${\boldsymbol{B}}$ is the local magnetic field), particles in the space-charge-limited flow are not accelerated up to relativistic energies and do not produce secondary pairs. Efficient particle acceleration, and, consequently, active pair formation, should thus be confined to the regions where ${{j}_{\parallel }}/{{j}_{{\rm GJ}}}\lt 0$ or ${{j}_{\parallel }}/{{j}_{{\rm GJ}}}\gt 1$. CB14 presented the first global aligned simulation with e^± discharges and concluded that pair formation in the outer magnetosphere is required to produce the solution similar to force-free but with strong charge separation, consistent with sub-GJ current flow in the polar cap region. This raised the question of whether the high-multiplicity state of the magnetosphere is ever achieved, particularly in more realistic oblique rotators.

The goal of this Letter, therefore, is to describe the magnetospheres of oblique pulsars. For this, we construct relativistic PIC simulations of the oblique pulsar with different prescriptions for populating the magnetosphere with plasma, focusing on solutions with realistic prescription for pair formation. The Letter is organized as follows. In Section 2 we describe our numerical method, in Section 3 we demonstrate the test solution with abundant neutral plasma supply, and in Section 4 we present our results on three-dimensional magnetosphere structure of aligned and oblique rotators with pair production.

To simulate the magnetosphere we use the 3D electromagnetic PIC code TRISTAN-MP (Spitkovsky et al. 2005). The neutron star can be described as a rotating conducting magnetized sphere, which generates electric fields corresponding to a quadrupolar surface charge. In the limit of zero work function, the accelerating components of electric fields extract charges from the surface and populate the magnetosphere with plasma. Below we ignore the complications of extracting ions from the surface and assume that extracted positive charges are positrons. In addition to the extraction of surface charges, we consider two scenarios for plasma supply. First, we extend the study of PS14 to oblique rotators by considering abundant neutral plasma injection in the whole magnetosphere. Second, we add a more realistic prescription, which approximates pair production by magnetic conversion of photons and two-photon collisions (see Section 4.1).

We use the linearity of Maxwell's equations to represent the electric field as a superposition: ${\boldsymbol{E}} ={{{\boldsymbol{E}} }_{{\bf vacuum}}}+{{{\boldsymbol{E}} }_{{\bf plasma}}}$, where ${{{\boldsymbol{E}} }_{{\bf vacuum}}}$ is the analytical field of the vacuum rotator (Michel & Li 1999), and ${{{\boldsymbol{E}} }_{{\bf plasma}}}$ is the field computed from PIC particles (Spitkovsky et al. 2002). The same decomposition is done for the magnetic field. Inside the conducting sphere, the plasma fields are forced to zero with a smoothing kernel, while the vacuum fields are set to corotation values.

We use a Cartesian grid of 1024³ cells, with stellar radius R_* = 35 cells, and light cylinder R_LC = 120 cells. Particles penetrating into the star by two cells are removed. The outer walls have radiation boundary conditions for fields and particles.

In this section we extend PS14 and build a force-free solution for the oblique rotator. Neutral pair plasma is injected with zero velocity at a fixed rate in every cell within a sphere of R = 2R_LC, if the magnetization in a cell exceeds a certain limit, $\sigma ={{B}^{2}}/(4\pi$ $n{{m}_{e}}{{c}^{2}})\gt 1000{{({{R}_{*}}/r)}^{3}}$. Here, B is the absolute value of the magnetic field, and n is the plasma density. Due to the limiter, the closed zone is not overloaded with plasma, and the flow is well-magnetized. The magnetization of our solution is σ ≈ 1000 near the pole (around 20 particles per cell; local skin depth is 1.1 cells), and σ ≈ 20 at the light cylinder (one particle per cell, skin depth is five cells).

The magnetic field in the plane defined by the magnetic moment ${\boldsymbol{ \mu }}$ and stellar angular velocity ${{{\bf \Omega }}_{*}}$ vectors for a 60° rotator is shown in Figure 1(c), where color represents the out-of-plane component of the magnetic field. In agreement with the MHD solution (Tchekhovskoy et al. 2013), poloidal field lines in this plane become purely radial beyond R_LC. The Y-point is located approximately at the light cylinder, and the current sheet oscillates around the equatorial plane. As shown in Figures 1(a) and (b), our kinetic solution has high multiplicity (≈10) in the bulk of the polar outflow, which carries sub-GJ current (see Figure 1(d)). The regions with current above the GJ value, favorable for efficient pair formation, are the current sheet and the rims of the polar cap, where volume return current flows.

Zoom In Zoom Out Reset image size

The field structure is quasi-stable in the corotating frame. We note that within the simulation box the role of drift-kink instability of the current sheet, pointed out for the aligned rotator by PS14 and studied by Cerutti et al. (2014), decreases with increasing obliquity, becoming negligible already for 60°. This is a consequence of the decrease of conduction current in the current sheet with increasing obliquity.

We ran simulations for obliquities 30°, 60°, and 90°, measuring the Poynting flux integrated over a sphere with r = R_LC. Combining with the aligned case (PS14), we conclude that the spin-down energy loss of the PIC solution is consistent with MHD studies:

$$\begin{eqnarray}&&L=\frac{{{\mu }^{2}}{\Omega }_{*}^{4}}{{{c}^{3}}}({{k}_{0}}+{{k}_{1}}{{{\rm sin} }^{2}}\chi ),\end{eqnarray} \tag{ 1 }$$

with coefficients k₀ = 1.0 ± 0.1, k₁ = 1.1 ± 0.1 (see Figure 5), and χ the inclination angle.

4.1. Pair Creation Process

4.2. Aligned Rotator

The magnetosphere of the aligned rotator with self-consistent pair formation is shown in Figure 2(c). Its magnetic field structure is close to the force-free solution: the field becomes mostly toroidal beyond R_LC, and poloidal field lines stretch and become radial. The opposite polarities of the toroidal field are supported by the equatorial current sheet. We observe moderate kinking of the sheet beyond 2R_LC, which does not affect the Y-point region during our simulation.

Zoom In Zoom Out Reset image size

The polar region is inactive to pair formation, as pointed out by Timokhin & Arons (2013) and CB14. This happens because the current flow is sub-GJ (see Figure 2(d)), thus ${{E}_{\parallel }}$ within the light cylinder is easily screened by the sub-relativistic electron outflow shown in Figure 2(a). However, pair production does take place in the return current layer and in the current sheet, which have high multiplicity of secondary pairs as shown in Figures 2(a) and (b). The acceleration of parent particles occurs due to strong ${{E}_{\parallel }}$ in the current sheet. The electric field points away from the star and positrons are accelerated outwards, while high-energy electrons stream toward the star. These electrons pass through the Y-point and penetrate into the inner magnetosphere by a few Larmor radii. The current sheet width scales with magnetic field strength at the light cylinder, staying of the order of the average Larmor radius of particles. Most of the open field lines do not cross the null surface, and the electron outflow on these lines extends to infinity. The last open field lines, however, go through the null surface. Since the electron outflow cannot cross the null surface, a vacuum gap opens above the current sheet (cf. CB14). The spin-down of this solution, measured as the flux of electromagnetic energy at the stellar surface, is $\approx 0.85{{\mu }^{2}}{\Omega }_{*}^{4}/{{c}^{3}}$.

We note that confining pair creation to the inner magnetosphere, $r\leqslant 2{{R}_{*}}$, or increasing the global pair creation threshold to values comparable with $\gamma \approx {{\sigma }_{{\rm LC}}}$, where σ_LC is the magnetization parameter at the light cylinder, leads to the suppression of return current, and the magnetosphere relaxes to the low multiplicity solution (Cerutti et al. 2014; CB14). In this state the magnetosphere consists of an equatorial disk of positrons and a polar dome of electrons. Positrons are pulled from the tip of the disk by strong electric field, E > B; however, the resulting current is not large enough to substantially modify the dipolar magnetic field structure.

4.3. Oblique Rotator

We found that magnetospheres of rotators with inclination angles χ ≲ 40° are qualitatively similar to the aligned solution, showing no pair creation activity in the polar cap outflow. The bulk outflow in these solutions is charge-separated, and pair production happens only in the current sheet and the return layer. These solutions relax to the oblique disk–dome structure if pair production in the outer magnetosphere is suppressed. This happens because the return current layer, where pair production is possible, cannot be formed by particles lifted from the surface, since the sign of the necessary current is opposite to the sign of available charges. As there is no current outflow, the oblique disk–dome solutions spin down at the vacuum rate.

Solutions with high obliquities, χ ≳ 40°, do show pair creation in the polar cap zone. There are regions on the polar cap with volume return current exceeding the GJ value (Timokhin & Arons 2013). Slices through the 3D magnetosphere with 60° inclination are shown in Figures 3(c) and (e); the 3D structure of the field lines is shown in Figure 3(e). The magnetospheric structure is similar to what we found in Section 3, with closed and open zones, Y-point, and the current sheet (see Figure 3(f) for the 3D current distribution). Pairs are produced in the region of volume return current and in the current sheet. The flow in the sheet is supported not only by the local pair creation, as in the case of the aligned rotator, but also by the reconnection-induced inflow of charges produced in the polar cap discharge. The plasma density is not uniform in the polar cap as seen from Figures 3(a) and (b).

Download figure:

Video Standard image High-resolution image

Part of the polar cap supports charge-separated electron outflow with $\alpha ={{j}_{\parallel }}/{{j}_{{\rm GJ}}}\approx 0.95\lt 1$. The Lorentz factor of this flow is ≈25, independent of the surface magnetic field. This is consistent with the estimate $\gamma \approx 2\alpha /(1-{{\alpha }^{2}})$ for the space-charge-limited flow (Beloborodov 2008; Timokhin & Arons 2013) in the case of sub-GJ current. The volume return current region shows active pair production. The polar discharge has quasi-periodic behavior on timescales that are short compared to the pulsar rotation period, with episodes of efficient pair production followed by quiet states. The three-dimensional study of the cascade intermittency is important for modeling the pulsar radio emission and will be discussed elsewhere. The multiplicity of the flow in the active region of the polar cap is ≈10, close to the limiting value. We performed a simulation with a twice larger limiting value and observed that the multiplicity of the flow increases by roughly the same factor. Though the whole polar cap does not show active pair formation, on average the pulsar wind is supplied with plasma of high multiplicity, concentrated in the equatorial current sheet wedge. The current sheet beyond the light cylinder should thus be the main source of γ-photons. As in the aligned case, mainly positrons are accelerated in the current sheet of the 60° rotator. We find that current sheet of the pair-forming solution shows more efficient kinking than the test solution from Section 4, though it does not strongly affect the global magnetospheric structure.

We find that the volume return current region in our simulations is bounded by the null surface that crosses the polar cap; this current is carried by positrons. However, in the limit of very small star ${{R}_{*}}/{{R}_{{\rm LC}}}\to 0$ the null surface is located outside the polar cap for inclinations that are not too close to 90°, and the outflowing current is supposed to be carried only by electrons. Thus, the pair-producing region with j > j_GJ may disappear.¹ We conclude that the precise obliquity angle for transition between active and inactive states of the polar cap may change for ${{R}_{*}}\ll {{R}_{{\rm LC}}}$. Although our present simulations are already physically relevant for millisecond pulsars with high ${{R}_{*}}/{{R}_{{\rm LC}}}$, we will investigate the smaller values of ${{R}_{*}}/{{R}_{{\rm LC}}}$ in future work.

The 90° rotator shows efficient pair production in the whole polar cap (Figure 4). In contrast to the aligned rotator, the current flow on the polar cap is anti-symmetric with respect to the equatorial plane. Current density exceeds j_GJ, necessitating the presence of e^± discharge. Pair formation also happens in the current sheet.

Zoom In Zoom Out Reset image size

In Figure 5 we present the spin-down of our oblique solutions with pair formation, measured at the stellar surface. In our solution for the aligned rotator with pair discharges, approximately 20% of the Poynting flux is dissipated within 2R_LC. The dissipated fraction of total energy losses decreases for larger inclination angles and reaches 3% for the orthogonal rotator. The dissipated energy is converted into particle energy in acceleration regions. The deviation of spin-down from solutions with abundant pair formation discussed in Section 3 is small and becomes negligible with increasing inclination. The origin of this deviation is due to the presence of vacuum-like regions inside the light cylinder that are not corotating with the star.

Zoom In Zoom Out Reset image size

In this Letter we presented the magnetospheric structure of oblique pulsars with pair formation. We found that pulsars with low obliquities produce active solutions only if pair formation in the outer magnetosphere is working. In these solutions the pairs are produced in the equatorial current sheet beyond the light cylinder and in the return current layer. If the radio emission is produced in the high-multiplicity flow at the polar cap, these pulsars should have mainly double-peaked profiles associated with return current layers. Pulsars with obliquity angles ≳40° show significant pair production activity in the polar regions. Their radio beams may a have multi-peaked structure. The spin-down of pair solutions is close to MHD model predictions.

Our prescription for pair formation is highly idealized. For simplicity we neglected the mean free path of photons and assumed the pair formation threshold to be independent of the local conditions. While this approximation is reliable for the polar cap discharge, the γ–γ process which can account for the activity in the current sheet is intrinsically non-local, which requires more careful treatment. Two processes may contribute to launching the e^± discharge in the current sheet: collisions of energetic synchrotron photons emitted by counter-streaming particle flows, and collisions of thermal X-ray photons from the neutron star surface, which may be up-scattered by energetic particles in the sheet, with synchrotron photons. Both mechanisms should operate most efficiently close to the Y-point, where magnetic field in the sheet is strongest, and the emitted synchrotron photons are the most energetic. We assumed efficient pair formation in the current sheet, which may not happen in pulsars with low magnetic field strength at the light cylinder. In low obliquity pulsars with no pair formation in the current sheet, the return current layer is suppressed, and such pulsars should not shine in the radio band, unless some physics beyond our model can reestablish pair formation.

A potentially important complication may arise due to the restricted ability of ions to produce pairs. Though ions can still emit curvature photons capable of producing pairs in γ–γ collisions with stellar photons, the resulting multiplicity of the flow due to this process is uncertain. This may affect the multiplicity of the polar outflow in the orthogonal rotator, where ions are extracted in the half of the polar cap.

The poloidal field at the stellar surface may be more complicated than the pure dipole field considered in this work due to the operation of Hall effect in the crust (e.g., Gourgouliatos & Cumming 2014). The resulting multipolar field can modify the operation of the polar discharge, producing more regions with efficient pair production.

The value of $\alpha ={{j}_{\parallel }}/{{j}_{{\rm GJ}}}$ that is essential for operation of polar discharge may be sensitive to general relativistic corrections. In particular, the value of GJ charge density is reduced due to the Lense–Thirring effect (Beskin 1990; Muslimov & Tsygan 1992). The increased value of α may lead to the ignition of polar discharge in low obliquity pulsars. This will be investigated in future work.

We thank Jonathan Arons, Vasily Beskin, Alexander Tchekhovskoy, and Andrey Timokhin for fruitful discussions. This research was supported by NASA grants NNX12AD01G and NNX13AO80G, Simons Foundation (grant 267233 to A.S.), and was facilitated by Max Planck/Princeton Center for Plasma Physics. The simulations presented in this article used computational resources supported by the PICSciE-OIT High Performance Computing Center and Visualization Laboratory.