Open Access
Issue
A&A
Volume 690, October 2024
Article Number A170
Number of page(s) 12
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202450238
Published online 08 October 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Magnetospheres forming around neutron stars and black holes, simply referred to as relativistic magnetospheres in the following, are involved in some of the most extreme high-energy astrophysical phenomena such as pulsars (Cerutti & Beloborodov 2017; Philippov & Kramer 2022), magnetar flares and outbursts (Thompson & Duncan 1995; Kaspi & Beloborodov 2017), relativistic jets from supermassive black holes (Blandford & Znajek 1977; Blandford et al. 2019), possibly fast radio bursts (Popov & Postnov 2010; Lyubarsky 2020), and precursor emission to compact object mergers (Crinquand et al. 2019; Most & Philippov 2020). Elaborating a quantitative and predictive model of relativistic magnetospheres turned out to be an outstanding endeavor at the forefront of modern physics, and our current understanding of such systems is at best incomplete. The fact that these regions concentrate large energy densities into various forms (i.e., gravitational, rotational, electromagnetic) and that the interplay between plasma physics, general relativity, and relativistic (quantum) electrodynamics dramatically increases the degree of complexity of the problem, which becomes analytically untractable.

In the quest for a comprehensive theory of relativistic magnetospheres, the aligned rotator has historically played the role of a Rosetta Stone in this field. Some of the most important developments were obtained in the framework of 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 and gas pressure, are irrelevant in this regime. This is a valid assumption almost everywhere within the magnetosphere where the field strength is high and the plasma density is low. In this limit, the structure of the magnetosphere is then a solution of the relativistic 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 of the morphology of relativistic magnetospheres in the limit where plasma is abundant and ultramagnetized everywhere. It also captures how the star spins down under the effect of magnetic torques well. However, the main shortcoming of this approach is that it cannot capture dissipation (see, however, Gruzinov 2008 and Li et al. 2012 for a dissipative force-free formulation), and more importantly nonthermal particle acceleration and radiation. It is therefore impossible to pin down the emitting regions and to constrain the model with observations without using ad hoc prescriptions, leading to large theoretical uncertainties (Bai & Spitkovsky 2010; Kalapotharakos et al. 2014; Cao & Yang 2019).

Observations of 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 20 TeV (Ansoldi et al. 2016; H. E. S. S. Collaboration 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. 2024; Hakobyan et al. 2023; Torres et al. 2024). Global PIC simulations of pulsar magnetospheres show robust features such as a sizeable dissipation rate of the outflowing Poynting flux via magnetic reconnection of open field lines beyond the light cylinder. This energy is efficiently channelled into nonthermal particle acceleration and energetic radiation, possibly explaining the observed 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 reproduced by local PIC models (Timokhin & Arons 2013; Philippov et al. 2020; Cruz et al. 2021) well.

In spite of these impressive developments, these results must be interpreted with caution because global PIC simulations are plagued with an unrealistically small separation between the scales where microphysics operates and the system size (i.e., 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 the distant future. In this work, we made 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 up physical parameters and to check whether salient features uncovered by pure kinetic models at smaller scales still hold strong. To this end, we develop a new hybrid numerical scheme to take advantage of the complementarity between the force-free and 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 of 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 discuss future applications and developments of this new hybrid approach in Section 6.

2. Scale separation

Capturing all relevant physical scales in pulsar magnetospheres is an outstanding challenge for simulations. In the following discussion, we consider, as our baseline, a millisecond pulsar with a rotation period of P = 1 ms and a surface dipolar magnetic field of B = 107 G, corresponding to the weakest observable gamma-ray pulsars of spindown power L ≈ 4.8 × 1033erg/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, r ∼ 106 cm, and the light-cylinder radius,

R LC = cP 2 π 5 r ( P 1 ms ) . $$ \begin{aligned} R_{\rm LC}=\frac{cP}{2\pi }\approx 5r_{\star } \left(\frac{P}{1\,\mathrm{ms}}\right). \end{aligned} $$(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 it 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 prep.), and thus it sets the largest relevant scale of the problem in this work.

While the macroscopic scales are well constrained, defining the microscopic scales is more complex because plasma physics and quantum electrodynamics (QED) effects must be taken into consideration. The electron plasma skin depth at the star surface sets the minimal relevant scale of the problem. Assuming that the plasma density is a multiple of the critical co-rotation density, n = κ n GJ = κ B / 2 π e R LC $ n_{\star}=\kappa n^{\star}_{\mathrm{GJ}}=\kappa B_{\star}/2\pi e R_{\mathrm{LC}} $ (Hones & Bergeson 1965; Goldreich & Julian 1969), where κ is the plasma multiplicity and e is the electron charge, yields

d e = γ m e c 2 4 π n e 2 2 ( γ 10 0 ) 1 / 2 ( κ 10 2 ) 1 / 2 ( P 1 ms ) 1 / 2 ( B 10 7 G ) 1 / 2 cm , $$ \begin{aligned} d^{\star }_{\rm e}&=\sqrt{\frac{\gamma m_{\rm e}c^2}{4\pi n_{\star } e^2}}\nonumber \\&\approx 2\left(\frac{\gamma }{10^0}\right)^{1/2}\left(\frac{\kappa }{10^2}\right)^{-1/2}\left(\frac{P}{1\mathrm{ms}}\right)^{1/2}\left(\frac{B_{\star }}{10^7\mathrm{G}}\right)^{-1/2} \mathrm{cm}, \end{aligned} $$(2)

where γ is the particle Lorentz factor and me is the electron mass; so, d e / r = 2 × 10 6 $ d^{\star}_{\mathrm{e}}/r_{\star}=2\times 10^{-6} $ for our reference pulsar case. The situation may still look hopeless at this stage. However, particles at the neutron star surface are not at rest. Instead, they are believed to move at ultra-relativistic speeds (γ ≫ 1), thereby reducing the above scale separation. A primary beam of particles is extracted and accelerated by the surface electric field induced by the rotation of the star (Sturrock 1971; Ruderman & Sutherland 1975). At best, the particles will undergo the full vacuum potential drop across the polar cap of the star. The size of the polar cap is defined between the magnetic axis and the first field line fully contained within the light cylinder. For an aligned dipole in vacuum, its angular size is sin θ pc = r / R LC $ \sin\theta_{\mathrm{pc}}=\sqrt{r_{\star}/R_{\mathrm{LC}}} $. Using the corotation surface electric field, E = −(Ω × rB/c, we can derive the potential drop across the polar cap as

Φ pc = μ Ω 2 c 2 , $$ \begin{aligned} \Phi _{\rm pc}=\frac{\mu \Omega ^2}{c^2}, \end{aligned} $$(3)

where Ω = 2π/P is the angular velocity of the star, and μ = Br3 is its magnetic moment. An electron experiencing the full potential drop will acquire a Lorentz factor given by

γ pc = e Φ pc m e c 2 2.6 × 10 8 ( B 10 7 G ) ( P 1 ms ) 2 . $$ \begin{aligned} \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\mathrm{G}}\right)\left(\frac{P}{1\mathrm{ms}}\right)^{-2}. \end{aligned} $$(4)

This calculation provides an upper limit; it is not a realistic estimate of the particle Lorentz factor in the bulk of the flow, because the presence of the plasma screens in large part the accelerating electric field (parallel to the magnetic field), and because of pair production. The main channel for pair production near the star surface is believed to be magnetic conversion, meaning that a 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)

χ ϵ b 0.1 , $$ \begin{aligned} \chi \equiv \epsilon b\gtrsim 0.1, \end{aligned} $$(5)

where ϵ = ℏν/mec2 is the dimensionless photon energy, and b = B / B QED $ b=\tilde{B}_{\perp}/B_{\mathrm{QED}} $ is the effective perpendicular magnetic field (see Eq. 17 below) normalized to the critical quantum field BQED = me2c3/ℏe ≈ 4.4 × 1013 G. If electrons radiate via synchrotron-curvature radiation as usually assumed, the critical photon energy is given by (Blumenthal & Gould 1970)

ϵ c = 3 2 b γ 2 . $$ \begin{aligned} \epsilon _{\rm c}=\frac{3}{2}b\gamma ^2. \end{aligned} $$(6)

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

γ th = 1 15 b 2 10 6 ( B 10 7 G ) 1 , $$ \begin{aligned} \gamma _{\rm th}=\sqrt{\frac{1}{15 b^2}}\approx 10^6\left(\frac{B_{\star }}{10^7\mathrm{G}}\right)^{-1}, \end{aligned} $$(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 = 1 20 b 2.2 × 10 5 ( B 10 7 G ) 1 . $$ \begin{aligned} \gamma _{\rm s}=\frac{\epsilon }{2}=\frac{1}{20b}\approx 2.2\times 10^5 \left(\frac{B_{\star }}{10^7\mathrm{G}}\right)^{-1}. \end{aligned} $$(8)

Assuming that pair creation is effective at producing a high-multiplicity plasma in the magnetosphere (κ ≫ 1), it sets the relevant scale of the particle Lorentz factor to be considered in estimating the plasma skin depth scale in Eq. (2), yielding d e s = d e ( γ s ) 10 3 cm = 10 3 r $ d^{\mathrm{s}}_{\mathrm{e}}=d^{\star}_{\mathrm{e}}(\gamma_{\mathrm{s}})\approx 10^3\mathrm{cm}=10^{-3}r_{\star} $ for our reference pulsar. One can realize from these order-of-magnitude calculations that a weaker magnetic field (B ≲ 107 G) would prevent pair formation as the threshold energy scale would be larger than the full potential drop.

At the light cylinder, a thin current sheet forms whose thickness, δ, is governed by the local plasma skin depth scale, meaning that δ d e LC $ \delta \sim d^{\mathrm{LC}}_{\mathrm{e}} $. 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 = B LC 2 4 π Γ LC κ n GJ LC m e c 2 = γ pc 2 Γ LC κ , $$ \begin{aligned} \gamma _{\rm LC}\sim \sigma _{\rm LC}=\frac{B^2_{\rm LC}}{4\pi \Gamma _{\rm LC}\kappa n^\mathrm{LC}_{\rm GJ}m_{\rm e}c^2}= \frac{\gamma _{\rm pc}}{2\Gamma _{\rm LC}\kappa }, \end{aligned} $$(9)

where BLC ∼ B(r/RLC)3 = μΩ3/c3, n GJ LC = Ω B LC / 2 π e c $ n^{\mathrm{LC}}_{\mathrm{GJ}}=\Omega B_{\mathrm{LC}}/2\pi e c $ is the Goldreich-Julian plasma density at the light cylinder and ΓLC is the pulsar wind bulk Lorentz factor at its base. Assuming that these energetic particles set the local skin depth scale, this yields

δ d e LC ( γ LC ) = γ LC m e c 2 4 π κ Γ LC n GJ LC e 2 = R LC 2 Γ LC κ , $$ \begin{aligned} \delta \sim d^\mathrm{LC}_{\rm e}(\gamma _{\rm LC}) = \sqrt{\frac{\gamma _{\rm LC}m_{\rm e}c^2}{4\pi \kappa \Gamma _{\rm LC} n^\mathrm{LC}_{\rm GJ}e^2}}=\frac{R_{\rm LC}}{2\Gamma _{\rm LC}\kappa }, \end{aligned} $$(10)

δ 2.4 × 10 4 ( Γ LC 10 0 ) 1 ( κ 10 2 ) 1 ( P 1 ms ) cm , $$ \begin{aligned} \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\mathrm{ms}}\right)\mathrm{cm}, \end{aligned} $$(11)

thus, δ/r ∼ 2 × 10−2 or δ/RLC ∼ 5 × 10−3 for these fiducial parameters. A similar estimate can be obtained using Ampère law (Coroniti 1990; Lyubarsky & Kirk 2001; Cerutti & Philippov 2017) (see also Arka & Dubus 2013; Uzdensky & Spitkovsky 2014 for other estimates).

Radiative cooling sets another energy scale in the problem at which the accelerating electric force is balanced by the radiation-reaction force, eE ∼ frad. The latter effectively acts as a continuous drag force (in the classical regime, i.e., γb ≪ 1), opposite to the particle’s direction of motion, whose magnitude in the ultrarelativistic regime (γ ≫ 1) can be estimated as f rad ( 2 / 3 ) r e 2 γ 2 B 2 $ f_{\mathrm{rad}}\approx(2/3)r^2_{\mathrm{e}}\gamma^2 \tilde{B}_{\perp}^2 $, where re is the classical radius of the electron. The fiducial radiation-reaction-limited electron Lorentz factor is then given by

γ rad = 3 e E 2 r e 2 B 2 3 × 10 4 ( E B ) 1 / 2 ( B 10 7 G ) 1 / 2 , $$ \begin{aligned} \gamma _{\rm rad}=\sqrt{\frac{3eE_{\parallel }}{2 r^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\mathrm{G}}\right)^{-1/2}, \end{aligned} $$(12)

or γ rad LC 3.3 × 10 5 $ \gamma^{\mathrm{LC}}_{\mathrm{rad}}\approx 3.3\times 10^5 $ at the light cylinder for the fiducial pulsar parameters, assuming E = B $ E_{\parallel}=\tilde{B}_{\perp} $ and B = B LC $ \tilde{B}_{\perp}=B_{\mathrm{LC}} $.

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 r:

d e s r 10 3 δ r 2.5 × 10 2 1 < R LC r = 5 , $$ \begin{aligned} \frac{d^\mathrm{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, \end{aligned} $$(13)

and in terms of energy scales normalized to γpc,

γ s γ pc 8.5 × 10 4 < γ rad LC γ pc 10 3 < γ th γ pc 4 × 10 3 1 . $$ \begin{aligned} \frac{\gamma _{\rm s}}{\gamma _{\rm pc}}\sim 8.5\times 10^{-4}< \frac{\gamma ^\mathrm{LC}_{\rm rad}}{\gamma _{\rm pc}}\sim 10^{-3} <\frac{\gamma _{\rm th}}{\gamma _{\rm pc}}\sim 4\times 10^{-3}\ll 1. \end{aligned} $$(14)

This exercise reveals that in such systems γs, γrad, and γth are very close together, meaning that the energy range of the particle spectrum may be quite narrow. This property is in agreement with a new important feature highlighted by the most recent analysis of the Fermi-LAT pulsars, which shows that the observed spectral energy distributions are narrower with decreasing spindown, nearly reaching a spectrum consistent with a monoenergetic particle population for the weakest pulsars with a spindown of L ∼ 1033 erg/s (Smith et al. 2023). Thus, even though the estimates presented in this section by no means represent a rigorous calculation, it provides an approximate assessment of the minimum computational needs to approach realistic systems using a kinetic approach. Our conclusion is that such a weak pulsar could in principle be captured with real scales using current computational facilities, provided that extreme resolutions are used in a 2D axisymmetric system, as, for example, in the recent works by Hu & Beloborodov (2022) and Bransgrove et al. (2023).

3. Hybrid force-free-PIC approach

In order to fulfill the minimum scale separation requirements highlighted in the previous section, a large number of grid cells must be used, and, hence, high computational power is needed. To alleviate some of this numerical cost, we developed a new hybrid approach specially designed for aligned magnetospheres, coupling time-dependent force-free electrodynamic and 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., κ ≫ 1). Using the expensive full PIC technique with a large number of particles per cell everywhere in the magnetosphere may seem excessive, as a numerically cheaper force-free description would be sufficient. In contrast, dissipative regions (e.g., electrostatic gaps, current sheets) should remain on microscopic scales (again, if κ ≫ 1), and so these are the scales for which the PIC approach is really needed. Therefore, it appears desirable to combine these two techniques, but a new scheme must be implemented to ensure that both approaches can coexist harmoniously.

In this work, we developed a new 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 – simply referred to as particles in the following– representing a very large number of physical particles with identical mass-to-charge ratios 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 u d t = q mc ( E + u × B γ ) + g , $$ \begin{aligned} \frac{\mathrm{d} \mathbf u }{\mathrm{d} t} = \frac{q}{mc}\left(\mathbf E +\frac{\mathbf u \times \mathbf B }{\gamma }\right)+\mathbf g , \end{aligned} $$(15)

where u = γv/c is the dimensionless particle four-momentum vector, and q is the particle electric charge, and where

g 2 3 r e 2 [ ( E + β × B ) × B + ( β · E ) E ] 2 3 r e 2 γ 2 B 2 β $$ \begin{aligned} \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 } \end{aligned} $$(16)

is the radiation reaction force following the Landau & Lifshitz (1971) formulation, and β = v/c, and where

B = ( E + β × B ) 2 ( β · E ) 2 $$ \begin{aligned} \tilde{B}_{\perp }=\sqrt{\left(\mathbf E +\boldsymbol{\beta }\times \mathbf B \right)^2-\left(\boldsymbol{\beta }\cdot \mathbf E \right)^2} \end{aligned} $$(17)

is the effective perpendicular magnetic field (Cerutti et al. 2016). Particle trajectories are evolved in time with a modified Boris pusher to incorporate the radiation reaction force to the Lorentz force. The current density carried by each particle is deposited onto the numerical grid where the electromagnetic fields are defined. Summing over all particles and over all species, the total current J reconstructed on the grid is then used to evolve the time-dependent Maxwell equations:

B t = c × E $$ \begin{aligned} \frac{\partial \mathbf B }{\partial t} = -c\,\boldsymbol{\nabla }\times \mathbf E \end{aligned} $$(18)

E t = c × B 4 π J , $$ \begin{aligned} \frac{\partial \mathbf E }{\partial t} = c\,\boldsymbol{\nabla }\times \mathbf B -4\pi \mathbf J , \end{aligned} $$(19)

allowing us 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).

thumbnail Fig. 1.

New hybrid numerical scheme proposed in this work during a time-integration cycle. This method is intended 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, J, using a blending function, f, which solely depends on the magnetic flux function, Ψ.

In general, each simulation cell is filled with a large number of particles, thus pushing the particles and depositing the currents largely dominate the overall computing time. These two steps typically take one order of magnitude longer than evolving the fields alone. This is where the benefits of force-free electrodynamics come into play. The force-free condition is given by

F μ ν I ν = 0 , $$ \begin{aligned} F_{\mu \nu }I^{\nu }=0, \end{aligned} $$(20)

where Fμν is the electromagnetic tensor and Iν is the four-current, which translates into the following conditions:

ρ E + J × B c = 0 , $$ \begin{aligned} \rho \mathbf E +\frac{\mathbf J \times \mathbf B }{c}=\mathbf 0 , \end{aligned} $$(21)

using the spatial components, where ρ =  ⋅ E/4π is the charge density; and

E · J = 0 , $$ \begin{aligned} \mathbf E \cdot \mathbf J =0, \end{aligned} $$(22)

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

E · B = 0 . $$ \begin{aligned} \mathbf E \cdot \mathbf B =0. \end{aligned} $$(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)

J = c 4 π · E ( E × B B 2 ) + c 4 π ( B · ( × B ) E · ( × E ) ) B B 2 , $$ \begin{aligned} \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}, \end{aligned} $$(24)

which can be understood as the sum of a current flowing across field lines, J, and a component flowing along field lines, J. Instead of pushing a large number of particles and depositing the current on the grid, the current is derived in a single step using Eq. (24), which saves on computing time and memory. Other than this, the Maxwell solver is unchanged whether the current is computed via the particles or via the fields. Thus, the coupling between both approaches is done via the current.

To implement the force-free module, we followed a similar numerical procedure as in Spitkovsky (2006), which consists in solely computing the perpendicular force-free current, J. The parallel component, J, prevents the formation of an electric-field component along the magnetic field at all times, E ⋅ B = 0, but computing J is a numerically delicate and cumbersome endeavor on a staggered mesh (see, however, Parfrey et al. 2012). Therefore, the electric field is renormalized at each time step to enforce the force-free condition

E E ( E · B ) B B 2 . $$ \begin{aligned} \mathbf E \leftarrow \mathbf E -\left(\mathbf E \cdot \mathbf B \right)\frac{\mathbf B }{B^2}. \end{aligned} $$(25)

The code also checks, at every time step, that the condition E2 < B2 is always satisfied; so, if it is not the case, the electric field is once again renormalized as

E B 2 E 2 E . $$ \begin{aligned} \mathbf E \leftarrow \sqrt{\frac{B^2}{E^2}}\mathbf E . \end{aligned} $$(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 that is also numerically convenient for tagging individual field lines. In ordinary spherical coordinates (r, θ), the magnetic flux function is given by

Ψ = 1 2 π B · d S = 0 θ B r r 2 sin θ d θ . $$ \begin{aligned} \Psi =\frac{1}{2\pi }\int \int \mathbf B \cdot d\mathbf S =\int _0^\theta B_\mathbf{r }r^2\sin \theta ^{\prime }d\theta ^{\prime }. \end{aligned} $$(27)

Once the boundary is set, a thin transition layer is implemented between neighboring field lines where the currents are blended to ensure a better continuity between the PIC and the force-free currents and to minimize numerical artifacts, such as spurious reflections or accumulation of charges. A simple linear interpolation scheme is sufficient for this purpose, so the full hybrid current is expressed as

J = ( 1 f ( Ψ ) ) J PIC + f ( Ψ ) J FFE , $$ \begin{aligned} \mathbf J =\Bigl (1-f\left(\Psi \right)\Bigr )\,\mathbf J _{\rm PIC}+f\left(\Psi \right)\,\mathbf J _{\rm FFE}, \end{aligned} $$(28)

where JPIC is the current density reconstructed from the particles, JFFE is the current from the force-free condition (Eq. 24), and f is the blending function. This hybrid current is then injected into the 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 < Ψ1, so that the blending function is as follows:

f ( Ψ ) = 1 , if Ψ < Ψ 0 , pure force-free, no particles , $$ \begin{aligned} f\left(\Psi \right)&=1,~\text{ if}~\Psi <\Psi _0,\rightarrow \text{ pure} \text{ force-free,} \text{ no} \text{ particles},\end{aligned} $$(29)

f ( Ψ ) = Ψ 1 Ψ Ψ 1 Ψ 0 , if Ψ 0 < Ψ < Ψ 1 , transition layer , $$ \begin{aligned} f\left(\Psi \right)&=\frac{\Psi _1-\Psi }{\Psi _1-\Psi _0},~\text{ if}~\Psi _0<\Psi <\Psi _1,\rightarrow \text{ transition} \text{ layer},\end{aligned} $$(30)

f ( Ψ ) = 0 , if Ψ > Ψ 1 , pure PIC . $$ \begin{aligned} f\left(\Psi \right)&=0,~\text{ if}~\Psi >\Psi _1,\rightarrow \text{ pure} \text{ PIC}. \end{aligned} $$(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 it dynamically evolves along with the magnetospheric activity. Figure 1 graphically summarizes the proposed numerical scheme over a time iteration.

4. Monopole

The hybrid scheme is first put to the test through modeling a 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 is free of dissipation – there is no gap or current sheet (as opposed to a split-monopole) – so spurious numerical effects in the scheme can be quantified. Moreover, it represents a good model for the asymptotic magnetospheric structure of pulsar open field lines. All things considered, it sets the best conditions to validate our new approach.

In this setup, the magnetic flux function is given by

Ψ M = Ψ ( 1 cos θ ) , $$ \begin{aligned} \Psi _{\rm M}=\Psi _{\star }\left(1-\cos \theta \right), \end{aligned} $$(32)

where Ψ = Br2. An exact solution to the Grad-Shafranov equation yields the following toroidal magnetic field component (Michel 1973):

B ϕ = Ψ R LC sin θ r , $$ \begin{aligned} B_{\phi }=-\frac{\Psi _{\star }}{R_{\rm LC}}\frac{\sin \theta }{r}, \end{aligned} $$(33)

if Ω ⋅ B > 0 in the northern hemisphere (and Ω ⋅ B < 0 in the southern hemisphere). In this solution, the electric field is purely along θ and matches the toroidal magnetic field, Eθ = Bϕ. The electric current is purely radial and is given by

J M = c 4 π × B = J GJ ( r r ) 2 cos θ e r , $$ \begin{aligned} \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 ~\mathbf e_{\rm r} , \end{aligned} $$(34)

where J GJ = Ω B / 2 π $ J^{\star}_{\mathrm{GJ}}=\Omega B_{\star}/2\pi $ corresponds to the fiducial Goldreich-Julian current, and the spindown power associated with the electromagnetic torque is

L M = c 4 π ( E × B ) · d S = c 2 1 1 r 2 B ϕ 2 d cos θ = 2 Ψ 2 Ω 2 3 c . $$ \begin{aligned} L_{\rm M}=\frac{c}{4\pi }\int \int \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}. \end{aligned} $$(35)

The goal of this first numerical experiment is to recover these properties using our hybrid scheme. To this end, we performed a set of hybrid simulations in 2D axisymmetric spherical coordinates and explored the effect of resolution on the results. The domain size explored here is composed of 10242, 20482, and 40962 cells. The grid spacing is logarithmic along the radial direction and constant along the θ direction. The domain size ranges from the star surface, rmin = r, up to rmax = (10/3)RLC, where RLC = 5r, and covers the full range in θ ∈ [0, π]. The magnetosphere is initialized in vacuum with a pure radial magnetic field, B = B(r/r)2. The surface magnetic field strength is fixed at B = 5 × 105 G, which is sufficiently intense in this experiment to achieve a quasi-force-free state in the PIC domain. The role of the field strength and scale separation is thoroughly investigated in the next section.

Corotation of the field lines with the star at the inner boundary is guaranteed by fixing the poloidal electric field to the ideal solution, E = −(Ω × rB/c. This procedure sets the magnetosphere into rotation. An absorbing layer starting at rabs = 3RLC progressively damps all outgoing electromagnetic waves and removes the particles, if there are any (Cerutti et al. 2015). In this experiment, the northern hemisphere is modeled in full PIC, while the southern hemisphere is captured with the force-free model. This choice is well suited for a comparison of both approaches thanks to the expected symmetry between both hemispheres; hence, any asymmetry or imbalance can be tracked easily. In this setup, we did not notice any difference with or without a transition layer for the two highest resolutions. Thus, we only show the runs performed with a sharp transition at θ = π/2, or Ψ = Ψ. In the PIC domain, electron-positron pairs are continuously injected with a high multiplicity, κinj = 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.

The top panel of Fig. 2 shows the state of the simulation for the highest numerical resolution, obtained after about two spin periods. A steady state is reached, and the solution is numerically stable. We recover the correct variations of the toroidal field and current densities. The latter follows the predicted cosine evolution in both hemispheres without apparent mismatch at the interface between both domains at the equator, except for the lowest resolution simulation where a small glitch is visible (middle panel of Fig. 2). This result demonstrates that both hemispheres are well balanced, meaning that there is no net accumulation of charges at the boundary or anywhere in the magnetosphere. The outgoing current in the 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

thumbnail Fig. 2.

Simulated monopolar magnetosphere using the hybrid force-free-PIC scheme. The division between the force-free and PIC domains is the equatorial plane (dashed magenta line in the top and middle panels). Top: Spatial distribution of the radial current density normalized to the fiducial current density, J GJ $ J^{\star}_{\mathrm{GJ}} $, and compensated by (r/r)2 for visualization purposes, for the 40962 cells run. Solid lines show the poloidal magnetic field lines. Middle: Cross-section of the current density profile at r = 1.5RLC 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), normalized to the monopole solution, LM (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.

ϵ = | L in L out | L out , $$ \begin{aligned} \epsilon =\frac{\left|L_{\rm in}-L_{\rm out}\right|}{L_{\rm out}}, \end{aligned} $$(36)

where Lin is the Poynting flux averaged in the r ∈ [r, 2r] interval, and Lout is averaged over the r ∈ [2RLC, 3RLC] interval (gray intervals in Fig. 2, bottom panel). We observe numerical convergence and a deviation smaller than ϵ < 1% for the 10242 run, down to ϵ ∼ 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 models of pulsar magnetospheres over a similar radial scale (e.g., Philippov & Spitkovsky 2014; Cerutti et al. 2015; Belyaev 2015; Cerutti et al. 2020; Hakobyan et al. 2023). Thus, the hybrid scheme does not add significantly more artificial dissipation to the simulated system, at least under the ideal condition of a force-free monopole.

5. Aligned dipole

The hybrid model is now applied to the aligned rotating dipole using 2D axisymmetric simulations, aiming at scaling relevant physical scales up to the fiducial weak gamma-ray pulsars introduced in Sect. 2.

5.1. Numerical setup

The properties of the simulation domain are nearly identical to the monopole configuration presented in the previous section (i.e., domain size, grid spacing and geometry, boundary conditions, pulsar spin period). Here, we focus on the highest resolution simulations that are composed of 81922 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 = μ sin 2 θ r . $$ \begin{aligned} \Psi _{\rm D}=\frac{\mu \sin ^2\theta }{r}. \end{aligned} $$(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 = 1836me. 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 γth, then a new pair is created with a Lorentz factor γs at the same position and along the same direction of motion as the parent particle. This procedure is similar to the heuristic models of pair production performed in previous studies (Philippov et al. 2015b; Chen et al. 2020; Cruz et al. 2021). It provides an effective way to inject pairs in the magnetosphere with the relevant energy scale, but this model should not be viewed as very realistic. Pair production in millisecond pulsars remains a puzzle that we do not wish to address in this work. In the 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. A natural choice is thus to limit the PIC domain to the equatorial and separatrix regions. The boundaries of the domain should be chosen with care. It is particularly important that the PIC domain extends from the star surface all the way to the absorbing boundary. If the chosen field line separating the PIC and the 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 defeat 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Ψpc. This region corresponds to the open field lines connected to the star polar caps.

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

  • Force-free domain II: 2.3Ψpc < Ψ < Ψmax. It is inside the closed field line regions deep inside the light cylinder.

In the above, Ψmax = μ/r is the maximum of the flux function that is located at the star’s surface in the equatorial plane, and Ψpc = μ/RLC 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-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 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 = 106 cm, and the light-cylinder radius is set to RLC = 5 × 106 cm, corresponding to a P = 1 ms spin period. The surface magnetic field of the reference production run is fixed at B = 107 G, giving BLC = μ/RLC3 ≈ 8 × 104 G at the light cylinder, and a force-free spindown power of L0 = μ2Ω4/c3 = 4.8 × 1033 erg/s. The magnetic field strength then sets the relevant energy scales given in Sect. 2 as γpc = 2.6 × 108, γth = 106, γs = 2.2 × 105, and γrad = 3 × 104. The Larmor radius and frequency scales are not resolved near the surface of the star, but they are irrelevant since the particles stream along magnetic field lines with negligible perpendicular momentum due to the strong radiative cooling. However, the plasma skin depth and frequency scales are resolved by d e s / Δ r 3 $ d^{\mathrm{s}}_{\mathrm{e}}/\Delta r\approx 3 $ and ω pe 1 / Δ t 8 $ \omega^{-1}_{\mathrm{pe}}/\Delta t\approx 8 $ at the surface where the conditions are numerically the most demanding, and they are well resolved elsewhere in the magnetosphere. The Larmor scales are also resolved in the outer parts of the magnetosphere, in particular in the equatorial current layer, where these scales become relevant. The Larmor radius of secondary pairs at the light cylinder, ρs, divided by the local grid spacing, RLCΔr/r, is ρs/(RLCΔr/r) = (γsmec2/eBLC)/(RLCΔr/r)≈3, and Larmor timescale ρs/(cΔt)≈36. The corresponding synchrotron cooling time is on the order of t syn LC 3 m e c / ( 2 r e 2 γ s B LC 2 ) 85 Δ t $ t^{\mathrm{LC}}_{\mathrm{syn}}\sim 3m_{\mathrm{e}}c/(2r^2_{\mathrm{e}}\gamma_{\mathrm{s}}B^2_{\mathrm{LC}})\approx 85\Delta t $.

To prevent spurious numerical effects near the star’s surface due to strong cooling in the integration of particle trajectories, any particle momentum transverse to the field is removed at each time step within an r-thick shell around the star. The particle pusher is unchanged everywhere else in the simulation. The radiation reaction force must, however, be artificially and momentarily reduced in the early stages of the simulation. The transient phase, when the magnetosphere is filling up with plasma, is quite violent and can generate extreme field strengths and gradients, leading to a crash in the simulation. The simulations are first evolved during one period where the strength of the radiation-reaction force is set to 1% of its nominal value, and it then slowly increases in time up to about 10%. In a last step, it is gradually increased up to the nominal strength at t ≳ 2P when a 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 feasible. Here, we scaled the magnetic field strength of the reference simulation down to B = 104 G, 105 G, and 106 G, but we kept the ratio between all energy scales the same as in the B = 107 G simulation, including the radiation-reaction-limited lepton Lorentz factor, γrad, given in Eq. (14). To keep the γpc/γrad ratio identical to that of the B = 107 G run, the strength of the radiation reaction force –or equivalently the Thomson cross-section– must be boosted by a factor of frad = 109 for B = 104 G, frad = 106 for B = 105, and frad = 103 for B = 106 G. Table 1 lists the numerical and physical parameters used in this study for the reference simulation.

Table 1.

Physical and numerical parameters for the reference weak gamma-ray pulsar simulation using the hybrid force-free-PIC scheme.

5.2. Results

Figure 3 shows the overall magnetospheric structure achieved by the hybrid model, once the initial transient has left the box and the full-scale radiation-reaction force is established at t ≈ 2.5P for the fiducial full-scale B = 107 G millisecond pulsar. The solution is consistent with previous global PIC studies. It is characterized by prominent separatrices and reconnecting current layers whose base stands near the light cylinder (the Y-point). It is important to note that this solution has not completely reached a steady state. During the initial transient phase, the Y-point forms well inside the light cylinder (around 0.6RLC), slowly migrates toward RLC, and 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. 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 us to bridge the pure force-free region to these vacuum gaps without noticeable numerical artifacts. We also verified that this feature is robust against the latitudinal extent of the PIC domain. While plasma depletion is expected in this regime, it is probably exaggerated here because of our pair production scheme, where the 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.

thumbnail Fig. 3.

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

Perhaps the most spectacular feature is the extreme narrowness of the current layer, the latter being barely visible in the global view (see the zoomed-in view in Figure 3). Near the light cylinder, the layer’s thickness is about δ ∼ 10−3RLC, which corresponds to the scale separation anticipated in Sect. 2, indicating that the simulation has achieved the desired scale separation. The pair multiplicity in the layer is κ ∼ 10 − 100, while the proton density is at most of order unity. The pairs being much more abundant than the protons, the layer’s thickness is controlled by the electronic scales. In contrast, the proton beam is much thicker, with a width of ∼0.1RLC near the light cylinder, and it is not significantly perturbed by the presence of plasmoids. The decoupling between the proton and the pair scales is amplified by radiative cooling, from which the protons are nearly immune.

Figure 3 also maps out the mean particle energies in the magnetosphere, demonstrating the crucial role of reconnection in particle acceleration. 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 these X-points present in the layer (see the inset in Figure 3) shows that the energy of the pairs peaks inside the layer, ℰpairs = γmec2 ≈ 1 TeV, and sharply decreases inside the nearest plasmoids on either side, down to ℰpairs ≲ 100 GeV. In this extreme synchrotron cooling regime, particles accelerate deep inside the layer where the effective magnetic-field strength is low, allowing them to accelerate beyond the radiation-reaction-limited energy, γ rad LC m e c 2 150 $ \gamma^{\mathrm{LC}}_{\mathrm{rad}}m_{\mathrm{e}}c^2\approx 150 $ GeV≪ℰpairs; however, as soon as they leave the layer and enter a plasmoid, they feel a sharp increase of the perpendicular magnetic field, they cool off catastrophically, and they radiate all their energy away in the form of synchrotron photons (Cerutti et al. 2013). Plasmoids effectively act as beam dumps, as clearly visible in Figure 3. Thus, other acceleration processes involving, for instance, plasmoid contraction or particle escape from plasmoids are irrelevant for the pairs in this environment, in contrast to 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 of ℰ ≈ 10 TeV, which is 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 × 130 TeV, where βrec ∼ 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 1 GeV to 1 TeV and peaks at about a few tens of gigaelectronvolt, below the radiation-reaction-limited energy ( γ rad LC m e c 2 150 $ \gamma^{\mathrm{LC}}_{\mathrm{rad}}m_{\mathrm{e}}c^2\approx 150 $ GeV). While there is no major difference between electrons and positrons at low energies, positrons significantly dominate over the electrons at the highest energies. Positrons clearly show an extra spectral component peaking at 300–400 GeV, which coincides with particles energized at X-points whose energy is determined by the magnetization parameter σLCmec2 ≈ 500 GeV (Eq. 9). This feature was already reported in previous works (Cerutti et al. 2015), and it seems to hold strong on more realistic scales. This asymmetry is due to the overall polarization of the magnetosphere, where a positive net charge is expected in the equatorial regions if Ω ⋅ μ > 0. We note that this effect is strongest for an aligned rotator and vanishes for an orthogonal rotator (Philippov & Spitkovsky 2018). As anticipated in Sect. 2, the pair spectrum is relatively narrow in energy band for a weak millisecond gamma-ray pulsar because of the proximity between the relevant energy scales. The proton spectrum is composed of a single component ranging from about 1 TeV up to 30 TeV. If we venture to fit a power law, the spectral index would be ∼ − 1.5, which is in qualitative agreement with local reconnection studies (Melzani et al. 2014; Werner et al. 2018). Although the exact spectral shape of all species must be considered with caution, given the oversimplified model used for pair production in this work, we believe that the energy scales, width, and budget are robust overall. Scaling down the magnetic field strength and keeping all energy scale separations the same with respect to the polar-cap potential drop, γpc, does not change the above picture. All magnetospheric features are recovered, and the particle maximum energies scale linearly with the field strength (Figure 5). This result suggests that the rescaling performed in previous studies is legitimate as long as the hierarchy of scales of the targeted system is respected.

thumbnail Fig. 4.

Total electron (solid blue), positron (solid red), and proton (dashed green) energy spectra for a B = 107 G and P = 1 ms pulsar. The relevant energy scales introduced in Sect. 2 are shown by vertical dotted lines, from the lowest to the highest energies: the radiation-reaction-limited energy at the light cylinder, γ rad LC m e c 2 $ \gamma^{\mathrm{LC}}_{\mathrm{rad}}m_{\mathrm{e}}c^2 $, the threshold energy for pair production, γthmec2, and the polar-cap potential drop, eΦpc.

thumbnail Fig. 5.

Maximum cut-off energy of the charged particle energy distributions, ℰmax (electrons: blue dots; positrons: red dots; protons: green dots), as a function of the surface magnetic field B. The evolution of the polar-cap potential drop is shown by the dashed line for comparison (Eq. 4).

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 1 MeV to 1 GeV, with a peak in the energy band of the Fermi-LAT above 100 MeV. It includes the emission from the equatorial current layer only, and thus synchrotron radiation largely dominates over curvature radiation because the Larmor radius is much smaller than the particle guiding center curvature radius (≳RLC). 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/α ≈ 160 MeV (Uzdensky et al. 2011) where α is the fine structure constant. Integrating over all frequencies, the synchrotron luminosity accounts for about 18% of the pulsar spindown. This sizeable fraction corresponds to the amount of dissipation of the Poynting flux measured in the magnetosphere within a few 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%), thus demonstrating that dissipation is by far physical and due to magnetic reconnection. This also means that nearly all of the dissipated power in the box is channelled into nonthermal radiation, which is possible only in the strong cooling regime. Indeed, about 2% of the spindown is channelled into the beam of protons while the pairs leave the box with ≲1% of the spindown power. The energy fraction carried by the pairs and the protons would increase at larger distances where dissipation proceeds and synchrotron cooling becomes inefficient. Integrated above the energy threshold of the Fermi-LAT, our fiducial pulsar shines about 4% of its total spindown power above 100 MeV range. Compared to the Fermi pulsar population, our solution is compatible with the typical reported gamma-ray efficiencies that range from 1–100% (Smith et al. 2023; see Figure 7).

thumbnail Fig. 6.

Total synchrotron spectral energy distribution emitted by the equatorial current sheet for B = 107 G and P = 1 ms (dashed black line). The electronic (positronic) contribution is shown by the blue (red) line. The vertical dotted line indicates the 100 MeV energy threshold of the Fermi-LAT for comparison.

thumbnail Fig. 7.

Pulsar gamma-ray efficiency above 100 MeV; that is, the gamma-ray luminosity over the pulsar spindown power (Lγ/L), as reported in the third Fermi-LAT catalog (gray stars; Smith et al. 2023), as a function of the pulsar spindown power (L). The properties of the pulsar simulated in this work (B = 107 G, P = 1 ms) are shown by the red dot. The error bars give a measure of the variations of the spindown and radiative efficiencies measured in a 0.5P time frame.

6. Conclusion

In this work, we developed a new hybrid scheme designed to blend the PIC and the force-free approaches in the context of relativistic magnetospheres. This effort was motivated by the need to increase the separation between the kinetic scales, where particle acceleration operates, and global scales. By this means, we increase the credibility of global PIC models in spite of their unrealistically low scale separations. Using these new capabilities and 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 us to circumvent the severe resolution constraints imposed by the physical conditions at the polar caps where the plasma is not as energetic as in the current layers, which would prevent a full PIC simulation to model such a system with the same numerical resolution as the hybrid run.

Our results recover all features reported by previous studies: a prominent reconnecting current layer efficiently accelerating pairs and ions. For the surface magnetic field, B = 107 G, simulated in this work, pairs are accelerated up to 1 TeV, while protons reach up to ≳10% of the full polar-cap potential drop, that is ≳10 TeV. With a luminosity of a few percent of the pulsar spindown power, the millisecond pulsar population could thus be a sizeable contributor to the Galactic 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% of the spindown power, including about 4% above 100 MeV 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 2023).

We also show that the rescaling procedure operated by previous PIC studies, in particular in the surface magnetic field strength, is valid as long as a realistic hierarchy of scales is respected with regard to the targeted system. Thus, the results from this simulation may not be extrapolated to larger field strength, as, for instance, the young pulsar population, because the scale separation would be different for a larger surface field strength. This work is intended to provide a proof of concept that magnetic reconnection within the equatorial current sheet is a viable scenario to explain the 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. Including a more realistic model for pair production could significantly improve the predictive power of the simulations, in particular regarding the pair multiplicities in the magnetosphere, the particle or photon spectral shapes and cut-off energies, as well as the radiative efficiencies. These last two quantities reported in this work are 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.

Acknowledgments

We thank Kyle Parfrey and Jens Mahlmann for their valuable comments and encouragements during the implementation of this new numerical scheme. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 863412). Computing resources were provided by TGCC under the allocation A0130407669 made by GENCI. We thank the International Space Science Institute (ISSI) for providing financial support for the organization of the meeting of ISSI Team No 459 led by I. Contopoulos and D. Kazanas where this project has been discussed. Software: Matplotlib (Hunter 2007), Numpy (Harris et al. 2020).

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460 [Google Scholar]
  2. Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17 [Google Scholar]
  3. Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Arka, I., & Dubus, G. 2013, A&A, 550, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Arons, J. 2003, ApJ, 589, 871 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bai, X.-N., & Spitkovsky, A. 2010, ApJ, 715, 1282 [Google Scholar]
  7. Belyaev, M. A. 2015, MNRAS, 449, 2759 [NASA ADS] [CrossRef] [Google Scholar]
  8. Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (Boca Raton: CRC Press) [CrossRef] [Google Scholar]
  9. Blandford, R. D. 2002, in Lighthouses of the Universe: The Most Luminous Celestial Objects and Their Use for Cosmology, eds. M. Gilfanov, R. Sunyeav, & E. Churazov, 381 [Google Scholar]
  10. Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
  12. Blasi, P., Epstein, R. I., & Olinto, A. V. 2000, ApJ, 533, L123 [NASA ADS] [CrossRef] [Google Scholar]
  13. Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
  14. Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cao, G., & Yang, X. 2019, ApJ, 874, 166 [Google Scholar]
  16. Cao, G., Zhang, L., & Sun, S. 2016, MNRAS, 455, 4267 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cerutti, B., & Beloborodov, A. M. 2017, Space Sci. Rev., 207, 111 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cerutti, B., & Philippov, A. A. 2017, A&A, 607, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [Google Scholar]
  20. Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [Google Scholar]
  21. Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
  22. Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, L22 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, ApJ, 889, 69 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chernoglazov, A., Hakobyan, H., & Philippov, A. 2023, ApJ, 959, 122 [NASA ADS] [CrossRef] [Google Scholar]
  26. Contopoulos, I. 2007, A&A, 472, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Contopoulos, I., & Stefanou, P. 2019, MNRAS, 487, 952 [CrossRef] [Google Scholar]
  28. Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351 [Google Scholar]
  29. Contopoulos, I., Pétri, J., & Stefanou, P. 2020, MNRAS, 491, 5579 [NASA ADS] [CrossRef] [Google Scholar]
  30. Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. N. 2024, MNRAS, 527, L127 [Google Scholar]
  31. Coroniti, F. V. 1990, ApJ, 349, 538 [Google Scholar]
  32. Crinquand, B., Cerutti, B., & Dubus, G. 2019, A&A, 622, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Cruz, F., Grismayer, T., Chen, A. Y., Spitkovsky, A., & Silva, L. O. 2021, ApJ, 919, L4 [NASA ADS] [CrossRef] [Google Scholar]
  34. Cruz, F., Grismayer, T., Chen, A. Y., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202347926 [Google Scholar]
  35. Erber, T. 1966, Rev. Mod. Phys., 38, 626 [NASA ADS] [CrossRef] [Google Scholar]
  36. Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118 [NASA ADS] [CrossRef] [Google Scholar]
  37. Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869 [Google Scholar]
  38. Gruzinov, A. 1999, ArXiv e-prints [arXiv: astro-ph/9902288] [Google Scholar]
  39. Gruzinov, A. 2008, ArXiv e-prints [arXiv:0802.1716] [Google Scholar]
  40. Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Gunn, J. E., & Ostriker, J. P. 1969, Phys. Rev. Lett., 22, 728 [NASA ADS] [CrossRef] [Google Scholar]
  42. H.E.S.S. Collaboration (Aharonian, F., et al.) 2023, Nat. Astron., 7, 1341 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, 943, 105 [NASA ADS] [CrossRef] [Google Scholar]
  44. Harding, A. K., & Lai, D. 2006, Rep. Prog. Phys., 69, 2631 [NASA ADS] [CrossRef] [Google Scholar]
  45. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  46. Hones, E. W., Jr, & Bergeson, J. E. 1965, J. Geophys. Res., 70, 4951 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hu, R., & Beloborodov, A. M. 2022, ApJ, 939, 42 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  49. Kalapotharakos, C., Harding, A. K., & Kazanas, D. 2014, ApJ, 793, 97 [Google Scholar]
  50. Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [NASA ADS] [CrossRef] [Google Scholar]
  51. Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261 [Google Scholar]
  52. Komissarov, S. S. 2002, MNRAS, 336, 759 [Google Scholar]
  53. Komissarov, S. S. 2006, MNRAS, 367, 19 [Google Scholar]
  54. Landau, L. D., & Lifshitz, E. M. 1971, in The Classical Theory of Fields, eds. L. D. Landau, & E. M. Lifshitz (New York: Pergamon Press) [Google Scholar]
  55. Li, J., Spitkovsky, A., & Tchekhovskoy, A. 2012, ApJ, 746, 60 [Google Scholar]
  56. Lyubarsky, Y. 2020, ApJ, 897, 1 [NASA ADS] [CrossRef] [Google Scholar]
  57. Lyubarsky, Y., & Kirk, J. G. 2001, ApJ, 547, 437 [NASA ADS] [CrossRef] [Google Scholar]
  58. McKinney, J. C. 2006, MNRAS, 368, L30 [NASA ADS] [CrossRef] [Google Scholar]
  59. Melzani, M., Walder, R., Folini, D., Winisdoerffer, C., & Favre, J. M. 2014, A&A, 570, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Michel, F. C. 1973, ApJ, 180, L133 [Google Scholar]
  61. Michel, F. C., & Tucker, W. H. 1969, Nature, 223, 277 [CrossRef] [Google Scholar]
  62. Most, E. R., & Philippov, A. A. 2020, ApJ, 893, L6 [NASA ADS] [CrossRef] [Google Scholar]
  63. Parfrey, K., Beloborodov, A. M., & Hui, L. 2012, MNRAS, 423, 1416 [CrossRef] [Google Scholar]
  64. Petropoulou, M., & Sironi, L. 2018, MNRAS, 481, 5687 [NASA ADS] [CrossRef] [Google Scholar]
  65. Philippov, A., & Kramer, M. 2022, ARA&A, 60, 495 [NASA ADS] [CrossRef] [Google Scholar]
  66. Philippov, A. A., & Spitkovsky, A. 2014, ApJ, 785, L33 [NASA ADS] [CrossRef] [Google Scholar]
  67. Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94 [NASA ADS] [CrossRef] [Google Scholar]
  68. Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJ, 815, L19 [NASA ADS] [CrossRef] [Google Scholar]
  69. Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015b, ApJ, 801, L19 [NASA ADS] [CrossRef] [Google Scholar]
  70. Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, Phys. Rev. Lett., 124, 245101 [Google Scholar]
  71. Popov, S. B., & Postnov, K. A. 2010, in Evolution of Cosmic Objects through their Physical Activity, eds. H. A. Harutyunian, A. M. Mickaelian, & Y. Terzian, 129 [Google Scholar]
  72. Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51 [Google Scholar]
  73. Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951 [Google Scholar]
  74. Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191 [NASA ADS] [CrossRef] [Google Scholar]
  75. Speiser, T. W. 1965, J. Geophys. Res., 70, 4219 [NASA ADS] [CrossRef] [Google Scholar]
  76. Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
  77. Sturrock, P. A. 1971, ApJ, 164, 529 [Google Scholar]
  78. Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255 [Google Scholar]
  79. Timokhin, A. N. 2006, MNRAS, 368, 1055 [Google Scholar]
  80. Timokhin, A. N., & Arons, J. 2013, MNRAS, 429, 20 [NASA ADS] [CrossRef] [Google Scholar]
  81. Torres, R., Grismayer, T., Cruz, F., Fonseca, R. A., & Silva, L. O. 2024, New A, 112, 102261 [NASA ADS] [CrossRef] [Google Scholar]
  82. Uzdensky, D. A. 2003, ApJ, 598, 446 [NASA ADS] [CrossRef] [Google Scholar]
  83. Uzdensky, D. A., & Spitkovsky, A. 2014, ApJ, 780, 3 [NASA ADS] [Google Scholar]
  84. Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJ, 737, L40 [NASA ADS] [CrossRef] [Google Scholar]
  85. Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [Google Scholar]
  86. Yee, K. 1966, IEEE Trans. Antennas Propagation, 14, 302 [CrossRef] [Google Scholar]
  87. Zhang, H., Sironi, L., & Giannios, D. 2021, ApJ, 922, 261 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Physical and numerical parameters for the reference weak gamma-ray pulsar simulation using the hybrid force-free-PIC scheme.

All Figures

thumbnail Fig. 1.

New hybrid numerical scheme proposed in this work during a time-integration cycle. This method is intended 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, J, using a blending function, f, which solely depends on the magnetic flux function, Ψ.

In the text
thumbnail Fig. 2.

Simulated monopolar magnetosphere using the hybrid force-free-PIC scheme. The division between the force-free and PIC domains is the equatorial plane (dashed magenta line in the top and middle panels). Top: Spatial distribution of the radial current density normalized to the fiducial current density, J GJ $ J^{\star}_{\mathrm{GJ}} $, and compensated by (r/r)2 for visualization purposes, for the 40962 cells run. Solid lines show the poloidal magnetic field lines. Middle: Cross-section of the current density profile at r = 1.5RLC 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), normalized to the monopole solution, LM (Eq. 35), for all three resolutions. Percentages correspond to the amount of numerical dissipation for each run, computed between the flux averaged in the gray areas.

In the text
thumbnail Fig. 3.

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

In the text
thumbnail Fig. 4.

Total electron (solid blue), positron (solid red), and proton (dashed green) energy spectra for a B = 107 G and P = 1 ms pulsar. The relevant energy scales introduced in Sect. 2 are shown by vertical dotted lines, from the lowest to the highest energies: the radiation-reaction-limited energy at the light cylinder, γ rad LC m e c 2 $ \gamma^{\mathrm{LC}}_{\mathrm{rad}}m_{\mathrm{e}}c^2 $, the threshold energy for pair production, γthmec2, and the polar-cap potential drop, eΦpc.

In the text
thumbnail Fig. 5.

Maximum cut-off energy of the charged particle energy distributions, ℰmax (electrons: blue dots; positrons: red dots; protons: green dots), as a function of the surface magnetic field B. The evolution of the polar-cap potential drop is shown by the dashed line for comparison (Eq. 4).

In the text
thumbnail Fig. 6.

Total synchrotron spectral energy distribution emitted by the equatorial current sheet for B = 107 G and P = 1 ms (dashed black line). The electronic (positronic) contribution is shown by the blue (red) line. The vertical dotted line indicates the 100 MeV energy threshold of the Fermi-LAT for comparison.

In the text
thumbnail Fig. 7.

Pulsar gamma-ray efficiency above 100 MeV; that is, the gamma-ray luminosity over the pulsar spindown power (Lγ/L), as reported in the third Fermi-LAT catalog (gray stars; Smith et al. 2023), as a function of the pulsar spindown power (L). The properties of the pulsar simulated in this work (B = 107 G, P = 1 ms) are shown by the red dot. The error bars give a measure of the variations of the spindown and radiative efficiencies measured in a 0.5P time frame.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.