Radio Emission and Electric Gaps in Pulsar Magnetospheres

The origin of pulsar radio emission is one of the old puzzles in theoretical astrophysics. In this Letter, we present a global kinetic plasma simulation that shows from first principles how and where radio emission can be produced in pulsar magnetospheres. We observe the self-consistent formation of electric gaps that periodically ignite electron-positron discharge. The gaps form above the polar cap and in the bulk return current. Discharge of the gaps excites electromagnetic modes, which share several features with the radio emission of real pulsars. We also observe the excitation of plasma waves and charge bunches by beam instabilities in the outer magnetosphere. Our numerical experiment demonstrates that global kinetic models can provide deep insight into the emission physics of pulsars and may help interpret their multiwavelength observations.


INTRODUCTION
Pulsars are rotating magnetized neutron stars which produce powerful beams of coherent radio emission.Despite an abundance of observational data, the mechanism generating radio waves in pulsar magnetospheres has remained elusive for more than fifty years.Pulsar magnetospheres are filled with highly magnetized collisionless electron-positron (e ± ) plasma which is produced in electric gaps -regions with voltage along magnetic field lines.Electric discharge in the gaps may be responsible for pulsar radio emission, and the only reliable way to solve this nonlinear problem is with a self-consistent numerical simulation.
In the last decade there has been a significant computational effort to model the magnetosphere structure and multi-wavelength emission from first-principles using the particle-in-cell (PIC) technique (Chen & Beloborodov 2014;Philippov et al. 2015;Philippov & Spitkovsky 2018).Global simulations unanimously display powerful gamma-ray emission from the outer magnetosphere (Gruzinov 2013;Chen & Beloborodov 2014; Corresponding author: Ashley Bransgrove abransgrove@princeton.eduPhilippov & Spitkovsky 2018), but radio waves were not observed.Recently Philippov et al. (2020) showed that radio waves can be generated directly by the e ± discharge above the polar-cap.The authors used local PIC simulations of the gap in Cartesian geometryit required high voltage and resolution which were not achieved in global simulations.
In this Letter we present a global kinetic plasma simulation of an axisymmetric magnetosphere.Our numerical experiment achieves voltage and spatial resolution sufficient to reveal new features, including locations of electric gaps where e ± discharge ignites and produces coherent radio waves.We observe beam instabilities in the magnetosphere for the first time.

METHOD
The kinetic plasma simulation is performed with the relativistic PIC code Pigeon (Hu et al. 2021;Hu & Beloborodov 2022), a descendant of Aperture (Chen & Beloborodov 2014), which solves Maxwell's equations with the equations of motion for charged particles and e ± creation.We include a general-relativistic correction to Faraday's law which describes the rotational framedragging effect near the star (Philippov & Spitkovsky 2018).We set the compactness of the star r s /r ⋆ = 0.5, where r s = 2GM/c 2 is the Schwarzschild radius and r ⋆ is the radius of the star.Hereafter lengths are given in units of r ⋆ , and times in units of r ⋆ /c.The e ± have mass m and charge ±e, while ions (protons) have charge e and mass m i = 10m.Electromagnetic fields are given in units of mc 2 e −1 r −1 ⋆ , number density in units of r −3 ⋆ , and charge density in units of er −3 ⋆ .The scales of the problem are set by the angular velocity Ω and the magnetic dipole moment µ of the pulsar.Rotation of the magnetized star induces a voltage Φ 0 ≈ µΩ 2 /c 2 capable of accelerating e ± to Lorentz factor γ 0 = eΦ 0 /mc 2 .For real pulsars γ 0 ∼ 10 10 .We scale it down to γ 0 = 2 × 10 4 , which allows us to resolve the plasma skin depth.Our value of γ 0 is twice as large as previous work (Hu & Beloborodov 2022), and our resolution is twice as high.We set Ω = 1/6, which implies the light-cylinder radius R LC = c/Ω = 6.The minimum charge density required to support the co-rotating magnetosphere ρ GJ = −Ω Ω Ω • B B B/(2πc) (Goldreich & Julian 1969) defines a characteristic number density n GJ = |ρ GJ |/e and a plasma frequency ω p = (4πn GJ e 2 /m) 1/2 .Similar to real pulsars, the simulation preserves the hierarchy of scales Ω ≪ ω p ≪ ω B , where ω B = eB/mc is the cyclotron frequency.
The gamma-ray emission, propagation, and e ± production is modelled using the Monte-Carlo method (Chen & Beloborodov 2014).The e ± emit gamma-rays of energy ϵ γ = 10mc 2 when they reach the threshold energy γ thr mc 2 .The threshold depends on the curvature of the field lines (Chen & Beloborodov 2014), and has typical value γ thr ∼ 100.The gamma-rays propagate with mean free path l before converting to secondary e ± with γ s mc 2 = ϵ γ /2.For r < 2, we set l to simulate conversion off the strong magnetic field near the stellar surface (γ − B channel).For r > 2, we set l to simulate collisions with soft target photons in the outer magnetosphere (γ − γ channel).The energy scales in our numerical experiment satisfy the correct hierarchy 1 ≪ γ s ≪ γ thr ≪ γ 0 .
The computational domain covers 1 ≤ r ≤ 30 and 0 ≤ θ ≤ π.The grid is uniform in log r and θ, and has resolution N r ×N θ = 8192×8192 cells.The plasma scale 2πc/ω p near the stellar surface is resolved by 25 grid cells.At the surface, electromagnetic fields satisfy a rotating conductor boundary condition, and we maintain a gravitationally bound electron-ion atmosphere with multiplicity M atm = n atm /n GJ = 10.At the outer boundary fields are damped, and particles are absorbed.

MAGNETOSPHERE STRUCTURE
The simulation begins with a non-rotating star with a vacuum dipole magnetosphere.We smoothly increase the angular velocity from zero to Ω during the first 5 r ⋆ /c, and observe that charges are lifted from the atmosphere and accelerated, triggering gamma-ray emission and e ± production which fills the magnetosphere with plasma.The magnetosphere reaches a quasi-steady state after the first revolution of the star (Fig. 1).Its global structure is similar to previous axisymmetric force-free and PIC models (Contopoulos et al. 1999;Gruzinov 2005;Parfrey et al. 2012;Chen & Beloborodov 2014).The higher voltage in our simulation leads to more efficient filling of the open field lines with plasma.A new feature is the e ± discharge in the outer magnetosphere outside the Y-shaped current sheet.This feature was not seen in previous PIC simulations with lower voltages.We also observe e ± discharge near the polar caps, an expected result of the relativistic frame-dragging effect (Muslimov & Tsygan 1992).The strongest dissipation and gamma-ray emission occur in the equatorial current sheet.

GAPS AND ELECTRIC DISCHARGE
The gaps (regions with non-zero E || = E E E •B B B/|B B B|) form because of a mismatch between the required parallel current j B = (c/4π)(∇ × B B B) • B B B/|B B B| and the maximum current j max = cρ GJ that could be supplied by the outward flow of a charge-separated plasma with the local charge density ρ = ∇ • E/4π.If the charged plasma under-supplies the required current (α ≡ j B /j max > 1), or supplies it with the wrong sign of charge (α < 0), Ampere's law guarantees the inductive growth of E || which can trigger the runaway production of e ± (Beloborodov 2008).The time-dependent gaps observed in the simulation differ from the gap models based on electrostatic considerations (Gauss' law), including the classical polar-gap (Ruderman & Sutherland 1975), slot-gap (Arons 1983;Muslimov & Harding 2004), and outer-gap (Cheng et al. 1986).
A time-dependent gap with α > 1 forms near the magnetic axis above the pulsar polar cap.This gap occurs because the general-relativistic frame-dragging effect reduces the apparent rotation of the star, and thus reduces ρ GJ near the surface, while leaving the required cur-rent j B unchanged (Muslimov & Tsygan 1992;Philippov et al. 2015;Belyaev & Parfrey 2016;Gralla et al. 2016).The resulting α > 1 renders the electron flow extracted from the star unable to supply j B , even if it moves at the speed of light.Thus, E || is induced, igniting e ± discharge.The discharge begins when electrons are lifted from the atmosphere by E || and accelerated to energies ∼ γ thr mc 2 .The electrons emit gamma-rays, which convert to e ± through the γ − B channel.The increasing density of pairs quickly screens E || , and the e ± shower proceeds outward.As the cloud of e ± leaves, E || grows again and the discharge cycle repeats.The discharge reaches maximum multiplicity M ≡ n/n GJ ∼ 10.The gap height is h ∼ 0.4r ⋆ , and the cycle time is ∼ h/c.
The magnetosphere prevents the accumulation of net electric charge on the neutron star by arranging positive return-currents, which compensate the negative charge flowing out of the polar-cap (Appendix A).Most of the return-current flows along the thin separatrix layer and is sustained by pair creation at the Y-point (Hu & Beloborodov 2022).We also observe a bulk return-current of macroscopic thickness on open field lines outside the separatrix layer which pass through the null-surface defined by ρ GJ = 0 [Fig. 2 top panel].The bulk return-current draws electrons from the pair plasma created around the light-cylinder.
The simulation reveals two time-dependent gaps on the field lines that conduct the bulk return-current: (i) an inner return-current gap of type α < 0 forms above the polar-cap where ρ GJ < 0 and j B > 0, and (ii) an outer return-current gap of type α > 1 forms outside the null-surface where ρ GJ > 0 and j B > 0. The presence of two gaps on the same field lines results in coupled discharge dynamics.The two interacting discharges reach e ± multiplicities M ∼ 5 and repeat on the timescale ∼ 2R LC /c.
The inner return-current gap has ρ GJ < 0, so here the charge density can only be supplied by electrons.Electrons flow toward the star in order to conduct j B > 0. When this flow dwindles, a gap opens inside the nullsurface (with E || directed away from the star) and expands inward until it reaches altitude r ∼ 2r ⋆ where pair production through the γ − B channel becomes possible.Then discharge develops, seeded by the inward electron flow accelerated in the gap.Pair creation quickly screens E ∥ , and the shower proceeds toward the star, spreading laterally due to the curvature of the magnetic field lines and the finite free paths of tangentially emitted gamma-rays.Some of the secondary positrons are reversed by E ∥ at the discharge onset; they escape outward and serve as seeds for the discharge in the outer return-current gap.
The outer return-current gap extends from the nullsurface toward R LC (Fig. 2, top).The gap occurs because the required current is positive, j B > 0, but positrons are not readily available from inside the nullsurface where ρ GJ < 0. This generates E ∥ directed away from the star.It pulls electrons into the gap from larger radii.The inward flowing electrons supply the correct sign of the current, but the wrong sign of charge (ρ GJ > 0), so they cannot screen the gap.Instead, they get accelerated by E ∥ , cross the null-surface and trigger the inner return-current discharge, which provides a source of positrons from inside the null-surface.These seed positrons get accelerated and trigger pair production in the outer return-current gap through the γ − γ channel.Thus, the two gaps (inner and outer) assist each other in repeating pair discharge.The waves are remarkably consistent with a mechanism previously seen in local-box simulations (Philippov et al. 2020;Cruz et al. 2021): pair production in the gap drives a variable current which couples to oscillations of δE E E and δB B B. Two properties of the observed δE E E show that it is the superluminal O-mode (Arons & Barnard 1986): 2], and δE ∥ /δE ⊥ is consistent with the superluminal O-mode polarization (Appendix B).Excitation occurs with frequencies ω comparable to the plasma frequency ω p = (4πne 2 /m) 1/2 , which is rapidly growing during the discharge (Tolman et al. 2022).The growing ω ∼ ω p and the moderate k (large length-scale of the discharge) lead to the excitation of the superluminal branch ω/k > c rather than subluminal Alfvén waves.The waves reach amplitude E ⋆ || ∼ γ s mc ω ⋆ /e, where ω ⋆ is the plasma frequency later in the discharge when screening reduces E || so that it can barely reverse e ± (Tolman et al. 2022).The emitted waves track adiabatically along the O-mode branch as they propagate through the decreasing plasma density, and should escape the magnetosphere as vacuum waves with k ⋆ = ω ⋆ /c (the escape is not captured in our simulation, because at large radii the grid cell exceeds k −1 ⋆ .)The wave emission is strongly modulated on the gap light-crossing timescale.
Our simulation also demonstrates the extraordinary (X) modes emitted from the equatorial plasmoid chain [Fig.2, insets (iii), (iv)] (Philippov et al. 2019;Lyubarsky 2019).These waves have δE E E ∥ (k k k × B B B) and show up only in δE E E ⊥ with no contribution to δE ∥ .Note that in axisymmetry k k k lies in the poloidal plane.Near the star the background magnetic field B B B is approximately poloidal.At large radii B B B is twisted out of the poloidal plane.

PLASMA INSTABILITIES
The simulation also reveals beam plasma instabilities, triggered at several sites in the outer magnetosphere.The triggering mechanism is non-local: the e ± wind becomes exposed to a gamma-ray beam emitted elsewhere in the magnetosphere, which injects an e ± beam with a different momentum (Fig. 3).Both streams are made of secondary e ± , and thus the instability growth rate is not strongly suppressed by ultra-high Lorentz factors typical for primary accelerated particles.The instability is electrostatic in nature, and controlled by plasma os- cillations along the magnetic field lines (rapid gyration of the particles increases their effective parallel inertia, but does not change the qualitative behaviour of the instability).Charge bunches form at each site where the instability occurs (Fig. 3).
The beam energy is deposited into growing plasma modes seen in Fig. 2, inset (iv).The modes grow in the outflowing wind over a characteristic distance d consistent with the growth rate of the relativistic beam instability d ∼ c/Γ (Appendix C) (Egorenkov et al. 1983).The measured wavelength of the modes is consistent with λ 0 ∼ 2πc/(ω p γ 1/2 ), indicating that the instability is first triggered on Cherenkov resonance (Egorenkov et al. 1983).The modes are inclined to the local magnetic field, and have transverse and longitudinal components [Fig.2, inset (iv)].The phase speed is necessarily subluminal.Therefore, we identify the excited waves as hybrid electrostatic-Alfvén modes (Arons & Barnard 1986;Rafat et al. 2019).The waves have amplitude E || ∼ 4π∆ρλ 0 , where ∆ρ ∼ ρ GJ is the charge density of the observed bunches (Fig. 3).The fate of the waves and radiation from the charge bunches should be investigated with detailed local simulations.

DISCUSSION
Our numerical experiment has achieved a voltage high enough to reveal the locations where the e ± discharge ignites and produces coherent radio waves.In addition to the well known polar-cap gap, the simulation has uncovered a pair of interacting gaps in the outer magnetosphere, which form around the null-surface in the bulk return-current.Furthermore, we observed waves and charge bunches excited by beam instabilities in the secondary e ± plasma -a result of multiple gamma-ray streams revealed by the global simulation.The high resolution is key to these findings, providing a glimpse of the origin of pulsar radio waves, a subtle phenomenon compared to pulsar gamma-ray emission.
The global kinetic simulation offers a physical framework to start interpreting the rich pulsar observations.In particular, the polar-cap discharge may be responsible for the so-called "core" radio emission, and the inner return-current discharge may produce "conal" radio emission.Their strong modulation by the electric discharge of the gaps may produce microstructure in the pulse profiles.Also interesting are the additional plasma instabilities found in the outer magnetosphere, which are special to energetic pulsars with active pair creation around the light cylinder.
However, significant work is needed to improve the simulations before they can be directly compared to individual pulsars.In particular, increasing multiplicity M may change the interaction of the two return-current gaps: instead of assisting each other with seed particles for discharge, the inner gap with a high M may flood the outer gap and switch off its radio emission.A more detailed implementation of gamma-ray emission and e ± production may change the beam instabilities and charge bunching.It has been suggested that charge bunches could produce coherent curvature emission in the radio band (Ruderman & Sutherland 1975;Gil et al. 2004) and the waves experience propagation effects (Barnard & Arons 1986;Lyubarsky 1993;Lyubarskii & Petrova 1998;Beskin & Philippov 2012), which needs further investigation.Finally, our simulation was limited to aligned rotators µ ∥ Ω.It will be essential to investigate radio emission in inclined rotators.Figure 4 summarizes the magnetospheric structure and the discharges that fill it with plasma, as observed in our simulation.The location of the discharges is controlled by the global magnetospheric currents, and occurs where α ≡ j B /(cρ GJ ) > 1 or α < 0. The discharge in the out-flowing current above the polar cap (Fig. 4, white region) is well known; it develops because here α is increased above unity by the frame dragging effect in the spacetime of the rotating star (Muslimov & Tsygan 1992).Discharge in the return-current zone is more subtle.Approximately 90% of the return current flows in the separatrix current sheet, and 10% flows in the grey region (Fig. 4), filling a significant volume with α < 0 near the star.In this "bulk" return-current region above the polar cap, we observe strong episodic e ± discharges.Similar discharges with α < 0 were previously simulated in one-dimensional (plane-parallel) models (Timokhin & Arons 2013).The discharge observed in our simulation, however, operates in the 'thin-tube' regime (the transverse size of the electric gap is smaller than its length along B B B) and does not admit a plane-parallel description.The discharge is seeded by the time-dependent plasma supply from larger radii, and the particles accelerate toward the star with secondary pair production, in repeating episodes.The separatrix return-current (Fig. 4, thick black curve) has α < 0 and is sustained by pair production at the Y-point as explained in (Hu & Beloborodov 2022).Pair production along and within the separatrix itself is inefficient because it is very thin and curved, so the emitted gamma-rays quickly travel outside its thickness before converting to e ± .As a result, the separatrix return-current does not display the cyclic e ± discharge episodes.

B. DISPERSION RELATION AND NORMAL MODES OF PULSAR PLASMA
Here we review the dispersion relation and normal modes of the e ± pair plasma expected to exist in pulsar magnetospheres.Near the neutron star surface the plasma is effectively infinitely magnetized.Therefore, charged particles only move parallel to the magnetic field, and the plasma supports a very limited set of normal modes (Arons & Barnard 1986).The dispersion relation is obtained from the Vlasov-Maxwell equations (Arons & Barnard 1986), where ω is the wave frequency, ω p the plasma frequency, k || and k ⊥ are the components of the wavevector parallel and perpendicular to B B B, g = ⟨γ −3 (1 − kv/ω) −2 ⟩, and ⟨...⟩ signifies an average over particle momenta.If the plasma is cold, g = 1 in the plasma rest frame.The dispersion relation describes three modes: The extraordinary (X) mode with ω = ck and electric field polarized in the k k k × B B B direction, and two branches of the ordinary (O) mode with electric field polarized in the k k k − B B B plane.In the dense discharge plasma (ck/ω p ≪ 1), the polarization of the two branches of the O-mode is described by where δE || is the wave electric field parallel to B B B, and δE ⊥ is the wave electric field perpendicular to B B B in the k k k − B B B plane (Arons & Barnard 1986).The Alfvén branch has phase speed ω/k < c, the superluminal O-mode mode ω/k > c, and the extraordinary mode has ω/k = c.For a given k k k, the three modes (extraordinary, Alfvén, and superluminal) are orthogonal.The Alfvén mode suffers from Landau damping, and cannot escape the magnetosphere without non-linear conversion to a different eigenmode (Arons & Barnard 1986).By contrast, as the superluminal mode propagates outward through the decreasing plasma density, it tracks adiabatically along the superluminal branch and becomes a freely propagating electromagnetic wave (Arons & Barnard 1986;Barnard & Arons 1986).During the propagation the wave frequency ω ⋆ is fixed, while k increases and approaches k ⋆ = ω ⋆ /c (the phase speed approaches c from above).The superluminal mode is attractive for pulsar radio emission because it does not suffer from Landau damping (ω/k > c), and it can escape the magnetosphere as a vacuum electromagnetic wave in the radio band (Philippov et al. 2020).

C. BEAM INSTABILITY IN PULSAR PLASMA
Here we review the relativistic beam instability in highly magnetized e ± plasma (for a detailed calculation see (Egorenkov et al. 1983)).Consider a uniform background of e ± with uniform density n p which streams relativistically along the magnetic field with Lorentz factor γ. The beam has uniform density n b < n p , and a Gaussian distribution of momentum parallel to the magnetic field with mean p b = γ b mc, and width ∆p b = ∆γ b mc.If the beam density is sufficiently high or ∆γ b is sufficiently small, the instability operates in the hydrodynamic regime with a growth rate Γ ≫ kc ∆γ b /γ 3 b .Then, the growth rate can be found by setting ∆γ b → 0. In particular, for purely longitudinal oscillations the dispersion relation has the form (Eq. B1) where the third term is due to the beam with ω b = 4πn b e 2 /m and v Egorenkov et al. 1983).Near Cherenkov resonance, the dispersion relation can be written in the form ω = kv b + ∆ω, with |∆ω| ≪ kv b .The growth rate is maximal when k approaches k 0 ≡ 2ω p ⟨γ⟩ 1/2 /c, and is given by The condition for the instability to operate in the hydrodynamic regime is then This condition may be satisisfied even when the momentum spread of the beam is large, ∆γ b ≳ γ b .Note that the instability also occurs at k ≪ k 0 , although at a slower rate (Egorenkov et al. 1983).The beam instability grows only the subluminal (Alfvén) mode; the superluminal mode cannot satisfy the Cherenkov resonance ω/k = v b .Beam instabilities in pulsar magnetospheres are usually discussed in the context of primary accelerated e ± with γ b ∼ 10 6 (beam) interacting with the secondary e ± with γ ∼ 10 2 -10 3 (background plasma) (Ruderman & Sutherland 1975).The growth rate (Eq.C4) is then strongly suppressed by the high γ b , resulting in Γ ≪ Ω, where Ω is the rotation rate of the pulsar.The plasma escapes the magnetosphere on the timescale r LC /c = Ω −1 , and so there is not enough time for the instability to develop.
By contrast, our simulation reveals locations where both the background plasma and the beam consist of secondary e ± : there is a wind of e ± produced near the star and flowing out to the light cylinder and e ± beams injected by the gamma-rays from the equatorial current sheet.The moderate Lorentz factors of the created e ± , γ b ∼ 10 2 -10 3 , imply a fast growth of the beam instability Γ > Ω.Indeed, using a typical γ ∼ 10 2 for the wind and γ b ∼ 10 3 for the injected beam, we find where M = n p /n GJ is the pair multiplicity, P = 2π/Ω, B LC is the magnetic field strength at the light-cylinder, and we set n b /n p ∼ 0.1.We conservatively normalized M to a modest value of 10 2 , and a higher M will only increase Γ/Ω.The beam instabilities require pair production by photon-photon collisions around the equatorial current sheet, which occurs in rapidly rotating pulsars; therefore, we normalized B LC to 10 6 G, as found e.g. in the Crab pulsar.Note also that for extremely energetic pulsars the e ± created around the current sheet may dominate the wind density and serve as the main plasma background; then, the wind emerging from the inner magnetosphere will play the role of the beam.The beam instabilities around the equatorial current sheet will lead to a new component of pulsar radio emission if there is significant conversion of the generated subluminal O-mode to another eigenmode, which can escape.

Figure 1 .
Figure 1.Poloidal cross section of the magnetosphere in the quasi-steady state.Green curves show poloidal magnetic field lines, and background color shows plasma density.The clouds of plasma above the polar cap form in e ± discharges.The absence of plasma (dark zone) above the separatrix shows the return-current gaps in the unscreened state; it fills with plasma every ∼ 2RLC/c during discharge episodes.Magnetic reconnection forms a dynamic plasmoid chain in the equatorial current sheet.

Figure 2 .
Figure 2. Excitation of waves in the pulsar magnetosphere.Color bar shows the components of the electromagnetic field multiplied by (r/r⋆) 2 .Top : E || in the quasi-steady state.The dashed white curve shows the null-surface ρGJ = 0 (the nullsurface breathes outside rLC due to plasmoid motion).Insets (i)-(iii): Electromagnetic modes excited during e ± discharges at the polar-cap, inner return-current, and outer return-current gaps.Inset (iv): Waves excited by a beam plasma instability.For insets (i) and (ii), δE ⊥ is calculated by projecting δE E E ⊥ onto a unit vector in the φ φ φ × B B B direction.
Fig. 2 shows the plasma waves excited during the polar-cap discharge [inset (i)], the inner return-current discharge [inset (ii)], and the outer return-current discharge [inset (iii)].The electromagnetic fields of the waves δB B B and δE E E are isolated by first subtracting the time average to remove the zero frequency (background)

Figure 3 .
Figure 3. Beam instabilities in the outer magnetosphere.Left: Parallel momentum (p || ) distributions of e ± (grey) and gamma-rays (cyan) in regions (a) and (b) in the right panel.For gamma-rays we show distributions of p || /2 (their momentum is split between the created e + and e − ).(a): Gamma-rays from the current-sheet (cyan peak) inject e ± (grey peak).(b): Gamma-rays from the separatrix (thin cyan peak) and current-sheet (wider cyan peak) inject different e ± beams.Right: Charge density ρ and charge bunches.
The authors thank A. Chen, R. Hui, A. Philippov, and M. Ruderman for useful discussions.Resources supporting this work were provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center.A.B. is supported by a PCTS fellowship and a Lyman Spitzer Jr. fellowship.A.M.B. acknowledges grant support by NSF AST 2009453, NASA 21-ATP21-0056, and Simons Foundation #446228.Y. L. is supported by NSF Grant AST-2009453.APPENDIX A. GLOBAL CURRENTS AND PAIR PRODUCTION

Figure 4 .
Figure 4. Global currents and e ± production in the aligned rotator.The white region above the polar cap is the zone of outflowing current, the shaded grey region indicates the bulk return-current, and the thick black curve shows the separatrix returncurrent layer.The hatched area indicates the closed zone.The null charge line (dashed) is defined by ρGJ = −Ω Ω Ω • B B B/(2πc) = 0; above the line ρGJ < 0, and below ρGJ > 0. Episodic e ± discharge occurs in the out-flowing current above the polar cap (α > 1), in the bulk return-current above the polar cap (α < 0), and in the bulk return-current outside the null line (α > 1).Beam instabilities occur outside the light cylinder where e ± are created by the intersecting gamma-ray beams with different directions.