Pulsar Radio Emission Mechanism: Radio Nanoshots as a Low-frequency Afterglow of Relativistic Magnetic Reconnection

The cause of pulsar radio emission has remained a mystery for nearly half a century. This is mainly because the radio emission is coherent, i.e., produced by microscopic magnetospheric collective plasma motions that are extremely difficult to calculate analytically. The recent discovery of coherent radiation from fast radio bursts (FRBs; Lorimer et al. 2007) underscores the role of pulsars as keys to understanding astrophysical coherent emission, and warrants a re-evaluation of the radio emission mechanisms with modern techniques.

Pulsar magnetospheres have recently been successfully studied numerically using a first-principle particle-in-cell (PIC) approach. PIC simulations helped uncover the conditions for the operation and intermittency of the polar cap discharge (Timokhin & Arons 2013; Chen & Beloborodov 2014; Philippov et al. 2015a), located a high-energy emission source at the current sheet beyond the light cylinder (LC; Cerutti et al. 2016), and explained the general morphology of gamma-ray profiles (Philippov & Spitkovsky 2018). They also provided clues about potential mechanisms of coherent radio emission. First, the non-stationary screening of the electric field by pair cascade at the polar cap could lead to the production of coherent electromagnetic (EM) waves (Timokhin & Arons 2013; Philippov et al. 2015a; A. Philippov et al. 2019, in preparation). Second, the return separatrix current layer bounding the closed-field zone supports multi-streaming particle distributions (Brambilla et al. 2018; Philippov & Spitkovsky 2018), which could be unstable to anomalous cyclotron instability (Lyutikov et al. 1999, A. Philippov & M. Lyutikov 2019, in preparation). Finally, magnetic reconnection in the current sheet beyond the LC proceeds in the highly non-stationary, plasmoid-dominated regime (Cerutti & Philippov 2017). Uzdensky & Spitkovsky (2014) suggested that rapid coherent motions of plasmoids in the hierarchical chain could contribute to radio emission, albeit without proposing a concrete physical mechanism.

The global magnetospheric structure suggests that beamed emission at the tip of the separatrix layer maps into the same caustic as the current-sheet emission (Bai & Spitkovsky 2010), so the three possible emission mechanisms described above lead to two potential sources of radio emission: emission produced in the non-stationary polar discharge, which we call Type I, and emission produced in the outer magnetosphere, which we call Type II. Type II emission should naturally coincide in phase with the magnetospheric gamma-ray emission, which also originates in the current layer (Lyubarskii 1996; Uzdensky & Spitkovsky 2014; Cerutti et al. 2016; Kalapotharakos et al. 2018). This emerging picture is motivated by the recent observations of the Crab pulsar, several other young energetic pulsars, and a few millisecond pulsars, which clearly show a radio emission source that is coincident with the location of the pulsed high-energy radiation source (Johnston & Romani 2004; Johnson et al. 2014; Eilek & Hankins 2016). In the Crab pulsar, this radio emission component (identified here with Type-II emission) is sometimes resolved into discrete bright pulses of nanosecond duration, "nanoshots" (Eilek & Hankins 2016). These nanoshots may form entirely the Crab's main pulse and low-frequency interpulse, and appear to come in packets making up microsecond-long bursts.

In this Letter we present a physical model of coherent Type II radio emission produced by merging plasmoids in the reconnecting current sheet beyond the pulsar's LC. In Section 2 we discuss the plasmoid instability in global PIC magnetospheric simulations. In Section 3 we present local PIC simulations of reconnection, which allow us to calculate the properties of coherent low-frequency EM radiation emitted by merging plasmoids. We summarize our model's observational implications and conclude in Section 4.

While this Letter was in preparation, the same study had been suggested independently by Y. Lyubarsky (Lyubarsky 2019). Our work complements his paper with simulations and shows that relativistic reconnection is indeed a powerful source of low-frequency waves, which for some pulsars are in the GHz band.

Numerical simulations of pulsar magnetospheres, starting with the pioneering work of Contopoulos et al. (1999), showed the formation of a prominent current sheet beyond the LC. Global 2D (Cerutti & Philippov 2017) and 3D PIC simulations (Philippov et al. 2015b; Cerutti et al. 2016; Philippov & Spitkovsky 2018) also confirmed that this current sheet is unstable to plasmoid instability. However, because of their limited scale separation between the current sheet's length, comparable to the LC radius, (${R}_{\mathrm{LC}}$), and its width, $\delta \sim {r}_{L}$ (here, r_L is the average Larmor radius of particles in the sheet), 3D simulations typically show only one to two plasmoids. In contrast, 2D simulations of the striped wind that are limited to the equatorial plane allow significantly greater scale separation and hence show more plasmoids (Cerutti & Philippov 2017).

To probe the expected number of plasmoids in the current sheet in the global pulsar magnetosphere as a function of plasma parameters, we simulate the structure of the striped wind using the setup of Cerutti & Philippov (2017). Our simulations are performed in spherical coordinates $(r,\theta ,\phi )$, and are restricted to the equatorial plane ($\theta =\pi /2$). The pulsar's LC is at $r={R}_{\mathrm{LC}}=3{R}_{* }$, where R_* is the stellar radius. We inject neutral plasma at the pulsar surface and give it a nonrelativistic kick along the magnetic field. We run this setup for several values of the plasma multiplicity, $\kappa =1,10,100$, at $r={R}_{* }$. Multiplicity is defined as $\kappa =n/{n}_{\mathrm{GJ}}$, where n is the injected plasma density and ${n}_{\mathrm{GJ}}={{\rm{\Omega }}}_{* }B/2\pi {ce}$ is the fiducial Goldreich–Julian plasma density (here, ${{\rm{\Omega }}}_{* }$ is the stellar angular velocity and B is the local magnetic field strength). We choose the magnetic field at the star to be sufficiently strong so that even the simulation with the largest injected multiplicity, $\kappa =100$, has a sufficiently high value of the magnetization parameter at the LC, ${\sigma }_{\mathrm{LC}}\equiv {B}_{\mathrm{LC}}^{2}/4\pi \kappa {n}_{\mathrm{GJ}}{m}_{e}{c}^{2}\approx 70$ (hereafter, ${B}_{\mathrm{LC}}$ is the magnetic field strength at the LC). Our simulations include the synchro-curvature radiation reaction on particle motion (see Cerutti & Philippov 2017 for details).

Simulations quickly settle to a solution with closed and open field lines, and a current sheet beyond the LC, separating regions of oppositely directed magnetic field. The typical structure of the plasma density is shown in Figure 1. We find that the current sheet width, δ, scales with the Larmor radius of particles heated up to energies $\langle \gamma \rangle \simeq {\sigma }_{\mathrm{LC}}$, so that $\delta \propto {r}_{L}\propto {m}_{e}{c}^{2}{\sigma }_{\mathrm{LC}}$/${{eB}}_{0}\propto {R}_{\mathrm{LC}}{\kappa }^{-1}$. For our κ = 1 simulation, the ratio $\delta /{R}_{\mathrm{LC}}$ is substantial, $\delta /{R}_{\mathrm{LC}}\approx 0.15$, and the sheet is stable to tearing. A simulation with $\kappa =10$ develops a current sheet with $\delta /{R}_{\mathrm{LC}}\approx 0.026$ and shows few plasmoids, while the κ = 100 simulation has $\delta /{R}_{\mathrm{LC}}\approx 0.0028$ and develops a vigorous plasmoid chain (see Figure 1). The zoom-in into the chain at $2{R}_{\mathrm{LC}}$ is shown in Figure 1(a). Overall, we find that the number of plasmoids within 2R_LC scales approximately as R_LC/δ, but we caution that our dynamic range is too small to predict an asymptotic scaling. However, it is clear that for real parameters of the Crab pulsar, this would correspond to very large numbers of plasmoids between one and two R_LC. We also observe plasmoid mergers at $r\lesssim (1-3){R}_{\mathrm{LC}}$, with the most frequent mergers of small plasmoids occurring close to the LC. As plasmoids grow when they propagate away from the LC, we observe mergers of bigger plasmoids at larger distances. We note that mergers of all sizes should be much more frequent for realistic parameters. Beyond r ≈ 4R_LC the flow expands nearly radially, and the plasmoid structure freezes out.

Zoom In Zoom Out Reset image size

Figure 1(b) shows the toroidal magnetic field component, exhibiting a wave-like structure (a light blue filament passing through the point with r ≈ 2.2R_LC, ϕ ≈ 3.1) coming from large plasmoids. As discussed below, these waves are fast waves emitted in the process of plasmoid mergers. We compute the light curve of this low-frequency radiation by measuring the Poynting flux at the simulation boundary, r = 6R_LC and ϕ = 0, and subtracting the "DC" components of the EM field that correspond to the bulk Poynting flux in the pulsar wind. The result is a sequence of bright short pulses with variable arrival times clustered around rotational phases 0.3 and 0.8 (see Figure 1(c)). The light curve averaged over a few pulsar rotations is shown as a black solid line; it is significantly narrower than the light curve of the incoherent high-energy synchrotron emission of particles accelerated in the current sheet (black dashed line, reproduced from Cerutti & Philippov 2017). As the synchrotron emission is in the gamma-ray band for young energetic pulsars (Uzdensky & Spitkovsky 2014; Hakobyan et al. 2018), we find that the peaks of coherent low-frequency emission coincide with the peaks of the gamma-ray light curve.⁵ Below we argue that the characteristic frequency of the low-frequency emission from plasmoid mergers falls in the GHz range for parameters of the Crab and other pulsars with high B_LC.

3.1. Numerical Setup

3.2. Results

Reconnection starts as the tearing instability breaks up the current layer into a chain of plasmoids. As these primary plasmoids move along the sheet, secondary current sheets between them become sufficiently elongated and also become unstable, yielding a hierarchical fragmented structure (see Figure 2(a) for the snapshot of magnetic field lines and plasma density). Over time, plasmoids grow and merge into larger ones. Current sheets perpendicular to the main current sheet form between sufficiently large merging plasmoids. These sheets also become unstable to secondary tearing and form secondary plasmoids as shown in the density map in Figure 2(b). Both small plasmoids in the main sheet and plasmoids in the secondary perpendicular current sheets can reach relativistic velocities. As they merge with each other or collide with the upstream magnetic field (as is the case for plasmoids in secondary perpendicular current sheets), coherent time-dependent currents are generated at the interface and launch a powerful EM wave that lasts until the merger is complete. These waves propagate in the cold upstream plasma across the magnetic field, and the wave electric field is perpendicular both to the background magnetic field and to the direction of propagation. This polarization corresponds to the extraordinary mode, which is the fast magnetosonic wave in the strongly magnetized pair plasma (Arons & Barnard 1986). As the plasma density goes to zero, these waves become vacuum EM waves. They are unaffected when they cross the buffer transition at x = 0.2L_x or x = 0.8L_x. The electric field of the waves is then recorded at the edge of the simulation domain, x = 0 and x = L_x.

Zoom In Zoom Out Reset image size

In Figure 2(c) we show a map of the x-component of the Poynting flux in a simulation snapshot, where we subtracted the average reconnection-driven inflowing Poynting flux, ${S}_{0}\approx 0.1{{cB}}_{0}^{2}/4\pi$. The circular spikes correspond to outward-propagating EM waves generated in plasmoid mergers and in collisions of plasmoids with the upstream field. Each individual spike is a half-wavelength oscillation, so that the spike duration τ and the associated frequency ν are related as τν ≃ 1.

Figure 3 shows the evolution of the electric field, normalized by the reconnection electric field ${E}_{0}=\langle {E}_{z}{\rangle }_{y}\simeq 0.1{B}_{0}$, at the x = 0 boundary of the box at y = L_y/2. Here, regions of E_z < 0 correspond to waves that carry energy away from the layer. In the bottom panel we show a zoom into a short time interval, when the waves from several different mergers reach the simulation boundary (only negative values of E_z shown for clarity). The pulse, e.g., a train of several spikes between t = 3600–3800ω_pt, comes from a merger of two large plasmoids, which produced several smaller plasmoids in the perpendicular sheet between them. Each individual merger of a small plasmoid with the upstream field generates a single electric-field spike. The overall pulse duration is τ ≈ 10/ν, where ν is the frequency associated with individual spikes. This relation is consistent with the observed duration of individual radio nanoshots from the Crab pulsar (Eilek & Hankins 2016).

Zoom In Zoom Out Reset image size

The duration distribution of single spikes and pulses⁶ is shown in Figure 4. It peaks around τ ∼ 10/ω_p, where ω_p is the background plasma frequency, and resembles a power law ${dN}/d\tau \approx {({\omega }_{p}\tau )}^{-\alpha }$ for larger durations, and α is between 1.5 and 2. The pulse-averaged spectral density of the EM radiation in 3D can be then estimated as ${S}_{\nu }d\nu \propto \left({ \mathcal E }/\tau \right)({dN}/d\tau )$ $d\tau \propto {\tau }^{2-\alpha }d\tau$, where we used τν ≃ 1, and estimated the energy of the individual pulse, ${ \mathcal E }$, as the magnetic energy stored in the corresponding merging plasmoid of size ℓ, which scales as the plasmoid's volume $\sim {{\ell }}^{3}\sim {(c\tau )}^{3}$. The obtained steep spectrum, ${S}_{\nu }\propto {\nu }^{-\beta }$, where $\beta =(4-\alpha )$ is between 2 and 2.5, implies that most of the energy is emitted as low-frequency radiation, produced by mergers of moderate-sized plasmoids. This spectrum is compatible with the observed β ≈ 3, for the main component of the Crab pulsar's radio profile (Moffett & Hankins 1999), although the exact value of β varies substantially between the pulses (Karuppusamy et al. 2010).

Zoom In Zoom Out Reset image size

While the occasional brightest spikes have instantaneous amplitudes of up to 3E₀, the electric field in typical pulses is about E₀, implying typical instantaneous Poynting fluxes of the EM waves as high as ${{cE}}_{0}^{2}/4\pi \approx 0.01{{cB}}_{0}^{2}/4\pi$. However, when averaged over the simulation duration T, the Poynting flux of the outgoing EM wave emission, measured at x = 0 boundary of the box at y = L_y/2, is much smaller, $(c/4\pi ){\int }_{0}^{T}{E}_{z}({B}_{y}-\langle {B}_{y}{\rangle }_{y}){dt}$/$T\approx {10}^{-4}{{cB}}_{0}^{2}/4\pi$. This implies that only a small fraction of the reconnection-dissipated power is transferred into coherent EM waves. Most of this coherent power is radiated at low frequencies, which may not necessarily fall in the observed GHz range in the case of pulsars. Thus, the conversion efficiency into observable radio waves could be even lower. We note that the average efficiency in our simulations is likely significantly exaggerated, because the periodic boundary conditions produce favorable conditions for plasmoid mergers. Mergers are still quite regular in simulations with outflowing boundaries (Sironi et al. 2016); however, it appears hard to quantify the average efficiency of the process in real pulsars with our idealized local calculations. Although the basic picture of the current sheet breaking into plasmoids as a result of tearing instability stays the same in three dimensions (Werner & Uzdensky 2017), our quantitative conclusions remain to be checked with future 3D simulations. Nonlinear interactions between the fast waves may further decrease the average efficiency of the escaping radiation (Lyubarsky 2019).

We now discuss implications of our results for observations of pulsar radio emission. In the Crab pulsar magnetosphere, the comoving plasma density and temperature inside the reconnecting current sheet are regulated by cross-layer pressure balance and the energy balance between reconnection-powered heating and synchrotron cooling, which is about 10¹³ cm⁻³ and 10 GeV, respectively (Lyubarskii 1996; Uzdensky & Spitkovsky 2014). Then, the typical thickness δ of elementary current layers is of order 1 m (Uzdensky & Spitkovsky 2014; Lyubarsky 2019), and the typical size of small plasmoids, about 10–100 times larger, is of order 30 m. This translates into a comoving duration of coherent radio-spikes of about 100 ns. After accounting for relativistic Doppler effect due to pulsar-wind bulk Lorentz factor Γ ∼ 10–100 (see, e.g., Lyubarskii 1996), this can explain the observed GHz range of pulsar radio emission (Uzdensky & Spitkovsky 2014; Lyubarsky 2019) and, in particular, ns-scale pulses (nanoshots). Pulsars with B_LC weaker than in the Crab have thicker current sheets and hence larger secondary plasmoids. Thus, coherent EM emission from plasmoid mergers in these systems may occur at frequencies lower than the traditional ≳100 MHz radio range, of interest for possible detections with the Low-Frequency Array and the Square Kilometre Array. This may explain why radio emission phase aligned with high-energy emission is only observed in pulsars with high B_LC (Johnston & Romani 2004).

It is remarkable that a nanosecond-long radio spike produced from a region just tens of meters across in a single plasmoid merger in the magnetosphere of the Crab pulsar 2 kpc away is so powerful that it can be detected on Earth. This is perhaps the smallest source of observable emission in all of astrophysics. Indeed, the amount of magnetic energy released in an individual merger event is about ${ \mathcal E }\sim ({B}_{\mathrm{LC}}^{2}/8\pi ){{\ell }}^{3}\sim {10}^{21}\,\mathrm{erg}$, where we took B_LC = 1 MG and ℓ = 30 m. If, say, 1% of this energy is emitted as a coherent EM pulse with a duration τ of a few ns-long wave-periods (i.e., τν ∼ few), one can estimate the observed spectral flux density (for a distance of d = 2 kpc) of order ${S}_{\mathrm{obs}}\sim 0.01{(\pi {d}^{2})}^{-1}({ \mathcal E }/\tau \nu ){{\rm{\Gamma }}}^{3}$ ∼ ${10}^{-25}{{\rm{\Gamma }}}^{3}$ $\mathrm{erg}/({\mathrm{cm}}^{2}\,{\rm{s}}\,\mathrm{Hz})$ ∼ $0.01\,{{\rm{\Gamma }}}^{3}\,\mathrm{Jy};$ taking Γ ∼ 30, we get S_obs ∼ 300 Jy, comparable to the observed values (Eilek & Hankins 2016, see also Lyubarsky 2019). Here we assumed that this emission is isotropic in the pulsar wind frame, which is the case in our idealized 2D reconnection simulation. The corresponding brightness temperature of 0.3 GHz radio emission is about ${T}_{B}\sim {S}_{\mathrm{obs}}{c}^{2}$/$(2{k}_{B}{\nu }^{2}\delta {\rm{\Omega }})$ ∼ ${10}^{38}\,{\rm{K}}\,{({\rm{\Gamma }}/30)}^{3}$, where $\delta {\rm{\Omega }}\sim \pi {{\ell }}^{2}/{d}^{2}$ is the solid angle covered by the colliding plasmoid, which is well in line with the observed values (Hankins et al. 2003).

To conclude, our model of the radio nanoshots explains their observed frequency, duration, phase alignment with the high-energy emission, extremely high instantaneous brightness, and low average efficiency. Our simulations do not address the sometimes strong circular polarization of nanoshots (Eilek & Hankins 2016). The observed circular polarization may arise due to propagation effects in the magnetospheric plasma (Lyubarskii & Petrova 1998; Beskin & Philippov 2012), which for high-B_LC pulsars imply the polarization formation radius located beyond the LC.

Given the right physical conditions, our results may have implications for coherent radio emission from other sources where relativistic reconnection may take place, e.g., magnetar magnetospheres and interacting magnetospheres of binary neutron stars before the merger. These possibilities and potential connections to FRBs will be investigated elsewhere.

We thank V. Beskin, Y. Lyubarsky, and M. Lyutikov for stimulating discussions. The Flatiron Institute is supported by the Simons Foundation. This research was supported by the NSF (grant AST-1411879), NASA ATP (grants NNX16AB28G, NNX17AK57G, and 80NSSC18K1099), CNES and the Université Grenoble Alpes (IDEX-IRS grant). D.U. acknowledges the hospitality of the Institute for Advanced Study and the Ambrose Monell Foundation for supporting his 2016–2017 sabbatical, when most of this work was done. We acknowledge PRACE for awarding us access to Curie at GENCI@CEA, France.