Coherent Emission from QED Cascades in Pulsar Polar Caps

Cascades of electron–positron pairs are a key source of plasma in pulsar magnetospheres. They result from a positive feedback loop that develops in vacuum gaps, regions of unscreened electric field in the magnetosphere. Polar caps have been proposed to host rotation-induced vacuum gaps (Sturrock 1971; Ruderman & Sutherland 1975), where TeV energy electrons and positrons emit gamma-ray photons via curvature radiation and these are absorbed in the local ∼10¹²G magnetic field to produce new pairs. Cascades cease when the fresh pair plasma screens the vacuum gap electric field. The plasma is then advected into the magnetosphere, the gap reopens and a new cascade begins (e.g., Timokhin 2010). The time-dependent dynamics of polar cap vacuum gaps has recently been proposed as a primary ingredient to explain the nature of pulsar radio emission (Beloborodov 2008; Timokhin & Arons 2012; Melrose et al. 2020; Philippov et al. 2020).

Kinetic plasma simulations are ideally suited to study the highly nonlinear interplay between time-dependent pair cascades and coherent plasma processes. Timokhin (2010) presented the first 1D particle-in-cell (PIC) simulations including curvature radiation and pair production from first principles. These simulations were performed in a frame corotating with the neutron star, which embeds the magnetic field twist imposed at macroscopic scales (Bai & Spitkovsky 2010) via a background current density j_m = α ρ_GJ c, where α is a constant of order unity, ρ_GJ is the Goldreich–Julian (Goldreich & Julian 1969) charge density, and c is the speed of light. They showed that cascades develop over regular short bursts followed by long quiet phases during which no pairs are produced. As the gap is locally screened, electric field oscillations are inductively driven due to collective kinetic-scale plasma motions (Levinson et al. 2005; Cruz et al. 2021). Similar simulations have been used by Timokhin & Arons (2012) and Timokhin & Harding (2015) to determine the spectra and multiplicity of the pair plasma created in cascades for a variety of initial conditions.

Philippov et al. (2020, hereafter PTS20) have recently resorted to an heuristic description of the emission and pair production processes to perform the first 2D Cartesian PIC simulations of pair cascades and identified a new process of coherent radiation emission. The key ingredient for this process is a finite angle between the pair production front and the background magnetic field, which cannot be captured in 1D. While PTS20 have been able to demonstrate this coherent emission mechanism using a simplified configuration, the shape of the pair production front and the spectrum and Poynting flux profile of the emitted radiation in realistic field geometries depend on the microscopic details of the emission and pair production cross sections. First-principles simulations of polar cap pair cascades are very challenging computationally, due to (i) the rapid and localized creation of a large number of particles and (ii) the large separation between gap and plasma kinetic scales. While the polar cap vacuum gap extends for ∼100 m (Ruderman & Sutherland 1975), the electron skin depth associated with the dense plasma produced in the cascade is ∼1 cm. This scale separation is a consequence of the large multiplicity of the cascade process that can be estimated as the ratio between the energies of primary and secondary particles, respectively ε_±/m_e c² ∼ 10⁷ and ${\varepsilon }_{\pm }^{{\prime} }/{m}_{e}{c}^{2}\sim {10}^{2}$, where m_e is the electron mass. In multidimensional configurations, the difficulties of simulating such scale disparity are heightened, and first-principles numerical models of the pulsar polar cap have not yet been possible.

In this Letter, we adopt a first-principles rescaling of the quantum electrodynamics (QED) processes responsible for gamma-ray and pair production processes in 2D axisymmetric PIC simulations of pulsar polar caps with a realistic field geometry. Using these simulations, we determine the multidimensional properties of pair cascades and their observational signatures.

We perform 2D QED–PIC simulations with the code OSIRIS (Fonseca et al. 2002, 2008) in axisymmetric cylindrical coordinates (R, z). The lower z boundary (z = 0) is a rotating conducting disk of radius R₀. The reference scale R₀ should be interpreted as the polar cap radius, and the z-axis as the magnetic (and rotation) axis of the neutron star. The rotation of the conducting disk is imposed by forcing the radial electric field at this boundary to be E_R (R, z = 0) = E₀ × (R/R₀) × g(R/R₀), where E₀ is a constant proportional to the angular velocity of the disk Ω and $g(x)=0.5\times (1-\tanh ((x-1)/0.2))$ is a smooth cutoff function at x ≃ 1. This boundary condition induces a unipolar electric field near the conductor, forcing plasma in this region into corotation with the disk. The upper z boundary is open for fields, such that any incident wave escapes freely from this boundary. Both z boundaries are open for particles.

The simulation domain is also permeated by the externally imposed dipolar magnetic field B _d with magnitude B₀ at R = z = 0. The center of the magnetic dipole is at (R, z) = (0, − R_*), with R_*/R₀ = 10. In analogy with a realistic polar cap, all magnetic field lines that cross the lower z boundary at a radius R < R₀ are open to infinity, and those that cross this boundary at R > R₀ are closed on the other hemisphere of the neutron star. We assume that closed field lines are filled with dense plasma that we model as an ideal conductor, i.e., we set E _∥ = ( E · B _d) B _d/∣ B _d∣² in this region to zero. This simulation setup forces a current to be driven along open field lines while providing a (virtual) return current along the edge of the conducting disk, mimicking the local conditions in pulsar polar caps.

The QED processes governing pair cascades in our simulations are photon emission via nonlinear Compton scattering (Erber 1966), the QED equivalent of curvature radiation for classical emission from ultrarelativistic particles (Kelner et al. 2015; Del Gaudio 2020), and multiphoton Breit–Wheeler pair production (Ritus 1985). In our simulations, these processes are included using Monte Carlo methods (Grismayer et al. 2016, 2017). In each time step, we first compute the quantum parameter χ_±,γ of each particle (subscripts ±, γ correspond to electrons, positrons, and photons, respectively), defined as ${\chi }_{\pm ,\gamma }=\sqrt{{({p}_{\mu }{F}^{\mu \nu })}^{2}}/({B}_{Q}{m}_{e}{c}^{2})$, where p^μ is the four-momentum of the particle, F^{μ ν} is the electromagnetic tensor, and B_Q ≃ 4.4 × 10¹³ G is the Schwinger field. We then rescale ${\chi }_{\pm }\to {\chi }_{\pm }^{{\prime} }\equiv {\zeta }_{\pm }{\chi }_{\pm }$ and ${\chi }_{\gamma }\to {\chi }_{\gamma }^{{\prime} }\equiv {\zeta }_{\gamma }{\chi }_{\gamma }$, with ζ_±,γ ≫ 1, and use ${\chi }_{\pm ,\gamma }^{\prime}$ to evaluate the probability of creating new particles. The rescaling is done to reduce the scale separation described above. A key aspect of our approach is that it retains the fundamental properties of the QED processes: (i) the spectrum of photons emitted via nonlinear Compton scattering is preserved (in particular its dependence on the energy and curvature of the trajectory of the emitting particle), and (ii) the probability of pair production critically depends on the angle between the photon propagation direction and the local electromagnetic field.

The surface electric field E₀ = ΩB₀ R_*/c and the rescaling parameters ζ_±,γ are adjusted such that a lepton is accelerated to the radiation reaction limited Lorentz factor (Daugherty & Harding 1982)

$$\begin{eqnarray}&&{\hat{\gamma }}_{\pm }={\left(\displaystyle \frac{3}{2{\alpha }_{\mathrm{fs}}}\displaystyle \frac{{E}_{0}}{{B}_{Q}}\right)}^{1/4}{\left(\displaystyle \frac{{\bar{\lambda }}_{{\rm{C}}}}{\rho }\right)}^{-1/2}{\zeta }_{\pm }^{-1/2}\end{eqnarray} \tag{ 1 }$$

over a distance

$$\begin{eqnarray}&&{{\ell }}_{a}\simeq {\hat{\gamma }}_{\pm }\left(\displaystyle \frac{{B}_{Q}}{{E}_{0}}\right){\bar{\lambda }}_{{\rm{C}}}\end{eqnarray} \tag{ 2 }$$

that we take to be ∼0.1R₀. In Equations (1) and (2), α_fs is the fine-structure constant, ${\bar{\lambda }}_{{\rm{C}}}$ is the reduced Compton wavelength, and $\rho \simeq {R}_{* }^{2}/{R}_{0}$ is the radius of curvature of the last open field line. We fix ζ_±,γ = 10³, eB₀ R₀/m_e c² = 10⁸, and B₀/B_Q = 0.1, and vary only E₀ such that ℓ_a /R₀ ≃ {0.2, 0.35, 0.5, 0.6}. These values were chosen for simulations to be computationally feasible yet comparable to real pulsars, in which ℓ_a /R₀ varies from ∼10⁻⁴ to ∼1 with increasing rotation period. The mean free path of photons emitted at the critical energy ${\varepsilon }_{c}=(3/2)({\bar{\lambda }}_{{\rm{C}}}/\rho ){\hat{\gamma }}_{\pm }^{3}{m}_{e}{c}^{2}{\zeta }_{\pm }$, defined as (e.g., Timokhin & Harding 2015)

$$\begin{eqnarray}&&{{\ell }}_{\gamma }\simeq \displaystyle \frac{8}{3}\displaystyle \frac{{B}_{Q}}{B}\displaystyle \frac{{m}_{e}{c}^{2}}{{\varepsilon }_{c}}\rho {\zeta }_{\gamma }^{-1},\end{eqnarray} \tag{ 3 }$$

is always ∼10⁻² R₀ for our choice of parameters. Thus, the ratio ℓ_γ /ℓ_a ∼ 0.1 is also within the range ∼10⁻²–1 expected in real systems. The electron skin depth associated with a density n_GJ = ∣ρ_GJ∣/e ≡ ΩB₀/2π ce is ∼5 × 10⁻³–10⁻² R₀. The simulation domain has a size L_R × L_z = (1.5R₀) × (2.5R₀) discretized in N_R × N_z = 6000 × 10,000 cells, and the time step is Δt = 10⁻⁴ R₀/c. The grid size was chosen to resolve the electron skin depth associated with the maximum density generated during pair cascades.

The rotating conducting disk is gradually spun up to its maximum angular velocity in the first 200 time steps of the simulations and kept constant thereafter. There is initially no plasma in the simulation domain, so this spin induces a vacuum corotation electric field in the open field line region: E_R increases and E_z decreases with R, whereas both components decay with z within a distance ∼R₀. We let the corotation field develop for a light-crossing time, L_z /c = 2.5 R₀/c, and then start injecting plasma in cells just above the z = 0 boundary for all R < R₀. At every time step, we inject one electron–positron pair per cell carrying a density n_inj = κ(E_∥/eR₀), where κ = 0.1 and E_∥ = ∣ E _∥∣. Pairs are injected at rest. We choose a value of κ ≪ 1 to ensure this injection provides only a seed plasma for the cascades and does not dominate the plasma outflow.

As pairs are injected and experience the vacuum electric field, only electrons are able to escape the surface of the conducting disk, whereas positrons are immediately reabsorbed at the boundary. In their acceleration along the magnetic field lines, electrons emit curvature photons, which then decay into pairs, triggering the pair cascade. The accelerating electric field E_∥ is then locally screened and the cascade stops. When the plasma flows away from the conducting disk, the vacuum gap reopens and the process restarts.

In Figure 1, we show the electron, positron, and photon densities and E_∥ for a simulation with ℓ_a /R₀ ≃ 0.2 when two pair production bursts are visible in the simulation domain. Due to horizontal gradients of the magnetic field curvature, photon emission and pair production do not occur uniformly across the open field lines. Instead, the cascade is triggered at R ≃ R₀/2 and z ≃ ℓ_a . For R ≳ R₀/2, E_∥ decreases, limiting acceleration and consequently photon emission, whereas for R/R₀ ≪ 1 the radius of curvature of magnetic field lines rapidly diverges and pair production is suppressed (Arons & Scharlemann 1979). For R ≲ R₀/2, the pair cascade develops at z < ℓ_a , i.e., the pair production front is inclined toward the magnetic axis.

As pair production bursts develop, inductive plasma waves are self-consistently excited (Levinson et al. 2005; Cruz et al. 2021). Due to the two-dimensional structure of the pair production front, these waves are not only emitted in the positive z direction, but also in all other directions—see Figure 1(d). The inclination of the wavevectors relative to B is not the result of cross-field particle motion, but rather of field-aligned currents coherently developed on adjacent field lines.

This inclination is a key ingredient in the self-consistent excitation of electromagnetic plasma modes. The wave generation process can be understood as follows (PTS20, Melrose et al. 2020): first, the nonuniformity of E_∥ across dipolar field lines gives rise to an oscillating azimuthal component of the magnetic field ${\tilde{B}}_{\phi }$ via $(\partial \tilde{{\boldsymbol{B}}}/\partial t)\sim -c{\rm{\nabla }}\times {{\boldsymbol{E}}}_{\parallel }$, where $\tilde{{\boldsymbol{B}}}={\tilde{B}}_{\phi }{{\boldsymbol{e}}}_{{\boldsymbol{\phi }}};$ then, an oscillatory electric field component $\tilde{{\boldsymbol{E}}}$ is excited via $(\partial \tilde{{\boldsymbol{E}}}/\partial t)\sim c{\rm{\nabla }}\times \tilde{{\boldsymbol{B}}}$. Gradients in $\tilde{{\boldsymbol{B}}}$ occur predominantly along the normal to the pair production front, which is also the direction of the wavevector k of these waves. In general, both k and $\tilde{{\bf{E}}}$ have components parallel and perpendicular to B _d. The requirement that the angle between the normal to the pair production front and B _d is finite is essential for this process to operate. Hence, modes with k ∥ B _d are never produced without accompanying modes with k ⊥ B _d. All our simulations show that these modes are excited at the local plasma frequency ω₀ (i.e., at the frequency of inductive E _∥ oscillations) and are linearly polarized. As pairs are produced, the local density increases and the wave frequency extends to ∼10–100ω₀. The properties of these modes are consistent with the superluminal O-modes ⁴ identified in previous works (Arons & Barnard 1986, PTS20).

A snapshot of the components of the Poynting vector ${\bf{S}}\,=(c/4\pi )\tilde{{\bf{E}}}\times \tilde{{\bf{B}}}$ is shown in Figure 2. The wave components $\tilde{{\bf{E}}}$ and $\tilde{{\bf{B}}}$ were calculated by subtracting local time-averaged E and B, respectively. Figure 2(a) shows rings of electromagnetic flux emitted around the most recently created plasma burst, centered at z/R₀ ≃ 0.4 and R/R₀ ≃ 0.5. Waves emitted upward carry a positive S_z , whereas waves emitted downward have initially a negative S_z that is then reversed after the waves are reflected at the z = 0 boundary. Electromagnetic waves are emitted in all directions. However, due to the inclination of the bursts relative to the magnetic field lines, a larger Poynting flux is generated on the flanks of the bursts. After some altitude, S_z thus exhibits a double-peaked structure, being larger close to the magnetic axis and the conducting field line boundary and smaller in the center of the open field line bundle.

Figure 2(b) shows that part of the electromagnetic flux rings convert into bands extended in the z direction (see, e.g., R/R₀ ≃ 0.5 and z/R₀ ≃ 1) over time. These bands result from the continuous reflection of the waves between the magnetic axis and the conducting region of closed field lines, and they drift across the open field line region over a timescale ∼R₀/c.

The electrodynamics of the two-dimensional pair cascades identified above for ℓ_a /R₀ ≃ 0.2 holds also for different values of ℓ_a /R₀. However, there are important effects to note: first, the wavelength of the E_∥ oscillations (and consequently of $\tilde{{\boldsymbol{E}}}$ and $\tilde{{\boldsymbol{B}}}$) decreases with decreasing ℓ_a, due to the larger multiplicity generated in the cascades in this regime; second, the pair production burst extends to smaller θ for smaller ℓ_a /R₀; third, the shape of the pair production burst is flatter in the radial direction for smaller ℓ_a /R₀, and more round for larger ℓ_a /R₀—see Figure 3. The latter effect is responsible for a more efficient generation of Poynting flux, since it allows for a larger portion of the pair production bursts to excite waves with k almost purely aligned with z. This is visible in Figure 4, where we show $\langle {S}_{r}\rangle \equiv \langle {S}_{z}\cos \theta +{S}_{R}\sin \theta \rangle$, averaged for all times after the initial vacuum transient at z/R₀ = 2 for simulations performed with different values of ℓ_a /R₀. The angle θ is normalized to ${\theta }_{0}\equiv \arcsin ({R}_{0}/{R}_{* })$, i.e., the analog of the polar cap angle in pulsar magnetospheres. For ℓ_a /R₀ ≃ 0.2, we observe the double-peaked structure in this profile described before (see peaks at θ ≃ 0 and θ ≃ θ₀). However, for large ℓ_a /R₀, the edge component is absent. We also observe that the peak value of 〈S_r 〉 on the magnetic axis decreases with ℓ_a /R₀, but is always ∼10⁻⁶–10⁻⁴ S₀, where ${S}_{0}={{ce}}^{2}{n}_{\mathrm{GJ}}^{2}{R}_{0}^{2}/16\pi$ is the total Poynting flux launched by the rotating conductor, a result consistent with the fraction of pulsar spin-down power observed in the radio (e.g., Lorimer & Kramer 2004).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this Letter, we have presented the first 2D simulations of pulsar polar cap pair cascades including the QED effects from first principles. Our results show that the gap dynamics has two significant two-dimensional features that could not be captured in one-dimensional simulations: (a) pair production is inhibited close to the magnetic axis, due to the null curvature of the magnetic field in this region, and (b) gradients of the magnetic field curvature across the polar cap induce an inclination between the normal to the pair production front and the background magnetic field.

We have observed that E_∥ oscillations inductively driven in the inclined pair production bursts can act as a source for coherent electromagnetic waves. Although the initial oscillations in E_∥ is longitudinal, the induced waves are oblique and electromagnetic in nature, which makes them distinct from the L-mode waves discussed in previous works (Rafat et al. 2019; Melrose et al. 2020). Moreover, the electromagnetic modes observed in our simulations are naturally produced in the oblique pair production fronts, and do not require the production of intermediary modes such as the L modes suggested by Melrose et al. (2020). The induced electromagnetic waves propagate in all directions, but the Poynting flux flows predominantly outward, away from the star. We have identified the coherent electromagnetic modes to be linearly polarized, superluminal O-modes, and verified that the resulting Poynting flux from this emission has two peaks: one on the magnetic axis and another on the edge of the open field line bundle. We interpret this as a consequence of the relative orientation between the normal to the pair production bursts and the background magnetic field. Thus, we expect this mechanism to operate also in magnetic field topologies more complex than a pure dipole (e.g., in multipolar magnetic fields), provided that the field curvature changes over a scale larger than the gap height ∼ℓ_a . In such topologies, the configuration of the open field lines, where cascades develop, may lead to complex-shaped pair production bursts and thus different emission power profiles.

We have also shown that a fraction of the emitted waves is continuously reflected at the last open field line boundary, giving rise to a drifting component in the Poynting flux on a timescale ∼1 μs. We note that this should not be confused with the observed radio drifting subpulses that occur on timescales ∼0.1–1 s (e.g., Weltevrede et al. 2007). The drifts identified in our simulations should be detectable at high time resolution and could be used to diagnose the height of the emission.

The simulations presented in this work adopt a very dilute, space charge-limited flow from a rotating conductor, and focus on the role of pair cascades in providing the current to screen the vacuum gap. In reality, the neutron star is expected to provide a larger current density ∼j_GJ = ρ_GJ c. In that case, the vacuum gaps are induced by the inability of the star to match the current density j > j_GJ, required by the global magnetosphere when general relativistic frame-dragging effects are taken into account (Gralla et al. 2016). It is not expected that pair cascades operate differently when these effects are included; however, given the same B₀ and Ω, the ratio ℓ_a /R₀ may differ slightly from what is presented in this work.

F.C., T.G., and L.O.S. are supported by the European Research Council (ERC-2015-AdG Grant 695088). F.C. is also supported by FCT (Portugal) (grant PD/BD/114307/2016) in the framework of the Advanced Program in Plasma Science and Engineering (APPLAuSE, FCT grant PD/00505/2012). A.C. is supported by NSF grants AST-1806084 and AST-1903335. A.S. is supported by NASA grant 80NSSC18K1099 and NSF grant PHY-1804048. We acknowledge PRACE for granting access to MareNostrum (Barcelona Supercomputing Center, Spain) and to TGCC Joliot Curie (CEA, France), where the simulations presented in this work were performed.