11institutetext: Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
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

Adrien Soudais \XeTeXLinkBox Benoît Cerutti\XeTeXLinkBox 1111    Ioannis Contopoulos\XeTeXLinkBox 22
(Received 04 April 2024 / Accepted 12 June 2024)
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, outflows

1 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 P=1𝑃1P=1italic_P = 1ms and a surface dipolar magnetic field B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG, corresponding to the weakest observable gamma-ray pulsars of spindown power L4.8×1033𝐿4.8superscript1033L\approx 4.8\times 10^{33}italic_L ≈ 4.8 × 10 start_POSTSUPERSCRIPT 33 end_POSTSUPERSCRIPTerg/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, r106similar-tosubscript𝑟superscript106r_{\star}\sim 10^{6}italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ∼ 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPTcm, and the light-cylinder radius,

RLC=cP2π5r(P1ms).subscript𝑅LC𝑐𝑃2𝜋5subscript𝑟𝑃1msR_{\rm LC}=\frac{cP}{2\pi}\approx 5r_{\star}\left(\frac{P}{1\rm{ms}}\right).italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT = divide start_ARG italic_c italic_P end_ARG start_ARG 2 italic_π end_ARG ≈ 5 italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ( divide start_ARG italic_P end_ARG start_ARG 1 roman_m roman_s end_ARG ) . (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, n=κnGJ=κB/2πeRLCsubscript𝑛𝜅subscriptsuperscript𝑛GJ𝜅subscript𝐵2𝜋𝑒subscript𝑅LCn_{\star}=\kappa n^{\star}_{\rm GJ}=\kappa B_{\star}/2\pi eR_{\rm LC}italic_n start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = italic_κ italic_n start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT = italic_κ italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / 2 italic_π italic_e italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT (Hones & Bergeson, 1965; Goldreich & Julian, 1969), where κ𝜅\kappaitalic_κ is the plasma multiplicity, and e𝑒eitalic_e is the electron charge, yields

desubscriptsuperscript𝑑e\displaystyle d^{\star}_{\rm e}italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT =\displaystyle== γmec24πne2𝛾subscript𝑚esuperscript𝑐24𝜋subscript𝑛superscript𝑒2\displaystyle\sqrt{\frac{\gamma m_{\rm e}c^{2}}{4\pi n_{\star}e^{2}}}square-root start_ARG divide start_ARG italic_γ italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 italic_π italic_n start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG (2)
\displaystyle\approx 2(γ100)1/2(κ102)1/2(P1ms)1/2(B107G)1/2cm,2superscript𝛾superscript10012superscript𝜅superscript10212superscript𝑃1ms12superscriptsubscript𝐵superscript107G12cm\displaystyle 2\left(\frac{\gamma}{10^{0}}\right)^{1/2}\left(\frac{\kappa}{10^% {2}}\right)^{-1/2}\left(\frac{P}{1\rm{ms}}\right)^{1/2}\left(\frac{B_{\star}}{% 10^{7}\rm{G}}\right)^{-1/2}\rm{cm},2 ( divide start_ARG italic_γ end_ARG start_ARG 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_κ end_ARG start_ARG 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_P end_ARG start_ARG 1 roman_m roman_s end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G end_ARG ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT roman_cm ,

where γ𝛾\gammaitalic_γ is the particle Lorentz factor and mesubscript𝑚em_{\rm e}italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT is the electron mass, so that de/r=2×106subscriptsuperscript𝑑esubscript𝑟2superscript106d^{\star}_{\rm e}/r_{\star}=2\times 10^{-6}italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 2 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 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 (γ1much-greater-than𝛾1\gamma\gg 1italic_γ ≫ 1), 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 sinθpc=r/RLCsubscript𝜃pcsubscript𝑟subscript𝑅LC\sin\theta_{\rm pc}=\sqrt{r_{\star}/R_{\rm LC}}roman_sin italic_θ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = square-root start_ARG italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG. Using the corotation surface electric field, 𝐄=(𝛀×𝐫)×𝐁/c𝐄𝛀𝐫𝐁𝑐\mathbf{E}=-(\boldsymbol{\Omega}\times\mathbf{r})\times\mathbf{B}/cbold_E = - ( bold_Ω × bold_r ) × bold_B / italic_c, we can derive the potential drop across the polar cap as

Φpc=μΩ2c2,subscriptΦpc𝜇superscriptΩ2superscript𝑐2\Phi_{\rm pc}=\frac{\mu\Omega^{2}}{c^{2}},roman_Φ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = divide start_ARG italic_μ roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (3)

where Ω=2π/PΩ2𝜋𝑃\Omega=2\pi/Proman_Ω = 2 italic_π / italic_P is the angular velocity of the star, and μ=Br3𝜇subscript𝐵subscriptsuperscript𝑟3\mu=B_{\star}r^{3}_{\star}italic_μ = italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT its magnetic moment. An electron experiencing the full potential drop will acquire a Lorentz factor given by

γpc=eΦpcmec22.6×108(B107G)(P1ms)2.subscript𝛾pc𝑒subscriptΦpcsubscript𝑚esuperscript𝑐22.6superscript108subscript𝐵superscript107Gsuperscript𝑃1ms2\gamma_{\rm pc}=\frac{e\Phi_{\rm pc}}{m_{\rm e}c^{2}}\approx 2.6\times 10^{8}% \left(\frac{B_{\star}}{10^{7}\rm{G}}\right)\left(\frac{P}{1\rm{ms}}\right)^{-2}.italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = divide start_ARG italic_e roman_Φ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT end_ARG start_ARG italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ≈ 2.6 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT ( divide start_ARG italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G end_ARG ) ( divide start_ARG italic_P end_ARG start_ARG 1 roman_m roman_s end_ARG ) start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT . (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)

χϵb0.1,𝜒italic-ϵ𝑏greater-than-or-equivalent-to0.1\chi\equiv\epsilon b\gtrsim 0.1,italic_χ ≡ italic_ϵ italic_b ≳ 0.1 , (5)

where ϵ=ν/mec2italic-ϵPlanck-constant-over-2-pi𝜈subscript𝑚esuperscript𝑐2\epsilon=\hbar\nu/m_{\rm e}c^{2}italic_ϵ = roman_ℏ italic_ν / italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT is the dimensionless photon energy, and b=B~/BQED𝑏subscript~𝐵perpendicular-tosubscript𝐵QEDb=\tilde{B}_{\perp}/B_{\rm QED}italic_b = over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT / italic_B start_POSTSUBSCRIPT roman_QED end_POSTSUBSCRIPT is the effective perpendicular magnetic field (see Eq. 17 below) normalized to the critical quantum field BQED=me2c3/e4.4×1013subscript𝐵QEDsubscriptsuperscript𝑚2esuperscript𝑐3Planck-constant-over-2-pi𝑒4.4superscript1013B_{\rm QED}=m^{2}_{\rm e}c^{3}/\hbar e\approx 4.4\times 10^{13}italic_B start_POSTSUBSCRIPT roman_QED end_POSTSUBSCRIPT = italic_m start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT / roman_ℏ italic_e ≈ 4.4 × 10 start_POSTSUPERSCRIPT 13 end_POSTSUPERSCRIPTG. If electrons radiate via synchrotron-curvature radiation as usually assumed, the critical photon energy is given by (Blumenthal & Gould, 1970)

ϵc=32bγ2.subscriptitalic-ϵc32𝑏superscript𝛾2\epsilon_{\rm c}=\frac{3}{2}b\gamma^{2}.italic_ϵ start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT = divide start_ARG 3 end_ARG start_ARG 2 end_ARG italic_b italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (6)

Combining Eq. (5) with Eq. (6) provides an estimate for the threshold electron Lorentz factor capable of producing new pairs,

γth=115b2106(B107G)1,subscript𝛾th115superscript𝑏2superscript106superscriptsubscript𝐵superscript107G1\gamma_{\rm th}=\sqrt{\frac{1}{15b^{2}}}\approx 10^{6}\left(\frac{B_{\star}}{1% 0^{7}\rm{G}}\right)^{-1},italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG 1 end_ARG start_ARG 15 italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG ≈ 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT ( divide start_ARG italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (7)

while the created pair shares the absorbed photon momentum, such that the Lorentz factor of the secondary pairs may be estimated as

γs=ϵ2=120b2.2×105(B107G)1.subscript𝛾sitalic-ϵ2120𝑏2.2superscript105superscriptsubscript𝐵superscript107G1\gamma_{\rm s}=\frac{\epsilon}{2}=\frac{1}{20b}\approx 2.2\times 10^{5}\left(% \frac{B_{\star}}{10^{7}\rm{G}}\right)^{-1}.italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = divide start_ARG italic_ϵ end_ARG start_ARG 2 end_ARG = divide start_ARG 1 end_ARG start_ARG 20 italic_b end_ARG ≈ 2.2 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT ( divide start_ARG italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (8)

Assuming that pair creation is effective at producing a high-multiplicity plasma in the magnetosphere (κ1much-greater-than𝜅1\kappa\gg 1italic_κ ≫ 1), it sets the relevant scale of the particle Lorentz factor to be considered in estimating the plasma skindepth scale in Eq. (2), yielding des=de(γs)103cm=103rsubscriptsuperscript𝑑sesubscriptsuperscript𝑑esubscript𝛾ssuperscript103cmsuperscript103subscriptrd^{\rm s}_{\rm e}=d^{\star}_{\rm e}(\gamma_{\rm s})\approx 10^{3}\rm{cm}=10^{-% 3}r_{\star}italic_d start_POSTSUPERSCRIPT roman_s end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT = italic_d start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT ( italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT ) ≈ 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT roman_cm = 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT roman_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT for our reference pulsar. One can realize from these order of magnitude calculations that a weaker magnetic field (B107less-than-or-similar-tosubscript𝐵superscript107B_{\star}\lesssim 10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ≲ 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG) 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, δ𝛿\deltaitalic_δ, is governed by the local plasma skindepth scale, meaning that δdeLCsimilar-to𝛿subscriptsuperscript𝑑LCe\delta\sim d^{\rm LC}_{\rm e}italic_δ ∼ italic_d start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT. 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,

γLCσLC=BLC24πΓLCκnGJLCmec2=γpc2ΓLCκ,similar-tosubscript𝛾LCsubscript𝜎LCsubscriptsuperscript𝐵2LC4𝜋subscriptΓLC𝜅subscriptsuperscript𝑛LCGJsubscript𝑚esuperscript𝑐2subscript𝛾pc2subscriptΓLC𝜅\gamma_{\rm LC}\sim\sigma_{\rm LC}=\frac{B^{2}_{\rm LC}}{4\pi\Gamma_{\rm LC}% \kappa n^{\rm LC}_{\rm GJ}m_{\rm e}c^{2}}=\frac{\gamma_{\rm pc}}{2\Gamma_{\rm LC% }\kappa},italic_γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ∼ italic_σ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT = divide start_ARG italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG start_ARG 4 italic_π roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_κ italic_n start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = divide start_ARG italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT end_ARG start_ARG 2 roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_κ end_ARG , (9)

where BLCB(r/RLC)3=μΩ3/c3similar-tosubscript𝐵LCsubscript𝐵superscriptsubscript𝑟subscript𝑅LC3𝜇superscriptΩ3superscript𝑐3B_{\rm LC}\sim B_{\star}(r_{\star}/R_{\rm LC})^{3}=\mu\Omega^{3}/c^{3}italic_B start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ∼ italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT = italic_μ roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT / italic_c start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, nGJLC=ΩBLC/2πecsubscriptsuperscript𝑛LCGJΩsubscript𝐵LC2𝜋𝑒𝑐n^{\rm LC}_{\rm GJ}=\Omega B_{\rm LC}/2\pi ecitalic_n start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT = roman_Ω italic_B start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT / 2 italic_π italic_e italic_c is the Goldreich-Julian plasma density at the light cylinder, and ΓLCsubscriptΓLC\Gamma_{\rm LC}roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT is the pulsar wind bulk Lorentz factor at its base. Assuming that these energetic particles set the local skindepth scale yields

δdeLC(γLC)=γLCmec24πκΓLCnGJLCe2=RLC2ΓLCκ,similar-to𝛿subscriptsuperscript𝑑LCesubscript𝛾LCsubscript𝛾LCsubscript𝑚esuperscript𝑐24𝜋𝜅subscriptΓLCsubscriptsuperscript𝑛LCGJsuperscript𝑒2subscript𝑅LC2subscriptΓLC𝜅\delta\sim d^{\rm LC}_{\rm e}(\gamma_{\rm LC})=\sqrt{\frac{\gamma_{\rm LC}m_{% \rm e}c^{2}}{4\pi\kappa\Gamma_{\rm LC}n^{\rm LC}_{\rm GJ}e^{2}}}=\frac{R_{\rm LC% }}{2\Gamma_{\rm LC}\kappa},italic_δ ∼ italic_d start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT ( italic_γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ) = square-root start_ARG divide start_ARG italic_γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 italic_π italic_κ roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_n start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG = divide start_ARG italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG start_ARG 2 roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_κ end_ARG , (10)
δ2.4×104(ΓLC100)1(κ102)1(P1ms)cm,similar-to𝛿2.4superscript104superscriptsubscriptΓLCsuperscript1001superscript𝜅superscript1021𝑃1mscm\delta\sim 2.4\times 10^{4}\left(\frac{\Gamma_{\rm LC}}{10^{0}}\right)^{-1}% \left(\frac{\kappa}{10^{2}}\right)^{-1}\left(\frac{P}{1\rm{ms}}\right)\rm{cm},italic_δ ∼ 2.4 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( divide start_ARG roman_Γ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( divide start_ARG italic_κ end_ARG start_ARG 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( divide start_ARG italic_P end_ARG start_ARG 1 roman_m roman_s end_ARG ) roman_cm , (11)

thus, δ/r2×102similar-to𝛿subscript𝑟2superscript102\delta/r_{\star}\sim 2\times 10^{-2}italic_δ / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ∼ 2 × 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT or δ/RLC5×103similar-to𝛿subscript𝑅LC5superscript103\delta/R_{\rm LC}\sim 5\times 10^{-3}italic_δ / italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ∼ 5 × 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT 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, eEfradsimilar-to𝑒subscript𝐸parallel-tosubscript𝑓radeE_{\parallel}\sim f_{\rm rad}italic_e italic_E start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT ∼ italic_f start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT. The latter effectively acts as a continuous drag force (in the classical regime, i.e., γb1much-less-than𝛾𝑏1\gamma b\ll 1italic_γ italic_b ≪ 1), opposite to the particle direction of motion, whose magnitude in the ultrarelativistic regime (γ1much-greater-than𝛾1\gamma\gg 1italic_γ ≫ 1) can be estimated as frad(2/3)re2γ2B~2subscript𝑓rad23subscriptsuperscript𝑟2esuperscript𝛾2superscriptsubscript~𝐵perpendicular-to2f_{\rm rad}\approx(2/3)r^{2}_{\rm e}\gamma^{2}\tilde{B}_{\perp}^{2}italic_f start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT ≈ ( 2 / 3 ) italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where resubscript𝑟er_{\rm e}italic_r start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT is the classical radius of the electron. The fiducial radiation-reaction-limited electron Lorentz factor is then given by

γrad=3eE2re2B~23×104(EB~)1/2(B~107G)1/2,subscript𝛾rad3𝑒subscript𝐸parallel-to2subscriptsuperscript𝑟2esubscriptsuperscript~𝐵2perpendicular-to3superscript104superscriptsubscript𝐸parallel-tosubscript~𝐵perpendicular-to12superscriptsubscript~𝐵perpendicular-tosuperscript107G12\gamma_{\rm rad}=\sqrt{\frac{3eE_{\parallel}}{2r^{2}_{\rm e}\tilde{B}^{2}_{% \perp}}}\approx 3\times 10^{4}\left(\frac{E_{\parallel}}{\tilde{B}_{\perp}}% \right)^{1/2}\left(\frac{\tilde{B}_{\perp}}{10^{7}\rm{G}}\right)^{-1/2},italic_γ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG 3 italic_e italic_E start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT over~ start_ARG italic_B end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT end_ARG end_ARG ≈ 3 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( divide start_ARG italic_E start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT end_ARG start_ARG over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G end_ARG ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT , (12)

or γradLC3.3×105subscriptsuperscript𝛾LCrad3.3superscript105\gamma^{\rm LC}_{\rm rad}\approx 3.3\times 10^{5}italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT ≈ 3.3 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT at the light cylinder for the fiducial pulsar parameters, assuming E=B~subscript𝐸parallel-tosubscript~𝐵perpendicular-toE_{\parallel}=\tilde{B}_{\perp}italic_E start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT = over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT and B~=BLCsubscript~𝐵perpendicular-tosubscript𝐵LC\tilde{B}_{\perp}=B_{\rm LC}over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT = italic_B start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT.

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 rsubscript𝑟r_{\star}italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT,

desr103δr2.5×1021<RLCr=5,similar-tosubscriptsuperscript𝑑sesubscript𝑟superscript103much-less-than𝛿subscript𝑟similar-to2.5superscript102much-less-than1subscript𝑅LCsubscript𝑟5\frac{d^{\rm s}_{\rm e}}{r_{\star}}\sim 10^{-3}\ll\frac{\delta}{r_{\star}}\sim 2% .5\times 10^{-2}\ll 1<\frac{R_{\rm LC}}{r_{\star}}=5,divide start_ARG italic_d start_POSTSUPERSCRIPT roman_s end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT end_ARG start_ARG italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG ∼ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT ≪ divide start_ARG italic_δ end_ARG start_ARG italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG ∼ 2.5 × 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT ≪ 1 < divide start_ARG italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG start_ARG italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG = 5 , (13)

and in terms of energy scales normalized to γpcsubscript𝛾pc\gamma_{\rm pc}italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT

γsγpc8.5×104<γradLCγpc103<γthγpc4×1031.similar-tosubscript𝛾ssubscript𝛾pc8.5superscript104subscriptsuperscript𝛾LCradsubscript𝛾pcsimilar-tosuperscript103subscript𝛾thsubscript𝛾pcsimilar-to4superscript103much-less-than1\frac{\gamma_{\rm s}}{\gamma_{\rm pc}}\sim 8.5\times 10^{-4}<\frac{\gamma^{\rm LC% }_{\rm rad}}{\gamma_{\rm pc}}\sim 10^{-3}<\frac{\gamma_{\rm th}}{\gamma_{\rm pc% }}\sim 4\times 10^{-3}\ll 1.divide start_ARG italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT end_ARG start_ARG italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT end_ARG ∼ 8.5 × 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT < divide start_ARG italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT end_ARG start_ARG italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT end_ARG ∼ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT < divide start_ARG italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT end_ARG start_ARG italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT end_ARG ∼ 4 × 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT ≪ 1 . (14)

This exercise reveals that in such system γssubscript𝛾s\gamma_{\rm s}italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT, γradsubscript𝛾rad\gamma_{\rm rad}italic_γ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT and γthsubscript𝛾th\gamma_{\rm th}italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT 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 L1033similar-to𝐿superscript1033L\sim 10^{33}italic_L ∼ 10 start_POSTSUPERSCRIPT 33 end_POSTSUPERSCRIPTerg/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., κ1much-greater-than𝜅1\kappa\gg 1italic_κ ≫ 1). 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 κ1much-greater-than𝜅1\kappa\gg 1italic_κ ≫ 1), 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

d𝐮dt=qmc(𝐄+𝐮×𝐁γ)+𝐠,d𝐮d𝑡𝑞𝑚𝑐𝐄𝐮𝐁𝛾𝐠\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\frac{q}{mc}\left(\mathbf{E}+\frac{% \mathbf{u}\times\mathbf{B}}{\gamma}\right)+\mathbf{g},divide start_ARG roman_d bold_u end_ARG start_ARG roman_d italic_t end_ARG = divide start_ARG italic_q end_ARG start_ARG italic_m italic_c end_ARG ( bold_E + divide start_ARG bold_u × bold_B end_ARG start_ARG italic_γ end_ARG ) + bold_g , (15)

where 𝐮=γ𝐯/c𝐮𝛾𝐯𝑐\mathbf{u}=\gamma\mathbf{v}/cbold_u = italic_γ bold_v / italic_c is the dimensionless particle 4-momentum vector, and q𝑞qitalic_q is the particle electric charge, and where

𝐠23re2[(𝐄+𝜷×𝐁)×𝐁+(𝜷𝐄)𝐄]23re2γ2B~2𝜷,𝐠23subscriptsuperscript𝑟2edelimited-[]𝐄𝜷𝐁𝐁𝜷𝐄𝐄23subscriptsuperscript𝑟2esuperscript𝛾2subscriptsuperscript~𝐵2perpendicular-to𝜷\mathbf{g}\approx\frac{2}{3}r^{2}_{\rm e}\left[\left(\mathbf{E}+\boldsymbol{% \beta}\times\mathbf{B}\right)\times\mathbf{B}+\left(\boldsymbol{\beta}\cdot% \mathbf{E}\right)\mathbf{E}\right]-\frac{2}{3}r^{2}_{\rm e}\gamma^{2}\tilde{B}% ^{2}_{\perp}\boldsymbol{\beta},bold_g ≈ divide start_ARG 2 end_ARG start_ARG 3 end_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT [ ( bold_E + bold_italic_β × bold_B ) × bold_B + ( bold_italic_β ⋅ bold_E ) bold_E ] - divide start_ARG 2 end_ARG start_ARG 3 end_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT over~ start_ARG italic_B end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT bold_italic_β , (16)

is the radiation reaction force following the Landau & Lifshitz (1971) formulation, and 𝜷=𝐯/c𝜷𝐯𝑐\boldsymbol{\beta}=\mathbf{v}/cbold_italic_β = bold_v / italic_c, and where

B~=(𝐄+𝜷×𝐁)2(𝜷𝐄)2,subscript~𝐵perpendicular-tosuperscript𝐄𝜷𝐁2superscript𝜷𝐄2\tilde{B}_{\perp}=\sqrt{\left(\mathbf{E}+\boldsymbol{\beta}\times\mathbf{B}% \right)^{2}-\left(\boldsymbol{\beta}\cdot\mathbf{E}\right)^{2}},over~ start_ARG italic_B end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT = square-root start_ARG ( bold_E + bold_italic_β × bold_B ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( bold_italic_β ⋅ bold_E ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (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 𝐉𝐉\mathbf{J}bold_J reconstructed on the grid is then used to evolve the time-dependent Maxwell equations,

𝐁t=c×𝐄𝐁𝑡𝑐bold-∇𝐄\frac{\partial\mathbf{B}}{\partial t}=-c\,\boldsymbol{\nabla}\times\mathbf{E}divide start_ARG ∂ bold_B end_ARG start_ARG ∂ italic_t end_ARG = - italic_c bold_∇ × bold_E (18)
𝐄t=c×𝐁4π𝐉,𝐄𝑡𝑐bold-∇𝐁4𝜋𝐉\frac{\partial\mathbf{E}}{\partial t}=c\,\boldsymbol{\nabla}\times\mathbf{B}-4% \pi\mathbf{J},divide start_ARG ∂ bold_E end_ARG start_ARG ∂ italic_t end_ARG = italic_c bold_∇ × bold_B - 4 italic_π bold_J , (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).

Refer to caption
Figure 1: This diagram presents the new hybrid numerical scheme proposed in this work during a time integration cycle. This method is meant to blend the force-free 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, 𝐉𝐉\mathbf{J}bold_J, using a blending function f𝑓fitalic_f which solely depends on the magnetic flux function ΨΨ\Psiroman_Ψ.

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

FμνIν=0,subscript𝐹𝜇𝜈superscript𝐼𝜈0F_{\mu\nu}I^{\nu}=0,italic_F start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT italic_I start_POSTSUPERSCRIPT italic_ν end_POSTSUPERSCRIPT = 0 , (20)

where Fμνsubscript𝐹𝜇𝜈F_{\mu\nu}italic_F start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT is the electromagnetic tensor and Iνsuperscript𝐼𝜈I^{\nu}italic_I start_POSTSUPERSCRIPT italic_ν end_POSTSUPERSCRIPT is the 4-current, which translates into the following conditions

ρ𝐄+𝐉×𝐁c=𝟎,𝜌𝐄𝐉𝐁𝑐0\rho\mathbf{E}+\frac{\mathbf{J}\times\mathbf{B}}{c}=\mathbf{0},italic_ρ bold_E + divide start_ARG bold_J × bold_B end_ARG start_ARG italic_c end_ARG = bold_0 , (21)

using the spatial components, where ρ=𝐄/4π𝜌bold-∇𝐄4𝜋\rho=\boldsymbol{\nabla}\cdot\mathbf{E}/4\piitalic_ρ = bold_∇ ⋅ bold_E / 4 italic_π is the charge density, and

𝐄𝐉=0,𝐄𝐉0\mathbf{E}\cdot\mathbf{J}=0,bold_E ⋅ bold_J = 0 , (22)

using the temporal component. Eq. (21) also leads to the condition

𝐄𝐁=0.𝐄𝐁0\mathbf{E}\cdot\mathbf{B}=0.bold_E ⋅ bold_B = 0 . (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)

𝐉=c4π𝐄(𝐄×𝐁B2)+c4π(𝐁(×𝐁)𝐄(×𝐄))𝐁B2,𝐉𝑐4𝜋bold-∇𝐄𝐄𝐁superscript𝐵2𝑐4𝜋𝐁bold-∇𝐁𝐄bold-∇𝐄𝐁superscript𝐵2\mathbf{J}=\frac{c}{4\pi}\boldsymbol{\nabla}\cdot\mathbf{E}\left(\frac{\mathbf% {E}\times\mathbf{B}}{B^{2}}\right)+\frac{c}{4\pi}\Bigl{(}\mathbf{B}\cdot\left(% \boldsymbol{\nabla}\times\mathbf{B}\right)-\mathbf{E}\cdot\left(\boldsymbol{% \nabla}\times\mathbf{E}\right)\Bigr{)}\frac{\mathbf{B}}{B^{2}},bold_J = divide start_ARG italic_c end_ARG start_ARG 4 italic_π end_ARG bold_∇ ⋅ bold_E ( divide start_ARG bold_E × bold_B end_ARG start_ARG italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) + divide start_ARG italic_c end_ARG start_ARG 4 italic_π end_ARG ( bold_B ⋅ ( bold_∇ × bold_B ) - bold_E ⋅ ( bold_∇ × bold_E ) ) divide start_ARG bold_B end_ARG start_ARG italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (24)

which can be understood as the sum of a current flowing across field lines, 𝐉subscript𝐉perpendicular-to\mathbf{J_{\perp}}bold_J start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT, and a component flowing along field lines, 𝐉subscript𝐉parallel-to\mathbf{J_{\parallel}}bold_J start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT. 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, 𝐉subscript𝐉perpendicular-to\mathbf{J}_{\perp}bold_J start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT. The parallel component, Jsubscript𝐽parallel-toJ_{\parallel}italic_J start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT, prevents the formation of an electric field component along the magnetic field at all times, 𝐄𝐁=0𝐄𝐁0\mathbf{E}\cdot\mathbf{B}=0bold_E ⋅ bold_B = 0, but computing Jsubscript𝐽parallel-toJ_{\parallel}italic_J start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT 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,

𝐄𝐄(𝐄𝐁)𝐁B2.𝐄𝐄𝐄𝐁𝐁superscript𝐵2\mathbf{E}\leftarrow\mathbf{E}-\left(\mathbf{E}\cdot\mathbf{B}\right)\frac{% \mathbf{B}}{B^{2}}.bold_E ← bold_E - ( bold_E ⋅ bold_B ) divide start_ARG bold_B end_ARG start_ARG italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (25)

The code also checks at every time step that the condition E2<B2superscript𝐸2superscript𝐵2E^{2}<B^{2}italic_E start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT < italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT is always satisfied, so that if it is not the case, the electric field is once again renormalized as

𝐄B2E2𝐄.𝐄superscript𝐵2superscript𝐸2𝐄\mathbf{E}\leftarrow\sqrt{\frac{B^{2}}{E^{2}}}\mathbf{E}.bold_E ← square-root start_ARG divide start_ARG italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_E start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG bold_E . (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 (r,θ𝑟𝜃r,\thetaitalic_r , italic_θ), the magnetic flux function is given by

Ψ=12π𝐁𝑑𝐒=0θB𝐫r2sinθdθ.Ψ12𝜋double-integral𝐁differential-d𝐒superscriptsubscript0𝜃subscript𝐵𝐫superscript𝑟2superscript𝜃𝑑superscript𝜃\Psi=\frac{1}{2\pi}\iint\mathbf{B}\cdot d\mathbf{S}=\int_{0}^{\theta}B_{% \mathbf{r}}r^{2}\sin\theta^{\prime}d\theta^{\prime}.roman_Ψ = divide start_ARG 1 end_ARG start_ARG 2 italic_π end_ARG ∬ bold_B ⋅ italic_d bold_S = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT italic_B start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_sin italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_d italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT . (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

𝐉=(1f(Ψ))𝐉PIC+f(Ψ)𝐉FFE,𝐉1𝑓Ψsubscript𝐉PIC𝑓Ψsubscript𝐉FFE\mathbf{J}=\Bigl{(}1-f\left(\Psi\right)\Bigr{)}\,\mathbf{J}_{\rm PIC}+f\left(% \Psi\right)\,\mathbf{J}_{\rm FFE},bold_J = ( 1 - italic_f ( roman_Ψ ) ) bold_J start_POSTSUBSCRIPT roman_PIC end_POSTSUBSCRIPT + italic_f ( roman_Ψ ) bold_J start_POSTSUBSCRIPT roman_FFE end_POSTSUBSCRIPT , (28)

where 𝐉PICsubscript𝐉PIC\mathbf{J}_{\rm PIC}bold_J start_POSTSUBSCRIPT roman_PIC end_POSTSUBSCRIPT is the current density reconstructed from the particles, 𝐉FFEsubscript𝐉FFE\mathbf{J}_{\rm FFE}bold_J start_POSTSUBSCRIPT roman_FFE end_POSTSUBSCRIPT is the current from the force-free condition (Eq. 24), and f𝑓fitalic_f 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, Ψ0<Ψ1subscriptΨ0subscriptΨ1\Psi_{0}<\Psi_{1}roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, so that the blending function is

f(Ψ)𝑓Ψ\displaystyle f\left(\Psi\right)italic_f ( roman_Ψ ) =\displaystyle== 1,ifΨ<Ψ0,pure force-free, no particles,\displaystyle 1,\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\ \Psi<\Psi_{% 0},\rightarrow\text{pure force-free, no particles},1 , if roman_Ψ < roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , → pure force-free, no particles , (29)
f(Ψ)𝑓Ψ\displaystyle f\left(\Psi\right)italic_f ( roman_Ψ ) =\displaystyle== Ψ1ΨΨ1Ψ0,ifΨ0<Ψ<Ψ1,transition layer,\displaystyle\frac{\Psi_{1}-\Psi}{\Psi_{1}-\Psi_{0}},\leavevmode\nobreak\ % \text{if}\leavevmode\nobreak\ \Psi_{0}<\Psi<\Psi_{1},\rightarrow\text{% transition layer},divide start_ARG roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - roman_Ψ end_ARG start_ARG roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG , if roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < roman_Ψ < roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , → transition layer , (30)
f(Ψ)𝑓Ψ\displaystyle f\left(\Psi\right)italic_f ( roman_Ψ ) =\displaystyle== 0,ifΨ>Ψ1,pure PIC.\displaystyle 0,\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\ \Psi>\Psi_{% 1},\rightarrow\text{pure PIC}.0 , if roman_Ψ > roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , → pure PIC . (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

ΨM=Ψ(1cosθ),subscriptΨMsubscriptΨ1𝜃\Psi_{\rm M}=\Psi_{\star}\left(1-\cos\theta\right),roman_Ψ start_POSTSUBSCRIPT roman_M end_POSTSUBSCRIPT = roman_Ψ start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ( 1 - roman_cos italic_θ ) , (32)

where Ψ=Br2subscriptΨsubscript𝐵subscriptsuperscript𝑟2\Psi_{\star}=B_{\star}r^{2}_{\star}roman_Ψ start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT. An exact solution to the Grad-Shafranov equation yields the toroidal magnetic field component, (Michel, 1973)

Bϕ=ΨRLCsinθr,subscript𝐵italic-ϕsubscriptΨsubscript𝑅LC𝜃𝑟B_{\phi}=-\frac{\Psi_{\star}}{R_{\rm LC}}\frac{\sin\theta}{r},italic_B start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = - divide start_ARG roman_Ψ start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT end_ARG divide start_ARG roman_sin italic_θ end_ARG start_ARG italic_r end_ARG , (33)

if 𝛀𝐁>0𝛀𝐁0\boldsymbol{\Omega}\cdot\mathbf{B}>0bold_Ω ⋅ bold_B > 0 in the northern hemisphere (and 𝛀𝐁<0𝛀𝐁0\boldsymbol{\Omega}\cdot\mathbf{B}<0bold_Ω ⋅ bold_B < 0 in the southern hemisphere). In this solution, the electric field is purely along θ𝜃\thetaitalic_θ and matches the toroidal magnetic field, Eθ=Bϕsubscript𝐸𝜃subscript𝐵italic-ϕE_{\theta}=B_{\phi}italic_E start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT = italic_B start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT. The electric current is purely radial and is given by

𝐉M=c4π×𝐁=JGJ(rr)2cosθ𝐞r,subscript𝐉M𝑐4𝜋bold-∇𝐁subscriptsuperscript𝐽GJsuperscriptsubscript𝑟𝑟2𝜃subscript𝐞r\mathbf{J_{\rm M}}=\frac{c}{4\pi}\boldsymbol{\nabla}\times\mathbf{B}=-J^{\star% }_{\rm GJ}\left(\frac{r_{\star}}{r}\right)^{2}\cos\theta\leavevmode\nobreak\ % \mathbf{e_{\rm r}},bold_J start_POSTSUBSCRIPT roman_M end_POSTSUBSCRIPT = divide start_ARG italic_c end_ARG start_ARG 4 italic_π end_ARG bold_∇ × bold_B = - italic_J start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT ( divide start_ARG italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG start_ARG italic_r end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_cos italic_θ bold_e start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT , (34)

where JGJ=ΩB/2πsubscriptsuperscript𝐽GJΩsubscript𝐵2𝜋J^{\star}_{\rm GJ}=\Omega B_{\star}/2\piitalic_J start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT = roman_Ω italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / 2 italic_π corresponds to the fiducial Goldreich-Julian current, and the spindown power associated to the electromagnetic torque is

LM=c4π(𝐄×𝐁)𝑑𝐒=c211r2Bϕ2dcosθ=2Ψ2Ω23c.subscript𝐿M𝑐4𝜋double-integral𝐄𝐁differential-d𝐒𝑐2superscriptsubscript11superscript𝑟2subscriptsuperscript𝐵2italic-ϕ𝑑𝜃2subscriptsuperscriptΨ2superscriptΩ23𝑐L_{\rm M}=\frac{c}{4\pi}\iint\left(\mathbf{E}\times\mathbf{B}\right)\cdot d% \mathbf{S}=\frac{c}{2}\int_{-1}^{1}r^{2}B^{2}_{\phi}d\cos\theta=\frac{2\Psi^{2% }_{\star}\Omega^{2}}{3c}.italic_L start_POSTSUBSCRIPT roman_M end_POSTSUBSCRIPT = divide start_ARG italic_c end_ARG start_ARG 4 italic_π end_ARG ∬ ( bold_E × bold_B ) ⋅ italic_d bold_S = divide start_ARG italic_c end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT - 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_d roman_cos italic_θ = divide start_ARG 2 roman_Ψ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 3 italic_c end_ARG . (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 10242superscript102421024^{2}1024 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, 20482superscript204822048^{2}2048 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT and 40962superscript409624096^{2}4096 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT cells. The grid spacing is logarithmic along the radial direction and constant along the θ𝜃\thetaitalic_θ direction. The domain size ranges from the star surface, rmin=rsubscript𝑟minsubscript𝑟r_{\rm min}=r_{\star}italic_r start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT = italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT, up to rmax=(10/3)RLCsubscript𝑟max103subscript𝑅LCr_{\rm max}=(10/3)R_{\rm LC}italic_r start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = ( 10 / 3 ) italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT, where RLC=5rsubscript𝑅LC5subscript𝑟R_{\rm LC}=5r_{\star}italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT = 5 italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT, and covers the full range in θ[0,π]𝜃0𝜋\theta\in[0,\pi]italic_θ ∈ [ 0 , italic_π ]. The magnetosphere is initialized in vacuum with a pure radial magnetic field, B=B(r/r)2𝐵subscript𝐵superscriptsubscript𝑟𝑟2B=B_{\star}(r_{\star}/r)^{2}italic_B = italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / italic_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. The surface magnetic field strength is fixed at B=5×105subscript𝐵5superscript105B_{\star}=5\times 10^{5}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 5 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPTG, 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 𝐄=(𝛀×𝐫)×𝐁/csubscript𝐄𝛀subscript𝐫subscript𝐁𝑐\mathbf{E_{\star}}=-(\boldsymbol{\Omega}\times\mathbf{r_{\star}})\times\mathbf% {B_{\star}}/cbold_E start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = - ( bold_Ω × bold_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ) × bold_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT / italic_c. This procedure sets the magnetosphere into rotation. An absorbing layer starting at rabs=3RLCsubscript𝑟abs3subscript𝑅LCr_{\rm abs}=3R_{\rm LC}italic_r start_POSTSUBSCRIPT roman_abs end_POSTSUBSCRIPT = 3 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT 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 θ=π/2𝜃𝜋2\theta=\pi/2italic_θ = italic_π / 2, or Ψ=ΨΨsubscriptΨ\Psi=\Psi_{\star}roman_Ψ = roman_Ψ start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT. In the PIC domain, electron-positron pairs are continuously injected with a high multiplicity, κinj=10subscript𝜅inj10\kappa_{\rm inj}=10italic_κ start_POSTSUBSCRIPT roman_inj end_POSTSUBSCRIPT = 10, 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.

Refer to caption
Refer to caption
Refer to caption
Figure 2: Simulated monopolar magnetosphere using the hybrid force-free-PIC scheme. The division between the force-free and the 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, JGJsubscriptsuperscript𝐽GJJ^{\star}_{\rm GJ}italic_J start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT, and compensated by (r/r)2superscript𝑟subscript𝑟2(r/r_{\star})^{2}( italic_r / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for visualization purposes, for the 40962superscript409624096^{2}4096 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT cells run. Solid lines show the poloidal magnetic field lines. Middle: Cross-section of the current density profile at r=1.5RLC𝑟1.5subscript𝑅LCr=1.5R_{\rm LC}italic_r = 1.5 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT as a function of the numerical resolution, and comparison with the exact profile (red dashed-dotted line, Eq. 34). Bottom: Radial profiles of the Poynting flux, L(r)𝐿𝑟L(r)italic_L ( italic_r ), normalized to the monopole solution, LMsubscript𝐿ML_{\rm M}italic_L start_POSTSUBSCRIPT roman_M end_POSTSUBSCRIPT (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.

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

ϵ=|LinLout|Loutitalic-ϵsubscript𝐿insubscript𝐿outsubscript𝐿out\epsilon=\frac{\left|L_{\rm in}-L_{\rm out}\right|}{L_{\rm out}}italic_ϵ = divide start_ARG | italic_L start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT - italic_L start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT | end_ARG start_ARG italic_L start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT end_ARG (36)

where Linsubscript𝐿inL_{\rm in}italic_L start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT is the Poynting flux averaged in the r[r,2r]𝑟subscript𝑟2subscript𝑟r\in[r_{\star},2r_{\star}]italic_r ∈ [ italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT , 2 italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ] interval, and Loutsubscript𝐿outL_{\rm out}italic_L start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT is averaged over the r[2RLC,3RLC]𝑟2subscript𝑅LC3subscript𝑅LCr\in[2R_{\rm LC},3R_{\rm LC}]italic_r ∈ [ 2 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT , 3 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ] interval (gray intervals in Fig. 2, bottom panel). We observe numerical convergence and a deviation smaller than ϵ<1%italic-ϵpercent1\epsilon<1\%italic_ϵ < 1 % for the 10242superscript102421024^{2}1024 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT run, down to ϵ0.15%similar-toitalic-ϵpercent0.15\epsilon\sim 0.15\%italic_ϵ ∼ 0.15 % 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

Physical parameters Values
Neutron star radius r=10kmsubscript𝑟10kmr_{\star}=10\;\mathrm{km}italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 roman_km
Spin period P=1ms𝑃1msP=1\,\mathrm{ms}italic_P = 1 roman_ms
Surface magnetic field B=107Gsubscript𝐵superscript107GB_{\star}=10^{7}\;\mathrm{G}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_G
Spindown power L0=4.8×1033erg/ssubscript𝐿04.8superscript1033ergsL_{0}=4.8\times 10^{33}\rm{erg/s}italic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 4.8 × 10 start_POSTSUPERSCRIPT 33 end_POSTSUPERSCRIPT roman_erg / roman_s
Polar cap energy γpc=2.6×108subscript𝛾pc2.6superscript108\gamma_{\rm pc}=2.6\times 10^{8}italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = 2.6 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT
Threshold energy γth=106subscript𝛾thsuperscript106\gamma_{\rm th}=10^{6}italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT
Secondary pairs energy γs=2.2×105subscript𝛾s2.2superscript105\gamma_{\rm s}=2.2\times 10^{5}italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 2.2 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT
Radiation-reaction energy γradLC=3.3×105subscriptsuperscript𝛾LCrad3.3superscript105\gamma^{\rm LC}_{\rm rad}=3.3\times 10^{5}italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = 3.3 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT
Electron plasma skindepth des=103cmsuperscriptsubscript𝑑essuperscript103cmd_{\rm e}^{\rm s}=10^{3}\mathrm{cm}italic_d start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_s end_POSTSUPERSCRIPT = 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT roman_cm
Plasma timescale ωpe1=3.3×108ssubscriptsuperscript𝜔1pe3.3superscript108s\omega^{-1}_{\rm pe}=3.3\times 10^{-8}\text{s}italic_ω start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_pe end_POSTSUBSCRIPT = 3.3 × 10 start_POSTSUPERSCRIPT - 8 end_POSTSUPERSCRIPT s
Synchrotron cooling time tsynLC=3.6×107ssubscriptsuperscript𝑡LCsyn3.6superscript107st^{\rm LC}_{\rm syn}=3.6\times 10^{-7}\text{s}italic_t start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_syn end_POSTSUBSCRIPT = 3.6 × 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT s
Numerical parameters Values
# grid cells (r,θ𝑟𝜃r,\thetaitalic_r , italic_θ) 8192×8192819281928192\times 81928192 × 8192
Proton to electron mass ratio mp/me=1836subscript𝑚psubscript𝑚e1836m_{\rm p}/m_{\rm e}=1836italic_m start_POSTSUBSCRIPT roman_p end_POSTSUBSCRIPT / italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT = 1836
Time step Δt=4.3×109Δ𝑡4.3superscript109\Delta t=4.3\times 10^{-9}roman_Δ italic_t = 4.3 × 10 start_POSTSUPERSCRIPT - 9 end_POSTSUPERSCRIPTs
Highest spatial resolution Δr=3.4×102Δ𝑟3.4superscript102\Delta r=3.4\times 10^{2}roman_Δ italic_r = 3.4 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPTcm
FFE domain I boundaries Ψmin=0,Ψ1=0.9Ψpcformulae-sequencesubscriptΨmin0subscriptΨ10.9subscriptΨpc\Psi_{\rm min}=0,\leavevmode\nobreak\ \Psi_{1}=0.9\Psi_{\rm pc}roman_Ψ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT = 0 , roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.9 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT
PIC domain boundaries Ψ0=0.85Ψpc,Ψ3=2.4Ψpcformulae-sequencesubscriptΨ00.85subscriptΨpcsubscriptΨ32.4subscriptΨpc\Psi_{0}=0.85\Psi_{\rm pc},\leavevmode\nobreak\ \Psi_{3}=2.4\Psi_{\rm pc}roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0.85 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT , roman_Ψ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = 2.4 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT
FFE domain II boundaries Ψ2=2.3Ψpc,Ψmax=μ/rformulae-sequencesubscriptΨ22.3subscriptΨpcsubscriptΨmax𝜇subscript𝑟\Psi_{2}=2.3\Psi_{\rm pc},\leavevmode\nobreak\ \Psi_{\rm max}=\mu/r_{\star}roman_Ψ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 2.3 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT , roman_Ψ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = italic_μ / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT
Table 1: Physical and numerical parameters for the reference weak gamma-ray pulsar simulation using the hybrid force-free-PIC scheme. Energy scales refer to electron Lorentz factors.
Refer to caption
Figure 3: Global hybrid force-free-PIC model of the aligned magnetosphere of a B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG surface magnetic field neutron star with a P=1𝑃1P=1italic_P = 1ms spin period at simulated time t=2.45𝑡2.45t=2.45italic_t = 2.45ms. Black solid lines show the poloidal magnetic field lines (i.e., isocontour of the magnetic flux function ΨΨ\Psiroman_Ψ), green and red solid lines show the transition from the force-free regions to the PIC region (i.e., Ψ0subscriptΨ0\Psi_{0}roman_Ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, Ψ1subscriptΨ1\Psi_{1}roman_Ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, Ψ2subscriptΨ2\Psi_{2}roman_Ψ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, Ψ3subscriptΨ3\Psi_{3}roman_Ψ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT, see Table 1). Panel a): Map of the pairs density normalized by nGJsubscriptsuperscript𝑛GJn^{\star}_{\rm GJ}italic_n start_POSTSUPERSCRIPT ⋆ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_GJ end_POSTSUBSCRIPT and compensated by (r/r)2superscript𝑟subscript𝑟2(r/r_{\star})^{2}( italic_r / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Panel b): Map of the mean pair energy, =γmec2𝛾subscript𝑚esuperscript𝑐2\mathcal{E}=\gamma m_{\rm e}c^{2}caligraphic_E = italic_γ italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Panels c) and d): same as panels a) and b) for the protons. The panels do not show the full domain but instead are focused on the PIC domain where the current sheet forms. A zoomed-in view on a reconnection site is shown in panel a). Notice how thin the reconnection layer is.

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 81922superscript819228192^{2}8192 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 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

ΨD=μsin2θr.subscriptΨD𝜇superscript2𝜃𝑟\Psi_{\rm D}=\frac{\mu\sin^{2}\theta}{r}.roman_Ψ start_POSTSUBSCRIPT roman_D end_POSTSUBSCRIPT = divide start_ARG italic_μ roman_sin start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_θ end_ARG start_ARG italic_r end_ARG . (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, mp=1836mesubscript𝑚p1836subscript𝑚em_{\rm p}=1836m_{\rm e}italic_m start_POSTSUBSCRIPT roman_p end_POSTSUBSCRIPT = 1836 italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT. 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 γthsubscript𝛾th\gamma_{\rm th}italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT, then a new pair is created with a Lorentz factor γssubscript𝛾s\gamma_{\rm s}italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT 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, 0<Ψ<0.9Ψpc0Ψ0.9subscriptΨpc0<\Psi<0.9\Psi_{\rm pc}0 < roman_Ψ < 0.9 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT. This region corresponds to the open field lines connected to the star polar caps.

  • PIC domain, 0.85Ψpc<Ψ<2.4Ψpc0.85subscriptΨpcΨ2.4subscriptΨpc0.85\Psi_{\rm pc}<\Psi<2.4\Psi_{\rm pc}0.85 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT < roman_Ψ < 2.4 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT. The domain encompasses the equatorial and separatrix current layers.

  • Force-free domain II, 2.3Ψpc<Ψ<Ψmax2.3subscriptΨpcΨsubscriptΨmax2.3\Psi_{\rm pc}<\Psi<\Psi_{\rm max}2.3 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT < roman_Ψ < roman_Ψ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. It is inside the closed field line regions deep inside the light cylinder.

In the above, Ψmax=μ/rsubscriptΨmax𝜇subscript𝑟\Psi_{\rm max}=\mu/r_{\star}roman_Ψ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = italic_μ / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT is the maximum of the flux function that is located at the star surface in the equatorial plane, and Ψpc=μ/RLCsubscriptΨpc𝜇subscript𝑅LC\Psi_{\rm pc}=\mu/R_{\rm LC}roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = italic_μ / italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT 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 ΔΨ=0.1ΨpcΔΨ0.1subscriptΨpc\Delta\Psi=0.1\Psi_{\rm pc}roman_Δ roman_Ψ = 0.1 roman_Ψ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT-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 r=106subscript𝑟superscript106r_{\star}=10^{6}italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPTcm, and the light-cylinder radius to RLC=5×106subscript𝑅LC5superscript106R_{\rm LC}=5\times 10^{6}italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT = 5 × 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPTcm corresponding to a P=1𝑃1P=1italic_P = 1ms spin period. The surface magnetic field of the reference production run is fixed at B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG, giving BLC=μ/RLC38×104subscript𝐵LC𝜇subscriptsuperscript𝑅3LC8superscript104B_{\rm LC}=\mu/R^{3}_{\rm LC}\approx 8\times 10^{4}italic_B start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT = italic_μ / italic_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ≈ 8 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPTG at the light cylinder, and a force-free spindown power L0=μ2Ω4/c3=4.8×1033subscript𝐿0superscript𝜇2superscriptΩ4superscript𝑐34.8superscript1033L_{0}=\mu^{2}\Omega^{4}/c^{3}=4.8\times 10^{33}italic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_μ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT / italic_c start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT = 4.8 × 10 start_POSTSUPERSCRIPT 33 end_POSTSUPERSCRIPTerg/s. The magnetic field strength then sets the relevant energy scales given in Sect. 2 as, γpc=2.6×108subscript𝛾pc2.6superscript108\gamma_{\rm pc}=2.6\times 10^{8}italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT = 2.6 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT, γth=106subscript𝛾thsuperscript106\gamma_{\rm th}=10^{6}italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT, γs=2.2×105subscript𝛾s2.2superscript105\gamma_{\rm s}=2.2\times 10^{5}italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 2.2 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT and γrad=3×104subscript𝛾rad3superscript104\gamma_{\rm rad}=3\times 10^{4}italic_γ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = 3 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT. 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 des/Δr3subscriptsuperscript𝑑seΔ𝑟3d^{\rm s}_{\rm e}/\Delta r\approx 3italic_d start_POSTSUPERSCRIPT roman_s end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT / roman_Δ italic_r ≈ 3 and ωpe1/Δt8subscriptsuperscript𝜔1peΔ𝑡8\omega^{-1}_{\rm pe}/\Delta t\approx 8italic_ω start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_pe end_POSTSUBSCRIPT / roman_Δ italic_t ≈ 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, ρssubscript𝜌s\rho_{\rm s}italic_ρ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT, divided by the local grid spacing, RLCΔr/rsubscript𝑅LCΔ𝑟subscript𝑟R_{\rm LC}\Delta r/r_{\star}italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT roman_Δ italic_r / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT, is ρs/(RLCΔr/r)=(γsmec2/eBLC)/(RLCΔr/r)3subscript𝜌ssubscript𝑅LCΔ𝑟subscript𝑟subscript𝛾ssubscript𝑚esuperscript𝑐2𝑒subscript𝐵LCsubscript𝑅LCΔ𝑟subscript𝑟3\rho_{\rm s}/(R_{\rm LC}\Delta r/r_{\star})=(\gamma_{\rm s}m_{\rm e}c^{2}/eB_{% \rm LC})/(R_{\rm LC}\Delta r/r_{\star})\approx 3italic_ρ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT / ( italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT roman_Δ italic_r / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ) = ( italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_e italic_B start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ) / ( italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT roman_Δ italic_r / italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT ) ≈ 3, and Larmor timescale ρs/(cΔt)36subscript𝜌s𝑐Δ𝑡36\rho_{\rm s}/(c\Delta t)\approx 36italic_ρ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT / ( italic_c roman_Δ italic_t ) ≈ 36. The corresponding synchrotron cooling time is of the order of tsynLC3mec/(2re2γsBLC2)85Δtsimilar-tosubscriptsuperscript𝑡LCsyn3subscript𝑚e𝑐2subscriptsuperscript𝑟2esubscript𝛾ssubscriptsuperscript𝐵2LC85Δ𝑡t^{\rm LC}_{\rm syn}\sim 3m_{\rm e}c/(2r^{2}_{\rm e}\gamma_{\rm s}B^{2}_{\rm LC% })\approx 85\Delta titalic_t start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_syn end_POSTSUBSCRIPT ∼ 3 italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c / ( 2 italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT italic_B start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT ) ≈ 85 roman_Δ italic_t.

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 rsubscript𝑟r_{\star}italic_r start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT-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 1%percent11\%1 % of its nominal value, and it is then slowly increasing in time up to about 10%percent1010\%10 %. In a last step, it is gradually increased up to the nominal strength at t2Pgreater-than-or-equivalent-to𝑡2𝑃t\gtrsim 2Pitalic_t ≳ 2 italic_P 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 B=104subscript𝐵superscript104B_{\star}=10^{4}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPTG, 105superscript10510^{5}10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPTG, and 106superscript10610^{6}10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPTG, but keep the ratio between all energy scales the same as the B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG simulation, including the radiation-reaction-limited lepton Lorentz factor γradsubscript𝛾rad\gamma_{\rm rad}italic_γ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT, given in Eq. (14). To keep the γpc/γradsubscript𝛾pcsubscript𝛾rad\gamma_{\rm pc}/\gamma_{\rm rad}italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT / italic_γ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT ratio identical to the B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG run, the strength of the radiation reaction force –or equivalently the Thomson cross section– must be boosted by a factor frad=109subscript𝑓radsuperscript109f_{\rm rad}=10^{9}italic_f start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 9 end_POSTSUPERSCRIPT for B=104subscript𝐵superscript104B_{\star}=10^{4}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPTG, frad=106subscript𝑓radsuperscript106f_{\rm rad}=10^{6}italic_f start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT for B=105subscript𝐵superscript105B_{\star}=10^{5}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT, and frad=103subscript𝑓radsuperscript103f_{\rm rad}=10^{3}italic_f start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT for B=106subscript𝐵superscript106B_{\star}=10^{6}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPTG. 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 t2.5P𝑡2.5𝑃t\approx 2.5Pitalic_t ≈ 2.5 italic_P for the fiducial full-scale B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG 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 0.6RLC0.6subscript𝑅LC0.6R_{\rm LC}0.6 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT), and slowly migrates towards RLCsubscript𝑅LCR_{\rm LC}italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT, 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 δ103RLCsimilar-to𝛿superscript103subscript𝑅LC\delta\sim 10^{-3}R_{\rm LC}italic_δ ∼ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT 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 κ10100similar-to𝜅10100\kappa\sim 10-100italic_κ ∼ 10 - 100, 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 0.1RLCsimilar-toabsent0.1subscript𝑅LC\sim 0.1R_{\rm LC}∼ 0.1 italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT 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.

Refer to caption
Figure 4: Total electron (solid blue), positron (solid red), and proton (dashed green) energy spectra for a B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG and P=1𝑃1P=1italic_P = 1ms pulsar. The relevant energy scales introduced in Sect. 2 are shown by vertical dotted lines, from the lowest to the highest energies: the radiation-reaction-limited energy at the light cylinder, γradLCmec2subscriptsuperscript𝛾LCradsubscript𝑚esuperscript𝑐2\gamma^{\rm LC}_{\rm rad}m_{\rm e}c^{2}italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, the threshold energy for pair production, γthmec2subscript𝛾thsubscript𝑚esuperscript𝑐2\gamma_{\rm th}m_{\rm e}c^{2}italic_γ start_POSTSUBSCRIPT roman_th end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, and the polar-cap potential drop, eΦpc𝑒subscriptΦpce\Phi_{\rm pc}italic_e roman_Φ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT.
Refer to caption
Figure 5: Maximum cut-off energy of the charged particle energy distributions, maxsubscriptmax\mathcal{E}_{\rm max}caligraphic_E start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT (electrons: blue dots, positrons: red dots, protons: green dots), as a function of the surface magnetic field Bsubscript𝐵B_{\star}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT. The evolution of the polar-cap potential drop is shown by the dashed line for comparison (Eq. 4).

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, pairs=γmec21subscriptpairs𝛾subscript𝑚esuperscript𝑐21\mathcal{E}_{\rm pairs}=\gamma m_{\rm e}c^{2}\approx 1caligraphic_E start_POSTSUBSCRIPT roman_pairs end_POSTSUBSCRIPT = italic_γ italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≈ 1TeV, and sharply decreases inside the nearest plasmoids on either side, down to pairs100less-than-or-similar-tosubscriptpairs100\mathcal{E}_{\rm pairs}\lesssim 100caligraphic_E start_POSTSUBSCRIPT roman_pairs end_POSTSUBSCRIPT ≲ 100GeV. 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, γradLCmec2150subscriptsuperscript𝛾LCradsubscript𝑚esuperscript𝑐2150\gamma^{\rm LC}_{\rm rad}m_{\rm e}c^{2}\approx 150italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≈ 150GeVpairsmuch-less-thanabsentsubscriptpairs\ll\mathcal{E}_{\rm pairs}≪ caligraphic_E start_POSTSUBSCRIPT roman_pairs end_POSTSUBSCRIPT; 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, 1010\mathcal{E}\approx 10caligraphic_E ≈ 10TeV, 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 maxβreceΦpcβrec×130similar-tosubscriptmaxsubscript𝛽rec𝑒subscriptΦpcsubscript𝛽rec130\mathcal{E}_{\rm max}\sim\beta_{\rm rec}e\Phi_{\rm pc}\approx\beta_{\rm rec}% \times 130caligraphic_E start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ∼ italic_β start_POSTSUBSCRIPT roman_rec end_POSTSUBSCRIPT italic_e roman_Φ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT ≈ italic_β start_POSTSUBSCRIPT roman_rec end_POSTSUBSCRIPT × 130TeV, where βrec0.1similar-tosubscript𝛽rec0.1\beta_{\rm rec}\sim 0.1italic_β start_POSTSUBSCRIPT roman_rec end_POSTSUBSCRIPT ∼ 0.1-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 (γradLCmec2150subscriptsuperscript𝛾LCradsubscript𝑚esuperscript𝑐2150\gamma^{\rm LC}_{\rm rad}m_{\rm e}c^{2}\approx 150italic_γ start_POSTSUPERSCRIPT roman_LC end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≈ 150GeV). 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 σLCmec2500subscript𝜎LCsubscript𝑚esuperscript𝑐2500\sigma_{\rm LC}m_{\rm e}c^{2}\approx 500italic_σ start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≈ 500GeV (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 𝛀𝝁>0𝛀𝝁0\boldsymbol{\Omega}\cdot\boldsymbol{\mu}>0bold_Ω ⋅ bold_italic_μ > 0. 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 1.5similar-toabsent1.5\sim-1.5∼ - 1.5, 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 γpcsubscript𝛾pc\gamma_{\rm pc}italic_γ start_POSTSUBSCRIPT roman_pc end_POSTSUBSCRIPT, 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.

Refer to caption
Figure 6: Total synchrotron spectral energy distribution emitted by the equatorial current sheet for B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG and P=1𝑃1P=1italic_P = 1ms (black dashed line). The electronic (positronic) contribution is shown by the blue (red) line. The vertical dotted line indicates the 100MeV energy threshold of the Fermi-LAT for comparison.
Refer to caption
Figure 7: Pulsar gamma-ray efficiency above 100MeV, that is the gamma-ray luminosity over the pulsar spindown power (Lγ/Lsubscript𝐿𝛾𝐿L_{\gamma}/Litalic_L start_POSTSUBSCRIPT italic_γ end_POSTSUBSCRIPT / italic_L), as reported in the third Fermi-LAT catalogue (gray stars, Smith et al. 2023), as a function of the pulsar spindown power (L𝐿Litalic_L). The properties of the pulsar simulated in this work (B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG, P=1𝑃1P=1italic_P = 1ms) 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.5P0.5𝑃0.5P0.5 italic_P time frame.

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 (RLCgreater-than-or-equivalent-toabsentsubscript𝑅LC\gtrsim R_{\rm LC}≳ italic_R start_POSTSUBSCRIPT roman_LC end_POSTSUBSCRIPT). 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 ϵrad(9/4)mec2/α160similar-tosubscriptitalic-ϵrad94subscript𝑚esuperscript𝑐2𝛼160\epsilon_{\rm rad}\sim(9/4)m_{\rm e}c^{2}/\alpha\approx 160italic_ϵ start_POSTSUBSCRIPT roman_rad end_POSTSUBSCRIPT ∼ ( 9 / 4 ) italic_m start_POSTSUBSCRIPT roman_e end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_α ≈ 160MeV (Uzdensky et al., 2011), where α𝛼\alphaitalic_α is the fine structure constant. Integrating over all frequencies, the synchrotron luminosity accounts for about 18%percent1818\%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 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 (0.1%less-than-or-similar-toabsentpercent0.1\lesssim 0.1\%≲ 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 non-thermal radiation, which is possible only in the strong cooling regime. Indeed, about 2%percent22\%2 % of the spindown is channelled into the beam of protons while the pairs leave the box with 1%less-than-or-similar-toabsentpercent1\lesssim 1\%≲ 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 Fermi-LAT, our fiducial pulsar shines about 4%percent44\%4 % of its total spindown power above  100100\>100100MeV 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 B=107subscript𝐵superscript107B_{\star}=10^{7}italic_B start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPTG simulated in this work, pairs are accelerated up to 1TeV, while protons reach up to 10%greater-than-or-equivalent-toabsentpercent10\gtrsim 10\%≳ 10 % of the full polar cap potential drop, that is 10greater-than-or-equivalent-toabsent10\gtrsim 10≳ 10TeV. 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 20%percent2020\%20 % of the spindown power, including about 4%percent44\%4 % 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