Electric Field Screening in Pair Discharges and Generation of Pulsar Radio Emission

Pulsars are rapidly rotating, highly magnetized neutron stars, most of which exhibit beamed radio emission. Several characteristics of this emission remain unexplained. Notably, the pulsar radio luminosity L_rad has a rough magnitude L_rad ∼ 10⁻⁷–10⁻² L_sd, with L_sd the spin-down luminosity, but has negligible dependence on L_sd (Lorimer & Kramer 2004; Szary et al. 2014). Furthermore, the radio spectrum is S_ω ∼ ω^−1.4±1.0, with ω the radiation frequency and S_ω the intensity emitted at that frequency (Bates et al. 2013).

Compelling observational evidence suggests that in most pulsars, radio emission is produced in the polar regions of the neutron star, where the star's rotation induces a strong electric field that extracts electrons from the neutron star surface and accelerates them to high energies (Sturrock 1971; Ruderman & Sutherland 1975). These electrons stream along curved magnetic field lines, radiating gamma-rays that are absorbed by the pulsar magnetic field, producing electron–positron pairs. The resulting cascade, which includes additional contributions from inverse Compton scattering and synchrotron radiation, leads to the rapid production of a pair plasma in the polar cap region, which screens out the initial inductive electric field (Tsai & Erber 1974; Cheng & Ruderman 1977; Levinson et al. 2005; Timokhin & Harding 2015; Cruz et al. 2021a). Some process occurring in the pair plasma created in this discharge is likely to produce the pulsar's radio emission (Melrose et al. 2021).

Recently, attention has turned toward the possibility that radio emission could be directly produced by the oscillating electric field that is created as the pair discharge screens the initial electric field (Beloborodov 2008; Timokhin 2010; Timokhin & Arons 2013; Melrose et al. 2020; Philippov et al. 2020). Specifically, direct numerical simulations (Timokhin 2010; Timokhin & Arons 2013) of polar cap pair discharges have demonstrated that pair creation is a nonstationary, repeating process in which the accelerating electric field is screened by freshly created pair plasma, which then leaves the polar cap, causing the field to emerge again. Philippov et al. (2020) showed that if pair creation is nonuniform across magnetic field lines, as it will be in any realistic pair cascade, a robust physical mechanism generates coherent radio emission during the screening process. In this Letter, we study the highly relativistic, collisionless, and nonlinear plasma physics governing the evolution of this pair discharge's electric field amplitude and show that this plasma physics may explain the enigmatic pulsar radio luminosity, contribute to the spectrum, and give insight into other aspects of pulsar physics. We verify our analytic models with direct kinetic simulations.

In particular, during an initial nonlinear phase of the discharge, the electric field experiences strong damping that is roughly exponential. This strong damping ceases when the quiver of newly injected particles in the electric field is small. The quantitative statement of the resulting condition on the electric field amplitude yields an expression roughly consistent with the observed magnitude of L_rad. After this point, the discharge experiences a linear phase, which only slightly decreases the emission amplitude but creates a relationship between the emission amplitude and frequency that is consistent with the observed radio spectrum.

This Letter begins in Section 3 by introducing the physical model and equations used to study the discharge plasma. Then, Section 4 describes the simulations of the discharge plasma, which will be used to check the analytical work that forms the bulk of the paper. Section 5 describes the first screening stage of the discharge; Section 6 the following stage of nonlinear, strongly damped oscillations; Section 7 the conditions that lead to the end of this stage's strong damping and set the overall emission amplitude; and Section 8 the following linear stage, which contributes to the pulsar spectrum. Finally, Section 9 shows how our model explains observations of pulsar luminosity and may also give insight into the spectrum and other aspects of pulsar behavior.

The discharge occurring in a pulsar polar cap has a complex and uncertain spatial structure that inhibits analytic study. In this section, we describe the current understanding of the realistic polar discharge structure and relate it to the more simplified plasma physics problem we study in this Letter. This will lead in Section 3 to the statement of the equations studied in this paper.

The actual polar cap discharge proceeds through the cyclic production of clouds of pair plasma (Timokhin 2010; Timokhin & Arons 2013; Timokhin & Harding 2015; Philippov et al. 2020; Cruz et al. 2021b). In this picture, an initially unshielded electric field in the pulsar's polar cap extracts electrons from the neutron star surface and accelerates them to significant Lorentz factors,

$\begin{eqnarray}&&{\gamma }_{b}\sim {10}^{7}.\end{eqnarray} \tag{ 1 }$

Here, we have defined

$\begin{eqnarray}&&\gamma \equiv \displaystyle \frac{1}{\sqrt{1-{\beta }^{2}}},\end{eqnarray} \tag{ 2 }$

with β ≡ v/c (v is the particle speed). These primary electrons radiate gamma-rays, which are eventually absorbed in the pulsar magnetic field, creating lower-energy electron–positron pairs with an average Lorentz factor (Timokhin & Harding 2015)

$\begin{eqnarray}&&{\gamma }_{l}\sim 10-{10}^{3}\end{eqnarray} \tag{ 3 }$

(here, the upper limit is more representative for pairs injected at earlier cascade stages, which are studied here) and a thermal spread T ∼ mc², with m the electron mass (Arendt & Eilek 2002; Timokhin & Arons 2013). Though these pairs are much lower energy than the primary electrons, they move at essentially the same speed, such that the primary and secondary pairs form one cloud of plasma moving away from the star. The secondary pairs interact strongly with the electric field, first shielding it, and then setting up plasma waves caused by the overshoot of the shielding process. The primary population is barely affected by these plasma waves, such that it maintains a constant γ_b ∼ 10⁷ and its ability to create curvature radiation. Thus, there is a nearly constant source of new secondary pairs being added to the oscillations occurring in the pair plasma cloud, and the oscillations damp in response.

Simulations of the two-dimensional polar cap region (Philippov et al. 2020) find that the created plasma waves are superluminal O-mode waves, which consist of an oscillating electric field with components along the background magnetic field and perpendicular to it, and an oscillating magnetic field perpendicular both to the background magnetic field and the oscillating electric field (Arons & Barnard 1986). In this Letter, we are concerned with understanding the physics governing the evolution of the amplitude of these waves.

The physical mechanism responsible for the generation of these waves works in at least two dimensions, requiring the nonuniformity of pair creation across magnetic field lines. However, the electric field driving the waves is the oscillating parallel electric field, which interacts directly with the pair plasma through a plasma current. This parallel electric field is determined by the dynamics of the discharge, i.e., by how the initial accelerating electric field is screened by injected pairs. Hence, wave evolution (at least in its initial stages) is determined by the evolution of the parallel electric field, which we can consider to be a one-dimensional problem. The waves observed in Philippov et al. (2020) have small spatial frequencies k, so to allow analytic progress, we can also consider only the temporal evolution of the electric field, i.e., we can ignore spatial variation along the magnetic field lines in both the initial electric field and in the particle injection rate. ⁵

Thus, in this Letter, we model the discharge plasma physics with a one-dimensional setup in which an initially strong electric field lies parallel to a magnetic field. Low-energy pairs with γ_l are added to this system uniformly in space at a chosen rate, representing the pairs created in the actual plasma cloud (which, unlike our simplified one-dimensional setup, has curved magnetic fields) by the higher-energy pairs. In response, the electric field is shielded and begins to oscillate with spatial frequency k = 0 (higher-k modes are not initially seeded because, unlike in Philippov et al. 2020, the setup does not have spatial structure).

The physics governing the evolution of this oscillating electric field is expected to be roughly equivalent to the physics governing the low-k O-mode oscillations observed in Philippov et al. (2020) and presumably also present in actual pulsar polar caps. (Indeed, as we discuss in Section 7, simulations of such an initially k = 0 setup eventually fragment into higher-k modes similar to those seen in two-dimensional systems.) However, a k = 0 setup is more tractable analytically and still gives significant insight into the physics of low-k electromagnetic modes controlled by their parallel dynamics.

In this section, we present the specific equations that we use to understand polar cap physics through the picture introduced previously.

The pulsar polar cap region can be represented by a uniform straight magnetic field B_⋆ and a uniform parallel initial electric field E₀, corresponding to the magnetic field in the pulsar polar region and the unshielded electric field created there by the pulsar's rotation. These fields initially exist in a vacuum. For an aligned rotator, the electric field E has initial strength (Deutsch 1955; Goldreich & Julian 1969)

$\begin{eqnarray}&&E\left(t=0\right)\equiv {E}_{0}=-\displaystyle \frac{{\rm{\Omega }}{R}_{\mathrm{pc}}}{c}{B}_{\star },\end{eqnarray} \tag{ 4 }$

with t the time, Ω the pulsar angular velocity, c the speed of light, and R_pc the polar cap radius, defined by

$\begin{eqnarray}&&{R}_{\mathrm{pc}}={R}_{\star }\sqrt{\displaystyle \frac{{R}_{\star }{\rm{\Omega }}}{c}},\end{eqnarray} \tag{ 5 }$

where R_⋆ is the neutron star radius, R_⋆ ≈ 1.2 × 10⁶ cm.

In an actual pulsar, the initial electric field E₀ extracts electrons from the neutron star into the polar cap vacuum and accelerates them to significant Lorentz factors (1), seeding a continual cascade of lower-energy pairs of average Lorentz factors (3). The high Lorentz factor γ_b strongly suppresses the interaction of the high-energy population with any plasma waves, so we ignore this population. Instead, we consider how Equation (4) evolves as lower-energy pairs at γ_l are added to the initial vacuum. The lower-energy population will damp the electric field; the coupled evolution of the electric field and the lower-energy pair population can be described by the relativistic Vlasov–Maxwell equations.

Specifically, the evolution of the electric field is governed by Ampère's law, which for a curl-free magnetic field reads

$\begin{eqnarray}&&{\partial }_{t}E+4\pi j=0,\end{eqnarray} \tag{ 6 }$

with j the plasma current. The plasma consists of two species s = + , −: electrons (with charge q₋ = −e) and positrons (with charge q₊ = +e). Both have mass m. Furthermore, we define

$\begin{eqnarray}&&u\equiv \gamma \beta \end{eqnarray} \tag{ 7 }$

as the particle momentum, equal to its component along the magnetic field. ⁶ The electrons and positrons are each described by their corresponding distributions,

$\begin{eqnarray}&&{f}_{s}(u,t),\end{eqnarray} \tag{ 8 }$

so that the plasma current in Equation (6) is given by

$\begin{eqnarray}&&j=c\sum _{s}{q}_{s}{\int }_{-\infty }^{\infty }\beta {f}_{s}\left(u,t\right){du},\end{eqnarray} \tag{ 9 }$

which couples Ampère's law to the Vlasov equation.

The evolution of ${f}_{s}\left(u,t\right)$ is given by the Vlasov equation with a source term representing the injection of lower-energy pairs:

$\begin{eqnarray}&&{\partial }_{t}{f}_{s}+\displaystyle \frac{{q}_{s}E\left(t\right)}{{mc}}{\partial }_{u}{f}_{s}=\displaystyle \frac{{ \mathcal S }\left(u-{u}_{l}\right)}{2},\end{eqnarray} \tag{ 10 }$

Here, u_l is the mean momentum of the newly injected lower-energy pairs,

$\begin{eqnarray}&&{u}_{l}\equiv {\gamma }_{l}{\beta }_{l},\end{eqnarray} \tag{ 11 }$

with β_l the value of β corresponding to γ_l . The form of ${ \mathcal S }$ characterizes the thermal spread of the newly injected pairs about their mean injection momentum; its integral gives the rate of increase of species density n_s due to the addition of new pairs, i.e., ${\partial }_{t}{n}_{s}={\int }_{-\infty }^{\infty }\left[{ \mathcal S }\left(u-{u}_{l}\right)/2\right]{du}$. For our analytical arguments, we will take a simplified source model with all pairs injected at one velocity, i.e.,

$\begin{eqnarray}&&{\partial }_{t}{f}_{s}+\displaystyle \frac{{q}_{s}E\left(t\right)}{{mc}}{\partial }_{u}{f}_{s}=\displaystyle \frac{S\delta \left(u-{u}_{l}\right)}{2},\end{eqnarray} \tag{ 12 }$

with S/2 giving the rate per species of pair creation. (The effect of thermal spread is considered in Appendix C.) For simplicity, we assume this rate to be constant in time.

Next, we apply a set of normalizations to Equations (6), (9), and (12). Let us define t₀ as the time it takes to change a pair's momentum by an amount Δu = 1 in the initial electric field, i.e.,

$\begin{eqnarray}&&{t}_{0}\equiv \displaystyle \frac{{mc}}{\left|{E}_{0}\right|e},\end{eqnarray} \tag{ 13 }$

with $\left|{E}_{0}\right|$ the absolute value of the initial electric field. Then, we normalize time to this value,

$\begin{eqnarray}&&\hat{t}\equiv \displaystyle \frac{t}{{t}_{0}},\end{eqnarray} \tag{ 14 }$

f_s to the density injected in that time,

$\begin{eqnarray}&&{\hat{f}}_{s}\equiv \displaystyle \frac{{f}_{s}}{{{St}}_{0}/2},\end{eqnarray} \tag{ 15 }$

and E to its initial strength,

$\begin{eqnarray}&&\hat{E}\equiv \displaystyle \frac{E}{\left|{E}_{0}\right|}.\end{eqnarray} \tag{ 16 }$

(Later, when we wish to refer to the amplitude of an oscillating electric field, we will use $\overline{\hat{E}}$.) Then, Equation (12) reduces to

$\begin{eqnarray}&&{\partial }_{\hat{t}}{\hat{f}}_{s}+\displaystyle \frac{{q}_{s}}{e}\hat{E}\left(\hat{t}\right){\partial }_{u}{\hat{f}}_{s}=\delta \left(u-{u}_{l}\right)\end{eqnarray} \tag{ 17 }$

and Equation (6) to

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}\left(\hat{t}\right)+\hat{j}=0,\end{eqnarray} \tag{ 18 }$

where Equation (9) becomes

$\begin{eqnarray}&&\hat{j}\equiv \displaystyle \frac{1}{\xi }\sum _{s}\displaystyle \frac{{q}_{s}}{e}{\int }_{-\infty }^{\infty }\beta {\hat{f}}_{s}{du}.\end{eqnarray} \tag{ 19 }$

Here, we have defined a parameter of key importance in determining the discharge dynamics:

$\begin{eqnarray}&&\xi \equiv \displaystyle \frac{2\left({E}_{0}^{2}/8\pi \right)}{{{mc}}^{2}\left({{St}}_{0}/2\right)},\end{eqnarray} \tag{ 20 }$

which relates the energy in the initial pulsar electric field to the rest-mass energy in the pairs injected during the time interval t₀. To find the appropriate value of this parameter to describe a pulsar, we introduce the Goldreich–Julian density (Goldreich & Julian 1969), ${n}_{\mathrm{GJ}}\equiv \left({\rm{\Omega }}{B}_{\star }\right)/\left(2\pi {ce}\right)$, and state that

$\begin{eqnarray}&&S\sim \displaystyle \frac{\lambda {n}_{\mathrm{GJ}}}{{R}_{\star }/c},\end{eqnarray} \tag{ 21 }$

with λ the multiplicity. Then, the value of ξ in a pulsar with magnetic field B_⋆, multiplicity λ, and period P is

$\begin{eqnarray}&&\xi \sim {10}^{12}\ {\left(\displaystyle \frac{{B}_{\star }}{{10}^{12}\,{\rm{G}}}\right)}^{2}\ \left(\displaystyle \frac{{10}^{5}}{\lambda }\right)\ {\left(\displaystyle \frac{1\,{\rm{s}}}{P}\right)}^{7/2},\end{eqnarray} \tag{ 22 }$

a large value.

To analyze the nonlinear and linear stages of the electric field oscillations, it will be useful to combine Equations (17), (18), and (19) to obtain another relationship for the electric field. Noting that the current, represented by the second term in Equation (18), receives no contribution from particles newly injected at u_l , we can trivially use Equation (17), without reference to the source term, to find the exact, nonlinear relationship,

$\begin{eqnarray}&&{\partial }_{\hat{t}}^{2}\hat{E}+{\hat{\omega }}^{2}\left(\hat{t}\right)\hat{E}\left(\hat{t}\right)=0,\end{eqnarray} \tag{ 23 }$

where

$\begin{eqnarray}&&{\hat{\omega }}^{2}\left(\hat{t}\right)\equiv \displaystyle \frac{{\hat{n}}_{+}}{\xi }{\left\langle \displaystyle \frac{1}{{\gamma }^{3}}\right\rangle }_{+}+\displaystyle \frac{{\hat{n}}_{-}}{\xi }{\left\langle \displaystyle \frac{1}{{\gamma }^{3}}\right\rangle }_{-};\end{eqnarray} \tag{ 24 }$

$\hat{\omega }$ is the normalized relativistic plasma frequency characterizing k = 0 oscillations (which may be nonlinear or linear). ⁷ The normalized density ${\hat{n}}_{s}$ is defined by

$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{du}{\hat{f}}_{s}\equiv {\hat{n}}_{s}=\hat{t}\end{eqnarray} \tag{ 25 }$

and the expected values ${\left\langle 1/{\gamma }^{3}\right\rangle }_{+}$, ${\left\langle 1/{\gamma }^{3}\right\rangle }_{-}$ by

$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }\displaystyle \frac{{du}}{{\gamma }^{3}}{\hat{f}}_{s}\equiv {\hat{n}}_{s}{\left\langle \displaystyle \frac{1}{{\gamma }^{3}}\right\rangle }_{s}.\end{eqnarray} \tag{ 26 }$

In deriving Equation (23), we have used the relation

$\begin{eqnarray}&&d\beta =\displaystyle \frac{{du}}{{\gamma }^{3}}\end{eqnarray} \tag{ 27 }$

and integrated by parts. The rest of the paper will consider the solutions to Equations (17), (18), and (23), which describe the damping of the electric field in the polar cap by the newly injected pairs and thus determine the radio emission characteristics. We emphasize that the relationships in this section are derived without linearization. They describe the oscillation of the plasma throughout all stages of the discharge, including the earliest stages, when the pair production rate is faster than the frequency $\hat{\omega }$.

The goal of the paper is the analytical study of the equations introduced previously, leading to relationships that can be used to characterize the radio emission produced in a pulsar polar cap discharge. To confirm these relationships, we compare them to kinetic particle-in-cell simulations with TRISTAN-MP v2 (Spitkovsky 2005; Hakobyan & Spitkovsky 2020).

We initialize a one-dimensional box with a uniform electric field and no plasma. At each subsequent time step, pair plasma of momentum u_l = 10 (at the lower end of the range described by Equation (3)) and small temperature T = 0.1 mc² (the effect of higher temperatures is considered in Appendix C) is injected into the box. The initial electric field amplitude and injection rate determine the value of ξ; as expected, the appropriately normalized simulation results depend only on the value of ξ and not on the electric field amplitude or source rate separately. (For convenience, after an initial test, only the electric field strength was varied to vary ξ.) We run a set of simulations with values of ξ ranging from ξ = 2.8 × 10⁶ to ξ = 4.4 × 10⁷.

After two time steps, there are 10 skin depths in the simulation box, with each resolved by 700 cells. The simulation is continued, with density increasing and ${\left\langle 1/{\gamma }^{3}\right\rangle }_{+,-}$ evolving, until the discharge is well past the transition to linearity; in all runs, the skin depth is resolved by at least 35 cells at the end of the simulation.

As pair creation begins in an unshielded electric field, the first event that occurs is the screening of the initial electric field, a process that has been considered previously in Levinson et al. (2005) and Cruz et al. (2021b). In this section, we examine the physics of this stage, leading to an expression for screening time, Equation (30), which agrees with previous results in the literature, and a new expression for the pair distribution function, Equation (31).

The strong electric field initially present in the polar cap rapidly accelerates the first newly injected lower-energy pairs to large speeds, such that β → 1 and in Equation (18) $\hat{j}$, defined in Equation (19), approaches

$\begin{eqnarray}&&\hat{j}=\displaystyle \frac{2}{\xi }{\int }_{-\infty }^{\infty }{\hat{f}}_{+}\left(u,\hat{t}\right){du}=\displaystyle \frac{2\hat{t}}{\xi },\end{eqnarray} \tag{ 28 }$

where in the second equality we have applied the definition in Equation (25). Integration of Equation (18) in time then gives

$\begin{eqnarray}&&\hat{E}\left(\hat{t}\right)=-1+\displaystyle \frac{{\hat{t}}^{2}}{\xi }.\end{eqnarray} \tag{ 29 }$

The initial electric field is fully screened when $\hat{E}=0$, which occurs at time

$\begin{eqnarray}&&{\hat{t}}_{\mathrm{screen}}={\xi }^{1/2}.\end{eqnarray} \tag{ 30 }$

As the electric field screens, particles are accelerated. Solving Equation (17) using Equation (29) shows the electron distribution function during the screening period is given by

$\begin{eqnarray}&&\begin{array}{l}{\hat{f}}_{-}(u,\hat{t})\\ =\ \left\{\begin{array}{ll}0, & u-{u}_{l}\lt 0\\ \displaystyle \frac{1}{1-{\hat{\tau }}_{\mathrm{source}}^{2}\left(u-{u}_{l},\hat{t}\right)/\xi }, & 0\lt u-{u}_{l}\lt \hat{t}-\displaystyle \frac{{\hat{t}}^{3}}{3\xi }\\ 0, & u-{u}_{l}\gt \hat{t}-\displaystyle \frac{{\hat{t}}^{3}}{3\xi },\end{array}\right.\end{array}\end{eqnarray} \tag{ 31 }$

where the parameter ${\hat{\tau }}_{\mathrm{source}}\left(u-{u}_{l},t\right)$ is defined by the cubic equation

$\begin{eqnarray}&&u-{u}_{l}=\left(\hat{t}-{\hat{\tau }}_{\mathrm{source}}\right)-\displaystyle \frac{1}{3\xi }\left({\hat{t}}^{3}-{\hat{\tau }}_{\mathrm{source}}^{3}\right).\end{eqnarray} \tag{ 32 }$

The positron distribution function has the same form, but consists of negative values of u − u_l . We plot the analytic electron distribution function at the time of screening ${\hat{t}}_{\mathrm{screen}}$ with the distribution function from a direct numerical simulation in Figure 1, finding good agreement. The largest disagreement with our analytic prediction occurs at low u − u_l , which corresponds to particles injected just before ${\hat{t}}_{\mathrm{screen}}$, when the electric field is weak and Equation (28) begins to break down. The maximum u obtained by the secondary pairs during the discharge is given by the maximum u at which Equation (36) has support at ${\hat{t}}_{\mathrm{screen}}$, i.e., $u-{u}_{l}=2\sqrt{\xi }/3$.

Zoom In Zoom Out Reset image size

At the time of screening, the electric field reaches zero, but the velocities obtained by the electrons and positrons continue to separate them, reestablishing an electric field in the plasma and beginning the oscillations capable of creating radio emission, which we study next, beginning in the next section with the study of the first, large oscillations.

The initial oscillations in the plasma following screening are, for large values of ξ, highly nonlinear and occur in a medium that changes quickly in comparison to the speed of the oscillation. As these oscillations proceed, they are damped by the continual creation of pairs, which are spun up into the wave oscillation, removing energy from the wave electric field. This section examines wave damping during this period.

The electric field amplitude in this stage evolves according to Equation (23), which was derived without making any assumption of linearity. The behavior of the electric field is thus specified by the behavior of the frequency ${\hat{\omega }}^{2}$, defined in Equation (24). Consideration of Equation (24) shows that ${\hat{\omega }}^{2}$ depends on density ${\hat{n}}_{s}$, which grows linearly as new pairs are added to the system, and on the values of ${\left\langle 1/{\gamma }^{3}\right\rangle }_{s}$, which evolve both due to the addition of new pairs at u = u_l to the system and due to the large displacements in particle momentum that occur in response to a nonlinear wave.

This evolution of ${\hat{\omega }}^{2}$ can be seen in the results of a simulation of the nonlinear stage, shown in the first panel of Figure 2. Here, ${\hat{\omega }}^{2}$ has large spikes, which represent times when the second derivative of the electric field increases and decreases dramatically in a short period of time. These spikes are caused by particles added at times when the electric field goes through zero. These newly added pairs are not initially accelerated, so they contribute to a build up of a large number of particles near u = u_l . Slightly later, one sign of pairs participating in this build up is dragged by the electric field through u = 0, where the built-up pairs contribute to a dramatic spike in $\left\langle 1/{\gamma }^{3}\right\rangle$. Afterwards, as the particles are accelerated to high u in the opposite direction, ${\hat{\omega }}^{2}$ decreases rapidly as pairs are accelerated to high γ. Crucially, the peak of the spike occurs slightly after the zero of the electric field.

Zoom In Zoom Out Reset image size

In Equation (23), a frequency described by this shifted spikes structure corresponds to efficient damping of the electric field. Let us first consider an elementary demonstration of this effect; later, we will build a more involved model of the damping for the specific case of the polar discharge.

Define time $\hat{t}={\hat{T}}_{h}$ as a point of zero-crossing for $\hat{E}$ that occurs at the end of the half-period h of the oscillation. We endeavor to determine how a spiked frequency such as that observed in our simulation will change the derivative and amplitude of $\hat{E}$ between period h = k and h = k + 1. Let us consider a simplified form of ${\hat{\omega }}^{2}$ containing only delta function spikes, offset by an amount Ξ from $\hat{t}={\hat{T}}_{k}$:

$\begin{eqnarray}&&{\hat{\omega }}^{2}=\delta \left(\hat{t}-{\hat{T}}_{k}-{\rm{\Xi }}\right).\end{eqnarray} \tag{ 33 }$

Plugging this expression into Equation (23) and integrating the resulting equation over a small region $\hat{t}-{\hat{T}}_{k}\in \left[-\kappa ,\kappa \right]$, where κ > Ξ, gives

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}+\kappa \right)-{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}-\kappa \right)=-\hat{E}\left({\hat{T}}_{k}+{\rm{\Xi }}\right).\end{eqnarray} \tag{ 34 }$

This equation represents a decrease in the magnitude of ${\partial }_{\hat{t}}\hat{E}$ between successive half-periods h = k and h = k + 1 of an electric field oscillation, causing the maximum amplitude of the electric field in the later half-period to be smaller, constituting damping of the electric field. The physical mechanism for this damping is the extraction of energy from the electric field to decelerate pairs in a way that causes a spike in the frequency.

This mechanism is at play in the nonlinear stage of the polar cap discharge. To obtain more specific results, we now build a more involved model; this model includes a study of the specific evolution of ${\hat{f}}_{-}$, yielding a more specific form for ${\hat{\omega }}^{2}$ than Equation (33) and, via Equation (23), a functional form for electric field damping, which we will present in Equation (48). The treatment is imprecise due to the nonlinear physics involved.

To develop a model for ${\hat{f}}_{-}$, let us call the amplitude of the electric field during the half-period h ${\overline{\hat{E}}}_{h};$ also, we consider the change in electric field amplitude and in $\hat{\omega }$ away from the zeros in $\hat{E}$ to be negligible. The electric field is approximately sinusoidal, ⁸ so that in the half-period h Equation (17) becomes

$\begin{eqnarray}&&{\partial }_{\hat{t}}{\hat{f}}_{-}+\left[{\overline{\hat{E}}}_{h}\cos \left(\hat{\omega }\hat{t}\right)\right]{\partial }_{u}{\hat{f}}_{-}=\delta \left(u-{u}_{l}\right).\end{eqnarray} \tag{ 35 }$

Solving Equation (35) gives

$\begin{eqnarray}&&{\hat{f}}_{-}={\hat{f}}_{-,\mathrm{new}}\left(u-{u}_{l},\hat{t}\right)+{\hat{f}}_{-,\mathrm{old}}\left[u\hat{\omega }-{\overline{\hat{E}}}_{h}\sin \left(\hat{\omega }\hat{t}\right)\right].\end{eqnarray} \tag{ 36 }$

Here, ${\hat{f}}_{-,\mathrm{old}}$, corresponding to the homogeneous solution of Equation (35), represents the evolution of particles already present at the beginning of the half-period h of the oscillation. The term ${\hat{f}}_{-,\mathrm{new}}\left(u,\hat{t}\right)$ gives the evolution of newly added particles and is given, in a half-period in which electrons are decelerated, by

$\begin{eqnarray}&&\begin{array}{l}{\hat{f}}_{-,\mathrm{new}}(u-{u}_{l},\hat{t})\\ =\,\left\{\begin{array}{ll}0, & u-{u}_{l}\gt {u}_{h}\left(\hat{t}\right)\\ \displaystyle \frac{1/{\overline{\hat{E}}}_{h}}{\sqrt{{g}_{+}\left(u-{u}_{l},\hat{t}\right){g}_{-}\left(u-{u}_{l},\hat{t}\right)}}, & {u}_{h}\left(\hat{t}\right)\lt u-{u}_{l}\lt 0\\ 0, & u-{u}_{l}\gt 0,\end{array}\right.\end{array}\end{eqnarray} \tag{ 37 }$

where

$\begin{eqnarray}&&{g}_{+}\left(u-{u}_{l},\hat{t}\right)=1+\displaystyle \frac{\left(u-{u}_{l}\right)\hat{\omega }}{{\overline{\hat{E}}}_{h}}-\sin \left(\hat{\omega }\hat{t}\right),\end{eqnarray} \tag{ 38 }$

$\begin{eqnarray}&&{g}_{-}\left(u-{u}_{l},\hat{t}\right)=1-\displaystyle \frac{\left(u-{u}_{l}\right)\hat{\omega }}{{\overline{\hat{E}}}_{h}}+\sin \left(\hat{\omega }\hat{t}\right),\end{eqnarray} \tag{ 39 }$

and

$\begin{eqnarray}&&{u}_{h}\left(\hat{t}\right)\equiv \displaystyle \frac{{\overline{\hat{E}}}_{h}\sin \left(\hat{\omega }\hat{t}\right)-{\overline{\hat{E}}}_{h}}{\hat{\omega }}\leqslant 0.\end{eqnarray} \tag{ 40 }$

The second panel of Figure 2 shows a comparison between Equation (37) and the distribution function at the second zero-crossing in the simulation illustrated in the first panel of Figure 2, indicating rough agreement. For convenience, the phase space in the simulation at this time is presented in the third panel of Figure 2.

Let us consider how this distribution function changes the electric field between the half-period h = k and h = k + 1; then, we will spread this change across a whole half-period to obtain a continuous expression for the damping.

Consideration of Equation (36) reveals that the frequency, Equation (24), represented in our elementary demonstration by Equation (33), is better represented near the transition between the half-period h = k and h = k + 1 by

$\begin{eqnarray}&&{\hat{\omega }}^{2}\left(\hat{t}\right)=\displaystyle \frac{1}{\xi }\left[g\left(\hat{t}\right)+A\delta \left(\displaystyle \frac{\hat{t}-{\hat{T}}_{k}-{\rm{\Xi }}}{{\hat{t}}_{b}}\right)\right].\end{eqnarray} \tag{ 41 }$

(Recall that ${\hat{T}}_{k}$ occurs at the end of half-period k.) Here, g reflects the slow change in the frequency due to the evolving nature of ${\hat{f}}_{-,\mathrm{old}}$ from one half-period to the next as the particles in each ${\hat{f}}_{-,\mathrm{new}}$ enter the old particles in the next half-period. This evolution includes a linear-in-time addition of density to the discharge and a slow increase in the overall value of $\left\langle 1/{\gamma }^{3}\right\rangle$. The second term models the spikes in $\left\langle 1/{\gamma }^{3}\right\rangle$ caused by ${\hat{f}}_{-,\mathrm{new}}$, which have characteristic height A and width ${\hat{t}}_{b}$. In Appendix B, we roughly estimate the scaling of these quantities with ξ:

$\begin{eqnarray}&&A\sim {\xi }^{0.2}\end{eqnarray} \tag{ 42 }$

and

$\begin{eqnarray}&&{\hat{t}}_{b}\sim {\xi }^{0.3}.\end{eqnarray} \tag{ 43 }$

Most damping is due to the delta function in Equation (41) via the physics described in our elementary demonstration. As in Equation (34), we integrate over a small region $\hat{t}-{\hat{T}}_{k}\,\in \left[-\kappa ,\kappa \right]$ to state (ignoring factors of order unity)

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}+\kappa \right)-{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}-\kappa \right)=-{\xi }^{-0.5}\hat{E}\left({\hat{T}}_{k}+{\rm{\Xi }}\right).\end{eqnarray} \tag{ 44 }$

The left side of the previous expression can be restated in terms of the frequency and the change in electric field amplitude:

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}+\kappa \right)-{\partial }_{\hat{t}}\hat{E}\left({\hat{T}}_{k}-\kappa \right)=\hat{\omega }{\overline{\hat{E}}}_{k+1}-\hat{\omega }{\overline{\hat{E}}}_{k}.\end{eqnarray} \tag{ 45 }$

(Recall that ${\overline{\hat{E}}}_{k+1}$ represents the wave amplitude in the next half-period of the oscillation.) Then, by representing the total change in amplitude as a rate multiplied by a time, we also state

$\begin{eqnarray}&&\hat{\omega }{\overline{\hat{E}}}_{k+1}-\hat{\omega }{\overline{\hat{E}}}_{k}=\hat{\omega }\left[\displaystyle \frac{{\rm{\Delta }}{\overline{\hat{E}}}_{h}}{{\rm{\Delta }}{\hat{t}}_{h}}{\rm{\Delta }}{\hat{t}}_{h}\right],\end{eqnarray} \tag{ 46 }$

with ${\rm{\Delta }}{\overline{\hat{E}}}_{h}$ representing the change in electric field amplitude from one half-period to the next, and ${\rm{\Delta }}{\hat{t}}_{h}$ the length of a half-period. In Equation (46), we can make the replacement ${\rm{\Delta }}{\overline{\hat{E}}}_{h}/{\rm{\Delta }}{\hat{t}}_{h}\to d\overline{\hat{E}}/d\hat{t}$, with $d\overline{\hat{E}}/d\hat{t}$ representing a continuous change in the amplitude of the electric field oscillations; furthermore, we note that ${\rm{\Delta }}{\hat{t}}_{h}\sim 1/\hat{\omega }$ (again ignoring order unity factors). The continuous analog of Equation (44), that is, an expression that spreads the localized damping caused by a spiked frequency across time, is then

$\begin{eqnarray}&&\displaystyle \frac{d\overline{\hat{E}}}{d\hat{t}}=-{\xi }^{-0.5}\overline{\hat{E}}.\end{eqnarray} \tag{ 47 }$

The evolution of the electric field amplitude during the nonlinear phase is then described by

$\begin{eqnarray}&&\overline{\hat{E}}\left(\hat{t}\right)={e}^{-\left({\xi }^{-0.5}\right)\hat{t}},\end{eqnarray} \tag{ 48 }$

where the initial electric field amplitude is $\overline{\hat{E}}=1$. The initial period of electric field evolution during the polar cap discharge is thus marked by strong, exponential damping of the amplitude of the electric field oscillations. This is caused by particles added during zero-crossings of the electric field having their momenta reversed and accelerated to high values in the opposite direction, sapping energy from the electric field. This damping grows weaker with higher values of ξ, which correspond to weaker pair damping.

To verify this model, we use our set of simulations with several values of ξ to fit an exponential decay $\overline{\hat{E}}\left(\hat{t}\right)={e}^{-a\hat{t}}$ to the peaks of the first three peaks of the electric field amplitude. An example of this fit is shown in the fourth panel of Figure 2. Then, we assemble the set of exponential decay fit constants a and fit these constants to a power-law function of ξ. The resulting fit is

$\begin{eqnarray}&&\overline{\hat{E}}\left(\hat{t}\right)={e}^{-\left(0.4{\xi }^{-0.5}\right)\hat{t}},\end{eqnarray} \tag{ 49 }$

in agreement with Equation (48). The fit for a compared with the values of a measured from simulation data can be seen in the final panel of Figure 2.

As the discharge proceeds, the initial stage of slow, large-amplitude oscillations, which experience the exponential damping described by Equation (48), will cease, and the discharge will proceed to a phase of smaller oscillations at a higher frequency. We call this later regime the "linear" stage, which we analyze in Section 8, revealing significantly weaker damping than in the nonlinear stage. ⁹ In this section, we discuss the conditions that cause the transition between the nonlinear stage and the linear stage.

In the nonlinear stage, damping is caused by sharp spikes in ${\hat{\omega }}^{2}$. These strong spikes will cease to exist when, in Equation (23),

$\begin{eqnarray}&&\displaystyle \frac{{\partial }_{\hat{t}}{\hat{\omega }}^{2}}{{\hat{\omega }}^{3}}\sim 1,\end{eqnarray} \tag{ 50 }$

where ${\hat{\omega }}^{2}$, defined in Equation (24), varies in time like Equation (41). The time variation in $\hat{\omega }$ comes from two sources: the change in the values of ${\hat{n}}_{s}$ and in the values of ${\left\langle 1/{\gamma }^{3}\right\rangle }_{s}$. The component of ${\partial }_{\hat{t}}{\hat{\omega }}^{2}$ coming from the change in ${\hat{n}}_{s}$ is already small at the time of screening and does not contribute to the spiking behavior. Instead, the spiking results from large displacements in the particle momentum as the wave electric field grows from zero to its full amplitude over the course of a half-period of oscillation, changing the values of ${\left\langle 1/{\gamma }^{3}\right\rangle }_{s}$.

These displacements become negligible roughly at the point where the displacement in particle u caused by the wave is insufficient to reverse the typical u of the injected particles, u_l . The variation in a particle's momentum, δ u, is found by considering the second term in Equation (36) or the maximum magnitude u − u_l with support in Equation (37),

$\begin{eqnarray}&&\delta u\sim \displaystyle \frac{\overline{\hat{E}}}{\hat{\omega }}.\end{eqnarray} \tag{ 51 }$

This quantity is a relativistic version of the quiver or jitter velocity familiar from classical plasma physics (Catto & Myra 1989). For this quantity to no longer be large enough to reverse u_l requires

$\begin{eqnarray}&&\displaystyle \frac{\overline{\hat{E}}}{\hat{\omega }{u}_{l}}\sim 1.\end{eqnarray} \tag{ 52 }$

After this point, spiking activity and the concomitant strong exponential damping considered in Section 6 cease. We will show that the subsequent discharge evolution is characterized by much weaker damping, such that the final radio emission amplitude is roughly determined by the electric field amplitude at the time Equations (50) and (52) are satisfied.

Additional simplification is possible when conditions stronger than Equations (50) and (52) hold, specifically,

$\begin{eqnarray}&&\displaystyle \frac{{\partial }_{\hat{t}}{\hat{\omega }}^{2}}{{\hat{\omega }}^{3}}\ll 1\end{eqnarray} \tag{ 53 }$

and

$\begin{eqnarray}&&\displaystyle \frac{\overline{\hat{E}}}{\hat{\omega }{u}_{l}}\ll 1,\end{eqnarray} \tag{ 54 }$

which state, respectively, that the timescale of background plasma variation is significantly longer than the period of oscillation and that the displacement in particle momentum caused by the wave is very small.

We can observe the onset of Equations (50) and (52), and a later progression to Equations (53) and (54), in our simulation. Figure 3 displays the values of ${\hat{\omega }}^{2}$, the rms value of the electric field in the simulation box divided by the frequency, $\sqrt{\left\langle {\hat{E}}^{2}\right\rangle }/\hat{\omega }$, and ${\partial }_{\hat{t}}{\hat{\omega }}^{2}/{\hat{\omega }}^{3}$ in two simulations. Considering the simulation with ξ = 2.8 × 10⁶ (a higher value of ξ, which displays similar trends, is shown in the right panel of Figure 3), we observe an initial stage where $\overline{\hat{E}}/\left(\hat{\omega }{u}_{l}\right)\gg 1$, the plasma frequency changes quickly, and the electric field damps rapidly. Once the simulation obeys Equation (52), around $\hat{t}=6000$, the strong damping ceases. The frequency still has significant variation, resulting from coherent displacements in the plasma distribution function caused by an electric field that, though not able to reverse newly injected particles, is able to cause noticeable momentum changes. This means that although Equation (50) roughly holds, Equation (53) is not satisfied. This situation is difficult to understand analytically but corresponds in our simulation to negligible damping.

Zoom In Zoom Out Reset image size

Later, around $\hat{t}=11,000$, our simulation of k = 0 electric field oscillations fragments into superluminal finite-k modes, which may be a result of a parametric-type instability like that considered in Cruz et al. (2021a). An illustration of this fragmentation is given in Figure 4. After this point, the value of ${\hat{\omega }}^{2}$, where here ${\hat{\omega }}^{2}$ is given by Equation (24) evaluated across the simulation box, which is proportional to, but not the same as, the frequency of a given k mode, does change slowly. This evolution occurs because the electric field is no longer able to cause coherent changes in the plasma distribution function. Conditions (53) and (54) hold, and a new method of analysis, to be explained in Section 8, becomes possible.

Zoom In Zoom Out Reset image size

We note here that we also ran simulations beginning with a finite-k mode. Specifically, we ran simulations starting from a sinusoidal electric field with one period in the simulation box (equivalent to a wavenumber k = 0.63 ${d}_{e}^{-1}$ with skin depth calculated from the plasma density two time steps into the simulation); the maximum initial electric field in these simulations is $\sqrt{2}$ larger than the uniform electric field in the normal simulations, such that the rms electric field amplitude is the same. Plasma is injected in the same way and at the same rate as in the k = 0 simulations. The discharge in the nonzero-k simulation quickly evolves into even higher-k oscillations because the initial field is more quickly screened near the zeros of the electric field. There is no later fragmentation, and the oscillations remain superluminal throughout the discharge.

The resulting damping of these simulations is shown in Figure 3 by the orange dashed line. The nonzero-k discharge exhibits a strong exponential damping phase, which ends when the rms amplitude of the electric field in the box approximately satisfies Equation (52). This damping is followed by a weak damping phase, without the intermediate stage of minimal damping and $\overline{\hat{E}}/\left({u}_{l}\hat{\omega }\right)\sim 1$ exhibited by the k = 0 simulations. However, the value of Equation (24) evaluated across the simulation box in finite-k simulations does not show the spiking behavior considered in Section 6, because the picture we advanced in that section will hold only locally.

In a plasma where Equations (53) and (54) are satisfied, the system enters a period in which the plasma wave is a small perturbation on the background and the frequency of the oscillation changes slowly, as quantified by Equation (53). In our simulations, this occurs after the fragmentation of the k = 0 mode into finite-k oscillations; however, because the observed oscillations still have low values of k, the evolution of their amplitude can be roughly analyzed in the same way as a k = 0 mode.

Specifically, the electric field differential equation, Equation (23), can be analyzed directly using textbook WKB techniques (Stix 1992), giving

$\begin{eqnarray}&&{\overline{\hat{E}}}^{2}\sim \displaystyle \frac{1}{\hat{\omega }},\end{eqnarray} \tag{ 55 }$

with the electric field oscillation frequency determined by the current value of Equation (24). In this stage, ${\hat{\omega }}^{2}$ changes due to increasing plasma density and evolving ${\left\langle 1/{\gamma }^{3}\right\rangle }_{+,-}$, which together yield an overall-increasing frequency, such that Equation (55) represents a slow damping of the electric field. (Ignoring the effect of changing ${\left\langle 1/{\gamma }^{3}\right\rangle }_{+,-}$, Equation (55) gives $\overline{\hat{E}}\sim 1/{\hat{t}}^{1/4}$, a dramatically weaker damping than Equation (48).)

Further physical insight into this damping can be obtained by arriving at Equation (55) via energy arguments. Using Equation (18), we can state

$\begin{eqnarray}&&\displaystyle \frac{1}{2}{\partial }_{\hat{t}}{\hat{E}}^{2}=-\hat{j}\hat{E},\end{eqnarray} \tag{ 56 }$

which relates the change in the energy in the electric field to the power exerted on the current by the electric field. Via the same method used to obtain Equation (23), we also have

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{j}={\hat{\omega }}^{2}\hat{E}\end{eqnarray} \tag{ 57 }$

with $\hat{\omega }$ defined in Equation (24). Using Equation (57) to substitute for $\hat{E}$ in Equation (56), we can arrive at

$\begin{eqnarray}&&{\partial }_{\hat{t}}\left(\displaystyle \frac{{\hat{E}}^{2}}{2}+\displaystyle \frac{{\hat{j}}^{2}}{2{\hat{\omega }}^{2}}\right)=-\displaystyle \frac{{\hat{j}}^{2}}{2{\hat{\omega }}^{2}}\left(\displaystyle \frac{1}{{\hat{\omega }}^{2}}\displaystyle \frac{\partial {\hat{\omega }}^{2}}{\partial \hat{t}}\right).\end{eqnarray} \tag{ 58 }$

The first term on the left side represents the wave energy in the electric field and the second the energy from the velocity of particles in the wave. When the wave frequency changes slowly, we can say these are roughly equal on average and replace Equation (58) with an expression in terms of the electric field amplitude:

$\begin{eqnarray}&&{\partial }_{\hat{t}}{\overline{\hat{E}}}^{2}=-\displaystyle \frac{{\overline{\hat{E}}}^{2}}{2}\left(\displaystyle \frac{1}{{\hat{\omega }}^{2}}\displaystyle \frac{\partial {\hat{\omega }}^{2}}{\partial \hat{t}}\right).\end{eqnarray} \tag{ 59 }$

From this, we can reproduce Equation (55).

The relationship in Equation (55), with a frequency calculated from the current value of Equation (24) across the whole simulation box, and the squared electric field strength given by its average value in the box, is compared to simulation in Figure 5, showing good agreement.

Zoom In Zoom Out Reset image size

In summary, the linear stage, which occurs late in the discharge, is marked by weak damping. The nature of this damping is such that the relationship between the electric field strength and frequency throughout the linear stage is given by Equation (55).

We have developed a one-dimensional model for the evolution of the wave amplitude during a pulsar polar cap discharge. In this model, the initial inductive electric field in the pulsar polar cap is screened by the creation of electron–positron pairs, setting up oscillations that continue to be damped by pair creation. This damping process has three phases: a screening phase, described in Section 5; a nonlinear phase, described in Section 6; and a linear phase, described in Section 8. The conditions for the transition between the nonlinear stage and the linear stage are considered in Section 7.

In a realistic, two-dimensional pulsar polar cap, the discharge waves are initially electromagnetic (Philippov et al. 2020); in addition, they will refract and may develop an even larger electromagnetic component as they also encounter cross-field gradients after their creation (Melrose et al. 2020). These electromagnetic waves eventually escape the plasma as Poynting flux. Our one-dimensional model does not consider the magnetic component of the plasma oscillations and cannot treat the escape process, which will be the focus of future work. However, the electromagnetic oscillations interact with the pair plasma through a parallel wave electric field and the resulting plasma current, and this interaction, in which the physics of our model is key, is responsible for the evolution of the amplitude of the radio emission. In this conclusion, we will thus use the physics understanding from our model to give new insight into multiple elements of pulsar observations.

To start, we develop an expression for the frequency of the pulsar's emission. We note that in the linear phase of the discharge, the pair distribution function approaches a Maxwell–Jüttner (Jüttner 1911) distribution drifting at u_l (at the large Lorentz factors given in Equation (3) characterizing the lower-energy pairs, γ_l ≈ u_l , and we will present results in terms of u_l ). For such a distribution, $\left\langle 1/{\gamma }^{3}\right\rangle$ depends on u_l and on the thermal spread of newly injected pairs, which is defined by the temperature, T. In particular, for pulsar conditions, where u_l ≫ 1, T ≳ mc², and T/(mc²) ≪ u_l , the plasma frequency, Equation (24), can be roughly represented by Rafat et al. (2019):

$\begin{eqnarray}&&\omega \approx \sqrt{\left(\displaystyle \frac{8\pi {n}_{+}{e}^{2}}{{{mu}}_{l}^{3}}\right)\left(\displaystyle \frac{2T}{{{mc}}^{2}}\right)}.\end{eqnarray} \tag{ 60 }$

The pulsar's radio emission is in the form of a superluminal O mode, which has a frequency ω_O above the plasma frequency. Pair multiplicity in polar cap cascades of active pulsars relatively quickly reaches values of ∼10³ and continues to grow until reaching the maximum multiplicity of a ∼few × 10⁵ (Timokhin & Harding 2015, 2019). Electromagnetic waves likely escape the plasma when it is characterized by a pair multiplicity between these values. ¹⁰ Evaluating Equation (60) numerically for a pair multiplicity at the bottom of this range, we get

$\begin{eqnarray}\begin{array}{rcl}{\nu }_{O} & \geqslant & 5\,\mathrm{MHz}{\left(\displaystyle \frac{{{\rm{n}}}_{+}/{{\rm{n}}}_{{GJ}}}{{10}^{3}}\right)}^{1/2}{\left(\displaystyle \frac{{{\rm{B}}}_{\star }}{{10}^{12}\,{\rm{G}}}\right)}^{1/2}\\ & & \times \,{\left(\displaystyle \frac{P}{1\,{\rm{s}}}\right)}^{-1/2}{\left(\displaystyle \frac{{u}_{l}}{{10}^{3}}\right)}^{-3/2}{\left(\displaystyle \frac{T}{{{mc}}^{2}}\right)}^{1/2},\end{array}\end{eqnarray} \tag{ 61 }$

with ν_O representing the frequency of O-mode waves. Pulsar radio emission has been observed at frequencies as low as 15 MHz (see, for example; Pilia et al. 2016), making Equation (61) possibly consistent with observations. O-mode oscillations with higher wavenumbers have higher frequencies; in addition, later in the discharge process, when the plasma multiplicity increases above ∼ 10³ and u_l drops down to ∼10, the frequency will increase. Both of these effects push Equation (61) toward the radio range.

To continue, we present a set of arguments leading to an expression for pulsar radio luminosity. The luminosity depends on the amplitude of the emission when it leaves the pulsar. As we described in Section 6, all strong damping of the plasma waves occurs in the nonlinear stage. Considering Equation (52), the end of the nonlinear phase occurs when

$\begin{eqnarray}&&\displaystyle \frac{e{\overline{E}}_{\star }}{{u}_{l}m{\omega }_{t}c}\sim 1,\end{eqnarray} \tag{ 62 }$

where ${\overline{E}}_{\star }$ is the electrostatic component of the wave amplitude at the point of transition between the nonlinear and linear stages, u_l ∼ 10³ is the injection momentum of the lower-energy pairs, and ω_t ∼ 10⁹ rad s⁻¹ defines the characteristic wave frequency at the end of the nonlinear phase. Evaluation gives ${\overline{E}}_{\star }={u}_{l}m{\omega }_{t}c/e\sim {10}^{5}\,{\rm{G}}$. The time to reach the end of the nonlinear phase, t_s , can be found in Equation (49), ${t}_{s}\sim 2\sqrt{\xi }\mathrm{ln}({E}_{0}/{\overline{E}}_{\star }){t}_{0}$. Estimating $\mathrm{ln}({E}_{0}/{\overline{E}}_{\star })\lesssim 10$ and using Equation (20), one obtains

$\begin{eqnarray}&&{t}_{s}\lesssim 0.5{\left(\displaystyle \frac{P}{1\,{\rm{s}}}\right)}^{1/4}\left(\displaystyle \frac{{R}_{\mathrm{pc}}}{c}\right).\end{eqnarray} \tag{ 63 }$

This time is shorter than the light-crossing time across the polar cap, R_pc/c, so the nonlinear phase of the discharge completes in all active pulsars.

As described in Section 8, the damping of the electric field is weak during the linear stage, so the rough overall radio luminosity is likely to be set at the end of the nonlinear phase. (The end of the nonlinear phase also corresponds to the lowest frequencies of emission, because u_l will decrease and multiplicity will increase in the linear phase. The pulsar spectrum declines with frequency, such that luminosity is dominantly determined by the low-frequency contributions occurring at the end of the nonlinear stage.) Emission across the pulsar polar cap gives a radio luminosity of ${L}_{\mathrm{rad}}\sim c\left({\overline{E}}_{\star }\sin \alpha {)}^{2}(\eta \pi {R}_{\mathrm{pc}}^{2}\right)/(4\pi )$, where α is the typical angle between the direction of wave propagation and the background magnetic field (see Philippov et al. 2020) and η characterizes the fraction of pair-producing field lines in the polar cap. Using Equation (62), one obtains

$\begin{eqnarray}&&\begin{array}{l}{L}_{\mathrm{rad}}\sim {10}^{28}\mathrm{erg}\,{{\rm{s}}}^{-1}\ {\left(\displaystyle \frac{{\omega }_{t}}{{10}^{9}\,\mathrm{rad}\,{{\rm{s}}}^{-1}}\right)}^{2}\ {\left(\displaystyle \frac{P}{1{\rm{s}}}\right)}^{-1}\\ \qquad \times \,{\left(\displaystyle \frac{{u}_{l}}{{10}^{3}}\right)}^{2}\eta {\sin }^{2}\alpha .\end{array}\end{eqnarray} \tag{ 64 }$

This luminosity is toward the lower end of the range of observed radio luminosities (Szary et al. 2014; Wu et al. 2020); higher values of luminosity may be explained by the escape of some radiation before the end of the nonlinear stage. The prediction in Equation (64) can also be expressed as a fraction of the pulsar's spin-down luminosity, ${L}_{\mathrm{sd}}={{cE}}_{0}^{2}\pi {R}_{\mathrm{pc}}^{2}/(4\pi )$ (Spitkovsky 2006):

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{L}_{\mathrm{rad}}}{{L}_{\mathrm{sd}}}\sim \displaystyle \frac{{\overline{E}}_{\star }^{2}}{{E}_{0}^{2}}(\eta {\sin }^{2}\alpha )\sim 2\,\times \,{10}^{-4}{\left(\displaystyle \frac{P}{1\,{\rm{s}}}\right)}^{3}\\ \quad \times \,{\left(\displaystyle \frac{{B}_{\star }}{{10}^{12}\,{\rm{G}}}\right)}^{-2}{\left(\displaystyle \frac{{\omega }_{t}}{{10}^{9}\,\mathrm{rad}\,{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{u}_{l}}{{10}^{3}}\right)}^{2}(\eta {\sin }^{2}\alpha ).\end{array}\end{eqnarray} \tag{ 65 }$

Observations show practically no dependence of the radio luminosity on pulsar parameters (Szary et al. 2014; Wu et al. 2020); however, Equation (64) shows a nonnegligible dependence on B_* and Ω. We should note that the uncertainties in our simplified estimate for ω_t and the numerical factors η and α may substantially reduce this dependence. For example, in young energetic pulsars with large Ω (and B_*) the pair creation front is expected to be less inclined with respect to the background magnetic field, which reduces α and hence the dependence of the radio luminosity on these parameters. We will refine these dependencies in our future multidimensional studies.

Next, although the overall radio luminosity is set at the transition between the nonlinear and linear phases, some modification in the amplitude does occur in the linear phase. We expect that radiation will escape the plasma throughout the linear stage (and possibly during the nonlinear stage as well) in a continuous process that occurs concurrently with the generation of the plasma waves. The relationship between amplitude and frequency in the linear stage, given by Equation (55), would thus be visible to an observer if, across the polar cap, the wave amplitude became unaffected by the plasma at a variety of different points in the linear damping process, when the value of the plasma frequency is different. This relationship roughly agrees with the observed spectrum of emission, S_ω ∼ ω^−1.4±1.0, with S_ω the intensity emitted at that frequency (Bates et al. 2013). This suggests that the pulsar radio spectrum may result in part from the linear damping physics of the pair discharge.

Finally, the distribution of particles right after the end of the initial screening, Equation (31), can be used as a model for the spectrum of particles flowing back from the discharge to the pulsar surface, allowing greater understanding of the X-ray hotspots observed by NICER (Gendreau et al. 2016; Salmi et al. 2020).

E.T. is indebted to a helpful comment on Math StackExchange from user "DinosaurEgg." This research made use of the "Tristan-MP v2" particle-in-cell code. We thank Peter Catto, Lev Arzamasskiy, Carolyn Raithel, and Fabio Cruz for helpful conversations. Hayk Hakobyan provided expert advice on code usage. E.T. was supported by the W.M. Keck Foundation Fund at the Institute for Advanced Study. A.P. was supported by NSF grant no. PHY-2010145. A.T. was supported by the grant 2019/35/B/ST9/03013 of the Polish National Science Centre.

Works that study pulsar polar cap discharges using a fluid treatment observe sawtooth-shaped oscillations, where electric field evolution is linear except at turning points (Levinson et al. 2005). There have been tentative suggestions in the literature that such sawteeth may be a generic feature of relativistic pair discharges, extending also to kinetic treatments (Cruz et al. 2021a). In fact, sawteeth do not occur in a kinetic treatment.

To understand why, let us first consider the origin of sawteeth in a fluid treatment. We can define a Lorentz factor associated with the electron and positron fluids as, with s = + , − again indicating species,

$\begin{eqnarray}&&{\gamma }_{s}\equiv \displaystyle \frac{1}{\sqrt{1-\tfrac{{v}_{s}^{2}}{{c}^{2}}}},\end{eqnarray} \tag{ A1 }$

with v_s the velocity of species s; furthermore, we define N_s as the density of fluid species s measured in the pulsar frame. The equations governing this system are (Levinson et al. 2005)

$\begin{eqnarray}&&m{\partial }_{t}\left({N}_{s}{\gamma }_{s}{v}_{s}\right)={N}_{s}{q}_{s}E,\end{eqnarray} \tag{ A2 }$

$\begin{eqnarray}&&{\partial }_{t}{N}_{s}=S/2,\end{eqnarray} \tag{ A3 }$

and

$\begin{eqnarray}&&{\partial }_{t}E+4\pi \sum _{s}{q}_{s}{N}_{s}{v}_{s}=0,\end{eqnarray} \tag{ A4 }$

with m again representing the pair mass, q_s the species charge, E the electric field, and S/2 the source rate. Under equivalent normalizations to those used in Section 3, these become

$\begin{eqnarray}&&{\partial }_{\hat{t}}\left({\hat{N}}_{s}{\gamma }_{s}{\beta }_{s}\right)={\hat{N}}_{s}\displaystyle \frac{{q}_{s}}{e}\hat{E},\end{eqnarray} \tag{ A5 }$

$\begin{eqnarray}&&{\partial }_{\hat{t}}{\hat{N}}_{s}=1,\end{eqnarray} \tag{ A6 }$

and

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}+\displaystyle \frac{1}{\xi }\sum _{s}\displaystyle \frac{{q}_{s}}{e}{\hat{N}}_{s}{\beta }_{s}=0.\end{eqnarray} \tag{ A7 }$

Except over short timescales when the pair fluid velocities reverse direction, the electron and positron fluids are highly relativistic, such that Equation (A7) approaches

$\begin{eqnarray}&&{\partial }_{\hat{t}}\hat{E}=-\displaystyle \frac{1}{\xi }\sum _{s}\displaystyle \frac{{q}_{s}}{e}{\hat{N}}_{s};\end{eqnarray} \tag{ A8 }$

with constant or slowly changing density, this represents a linear or nearly linear change in the magnetic field, giving the sawtooth form seen in the first panel of Figure 1 of Levinson et al. (2005).

A kinetic description of the polar cap oscillations, presented in Section 3, differs fundamentally from a fluid one in that it takes into account that a plasma is composed of particles of many different velocities, unlike a fluid, which has a single velocity. Thus, the replacement presented in Equation (A8) is not possible, and the electric field evolution is sinusoidal, with a frequency given by Equation (24). We note that to obtain a fluid-like result, i.e., one with linear oscillations of ${\partial }_{\hat{t}}^{2}\hat{E}=0$, we could take the limit ${\left\langle 1/{\gamma }^{3}\right\rangle }_{+,-}\to 0$ in the frequency given by Equation (24), which corresponds to the artificial situation of a plasma composed just of particles moving with c, which bears some resemblance to a fluid moving at c.

As a further illustration of the difference between kinetic and fluid oscillations, we present in Figure 6 a comparison between the first kinetic bounce shown in Figure 2 and a calculation of the fluid analog evolution of the electric field that would occur if Equation (A8) governed the system. The difference in the two lines is significant.

Zoom In Zoom Out Reset image size

In this appendix, we develop a rough model for A and ${\hat{t}}_{b}$. To do this, we consider the behavior of Equation (37) near the beginning of a period in which electrons are decelerated, $\hat{t}=\tfrac{\pi }{2\hat{\omega }}$. Specifically, consider that the contribution of new particles to the frequency is

$\begin{eqnarray}\begin{array}{rcl}{\hat{\omega }}_{\mathrm{new}}^{2}\left(\hat{t}\right) & = & {\left(\displaystyle \frac{{\hat{n}}_{+}}{\xi }{\left\langle \displaystyle \frac{1}{{\gamma }^{3}}\right\rangle }_{+}\right)}_{\mathrm{new}}+{\left(\displaystyle \frac{{\hat{n}}_{-}}{\xi }{\left\langle \displaystyle \frac{1}{{\gamma }^{3}}\right\rangle }_{-}\right)}_{\mathrm{new}}\\ & \approx & \displaystyle \frac{1}{\xi }{\int }_{{u}_{\max }\left(\hat{t}\right)}^{{u}_{l}}{du}{\hat{f}}_{-,\mathrm{new}}\left(u,\hat{t}\right){\left(1+{u}^{2}\right)}^{-3/2},\end{array}\end{eqnarray} \tag{ B1 }$

with ${\hat{f}}_{-,\mathrm{new}}$ defined in Equation (37) and ${u}_{\max }\left(\hat{t}\right)$ defined in Equation (40). The frequency is dominated by contributions from the electrons because positrons are accelerated to even higher u than their injection momentum u_l , increasing their Lorentz factor and decreasing their contribution to the frequency. The value of A represents the height of the spike in $\xi {\hat{\omega }}^{2}$ that occurs just past $\hat{t}=\tfrac{\pi }{2\hat{\omega }}$, at the time where particles added near the zero of the electric field are pulled through u = 0. The value of $A{\hat{t}}_{b}$ represents the area under this spike. Specifically, we can write

$\begin{eqnarray}&&A\sim \xi \mathop{\max }\limits_{\tfrac{\pi }{2\hat{\omega }}\lt \hat{t}\lt \displaystyle \frac{3\pi }{2\hat{\omega }}}\left[{\hat{\omega }}_{\mathrm{new}}^{2}\left(\hat{t}\right)\right].\end{eqnarray} \tag{ B2 }$

and

$\begin{eqnarray}&&A{\hat{t}}_{b}\sim \xi {\int }_{\pi /2\hat{\omega }}^{3\pi /2\hat{\omega }}{\hat{\omega }}_{\mathrm{new}}^{2}\left(\hat{t}\right)d\hat{t}.\end{eqnarray} \tag{ B3 }$

We numerically evaluate these quantities for ${\hat{E}}_{n}\sim 1$ at several values of $\hat{\omega }$ and fit the resulting data, giving

$\begin{eqnarray}&&A\sim {\hat{\omega }}^{-0.5}\end{eqnarray} \tag{ B4 }$

and

$\begin{eqnarray}&&{\hat{t}}_{b}\sim {\hat{\omega }}^{-0.7}.\end{eqnarray} \tag{ B5 }$

We note from Equation (41) that at times of zero-crossing ${\hat{\omega }}^{2}\sim A/\xi$ and use Equation (B4) to write

$\begin{eqnarray}&&{\hat{\omega }}^{2}\sim \displaystyle \frac{1}{\xi {\hat{\omega }}^{0.5}},\end{eqnarray} \tag{ B6 }$

whence,

$\begin{eqnarray}&&A\sim {\xi }^{0.2}\end{eqnarray} \tag{ B7 }$

and

$\begin{eqnarray}&&{\hat{t}}_{b}\sim {\xi }^{0.3}.\end{eqnarray} \tag{ B8 }$

Simulation confirms that the height and width of the spikes increase slowly with ξ.

All analysis after Equation (12) assumes negligible injection temperature; our simulations use a small injection temperature of T = 0.1 mc², as described in Section 4. However, as noted in Section 3, newly injected pairs in fact have a thermal spread of roughly T ≳ mc². To examine the effect of this thermal spread, we run two more cases of our ξ = 2.8 × 10⁶ simulation with injection temperatures of T = 1 mc² and T = 10 mc². The effect of finite injection temperature on these simulations can be seen in Figure 7. The left figure, which depicts the quantity also shown in Figure 5, shows that higher-temperature simulations have a lower frequency and a higher saturation amplitude. In addition, the fragmentation into finite-k modes occurs later in the higher-temperature simulation. The right plot in Figure 7 displays the spatial profile of $\hat{E}$ at a time late in the discharge, after the simulation fragments into finite-k modes, and should be compared to the quantities in Figure 4, showing that higher-temperature simulations have slightly larger wavelength modes.

Zoom In Zoom Out Reset image size