Magnetic Energy Dissipation and γ-Ray Emission in Energetic Pulsars

We employ radiative PIC algorithm implemented in the Tristan-MP v2 code designed by Hakobyan & Spitkovsky (2020) to simulate the 3D dynamics of the entire magnetosphere. Synchrotron radiation drag force is included in the equations of motion for particles in the Landau–Lifshitz form. The magnitude of the force is artificially enhanced with respect to the Lorentz force, so the most energetic particles lose a substantial amount of energy on the gyration timescales. A similar approach for modeling synchrotron drag has been used in previous works (e.g., Cerutti et al. 2016).

Our simulations start with an empty Cartesian domain of the size $\sim {(5{R}_{\mathrm{LC}})}^{3}$, where R_LC = c/Ω is the light cylinder of the magnetosphere, and Ω is the spin frequency of the neutron star. The star itself is modeled as a perfectly conducting rotating sphere in the center of the domain (similar to Philippov et al. 2015). A dipolar magnetic field is imposed as a boundary condition near the surface of the sphere, as well as the initial condition in the whole domain. The angle between the magnetic axis and the rotation axis is further denoted by χ. Near outer boundaries, all of the electromagnetic field components are damped to zero, and the particles are allowed to leave. Particle distribution is sampled on average by ∼10 macroparticles per grid cell (PPCs), with fewer PPCs farther from the star. ⁶ To mitigate the numerical artifacts from finite PPCs, we employ digital filtering of the deposited currents with a stencil of size 2. We also tested the setup with a lower-resolution simulation and 10 times more PPCs, arriving at the same results.

To fill the magnetosphere with plasma we mimic the pair-production process near the surface of the star by injecting pair plasma in a small spherical shell of size Δr at a rate proportional to the local GJ density: Δn(θ, ϕ)/Δt = f_inj∣ρ_GJ(θ, ϕ)/e∣(c/Δr). The dimensionless parameter f_inj controls the injection multiplicity. The exact value of this number is unimportant, as long as enough plasma is injected to establish the FF solution (we typically set f_inj = 1). We also give a marginal initial velocity to the newly injected particles along the local magnetic field lines (typically a Lorentz factor of γ ≈ 2). ⁷

To be able to resolve the plasma skin depth, d_e , everywhere except for the surface of the star, we greatly reduce the scale separation compared to realistic pulsars. The radius of the star is resolved by 75Δx (we will denote Δx as the size of the grid cell), while the size of the light cylinder is R_LC ∼ 440Δx (or ∼6 times the radius of the star). The strength of the magnetic field at the surface, B_*, which also rescales the Goldreich–Julian density, is chosen in such a way that ${d}_{e}^{\mathrm{LC}}\sim \mathrm{few}\,{\rm{\Delta }}x$, where ${d}_{e}^{\mathrm{LC}}$ is the plasma skin depth near the light cylinder. The scale separation between macroscopic and microscopic (plasma-kinetic) length scales is, thus, at most ∼200 in our highest-resolution simulation (see almost 8 orders of magnitude scale separation in realistic pulsars). Large separation of scales ensures that the growth of microscopic plasma instabilities that develop at kinetic timescales (inverse plasma frequency, ${\omega }_{\mathrm{pe}}^{-1}$) is much faster than the dynamical timescale characterized by the rotation period, P. Localized 2D simulations (see, e.g., Werner et al. 2016) show that at scale separation ≳100, kinetic instabilities establish an asymptotic regime, which justifies our choice of ${R}_{\mathrm{LC}}/{d}_{e}^{\mathrm{LC}}\sim {\omega }_{\mathrm{pe}}^{\mathrm{LC}}/{\rm{\Omega }}\sim 200$. The choice of these parameters yields the total size of the simulation domain (2200Δx)³.

In regions where the magnetic field is strong (close to the surface of the star), we significantly underresolve particle gyration timescale (ω_B Δt/γ ≫ 1, where ω_B = ∣e∣B/m_e c, and γ ≈ few). To avoid associated numerical errors, we employ a coupled guiding-center/Boris algorithm for solving particle equations of motion (Bacchini et al. 2020). This approach allows us to ignore gyrations of most of the particles in the bulk of the magnetosphere where the gyration is underresolved, reducing their motion to that of their guiding centers. At the same time, switching to conventional Boris pusher in regions where E/B > 0.95 allows us to recover the full equation of motion for high-energy particles in regions with vanishing magnetic field (i.e., current sheets). ⁸ By ignoring the synchrotron drag force for particles in the guiding center regime, we can also avoid having numerical errors in the strong cooling regime, when the synchrotron cooling algorithm heavily relies on the resolution of gyration orbit.

Overall, the three dimensionless parameters that we fix independently in our simulations are:

• 1.  
  ${R}_{\mathrm{LC}}/{d}_{e}^{\mathrm{LC}}\sim 100\mbox{--}200$: the ratio of the size of the light cylinder and the plasma skin depth at a corresponding GJ plasma density. We refer to this as the scale separation of our simulation;
• 2.  
  R_*/Δx ≈ 75: the number of simulation cells per stellar radius, which we call the resolution of our simulation; and
• 3.  
  R_LC/R_* ≈ 6: the size of the light cylinder, w.r.t. the size of the star.

In addition we also have the obliquity angle, χ, which we vary from 0°–90°, and the cooling strength, which we detail below.

In all of our simulations, after a brief transient lasting for about one rotation period of the star, P, a steady-state structure is established that closely resembles the FF solution. The snapshot of the steady-state for the neutron star with an inclination angle χ = 30° is outlined show in Figure 1. Plasma close to the surface corotates with the neutron star, flowing along almost poloidal magnetic field lines that start and end at the surface of the star. This corotation extends up to R_LC = c/Ω from the rotation axis (Ω being the rotation frequency of the neutron star), the light cylinder. Outside the light cylinder, plasma can no longer corotate with the star, and pairs slide along spiral magnetic field lines beyond the R_LC, moving almost radially outward. While in the inner magnetosphere (inside the light cylinder), the magnetic field is almost purely poloidal, in the outer magnetosphere it has a toroidal components.

Zoom In Zoom Out Reset image size

The rotation of the magnetosphere imposes a poloidal electric field in the wind zone, E_θ ∼ B_ϕ . As a result, the pulsar radiates electromagnetic energy in the form of a radial Poynting flux: $(4\pi /c){\boldsymbol{S}}={\boldsymbol{E}}\times {\boldsymbol{B}}={E}_{\theta }{B}_{\phi }\hat{{\boldsymbol{r}}}$. This flux, integrated over a sphere that encloses the light cylinder, is the spin-down energy of the neutron star, often denoted as $\dot{E}$, which for an aligned rotator reads:

$$\begin{eqnarray}\begin{array}{rcl}\dot{E}\equiv {L}_{0} & = & {\unicode{∯} }_{r={R}_{\mathrm{LC}}}{\boldsymbol{S}}\cdot d{\boldsymbol{a}}\\ & = & 2\pi \displaystyle \frac{{B}_{* }^{2}{R}_{* }^{3}}{P}{\left(\displaystyle \frac{{R}_{* }}{{R}_{\mathrm{LC}}}\right)}^{3},\end{array}\end{eqnarray} \tag{ 1 }$$

where B_* is the magnetic field strength near the stellar surface at the equator.

Opposing field lines from the northern and southern hemispheres of the neutron star are divided in the outer magnetosphere by a current sheet. The structure of the equatorial current sheet is more visible in the simulation with no obliquity, shown in Figures 2(a), (b) (hereafter we refer to this simulation as R75_ang0, since χ = 0°, and R_* = 75Δx). ⁹ Cartesian coordinates are normalized to the R_LC, and (0, 0, 0) is the middle of the box, where the center of the neutron star is. Figure 2(a) shows the slice of the simulation in the x = 0 plane, while Figure 2(b) shows the slice in the z = 0 plane (in this simulation, χ = 0°). Color indicates the plasma number density compensated by the cylindrical radius squared, ${{nr}}_{{xy}}^{2};$ its value is normalized by the corresponding GJ density at the surface of the star times the radius of the star squared, ${n}_{\mathrm{GJ}}^{* }{R}_{* }^{2}$.

Zoom In Zoom Out Reset image size

From Figure 2 it is evident that the total plasma density close to the current sheet is few to ten times larger than the local GJ density, which means the magnetosphere has enough plasma to screen the accelerating electric field almost everywhere. If the parallel electric field were screened everywhere, magnetic energy dissipation in the magnetosphere would not be possible, and the integral in Equation (1) would be constant with r, $L(r)\equiv {\unicode{∯} }_{r}{\boldsymbol{S}}\cdot d{\boldsymbol{a}}=\mathrm{const}$, yielding no γ-ray emission.

The accelerating electric field can survive in microscopic subregions of the equatorial current sheet. The characteristic scale of these regions does not exceed a few-to-tens of plasma skin depths, d_e = c/ω_pe, where ω_pe is the plasma frequency for relativistic electron-positron pairs in the current sheet. For the Crab pulsar, this parameter close to the light cylinder is of the order of a few centimeters. The light cylinder, on the other hand, is almost 8 orders of magnitude larger (for the Crab, it is ∼1500 km). To properly model the dynamics of these subregions, first-principles PIC simulations are necessary.

The equatorial current sheet highlighted in Figure 2 is not laminar even in the aligned case, χ = 0°. Rather, it is prone to kinetic instabilities that develop on microscopic timescales comparable to ${\omega }_{\mathrm{pe}}^{-1}\ll {{\rm{\Omega }}}^{-1}$. Drift-kink instability displaces the current sheet in the z-direction, which results in the undulation of the sheet at a certain saturated amplitude (shown in the poloidal slice of Figure 2(b); also see Philippov & Spitkovsky 2014; Cerutti et al. 2015). The wavevector of this perturbation is perpendicular to the local magnetic field and is parallel to the local current.

More importantly, the current sheet also experiences tearing instability due to relativistic reconnection of magnetic field lines from the upstream (see, e.g., Cerutti & Philippov 2017; Philippov & Spitkovsky 2018; Cerutti et al. 2020). As a result of this process, magnetic islands ("plasmoids") are formed containing hot plasma energized in the reconnection process. In-between the plasmoids, there are magnetic nulls—regions where the magnetic field lines tear, and energy dissipation happens, the "x-points."

Figure 3(c) shows a 3D rendering of the plasma density as well as two slices of the same quantity: one along the magnetic field lines (Figure 3(a) indicated with a red surface in panel 3(c)), and the other one—perpendicular to field lines, along the direction of the current (Figure 3(b) indicated with a blue surface). In Figure 3(b) one can see the undulation of the current sheet due to drift-kink instability. The reconnection of magnetic field lines and tearing instability, on the other hand, are better visible in Figure 3(a); overdense regions in the sheet are the plasmoids, which, in 3D, look like tubes elongated almost radially (Figure 3(c)).

The dynamics of the reconnecting current sheet in Figure 3(a) is very similar to a 2D Harris sheet, except for the fact that the upstream plasma moves along the magnetic field lines with bulk Lorentz factor ${\rm{\Gamma }}\sim { \mathcal O }(1)$ (in reality, this is expected to be ${\rm{\Gamma }}\sim { \mathcal O }(100);$ see, e.g., discussion by Cerutti et al. 2020). ¹⁰ In Figures 4(a) and (b) we take the slice along the upstream magnetic field lines (similar to Figure 3(a)) where we plot quantities crucial for understanding the energy dissipation in the current sheet: the vertical component of the E × B drift velocity in units of c, and the energy dissipation rate, j · E . In the upstream (above and below the current sheet), the plasma motion is strictly constrained by magnetic field lines, which coincide with the plane of Figures 4(a) and (b). Cold plasma and opposing magnetic field lines are brought together by an $\left({\boldsymbol{E}}\times {\boldsymbol{B}}\right)$ drift (shown in Figure 4(a)); the current sheet then becomes unstable and tears producing plasmoids (shown with magenta contours). The dimensionless rate of inflow is set by the reconnection; this value can be directly measured from the simulations. For our simulations, this value is close to

$$\begin{eqnarray}&&{\beta }_{\mathrm{rec}}\equiv \displaystyle \frac{{\left({\boldsymbol{E}}\times {\boldsymbol{B}}\right)}_{\mathrm{in}}}{| {\boldsymbol{B}}{| }^{2}}\approx 0.1-0.2,\end{eqnarray} \tag{ 2 }$$

where the subscript "in" corresponds to the velocity component directed into the current sheet. This value is consistent with isolated 2D/3D simulations of the Harris sheet, and has been demonstrated to be independent of any numerical or physical parameters (such as the number of macroparticles per cell, the resolution, or the extent of the box; see, e.g., Werner et al. 2018).

Zoom In Zoom Out Reset image size

Reconnecting magnetic field in the current sheet generates a nonideal electric field in the direction of ∇ × B , which is the direction of the current in Figure 2(a) (in Figures 3(a) and 4(a) and (b), this corresponds to the out-of-plane direction). The magnitude of the nonideal electric field is of the order of E_rec ∼ β_rec B_up, where B_up is the magnetic field strength in the upstream. Because of the emergence of j · E (shown in Figure 4(b)), some of the Poynting flux radiated by the star (Equation (1)) dissipates in the magnetosphere. As seen in Figure 4(b), dissipation is concentrated exclusively in the thin current sheet outside the light cylinder, where the energy of the electromagnetic field is converted into the kinetic energy of plasma through magnetic reconnection. According to Poynting's theorem:

$$\begin{eqnarray}L(r)-{L}_{0}\equiv {\unicode{∯} }_{{R}_{\mathrm{LC}}}^{r}{\boldsymbol{S}}d{\boldsymbol{a}}=-\iiint {\boldsymbol{j}}\cdot {\boldsymbol{E}}{d}^{3}{\boldsymbol{r}}.\end{eqnarray} \tag{ 3 }$$

This means that L(r) is constant within the light cylinder, r < R_LC (where j · E = 0), but decays with radius for r > R_LC. Figure 5 shows L(r) from one of our simulations with χ = 0° computed using the flux of the Poynting vector (blue line), and the volume integral of j · E (red line). The orange band corresponds to an analytical model described further in this section.

Zoom In Zoom Out Reset image size

To understand the slope of L(r) beyond the light cylinder, it is useful to build a simple analytical model of the Poynting-flux dissipation in the current sheet. For χ = 90° rotator this model has been studied by Cerutti et al. (2020); here we focus on χ = 0°. As mentioned earlier, the dissipation of Poynting flux in our simulations is caused by the work done by the reconnection electric field, j · E , as described by Equation (3). The strength of the electric field generated due to magnetic reconnection at a given distance r from the star is E_rec(r) ∼ β_rec B_up(r). Here, the reconnecting magnetic field is a combination of poloidal and toroidal components above and below the current layer, ${B}_{\mathrm{up}}(r)=\sqrt{{B}_{r}^{2}+{B}_{\phi }^{2}}$, the value of which beyond the light cylinder can be approximated as ${B}_{\mathrm{up}}(r)={B}_{\mathrm{LC}}({R}_{\mathrm{LC}}/r)\sqrt{1+{\left({R}_{\mathrm{LC}}/r\right)}^{2}}$ (rotating monopole; see, e.g., Deutsch 1955 and Michel 1973). Here and further, ${B}_{\mathrm{LC}}\equiv {B}_{* }{\left({R}_{* }/{R}_{\mathrm{LC}}\right)}^{3}$. Since (4π/c) j = ∇ × B , the current generated during reconnection can be estimated as j ∼ cB_up/(2π δ_cs), where δ_cs is the characteristic width of the current sheet. Thus, the volume integral in Equation (3) reduces to:

$$\begin{eqnarray}&&\begin{array}{l}L(r)-{L}_{0}=-2\pi {\int }_{{R}_{\mathrm{LC}}}^{r}{dr}\,{r}^{2}{\int }_{\pi /2-{\rm{\Delta }}\theta /2}^{\pi /2+{\rm{\Delta }}\theta /2}d\theta \,\sin \theta \\ \times \,\left(\displaystyle \frac{{\beta }_{\mathrm{rec}}c}{2\pi {\delta }_{\mathrm{cs}}}{B}_{\mathrm{up}}^{2}\right).\end{array}\end{eqnarray} \tag{ 4 }$$

The integral over the polar angle, θ, at each r is accumulated in a small region near the equator (θ = π/2) of angular size Δθ ∼ δ_cs/r (which is the current sheet highlighted in Figure 4(b)). ¹¹ Integrating Equation (4) then yields

$$\begin{eqnarray}&&\displaystyle \frac{L(r)}{{L}_{0}}=1-{\beta }_{\mathrm{rec}}\left(\mathrm{ln}\displaystyle \frac{r}{{R}_{\mathrm{LC}}}+\displaystyle \frac{1}{2}\left(1-{\left[\displaystyle \frac{r}{{R}_{\mathrm{LC}}}\right]}^{-2}\right)\right),\end{eqnarray} \tag{ 5 }$$

where we also used ${L}_{0}\equiv {B}_{\mathrm{LC}}^{2}{R}_{\mathrm{LC}}^{2}c$. The logarithmic term in this expression corresponds to the dissipation of toroidal field, B_ϕ , while the second term corresponds to the dissipation of poloidal field, B_r . At large distances, r ≫ R_LC, the first term prevails, as the field is almost purely toroidal. However, for the small region, R_LC < r < 2R_LC, where γ-ray emission likely originates, contributions from both terms are important. In Figure 5 we plot the total Poynting flux, L, across a spherical shell of radius r (blue line). The orange band in Figure 5 corresponds to Equation (5) with β_rec varying from 0.09–0.12 (this number can be directly measured from Figure 4(a)). From the plot, it is evident that about 10% of the radiated Poynting flux, L₀, is dissipated in the outer magnetospheric current sheet within the first R_LC. The magnetic field strength weakens with distance, and particles can only radiate (via synchrotron) in the gigaelectronvolt range within the first one-to-few R_LC. Thus, the measured γ-ray luminosity, L_γ , should be directly correlated with the energy dissipated in this narrow range of radii. ¹² The precise value of the dissipated energy, as we have demonstrated above, is only determined by the reconnection rate at plasma microscales.

As we demonstrate in the next section, in our simulations, the plasma in the upstream part of the current sheet moves along the magnetic field lines with relativistic velocities of Γ ∼ 10. Naively, one would expect that since the rate of reconnection is universally defined in the proper frame of the upstream, in the lab frame the measured rate should be Γ times slower. Nevertheless, we do not observe any significant change in β_rec ∼ 0.1 (similar observations have been made by Cerutti et al. 2020). To check the consistency of this result in a more controlled way, we performed local 2D simulations of the current sheet, where the upstream plasma was initialized with a relativistic velocity component along the unreconnected magnetic field. From these simulations, we find that the global reconnection rate in a given reference frame is controlled by the strength of the electric field in the x-points, which are at rest. Despite the motion of the upstream, the reconnection rate remains constant, as long as the outflow velocities from the x-points remain relativistic. This condition is satisfied as long as the characteristic Alfvén four-velocity in the current sheet is much larger than the square root of the bulk four-velocity of the upstream, Γ: ${u}_{{ \mathcal A }}\sim \sqrt{\sigma }\gg \sqrt{{\rm{\Gamma }}}$. ¹³ For realistic pulsar parameters, the particle four-velocity, Γ, which (close to the light cylinder) is mostly along the magnetic field lines, is of the order of 10–100, while the magnetization is close to 10⁵–10⁶. This means that the reconnection is unlikely to slow down because of the bulk motion of the upstream.

For other values of χ, our results are consistent with earlier works (PIC: Philippov et al. 2015; and MHD: Tchekhovskoy et al. 2013). Namely, the spin-down luminosity near the light cylinder $\oint {\boldsymbol{S}}d{\boldsymbol{a}}\approx {L}_{0}(1+{\sin }^{2}\chi )$, and for larger values of χ, we see inhibited dissipation in the current sheet, as the jump of the magnetic field across the current sheet is balanced primarily by the displacement current (as opposed to the conduction current) for larger inclinations (see, e.g., Philippov & Spitkovsky 2018).

Some of the previous 2D axisymmetric simulations of aligned pulsars (Belyaev 2015; Hu & Beloborodov 2022) reported the presence of excess dissipation near the Y-point, at the level of a few percent of L₀, with a lower dissipation rate at r > R_LC. While the reduced dissipation rate is caused by the fact that in 2D simulations only one of the magnetic field components can be reconnected (only the poloidal one), the anomalous dissipation near the Y-point has unclear origins. In our 3D simulations, we also see signs of excess dissipation in the Y-point, but only in the aligned case, and only when the skin depth near the light cylinder is marginally resolved (dashed blue line in Figure 5). ¹⁴

The electromagnetic energy dissipated during the reconnection in the current sheet is deposited into plasma particles. In the next section, we focus on particle energization in the current sheet and its consequences for the observed high-energy emission in γ-ray pulsars.

In the current sheets of pulsar magnetospheres, the synchrotron cooling timescale for particles is always much shorter than the system-crossing timescale, characterized by the rotation period. However, it can be comparable to the acceleration time of the highest-energy particles, meaning that the relative importance of cooling can be significant at microscopic timescales. The efficiency of cooling can be conveniently quantified with a dimensionless number γ_rad. This number corresponds to the Lorentz factor of particles for which the synchrotron drag force in the given background magnetic field B is comparable to the Lorentz force from the reconnection-driven electric field E ∼ β_rec B. This condition reads: $| e| {\beta }_{\mathrm{rec}}B=(4/3){\sigma }_{T}{B}^{2}{\gamma }_{\mathrm{rad}}^{2}$, where σ_T is the Thomson cross section. ¹⁵ Particles with γ ≫ γ_rad will lose their energies while being exposed to a perpendicular magnetic field component much faster than they are able to accelerate; the opposite is true for particles with γ ≪ γ_rad. Notice that strictly speaking this is not applicable to particles at x-points, since the magnetic field vanishes there.

The main parameter that determines the characteristic maximum energy to which particles accelerate during reconnection is the magnetization of the upstream (unreconnected) plasma, σ. For a given background magnetic field, B, and number density of inflowing plasma, n, this dimensionless parameter is equal to twice the ratio of the magnetic energy density and the rest mass energy density of plasma:

$$\begin{eqnarray}&&\sigma \equiv \displaystyle \frac{{B}^{2}/4\pi }{{{nm}}_{e}{c}^{2}},\end{eqnarray} \tag{ 6 }$$

where B and n are measured in the proper frame, where plasma is at rest. In Figure 6 we show the magnetization parameter in the upstream of the reconnecting equatorial current sheet in a poloidal slice. In our typical simulations, magnetization of the upstream, σ^LC, is close to 500–1000, as shown with a vertical slice in Figure 6(b). This parameter, which in relativistic reconnection is always ≫1, determines the characteristic maximum energy to which particles can be accelerated in a single x-point (or x-line in three dimensions; see Figure 4(b)). If secondary reacceleration is prohibited, the nonthermal distribution function of accelerated particles extends at most to energies of a few times σ m_e c². Energized particles, however, may either then reenter the current sheet or become trapped in plasmoids and get accelerated again via secondary processes to even higher energies (Petropoulou & Sironi 2018; Hakobyan et al. 2021; Zhang et al. 2021). However, as we demonstrate below, in pulsar magnetospheres with realistic parameters, this process is suppressed due to synchrotron losses.

Zoom In Zoom Out Reset image size

Later in this section we will refer to the case when γ_rad < σ as the strong cooling regime, while the opposite case will be referred to as the weak cooling regime. The value of γ_rad (close to the light cylinder) only depends on the magnetic field strength: ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\approx {10}^{5}{\left({B}_{\mathrm{LC}}/{10}^{5}\,{\rm{G}}\right)}^{-1/2}$, which we know rather reliably for pulsars by extrapolating the magnetic field strength at the surface. The magnetization, σ^LC, on the other hand, depends on plasma density near the light cylinder, which is harder to constrain. However, we can estimate that approximately, by assuming that the cutoff frequency of the γ-ray photons corresponds to the synchrotron peak energy for the highest-energy particles with E_e ∼ 4σ^LC m_e c², ${\hslash }{\omega }_{B}^{\mathrm{LC}}{\left(4{\sigma }^{\mathrm{LC}}\right)}^{2}\approx {E}_{\mathrm{cut}}$. This provides a rough empirical estimate for σ^LC, which for the population of γ-ray pulsars varies between 10⁵ and 10⁶. Consequently, the reconnection in the current sheets of pulsars proceeds in the marginally radiative regime, where ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\sim {\sigma }^{\mathrm{LC}}$. In the next section we quantitatively investigate how the electromagnetic energy dissipation rate, as well as the particle distribution and synchrotron spectra, change depending on the ratio ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}/{\sigma }^{\mathrm{LC}}$.

In our simulations, we keep σ ≫ 1 in the outer magnetosphere (to ensure the reconnection proceeds in the ultrarelativistic regime), and vary the ratio γ_rad/σ. Following the discussion earlier, this ratio determines the relative importance of the effects of synchrotron cooling on reconnection dynamics, and, as we demonstrate later, results in qualitatively different high-energy radiation spectra. In Figure 7 we show snapshots from simulations with three different inclination angles, χ = 0°, 20°, and 60° (three rows), comparing two extremes of the cooling efficiency in each case: no cooling (left column, panels (a), (c), and (e)), and strong cooling ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}/{\sigma }^{\mathrm{LC}}\approx 1/15$ (right column, panels (b), (d), and (f)). In the same snapshot we also plot the integrated Poynting flux, L, as a function of radius with the x-axis corresponding to the radius in R_LC (shared with the x-axis of respective snapshots). We find that the general structure of the current sheet is largely unaffected by the strength of the synchrotron cooling. At low inclination angles, χ ≲ 20°, the current sheet is more intermittent with stronger cooling (compare Figure 7(a) and Figure 7(b)), but at higher inclinations, the difference is negligible. The rates of magnetic reconnection and, as a result, the Poynting-flux dissipation curves, L(r), are only marginally affected by the cooling strength for all values of χ.

Zoom In Zoom Out Reset image size

We further closely inspect two simulations with χ = 20° and χ = 60° (these simulations will be referred to as R75_ang20, and R75_ang60; all of the other parameters except for the inclination angle are the same as in R75_ang0). In the first series of runs, we keep σ^LC ≈ 500–1000 and vary ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$ to capture the following regimes: ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}/{\sigma }^{\mathrm{LC}}=1/15,\,1/3,\,2,\,\infty$ (where γ_rad = ∞ corresponds to a simulation without synchrotron cooling; for future reference, we denote these simulations as R75_ang20_gr1o15, R75_ang60_gr1o15, ..., R75_ang60_grINF).

Particle distribution functions in the current sheet are also very similar for all of the values of the cooling strength, as shown in Figures 8(a) and (c). In all of the studied cases, particles in the layer are able to accelerate to γ ∼ σ^LC in x-points, forming a power law of f ∝ γ⁻¹ − γ⁻² (for larger values of χ = 60°, the spectrum steepens from γ⁻¹ to about γ^−1.3). Only in the most strongly cooled simulation do we see that the cutoff of particle distribution is slightly shifted toward lower energies. This marginal difference is likely an effect of comparably small-scale separation in our simulations, and will most likely be unnoticeable for realistic systems. In weakly cooled simulations above γ > σ^LC, we see a smooth transition to γ^−2.5 − γ⁻³ and, eventually, to an almost exponential cutoff. Naively, one would expect that in the reconnection process with strong cooling, particles cannot gain energies higher than γ_rad m_e c²; however, in reality, the cooling becomes faster than the acceleration for these particles. As mentioned above, this is not necessarily true, as the acceleration takes place in the region of the current sheet where there is virtually no cooling (this has also been demonstrated in 2D simulations by Cerutti et al. 2014; Kagan et al. 2016; Hakobyan et al. 2019, etc.). On the other hand, in the strong cooling regime, once particles leave the accelerating regions (either back to the upstream or into plasmoids), they very quickly lose their energies without a chance to get reaccelerated again to energies larger than a few σ m_e c². As a result, the nonthermal distribution of particles, even in the case of a very strong cooling, is determined by the acceleration at x-points and extends up to γ ∼ few σ (see, e.g., Sironi 2022).

Zoom In Zoom Out Reset image size

In Figures 8(b) and (d), we show the spectra of emitted synchrotron photons for both χ = 20° and χ = 60° simulations for different cooling strengths. In all cases, we see a recurring pattern in photon spectra: a rise ν F_ν ∝ ν, a transition with a peak and a decay at higher energies. The power-law index for the distribution of photons is roughly consistent with the power-law index for the distribution of particles, i.e., f ∝ γ^−p leads to ν F_ν ∝ ν^−(p−3)/2.

Since the majority of the particles are unable to retain energies γ > γ_rad (and since σ > γ_rad), the spectral peak in photon energies, E_p, roughly corresponds to the synchrotron emission of particles with γ ∼ γ_rad: ${E}_{{\rm{p}}}\approx {\hslash }{\omega }_{B}^{\mathrm{LC}}{({\gamma }_{\mathrm{rad}}^{\mathrm{LC}})}^{2}$. This is evident from a comparison of the color bars in Figure 8(a), (b): we see a shift of the emission peak toward lower energies for simulations with stronger cooling. The spectrum, nevertheless, extends farther, as particles are still able to accelerate to energies >γ_rad in the x-points, albeit rapidly radiating away. In the case of moderate-to-weak cooling (σ < γ_rad), the majority of particles is accelerated to γ ∼ σ, and the peak is controlled by the value of σ: ${E}_{{\rm{p}}}={\hslash }{\omega }_{B}^{\mathrm{LC}}{\left({\sigma }_{}^{\mathrm{LC}}\right)}^{2}$ (see green lines in Figure 8(b), (d)).

Note that we do not capture significant intermittency at higher energies for the strong cooling simulation, which was prominent in 2D localized simulations by Hakobyan et al. (2019). Because of this, simulations with stronger cooling also have a smaller spectral cutoff energy together with a smaller spectral peak, despite the fact that the cutoff in particle distribution is unaffected. In two dimensions, the high-energy intermittency, which was responsible for extending the cutoff to values close to ℏ ω_B σ², was primarily caused by merger events of plasmoids of different sizes varying from just a few to hundreds of plasma skin depths. In our global 3D simulations, plasmoids are at most several skin depths in size, and there is only a handful of them in the current sheet. This striking difference in the separation of scales between local 2D and global 3D simulations, as well as the overall complexity in the structure of the current sheet in three dimensions, may explain the observed difference. Future large-scale localized 3D simulations of magnetic reconnection in the high-σ strong cooling regime will be able to clarify this conflict and provide a more definitive picture.

To test how particle acceleration and radiation spectra change with the value of magnetization, σ^LC, we artificially load the current sheet with extra plasma to alter the effective value of σ in the current layer (see Figure 9). For each pair injected from the surface, we initialize M − 1 additional pairs at rest in the cyan region shown in Figure 9. ¹⁶ All of the other parameters, such as ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\approx {10}^{3}$ and V_pc ≈ 5000 m_e c², are kept constant. These "secondary" pairs rapidly catch up with the bulk E × B outflow, and their ultimate effect is the decrease of the effective magnetization ∼M times (as shown in the left panels of Figure 9).

Zoom In Zoom Out Reset image size

In Figure 10 we show corresponding particle distributions in the current layer, as well as the emerging radiation spectra for all of these simulations. The blue line shows the fiducial case without additional injection. In that simulation the magnetization is σ^LC ∼ 10³, and particles in the current layer form an extended power law of f ∝ γ⁻¹ extending to γ ∼ σ^LC. In this moderate cooling regime, the peak of the radiation spectrum corresponds roughly to the synchrotron peak of the particles with $\gamma \sim {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\sim {10}^{3}$. With increasing pair density, M, magnetization at the light cylinder drops proportionally (as shown with the red and green lines in the particle spectra of Figure 10). In the cases of M = 2 and M = 20, the magnetization values near the Y-point are, respectively, σ^LC ≈ 500 and ≈50. Particles form a clear hard power-law distribution up to γ ∼ σ^LC with a steepening to f ∝ γ⁻² − γ⁻³ at higher energies, γ > σ^LC, followed by a cutoff. In these cases, pairs are still able to accelerate to energies >σ^LC m_e c² since the cooling is weak (${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\sim {10}^{3}\lt {\sigma }^{\mathrm{LC}}$). In the simulation shown with a dashed green line, we inhibit this secondary acceleration by enhancing the synchrotron cooling strength (decreasing ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$ to ≈200). These simulations clearly demonstrate that the particle acceleration potential of the pulsar outer magnetosphere is indeed determined by the local magnetization parameter, ${\gamma }_{\max }\sim {\sigma }^{\mathrm{LC}}$, rather than the electric potential drop near the polar cap, V_pc/m_e c², defined in equation Equation (A1). In pulsars these values are correlated, with σ^LC being significantly smaller:

$$\begin{eqnarray}&&{\sigma }^{\mathrm{LC}}\approx \displaystyle \frac{1}{2}{\left(\displaystyle \frac{{n}^{\mathrm{LC}}}{{n}_{\mathrm{GJ}}^{\mathrm{LC}}}\right)}^{-1}\displaystyle \frac{{V}_{\mathrm{pc}}}{{m}_{e}{c}^{2}}\ll \displaystyle \frac{{V}_{\mathrm{pc}}}{{m}_{e}{c}^{2}},\end{eqnarray} \tag{ 7 }$$

where ${n}^{\mathrm{LC}}/{n}_{\mathrm{GJ}}^{\mathrm{LC}}\gg 1$ is the plasma multiplicity near the light cylinder.

Zoom In Zoom Out Reset image size

Photon spectra in these runs are consistent with our earlier predictions (see Section 4.2). In the first run (blue line) the spectral peak is set by particles with energies $\gamma \sim {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\sim {\sigma }^{\mathrm{LC}}$. When we decrease σ^LC (while keeping the cooling strength, ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$, constant) we see that the peak is only marginally sensitive to σ^LC (which drops from 10³ in M = 1, to 200 in M = 2, and to 50 in M = 20), and is controlled by particles with $\gamma \sim {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$ (hence the similarity of the red and green curves). When we decrease ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$ to 200 for the σ^LC ∼ 50 run (green dashed line), the peak is shifted toward lower energies (despite the fact that particle distributions are only marginally affected).

In this work we closely inspect 3D PIC simulations of neutron star magnetospheres with the strong synchrotron radiation-reaction force included. Our primary focus is the energy dissipation and the plasma dynamics in the reconnecting current sheet. To be able to properly capture the kinetic physics governing the energy extraction in the current layer, we push the separation of macro-to-microscales in our simulations to ∼200.

We show that the rate of plasma inflow into the reconnecting equatorial current sheet, which is controlled by plasma kinetics at the scales of microscopic x-lines, determines the dissipation rate in the entire magnetosphere. This puts a stringent constraint on how much electromagnetic energy can be dissipated (and, ultimately, radiated) in the current layer within a few light cylinders. From simulations we find that this dissipation rate is almost insensitive to both the bulk motion of the background unreconnected plasma, and the synchrotron cooling strength, as long as the magnetization is high: σ^LC ≫ 1. The fraction of the dissipated energy within the region between R_LC < r < 2R_LC in the outer magnetosphere varies between 1% and 10%, with the exact value depending only on the pulsar inclination angle.

Due to the fast synchrotron cooling timescale, reconnection in the equatorial current layer proceeds in the radiatively efficient regime. A large fraction of the dissipated electromagnetic energy is radiated away by the accelerated pairs. Since cooling is always faster than the dynamical time of the current sheet, the amount of radiated energy is insensitive to the cooling strength and only depends on the amount of magnetic energy dissipated.

We confirm that the particle distribution is almost unaffected by synchrotron cooling, with particles in the current sheet forming a hard power-law spectrum with steep cutoff at energies of ∼σ^LC m_e c². Photon spectra, on the other hand, are evidently different. In the marginal cooling regime, ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\sim {\sigma }^{\mathrm{LC}}$, both the peak and the high-energy cutoff of the emission are controlled by the magnetization parameter, ${E}_{{\rm{p}}}\approx {E}_{\mathrm{cut}}\approx {\hslash }{\omega }_{B}^{\mathrm{LC}}{\left(\mathrm{few}\cdot {\sigma }^{\mathrm{LC}}\right)}^{2}$. In the magnetospheres where synchrotron cooling is strong, ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\lt {\sigma }^{\mathrm{LC}}$, plasmoids that carry the bulk of the plasma and contribute to the high-energy emission the most quickly cooled down to the characteristic temperatures of ${T}_{e}\sim {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}{m}_{e}{c}^{2}$ (see, e.g., Uzdensky & Spitkovsky 2014). Our simulations suggest that the peaks of the emission are thus controlled by the synchrotron emission of particles with energies comparable to T_e : ${E}_{{\rm{p}}}\approx {\hslash }{\omega }_{B}^{\mathrm{LC}}{\left({T}_{e}/{m}_{e}{c}^{2}\right)}^{2}\sim {\hslash }{\omega }_{B}^{\mathrm{LC}}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\right)}^{2}\approx 20\,\mathrm{MeV}$.

Observations by the Fermi telescope during the last decade have shown that some of the pulsars with $\dot{E}\gtrsim {10}^{34}$ erg s⁻¹ emit in γ-rays with characteristic luminosities ranging between 0.1% and 10% of the spin-down power. Our 3D PIC simulations provide strong evidence that this emission is powered by magnetic reconnection in the equatorial current sheets of pulsar magnetospheres. Reconnection occurs in the nonlinear stage of the tearing instability, which develops in the equatorial current sheet on microscopic timescales, ${\omega }_{{\rm{t}}}^{-1}\approx {\rm{\Gamma }}{w}_{\mathrm{cs}}/c$, where w_cs is the characteristic width of the sheet, and ${\rm{\Gamma }}\approx \sqrt{{R}_{\mathrm{LC}}/{R}_{* }}$ is the bulk Lorentz factor of the upstream plasma in the wind. ¹⁷ The width of the current sheet, w_cs, is defined by the pressure equilibrium, and is set by the Larmor radii of particles with characteristic Lorentz factors of $\gamma \approx ({T}_{e}/{m}_{e}{c}^{2})\approx {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$: ${w}_{\mathrm{cs}}\approx {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}{m}_{e}{c}^{2}/(| e| {B}_{\mathrm{LC}})$, where the ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$ parameter is set by radiation-reaction force. In this paper we explicitly assumed that the plasma microphysics is disentangled from the large-scale structure of the pulsar magnetosphere by taking a large enough separation of micro-to-macro scales. For this assumption to be valid, it is necessary for the tearing instability to develop much faster than the dynamical time of the magnetosphere—the pulsar rotation period, Ω⁻¹. The ratio of these timescales is always large for γ-ray pulsars, ${\omega }_{{\rm{t}}}/{\rm{\Omega }}\approx {10}^{3}{\dot{{ \mathcal E }}}_{36}/\sqrt{{{ \mathcal B }}_{12}}$, meaning that at least for the young pulsars with $\dot{E}\gtrsim {10}^{35}$ erg s⁻¹, the reconnection develops almost instantly in the current sheet near the Y-point. ¹⁸

The amount of energy dissipation within a few light cylinders beyond the Y-point, where the high-energy radiation originates, is determined by the reconnection rate, β_rec ≈ 0.1, and the pulsar inclination angle. Importantly, we find that this rate is not suppressed by the presence of the bulk flow of the upstream plasma along the magnetic field lines, as long as the Lorentz factor of the Alfvén velocity is larger than that of the bulk flow, $\sqrt{{\sigma }^{\mathrm{LC}}}\gg {\rm{\Gamma }}\approx \sqrt{{R}_{\mathrm{LC}}/{R}_{* }}$. ¹⁹ In young pulsars ${\rm{\Gamma }}\approx 20\,{{ \mathcal B }}_{12}^{1/4}/{\dot{{ \mathcal E }}}_{36}^{1/8}$, whereas $\sqrt{{\sigma }^{\mathrm{LC}}}\approx 500\,{\dot{{ \mathcal E }}}_{36}^{1/4}{{ \mathcal K }}_{4}^{-1/2}$, implying that the two are separated by at least an order of magnitude. However, even when $\sqrt{{\sigma }^{\mathrm{LC}}}\sim {\rm{\Gamma }}$ (which might be the case for the older pulsars with $\dot{E}\approx {10}^{34}$ erg s⁻¹), our findings indicate that the rate is only marginally slower (by a factor of 2 at most), which is virtually unnoticeable for the purposes of observations.

Our simulations demonstrate that the dissipated power is further radiated away via synchrotron radiation, which we identify as the main emission mechanism in our simulations. Curvature radiation in this context is negligible (contrary to what has been reported by, e.g., Kalapotharakos et al. 2018), as the dynamically strong cooling prevents particles from accelerating to energies beyond σ^LC m_e c², and the Larmor radii of the most energetic particles, ρ_L , are limited to microscopic plasma scales, ${R}_{\mathrm{LC}}/{\rho }_{L}\approx 2\cdot {10}^{4}\,{{ \mathcal K }}_{4}\gg 1$. For the Vela-like pulsars, where the discharge near the polar cap is the main source of pairs, and there is no additional pair production in the current sheet, most of the dissipated energy is emitted directly in the form of the highest-energy synchrotron photons around one to a few gigaelectronvolts, leading to characteristically high values of ${L}_{\gamma }/\dot{E}\approx 1 \% \mbox{--}10 \%$. The youngest Crab-like pulsars, on the other hand, produce a sufficient amount of high-energy photons, so that pair creation near the light cylinder becomes an important source of additional pair plasma (Lyubarskii 1996; Hakobyan et al. 2019). To estimate the fraction of the dissipated magnetic energy that goes into these secondary pairs, we need to evaluate the effective optical depth for the γ-ray photons, τ_{γ x} = f_{γ x} σ_T n_x R_LC, where f_{γ x} ≈ 0.25 (given by the maximum cross section of the pair-production process: ${\sigma }_{\gamma \gamma }^{(\max )}\approx 0.25{\sigma }_{T};$ see, e.g., Akhiezer & Berestetskij 1965), and n_x is the number density of X-ray photons, the presence of which is required for the gigaelectronvolt photons to pair produce. ²⁰ The effective number density of the low-energy photons can be estimated as ${n}_{x}\approx ({L}_{x}/c{\varepsilon }_{x})/(4\pi {R}_{\mathrm{LC}}^{2})$, where L_x is the luminosity in X-rays, and ε_x is the characteristic energy of X-ray photons. The pair-production cross section is the largest for the photons that satisfy ${\varepsilon }_{x}{\varepsilon }_{\gamma }\approx {({m}_{e}{c}^{2})}^{2}$ (here, ε_γ ≈ GeV). Substituting all of the values, and taking ${L}_{x}/\dot{E}\sim 1 \%$, which is the case for the youngest pulsars (Marelli et al. 2011), we finally find: ${\tau }_{\gamma x}\approx {10}^{-2}{\dot{{ \mathcal E }}}_{36}^{5/4}/\sqrt{{{ \mathcal B }}_{12}}$. For the Crab, $\dot{E}\approx 4\cdot {10}^{38}$ erg s⁻¹, and we find that τ_{γ x} ≈ 1. The relatively large value for the optical depth means that the nonnegligible fraction of the dissipated energy is deposited into secondary pairs produced in the reconnection upstream. Characteristic energies of the produced pairs can be estimated as ${\gamma }_{\sec }\approx {E}_{\gamma }/({m}_{e}{c}^{2})$, where E_γ ≈ 0.1–1 GeV is the characteristic energy of the pair-producing photons. For the Crab, ${\gamma }_{\sec }\approx {10}^{3}$, and the synchrotron emission of these secondary pairs corresponds to energies ${E}_{\sec }\approx {\hslash }{\gamma }_{\sec }^{2}| e| {B}_{\mathrm{LC}}/({m}_{e}c)\approx 1\,\mathrm{keV}$.

The maximum energy to which particles accelerate during the reconnection process is controlled by the plasma magnetization parameter near the light cylinder, few-to-ten σ^LC, which in turn determines the cutoff energy in the observed emission. The synchrotron cutoff can be expressed as ${E}_{\mathrm{cut}}\approx (3/2){\hslash }{\gamma }_{\mathrm{cut}}^{2}| e| {\tilde{B}}_{\perp }/({m}_{e}c)$, where ${\tilde{B}}_{\perp }^{2}={\left({\boldsymbol{E}}+{\boldsymbol{\beta }}\times {\boldsymbol{B}}\right)}^{2}-{\left({\boldsymbol{\beta }}\cdot {\boldsymbol{E}}\right)}^{2}$ (see, e.g., Cerutti et al. 2016), and γ_cut ≈ few · σ^LC. For reconnection in the strong cooling regime, most of the dissipated energy is radiated away via emission from the highest-energy pairs. The dissipation power is equal to ∣e∣E_rec c, where E_rec ≈ β_rec B_LC, while the synchrotron radiation power is ${\sigma }_{T}c{\gamma }_{\mathrm{cut}}^{2}{\tilde{B}}_{\perp }^{2}/(4\pi )$. From the definition of ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$, we may then find: ${\gamma }_{\mathrm{cut}}{\tilde{B}}_{\perp }\approx {\gamma }_{\mathrm{rad}}^{\mathrm{LC}}{B}_{\mathrm{LC}}$, where B_LC is the upstream magnetic field strength. Substituting this into the expression for the cutoff energy yields: ${E}_{\mathrm{cut}}\approx (3/2){\hslash }{\left({\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\right)}^{2}| e| {B}_{\mathrm{LC}}/({m}_{e}c)\left({\gamma }_{\mathrm{cut}}/{\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\right)$ $\approx (20\,\mathrm{MeV})\left({\gamma }_{\mathrm{cut}}/{\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\right)$. The cutoff energy is then uniquely set by the ratio ${\gamma }_{\mathrm{cut}}/{\gamma }_{\mathrm{rad}}^{\mathrm{LC}}$, resulting in ${E}_{\mathrm{cut}}\approx (200\,\mathrm{MeV}){\dot{{ \mathcal E }}}_{36}^{7/8}{{ \mathcal B }}_{12}^{-1/4}/{{ \mathcal K }}_{4}$. For Vela-like pulsars, κ ≈ 10³ − 10⁴ (Timokhin & Harding 2015), $\dot{E}\approx {10}^{36}-{10}^{37}$ erg s⁻¹, resulting in E_cut ≈ 1–10 GeV. ²¹

For the most energetic Crab-like pulsars with a strong magnetic field near the light cylinder (B_LC ∼ 4 × 10⁶ G for Crab), pair creation near the Y-point and beyond is important. The amount of produced pairs is proportional to the optical depth, τ_{γ x} , estimated earlier, and the number of gigaelectronvolt photons. The production rate can be estimated as ${\dot{N}}_{\pm }\approx {\tau }_{\gamma x}{L}_{\gamma }/{\varepsilon }_{\gamma }$, where ε_γ ≈ GeV. To obtain the multiplicity κ, we compare this with the GJ flux, ${\dot{N}}_{\mathrm{GJ}}\approx \sqrt{c\dot{E}}/| e|$, where we assumed that the ${\dot{N}}_{\mathrm{GJ}}e$ is equal to the GJ current: $\sqrt{c\dot{E}}$. We then find the multiplicity of pair production near the light cylinder ${\kappa }_{\mathrm{LC}}={\dot{N}}_{\pm }/{\dot{N}}_{\mathrm{GJ}}\approx 40\,{\dot{{ \mathcal E }}}_{36}^{7/4}{{ \mathcal B }}_{12}^{-1/2}({L}_{x}/1 \% \dot{E})({L}_{\gamma }/0.1 \% \dot{E})$. For the Crab pulsar, ${L}_{x}/\dot{E}\approx 1 \%$, and ${L}_{\gamma }/\dot{E}\approx 0.1 \%$ yielding κ_LC ≈ 10⁶. For the Vela, on the other hand, ${L}_{x}/\dot{E}\approx 0.01 \%$, and ${L}_{\gamma }/\dot{E}\approx 1 \%$ yielding κ_LC ≈ 3 (similar to what has been obtained by Lyubarskii (1996). ²² Substituting κ ≈ κ_LC for the Crab to the relation for the cutoff energy, we then find E_cut ≈ 0.1–1 GeV.

To reiterate, while the amount of the dissipated magnetic energy is uniquely set by the microphysics of the reconnection and the inclination angle, the luminosity of the observed emission can be somewhat different between the high-$\dot{E}$ Crab-like, and the low-$\dot{E}$ Vela-like pulsars. The cutoff energies are dictated by the highest-energy pairs, and in both cases, the numbers are close to the observed one to a few gigaelectronvolts. Spectral peaks, on the other hand, are different. In the case of Vela-like pulsars, most of the dissipated Poynting flux is radiated away via synchrotron at energies close to the cutoff around 1–10 GeV. For the Crab-like pulsars, however, secondary pair production near the light cylinder is important, and a large fraction of the dissipated energy goes into pairs that further re-radiate via synchrotron at a range of energies spanning all the way to kiloelectronvolts, resulting in a much broader spectrum. This yields a smaller γ-ray luminosity in the Fermi band (between 0.1 and 100 GeV), which is clearly observed by Abdo et al. (2013), who reported $\dot{E}\propto {L}_{\gamma }^{1/2}$ beyond $\dot{E}\gtrsim {10}^{36}$ erg s⁻¹.

In our simulations we did not include self-consistent pair production (neither the γ B → e^± process near the polar cap, nor the γ γ → e^±, Breit-Wheeler, process near the light cylinder), instead choosing to provide a constant inflow of plasma from the surface of the star. Our reasonable assumption was that as long as enough plasma is provided to the magnetosphere (n ≳ n_GJ), its general structure, as well as the dynamics of the current layer would not depend on how exactly the plasma was supplied. However, as mentioned earlier, for the youngest Crab-like pulsars, the secondary pairs produced by the two-photon pair production in the outer magnetosphere can carry a large fraction of the dissipated energy, re-radiating it at optical to X-ray frequencies. Additionally, as was shown in 2D simulations of an isolated current sheet (Hakobyan et al. 2019), and as mentioned in Section 4.3 of this paper, these secondary pairs inhibit the acceleration process during reconnection, thus controlling the extent of the γ-ray signal. To understand this process in more details, and also study the intermittency, the light curves, and the polarization characteristics of the distinct X-ray signal these pairs produce, global 3D simulations with proper two-photon pair-production physics are necessary. The Tristan-MP v2 code (Hakobyan & Spitkovsky 2020), used in this paper, supports the capabilities of tracking individual photons and modeling their interaction pair by pair. In the future we plan to apply this capability, coupled with the hybrid guiding center particle pusher, to model the pair-production process and its dynamics in the global context of the magnetosphere more realistically. While we think our main conclusions outlined in Section 5.2 will still hold, these more robust simulations will allow us to more self-consistently predict the optical/X-ray luminosity for the young pulsars (which we simply postulated empirically for the purposes of our estimations), to understand the long-term evolution of the secondary pairs, as well as their effect on the structure of the pulsar wind.