Issue 
A&A
Volume 690, October 2024



Article Number  A170  
Number of page(s)  12  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202450238  
Published online  08 October 2024 
Scaling up global kinetic models of pulsar magnetospheres using a hybrid forcefreePIC numerical approach
^{1}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
^{2}
Research Center for Astronomy and Applied Mathematics, Academy of Athens, GR 11527 Athens, Greece
Received:
4
April
2024
Accepted:
12
June
2024
Context. The particleincell approach has proven effective in modeling neutronstar and blackhole magnetospheres from first principles, but global simulations are plagued with an unrealistically small separation between the scales where microphysics operates and the systemsize scales due to limited numerical resources. A legitimate concern is whether the scale separation achieved to date is large enough for results to be safely extrapolated to realistic scales.
Aims. In this work, our aim is to explore the effect of scaling up physical parameters 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 highenergy radiation emitted beyond the light cylinder.
Methods. To reach this objective, we developed a new hybrid numerical scheme coupling the ideal forcefree and the particleincell methods to optimize the numerical cost of global models. We propose a domain decomposition of the magnetosphere based on the magneticfield topology using the flux function. The forcefree model is enforced along open field lines while the particleincell 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 highresolution axisymmetric simulations. Magnetospheric features reported by previous kinetic models are recovered, and strong synchrotron radiation above 100MeV consistent with the FermiLAT gammaray pulsar population is successfully reproduced.
Conclusions.This work further consolidates the shiningreconnecting current sheet scenario as the origin of the gammaray emission in pulsars, as well as firmly establishing pulsar magnetospheres as at least teraelectronvolt particle accelerators.
Key words: acceleration of particles / magnetic reconnection / radiation mechanisms: nonthermal / methods: numerical / pulsars: general / stars: winds / outflows
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Magnetospheres forming around neutron stars and black holes, simply referred to as relativistic magnetospheres in the following, are involved in some of the most extreme highenergy astrophysical phenomena such as pulsars (Cerutti & Beloborodov 2017; Philippov & Kramer 2022), magnetar flares and outbursts (Thompson & Duncan 1995; Kaspi & Beloborodov 2017), relativistic jets from supermassive black holes (Blandford & Znajek 1977; Blandford et al. 2019), 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 that the interplay between plasma physics, general relativity, and relativistic (quantum) electrodynamics dramatically increases the degree of complexity of the problem, which 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 forcefree 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 and 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 GradShafranov 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 forcefree 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 lightcylinder radius, where corotation with the star is superluminal (Goldreich & Julian 1969; Michel & Tucker 1969). Thereby, the forcefree framework provides a realistic description of the morphology of relativistic magnetospheres in the limit where plasma is abundant and ultramagnetized everywhere. It also captures how the star spins down under the effect of magnetic torques well. 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 forcefree formulation), and more importantly nonthermal 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 rotationpowered pulsars in the gammaray 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 20 TeV (Ansoldi et al. 2016; H. E. S. S. Collaboration 2023), which represents a sizeable fraction of the full polarcap 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 radiationreaction force, and general relativistic effects. In this respect, the particleincell (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. 2024; 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 nonthermal particle acceleration and energetic radiation, possibly explaining the observed highenergy radiation (Cerutti et al. 2016; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018). At the base of the open field lines near the star, the sparkgap dynamics and the putative radio emission are also reproduced by local PIC models (Timokhin & Arons 2013; Philippov et al. 2020; Cruz et al. 2021) well.
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., lightcylinder) 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 polarcap physics and the lightcylinder reconnection physics with full scales in the same simulation seems impossible, even in the distant future. In this work, we made the choice to sacrifice the polarcap microphysics to focus our numerical resources on the lightcylinder 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 up physical parameters 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 forcefree and PIC approaches. We propose a domain decomposition based on the magneticfield topology of the system, where open field lines are described by an ideal forcefree 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 of reaching realistic scales of a weak millisecond pulsar, just above the limit of the observed FermiLAT 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 forcefreePIC numerical scheme. This approach is first tested against the monopole solution (Sect. 4). It is then applied to the aligned dipole with highresolution 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 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 consider, as our baseline, a millisecond pulsar with a rotation period of P = 1 ms and a surface dipolar magnetic field of B_{⋆} = 10^{7} G, corresponding to the weakest observable gammaray pulsars of spindown power L ≈ 4.8 × 10^{33}erg/s as reported by the FermiLAT (Abdo et al. 2010, 2013; Smith et al. 2023). It also represents a realistic gammaray pulsar with the smallest scale separation that we strive to capture in this work as a proof of principles (see Sect. 5).
The systemsize macroscopic scales are set by the neutron star radius, r_{⋆} ∼ 10^{6} cm, and the lightcylinder radius,
$$\begin{array}{c}\hfill {R}_{\mathrm{LC}}=\frac{\mathit{cP}}{2\pi}\approx 5{r}_{\star}\left(\frac{P}{1\phantom{\rule{0.166667em}{0ex}}\mathrm{ms}}\right).\end{array}$$(1)
The lightcylinder 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 “Ypoint” (e.g., Uzdensky 2003; Contopoulos et al. 2024), and it proceeds at all radii beyond this point. Energetic synchrotron emission is produced within a few lightcylinder radii because of the steep decay of the magneticfield strength with radius (Cerutti et al. 2016; Cerutti et al. in prep.), and thus it sets the largest relevant scale of the problem in this work.
While the macroscopic scales are well constrained, defining the microscopic scales is more complex because plasma physics and quantum electrodynamics (QED) effects must be taken into consideration. The electron plasma skin depth at the star surface sets the minimal relevant scale of the problem. Assuming that the plasma density is a multiple of the critical corotation density, ${n}_{\star}=\kappa {n}_{\mathrm{GJ}}^{\star}=\kappa {B}_{\star}/2\pi e{R}_{\mathrm{LC}}$ (Hones & Bergeson 1965; Goldreich & Julian 1969), where κ is the plasma multiplicity and e is the electron charge, yields
$$\begin{array}{cc}\hfill {d}_{\mathrm{e}}^{\star}& =\sqrt{\frac{\gamma {m}_{\mathrm{e}}{c}^{2}}{4\pi {n}_{\star}{e}^{2}}}\hfill \\ \hfill & \approx 2{\left(\frac{\gamma}{{10}^{0}}\right)}^{1/2}{\left(\frac{\kappa}{{10}^{2}}\right)}^{1/2}{\left(\frac{P}{1\mathrm{ms}}\right)}^{1/2}{\left(\frac{{B}_{\star}}{{10}^{7}\mathrm{G}}\right)}^{1/2}\mathrm{cm},\hfill \end{array}$$(2)
where γ is the particle Lorentz factor and m_{e} is the electron mass; so, ${d}_{\mathrm{e}}^{\star}/{r}_{\star}=2\times {10}^{6}$ 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 ultrarelativistic speeds (γ ≫ 1), thereby 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 $sin{\theta}_{\mathrm{pc}}=\sqrt{{r}_{\star}/{R}_{\mathrm{LC}}}$. Using the corotation surface electric field, E = −(Ω × r)×B/c, we can derive the potential drop across the polar cap as
$$\begin{array}{c}\hfill {\mathrm{\Phi}}_{\mathrm{pc}}=\frac{\mu {\mathrm{\Omega}}^{2}}{{c}^{2}},\end{array}$$(3)
where Ω = 2π/P is the angular velocity of the star, and μ = B_{⋆}r_{⋆}^{3} is its magnetic moment. An electron experiencing the full potential drop will acquire a Lorentz factor given by
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{pc}}=\frac{e{\mathrm{\Phi}}_{\mathrm{pc}}}{{m}_{\mathrm{e}}{c}^{2}}\approx 2.6\times {10}^{8}\left(\frac{{B}_{\star}}{{10}^{7}\mathrm{G}}\right){\left(\frac{P}{1\mathrm{ms}}\right)}^{2}.\end{array}$$(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 large 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 magnetic conversion, meaning that a gammaray 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)
$$\begin{array}{c}\hfill \chi \equiv \u03f5b\gtrsim 0.1,\end{array}$$(5)
where ϵ = ℏν/m_{e}c^{2} is the dimensionless photon energy, and $b={\stackrel{\sim}{B}}_{\perp}/{B}_{\mathrm{QED}}$ is the effective perpendicular magnetic field (see Eq. 17 below) normalized to the critical quantum field B_{QED} = m_{e}^{2}c^{3}/ℏe ≈ 4.4 × 10^{13} G. If electrons radiate via synchrotroncurvature radiation as usually assumed, the critical photon energy is given by (Blumenthal & Gould 1970)
$$\begin{array}{c}\hfill {\u03f5}_{\mathrm{c}}=\frac{3}{2}b{\gamma}^{2}.\end{array}$$(6)
Combining Eq. (5) with Eq. (6) provides an estimate for the threshold electron Lorentz factor capable of producing new pairs,
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{th}}=\sqrt{\frac{1}{15{b}^{2}}}\approx {10}^{6}{\left(\frac{{B}_{\star}}{{10}^{7}\mathrm{G}}\right)}^{1},\end{array}$$(7)
while the created pair shares the absorbed photon momentum, such that the Lorentz factor of the secondary pairs may be estimated as
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{s}}=\frac{\u03f5}{2}=\frac{1}{20b}\approx 2.2\times {10}^{5}{\left(\frac{{B}_{\star}}{{10}^{7}\mathrm{G}}\right)}^{1}.\end{array}$$(8)
Assuming that pair creation is effective at producing a highmultiplicity plasma in the magnetosphere (κ ≫ 1), it sets the relevant scale of the particle Lorentz factor to be considered in estimating the plasma skin depth scale in Eq. (2), yielding ${d}_{\mathrm{e}}^{\mathrm{s}}={d}_{\mathrm{e}}^{\star}({\gamma}_{\mathrm{s}})\approx {10}^{3}\mathrm{cm}={10}^{3}{r}_{\star}$ for our reference pulsar. One can realize from these orderofmagnitude calculations that a weaker magnetic field (B_{⋆} ≲ 10^{7} 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 skin depth scale, meaning that $\delta \sim {d}_{\mathrm{e}}^{\mathrm{LC}}$. 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,
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{LC}}\sim {\sigma}_{\mathrm{LC}}=\frac{{B}_{\mathrm{LC}}^{2}}{4\pi {\mathrm{\Gamma}}_{\mathrm{LC}}\kappa {n}_{\mathrm{GJ}}^{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}}=\frac{{\gamma}_{\mathrm{pc}}}{2{\mathrm{\Gamma}}_{\mathrm{LC}}\kappa},\end{array}$$(9)
where B_{LC} ∼ B_{⋆}(r_{⋆}/R_{LC})^{3} = μΩ^{3}/c^{3}, ${n}_{\mathrm{GJ}}^{\mathrm{LC}}=\mathrm{\Omega}{B}_{\mathrm{LC}}/2\pi ec$ is the GoldreichJulian plasma density at the light cylinder and Γ_{LC} is the pulsar wind bulk Lorentz factor at its base. Assuming that these energetic particles set the local skin depth scale, this yields
$$\begin{array}{c}\hfill \delta \sim {d}_{\mathrm{e}}^{\mathrm{LC}}({\gamma}_{\mathrm{LC}})=\sqrt{\frac{{\gamma}_{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}}{4\pi \kappa {\mathrm{\Gamma}}_{\mathrm{LC}}{n}_{\mathrm{GJ}}^{\mathrm{LC}}{e}^{2}}}=\frac{{R}_{\mathrm{LC}}}{2{\mathrm{\Gamma}}_{\mathrm{LC}}\kappa},\end{array}$$(10)
$$\begin{array}{c}\hfill \delta \sim 2.4\times {10}^{4}{\left(\frac{{\mathrm{\Gamma}}_{\mathrm{LC}}}{{10}^{0}}\right)}^{1}{\left(\frac{\kappa}{{10}^{2}}\right)}^{1}\left(\frac{P}{1\mathrm{ms}}\right)\mathrm{cm},\end{array}$$(11)
thus, δ/r_{⋆} ∼ 2 × 10^{−2} or δ/R_{LC} ∼ 5 × 10^{−3} for these fiducial parameters. A similar estimate can be obtained using 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 radiationreaction force, eE_{∥} ∼ f_{rad}. The latter effectively acts as a continuous drag force (in the classical regime, i.e., γb ≪ 1), opposite to the particle’s direction of motion, whose magnitude in the ultrarelativistic regime (γ ≫ 1) can be estimated as ${f}_{\mathrm{rad}}\approx (2/3){r}_{\mathrm{e}}^{2}{\gamma}^{2}{\stackrel{\sim}{B}}_{\perp}^{2}$, where r_{e} is the classical radius of the electron. The fiducial radiationreactionlimited electron Lorentz factor is then given by
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{rad}}=\sqrt{\frac{3e{E}_{\Vert}}{2{r}_{\mathrm{e}}^{2}{\stackrel{\sim}{B}}_{\perp}^{2}}}\approx 3\times {10}^{4}{\left(\frac{{E}_{\Vert}}{{\stackrel{\sim}{B}}_{\perp}}\right)}^{1/2}{\left(\frac{{\stackrel{\sim}{B}}_{\perp}}{{10}^{7}\mathrm{G}}\right)}^{1/2},\end{array}$$(12)
or ${\gamma}_{\mathrm{rad}}^{\mathrm{LC}}\approx 3.3\times {10}^{5}$ at the light cylinder for the fiducial pulsar parameters, assuming ${E}_{\Vert}={\stackrel{\sim}{B}}_{\perp}$ and ${\stackrel{\sim}{B}}_{\perp}={B}_{\mathrm{LC}}$.
In summary, modeling the weakest gammaray pulsar must satisfy at the very least the following scale separation, in terms of spatial quantities normalized to r_{⋆}:
$$\begin{array}{c}\hfill \frac{{d}_{\mathrm{e}}^{\mathrm{s}}}{{r}_{\star}}\sim {10}^{3}\ll \frac{\delta}{{r}_{\star}}\sim 2.5\times {10}^{2}\ll 1<\frac{{R}_{\mathrm{LC}}}{{r}_{\star}}=5,\end{array}$$(13)
and in terms of energy scales normalized to γ_{pc},
$$\begin{array}{c}\hfill \frac{{\gamma}_{\mathrm{s}}}{{\gamma}_{\mathrm{pc}}}\sim 8.5\times {10}^{4}<\frac{{\gamma}_{\mathrm{rad}}^{\mathrm{LC}}}{{\gamma}_{\mathrm{pc}}}\sim {10}^{3}<\frac{{\gamma}_{\mathrm{th}}}{{\gamma}_{\mathrm{pc}}}\sim 4\times {10}^{3}\ll 1.\end{array}$$(14)
This exercise reveals that in such systems γ_{s}, γ_{rad}, and γ_{th} 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 FermiLAT pulsars, which shows that the observed spectral energy distributions are narrower with decreasing spindown, nearly reaching a spectrum consistent with a monoenergetic particle population for the weakest pulsars with a spindown of L ∼ 10^{33} erg/s (Smith et al. 2023). Thus, even though the estimates presented in this section by no means represent a rigorous calculation, it provides an approximate assessment of the minimum computational needs to approach realistic systems using a kinetic approach. Our conclusion is that such a weak pulsar could in principle be captured with real scales using current computational facilities, provided that extreme resolutions are used in a 2D axisymmetric system, as, for example, in the recent works by Hu & Beloborodov (2022) and Bransgrove et al. (2023).
3. Hybrid forcefreePIC approach
In order to fulfill the minimum scale separation requirements highlighted in the previous section, a large number of grid cells must be used, and, hence, high computational power is needed. To alleviate some of this numerical cost, we developed a new hybrid approach specially designed for aligned magnetospheres, coupling timedependent forcefree electrodynamic and PIC approaches in the same simulation. This idea relies on the fact that gammaray pulsar magnetospheres are mostly very close to a dissipationless forcefree state because of efficient pair creation (i.e., κ ≫ 1). Using the expensive full PIC technique with a large number of particles per cell everywhere in the magnetosphere may seem excessive, as a numerically cheaper forcefree description would be sufficient. In contrast, dissipative regions (e.g., electrostatic gaps, current sheets) should remain on microscopic scales (again, if κ ≫ 1), and so these are the scales for 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.
In this work, we developed a new forcefree 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 pointlike superparticles – simply referred to as particles in the following– representing a very large number of physical particles with identical masstocharge ratios following the same path in phase space (Birdsall & Langdon 1991). Their evolution is governed by the Lorentz force, and the radiationreaction force in this pulsarspecific context, as
$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\frac{q}{\mathit{mc}}(\mathbf{E}+\frac{\mathbf{u}\times \mathbf{B}}{\gamma})+\mathbf{g},\end{array}$$(15)
where u = γv/c is the dimensionless particle fourmomentum vector, and q is the particle electric charge, and where
$$\begin{array}{c}\hfill \mathbf{g}\approx \frac{2}{3}{r}_{\mathrm{e}}^{2}[(\mathbf{E}+\mathit{\beta}\times \mathbf{B})\times \mathbf{B}+(\mathit{\beta}\xb7\mathbf{E})\mathbf{E}]\frac{2}{3}{r}_{\mathrm{e}}^{2}{\gamma}^{2}{\stackrel{\sim}{B}}_{\perp}^{2}\mathit{\beta}\end{array}$$(16)
is the radiation reaction force following the Landau & Lifshitz (1971) formulation, and β = v/c, and where
$$\begin{array}{c}\hfill {\stackrel{\sim}{B}}_{\perp}=\sqrt{{(\mathbf{E}+\mathit{\beta}\times \mathbf{B})}^{2}{(\mathit{\beta}\xb7\mathbf{E})}^{2}}\end{array}$$(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 J reconstructed on the grid is then used to evolve the timedependent Maxwell equations:
$$\begin{array}{c}\hfill \frac{\partial \mathbf{B}}{\partial t}=c\phantom{\rule{0.166667em}{0ex}}\mathbf{\nabla}\times \mathbf{E}\end{array}$$(18)
$$\begin{array}{c}\hfill \frac{\partial \mathbf{E}}{\partial t}=c\phantom{\rule{0.166667em}{0ex}}\mathbf{\nabla}\times \mathbf{B}4\pi \mathbf{J},\end{array}$$(19)
allowing us to close the PIC loop performed at each time step (Fig. 1). The field integration step follows a standard finitedifference timedomain method that involves a staggered mesh (Yee 1966).
Fig. 1. New hybrid numerical scheme proposed in this work during a timeintegration cycle. This method is intended to blend the forcefree electrodynamic (FFE) and the full kinetic (PIC) approaches within the same numerical framework using a domain decomposition. The coupling is achieved via the electric current densities, J, using a blending function, f, which solely depends on the magnetic flux function, Ψ. 
In general, each simulation cell is filled with a large number of particles, 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 forcefree electrodynamics come into play. The forcefree condition is given by
$$\begin{array}{c}\hfill {F}_{\mu \nu}{I}^{\nu}=0,\end{array}$$(20)
where F_{μν} is the electromagnetic tensor and I^{ν} is the fourcurrent, which translates into the following conditions:
$$\begin{array}{c}\hfill \rho \mathbf{E}+\frac{\mathbf{J}\times \mathbf{B}}{c}=\mathbf{0},\end{array}$$(21)
using the spatial components, where ρ = ∇ ⋅ E/4π is the charge density; and
$$\begin{array}{c}\hfill \mathbf{E}\xb7\mathbf{J}=0,\end{array}$$(22)
using the temporal component. Eq. (21) also leads to the condition
$$\begin{array}{c}\hfill \mathbf{E}\xb7\mathbf{B}=0.\end{array}$$(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)
$$\begin{array}{c}\hfill \mathbf{J}=\frac{c}{4\pi}\mathbf{\nabla}\xb7\mathbf{E}\left(\frac{\mathbf{E}\times \mathbf{B}}{{B}^{2}}\right)+\frac{c}{4\pi}(\mathbf{B}\xb7(\mathbf{\nabla}\times \mathbf{B})\mathbf{E}\xb7(\mathbf{\nabla}\times \mathbf{E}))\frac{\mathbf{B}}{{B}^{2}},\end{array}$$(24)
which can be understood as the sum of a current flowing across field lines, J_{⊥}, and a component flowing along field lines, J_{∥}. 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 forcefree module, we followed a similar numerical procedure as in Spitkovsky (2006), which consists in solely computing the perpendicular forcefree current, J_{⊥}. The parallel component, J_{∥}, prevents the formation of an electricfield component along the magnetic field at all times, E ⋅ B = 0, but computing J_{∥} 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 forcefree condition
$$\begin{array}{c}\hfill \mathbf{E}\leftarrow \mathbf{E}(\mathbf{E}\xb7\mathbf{B})\frac{\mathbf{B}}{{B}^{2}}.\end{array}$$(25)
The code also checks, at every time step, that the condition E^{2} < B^{2} is always satisfied; so, if it is not the case, the electric field is once again renormalized as
$$\begin{array}{c}\hfill \mathbf{E}\leftarrow \sqrt{\frac{{B}^{2}}{{E}^{2}}}\mathbf{E}.\end{array}$$(26)
Computing the forcefree 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 secondorder 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 forcefree 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 that is also numerically convenient for tagging individual field lines. In ordinary spherical coordinates (r, θ), the magnetic flux function is given by
$$\begin{array}{c}\hfill \mathrm{\Psi}=\frac{1}{2\pi}\int \int \mathbf{B}\xb7d\mathbf{S}={\int}_{0}^{\theta}{B}_{\mathbf{r}}{r}^{2}sin{\theta}^{\prime}d{\theta}^{\prime}.\end{array}$$(27)
Once the boundary is set, a thin transition layer is implemented between neighboring field lines where the currents are blended to ensure a better continuity between the PIC and the forcefree currents and to minimize numerical artifacts, such as spurious reflections or accumulation of charges. A simple linear interpolation scheme is sufficient for this purpose, so the full hybrid current is expressed as
$$\begin{array}{c}\hfill \mathbf{J}=(1f\left(\mathrm{\Psi}\right))\phantom{\rule{0.166667em}{0ex}}{\mathbf{J}}_{\mathrm{PIC}}+f\left(\mathrm{\Psi}\right)\phantom{\rule{0.166667em}{0ex}}{\mathbf{J}}_{\mathrm{FFE}},\end{array}$$(28)
where J_{PIC} is the current density reconstructed from the particles, J_{FFE} is the current from the forcefree condition (Eq. 24), and f is the blending function. This hybrid current is then injected into the MaxwellAmpè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 forcefree domain). Considering a single transition layer in the domain, two distinct fluxes must be chosen, Ψ_{0} < Ψ_{1}, so that the blending function is as follows:
$$\begin{array}{cc}\hfill f\left(\mathrm{\Psi}\right)& =1,\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{0.333333em}{0ex}}\text{if}\phantom{\rule{3.33333pt}{0ex}}\mathrm{\Psi}<{\mathrm{\Psi}}_{0},\to \phantom{\rule{0.333333em}{0ex}}\text{pure forcefree, no particles},\hfill \end{array}$$(29)
$$\begin{array}{cc}\hfill f\left(\mathrm{\Psi}\right)& =\frac{{\mathrm{\Psi}}_{1}\mathrm{\Psi}}{{\mathrm{\Psi}}_{1}{\mathrm{\Psi}}_{0}},\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{0.333333em}{0ex}}\text{if}\phantom{\rule{3.33333pt}{0ex}}{\mathrm{\Psi}}_{0}<\mathrm{\Psi}<{\mathrm{\Psi}}_{1},\to \phantom{\rule{0.333333em}{0ex}}\text{transition layer},\hfill \end{array}$$(30)
$$\begin{array}{cc}\hfill f\left(\mathrm{\Psi}\right)& =0,\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{0.333333em}{0ex}}\text{if}\phantom{\rule{3.33333pt}{0ex}}\mathrm{\Psi}>{\mathrm{\Psi}}_{1},\to \phantom{\rule{0.333333em}{0ex}}\text{pure PIC}.\hfill \end{array}$$(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 magneticfield topology in time in the rest of the domain. In this sense, the domain decomposition proposed in this work is not static, but it 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 through modeling a steadystate forcefree monopolar magnetosphere (i.e., a single magnetic monopole). Although unphysical, studying this configuration is of great interest: it admits an exact analytical solution and is free of dissipation – there is no gap or current sheet (as opposed to a splitmonopole) – so spurious numerical effects in the scheme can be quantified. Moreover, 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
$$\begin{array}{c}\hfill {\mathrm{\Psi}}_{\mathrm{M}}={\mathrm{\Psi}}_{\star}(1cos\theta ),\end{array}$$(32)
where Ψ_{⋆} = B_{⋆}r_{⋆}^{2}. An exact solution to the GradShafranov equation yields the following toroidal magnetic field component (Michel 1973):
$$\begin{array}{c}\hfill {B}_{\varphi}=\frac{{\mathrm{\Psi}}_{\star}}{{R}_{\mathrm{LC}}}\frac{sin\theta}{r},\end{array}$$(33)
if Ω ⋅ B > 0 in the northern hemisphere (and Ω ⋅ B < 0 in the southern hemisphere). In this solution, the electric field is purely along θ and matches the toroidal magnetic field, E_{θ} = B_{ϕ}. The electric current is purely radial and is given by
$$\begin{array}{c}\hfill {\mathbf{J}}_{\mathrm{M}}=\frac{c}{4\pi}\mathbf{\nabla}\times \mathbf{B}={J}_{\mathrm{GJ}}^{\star}{\left(\frac{{r}_{\star}}{r}\right)}^{2}cos\theta \phantom{\rule{3.33333pt}{0ex}}{\mathbf{e}}_{\mathrm{r}},\end{array}$$(34)
where ${J}_{\mathrm{GJ}}^{\star}=\mathrm{\Omega}{B}_{\star}/2\pi $ corresponds to the fiducial GoldreichJulian current, and the spindown power associated with the electromagnetic torque is
$$\begin{array}{c}\hfill {L}_{\mathrm{M}}=\frac{c}{4\pi}\int \int (\mathbf{E}\times \mathbf{B})\xb7d\mathbf{S}=\frac{c}{2}{\int}_{1}^{1}{r}^{2}{B}_{\varphi}^{2}dcos\theta =\frac{2{\mathrm{\Psi}}_{\star}^{2}{\mathrm{\Omega}}^{2}}{3c}.\end{array}$$(35)
The goal of this first numerical experiment is to recover these properties using our hybrid scheme. To this end, we performed a set of hybrid simulations in 2D axisymmetric spherical coordinates and explored the effect of resolution on the results. The domain size explored here is composed of 1024^{2}, 2048^{2}, and 4096^{2} cells. The grid spacing is logarithmic along the radial direction and constant along the θ direction. The domain size ranges from the star surface, r_{min} = r_{⋆}, up to r_{max} = (10/3)R_{LC}, where R_{LC} = 5r_{⋆}, and covers the full range in θ ∈ [0, π]. The magnetosphere is initialized in vacuum with a pure radial magnetic field, B = B_{⋆}(r_{⋆}/r)^{2}. The surface magnetic field strength is fixed at B_{⋆} = 5 × 10^{5} G, which is sufficiently intense in this experiment to achieve a quasiforcefree 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, E_{⋆} = −(Ω × r_{⋆})×B_{⋆}/c. This procedure sets the magnetosphere into rotation. An absorbing layer starting at r_{abs} = 3R_{LC} progressively damps all outgoing electromagnetic waves and removes the particles, if there are any (Cerutti et al. 2015). In this experiment, the northern hemisphere is modeled in full PIC, while the southern hemisphere is captured with the forcefree model. This choice is well suited for a comparison of both approaches thanks to the expected symmetry between both hemispheres; 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. Thus, we only show the runs performed with a sharp transition at θ = π/2, or Ψ = Ψ_{⋆}. In the PIC domain, electronpositron pairs are continuously injected with a high multiplicity, κ_{inj} = 10, at the star surface to ensure a forcefree 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.
The top panel of Fig. 2 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 of 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 forcefree 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 forcefree 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
Fig. 2. Simulated monopolar magnetosphere using the hybrid forcefreePIC scheme. The division between the forcefree and PIC domains is the equatorial plane (dashed magenta line in the top and middle panels). Top: Spatial distribution of the radial current density normalized to the fiducial current density, ${J}_{\mathrm{GJ}}^{\star}$, and compensated by (r/r_{⋆})^{2} for visualization purposes, for the 4096^{2} cells run. Solid lines show the poloidal magnetic field lines. Middle: Crosssection of the current density profile at r = 1.5R_{LC} as a function of the numerical resolution, and comparison with the exact profile (red dasheddotted line, Eq. 34). Bottom: Radial profiles of the Poynting flux, L(r), normalized to the monopole solution, L_{M} (Eq. 35), for all three resolutions. Percentages correspond to the amount of numerical dissipation for each run, computed between the flux averaged in the gray areas. 
$$\begin{array}{c}\hfill \u03f5=\frac{{L}_{\mathrm{in}}{L}_{\mathrm{out}}}{{L}_{\mathrm{out}}},\end{array}$$(36)
where L_{in} is the Poynting flux averaged in the r ∈ [r_{⋆}, 2r_{⋆}] interval, and L_{out} is averaged over the r ∈ [2R_{LC}, 3R_{LC}] interval (gray intervals in Fig. 2, bottom panel). We observe numerical convergence and a deviation smaller than ϵ < 1% for the 1024^{2} run, down to ϵ ∼ 0.15% for the highest resolution. We note that numerical dissipation mainly originates from the forcefree 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 models 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 forcefree 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 gammaray pulsars introduced in Sect. 2.
5.1. Numerical setup
The properties of the simulation domain are nearly identical to 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 8192^{2} 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
$$\begin{array}{c}\hfill {\mathrm{\Psi}}_{\mathrm{D}}=\frac{\mu {sin}^{2}\theta}{r}.\end{array}$$(37)
There is no exact analytical solution in the forcefree limit for an aligned dipole, and thus results are compared with previous forcefree 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, m_{p} = 1836m_{e}. This processes can be faithfully reproduced by maintaining the GoldreichJulian 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 pairproduction 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 γ_{th}, then a new pair is created with a Lorentz factor γ_{s} 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 forcefree domain, there is no other condition to consider apart from corotation at the neutron star surface.
Our choice for an efficient forcefreePIC domain decomposition in this setup is educated from the magnetic topology revealed by previous studies. For a plasmafilled magnetosphere, the main loci of dissipation are the equatorial current sheet and at its base around the Ypoint 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 gammaray emission in observed pulsars that we seek to reproduce in this work. A natural choice is thus 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 forcefree 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 defeat the purpose to choose a wide PIC domain. We found that a good compromise is to use the following magnetic flux function boundaries:

Forcefree domain I: 0 < Ψ < 0.9Ψ_{pc}. This region corresponds to the open field lines connected to the star polar caps.

PIC domain: 0.85Ψ_{pc} < Ψ < 2.4Ψ_{pc}. The domain encompasses the equatorial and separatrix current layers.

Forcefree domain II: 2.3Ψ_{pc} < Ψ < Ψ_{max}. It is inside the closed field line regions deep inside the light cylinder.
In the above, Ψ_{max} = μ/r_{⋆} is the maximum of the flux function that is located at the star’s surface in the equatorial plane, and Ψ_{pc} = μ/R_{LC} 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 forcefree dipole. There is a ΔΨ = 0.1Ψ_{pc}thick transition layer between the PIC and forcefree domains. However, during the transient phase, we found that the transition layer between the forcefree domain I and the PIC domain can sometimes partially prevent field line opening, so 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 r_{⋆} = 10^{6} cm, and the lightcylinder radius is set to R_{LC} = 5 × 10^{6} cm, corresponding to a P = 1 ms spin period. The surface magnetic field of the reference production run is fixed at B_{⋆} = 10^{7} G, giving B_{LC} = μ/R_{LC}^{3} ≈ 8 × 10^{4} G at the light cylinder, and a forcefree spindown power of L_{0} = μ^{2}Ω^{4}/c^{3} = 4.8 × 10^{33} erg/s. The magnetic field strength then sets the relevant energy scales given in Sect. 2 as γ_{pc} = 2.6 × 10^{8}, γ_{th} = 10^{6}, γ_{s} = 2.2 × 10^{5}, and γ_{rad} = 3 × 10^{4}. 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 skin depth and frequency scales are resolved by ${d}_{\mathrm{e}}^{\mathrm{s}}/\mathrm{\Delta}r\approx 3$ and ${\omega}_{\mathrm{pe}}^{1}/\mathrm{\Delta}t\approx 8$ 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, ρ_{s}, divided by the local grid spacing, R_{LC}Δr/r_{⋆}, is ρ_{s}/(R_{LC}Δr/r_{⋆}) = (γ_{s}m_{e}c^{2}/eB_{LC})/(R_{LC}Δr/r_{⋆})≈3, and Larmor timescale ρ_{s}/(cΔt)≈36. The corresponding synchrotron cooling time is on the order of ${t}_{\mathrm{syn}}^{\mathrm{LC}}\sim 3{m}_{\mathrm{e}}c/(2{r}_{\mathrm{e}}^{2}{\gamma}_{\mathrm{s}}{B}_{\mathrm{LC}}^{2})\approx 85\mathrm{\Delta}t$.
To prevent spurious numerical effects near the star’s 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 an r_{⋆}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 radiationreaction force is set to 1% of its nominal value, and it then slowly increases in time up to about 10%. In a last step, it is gradually increased up to the nominal strength at t ≳ 2P when a quasisteady 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 feasible. Here, we scaled the magnetic field strength of the reference simulation down to B_{⋆} = 10^{4} G, 10^{5} G, and 10^{6} G, but we kept the ratio between all energy scales the same as in the B_{⋆} = 10^{7} G simulation, including the radiationreactionlimited lepton Lorentz factor, γ_{rad}, given in Eq. (14). To keep the γ_{pc}/γ_{rad} ratio identical to that of the B_{⋆} = 10^{7} G run, the strength of the radiation reaction force –or equivalently the Thomson crosssection– must be boosted by a factor of f_{rad} = 10^{9} for B_{⋆} = 10^{4} G, f_{rad} = 10^{6} for B_{⋆} = 10^{5}, and f_{rad} = 10^{3} for B_{⋆} = 10^{6} G. Table 1 lists the numerical and physical parameters used in this study for the reference simulation.
Physical and numerical parameters for the reference weak gammaray pulsar simulation using the hybrid forcefreePIC scheme.
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 fullscale radiationreaction force is established at t ≈ 2.5P for the fiducial fullscale B_{⋆} = 10^{7} 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 Ypoint). It is important to note that this solution has not completely reached a steady state. During the initial transient phase, the Ypoint forms well inside the light cylinder (around 0.6R_{LC}), slowly migrates toward R_{LC}, and continues to do so at this stage. This migration is still visible in Figure 3 in the form of density striations inside the Ypoint, 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. 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 layerthickness scales to stellarradii scales for the largest ones. The layer is also slightly kinked due to currentdriven 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 us to bridge the pure forcefree region to these vacuum gaps without noticeable numerical artifacts. We 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 pairproducing photons have virtually zero meanfree path. However, we believe that this has no dramatic effect on the highenergy output of the magnetosphere described below.
Fig. 3. Global hybrid forcefreePIC model of the aligned magnetosphere of a B_{⋆} = 10^{7} G surface magnetic field neutron star with a P = 1 ms spin period at simulated time, t = 2.45 ms. Solid black lines show the poloidal magnetic field lines (i.e., isocontour of the magnetic flux function Ψ), solid green and red lines show the transition from the forcefree regions to the PIC region (i.e., Ψ_{0}, Ψ_{1}, Ψ_{2}, and Ψ_{3}; see Table 1). Panel a: Map of the pairs density normalized by ${n}_{\mathrm{GJ}}^{\star}$ and compensated by (r/r_{⋆})^{2}. Panel b: Map of the mean pair energy, ℰ = γm_{e}c^{2}. Panels c and d: Same as panels (a) and (b), but for the protons. The panels do not show the full domain, but instead they are focused on the PIC domain where the current sheet forms. A zoomedin view of a reconnection site is shown in panel (a). Notice how thin the reconnection layer is. 
Perhaps the most spectacular feature is the extreme narrowness of the current layer, the latter being barely visible in the global view (see the zoomedin view in Figure 3). Near the light cylinder, the layer’s thickness is about δ ∼ 10^{−3}R_{LC}, 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 κ ∼ 10 − 100, while the proton density is at most of order unity. The pairs being much more abundant than the protons, the layer’s thickness is controlled by the electronic scales. In contrast, the proton beam is much thicker, with a width of ∼0.1R_{LC} 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, from which the protons are nearly immune.
Figure 3 also maps out the mean particle energies in the magnetosphere, demonstrating the crucial role of reconnection in particle acceleration. Xpoints, 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 counterstreaming beams at Xpoints. Positrons (and protons) move radially outwards along with the wind and plasmoids, while electrons tend to precipitate inwards. Counterstreaming 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 zoomedin view on one of these Xpoints present in the layer (see the inset in Figure 3) shows that the energy of the pairs peaks inside the layer, ℰ_{pairs} = γm_{e}c^{2} ≈ 1 TeV, and sharply decreases inside the nearest plasmoids on either side, down to ℰ_{pairs} ≲ 100 GeV. In this extreme synchrotron cooling regime, particles accelerate deep inside the layer where the effective magneticfield strength is low, allowing them to accelerate beyond the radiationreactionlimited energy, ${\gamma}_{\mathrm{rad}}^{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}\approx 150$ GeV≪ℰ_{pairs}; however, 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 they 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 nonradiative reconnection (Petropoulou & Sironi 2018; Zhang et al. 2021). However, plasmoids, Xpoints, 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 of ℰ ≈ 10 TeV, which is close to the full polarcap potential drop (Philippov & Spitkovsky 2018; Guépin et al. 2020). This energy scale is consistent with systemsizelimited acceleration achievable with reconnection, given by ℰ_{max} ∼ β_{rec}eΦ_{pc} ≈ β_{rec} × 130 TeV, where β_{rec} ∼ 0.1 − 0.2 is the reconnection rate.
Figure 4 shows the energy distribution for all the chargedparticle species. The pair spectrum ranges from about 1 GeV to 1 TeV and peaks at about a few tens of gigaelectronvolt, below the radiationreactionlimited energy (${\gamma}_{\mathrm{rad}}^{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}\approx 150$ 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 at 300–400 GeV, which coincides with particles energized at Xpoints whose energy is determined by the magnetization parameter σ_{LC}m_{e}c^{2} ≈ 500 GeV (Eq. 9). This feature was already reported in previous works (Cerutti et al. 2015), and it seems to hold strong on 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 Ω ⋅ μ > 0. We 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 gammaray pulsar because of the proximity between the relevant energy scales. The proton spectrum is composed of a single component ranging from about 1 TeV up to 30 TeV. If we venture to fit a power law, the spectral index would be ∼ − 1.5, which is 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 polarcap potential drop, γ_{pc}, 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.
Fig. 4. Total electron (solid blue), positron (solid red), and proton (dashed green) energy spectra for a B_{⋆} = 10^{7} G and P = 1 ms pulsar. The relevant energy scales introduced in Sect. 2 are shown by vertical dotted lines, from the lowest to the highest energies: the radiationreactionlimited energy at the light cylinder, ${\gamma}_{\mathrm{rad}}^{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}$, the threshold energy for pair production, γ_{th}m_{e}c^{2}, and the polarcap potential drop, eΦ_{pc}. 
Fig. 5. Maximum cutoff energy of the charged particle energy distributions, ℰ_{max} (electrons: blue dots; positrons: red dots; protons: green dots), as a function of the surface magnetic field B_{⋆}. The evolution of the polarcap potential drop is shown by the dashed line for comparison (Eq. 4). 
The narrowness of the pair energy spectrum translates into a narrowband synchrotron energy distribution. The total synchrotron spectrum reported in Figure 6 ranges from about 1 MeV to 1 GeV, with a peak in the energy band of the FermiLAT above 100 MeV. 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 (≳R_{LC}). 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 Xpoints beyond the synchrotron burnoff limit, given by ϵ_{rad} ∼ (9/4)m_{e}c^{2}/α ≈ 160 MeV (Uzdensky et al. 2011) where α is the fine structure constant. Integrating over all frequencies, the synchrotron luminosity accounts for about 18% of the pulsar spindown. This sizeable fraction corresponds to the amount of dissipation of the Poynting flux measured in the magnetosphere within a few lightcylinder 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 (≲0.1%), 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 nonthermal radiation, which is possible only in the strong cooling regime. Indeed, about 2% of the spindown is channelled into the beam of protons while the pairs leave the box with ≲1% 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 FermiLAT, our fiducial pulsar shines about 4% of its total spindown power above 100 MeV range. Compared to the Fermi pulsar population, our solution is compatible with the typical reported gammaray efficiencies that range from 1–100% (Smith et al. 2023; see Figure 7).
Fig. 6. Total synchrotron spectral energy distribution emitted by the equatorial current sheet for B_{⋆} = 10^{7} G and P = 1 ms (dashed black line). The electronic (positronic) contribution is shown by the blue (red) line. The vertical dotted line indicates the 100 MeV energy threshold of the FermiLAT for comparison. 
Fig. 7. Pulsar gammaray efficiency above 100 MeV; that is, the gammaray luminosity over the pulsar spindown power (L_{γ}/L), as reported in the third FermiLAT catalog (gray stars; Smith et al. 2023), as a function of the pulsar spindown power (L). The properties of the pulsar simulated in this work (B_{⋆} = 10^{7} G, P = 1 ms) are shown by the red dot. The error bars give a measure of the variations of the spindown and radiative efficiencies measured in a 0.5P time frame. 
6. Conclusion
In this work, we developed a new hybrid scheme designed to blend the PIC and the forcefree approaches in the context of relativistic magnetospheres. This effort was motivated by the need to increase the separation between the kinetic scales, where particle acceleration operates, and global scales. By this means, we increase the credibility of global PIC models in spite of their unrealistically low scale separations. Using these new capabilities and highresolution simulations, we focused our numerical resources on the modeling of a weak millisecond gammaray 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 us 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 a 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, B_{⋆} = 10^{7} G, simulated in this work, pairs are accelerated up to 1 TeV, while protons reach up to ≳10% of the full polarcap potential drop, that is ≳10 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 cosmicray 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 20% of the spindown power, including about 4% above 100 MeV gamma rays, a number compatible with the FermiLAT 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 2023).
We also show 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 intended to provide a proof of concept that magnetic reconnection within the equatorial current sheet is a viable scenario to explain the gammaray emission, without scaling down the physical parameters, at least in the lowluminosity range of the gammaray pulsar population. In this respect, our fiducial simulation has successfully achieved this objective, although much remains to be done in the future. 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 or photon spectral shapes and cutoff energies, as well as the radiative efficiencies. These last two quantities reported in this work are low compared to the bulk of the FermiLAT millisecond pulsar population. Longer integration times to establish a more robust steady state, in particular inside the light cylinder (separatrices and Ypoint), 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.
Acknowledgments
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, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460 [Google Scholar]
 Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17 [Google Scholar]
 Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arka, I., & Dubus, G. 2013, A&A, 550, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arons, J. 2003, ApJ, 589, 871 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Spitkovsky, A. 2010, ApJ, 715, 1282 [Google Scholar]
 Belyaev, M. A. 2015, MNRAS, 449, 2759 [NASA ADS] [CrossRef] [Google Scholar]
 Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (Boca Raton: CRC Press) [CrossRef] [Google Scholar]
 Blandford, R. D. 2002, in Lighthouses of the Universe: The Most Luminous Celestial Objects and Their Use for Cosmology, eds. M. Gilfanov, R. Sunyeav, & E. Churazov, 381 [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Blasi, P., Epstein, R. I., & Olinto, A. V. 2000, ApJ, 533, L123 [NASA ADS] [CrossRef] [Google Scholar]
 Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
 Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Cao, G., & Yang, X. 2019, ApJ, 874, 166 [Google Scholar]
 Cao, G., Zhang, L., & Sun, S. 2016, MNRAS, 455, 4267 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., & Beloborodov, A. M. 2017, Space Sci. Rev., 207, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., & Philippov, A. A. 2017, A&A, 607, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [Google Scholar]
 Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, ApJ, 889, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Chernoglazov, A., Hakobyan, H., & Philippov, A. 2023, ApJ, 959, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, I. 2007, A&A, 472, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Contopoulos, I., & Stefanou, P. 2019, MNRAS, 487, 952 [CrossRef] [Google Scholar]
 Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351 [Google Scholar]
 Contopoulos, I., Pétri, J., & Stefanou, P. 2020, MNRAS, 491, 5579 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. N. 2024, MNRAS, 527, L127 [Google Scholar]
 Coroniti, F. V. 1990, ApJ, 349, 538 [Google Scholar]
 Crinquand, B., Cerutti, B., & Dubus, G. 2019, A&A, 622, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cruz, F., Grismayer, T., Chen, A. Y., Spitkovsky, A., & Silva, L. O. 2021, ApJ, 919, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, F., Grismayer, T., Chen, A. Y., et al. 2024, A&A, in press, https://doi.org/10.1051/00046361/202347926 [Google Scholar]
 Erber, T. 1966, Rev. Mod. Phys., 38, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869 [Google Scholar]
 Gruzinov, A. 1999, ArXiv eprints [arXiv: astroph/9902288] [Google Scholar]
 Gruzinov, A. 2008, ArXiv eprints [arXiv:0802.1716] [Google Scholar]
 Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gunn, J. E., & Ostriker, J. P. 1969, Phys. Rev. Lett., 22, 728 [NASA ADS] [CrossRef] [Google Scholar]
 H.E.S.S. Collaboration (Aharonian, F., et al.) 2023, Nat. Astron., 7, 1341 [NASA ADS] [CrossRef] [Google Scholar]
 Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, 943, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Harding, A. K., & Lai, D. 2006, Rep. Prog. Phys., 69, 2631 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Hones, E. W., Jr, & Bergeson, J. E. 1965, J. Geophys. Res., 70, 4951 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, R., & Beloborodov, A. M. 2022, ApJ, 939, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Kalapotharakos, C., Harding, A. K., & Kazanas, D. 2014, ApJ, 793, 97 [Google Scholar]
 Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261 [Google Scholar]
 Komissarov, S. S. 2002, MNRAS, 336, 759 [Google Scholar]
 Komissarov, S. S. 2006, MNRAS, 367, 19 [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1971, in The Classical Theory of Fields, eds. L. D. Landau, & E. M. Lifshitz (New York: Pergamon Press) [Google Scholar]
 Li, J., Spitkovsky, A., & Tchekhovskoy, A. 2012, ApJ, 746, 60 [Google Scholar]
 Lyubarsky, Y. 2020, ApJ, 897, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lyubarsky, Y., & Kirk, J. G. 2001, ApJ, 547, 437 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C. 2006, MNRAS, 368, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Melzani, M., Walder, R., Folini, D., Winisdoerffer, C., & Favre, J. M. 2014, A&A, 570, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michel, F. C. 1973, ApJ, 180, L133 [Google Scholar]
 Michel, F. C., & Tucker, W. H. 1969, Nature, 223, 277 [CrossRef] [Google Scholar]
 Most, E. R., & Philippov, A. A. 2020, ApJ, 893, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Parfrey, K., Beloborodov, A. M., & Hui, L. 2012, MNRAS, 423, 1416 [CrossRef] [Google Scholar]
 Petropoulou, M., & Sironi, L. 2018, MNRAS, 481, 5687 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A., & Kramer, M. 2022, ARA&A, 60, 495 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., & Spitkovsky, A. 2014, ApJ, 785, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJ, 815, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015b, ApJ, 801, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, Phys. Rev. Lett., 124, 245101 [Google Scholar]
 Popov, S. B., & Postnov, K. A. 2010, in Evolution of Cosmic Objects through their Physical Activity, eds. H. A. Harutyunian, A. M. Mickaelian, & Y. Terzian, 129 [Google Scholar]
 Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51 [Google Scholar]
 Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951 [Google Scholar]
 Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Speiser, T. W. 1965, J. Geophys. Res., 70, 4219 [NASA ADS] [CrossRef] [Google Scholar]
 Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
 Sturrock, P. A. 1971, ApJ, 164, 529 [Google Scholar]
 Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255 [Google Scholar]
 Timokhin, A. N. 2006, MNRAS, 368, 1055 [Google Scholar]
 Timokhin, A. N., & Arons, J. 2013, MNRAS, 429, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, R., Grismayer, T., Cruz, F., Fonseca, R. A., & Silva, L. O. 2024, New A, 112, 102261 [NASA ADS] [CrossRef] [Google Scholar]
 Uzdensky, D. A. 2003, ApJ, 598, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Uzdensky, D. A., & Spitkovsky, A. 2014, ApJ, 780, 3 [NASA ADS] [Google Scholar]
 Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJ, 737, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [Google Scholar]
 Yee, K. 1966, IEEE Trans. Antennas Propagation, 14, 302 [CrossRef] [Google Scholar]
 Zhang, H., Sironi, L., & Giannios, D. 2021, ApJ, 922, 261 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Physical and numerical parameters for the reference weak gammaray pulsar simulation using the hybrid forcefreePIC scheme.
All Figures
Fig. 1. New hybrid numerical scheme proposed in this work during a timeintegration cycle. This method is intended to blend the forcefree electrodynamic (FFE) and the full kinetic (PIC) approaches within the same numerical framework using a domain decomposition. The coupling is achieved via the electric current densities, J, using a blending function, f, which solely depends on the magnetic flux function, Ψ. 

In the text 
Fig. 2. Simulated monopolar magnetosphere using the hybrid forcefreePIC scheme. The division between the forcefree and PIC domains is the equatorial plane (dashed magenta line in the top and middle panels). Top: Spatial distribution of the radial current density normalized to the fiducial current density, ${J}_{\mathrm{GJ}}^{\star}$, and compensated by (r/r_{⋆})^{2} for visualization purposes, for the 4096^{2} cells run. Solid lines show the poloidal magnetic field lines. Middle: Crosssection of the current density profile at r = 1.5R_{LC} as a function of the numerical resolution, and comparison with the exact profile (red dasheddotted line, Eq. 34). Bottom: Radial profiles of the Poynting flux, L(r), normalized to the monopole solution, L_{M} (Eq. 35), for all three resolutions. Percentages correspond to the amount of numerical dissipation for each run, computed between the flux averaged in the gray areas. 

In the text 
Fig. 3. Global hybrid forcefreePIC model of the aligned magnetosphere of a B_{⋆} = 10^{7} G surface magnetic field neutron star with a P = 1 ms spin period at simulated time, t = 2.45 ms. Solid black lines show the poloidal magnetic field lines (i.e., isocontour of the magnetic flux function Ψ), solid green and red lines show the transition from the forcefree regions to the PIC region (i.e., Ψ_{0}, Ψ_{1}, Ψ_{2}, and Ψ_{3}; see Table 1). Panel a: Map of the pairs density normalized by ${n}_{\mathrm{GJ}}^{\star}$ and compensated by (r/r_{⋆})^{2}. Panel b: Map of the mean pair energy, ℰ = γm_{e}c^{2}. Panels c and d: Same as panels (a) and (b), but for the protons. The panels do not show the full domain, but instead they are focused on the PIC domain where the current sheet forms. A zoomedin view of a reconnection site is shown in panel (a). Notice how thin the reconnection layer is. 

In the text 
Fig. 4. Total electron (solid blue), positron (solid red), and proton (dashed green) energy spectra for a B_{⋆} = 10^{7} G and P = 1 ms pulsar. The relevant energy scales introduced in Sect. 2 are shown by vertical dotted lines, from the lowest to the highest energies: the radiationreactionlimited energy at the light cylinder, ${\gamma}_{\mathrm{rad}}^{\mathrm{LC}}{m}_{\mathrm{e}}{c}^{2}$, the threshold energy for pair production, γ_{th}m_{e}c^{2}, and the polarcap potential drop, eΦ_{pc}. 

In the text 
Fig. 5. Maximum cutoff energy of the charged particle energy distributions, ℰ_{max} (electrons: blue dots; positrons: red dots; protons: green dots), as a function of the surface magnetic field B_{⋆}. The evolution of the polarcap potential drop is shown by the dashed line for comparison (Eq. 4). 

In the text 
Fig. 6. Total synchrotron spectral energy distribution emitted by the equatorial current sheet for B_{⋆} = 10^{7} G and P = 1 ms (dashed black line). The electronic (positronic) contribution is shown by the blue (red) line. The vertical dotted line indicates the 100 MeV energy threshold of the FermiLAT for comparison. 

In the text 
Fig. 7. Pulsar gammaray efficiency above 100 MeV; that is, the gammaray luminosity over the pulsar spindown power (L_{γ}/L), as reported in the third FermiLAT catalog (gray stars; Smith et al. 2023), as a function of the pulsar spindown power (L). The properties of the pulsar simulated in this work (B_{⋆} = 10^{7} G, P = 1 ms) are shown by the red dot. The error bars give a measure of the variations of the spindown and radiative efficiencies measured in a 0.5P time frame. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.