Magnetic energy dissipation and gamma-ray emission in energetic pulsars

Some of the most energetic pulsars exhibit rotation-modulated gamma-ray emission in the $0.1$ to $100$ GeV band. The luminosity of this emission is typically $0.1\text{-}10\%$ of the pulsar spin-down power (gamma-ray efficiency), implying that a significant fraction of the available electromagnetic energy is dissipated in the magnetosphere and reradiated as high-energy photons. To investigate this phenomenon we model a pulsar magnetosphere using 3D particle-in-cell simulations with strong synchrotron cooling. We particularly focus on the dynamics of the equatorial current sheet where magnetic reconnection and energy dissipation take place. Our simulations demonstrate that a fraction of the spin-down power dissipated in the magnetospheric current sheet is controlled by the rate of magnetic reconnection at microphysical plasma scales and only depends on the pulsar inclination angle. We demonstrate that the maximum energy and the distribution function of accelerated pairs is controlled by the available magnetic energy per particle near the current sheet, the magnetization parameter. The shape and the extent of the plasma distribution is imprinted in the observed synchrotron emission, in particular, in the peak and the cutoff of the observed spectrum. We study how the strength of synchrotron cooling affects the observed variety of spectral shapes. Our conclusions naturally explain why pulsars with higher spin-down power have wider spectral shapes and, as a result, lower gamma-ray efficiency.


INTRODUCTION
In γ-ray pulsars a significant fraction of the spin-down power (between 0.1% and 10%) is converted into high energy photons (Abdo et al. 2013). This suggests that somewhere in the magnetosphere the Poynting flux radiated by the rotating star is efficiently dissipated and converted into the energy of pairs populating the magnetosphere and, ultimately, to γ-rays. The equatorial current sheet beyond the light cylinder, where magnetic reconnection takes place, is a likely candidate where this energy extraction can take place. The emergence of the current sheet in plasma-filled magnetospheres of neutron stars was predicted in the force-free formulation (Contopoulos et al. 1999;Gruzinov 2005;Timokhin 2006), while kinetic plasma simulations allowed to model them self-consistently, capturing the process of magnetic reconnection (Philippov & Spitkovsky 2014;Chen & Beloborodov 2014;Cerutti et al. 2015;Belyaev 2015;Kalapotharakos et al. 2018;Cerutti et al. 2020). 2D simulations of axisymmetric magnetospheres were able to achieve significant separation of scales between the macroscopic extent of the current layer and the micro-scopic plasma-kinetic scale. However, the dynamics of a 2D reconnecting current sheet may be different from that in 3D. Current layers in 3D are prone to kinetic instabilities, such as the kink instability, which may disrupt the layer, interfering with the magnetic reconnection and tearing instability, and potentially leading to the suppresion of dissipation rate (Guo et al. 2021;Werner & Uzdensky 2021;Zhang et al. 2021). Neutron stars which have magnetic axes misaligned with respect to their rotation axes have intrinsically nonaxisymmetric magnetospheres and, as a result, more complex structures of current layers which can only be studied in 3D.
In global 3D simulations large scale separation is numerically challenging and expensive, and yet critically necessary to resolve the hierarchical chain of plasmoids that occur in the magnetic reconnection and mediate the particle acceleration process. As a result, while previous 3D simulations properly capture the general structure of the magnetosphere, the complex dynamics of the current layer is still largely unexplored due to the rather limited scale separation in prior works. In this work we make an attempt to properly capture the microphysics of the current sheet in 3D PIC simulations by extending the separation of macro-to-micro scales to 100. We achieve this by employing a hybrid approach for particle pushing in our simulations, where the motion of particles in highly magnetized regions is reduced to that of their guiding centers, while in the accelerating regions the full motion is recovered (Bacchini et al. 2020). This approach also enables us to include a strong selfconsistent synchrotron radiation-reaction force acting on the subsect of particles for which gyration is resolved.
The goal of the present work is to quantify the amount of magnetic energy dissipation, determine what plasma parameters control the average and the maximum energy that particles gain during the acceleration process, and what role does the strong synchrotron cooling, present in most of the energetic pulsars, play. We begin in section 2 with a brief introduction to the 3D structure of the magnetosphere, and the numerical setup we employ. We then describe the energy dissipation process in the reconnecting current layer and demonstrate that the rate of this process controls the overall dissipation rate in the magnetosphere (section 3). In section 4 we study the particle acceleration and high-energy radiation in the regimes of dynamically strong and weak synchrotron cooling. Section 5 summarizes our main findings, putting them in the context of Fermi observations. In particular, we show how the interplay of radiative losses and particle acceleration in magnetic reconnection explains the observed diversity of γ-ray spectra.

PLASMA-FILLED PULSAR MAGNETOSPHERES
Pair cascade near the neutron star surface has long been thought to populate the neutron star magnetospheres with pair plasma (Sturrock 1971;Ruderman & Sutherland 1975). When the magnetosphere has enough plasma supply, n e ρ GJ /|e|, its structure relaxes to the force-free (FF) solution, where the electric field component parallel to the magnetic field vanishes everywhere, E · B = 0. Charge density necessary to sustain FF solution is called the Goldreich-Julian (GJ) density, ρ GJ = Ω · B/(2πc) (Goldreich & Julian 1969). For the most energetic pulsars pair-producing cascade is highly efficient and can populate the magnetosphere with pair multiplicities up to n e |e|/ρ GJ ∼ 10 4 (Timokhin & Harding 2015. In this work we only consider neutron star magnetospheres with the abundant pair plasma supply, where deviations from the FF solution are at the scale of the local plasma skin-depth, d e .

Numerical setup
We employ radiative particle-in-cell (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 time scales. Similar approach for modeling synchrotron drag has been used in the previous works (e.g., Cerutti et al. 2016).
Our simulations start with an empty Cartesian domain of the size ∼ (5R 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). 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 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 (PPC), with fewer PPC further from the star. 1 To mitigate the numerical artifacts from finite PPC we employ digital filtering of the deposited currents with a stencil of size 2. We also tested the setup with smaller resolution simulation and 10 times more PPC, arriving to 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). 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 force-free 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). 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 LC e ∼ few ∆x, where d LC e 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 (cf. 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, ω −1 pe ) 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 LC /d LC e ∼ ω LC pe /Ω ∼ 200. The choice of these parameters yields the total size of the simulation domain (2200∆x) 3 .
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 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 to recover the full equation of motion for high-energy particles in regions with vanishing magnetic field (i.e., current sheets). 3 By 2 It is important to mention that our technique of providing the magnetosphere with fresh plasma is different from the more selfconsistent methods used in previous works (Chen & Beloborodov 2014;Philippov et al. 2015). To simplify the simulation as well as to reduce the number of input parameters, we do not model the full pair production process. However, as was demonstrated in the earlier works Belyaev 2015), the structure of the outer magnetosphere is unaffected by the pair injection prescription near the polar cap. 3 More details on how the coupled particle pusher algorithm works can be found in the appendix B.
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: ∼ 100-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; • R * /∆x ≈ 75: the number of simulation cells per stellar radius, which we call the resolution of our simulation; • 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 • to 90 • , and the cooling strength, which we describe in details further.

Steady-state numerical solution
In all of our simulations, after a brief transient lasting for about 1 rotation period of the star, P , a steady-state structure is established which 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.
The rotation of the magnetosphere imposes a poloidal electric field in the wind zone, E θ ∼ B φ . As a result, pulsar radiates electromagnetic energy in the form of a radial Poynting flux: (4π/c)S = E × B = E θ B φr . This flux, integrated over a sphere that encloses the light cylinder, is the spin-down energy of the neutron star, often denoted asĖ, which for an aligned rotator reads: Poloidal slice of the 3D plasma density normalized to n * GJ R 2 * /r 2 from a simulation of an inclined rotator with χ = 30 • (here n * GJ is the GJ density near the pole of the star). The surface of the neutron star and its rotation axis are shown in blue. Magnetic field lines traced from the surface of the star, as well as the direction of the magnetic moment, are shown in green. The light cylinder, RLC, shown with white dashed lines, separates the inner magnetosphere from the outer magnetosphere and the wind. Current sheet originating near the Y-point and spreading into the outer magnetosphere is clearly visible.
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 better visible in the simulation with no obliquity shown in Figure 2a,b (hereafter we refer to this simulation as R75 ang0, since χ = 0 • , and R * = 75∆x). 4 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 2a shows the slice of the simulation in the x = 0 plane, while Figure 2b 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 2 xy ; its value is normalized by the corresponding GJ density at the surface of the star times the radius of the star squared, n * GJ R 2 * . 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 Eq. (1) would be constant with r, L(r) ≡ ‚ r S · da = const, yielding no γ-ray emission. The accelerating electric field can survive in microscopic sub-regions of the equatorial current sheet. The characteristic scale of these regions does not exceed fewto-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 kilometers). To properly model the dynamics of these sub-regions first-principles PIC simulations are necessary.

DISSIPATION OF SPIN-DOWN ENERGY
The structure of the magnetosphere in PIC, as was demonstrated in the previous section, is very similar to that in force-free. This should come as no surprise, since in our simulations E is perfectly screened with abundant plasma, and magnetized particles follow field lines everywhere except for the equatorial current sheet where the magnetic field vanishes. The key difference from FF solutions is the dynamics of the current sheet. In ideal FF the energy dissipation in the current sheet is due to the finite resolution of the grid, and is thus numerical in its nature. PIC algorithm, on the other hand, enables us to "resolve" this dissipation on plasma kinetic scales, by capturing the reconnection of magnetic field lines from first principles.

Kinetic instabilities in the current sheet
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 ω −1 pe Ω −1 . Driftkink instability displaces the current sheet in z-direction, which results in the undulation of the sheet at a certain saturated amplitude (shown in the poloidal slice of Figure 2b; also see Philippov & Spitkovsky 2014;Cerutti et al. 2015). The wave-vector 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 3c shows a 3D rendering of the plasma density as well as two slices of the same quantity: one along the magnetic field lines ( Figure 3a indicated with a red surface in panel 3c), and the other one -perpendicular to field lines, along the direction of the current ( Figure 3b indicated with a blue surface). In Figure 3b 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 3a; overdense regions in the sheet are the plasmoids, which in 3D look like tubes elongated almost radially ( Figure 3c).
The dynamics of the reconnecting current sheet in slice 3a 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 Γ ∼ O(1) (in reality this is expected to be Γ ∼ O(100); see, e.g., discussion by Cerutti et al. 2020). 5 In Figure 4a,  . 3D rendering of the plasma density (compensated by r 2 ) from a simulation of an aligned rotator (R75 ang0). 3D rendering (panel c) is accompanied by two slices (a and b, also visible in 3D as red and blue surfaces). One of the slices (a) is orthogonal to the equator and curves along the upstream magnetic field lines, while the second one (b) is along the direction of the current perpendicular to the magnetic field. Overdense elongated tubes in panel (c) are the 3D flux ropes, the plasmoids, which are produced as a result of the non-linear tearing instability. In a 2D slice (a) the dynamics of both the current sheet and the plasmoids look very similar to those in isolated current sheet simulations. In panel (b) the drift-kink instability is visible, which deforms the current sheet in the direction perpendicular to the direction of the tearing instability. An animated version of this figure is available at the following link: https://youtu.be/-YXJ4yTlhWw. 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 the Figure 4a,b. Cold plasma and opposing magnetic field lines are brought together by an (E × B) drift (shown in Figure 4a); 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 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).

Magnetic energy dissipation via reconnection
Reconnecting magnetic field in the current sheet generates a non-ideal electric field in the direction of ∇×B, which is the direction of the current in Figure 2a Figure 4a,b this corresponds to the out-ofplane direction). The magnitude of the non-ideal 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 4b) some of the Poynting flux radiated by the star (Eq. 1) dissipates in the magnetosphere. As seen in Figure 4b 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 . . Figure 4. Two panels show different quantities in the same slice as Figure 3a. (a) shows the E × B inflow rate into the current sheet; this corresponds to the rate at which the reconnection of magnetic field lines occur, roughly, βrec ∼ 0.1. In panel (b) we plot the work done by the electric field (compensated by r 3 xy ), which also indicates the electromagnetic energy dissipation rate. The latter is localized in the thin equatorial current sheet. On both panels we overplot the plasmoids in magenta (their boundaries are selected as contours of large overdensities clearly visible in Figure 3a). of plasma through magnetic reconnection. According to Poynting's theorem: 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). Orange band corresponds to an analytical model described further in this section.
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 Eq. (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 up (r) = B 2 r + B 2 φ , the value of which beyond the light cylinder can be approximated as B up (r) = B LC (R LC /r) 1 + (R LC /r) 2 (rotating monopole; see, e.g., Deutsch 1955, andMichel 1973). Here and further, B LC ≡ B * (R * /R LC ) 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 Eq. (3) reduces to: (4) 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 4b). 6 Integrating Eq. (4) then yields where we also used L 0 ≡ B 2 LC R 2 LC 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 Fig. 5 we plot the total Poynting flux, L, across a spherical shell of radius r (blue line). Orange band in Figure 5 corresponds to Eq. (5) with β rec varying from 0.09 to 0.12 (this number can be directly measured from Figure 4a). From the plot it is evident that about 10% of the radiated Poynting flux, L 0 , 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 GeV 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. 7 The precise value of the dissipated energy, as we have demon-strated above, is only determined by the reconnection rate at plasma microscales.
5∆x. The red line shows the j · E, which accounts for the dissipation rate of the electromagnetic energy. Orange band is the theoretical fit, Eq. (5), with βrec ≈ 0.09-0.12. L0 is computed from Eq. (1). A slight increase of L with respect to L0 in the inner magnetosphere is due to the fact that the Y-point in our simulation is slightly inside the r/RLC = 1; this leads to slightly more open field lines and thus more Poynting flux. Weak sycnhrotron cooling of the particles is present in both of these simulations; however, as we show further, the effect of cooling on the reconnection rate in this high magnetization regime is negligible.
As we demonstrate in the next section, in our simulations the plasma in the upstream 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, reconnection rate remains constant, as long as the outflow velocities from the x-points remain rela-tivistic. 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 A ∼ √ σ √ Γ. 8 For realistic pulsar parameters the particle four-velocity, Γ, which close to the light cylinder is mostly along the magnetic fieldlines, is of the order of 10-100, while the magnetization is close to 10 5 -10 6 . 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, andMHD: Tchekhovskoy et al. 2013). Namely, the spin-down luminosity near the light cylinder¸Sda ≈ L 0 (1 + sin 2 χ), 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 2021) reported the presence of excess dissipation near the Ypoint, at the level of a few percent of L 0 , 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). 9 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.

γ-RAY EMISSION
Relativistic magnetic reconnection is known to produce non-thermal particle population (Sironi & Spitkovsky 2014;Guo et al. 2014;Werner et al. 2016); reconnection in the current sheets of neutron star magnetospheres is no exception. These energized particles are believed to radiate via synchrotron mechanism, pro-8 More details of our findings together with some discussion can be found in the appendix C. 9 We see this at resolution of (860∆x) 3 , where d LC e ∼ 0.5∆x. For higher resolution of (2200x) 3 considered, e.g., in Fig. 4, and d LC e ∼ 5∆x we see no signs of this effect.
ducing the pulsed high-energy emission observed in γray pulsars Philippov & Spitkovsky 2018). However, the synchrotron drag, the strength of which depends on the value of the magnetic field outside the current sheet, can inhibit the acceleration process, greatly affecting the outgoing photon spectrum.

Dynamical importance of synchrotron drag
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 highestenergy 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|β rec B = (4/3)σ T B 2 γ 2 rad , where σ T is the Thomson cross-section. 10 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: with B and n being measured in the proper frame where plasma is at rest. In Figure 6 we show the magnetization parameter in the upstream of the reconnecting 10 In terms of PIC simulations, defining dimensionless γ rad is equivalent to upscaling the classical electron radius, re, to boost the relative efficiency of synchrotron losses. Also note that in this rough definition the magnetic field is considered to be exactly perpendicular to the motion of a particle. In our simulations, as well as in reality, particles flying at small pitch angles w.r.t. the magnetic field will be cooled less efficiently. Thus, γ rad can be thought of as just a proxy for the average cooling efficiency for a population of isotropically distributed particles.
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 6b. This parameter, which in relativist reconnection is always 1, determines the characteristic maximum energy to which particles can be accelerated in a single xpoint (or x-line in 3D, see Figure 4b). If secondary reacceleration is prohibited, the non-thermal distribution function of accelerated particles extends at most to energies of a few times σm e c 2 . 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 supressed due to synchrotron losses.
Further 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: γ LC rad ≈ 10 5 B LC /10 5 G −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 2 : ω LC B (4σ LC ) 2 ≈ E cut . This provides a rough empirical estimate for σ LC , which for the population of γ-ray pulsars varies between 10 5 and 10 6 . Consequently, the reconnection in the current sheets of pulsars proceeds in the marginally radiative regime, where γ LC rad ∼ σ 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 γ LC rad /σ LC .

Particle acceleration and synchrotron spectra in the radiatively cooled magnetosphere
In our simulations we keep σ 1 in the outer magnetosphere (to ensure the reconnection proceeds in the ultra-relativistic 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 • , 60 • (three rows), comparing two extremes of the cooling effi- ciency in each case: no cooling (left column, panels a, c, e), and strong cooling γ LC rad /σ LC ≈ 1/15 (right column, panels b, d, 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 7a and Figure 7b), but at higher inclinations the difference is negligible. The rate 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 χ.
Particle distribution functions in the current sheet are also very similar for all the values of the cooling strength, as shown in Figures 8a,c. In all of the studied cases par-ticles in the layer are able to accelerate to γ ∼ σ LC in x-points, forming a power law of f ∝ γ −1 -γ −2 (for larger values of χ = 60 • the spectrum steepens from γ −1 to about γ −1.3 ). Only in the most strongly cooled simulation we see that the cutoff of particle distribution is slightly shifted towards lower energies. This marginal difference is likely an effect of comparably small scale separation in our simulations, and will most likely be unnoticable for realistic systems. In weakly cooled simulations above γ > σ LC we see a smooth transition to γ −2.5 -γ −3 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 2 , as for these particles the cooling becomes faster than the acceleration. 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. 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 2 . As a result, the non-thermal 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).
In Figures 8b,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. 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 : This is evident from a comparison of the colorbars in 8a,b: we see a shift of the emission peak towards lower energies for simulations with stronger cooling. The spectrum, nevertheless, extends further, 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 p = ω LC B (σ LC ) 2 (see green lines in Figures 8b,d). . Density snapshots from late stages (more than 2 rotations) of six different simulations with varying obliquity angle, χ, and synchrotron cooling efficiency, γ LC rad /σ LC . All of the slices are in the poloidal plane. Rows correspond to, respectively, χ = 0 • (a, b), 20 • (c, d), and 60 • (e, f). Simulations in the left column (a, c, e) have synchrotron cooling turned off, while the ones on the right (b, d, f) have very strong synchrotron cooling, γ LC rad /σ LC ≈ 1/15. Each panel also contains 1D plots of L(r)/L0 (white line; L0 is the flux measured at r = 0.75RLC) with horizontal axes, r, having the same scale as the x-axes of snapshots. Green bands correspond to the dissipation due to reconnetion estimated by the Eq. (5) with the rate varying between 0.8 and 0.15.
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 smaller spectral cutoff energy together with a smaller spectral peak, despite the fact that the cutoff in particle distribution is unaffected. In 2D the high energy intermittency, which was responsible for extending the cutoff to values close to ω B σ 2 , 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 3D, 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.

Effects of additional pair-loading
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 Figure 8. (a, c) particle and (b, d) photon spectra for the R75 ang20, and R75 ang60 simulations in different synchrotron cooling regimes. Effective σ at the light cylinder is marked with yellow stripes. Three smaller bars of different colors in (a) and (c) indicate the effective γ rad . Colorbar at the top of both panels puts particle energies into correspondence with synchrotron peak energies, E ∝ γ 2 BLC. Photon energies are normalized to E0 ∝ σ LC 2 BLC. While particle spectra look almost identical (except for the strongest cooled case), peaks of photons are shifted to smaller energies for smaller γ rad /σ. Only particles in the current layer are accounted for. Simulations without synchrotron cooling (black lines) are not shown in panels (b) and (d), since we do not collect photons in that case.
initialize M − 1 additional pairs at rest in the cyan region shown in Figure 9. 11 All the other parameters, such as γ LC rad ≈ 10 3 and V pc ≈ 5000 m e c 2 , 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).
In Figure 10 we show corresponding particle distributions in the current layer, as well as the emerging ra-11 For computational purposes we limit injection only to fieldlines close to the equator that participate in reconnection. Since particles are well magnetized and any transport in the transverse direction to the magnetic field is suppressed, the fieldlines at higher altitudes have little-to-no effect on the reconnection process.
diation spectra for all these simulations. The blue line shows the fiducial case without additional injection. In that simulation the magnetization is σ LC ∼ 10 3 , and particles in the current layer form an extended powerlaw of f ∝ γ −1 extending to γ ∼ σ LC . In this moderate cooling regime the peak of radiation spectrum corresponds roughly to the synchrotron peak of the particles with γ ∼ γ LC rad ∼ 10 3 . With increasing pair density, M , magnetization at the light cylinder drops proportionally (as shown with red and green lines in the particle spectra of Figure 10). In 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 ∝ γ −2 -γ −3 at higher energies, γ > σ LC , followed Figure 9. Plasma density, compensated for the falloff with the distance, nr 2 /(n * GJ R 2 * ), for three simulations R75 ang0 with different pair loading rates. Additional injection region is highlighted with cyan. Newly created particles are initialized with zero velocity. From top to bottom, M = 1 (no extra injection, all particles originate at the surface), M = 2 (one particle is injected per each surface-injected particle), M = 20. To the left of each panel we plot the magnetization of the steady-state solution for each case as a function of vertical coordinate (corresponding slice is shown with a white dashed line). by a cutoff. In these cases pairs are still able to accelerate to energies > σ LC m e c 2 since the cooling is weak (γ LC rad ∼ 10 3 < σ LC ). In the simulation shown with a dashed green line we inhibit this secondary acceleration by enhancing the synchrotron cooling strength (decreasing γ LC rad to ≈ 200). These simulations clearly demonstrate that the particle acceleration potential of the pulsar outer magnetosphere is indeed determined by the local magnetization parameter, γ max ∼ σ LC , rather than the electric potential drop near the polar cap, V pc /m e c 2 , defined in equation Eq. (A1). In pulsars these values are correlated, with σ LC being significantly smaller: where n LC /n LC GJ 1 is the plasma multiplicity near the light cylinder.
Photon spectra in these runs are consistent with our earlier predictions (see Sec. 4.2). In the first run (blue line) the spectral peak is set by particles with energies γ ∼ γ LC rad ∼ σ LC . When we decrease σ LC (while keeping the cooling strength, γ LC rad , constant) we see that the peak is only marginally sensitive to σ LC (which drops from 10 3 in M = 1, to 200 in M = 2, and to 50 in M = 20), and is controlled by particles with γ ∼ γ LC rad (hence the similarity of the red and green curves). When we decrease γ LC rad to 200 for the σ LC ∼ 50 run (green dashed line), the peak is shifted towards lower energies (despite the fact that particle distributions are only marginally affected).

Summary
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-micro scales 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-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. 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 2 . Photon spectra, on the other hand, are evidently different. In the marginal cooling regime, γ LC rad ∼ σ LC , both the peak and the high-energy cutoff of the emission are controlled by the magnetization parameter, E p ≈ E cut ≈ ω LC B (few · σ LC ) 2 . In the magnetospheres where synchrotron cooling is strong, γ LC rad < σ LC , plasmoids that carry the bulk of the plasma and contribute to the high energy emission the most quickly cool down to the characteristic temperatures of T e ∼ γ LC rad 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 :

Observational Implications
Observations by the Fermi telescope during the last decade have shown that some of the pulsars withĖ 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 non-linear stage of the tearing instability which develops in the equatorial current sheet on microscopic timescales, ω −1 t ≈ Γw cs /c, where w cs is the characteristic width of the sheet, and Γ ≈ R LC /R * is the bulk Lorentz factor of the upstream plasma in the wind. 12 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 γ ≈ (T e /m e c 2 ) ≈ γ LC rad : w cs ≈ γ LC rad m e c 2 /(|e|B LC ), where the γ LC rad 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, Ω −1 . The ratio of these timescales is always large for γ-ray pulsars, ω t /Ω ≈ 10 3Ė 36 / √ B 12 , meaning that at least for the young pulsars withĖ 10 35 erg/s the reconnection develops almost instantly in the current sheet near the Y-point. 13 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, √ σ LC Γ ≈ R LC /R * . 14 In young pulsars Γ ≈ 20 B , implying that the two are separated by at least an order of magnitude. However, even when √ σ LC ∼ Γ (which might be the case for the older pulsars witḣ E ≈ 10 34 erg/s) our findings indicate that the rate is only marginally slower (by a factor of two 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 2 , and the Larmor radii of the most energetic particles, ρ L , are limited to microscopic plasma scales, R LC /ρ L ≈ 2 · 10 4 K 4 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 13 Here and further in this section we normalize all the quantities to fiducial values ofĖ 36 ≡Ė/(10 36 erg/s), B 12 ≡ B * /(10 12 G), and K 4 ≡ κ/10 4 (plasma multiplicity w.r.t. the local Goldreich-Julian value). For reference,Ė varies between 10 32 and 10 38 erg/s (Vela: 7 × 10 36 erg/s, Crab: 4 × 10 38 erg/s), B * weakly changes between 10 12 -10 13 G, and κ ≈ 10 2 -10 7 (with 10 4 being the case only for a few youngest pulsars, e.g., the Crab). In practice κ is determined by the dynamics of the polar cap discharge, and the pair-production efficiency near the light cylinder, both of which indirectly depend onĖ and B * (or P andṖ ). But for the purposes of this discussion we omit this dependency. 14 The value of Γ, which is the Lorentz factor of the bulk motion of pairs flowing from the polar cap to the outer magnetosphere, is set by the last curvature photons to pair produce before escaping the polar cap discharge region. The energy of pairs produced by these photons is determined by the angle, ψ, between the photon momentum and the magnetic field: γ ± ≈ 1/ψ (this is valid in the limit where the energies of photons are large: mec 2 , and the angle ψ 1; see, e.g., Timokhin 2010). Since the discharge operates up to ≈ R * from the surface of the star, the last photons to pair produce will have ψ ≈ R * /Rc, where Rc is the radius of curvature of the magnetic field lines Rc ≈ R * R LC /R * . From these we finally find: Γ = γ (last) ± ≈ R LC /R * . is emitted directly in the form of the highest-energy synchrotron photons around one-to-few GeV, leading to characteristically high values of L γ /Ė ≈ 1%-10%. The youngest Crab-like pulsars, on the other hand, produce 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: σ (max) γγ ≈ 0.25σ T ; see, e.g., Akhiezer & Berestetskij 1985), and n x is the number density of X-ray photons, the presence of which is required for the GeV photons to pair produce. 15 The effective number density of the low-energy photons can be estimated as n x ≈ (L x /cε x )/(4πR 2 LC ), where L x is the luminosity in X-rays, and ε x is the characteristic energy of X-ray photons. Pair production cross section is the largest for the photons that satisfy ε x ε γ ≈ (m e c 2 ) 2 (here, ε γ ≈ GeV). Substituting all the values, and taking L x /Ė ∼ 1%, which is the case for the youngest pulsars (Marelli et al. 2011), we finally find: τ γx ≈ 10 −2Ė 5/4 36 / √ B 12 . For the Crab,Ė ≈ 4 · 10 38 erg/s, and we find that τ γx ≈ 1. The relatively large value for the optical depth means that the non-negligible 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 γ sec ≈ E γ /(m e c 2 ), where E γ ≈ 0.1-1 GeV is the charactertic energy of the pair-producing photons. For the Crab, γ sec ≈ 10 3 and the synchrotron emission of these secondary pairs corresponds to energies E sec ≈ γ 2 sec |e|B LC /(m e c) ≈ 1 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 cut ≈ (3/2) γ 2 cut |e|B ⊥ /(m e c), 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 highestenergy pairs. The dissipation power is equal to |e|E rec c, where E rec ≈ β rec B LC , while the synchrotron radiation power is σ T cγ 2 cutB 2 ⊥ /(4π). From the definition of 15 See also similar arguments in the context of the M87* by Ripperda et al. 2022. γ LC rad we may then find: γ cutB⊥ ≈ γ LC rad B LC , where B LC is the upstream magnetic field strength. Substituting this into the expression for the cutoff energy yields: The cutoff energy is then uniquely set by the ratio γ cut /γ LC rad , resulting in E cut ≈ (200 MeV)Ė 7/8 36 B −1/4 12 /K 4 . For Vela-like pulsars, κ ≈ 10 3 -10 4 (Timokhin & Harding 2015),Ė ≈ 10 36 -10 37 erg/s, resulting in E cut ≈ 1-10 GeV. 16 For the most energetic Crab-like pulsars with a strong magnetic field near the light cylinder (B LC ∼ 4 · 10 6 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 GeV photons. The production rate can be estimated asṄ ± ≈ τ γx L γ /ε γ , where ε γ ≈ GeV. To obtain the multiplicity κ we compare this with the GJ flux,Ṅ GJ ≈ cĖ/|e|, where we assumed that theṄ GJ e is equal to the GJ current: cĖ. We then find the multiplicity of pair-production near the light cylinder For the Crab pulsar, L x /Ė ≈ 1%, and L γ /Ė ≈ 0.1% yielding κ LC ≈ 10 6 . For the Vela, on the other hand, L x /Ė ≈ 0.01%, and L γ /Ė ≈ 1% yielding κ LC ≈ 3 (similar to what has been obtained by Lyubarskii (1996). 17 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-Ė Crab-like, and the low-Ė Vela-like pulsars. The cutoff energies are dictated by the highestenergy pairs, and in both cases the numbers are close to the observed one-to-few GeV. 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 which further re-radiate via synchrotron at a range of energies spanning all the way to keV, 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 where authors reportĖ ∝ L 1/2 γ beyondĖ 10 36 erg/s.

Future work
In our simulations we did not include self-consistent pair production, either the γB → e ± process near the polar cap, or 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, reradiating 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 gamma-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. 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 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. In this Appendix we present some dimensionless relations both in the context of astrophysical pulsars, and our simulations. We discuss how each of these parameters is up/down scaled in our global simulations relative to realistic values.
Consider a rotating perfectly conducting neutron star with a radius R * and a magnetic moment µ = B * R 3 * , where B * is the magnetic field strength near the equator. For simplicity let us consider χ = 0 • , i.e., the magnetic moment is aligned with the rotation axis. If the magnetosphere is filled with enough plasma to provide the necessary charge density for screening the parallel electric field, then magnetic field lines that originate inside a circular region of radius R pc around the magnetic axis will be open, i.e., will continue to infinity. This region, the polar cap, has a radius of R pc = R * R * /R LC , where R LC = c/Ω is the light cylinder of the neutron star. Electric field generated at the surface of the star across the polar cap generates a potential drop, V pc , which can be written in a dimensionless form where ω * B = |e|B * /m e c is the nominal gyrofrequency at the surface of the star, and P is the rotation period of the star.
Defining the plasma multiplicity close to the light cylinder, λ, as n LC = λΩB LC /2πc|e|, the scale separation (ratio of the light cylinder to the plasma skin depth) in the magnetospheric current layer is given by the following expression: V pc m e c 2 1/2 ≈ 7 × 10 6 λ 10 4 1/2 B * 10 12 G 1/2 P 100 ms −1 . (A2) The plasma magnetization near the light cylinder is The duration, resolution and other parameters of our simulations are chosen in such a way as to satisfy several requirements, while still being computationally feasible. These requirements are the following: 1. T 3P : the durations of our simulations are typically a few rotation periods to ensure the initial transient state has passed and a steady-state is reached; 2. L 3 (5R LC ) 3 : size of the domain has to be large enough to fit at least 1R LC of the current sheet; 3. R * ∆x: stellar surface needs to be well resolved; 4. V pc /m e c 2 1: potential drop at the polar cap needs to be very large; 5. λ > 1: to ensure enough plasma is provided to screen the accelerating electric field, the plasma multiplicity has to be at least a few; 6. d * e ∼ ∆x: the plasma skin depth at the surface of the star has to be marginally resolved (or at least not strongly underresolved); 7. d LC e ∆x: the plasma skin depth near the magnetospheric current sheet has to be resolved to properly capture the reconnection process; 8. R LC /d LC e 100: the scale separation in the current sheet has to be large; 9. σ LC 1: to ensure the reconnection is in high-σ regime (as shown in Eq. (A2) this is guaranteed from the first and second conditions); 10. ω LC B ∆t −1 : to properly capture the synchrotron cooling and gyration of particles in the current layer, their corresponding gyration periods have to be resolved.

B. GUIDING CENTER APPROXIMATION
In all of our simulations of pulsar magnetospheres we employ a hybrid particle mover algorithm (Bacchini et al. 2020). The algorithm allows to switch between two numerical schemes for solving the equation of motion of particles depending on the electric and magnetic field values at particle location, as well as the momentum of the particle. In the GCA (guiding center approximation) regime the motion of the particle is reduced to the motion of its guiding center, ignoring the curvature and ∇B drift terms which are much smaller than the E × B term. The equation of motion in this regime can be written in the following form: Here v E is the three-velocity of the frame where E and B are parallel, with w E /c = E ×B/(E 2 +B 2 ). Here we implicitly decomposed the velocity of particle into three different components: where u ⊥g is the particle velocity component responsible for its gyration. Notice, that we enforce the magnetic moment of the particle to be zero in the GCA regime. This assumption has physical justification: in pulsar magnetospheres bulk of the particles essentially reside at the zeroth Landau level, radiating their perpendicular momentum via synchrotron mechanism almost instantly; setting µ = 0 when particle switches from the conventional algorithm to GCA is similar to that process. GCA algorithm allows to significantly underresolve cold particle gyration close to the stellar surface as well as in the most of the magnetosphere, while still being able to model the strong synchrotron cooling for the high-energy particles treated with the conventional pusher. In order to properly capture the regions where important kinetic plasma physics occurs (such as reconnection in the current layer), while still retaining the GCA approach for the bulk of the particles in the magnetosphere, we employ simple criteria for switching between the conventional Boris and GCA. If either of the following two conditions is satisfied, the full equation of motion of the particle is solved with the conventional algorithm: Here f E and f B are dimensionless numbers which we typically choose to be 1 and 0.95 correspondingly. The first condition ensures that only particles with unresolved gyroradii are reduced to their corresponding guiding centers. The second condition ensures that the GCA is not applied in the regions of vanishing magnetic field where particle gyration is ill-defined.

C. DEPENDENCE OF THE RECONNECTION RATE ON THE UPSTREAM MOTION
To test how the rate of magnetic reconnection depends on the bulk velocity of upstream plasma along the layer we perform 2D localized simulations. We start with a Harris equilibrium defined in the frame moving with a Lorentz factor Γ along the current layer with respect to the lab frame. In the first set of simulations we keep σ up /Γ = B 2 /(4πρ up Γ) = 20 constant (ρ up is the upstream density in the lab frame) and vary Γ (meaning that we also vary σ up ). In the second set, we keep σ up = 100 (as measured in the lab frame), while, again, varying Γ. After about a light crossing time of the box, when the reconnection proceeds in the plasmoid-dominated stage (but before it shuts down because of the finite size of the box) we can measure the effective reconnection rate in two ways. We can either directly measure the upstream drift velocity towards the current layer, β in = (E × B) in /B 2 , or we can compute how fast the magnetic energy is dissipated by evaluating the slope of the dE B /dt curve, where E B =´V d 3 rB 2 /8π is the total magnetic energy contained in the box. Both of these methods unsurprisingly provide the same answer, so below we only report the direct measurement of β rec .
In Figure 11 (left) we plot our measurements for the inflow velocity, β in = β rec as a function of the bulk Lorentz factor of the upstream, Γ, for the two simulation series mentioned above. Red curve corresponds to the case when σ up /Γ = 20 is kept constant; in this series we see almost no variation in the inflow velocity, the value of which lies between 0.8 and 0.14. Notice, that keeping σ up /Γ fixed means that we scale σ up (in the lab frame) up with increasing Γ. In the case when σ up is kept fixed we see a moderate decrease in β rec .
To understand this effect, let us consider separately the velocities of relativistic outflows from the x-point (as shown in the right panel of Figure 11). In the frame co-moving with the upstream plasma the velocities of these outflows should be symmetric w.r.t. the y-axis. Or in other words, u ± ≈ u A = √ σ , where dashed quantities are measured in the frame co-moving with upstream plasma. Here σ ≡ (B ) 2 /(4πρ e c 2 ); "±" indicate the direction of motion along and opposite to the current layer (the y-axis). Boosting back to the lab frame we find: u + ≈ u + Γ, and u − ≈ u − /Γ, assuming the "upstream" frame moves along the y-axis with a Lorentz factor of Γ 1. Cold upstream magnetization (not including bulk Lorentz factor) in the lab frame, σ ≡ B 2 /(4πρ e c 2 ), can be expressed as σ = σ /Γ, since B up = B y does not transform, and ρ e = Γρ e . Then the outflow velocities in the lab frame can be written as: u + ≈ u A Γ 3/2 , and u − ≈ u A Γ −1/2 , where u A = √ σ is the relativistic Alfvén velocity in the lab frame. This can be also seen in the right panel of Figure 11: outflow velocities along the y-axis measured in the lab frame are clearly asymmetric. When Γ 1/2 becomes comparable to u A , the outflow velocity in the direction opposite to the upstream boost, u − , may become non-relativistic. In this case it is natural to expect the reconnection rate, defined in our case as β rec ≡ v in /v out ≈ v in /c to drop, since the effective outflow three-velocity, v out , can now be smaller than c.
Arguments provided in this section are not exhaustive, and should rather be taken as possible explanations of empirical facts, observed in numerous simulations. Further investigation is definitely necessary to better understand the nature of the mechanism that controls the reconnection rate in these cases.

D. ENERGY DISTRIBUTION OF PARTICLES IN THE DIFFERENT PARTS OF THE SIMULATION DOMAIN
For the simulation R75 ang0 magnetization parameter in the upstream, shown in Figure 6, is close to σ LC ∼ 10 3 (as shown with the line-out plot on the right) and drops to zero inside the current sheet, where the magnetic field vanishes. This means that the characteristic energies to which particles can get accelerated in the current sheet are comparable to σm e c 2 (as will be demonstrated shortly).
Let us consider distribution functions for e ± in different regions of our simulation R75 ang0, as shown in Figure 12. On a poloidal slice in panel 12a we show the bulk Lorentz factor of the plasma, Γ = 1 + U 2 /c 2 , which is computed Initially the upstream plasma is pushed in the +y direction. In the current sheet, nevertheless, relativistic outflows move both directions, with the x-points in-between being static.
using the bulk four-velocity, U =´uf (u)d 3 r. Panels 12b, 12c, and 12d show distribution functions of the electrons and positrons in three different regions: in the separatrix (last closed field line), in the current sheet, and in the upstream correspondingly. Upstream plasma is relatively cold; electrons and positron flow outwards along the magnetic field lines, gaining characteristic bulk Lorentz factors of Γ ∼ 10 ( Figure 12d). On the separatrix (Figure 12b) both electrons and positrons from the surface are marginally accelerated by an unscreened electric field, gaining energies close to γ ∼ 10-100. This region also hosts hot electrons returning from the Y-point, which is evident from the excess of electrons at Lorentz factors of a few 10 2 shown in Figure 12b. These two regions, the upstream and the separatrix, act as an "intake" and an "exhaust" for the current sheet. It is important to note that particles in these regions are generally exposed to a substantial orthogonal magnetic field component. This means that high-energy particles from the current sheet cannot exist there for timescales longer than the short cooling time. Thus, the main role in shaping the observed γ-ray emission is played by the current sheet, where constant release of magnetic energy sustainably accelerates particles in x-lines, where the magnetic field strength is zero (the cooling is non-existent).
In the current sheet (Figure 12c) we see a substantial non-thermal particle population, extending to γ ∼ σ ∼ 500-1000. The bulk motion of the current sheet, on the other hand, has a Lorentz factor of Γ ∼ 10-100, consistent with characteristic velocities of relativistic flows along the current sheet during the reconnection, Γ ∼ √ σ (which corresponds to the Alfvén speed). 19 Notice also, that the distributions of electrons and positrons are slightly different: positrons extend to slightly higher energies. This is due to the fact that the electric field generated during reconnection "pushes" the electrons opposite to the global E × B drift, which is charge-invariant and directed radially outward (in Figure 3a the reconnection electric field is pointing in the out-of-plane direction and has both toroidal and radial bulk Γ factor of e ± e ± distributions e + + e − e + e − Figure 12. Distribution functions of electrons and positrons measured in different locations of our simulation shown on an azimuthally averaged poloidal slice. Synchrotron cooling in this particular simulation (R75 ang0 gr2) is weak. Bulk Lorentz factor is shown in panel (a). Panels (b), (c), and (d) show distribution functions in the separatrix, the current layer, and the upstream, respectively. Plasma in the upstream is typically cold with a bulk outflow Lorentz factor of a few. In the current sheet the bulk motion of plasma is dictated by the dynamics of reconnection, and the characteristic outflow velocities are comparable to the Alfvén speed with a Lorentz factor ∼ √ σ. As expected, the highest energy particles are produced in the current layer where magnetic reconnection occurs. components). This difference in the electrons and positrons has been observed in earlier simulations (see, e.g., Cerutti et al. 2015); however, it almost vanishes for large obliquity angles, χ, which we also see in our simulations.