11email: adrien.soudais@univ-grenoble-alpes.fr 22institutetext: Research Center for Astronomy and Applied Mathematics, Academy of Athens, GR 11527 Athens, Greece
Scaling up global kinetic models of pulsar magnetospheres using a hybrid force-free-PIC numerical approach
Abstract
Context. The particle-in-cell approach has proven effective at modeling neutron star and black hole magnetospheres from first principles, but global simulations are plagued with an unrealistically small separation between the scales where microphysics operates and the system-size scales due to limited numerical resources. A legitimate concern is whether the scale separation currently achieved is large enough, such that results can be safely extrapolated to realistic scales.
Aims. In this work, our aim is to explore the effect of scaling physical parameters up, and to check whether salient features uncovered by pure kinetic models at smaller scales are still valid, with a special emphasis on particle acceleration and high-energy radiation emitted beyond the light cylinder.
Methods. To reach this objective, we develop a new hybrid numerical scheme coupling the ideal force-free and the particle-in-cell methods, to optimize the numerical cost of global models. We propose a domain decomposition of the magnetosphere based on the magnetic field topology using the flux function. The force-free model is enforced along open field lines while the particle-in-cell model is restricted to the reconnecting field line region.
Results. As a proof of concept, this new hybrid model is applied to simulate a weak millisecond pulsar magnetosphere with realistic scales using high-resolution axisymmetric simulations. Magnetospheric features reported by previous kinetic models are recovered, and strong synchrotron radiation above 100MeV consistent with the Fermi-LAT gamma-ray pulsar population is successfully reproduced.
Conclusions. This work further consolidates the shining reconnecting current sheet scenario as the origin of the gamma-ray emission in pulsars, as well as firmly establishes pulsar magnetospheres as at least TeV particle accelerators.
Key Words.:
acceleration of particles – magnetic reconnection – radiation mechanisms: non-thermal – methods: numerical – pulsars: general – stars: winds, outflows1 Introduction
Magnetospheres forming around neutron stars and black holes, or simply referred to as relativistic magnetospheres in the following, are involved in some of the most extreme high-energy astrophysical phenomena, such as pulsars (Cerutti & Beloborodov, 2017; Philippov & Kramer, 2022), magnetars flares and outbursts (Thompson & Duncan, 1995; Kaspi & Beloborodov, 2017), relativistic jets from supermassive black holes (Blandford & Znajek, 1977; Blandford et al., 2019), and possibly fast radio bursts (Popov & Postnov, 2010; Lyubarsky, 2020) and precursor emission to compact object mergers (Crinquand et al., 2019; Most & Philippov, 2020). Elaborating a quantitative and predictive model of relativistic magnetospheres turned out to be an outstanding endeavor at the forefront of modern physics, and our current understanding of such systems is at best incomplete. The fact that these regions concentrate large energy densities into various forms (i.e., gravitational, rotational, electromagnetic), and the interplay between plasma physics, general relativity and relativistic (quantum) electrodynamics dramatically increases the degree of complexity of the problem that becomes analytically untractable.
In the quest for a comprehensive theory of relativistic magnetospheres, the aligned rotator has historically played the role of a Rosetta stone in this field. Some of the most important developments were obtained in the framework of force-free electrodynamics, a degenerate form of relativistic magnetohydrodynamic (MHD) in which the evolution of the plasma is solely governed by the Lorentz force (Blandford, 2002; Komissarov, 2002). Thus, all of the other forces including inertia or gas pressure are irrelevant in this regime. This is a valid assumption almost everywhere within the magnetosphere where the field strength is high and the plasma density is low. In this limit, the structure of the magnetosphere is then a solution of the relativistic Grad-Shafranov equation (Scharlemann & Wagoner, 1973; Michel, 1973), and an exact solution can be derived for a magnetic monopole (Michel, 1973; Blandford & Znajek, 1977). Extending these results to a more complex magnetic field topology, even to an aligned magnetic dipole, has never been achieved by analytical means.
Instead, an approximate numerical model must be used. The force-free solution of the aligned dipole recovers essential magnetospheric features (Contopoulos et al., 1999; Komissarov, 2006; Spitkovsky, 2006; Timokhin, 2006; McKinney, 2006; Parfrey et al., 2012; Cao et al., 2016), in particular the overall hybrid magnetic topology composed of closed and open field lines and the formation of an equatorial current layer separating both magnetic hemispheres beyond the light-cylinder radius, where co-rotation with the star is superluminal (Goldreich & Julian, 1969; Michel & Tucker, 1969). Thereby, the force-free framework provides a realistic description for the morphology of relativistic magnetospheres in the limit where plasma is abundant and ultramagnetized everywhere. It also captures well how the star spins down under the effect of magnetic torques. However, the main shortcoming of this approach is that it cannot capture dissipation (see, however Gruzinov 2008 and Li et al. 2012 for a dissipative force-free formulation), and more importantly non-thermal particle acceleration and radiation. It is therefore impossible to pin down the emitting regions and to constrain the model with observations without using ad-hoc prescriptions, leading to large theoretical uncertainties (Bai & Spitkovsky, 2010; Kalapotharakos et al., 2014; Cao & Yang, 2019).
Observations of rotation-powered pulsars in the gamma-ray range demonstrate that the magnetosphere accelerates particles with a disconcerting efficiency (Abdo et al., 2010, 2013; Smith et al., 2023). The recent detections of TeV pulsed emission in the Crab and Vela pulsars indicate that electrons are accelerated at least up to 20TeV (Ansoldi et al., 2016; H. E. S. S. Collaboration et al., 2023), which represents a sizeable fraction of the full polar-cap potential drop. Acceleration of hadrons in the magnetosphere is also an open issue of great interest for understanding the origin of cosmic rays (e.g., Gunn & Ostriker 1969; Blasi et al. 2000; Arons 2003; Fang et al. 2012). To address these outstanding questions and to connect the physical model to observations, an ab-initio plasma model is required, taking into account basic plasma processes at microscopic (i.e., kinetic) scales, pair production, the radiation-reaction force, and general relativistic effects. In this respect, the particle-in-cell (PIC) approach has proven effective at integrating all of the above physical ingredients necessary to model relativistic magnetospheres from first principles (Philippov & Spitkovsky, 2014; Chen & Beloborodov, 2014; Cerutti et al., 2015; Belyaev, 2015; Philippov et al., 2015b; Guépin et al., 2020; Chen et al., 2020; Hu & Beloborodov, 2022; Bransgrove et al., 2023; Cruz et al., 2023; Hakobyan et al., 2023; Torres et al., 2024). Global PIC simulations of pulsar magnetospheres show robust features such as a sizeable dissipation rate of the outflowing Poynting flux via magnetic reconnection of open field lines beyond the light cylinder. This energy is efficiently channelled into non-thermal particle acceleration and energetic radiation, possibly explaining the observed high-energy radiation (Cerutti et al., 2016; Philippov & Spitkovsky, 2018; Kalapotharakos et al., 2018). At the base of the open field lines near the star, the spark-gap dynamics and the putative radio emission are also well reproduced by local PIC models (Timokhin & Arons, 2013; Philippov et al., 2020; Cruz et al., 2021).
In spite of these impressive developments, these results must be interpreted with caution because global PIC simulations are plagued with an unrealistically small separation between the scales where microphysics operates, and the system-size (i.e., light-cylinder) scales. A legitimate concern is whether the scale separation achieved by current global PIC simulations is large enough such that results can be safely extrapolated to realistic scales. Thus, capturing the polar-cap physics and the light-cylinder reconnection physics with full scales in the same simulation seems impossible, even in a far future. In this work, we make the choice to sacrifice the polar-cap microphysics to focus our numerical resources on the light-cylinder scales, where most of the dissipation and particle acceleration take place (Cerutti et al., 2020; Hakobyan et al., 2023). Our main objective is to scale physical parameters up, and to check whether salient features uncovered by pure kinetic models at smaller scales still hold strong. To this end, we develop a new hybrid numerical scheme to take advantage of the complementarity between the force-free and the PIC approaches. We propose a domain decomposition based on the magnetic field topology of the system, where open field lines are described by an ideal force-free model while the reconnecting region is described by a PIC model. As a proof of concept, we apply this method to the aligned pulsar magnetosphere with the aim at reaching realistic scales of a weak millisecond pulsar, just above the limit of the observed Fermi-LAT pulsar population. In the next section (Sect. 2), we present the relevant physical scales needed to model a pulsar magnetosphere. In Section 3, we describe the new hybrid force-free-PIC numerical scheme. This approach is first tested against the monopole solution (Sect. 4). It is then applied to the aligned dipole with high-resolution simulations, where we explore the effect of scaling up physical parameters, from previously achieved scales in full PIC simulations, to more realistic scales (Sect. 5). We summarize this work and we discuss future applications and developments of this new hybrid approach in Section 6.
2 Scale separation
Capturing all relevant physical scales in pulsar magnetospheres is an outstanding challenge for simulations. In the following discussion, we will consider, as our baseline, a millisecond pulsar with a rotation period ms and a surface dipolar magnetic field G, corresponding to the weakest observable gamma-ray pulsars of spindown power erg/s as reported by the Fermi-LAT (Abdo et al., 2010, 2013; Smith et al., 2023). It also represents a realistic gamma-ray pulsar with the smallest scale separation that we strive to capture in this work as a proof of principles (see Sect. 5).
The system-size macroscopic scales are set by the neutron star radius, cm, and the light-cylinder radius,
(1) |
The light-cylinder radius sets the scale at which the plasma streams away from the star along open fields in the form of a relativistic wind. Magnetic reconnection starts operating at the base of the wind current layer, known as the “Y-point” (e.g., Uzdensky 2003; Contopoulos et al. 2024), and proceeds at all radii beyond this point. Energetic synchrotron emission is produced within a few light-cylinder radii because of the steep decay of the magnetic field strength with radius (Cerutti et al. 2016; Cerutti et al. in preparation), and thus sets the largest relevant scale of the problem in this work.
While the macroscopic scales are well constrained, defining the microscopic scales is more involved because plasma physics and QED effects must be taken into consideration. The electronic plasma skindepth at the star surface sets the minimal relevant scale of the problem. Assuming that the plasma density is a multiple of the critical co-rotation density, (Hones & Bergeson, 1965; Goldreich & Julian, 1969), where is the plasma multiplicity, and is the electron charge, yields
(2) | |||||
where is the particle Lorentz factor and is the electron mass, so that for our reference pulsar case. The situation may still look hopeless at this stage. However, particles at the neutron star surface are not at rest. Instead, they are believed to move at ultra-relativistic speeds (), hereby reducing the above scale separation. A primary beam of particles is extracted and accelerated by the surface electric field induced by the rotation of the star (Sturrock, 1971; Ruderman & Sutherland, 1975). At best, the particles will undergo the full vacuum potential drop across the polar cap of the star. The size of the polar cap is defined between the magnetic axis and the first field line fully contained within the light cylinder. For an aligned dipole in vacuum, its angular size is . Using the corotation surface electric field, , we can derive the potential drop across the polar cap as
(3) |
where is the angular velocity of the star, and its magnetic moment. An electron experiencing the full potential drop will acquire a Lorentz factor given by
(4) |
This calculation provides an upper limit, it is not a realistic estimate of the particle Lorentz factor in the bulk of the flow, because the presence of the plasma screens in great part the accelerating electric field (parallel to the magnetic field), and because of pair production. The main channel for pair production near the star surface is believed to be by magnetic conversion, meaning that a gamma-ray photon is annihilated by virtual photons from the magnetic field to produce a pair. Quantum electrodynamics predicts that pair production occurs if (Erber, 1966; Harding & Lai, 2006)
(5) |
where is the dimensionless photon energy, and is the effective perpendicular magnetic field (see Eq. 17 below) normalized to the critical quantum field G. If electrons radiate via synchrotron-curvature radiation as usually assumed, the critical photon energy is given by (Blumenthal & Gould, 1970)
(6) |
Combining Eq. (5) with Eq. (6) provides an estimate for the threshold electron Lorentz factor capable of producing new pairs,
(7) |
while the created pair shares the absorbed photon momentum, such that the Lorentz factor of the secondary pairs may be estimated as
(8) |
Assuming that pair creation is effective at producing a high-multiplicity plasma in the magnetosphere (), it sets the relevant scale of the particle Lorentz factor to be considered in estimating the plasma skindepth scale in Eq. (2), yielding for our reference pulsar. One can realize from these order of magnitude calculations that a weaker magnetic field (G) would prevent pair formation as the threshold energy scale would be larger than the full potential drop.
At the light cylinder, a thin current sheet forms whose thickness, , is governed by the local plasma skindepth scale, meaning that . At this stage, efficient particle acceleration proceeds via magnetic reconnection, increasing on average the particle Lorentz factor to a scale governed by the plasma magnetization parameter,
(9) |
where , is the Goldreich-Julian plasma density at the light cylinder, and is the pulsar wind bulk Lorentz factor at its base. Assuming that these energetic particles set the local skindepth scale yields
(10) |
(11) |
thus, or for these fiducial parameters. A similar estimate can be obtained from Ampère law (Coroniti, 1990; Lyubarsky & Kirk, 2001; Cerutti & Philippov, 2017) (see also, Arka & Dubus 2013; Uzdensky & Spitkovsky 2014, for other estimates).
Radiative cooling sets another energy scale in the problem at which the accelerating electric force is balanced by the radiation-reaction force, . The latter effectively acts as a continuous drag force (in the classical regime, i.e., ), opposite to the particle direction of motion, whose magnitude in the ultrarelativistic regime () can be estimated as , where is the classical radius of the electron. The fiducial radiation-reaction-limited electron Lorentz factor is then given by
(12) |
or at the light cylinder for the fiducial pulsar parameters, assuming and .
In summary, modeling the weakest gamma-ray pulsar must satisfy at the very least the following scale separation, in terms of spatial quantities normalized to ,
(13) |
and in terms of energy scales normalized to
(14) |
This exercise reveals that in such system , and are very close together meaning that the energy range of the particle spectrum may be quite narrow. This property is in agreement with a new important feature highlighted by the most recent analysis of the Fermi-LAT pulsars, which shows that the observed spectral energy distributions is narrower with decreasing spindown, nearly reaching a spectrum consistent with a monoenergetic particle population for the weakest pulsars of spindown erg/s (Smith et al., 2023). Hereby, even though the estimates presented in this section represent by no means a rigorous calculation, it provides an approximate assessment of the minimum computational needs to approach a realistic systems using a kinetic approach. Our conclusion is that such weak pulsar could in principle be captured with full scales with current computational facilities, provided that extreme resolutions are used in a 2D axisymmetric system, as for instance in the recent works by Hu & Beloborodov (2022); Bransgrove et al. (2023).
3 Hybrid force-free-PIC approach
In order to fulfil the minimum scale separation requirements highlighted in the previous section, a large number of grid cells must be used, and hence a large computational power is needed. To alleviate some of this numerical cost, we develop a new hybrid approach specially designed for aligned magnetospheres, coupling time-dependent force-free electrodynamic and the PIC approaches in the same simulation. This idea relies on the fact that gamma-ray pulsar magnetospheres are mostly very close to a dissipationless force-free state because of efficient pair creation (i.e., ). Using the expensive full PIC technique with a large number of particles per cell everywhere in the magnetosphere may appear as excessive, while a numerically cheaper force-free description would be sufficient. In contrast, dissipative regions (e.g., electrostatic gaps, current sheets) should remain on microscopic scales (again, if ), and so are the scales at which the PIC approach is really needed. Therefore, it appears desirable to combine these two techniques, but a new scheme must be implemented to ensure that both approaches can coexist harmoniously together.
In this work, we develop a new force-free module implemented in the Zeltron code, an explicit relativistic PIC code (Cerutti et al., 2013) used in the past for global full PIC simulations of pulsar magnetospheres (e.g., Cerutti et al. 2015). In the PIC approach, the plasma is modeled by point-like superparticles –or simply referred as particles in the following– representing a very large number of physical particles with identical mass-to-charge ratio following the same path in phase space (Birdsall & Langdon, 1991). Their evolution is governed by the Lorentz force, and the radiation-reaction force in this pulsar-specific context, as
(15) |
where is the dimensionless particle 4-momentum vector, and is the particle electric charge, and where
(16) |
is the radiation reaction force following the Landau & Lifshitz (1971) formulation, and , and where
(17) |
is the effective perpendicular magnetic field (Cerutti et al., 2016). Particle trajectories are evolved in time with a modified Boris pusher to incorporate the radiation reaction force to the Lorentz force. The current density carried by each particle is deposited onto the numerical grid where the electromagnetic fields are defined. Summing over all particles and over all species, the total current reconstructed on the grid is then used to evolve the time-dependent Maxwell equations,
(18) |
(19) |
allowing to close the PIC loop performed at each time step (Fig. 1). The field integration step follows a standard finite-difference time-domain method that involves a staggered mesh (Yee, 1966).
In general, each simulation cell is filled with a large number of particles, and thus pushing the particles and depositing the currents largely dominate the overall computing time. These two steps typically take one order of magnitude longer than evolving the fields alone. This is where the benefits of force-free electrodynamics come into play. The force-free condition is given by
(20) |
where is the electromagnetic tensor and is the 4-current, which translates into the following conditions
(21) |
using the spatial components, where is the charge density, and
(22) |
using the temporal component. Eq. (21) also leads to the condition
(23) |
Combining Eqs. (21)-(23) with Maxwell equations Eqs. (18)-(19), one can express the total current density solely from the fields in the form of an effective Ohm’s law, as (Gruzinov, 1999; Blandford, 2002)
(24) |
which can be understood as the sum of a current flowing across field lines, , and a component flowing along field lines, . Instead of pushing a large number of particles and depositing the current on the grid, the current is derived in a single step using Eq. (24), which saves on computing time and memory. Other than this, the Maxwell solver is unchanged whether the current is computed via the particles or via the fields. Thus, the coupling between both approaches is done via the current.
To implement the force-free module, we follow a similar numerical procedure as in Spitkovsky (2006), which consists in solely computing the perpendicular force-free current, . The parallel component, , prevents the formation of an electric field component along the magnetic field at all times, , but computing is a numerically delicate and cumbersome endeavor on a staggered mesh (see, however, Parfrey et al. 2012). Therefore, the electric field is renormalized at each time step to enforce the force-free condition,
(25) |
The code also checks at every time step that the condition is always satisfied, so that if it is not the case, the electric field is once again renormalized as
(26) |
Computing the force-free current and renormalizing the electric field directly on the staggered mesh is delicate, because field components defined at different locations on the cell and times are mixed, leading to a disruption of the second-order accuracy of the scheme. Instead, we interpolate all quantities into the same grid point to compute the current and the electric field, and then interpolate these back onto the staggered grid. The interpolation scheme is based on an arithmetic average between neighboring points.
To blend the two approaches, the computing box is decomposed into pure force-free and pure PIC domains. In axisymmetric systems, using the magnetic flux function is a natural choice to make this division. This scalar field is constant along poloidal field lines which are themselves equipotential surfaces, it is therefore a physically motivated criterion which is also numerically convenient to tag individual field lines. In ordinary spherical coordinates (), the magnetic flux function is given by
(27) |
Once the boundary is set, a thin transition layer is implemented between neighbouring field lines where the currents are blended to ensure a better continuity between the PIC and the force-free currents and to minimize numerical artefacts, such as spurious reflections or accumulation of charges. A simple linear interpolation scheme is sufficient for this purpose, so that the full hybrid current is expressed as
(28) |
where is the current density reconstructed from the particles, is the current from the force-free condition (Eq. 24), and is the blending function. This hybrid current is then injected into Maxwell-Ampère law (Eq. 19) to evolve the field in the whole simulation domain. The formulation of Eq. (28) implies that both schemes are computing the current in the transition layer (i.e., particles are present everywhere within this layer, but they are of course absent in the pure force-free domain). Considering a single transition layer in the domain, two distinct fluxes must be chosen, , so that the blending function is
(29) | |||||
(30) | |||||
(31) |
The magnetic flux function is computed at each time step so that the blending function is updated accordingly. In the context of pulsar magnetospheres where the magnetic field lines are frozen into the neutron star surface, this decomposition is fixed on the neutron star surface and then follows the evolution of the magnetic field topology in time in the rest of the domain. In this sense, the domain decomposition proposed in this work is not static but dynamically evolves along with the magnetospheric activity. Figure 1 graphically summarizes the proposed numerical scheme over a time iteration.
4 Monopole
The hybrid scheme is first put to the test on the modeling of a steady-state force-free monopolar magnetosphere (i.e., a single magnetic monopole). Although unphysical, studying this configuration is of great interest: it admits an exact analytical solution and it is free of dissipation –there is no gap or current sheet (as opposed to a split-monopole)– so that spurious numerical effects in the scheme can be quantified, and it represents a good model for the asymptotic magnetospheric structure of pulsar open field lines. All things considered, it sets the best conditions to validate our new approach.
In this setup, the magnetic flux function is given by
(32) |
where . An exact solution to the Grad-Shafranov equation yields the toroidal magnetic field component, (Michel, 1973)
(33) |
if in the northern hemisphere (and in the southern hemisphere). In this solution, the electric field is purely along and matches the toroidal magnetic field, . The electric current is purely radial and is given by
(34) |
where corresponds to the fiducial Goldreich-Julian current, and the spindown power associated to the electromagnetic torque is
(35) |
The goal of this first numerical experiment is to recover these properties using our hybrid scheme. To this end, we perform a set of hybrid simulations in 2D axisymmetric spherical coordinates and explore the effect of resolution on the results. The domain size explored here is composed of , and cells. The grid spacing is logarithmic along the radial direction and constant along the direction. The domain size ranges from the star surface, , up to , where , and covers the full range in . The magnetosphere is initialized in vacuum with a pure radial magnetic field, . The surface magnetic field strength is fixed at G, which is sufficiently intense in this experiment to achieve a quasi force-free state in the PIC domain. The role of the field strength and scale separation is thoroughly investigated in the next section.
Corotation of the field lines with the star at the inner boundary is guaranteed by fixing the poloidal electric field to the ideal solution . This procedure sets the magnetosphere into rotation. An absorbing layer starting at progressively damps all outgoing electromagnetic waves and removes the particles, if any (Cerutti et al., 2015). In this experiment, the northern hemisphere is modeled in full PIC, while the southern hemisphere is captured with the force-free model. This choice is well suited for a comparison between both approaches thanks to the expected symmetry between both hemispheres, and hence, any asymmetry or imbalance can be tracked easily. In this setup, we did not notice any difference with or without a transition layer for the two highest resolutions, so that here, we only show the runs performed with a sharp transition at , or . In the PIC domain, electron-positron pairs are continuously injected with a high multiplicity, , at the star surface to ensure a force-free state, and pair creation is not allowed elsewhere in the domain. The pairs are injected within the first layer of cells above the surface, and there are no ions at this stage.
Figure 2, top panel, shows the state of the simulation for the highest numerical resolution, obtained after about two spin periods. A steady state is reached and the solution is numerically stable. We recover the correct variations of the toroidal field and current densities. The latter follows the predicted cosine evolution in both hemispheres without apparent mismatch at the interface between both domains at the equator, except for the lowest resolution simulation where a small glitch is visible (middle panel in Fig. 2). This result demonstrates that both hemispheres are well balanced, meaning that there is no net accumulation of charges at the boundary or anywhere in the magnetosphere, the outgoing current in the force-free region is perfectly matched by the incoming current in the PIC region, such that there is no net current overall. The only noticeable difference in the current profile is the high level of noise in the PIC domain as opposed to the force-free region. It is in part due to the shot noise associated with the finite number of particles, but also to fast plasma oscillations. The bottom panel in Figure 2 zooms in on the radial evolution of the Poynting flux. The solution predicts that this quantity should be perfectly conserved at all radii (see, Eq. 35), so that any deviation offers a way to measure the amount of numerical dissipation as a function of resolution. The rate of dissipation is computed as
(36) |
where is the Poynting flux averaged in the interval, and is averaged over the interval (gray intervals in Fig. 2, bottom panel). We observe numerical convergence and a deviation smaller than for the run, down to for the highest resolution. We note that numerical dissipation mainly originates from the force-free hemisphere. This rate is about two orders of magnitude smaller than the amount of physical dissipation reported in the next section and in previous global kinetic model of pulsar magnetospheres over a similar radial scale (e.g., Philippov & Spitkovsky 2014; Cerutti et al. 2015; Belyaev 2015; Cerutti et al. 2020; Hakobyan et al. 2023). Thus, the hybrid scheme does not add significantly more artificial dissipation to the simulated system, at least under the ideal condition of a force-free monopole.
5 Aligned dipole
The hybrid model is now applied to the aligned rotating dipole using 2D axisymmetric simulations, aiming at scaling relevant physical scales up to the fiducial weak gamma-ray pulsars introduced in Sect. 2.
5.1 Numerical setup
|
The properties of the simulation domain are nearly identical as for the monopole configuration presented in the previous section (i.e., domain size, grid spacing and geometry, boundary conditions, pulsar spin period). Here, we focus on the highest resolution simulations that are composed of cells. As for the monopolar magnetosphere, the box is initially empty of plasma and a dipolar field fills the full domain, and corotation is enforced at the inner boundary by setting the ideal electric field. The initial magnetic flux function is
(37) |
There is no exact analytical solution in the force-free limit for an aligned dipole, and thus results are compared with previous force-free and full PIC simulations. In the PIC domain, plasma injection is done differently than in the monopolar case. First, the primary beam of particles is generated by the extraction from the crust of the star of electrons and protons with a realistic mass ratio, . This processes can be faithfully reproduced by maintaining the Goldreich-Julian plasma density at the inner boundary and at every time step. The plasma is injected in corotation with the star, but with no momentum along the poloidal field lines. The surface electric field induced by the stellar rotation then pulls and accelerates charges in the magnetosphere where pair creation can then operate. The pair production process is highly simplified, it is based on a constant energy threshold criterion of the parent electron or positron everywhere in the simulation, including in the wind zone. If the lepton Lorentz factor is larger than , then a new pair is created with a Lorentz factor at the same position and along the same direction of motion as the parent particle. This procedure is similar to the heuristic models of pair production performed in previous studies (Philippov et al., 2015b; Chen et al., 2020; Cruz et al., 2021). It provides an effective way to inject pairs in the magnetosphere with the relevant energy scale, but this model should not be viewed as very realistic. Pair production in millisecond pulsars remains a puzzle that we do not wish to address in this work. In the force-free domain, there is no other condition to consider apart from corotation at the neutron star surface.
Our choice for an efficient force-free-PIC domain decomposition in this setup is educated from the magnetic topology revealed by previous studies. For a plasma-filled magnetosphere, the main loci of dissipation are the equatorial current sheet and at its base around the Y-point region near the light cylinder, and to a lesser extent the separatrix current layers forming between the closed and open field line boundary that connect the equatorial current sheet to the star polar caps. These are also the sites for the putative gamma-ray emission in observed pulsars that we seek to reproduce in this work. Hereby, a natural choice is to limit the PIC domain to the equatorial and separatrix regions. The boundaries of the domain should be chosen with care. It is particularly important that the PIC domain extends from the star surface all the way to the absorbing boundary. If the chosen field line separating the PIC and the force-free domains reconnects inside the box, particles can no longer escape and the current layer is abruptly disrupted, leading to artificial dissipation. On the other hand, it would fail the purpose to choose a wide PIC domain. We found that a good compromise is to use the following magnetic flux function boundaries:
-
•
Force-free domain I, . This region corresponds to the open field lines connected to the star polar caps.
-
•
PIC domain, . The domain encompasses the equatorial and separatrix current layers.
-
•
Force-free domain II, . It is inside the closed field line regions deep inside the light cylinder.
In the above, is the maximum of the flux function that is located at the star surface in the equatorial plane, and is the magnetic flux crossing the light cylinder for a dipole in vacuum. This latter quantity is a good proxy for the amount of flux of open magnetic field lines even for the force-free dipole. There is a -thick transition layer between the PIC and force-free domains. However, during the transient phase, we found that the transition layer between the force-free domain I and the PIC domain can sometimes partially prevent field line opening, so that it is best to set a sharp transition in this phase. This transition layer is enforced once the initial transient has passed.
The star radius is fixed at cm, and the light-cylinder radius to cm corresponding to a ms spin period. The surface magnetic field of the reference production run is fixed at G, giving G at the light cylinder, and a force-free spindown power erg/s. The magnetic field strength then sets the relevant energy scales given in Sect. 2 as, , , and . The Larmor radius and frequency scales are not resolved near the surface of the star, but they are irrelevant since the particles stream along magnetic field lines with negligible perpendicular momentum due to the strong radiative cooling. However, the plasma skindepth and frequency scales are resolved by and at the surface where the conditions are numerically the most demanding, and they are well resolved elsewhere in the magnetosphere. The Larmor scales are also resolved in the outer parts of the magnetosphere, in particular in the equatorial current layer where these scales become relevant. The Larmor radius of secondary pairs at the light cylinder, , divided by the local grid spacing, , is , and Larmor timescale . The corresponding synchrotron cooling time is of the order of .
To prevent spurious numerical effects near the star surface due to strong cooling in the integration of particle trajectories, any particle momentum transverse to the field is removed at each time step within a -thick shell around the star. The particle pusher is unchanged everywhere else in the simulation. The radiation reaction force must however be artificially and momentarily reduced in the early stages of the simulation. The transient phase when the magnetosphere is filling up with plasma is quite violent and can generate extreme field strengths and gradients, leading to a crash in the simulation. The simulations are first evolved during one period where the strength of the radiation-reaction force is set to of its nominal value, and it is then slowly increasing in time up to about . In a last step, it is gradually increased up to the nominal strength at when a quasi-steady state can be physically exploited.
An important objective of this work is to study the effect of a rescaling of the magnetic field strength. Indeed, previous PIC simulations used unrealistically low field strength but physical quantities were rescaled to make global simulations doable. Here, we scale the magnetic field strength of the reference simulation down to G, G, and G, but keep the ratio between all energy scales the same as the G simulation, including the radiation-reaction-limited lepton Lorentz factor , given in Eq. (14). To keep the ratio identical to the G run, the strength of the radiation reaction force –or equivalently the Thomson cross section– must be boosted by a factor for G, for , and for G. Table 1 lists the numerical and physical parameters used in this study for the reference simulation.
5.2 Results
Figure 3 shows the overall magnetospheric structure achieved by the hybrid model, once the initial transient has left the box and the full scale radiation-reaction force is established at for the fiducial full-scale G millisecond pulsar. The solution is consistent with previous global PIC studies, it is characterized by prominent separatrices and reconnecting current layers whose base stands near the light cylinder (the Y-point). It is important to note that this solution has not completely reached a steady state. During the initial transient phase, the Y-point forms well inside the light cylinder (around ), and slowly migrates towards , and still continues to do so at this stage. This migration is still visible in Figure 3 in the form of density striations inside the Y-point, a feature that should vanish if a perfect steady state is achieved. The separatrix current layers are also quite thick and have not fully relaxed to the expected thin rim along the last closed field lines, and thus the final state may depart from what is shown here. Nevertheless, the magnetospheric structure in the wind zone where the equatorial current forms, which is the region of prime interest for what follows, has reached the desired state. The reconnecting layer is filled with a chain of plasmoids of various sizes, ranging from layer-thickness scales to stellar-radii scales for the largest ones. The layer is also slightly kinked due to current-driven instabilities as also reported in previous studies (e.g., see Cerutti et al. 2015; Belyaev 2015). Large vacuum gaps form outside the layer, a feature also already reported in previous studies (Chen & Beloborodov, 2014; Philippov et al., 2015a). The transition layer in the blending function (Eq. 30) allows to bridge the pure force-free region to these vacuum gaps without noticeable numerical artefacts. We have also verified that this feature is robust against the latitudinal extent of the PIC domain. While plasma depletion is expected in this regime, it is probably exaggerated here because of our pair production scheme, where the pair producing photons have virtually zero mean-free path. However, we believe that this has no dramatic effect on the high-energy output of the magnetosphere described below.
Perhaps the most spectacular feature is the extreme narrowness of the current layer, the latter being barely visible in the global view (see the zoomed-in view in Figure 3). Near the light cylinder, the layer thickness is about which corresponds to the scale separation anticipated in Sect. 2, indicating that the simulation has achieved the desired scale separation. The pair multiplicity in the layer is , while the proton density is at most of the order unity. The pairs being much more abundant than the protons, the layer thickness is controlled by the electronic scales. In contrast, the proton beam is much thicker, with a width near the light-cylinder, and it is not significantly perturbed by the presence of plasmoids. The decoupling between the proton and the pair scales is amplified by radiative cooling, for which the protons are nearly immune from.
Figure 3 also maps out the mean particle energies in the magnetosphere, demonstrating the crucial role of reconnection in particle acceleration. X-points, where the field reconnects in-between two plasmoids, are the main acceleration sites for the pairs. An examination of the bulk velocities reveals that electrons and positrons form two counter-streaming beams at X-points. Positrons (and protons) move radially outwards along with the wind and plasmoids, while electrons tend to precipitate inwards. Counter-streaming is needed to carry the magnetospheric return current (Cerutti et al., 2015; Contopoulos et al., 2020); it is strongest near the light cylinder and tends to diminish with distance. This effect could also be reduced by a larger pair multiplicity in the current sheet. The zoomed-in view on one of this X-point present in the layer (see the inset in Figure 3) shows that the energy of the pairs peaks inside the layer, TeV, and sharply decreases inside the nearest plasmoids on either side, down to GeV. In this extreme synchrotron cooling regime, particles accelerate deep inside the layer where the effective magnetic field strength is low, allowing them to accelerate beyond the radiation-reaction-limited energy, GeV; but as soon as they leave the layer and enter a plasmoid, they feel a sharp increase of the perpendicular magnetic field, they cool off catastrophically and radiate all their energy away in the form of synchrotron photons (Cerutti et al., 2013). Plasmoids effectively act as beam dumps, as clearly visible in Figure 3. Thus, other acceleration processes involving, for instance, plasmoid contraction or particle escape from plasmoids are irrelevant for the pairs in this environment, in contrast to non-radiative reconnection (Petropoulou & Sironi, 2018; Zhang et al., 2021). However, plasmoids, X-points, and radiative cooling are unimportant for ions. From their perspective, the layer is a sharp featureless magnetic null along which they follow relativistic Speiser orbits (Speiser, 1965; Contopoulos, 2007; Cerutti et al., 2013; Contopoulos & Stefanou, 2019; Zhang et al., 2021; Chernoglazov et al., 2023), allowing them to reach extremely-high energies, TeV, close to the full polar-cap potential drop (Philippov & Spitkovsky, 2018; Guépin et al., 2020). This energy scale is consistent with system-size-limited acceleration achievable with reconnection, given by TeV, where -0.2 is the reconnection rate.
Figure 4 shows the energy distribution for all the charged particle species. The pair spectrum ranges from about 1GeV to 1TeV and peaks at about a few tens of GeV, below the radiation-reaction limited energy (GeV). While there is no major difference between electrons and positrons at low energies, positrons significantly dominate over the electrons at the highest energies. Positrons clearly show an extra spectral component peaking 300-400GeV, which coincides with particles energized at X-points whose energy is determined by the magnetization parameter GeV (Eq. 9). This feature was already reported in previous works (Cerutti et al., 2015), and seems to hold strong at more realistic scales. This asymmetry is due to the overall polarization of the magnetosphere where a positive net charge is expected in the equatorial regions, if . Note that this effect is strongest for an aligned rotator and vanishes for an orthogonal rotator (Philippov & Spitkovsky, 2018). As anticipated in Sect. 2, the pair spectrum is relatively narrow in energy band for a weak millisecond gamma-ray pulsar because of the proximity between the relevant energy scales. The proton spectrum is composed of a single component ranging from about 1TeV up to 30TeV. If we venture to fit a power law, the spectral index would be , in qualitative agreement with local reconnection studies (Melzani et al., 2014; Werner et al., 2018). Although the exact spectral shape of all species must be considered with caution, given the oversimplified model used for pair production in this work, we believe that the energy scales, width, and budget are robust overall. Scaling down the magnetic field strength, and keeping all energy scale separations the same with respect to the polar-cap potential drop , does not change the above picture. All magnetospheric features are recovered and the particle maximum energies scale linearly with the field strength (Figure 5). This result suggests that the rescaling performed in previous studies is legitimate as long as the hierarchy of scales of the targeted system is respected.
The narrowness of the pair energy spectrum translates into a narrow-band synchrotron energy distribution. The total synchrotron spectrum reported in Figure 6 ranges from about 1MeV to 1GeV, with a peak in the energy band of the Fermi-LAT above 100MeV. It includes the emission from the equatorial current layer only, and thus synchrotron radiation largely dominates over curvature radiation because the Larmor radius is much smaller than the particle guiding center curvature radius (). Given that the emission is strongly beamed, it corresponds to the spectrum an observer would see if looking in the equatorial plane of the system. The GeV component corresponds to the particles (mostly positrons) accelerated at X-points beyond the synchrotron burnoff limit, given by MeV (Uzdensky et al., 2011), where is the fine structure constant. Integrating over all frequencies, the synchrotron luminosity accounts for about of the pulsar spindown. This sizeable fraction corresponds to the amount of dissipation of the Poynting flux measured in the magnetosphere within a few light-cylinder radii away from the base of the layer that is controlled by the reconnection rate (Cerutti et al., 2020; Hakobyan et al., 2023). It is way larger than the numerical dissipation reported in Sect. 4 (), thus demonstrating that dissipation is by far physical and due to magnetic reconnection. This also means that nearly all of the dissipated power in the box is channelled into non-thermal radiation, which is possible only in the strong cooling regime. Indeed, about of the spindown is channelled into the beam of protons while the pairs leave the box with of the spindown power. The energy fraction carried by the pairs and the protons would increase at larger distances where dissipation proceeds and synchrotron cooling becomes inefficient. Integrated above the energy threshold of the Fermi-LAT, our fiducial pulsar shines about of its total spindown power above MeV range. Compared to the Fermi pulsar population, our solution is compatible with the typical reported gamma-ray efficiencies that range from 1-100% (Smith et al. 2023, see Figure 7).
6 Conclusion
In this work, we have developed a new hybrid scheme designed to blend the PIC and the force-free approaches in the context of relativistic magnetospheres. This effort is motivated by the need to increase the separation between the kinetic scales where particle acceleration operates and global scales, and by this means add more credibility to global PIC models in spite of their unrealistically low scale separations. Using these new capabilities and high-resolution simulations, we focused our numerical resources on the modeling of a weak millisecond gamma-ray pulsar using realistic (i.e., without artificial rescaling) energy scales that are relevant in the equatorial current sheet, while sacrificing the polar cap physics near the star surface. This hybrid approach is about four times faster than a full PIC simulation of the same problem with identical numerical parameters. Perhaps more importantly, it allows to circumvent the severe resolution constraints imposed by the physical conditions at the polar caps where the plasma is not as energetic as in the current layers, which would prevent a full PIC simulation to model such system with the same numerical resolution as the hybrid run.
Our results recover all features reported by previous studies: a prominent reconnecting current layer efficiently accelerating pairs and ions. For the surface magnetic field G simulated in this work, pairs are accelerated up to 1TeV, while protons reach up to of the full polar cap potential drop, that is TeV. With a luminosity of a few percent of the pulsar spindown power, the millisecond pulsar population could thus be a sizeable contributor to the Galactic cosmic-ray population (Philippov & Spitkovsky, 2018; Guépin et al., 2020). The extreme cooling rate of pairs leads to a nearly complete conversion of the pair energy into synchrotron radiation, including above the burnoff limit in the GeV range. The bolometric radiative efficiency reaches about of the spindown power, including about above 100MeV gamma rays, a number compatible with the Fermi-LAT pulsar population (Smith et al., 2023). The fact that pairs can be accelerated at least up to TeV energies, even for a weak pulsar, is encouraging for explaining the recent detections of pulsed TeV emission in the Crab and the Vela pulsars (Ansoldi et al., 2016; H. E. S. S. Collaboration et al., 2023).
We have also shown that the rescaling procedure operated by previous PIC studies, in particular in the surface magnetic field strength, is valid as long as a realistic hierarchy of scales is respected with regard to the targeted system. Thus, the results from this simulation may not be extrapolated to larger field strength, as for instance the young pulsar population, because the scale separation would be different for a larger surface field strength. This work is meant to provide a proof of concept that magnetic reconnection within the equatorial current sheet is a viable scenario to explain the gamma-ray emission, without scaling down the physical parameters, at least in the low-luminosity range of the gamma-ray pulsar population. In this respect, our fiducial simulation has successfully achieved this objective, although much remains to be done in the future. In particular, including a more realistic model for pair production could significantly improve the predictive power of the simulations, in particular regarding the pair multiplicities in the magnetosphere, the particle/photon spectral shapes and cutoff energies, as well as the radiative efficiencies. These last two quantities reported in this work being low compared to the bulk of the Fermi-LAT millisecond pulsar population. Longer integration times to establish a more robust steady state, in particular inside the light cylinder (separatrices and Y-point), would be desirable. The hybrid method presented in this paper provides interesting new avenues to model relativistic magnetospheres, including those forming around black holes. More work is needed to apply this method to a more diverse array of astrophysical problems, and to generalize it into a full 3D setup.
Acknowledgements.
We thank Kyle Parfrey and Jens Mahlmann for their valuable comments and encouragements during the implementation of this new numerical scheme. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 863412). Computing resources were provided by TGCC under the allocation A0130407669 made by GENCI. We thank the International Space Science Institute (ISSI) for providing financial support for the organization of the meeting of ISSI Team No 459 led by I. Contopoulos and D. Kazanas where this project has been discussed. Software: Matplotlib (Hunter 2007), Numpy (Harris et al. 2020).References
- Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460
- Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17
- Ansoldi et al. (2016) Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133
- Arka & Dubus (2013) Arka, I. & Dubus, G. 2013, A&A, 550, A101
- Arons (2003) Arons, J. 2003, ApJ, 589, 871
- Bai & Spitkovsky (2010) Bai, X.-N. & Spitkovsky, A. 2010, ApJ, 715, 1282
- Belyaev (2015) Belyaev, M. A. 2015, MNRAS, 449, 2759
- Birdsall & Langdon (1991) Birdsall, C. K. & Langdon, A. B. 1991, Plasma Physics via Computer Simulation
- Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467
- Blandford (2002) Blandford, R. D. 2002, in Lighthouses of the Universe: The Most Luminous Celestial Objects and Their Use for Cosmology, ed. M. Gilfanov, R. Sunyeav, & E. Churazov, 381
- Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
- Blasi et al. (2000) Blasi, P., Epstein, R. I., & Olinto, A. V. 2000, ApJ, 533, L123
- Blumenthal & Gould (1970) Blumenthal, G. R. & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237
- Bransgrove et al. (2023) Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9
- Cao & Yang (2019) Cao, G. & Yang, X. 2019, ApJ, 874, 166
- Cao et al. (2016) Cao, G., Zhang, L., & Sun, S. 2016, MNRAS, 455, 4267
- Cerutti & Beloborodov (2017) Cerutti, B. & Beloborodov, A. M. 2017, Space Sci. Rev., 207, 111
- Cerutti et al. (2015) Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606
- Cerutti & Philippov (2017) Cerutti, B. & Philippov, A. A. 2017, A&A, 607, A134
- Cerutti et al. (2020) Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204
- Cerutti et al. (2016) Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401
- Cerutti et al. (2013) Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147
- Chen & Beloborodov (2014) Chen, A. Y. & Beloborodov, A. M. 2014, ApJ, 795, L22
- Chen et al. (2020) Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, ApJ, 889, 69
- Chernoglazov et al. (2023) Chernoglazov, A., Hakobyan, H., & Philippov, A. 2023, ApJ, 959, 122
- Contopoulos (2007) Contopoulos, I. 2007, A&A, 472, 219
- Contopoulos et al. (1999) Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351
- Contopoulos et al. (2024) Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. N. 2024, MNRAS, 527, L127
- Contopoulos et al. (2020) Contopoulos, I., Pétri, J., & Stefanou, P. 2020, MNRAS, 491, 5579
- Contopoulos & Stefanou (2019) Contopoulos, I. & Stefanou, P. 2019, MNRAS, 487, 952
- Coroniti (1990) Coroniti, F. V. 1990, ApJ, 349, 538
- Crinquand et al. (2019) Crinquand, B., Cerutti, B., & Dubus, G. 2019, A&A, 622, A161
- Cruz et al. (2023) Cruz, F., Grismayer, T., Chen, A. Y., et al. 2023, arXiv e-prints, arXiv:2309.04834
- Cruz et al. (2021) Cruz, F., Grismayer, T., Chen, A. Y., Spitkovsky, A., & Silva, L. O. 2021, ApJ, 919, L4
- Erber (1966) Erber, T. 1966, Reviews of Modern Physics, 38, 626
- Fang et al. (2012) Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118
- Goldreich & Julian (1969) Goldreich, P. & Julian, W. H. 1969, ApJ, 157, 869
- Gruzinov (1999) Gruzinov, A. 1999, arXiv e-prints, astro
- Gruzinov (2008) Gruzinov, A. 2008, arXiv e-prints, arXiv:0802.1716
- Guépin et al. (2020) Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138
- Gunn & Ostriker (1969) Gunn, J. E. & Ostriker, J. P. 1969, Phys. Rev. Lett., 22, 728
- H. E. S. S. Collaboration et al. (2023) H. E. S. S. Collaboration, Aharonian, F., Ait Benkhali, F., et al. 2023, Nature Astronomy, 7, 1341
- Hakobyan et al. (2023) Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, 943, 105
- Harding & Lai (2006) Harding, A. K. & Lai, D. 2006, Reports on Progress in Physics, 69, 2631
- Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
- Hones & Bergeson (1965) Hones, Edward W., J. & Bergeson, J. E. 1965, J. Geophys. Res., 70, 4951
- Hu & Beloborodov (2022) Hu, R. & Beloborodov, A. M. 2022, ApJ, 939, 42
- Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
- Kalapotharakos et al. (2018) Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44
- Kalapotharakos et al. (2014) Kalapotharakos, C., Harding, A. K., & Kazanas, D. 2014, ApJ, 793, 97
- Kaspi & Beloborodov (2017) Kaspi, V. M. & Beloborodov, A. M. 2017, ARA&A, 55, 261
- Komissarov (2002) Komissarov, S. S. 2002, MNRAS, 336, 759
- Komissarov (2006) Komissarov, S. S. 2006, MNRAS, 367, 19
- Landau & Lifshitz (1971) Landau, L. D. & Lifshitz, E. M. 1971, The classical theory of fields, ed. Landau, L. D. & Lifshitz, E. M.
- Li et al. (2012) Li, J., Spitkovsky, A., & Tchekhovskoy, A. 2012, ApJ, 746, 60
- Lyubarsky (2020) Lyubarsky, Y. 2020, ApJ, 897, 1
- Lyubarsky & Kirk (2001) Lyubarsky, Y. & Kirk, J. G. 2001, ApJ, 547, 437
- McKinney (2006) McKinney, J. C. 2006, MNRAS, 368, L30
- Melzani et al. (2014) Melzani, M., Walder, R., Folini, D., Winisdoerffer, C., & Favre, J. M. 2014, A&A, 570, A112
- Michel (1973) Michel, F. C. 1973, ApJ, 180, L133
- Michel & Tucker (1969) Michel, F. C. & Tucker, W. H. 1969, Nature, 223, 277
- Most & Philippov (2020) Most, E. R. & Philippov, A. A. 2020, ApJ, 893, L6
- Parfrey et al. (2012) Parfrey, K., Beloborodov, A. M., & Hui, L. 2012, MNRAS, 423, 1416
- Petropoulou & Sironi (2018) Petropoulou, M. & Sironi, L. 2018, MNRAS, 481, 5687
- Philippov & Kramer (2022) Philippov, A. & Kramer, M. 2022, ARA&A, 60, 495
- Philippov et al. (2020) Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, Phys. Rev. Lett., 124, 245101
- Philippov et al. (2015a) Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJ, 815, L19
- Philippov & Spitkovsky (2014) Philippov, A. A. & Spitkovsky, A. 2014, ApJ, 785, L33
- Philippov & Spitkovsky (2018) Philippov, A. A. & Spitkovsky, A. 2018, ApJ, 855, 94
- Philippov et al. (2015b) Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015b, ApJ, 801, L19
- Popov & Postnov (2010) Popov, S. B. & Postnov, K. A. 2010, in Evolution of Cosmic Objects through their Physical Activity, ed. H. A. Harutyunian, A. M. Mickaelian, & Y. Terzian, 129–132
- Ruderman & Sutherland (1975) Ruderman, M. A. & Sutherland, P. G. 1975, ApJ, 196, 51
- Scharlemann & Wagoner (1973) Scharlemann, E. T. & Wagoner, R. V. 1973, ApJ, 182, 951
- Smith et al. (2023) Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191
- Speiser (1965) Speiser, T. W. 1965, J. Geophys. Res., 70, 4219
- Spitkovsky (2006) Spitkovsky, A. 2006, ApJ, 648, L51
- Sturrock (1971) Sturrock, P. A. 1971, ApJ, 164, 529
- Thompson & Duncan (1995) Thompson, C. & Duncan, R. C. 1995, MNRAS, 275, 255
- Timokhin (2006) Timokhin, A. N. 2006, MNRAS, 368, 1055
- Timokhin & Arons (2013) Timokhin, A. N. & Arons, J. 2013, MNRAS, 429, 20
- Torres et al. (2024) Torres, R., Grismayer, T., Cruz, F., Fonseca, R. A., & Silva, L. O. 2024, arXiv e-prints, arXiv:2401.02908
- Uzdensky (2003) Uzdensky, D. A. 2003, ApJ, 598, 446
- Uzdensky et al. (2011) Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJ, 737, L40
- Uzdensky & Spitkovsky (2014) Uzdensky, D. A. & Spitkovsky, A. 2014, ApJ, 780, 3
- Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840
- Yee (1966) Yee, K. 1966, IEEE Transactions on Antennas and Propagation, 14, 302
- Zhang et al. (2021) Zhang, H., Sironi, L., & Giannios, D. 2021, ApJ, 922, 261