Issue 
A&A
Volume 690, October 2024



Article Number  A229  
Number of page(s)  18  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202347926  
Published online  10 October 2024 
Particleincell simulations of pulsar magnetospheres: Transition between electrosphere and forcefree regimes
^{1}
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049001 Lisboa, Portugal
^{2}
Physics Department and McDonnell Center for the Space Sciences, Washington University in St. Louis, St. Louis, MO 63130, USA
^{3}
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
^{4}
DCTI/ISCTE Instituto Universitário de Lisboa, 1649026 Lisboa, Portugal
Received:
9
September
2023
Accepted:
19
June
2024
Aims. Global particleincell (PIC) simulations of pulsar magnetospheres are performed with volume, surface, and pairproductionbased plasma injection schemes to systematically investigate the transition between electrosphere and forcefree pulsar magnetospheric regimes.
Methods. We present a new extension of the PIC code OSIRIS that can be used to model pulsar magnetospheres with a twodimensional axisymmetric spherical grid. The subalgorithms of the code and thorough benchmarks are presented in detail, including a new firstorder current deposition scheme that conserves charge to machine precision.
Results. We show that all plasma injection schemes produce a range of magnetospheric regimes. Active solutions can be obtained with surface and volume injection schemes when using artificially large plasmainjection rates, and with pairproductionbased plasma injection for sufficiently large separation between kinematic and pairproduction energy scales.
Key words: relativistic processes / methods: numerical / pulsars: general
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Over the last decade, global kinetic simulations have proven to be essential tools for understanding the electrodynamics of pulsar magnetospheres. They have been used to study the organization of plasma currents in the vicinity of the neutron star (Philippov et al. 2015a; Chen 2017; Kalapotharakos et al. 2018) and the acceleration of leptons (Chen & Beloborodov 2014; Belyaev 2015a; Cerutti et al. 2015; Philippov & Spitkovsky 2014; Philippov et al. 2015b; Brambilla et al. 2018) and ions (Guépin et al. 2020) in the current sheets that develop beyond the light cylinder, leading to gammaray emission consistent with observations.
Particleincell (PIC) (Dawson 1962, 1983; Hockney & Eastwood 1988; Birdsall & Langdon 1991) has been the main methodology used in global kinetic simulations of pulsar magnetospheres. PIC simulations reproduce the kinetic plasma phenomena relevant in pulsars with high fidelity, including the evolution of highly nonthermal particle distributions and kineticscale fluctuations (Touati et al. 2022). Recent extensions of the PIC method have also allowed the inclusion of quantum electrodynamical effects such as pair production (Grismayer et al. 2016, 2017) and general relativity corrections (Philippov et al. 2015a; Torres et al. 2024), which are also relevant in pulsars.
Depending on the details of the injection and/or pairproduction model, the global asymptotic magnetospheric topology varies quite significantly: in some cases, the system autoregulates to a fully chargeseparated configuration (also called electrosphere) that does not produce a Poynting flux, whereas in other cases the magnetosphere converges to a forcefree regime (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Cerutti et al. 2015; Guépin et al. 2020; Hakobyan et al. 2023). While this range of solutions has been identified in several works, a systematic study has not been performed to compare volume, surface, and pairproductionbased injection schemes.
In this work, we performed twodimensional axisymmetric global simulations of pulsar magnetospheres with three different pairinjection schemes: (i) over large volumes of the magnetosphere; (ii) from the stellar surface only; and (iii) using a prescription model for pair production. We use these simulations to systematically characterize the obtained magnetospheric solutions as a function of the injection and/or pairproduction model parameters. We show that all plasma sources produce near forcefree solutions in the regime of large plasma supply and inactive electrosphere solutions with small plasma supply. All plasma sources also allow a transitional regime with subforcefree surface Poynting flux and wide equatorial current sheets.
The simulations presented in this work are performed with a recent extension of the PIC code OSIRIS (Fonseca et al. 2002, 2008), which was developed for magnetospheric models of compact objects and is presented in this work for completeness.
This paper is organized as follows. In Sect. 2, we describe the set of numerical techniques used to generalize the PIC method to perform twodimensional axisymmetric global kinetic simulations of pulsar magnetospheres with OSIRIS: the adopted discretization of the spatial domain is presented in Sect. 2.1, and the numerical schemes used to advance the field and particle equations and the corresponding boundary conditions are detailed in Sects. 2.2 and 2.3. A new chargeconserving currentdeposition scheme is presented in Sect. 2.4, and the typical scales and normalizations adopted in the code are presented in Sect. 2.5. In Sect. 3, we present simulations with volumebased (Sect. 3.1), surfacebased (Sect. 3.2), and pairproductionbased (Sect. 3.3) plasma injection. Our conclusions are presented in Sect. 4.
2. Numerical tool
2.1. Discretization and spatial grid
The numerical tool presented in this work is designed to model the global plasma environment surrounding neutron stars; that is, the spatial volume between the stellar surface and a few light cylinder radii above it, corresponding to the cylinder radius where the corotation velocity matches the speed of light, R_{LC} ≡ c/Ω. We describe this system in spherical coordinates, with the radial coordinate r measured from the center of the neutron star and the polar angle θ measured from the star’s rotation axis Ω. We assume that Ω is either parallel or antiparallel to the star’s magnetic axis μ, such that we can assume axisymmetry about Ω and derivatives with respect to the azimuthal angle ϕ can be dropped, ∂/∂ϕ = 0.
Similarly to Chen & Beloborodov (2014), Cerutti et al. (2015), we discretize the simulation domain r ∈ [r_{min}, r_{max}], θ ∈ [0, π] in a grid with N_{r} × N_{θ} cells. We adopt a regular grid spacing in θ, Δθ = π/N_{θ}, and in log r. The latter choice allows for a grid spacing that monotonically increases with r. In pulsar magnetosphere simulations, this choice favors the resolution of shorter spatial scales close to the stellar surface, where denser plasmas are expected, and relaxes it far from the neutron star, where it is less needed. The discretization in the radial direction can be formally written as
$$\begin{array}{c}\hfill log{r}_{n}=log{r}_{\mathrm{min}}+(n1)\mathrm{\Delta}\phantom{\rule{4pt}{0ex}},\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}n=1,2,...,{N}_{r}+1\phantom{\rule{4pt}{0ex}},\end{array}$$(1)
with Δ ≡ log(r_{max}/r_{min})/N_{r}. Equation (1) can be manipulated to write the useful relation r_{n} = r_{min}δ^{n − 1}, where δ ≡ (r_{max}/r_{min})^{1/Nr} is a parameter that combines all properties of the radial axis.
A schematic representation of the grid used to discretize a typical simulation domain in illustrated in Fig. 1a. The edges of grid cells are shown in black lines, and domain boundaries are highlighted in blue and dark red. The lower radial boundary coincides with the stellar surface, r_{min} = r_{*}, whereas the upper radial boundary is at r_{max}∼ tens of r_{*}, and acts as an open boundary. The θ = 0, π boundaries enforce axisymmetry, effectively serving as reflecting boundaries. More details about these boundaries are provided in Sects. 2.2 and 2.3.
Due to the large disparity between kinetic and system scales in pulsars, PIC simulations typically employ a phenomenological description of the pairproduction processes responsible for filling the pulsar magnetosphere. This description can be as simple as injecting plasma into a significant fraction of the simulation domain (Philippov & Spitkovsky 2014; Belyaev 2015a; Kalapotharakos et al. 2018; Brambilla et al. 2018), limiting this injection to a small volume close to the stellar surface (Cerutti et al. 2015; Hakobyan et al. 2023), or even considering heuristic pairproduction models (Chen & Beloborodov 2014; Philippov et al. 2015a,b; Chen et al. 2020; Guépin et al. 2020; Bransgrove et al. 2023).
Fig. 1. Schematic representation of a spherical PIC grid: Panel (a) shows the grid layout and identifies the coordinate system and boundary types, while panel (b) shows the edges of the grid cell where each field component is defined. 
In Fig. 1b, we show a schematic representation of a typical grid cell, which we label with indices (i, j) in the radial and polar directions, respectively. Cell boundaries are drawn in solid black lines, and auxiliary lines are drawn in dashed black lines. The positions where the electric and magnetic field components are defined are indicated in dark red and blue. Half integer indices i + 1/2 and j + 1/2 indicate positions defined as r_{i + 1/2} ≡ (r_{i} + r_{i + 1})/2 and θ_{j + 1/2} ≡ (θ_{j} + θ_{j + 1})/2, respectively. The grid illustrated in Fig. 1 presents two key differences with respect to a typical Cartesian grid: a) its cells have curvilinear boundaries and b) their shape and volume change across the grid. These conditions make each step of the PIC method in spherical coordinates more challenging, requiring conversions between coordinate systems in the particle pusher and adjustments in the current deposition scheme to accomodate particle shrinking/expansion in each time step. We explore these challenges and workarounds in Sects. 2.2, 2.3 and 2.4.
2.2. Electromagnetic field solver
Electric and magnetic field components are defined in the edges of the staggered grid cells indicated in Fig. 1b. This definition is analogous to that used in traditional Cartesian grids, and allows the use of the Yee algorithm (Yee 1966) to advance the electric and magnetic field in time via Maxwell’s equations,
$$\begin{array}{cc}\hfill {\mathbf{B}}^{n+1/2}& ={\mathbf{B}}^{n1/2}=c\mathrm{\Delta}t{(\mathrm{\nabla}\times \mathbf{E})}^{n}\phantom{\rule{4pt}{0ex}},\hfill \\ \hfill {\mathbf{E}}^{n+1}& ={\mathbf{E}}^{n}+c\mathrm{\Delta}t{(\mathrm{\nabla}\times \mathbf{B})}^{n+1/2}4\pi \mathrm{\Delta}t{\mathbf{j}}^{n+1/2}\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(2)
where quantities with integer/half integer superscripts are defined in integer/half integer times and Δt is the time step.
Here we adopt the same methodology as Cerutti et al. (2015), Belyaev (2015b) and use an integral form of Maxwell’s equations that avoids divergences on the polar boundaries. This integral form is obtained using Stokes’ theorem to evaluate the curl of electric and magnetic fields in a given cell as
$$\begin{array}{c}\hfill {(\mathrm{\nabla}\times \mathbf{E})}_{\mathrm{cell}}=({\displaystyle {\oint}_{{\mathcal{C}}_{\mathrm{cell}}}}\mathbf{E}\xb7\mathrm{d}{\mathcal{C}}_{\mathrm{cell}})/{\mathcal{S}}_{\mathrm{cell}}\phantom{\rule{4pt}{0ex}},\end{array}$$(3)
where 𝒞_{cell} is the contour defining the edge of that cell, 𝒮_{cell} is the corresponding area, and the closed integral and dot product have the usual definition of Stokes’ theorem. The cell label and corresponding integrations in Eq. (4) change according to the field component under consideration. For instance, we can write the radial component of ∇ × E as
$$\begin{array}{c}\hfill {{(\mathrm{\nabla}\times \mathbf{E})}_{r}}_{(i,j+1/2)}=\frac{sin{\theta}_{j+1}{{E}_{\varphi}}_{(i,j+1)}sin{\theta}_{j}{{E}_{\varphi}}_{(i,j)}}{{r}_{i}(cos{\theta}_{j}cos{\theta}_{j+1})}\phantom{\rule{4pt}{0ex}}.\end{array}$$(4)
This expression is derived by noting that, according to Eq. (2), (∇×E)_{r} should be defined in the same position as B_{r}, that is, at cell indices (i, j + 1/2). This defines the integration surface relevant to Stokes’ theorem as r = r_{i}, θ ∈ [θ_{j}, θ_{j + 1}]. The numerator and denominator in Eq. (4) then read respectively 2π(r_{i}sin θ_{j + 1}E_{ϕ}_{(i, j + 1)} − r_{i}sin θ_{j}E_{ϕ}_{(i, j)}) and 2πr_{i}^{2}(cosθ_{j} − cosθ_{j + 1}), where the 2π factor comes from the integration along ϕ. A similar calculation can be performed for all other components (Cerutti et al. 2015).
We note that at the simulation boundaries (i = {1, N_{r} + 1}, j = {1, N_{θ} + 1}), the integration regions are adapted to fit within the domain. For example, the θ integration is changed to θ ∈ [0, θ_{1 + 1/2}] and θ ∈ [θ_{Nθ + 1/2}, π] at the θ = 0 and θ = π boundaries, respectively. We also apply special rules to the field components at the boundaries; for example, in the polar boundaries we enforce the axisymmetry conditions E_{ϕ}_{(i, 1)} = E_{ϕ}_{(i, Nθ + 1)} = 0 and B_{θ}_{(i + 1/2, 1)} = B_{θ}_{(i + 1/2, Nθ + 1)} = 0. The inner radial boundary acts generally as a rotating conductor mimicking the stellar surface, whereas the outer boundary acts as a firstorder standard Mur open boundary condition (Mur 1981), that is, a perfect absorber of perturbations propagating perpendicularly to the boundary, thus allowing the system to converge to a quasistationary solution. For the sole purpose of code benchmarking, we have also implemented static conductor boundary conditions for both inner and outer radial boundaries, which enforce tangent (normal) electric (magnetic) field components to be null, E_{ϕ}_{(1, j)} = E_{ϕ}_{(Nr + 1, j)} = 0, E_{θ}_{(1, j + 1/2)} = E_{θ}_{(Nr + 1, j + 1/2)} = 0 and B_{r}_{(1, j + 1/2)} = B_{r}_{(Nr + 1, j + 1/2)} = 0.
We benchmarked our field solver implementation by studying stationary electromagnetic TM modes between two spherical static conductors (Jackson 1975). We verified that the solution obtained numerically is in excellent agreement with the analytical solution of Maxwell’s equations for these modes, as well as with the detailed discussion about a similar solver in Belyaev (2015b).
2.3. Particle pusher
Particle position and momentum components are updated in Cartesian coordinates with either the Boris (Boris 1970; Birdsall & Langdon 1991) or Vay (Vay 2008) pushers, although other pushers are also compatible with the remaining modified subalgorithms of PIC presented in this work. In each time step, a particle push is done as follows: first, the electric and magnetic fields are interpolated from the edges of the corresponding grid cell to the particle position x_{p}^{n} ≡ (r_{p}, θ_{p}), an operation that we write schematically as (E_{(i, j)}^{n}, B_{(i, j)}^{n})→(E_{p}^{n}, B_{p}^{n}). This interpolation is done using a area/volume weighting scheme. For example, the toroidal component of the electric field can be written as
$$\begin{array}{c}\hfill {{\mathbf{E}}_{\varphi}}_{p}={\displaystyle \sum _{{i}^{\prime}=i,i+1}}{\displaystyle \sum _{{j}^{\prime}=j,j+1}}{{f}_{r}}_{{i}^{\prime}}{{f}_{\theta}}_{{j}^{\prime}}{{\mathbf{E}}_{\varphi}}_{({i}^{\prime},{j}^{\prime})}\phantom{\rule{4pt}{0ex}},\end{array}$$(5)
with
$$\begin{array}{cc}\hfill {{f}_{r}}_{i}& =1{{f}_{r}}_{i+1}=\frac{{r}_{p}^{3}{r}_{i}^{3}}{{r}_{i+1}^{3}{r}_{i}^{3}}\phantom{\rule{4pt}{0ex}},\hfill \\ \hfill {{f}_{\theta}}_{j}& =1{{f}_{\theta}}_{j+1}=\frac{cos{\theta}_{j}cos{\theta}_{p}}{cos{\theta}_{j}cos{\theta}_{j+1}}\phantom{\rule{4pt}{0ex}}.\hfill \end{array}$$(6)
After the interpolation, the field components are converted from spherical to Cartesian coordinates, (E_{p}^{n}, B_{p}^{n})→(E_{p, C}^{n}, B_{p, C}^{n}), a calculation that depends on the particle position at time t^{n}, x^{n}. Finally, the particle momentum and position are updated in time, u^{n − 1/2} ≡ p^{n − 1/2}/m_{e}c → u^{n + 1/2} ≡ p^{n + 1/2}/m_{e}c and x^{n} → x^{n + 1} respectively. Choosing to advance position and momentum components in Cartesian coordinates guarantees that we are solving the simplest possible equations of motion and also allows for an easy integration with other modules in OSIRIS, such as those accounting for classical radiation reaction losses (Vranic et al. 2016) and QED effects (Grismayer et al. 2016, 2017), not explored in the current study. Other works seem to point out that the radiation reaction does not significantly modify the global magnetospheric structure, leading only to thinner current sheets (Uzdensky & McKinney 2011; Philippov et al. 2015a; Hakobyan et al. 2019, 2023; Schoeffler et al. 2023). This effect may be important for the hot population near the light cylinder radius, which in the current study may possess unphysical Larmor radii. Further investigation on the effect of synchrotron cooling on the magnetospheric state is deferred to future work.
We note that advancing the particle position in (x, y, z) does not introduce any asymmetry in the azimuthal direction ϕ; in fact, each macroparticle in our simulation represents a charged ring with azimuthal symmetry and ϕ is never used throughout the rest of the numerical scheme. In summary, the simulation stores the particle’s 3D Cartesian coordinates and velocities, and the 2D spherical poloidal coordinates,
$$\begin{array}{cc}& {\mathbf{x}}_{p}\equiv ({r}_{p},{\theta}_{p},{x}_{p},{y}_{p},{z}_{p}),\hfill \end{array}$$(7)
$$\begin{array}{cc}\hfill & {\mathbf{u}}_{p}\equiv ({u}_{x},{u}_{y},{u}_{z}).\hfill \end{array}$$(8)
We tested our implementation of the particle pushers in a large set of background electric and/or magnetic field configurations. In Fig. 2, we show results from a relevant subset of these configurations, namely a particle moving in (a) a uniform azimuthal magnetic field, (b) crossed constant magnetic and electric fields, and (c) the timevarying electric and magnetic field components of the TM modes described in the electromagnetic field solver benchmark presented in Sect. 2.2. For all these cases, we show a comparison between the solutions obtained with the Boris pusher and analytical or other numerical solutions. We obtain an excellent agreement between the results of the Boris pusher and the reference analytical/numerical curves. Solutions obtained with the Vay pusher show a similar agreement with the reference curves. In Fig. 2a2, we represent the temporal evolution of the particle energy for over ∼1000 periods, showing that it is conserved to machine precision. We note that in all these benchmarks, the only electromagnetic fields were those either imposed externally or calculated with the field solver; that is, they do not include the fields selfconsistently created due to particle motion via plasma currents.
Fig. 2. Particle pusher benchmarks corresponding to particle motions in (a12) a uniform azimuthal magnetic field, (b12) crossed constant magnetic and electric fields, and (c12) the timevarying electric and magnetic field components of TM modes. 
2.4. Current deposition
A current deposition algorithm computes the current density j on the edges of grid cells as the positions and momenta of particles are updated. A trivial choice is to compute this current as the sum over the macroparticles of the product of their charge density and instantaneous velocity. However, such algorithm in general does not satisfy the continuity equation (Villasenor & Buneman 1992; Esirkepov 2001),
$$\begin{array}{c}\hfill \frac{\partial \rho}{\partial t}+\mathrm{\nabla}\xb7\mathbf{j}=0\phantom{\rule{4pt}{0ex}},\end{array}$$(9)
where ρ is the total plasma density. Solving Eq. (9) ensures also that Gauss’ law, written as
$$\begin{array}{c}\hfill \mathrm{\nabla}\xb7\mathbf{E}=4\pi \rho \phantom{\rule{4pt}{0ex}},\end{array}$$(10)
is satisfied. Finding a current deposition algorithm that satisfies Eq. (9), and consequently Eq. (10) (i.e., a chargeconserving current deposition algorithm) is one of the key challenges in PIC codes. For Cartesian grids, there is a well established method for any interpolation order proposed in Esirkepov (2001). However, for nonuniform spherical grids, this challenge is more substantial, as grid cells (and particle shapes, which we define below) change across the grid. Other codes that adopt such grids (Cerutti et al. 2015; Philippov et al. 2015a) usually do not include chargeconserving current deposition algorithms, and adopt instead numerical schemes to enforce the validity of Eq. (10), for example Poisson solvers.
Here, we propose a new current deposition scheme that conserves charge to machine precision in the nonuniform grid defined in section 2.1. We start by defining the volume occupied by a macroparticle centered at (r_{p}, θ_{p}). The function that defines this volume is usually called the particle shape, S(r, θ, r_{p}, θ_{p}). Before writing the exact form of S, let us define some of its important properties, which we illustrate schematically in Fig. 3. First, the particle shape should only coincide with the shape of the cell in which its center is located, labeled with indices (i, j), when and only when (r_{p}, θ_{p}) = (r_{i + 1/2}, θ_{j + 1/2}). Since the grid spacing in the radial direction is a function of r, the particle width in this direction should also be a function of r_{p}, Δr ≡ Δr(r_{p}). Furthermore, the charge density associated with each macroparticle should also be a function of r_{p}. More specifically, the charge density should decrease with r_{p} to compensate the corresponding increase in volume of the macroparticle, such that its total charge remains constant.
Fig. 3. Schematic representation of (a) the spherical particle shape and (b) the variation of its flattop density value with the radial coordinate. The blue shaded region in (a) represents the particle shape and identifies its widths in the radial and polar directions. 
Defining the number of real particles in a macroparticle as N_{p}, we formally wish to find a waterbaglike particle number density n(r) such that
$$\begin{array}{c}\hfill {\displaystyle {\int}_{{V}_{i}}}n({r}_{i+1/2})\phantom{\rule{4pt}{0ex}}\mathrm{d}{V}_{i}={\displaystyle {\int}_{{V}_{{i}^{\prime}}}}n({r}_{{i}^{\prime}+1/2})\phantom{\rule{4pt}{0ex}}\mathrm{d}{V}_{{i}^{\prime}}={N}_{p}\phantom{\rule{4pt}{0ex}},\end{array}$$(11)
where V_{i, i′} are the volumes of cells with radial labels i, i′ (see Figure 3b)). For simplicity, we assume that the particle density is only a function of r, and generalize it later to include the natural dependence in θ as well. Assuming that n(r_{i + 1/2}) is constant within cell i, we can solve Eq. (11) to obtain
$$\begin{array}{c}\hfill n({r}_{i+1/2})=\frac{3{N}_{p}}{4\pi}\frac{1}{{r}_{i+1}^{3}{r}_{i}^{3}}=\frac{3{N}_{p}}{32\pi}\frac{{(\delta +1)}^{3}}{{\delta}^{3}1}\frac{1}{{r}_{i+1/2}^{3}}\phantom{\rule{4pt}{0ex}},\end{array}$$(12)
where we have used the relation r_{i + 1/2} = r_{i}(1 + δ)/2 = r_{i + 1}(1 + δ^{−1})/2. We note that Eq. (12) defines n(r) for any r_{i + 1/2}, but not for r ≠ r_{i + 1/2}. We choose to take the continuous limit of n(r_{i + 1/2}) for an arbitrary radius, replacing r_{i + 1/2} for an arbitrary r_{p},
$$\begin{array}{c}\hfill n({r}_{p})=\frac{3{N}_{p}}{32\pi}\frac{{(\delta +1)}^{3}}{{\delta}^{3}1}\frac{1}{{r}_{p}^{3}}\phantom{\rule{4pt}{0ex}}.\end{array}$$(13)
Eq. (13) ensures that n(r) satisfies exactly Eq. (11) when r_{p} = r_{i + 1/2} and that the particle shape is a smooth function of r_{p}. The particle width Δr(r_{p}) is determined in a similar manner; first, we express the grid spacing in terms of r_{i + 1/2}, Δr_{i} = r_{i + 1} − r_{i} = 2r_{i + 1/2}(δ − 1)/(δ + 1), and we extend this definition to an arbitrary radius r_{p},
$$\begin{array}{c}\hfill \mathrm{\Delta}r({r}_{p})=2{r}_{p}\frac{\delta 1}{\delta +1}\phantom{\rule{4pt}{0ex}}.\end{array}$$(14)
This quantity is represented for a typical grid in Fig. 4a, together with the grid spacing Δr_{i}. As expected, both quantities match exactly when r = r_{i + 1/2}, and Δr is a smooth function of r. Equations (13) and (14) ensure that the conservation law expressed in Eq. (11) can be extended to any radius, which is shown in Fig. 4b.
Fig. 4. Particle shape properties. Panel (a) shows the radial width and panel (b) the density and real particle number. 
The general particle shape S can be inferred from this discussion, and in particular from Eq. (13). It reads
$$\begin{array}{cc}\hfill S(r,\theta ,{r}_{p},{\theta}_{p})& =\frac{3}{16\pi}\frac{{(\delta +1)}^{3}}{{\delta}^{3}1}\frac{1}{{r}_{p}^{3}}{b}_{0}\left(\frac{r{r}_{p}}{\mathrm{\Delta}r({r}_{p})}\right)\times \hfill \\ \hfill & \times \frac{1}{cos({\theta}_{p}\mathrm{\Delta}\theta /2)cos({\theta}_{p}+\mathrm{\Delta}\theta /2)}{b}_{0}\left(\frac{\theta {\theta}_{p}}{\mathrm{\Delta}\theta}\right)\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(15)
where b_{0}(x) is the zeroth order bspline function, defined as b_{0}(x) = 1 if x< 0.5 and 0 otherwise. We note that Eq. (15) generalizes the particle shape to a twodimensional (r, θ) grid, hence the cos(θ_{p} ± Δθ) terms resulting from the integral in Eq. (11). With the shape function in Eq. (15), we can compute the charge density at any point (r, θ) due to the presence of a macroparticle with N_{p} real particles of charge q_{p} and coordinates (r_{p}, θ_{p}) as ρ_{p}(r, θ, r_{p}, θ_{p}) = q_{p}N_{p}S(r, θ, r_{p}, θ_{p}). The charge density at cell edges is defined resorting to the area/volume weighting technique described in Sect. 2.3, and can be formally derived as
$$\begin{array}{cc}\hfill {\rho}_{(i,j)}({r}_{p},{\theta}_{p})& =\frac{{\int}_{{V}_{i,j}}{\rho}_{s}(r,\theta ,{r}_{p},{\theta}_{p})\phantom{\rule{4pt}{0ex}}\mathrm{d}{V}_{i,j}}{{V}_{i,j}}\hfill \\ \hfill & =\frac{{q}_{p}{N}_{p}\frac{3}{16\pi}\frac{{(\delta +1)}^{3}}{{\delta}^{3}1}}{({r}_{i+1/2}^{3}{r}_{i1/2}^{3})(cos{\theta}_{j1/2}cos{\theta}_{j+1/2})}\times \hfill \\ \hfill & \times [\frac{{r}_{>}^{3}{r}_{<}^{3}}{{r}_{p}^{3}}][\frac{cos({\theta}_{p}\mathrm{\Delta}\theta /2)cos{\theta}_{j+1/2}}{cos({\theta}_{p}\mathrm{\Delta}\theta /2)cos({\theta}_{p}+\mathrm{\Delta}\theta /2)}]\phantom{\rule{4pt}{0ex}}.\hfill \end{array}$$(16)
We note that the special integration limits r_{>} = min(r_{p} + Δr(r_{p})/2, r_{i + 1/2}) and r_{<} = max(r_{p} − Δr(r_{p})/2, r_{i − 1/2}) result from the subtlety that the particle radial width is a function of the particle radial coordinate, r_{p}. The expressions in square brackets are often referred to as the weighting functions in PIC current deposition algorithms. We note also that Eq. (16) and the poloidal weighting functions defined in Eq. (6) constitute the same operation, with the volume integral centered in different grid positions, (r_{p}, θ_{p}) and (r_{i}, θ_{j}), respectively. Excluding the volume terms of the particle, given in Eq. (15), this point can be made explicit through
$$\begin{array}{c}\hfill {f}_{\mathit{ri}}=\frac{{r}_{p}^{3}{r}_{i}^{3}}{{r}_{i+1}^{3}{r}_{i}^{3}}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\mathrm{and}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\phantom{\rule{0.277778em}{0ex}}{\rho}_{(i,j)}({r}_{p},{\theta}_{p})\propto \frac{{r}_{>}^{3}{r}_{<}^{3}}{{r}_{i+1/2}^{3}{r}_{i1/2}^{3}},\end{array}$$(17)
for the radial components, which shows that current deposition and field interpolation schemes have the same order.
The particle shape in Eq. (15) and the deposition rule in Eq. (16) are the key ingredients in our chargeconserving current deposition scheme. This scheme is inspired by the seminal work of Villasenor & Buneman (1992) (hereafter VB92), which presented a scheme that predecessed the widely used method of Esirkepov (2001) for PIC current deposition in Cartesian grids. The VB92 method is schematically represented in Fig. 5a. Villasenor & Buneman (1992) proposed that the current density j should be computed directly by inverting the continuity equation, thus enforcing by construction that it is satisfied. In practice, when a particle is pushed in time from a position x^{n} to a position x^{n + 1}, part of its shape crosses the boundaries over which the current density is defined in the Cartesian PIC grid. These boundaries, and the exact locations where each of the components of j are defined, are shown in Fig. 5a in green and red lines and arrows, respectively. VB92 recognized that we can simply compute the different current density components by evaluating the fraction of charge density carried by each macroparticle that crosses the boundaries identified in green and red. For a Cartesian grid, this fraction can be computed geometrically as the ratio between the areas A_{green} and A_{red} and the total area corresponding to the particle shape, A_{total}. This calculation is simple in Cartesian grids because the particle shape does not change across the grid, which allows us to label which parts of the colored area at x > x_{i + 1/2} and y > y_{j + 1/2} crossed each of the green or red lines. In a spherical grid, this condition is not met, and the calculation becomes more involved.
Fig. 5. Schematic representation of the currentdeposition algorithm in (a) Cartesian and (b) spherical coordinates (see text for details). 
A schematic representation of the method equivalent to VB92 in a spherical grid is shown in Fig. 5b, where same rationale described above is easily applied except for the determination of the area identified with A_{?}. Because the particle expands during its motion from x^{n} to x^{n + 1}, it is not trivial to determine which fraction of A_{?} should be combined with A_{green} (A_{red}) to compute j_{r}_{(i + 1/2, j)} (j_{θ}_{(i, j + 1/2)}). We circumvent this issue by generalizing the geometrical interpretation of ∇ ⋅ j proposed by VB92. They suggested that the total current divergence can be split as ∇ ⋅ j = (∇⋅j)_{x} + (∇⋅j)_{y} in a Cartesian grid, with (∇⋅j)_{y} ∝ A_{green}/A_{total} and (∇⋅j)_{x} ∝ A_{red}/A_{total}, and that these terms could be computed directly by evaluating −∂ρ_{(i, j)}/∂t assuming that the particle moves purely along the corresponding direction at an average position along the orthogonal direction. Formally, this is expressed as
$$\begin{array}{cc}\hfill {{(\mathrm{\nabla}\xb7\mathbf{j})}_{x}}_{(i,j)}& ={\frac{\partial {\rho}_{(i,j)}}{\partial t}}_{{x}^{n},\overline{y}}^{{x}^{n+1},\overline{y}}=\frac{{\rho}_{(i,j)}({x}^{n+1},\overline{y}){\rho}_{(i,j)}({x}^{n},\overline{y})}{\mathrm{\Delta}t}\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(18)
$$\begin{array}{cc}\hfill {{(\mathrm{\nabla}\xb7\mathbf{j})}_{y}}_{(i,j)}& ={\frac{\partial {\rho}_{(i,j)}}{\partial t}}_{\overline{x},{y}^{n}}^{\overline{x},{y}^{n+1}}=\frac{{\rho}_{(i,j)}(\overline{x},{y}^{n+1}){\rho}_{(i,j)}(\overline{x},{y}^{n})}{\mathrm{\Delta}t}\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(19)
where $\overline{x}=({x}^{n+1}+{x}^{n})/2$ and $\overline{\mathit{y}}=({\mathit{y}}^{n+1}+{\mathit{y}}^{n})/2$. From Eqs. (18) and (19), we can express the divergence operators using finite differences and obtain j_{x}_{(i + 1/2, j)} and j_{y}_{(i, j + 1/2)}. This approach can be generalized to spherical coordinates, as we can write ∇ ⋅ j = (∇⋅j)_{r} + (∇⋅j)_{θ}. However, because the particle shape changes continuously in the radial direction, (∇⋅j)_{θ} cannot be computed assuming that the particle moves purely along the polar direction with $\overline{r}=({r}^{n+1}+{r}^{n})/2$. Instead, we proceed as follows: first, we compute ∇ ⋅ j and (∇⋅j)_{r} using
$$\begin{array}{cc}\hfill {(\mathrm{\nabla}\xb7\mathbf{j})}_{(i,j)}& ={\frac{\partial {\rho}_{(i,j)}}{\partial t}}_{{r}^{n},{\theta}^{n}}^{{r}^{n+1},{\theta}^{n+1}}=\frac{{\rho}_{(i,j)}({r}^{n+1},{\theta}^{n+1}){\rho}_{(i,j)}({r}^{n},{\theta}^{n})}{\mathrm{\Delta}t}\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(20)
$$\begin{array}{cc}\hfill {{(\mathrm{\nabla}\xb7\mathbf{j})}_{r}}_{(i,j)}& ={\frac{\partial {\rho}_{(i,j)}}{\partial t}}_{{r}^{n},\overline{\theta}}^{{r}^{n+1},\overline{\theta}}=\frac{{\rho}_{(i,j)}({r}^{n+1},\overline{\theta}){\rho}_{(i,j)}({r}^{n},\overline{\theta})}{\mathrm{\Delta}t}\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(21)
where $\overline{\theta}=({\theta}^{n+1}+{\theta}^{n})/2$. Then, we compute (∇⋅j)_{θ} = ∇ ⋅ j − (∇⋅j)_{r}. Finally, we invert the nabla operators,
$$\begin{array}{cc}\hfill {{(\mathrm{\nabla}\xb7\mathbf{j})}_{r}}_{(i,j)}& =3[\frac{{r}_{i+1/2}^{2}{{j}_{r}}_{(i+1/2,j)}{r}_{i1/2}^{2}{{j}_{r}}_{(i1/2,j)}}{{r}_{i+1/2}^{3}{r}_{i1/2}^{3}}]\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(22)
$$\begin{array}{cc}\hfill {{(\mathrm{\nabla}\xb7\mathbf{j})}_{\theta}}_{(i,j)}& =\frac{3}{2}\frac{{r}_{i+1/2}^{2}{r}_{i1/2}^{2}}{{r}_{i+1/2}^{3}{r}_{i1/2}^{3}}\times \hfill \\ \hfill & \times [\frac{sin{\theta}_{j+1/2}{{j}_{\theta}}_{(i,j+1/2)}sin{\theta}_{j1/2}{{j}_{\theta}}_{(i,j1/2)}}{cos{\theta}_{j1/2}cos{\theta}_{j+1/2}}]\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(23)
to find the current components. The inversion of (∇⋅j)_{θ}_{(i, j)} is simple, because the second term in the square brackets of Eq. (23) is always zero given that the particle motion is restricted to cell (i, j). The same is applicable to the inversion of (∇⋅j)_{r}_{(i, j)} for most particle positions in cell (i, j); however, due to the fact that the particle expands with r_{p}, it can deposit current at the grid position (i − 1/2, j) when r_{p} is close to r_{i}. When this happens, we determine (∇⋅j)_{r}_{(i − 1, j)} using Eq. (21), invert the corresponding operator to obtain j_{r}_{(i − 1/2, j)} and use it to solve for j_{r}_{(i + 1/2, j)} in Eq. (22). When particles cross two cells from x^{n} to x^{n + 1}, we split their trajectory such that each split is within a single cell, and apply the method described before to each trajectory split. The same strategy is applied in the algorithms proposed in Villasenor & Buneman (1992) and Esirkepov (2001). This method does not impose any restriction on the azimuthal current component, which we take to be simply j_{ϕ}_{(i, j)} = ρ_{(i, j)}v_{ϕ}, where v_{ϕ} is the macroparticle velocity in the azimuthal direction. Also, we adapt the integration regions as described in Sect. 2.2 to correct the evaluation of Eqs. (16) and (23) on the poloidal axes.
Finally, we note that Eqs. (16) and (22)(23) can also be derived by applying the algorithms in Villasenor & Buneman (1992) or Esirkepov (2001) (in firstorder) in a Cartesian logical space with the spherical coordinates metric; that is, applying a coordinate transformation to $(\stackrel{\sim}{r},\theta )\equiv (logr,\theta )$. In this Cartesianlike logical space, the grid becomes uniform and the cell center is located at r_{i + 1/2} = exp((r_{i} + r_{i + 1})/2), thus rendering the particle shape function S asymmetric in polar spherical coordinates. These special coordinates are less intuitive from a physical perspective but solve the particle shrinking/expansion problem described in this section, not requiring the special integration limit verification. A similar approach using logical coordinates was taken in (Belyaev 2015b). The method introduced in this paper is generic and conserves charge to machine precision even when adopting nonuniform grids, as we shall now demonstrate.
We benchmarked the current deposition method presented here by initializing particles all over the simulation domain with a random velocity (fixed number of particles per cell), depositing their current over a time step Δt and evaluating
$$\begin{array}{cc}\hfill {\mathrm{\Delta}}_{\mathrm{Continuity}}& =\frac{\mathrm{\Delta}t}{{\rho}_{(i,j)}}(\frac{\partial {\rho}_{(i,j)}}{\partial t}+{(\mathrm{\nabla}\xb7\mathbf{j})}_{(i,j)})\phantom{\rule{4pt}{0ex}},\hfill \end{array}$$(24)
$$\begin{array}{cc}\hfill {\mathrm{\Delta}}_{\mathrm{Gauss}}& =\frac{1}{{\rho}_{(i,j)}}({(\mathrm{\nabla}\xb7\mathbf{E})}_{(i,j)}4\pi {\rho}_{(i,j)})\phantom{\rule{4pt}{0ex}}.\hfill \end{array}$$(25)
Both Δ_{Continuity} and Δ_{Gauss} should be zero if the continuity equation and Gauss’ law are satisfied. Figure 6 shows that these quantities are both of the order of machine precision, 10^{−15} − 10^{−11}. The value of both Δ_{Continuity} and Δ_{Gauss} tends to be larger when closer to the star. Particles crossing cell boundaries lead to the accumulation of roundoff errors. Closer to the star the resolution is higher, and thus there are more crossing operations, which justifies the radial error dependence. The shrinkage and expansion of particles due to the use of a nonuniform grid with symmetric particle shape function S aggravates this error, as one single particle may deposit charge and current in three consecutive cells (in a single direction). The choice of logical coordinates reduces this number to two consecutive cells, thus performing better. Nevertheless, we have verified that the accuracy of the method is maintained over multiple time steps for a nonuniform grid by ensuring that the evolution of the grid integrals of Δ_{Continuity} and Δ_{Gauss} remain at machine precision level, for both single and multiparticle tests.
Fig. 6. Current deposition benchmarks, showing that both (a) the continuity equation and (b) Gauss’ law are satisfied to machine precision. 
This current deposition method thus accurately conserves charge, avoiding the need for other correcting algorithms. It is also inexpensive, since most factors in Eqs. (20)–(23) can be precomputed and reused throughout a simulation. Also, this method can be generalized for curved spacetime configurations (Torres et al. 2024).
2.5. Typical scales and normalizations
In the benchmarks presented above, the normalization units of distances, times, and fields varied according to what best suits the respective tests. However, for pulsar magnetosphere simulations, we adopt a common normalization that we introduce here. We choose to normalize distances to the stellar radius r_{*} and times to r_{*}/c. Electric and magnetic fields are normalized to m_{e}c^{2}/er_{*}, however we typically represent them in units of en_{GJ}r_{*}, where n_{GJ} = ΩB_{*}/2πec is the surface GoldreichJulian (GJ) (Goldreich & Julian 1969) particle number density. The GJ density also defines a typical frequency ${\omega}_{p,\mathrm{GJ}}=\sqrt{4\pi {e}^{2}{n}_{\mathrm{GJ}}/{m}_{e}}$ and an electron skin depth d_{e, GJ} = c/ω_{p, GJ}. The time step and grid spacing are chosen to resolve these temporal and spatial scales, respectively.
In pulsar magnetosphere simulations, the main parameter responsible for setting the typical temporal, spatial and energy scales is the normalized value of the surface magnetic field, B_{*}(er_{*}/m_{e}c^{2}). For realistic parameters, B_{*} ≃ 10^{12} G and r_{*} ≃ 10 km, we have B_{*}(er_{*}/m_{e}c^{2})∼10^{15}. Global simulations are not feasible with such values, since they would have to resolve scales of the order of ∼tens of r_{*} down to d_{e, GJ} ∼ 10^{−7}r_{*}. For this reason, we use more modest values of B_{*}(er_{*}/m_{e}c^{2})∼10^{3} − 10^{6}, such that we respect the ordering in these objects, Ω ≪ ω_{p, GJ} ≪ ω_{c}, where ω_{c} = eB_{*}/m_{e}c is the cyclotron frequency associated with a field magnitude B_{*}.
3. Global simulations of pulsar magnetospheres
In this Section, we present global PIC simulations of pulsar magnetospheres obtained with the OSIRIS framework (Fonseca et al. 2002, 2008). We start by allowing electronpositron pairs to be artificially and abundantly injected in our simulations, and then make increasingly realistic assumptions about the plasma supply processes, in particular regarding the regions of space where pair cascades operate, and the separation between kinetic and system scales.
All simulations presented here have a similar initial configuration: the system starts in vacuum and with an initial dipolar magnetic field of polar surface magnitude B_{*}; that is, with components B_{r}(r, θ) = B_{*}(r_{*}/r)^{3}cosθ and B_{θ}(r, θ) = (1/2)B_{*}(r_{*}/r)^{3}sin θ. The inner radial boundary is treated as a rotating conductor of angular velocity $\mathbf{\Omega}=\mathrm{\Omega}\widehat{\mathbf{z}}$; at the surface of the neutron star, we impose the corotation electric field E = −(v_{rot} × B)/c, with ${\mathbf{v}}_{\mathrm{rot}}=\mathbf{\Omega}\times ({r}_{\ast}\widehat{\mathbf{r}})$. In all simulations, we consider the stellar rotation frequency to be initially zero and increase it linearly over a time t_{rise}c/r_{*} = 1.5 to Ωr_{*}/c = 0.125. For times t > t_{rise}, the stellar frequency is kept constant. The stellar period is T = 2π/Ω = 50 r_{*}/c and the lightcylinder radius is R_{LC}/r_{*} = 8. All simulations use also r_{min}/r_{*} = 1 and r_{max}/r_{*} = 20, such that the plasma dynamics can be captured up to r/R_{LC} > 2. The value of B_{*} is chosen to satisfy the ordering Ω ≪ ω_{p, GJ} ≪ ω_{c} described in Sect. 2.5 while maintaining simulations numerically feasible. This choice and others regarding for example grid resolution vary according to the injection scheme and parameter regime under study, and are detailed alongside the corresponding simulations.
3.1. Volume injection
In this section, we inject plasma everywhere in the simulation domain where the local electric field component parallel to the magnetic field satisfies the condition E_{∥}c/r_{*}ΩB_{*} > k_{lim}, where k_{lim} is a constant. Similar injection criteria have been used in Belyaev (2015a), whereas in Philippov & Spitkovsky (2014), Kalapotharakos et al. (2018), Brambilla et al. (2018) plasma is only injected if the local magnetization is also above a given threshold. Physically, this injection scheme is equivalent to assuming that electronpositron pair cascades may develop wherever E_{∥} is sufficiently large, as it neglects any role of the local magnetic field magnitude or curvature. Since all fields (and in particular E_{∥}) decay with r, the choice of k_{lim} can also be interpreted as a spatial limitation to the plasma supply: infinitely small values of k_{lim} allow plasma to be injected up to r ≫ r_{*}, whereas k_{lim} ∼ 1 restricts the plasma supply to radii r ∼ r_{*}. A macroelectronpositron pair with particle weight w_{vol} = k_{vol}V_{p}E_{∥}/er_{*} and carrying a number density n_{vol} = w_{vol}/V_{p}, with k_{vol} = 0.2 and V_{p} being the macroparticle’s volume, is injected at rest in each cell and time step in which the injection condition is met. The choice of k_{vol} is such that a few macroparticles are required to supply the charge density that screens E_{∥} and stops the injection. We can also interpret k_{vol} as a parameter proportional to the local GJ density, since E_{∥}/er_{*} ∼ n_{GJ}. In all the simulations presented in this section, B_{*}er_{*}/m_{e}c^{2} = 8 × 10^{3}, N_{r} × N_{θ} = 1000^{2} and Δtc/r_{*} = 10^{−3}. In these conditions, c/ω_{p, GJ}r_{*} ≃ 0.022, whereas the minimum grid spacing is min(Δr_{i})/r_{*} ≃ 0.003.
In Fig. 7, we present an overview of the quasisteadystate solution obtained with k_{lim} = 0.005. This solution is achieved after a time ∼25 r_{*}/c ∼ T/2^{1}. In the first half stellar period, the simulation undergoes a transient stage in which the vacuum corotation fields are established and plasma is created. The solution presented in Fig. 7 resembles the canonical forcefree regime of pulsar magnetospheres: the magnetosphere is divided in two regions permeated by closed and open magnetic field lines (shown in white/black solid lines in all panels), with the last closed field line crossing the equatorial plane at the lightcylinder radius (shown in a white/black dashed vertical line in all panels). The open and closed field line regions are respectively negatively and positively charged, even if electrons and positrons exist in both regions – see Fig. 7ac, showing the electron and positron number density and the total charge density, respectively. As shown in Fig. 7d, a negative radial current density j_{r} (blue) is conducted from the polar regions and along the open field lines, which is compensated by return current layers (red) established on the last closed field line. The return current layers are connected with each other at a distance r ≃ R_{LC} on the equatorial plane, where the poloidal magnetic field lines resemble a Y shape. A radial current density layer extends along the equatorial plane to large distances, supporting a strong gradient in the toroidal magnetic field component B_{ϕ}, illustrated in Fig. 7e. The poloidal magnetic field lines have also opposite polarity in opposite sides of this equatorial current layer, and reconnect sporadically, leading to the formation of outflowing plasmoids – see the large density structures at r/r_{*} ≃ 12 in Fig. 7ab. The plasma supply in this simulation is large enough such that E_{∥} is effectively screened in the whole simulation domain, as shown in Fig. 7f, and thus lies well within the assumptions of the forcefree regime for pulsar magnetospheres.
Fig. 7. Forcefree magnetosphere obtained with volume injection. Panels (af) show the electron and positron density, total charge density, radial current density, azimuthal magnetic field, and electric field component parallel to the local magnetic field, respectively. Quantities are multiplied by powers of r to enhance the features at large radii. Black and white solid lines represent magnetic field lines, and vertical dashed lines show the location of the light cylinder. 
The quasisteadystate shown in Fig. 7 is sustained via intermittent injection, mainly along the return current layers. In these regions, E_{∥} is less efficiently screened, leading to the injection of plasma which, in turn, screens the field as it flows along the return current layers. As we shall demonstrate, this intermittency has a period of ≃0.3 − 0.5 T, and it may play a significant role in the temporal evolution of the magnetospheric state. However, for k_{lim} = 0.005 the solution never deviates significantly from the forcefree regime.
In order to demonstrate how the magnetospheric solution changes with k_{lim}, in Fig. 8 we compare the total charge density of the solutions obtained with k_{lim} = {0.005, 0.01, 0.1}. We reiterate the fact that k_{lim} is the minimum value of E_{∥}c/r_{*}ΩB_{*} for which we inject plasma. It is clear that the forcefree regime is only observed for k_{lim} = 0.005. For k_{lim} = 0.01, the equatorial current sheet (positively charged region at r ≳ R_{LC}) is wide and the return current layers are not positively charged everywhere, and for k_{lim} = 0.1 the solution does not even produce an outflow. In fact, by increasing k_{lim}, we are limiting the plasma supply to regions closer and closer to the stellar surface. This can be understood by noting that this parameter compares the local E_{∥} with the reference value ΩB_{*}r_{*}/c (i.e., the surface magnitude of the electric field in vacuum). Since the typical magnitude of E_{∥} decreases with r, increasing k_{lim} limits plasma injection to smaller radii. In the k_{lim} = 0.01 run, this supply occurs only up to radii r/r_{*} ≃ 3, and the solution shows the same intermittency observed for k_{lim} = 0.005. However, the injection stage is not as efficient in this case, and the equatorial outflow is not dense enough to produce a thin current sheet. For k_{lim} = 0.1, only regions close to the surface can initially fulfil the injection criteria, and no plasma is supplied to large radii. The system relaxes in this case to a fully chargeseparated configuration, with only electrons (positrons) in the poles (equatorial region). This solution is often denominated as the diskdome or electrosphere solution (Jackson 1976; KrausePolstorff & Michel 1985). In the charged regions, the electric field is screened, injection ceases and no plasma outflows are formed.
Fig. 8. Magnetospheric solutions obtained with volume injection. The panels show the total charge density after a stellar rotation period. 
An important property of the magnetospheric solution is the integrated Poynting flux L(r), defined as
$$\begin{array}{c}\hfill L(r)=\frac{c}{2}{\displaystyle {\int}_{0}^{\pi}}{(\mathbf{E}\times \mathbf{B})}_{r}\phantom{\rule{4pt}{0ex}}{r}^{2}sin\theta \mathrm{d}\theta \phantom{\rule{4pt}{0ex}}.\end{array}$$(26)
Figure 9 shows L(r) as a function of time for the three simulations described before. This quantity is normalized to the theoretical value of the spindown luminosity, L_{0} = μ^{2}Ω^{4}/c^{3}, with μ = B_{*}r_{*}^{3}. We observe a large spindown at early times for all simulations, which is a consequence of the initial transient stage. After this transient, the k_{lim} = 0.1 simulation converges to a surface Poynting flux L_{*}/L_{0} ≪ 1, which is a consequence of the inactivity of diskdome solution. On the contrary, the simulations with lower k_{lim} have L_{*}/L_{0} ∼ 1. The Poynting flux remains approximately constant within the lightcylinder for these runs, and decays with r for r > R_{LC}, which is a signature of the conversion from magnetic to kinetic energy due to magnetic reconnection in the equatorial plane. The surface Poynting flux shows variations of periodicity 0.3 − 0.5 T, which are correlated with the intermittency of the solution identified above in this section. The timeaveraged radial dependence of the luminosity ⟨L⟩ after a stellar period and the temporal dependence of L_{*} is shown in Fig. 10.
Fig. 9. Poynting flux in simulations with volume injection. Values are normalized to the theoretical value L_{0} = μ^{2}Ω^{4}/c^{3}. 
Fig. 10. Radial and temporal dependencies of Poynting flux in simulations with volume injection. Panel (a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and panel (b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in (a) and (b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 
The simulations presented in this section show that the efficiency of the plasma supply critically determines the global structure of the pulsar magnetosphere. Independently of the injection criteria variants, the volumetric prescription of plasma allows for efficient screening of all magnetospheric gaps, thus approaching the forcefree solution everywhere, as documented in similar works (Philippov & Spitkovsky 2014; Belyaev 2015a; Kalapotharakos et al. 2018; Brambilla et al. 2018). In the present study, we show that k_{lim} effectively controls the radial extent of the plasma supply, demonstrating that quasiforcefree solutions can only be achieved when the outer gap is, at least, partially screened. This fact is physically very relevant, as it determines the amount of fully open magnetic field lines contributing towards the stellar spindown power, reflected in the observed luminosity and also in the poloidal extent of the polar cap. Here, we show that pulsars with quasiforcefree magnetospheres (e.g., k_{lim} = 0.005 simulation) must be observed with ⟨L⟩∼L_{0} and the radio beam poloidal extent, if resolved, must be close to the polar cap angle provided by forcefree models, ${\theta}_{\mathrm{PC}}\sim \sqrt{{r}_{\ast}/{R}_{\mathrm{LC}}}$. This is in contrast with subforcefree magnetospheres (e.g., k_{lim} = 0.01 simulation), resembling the weakpulsar solution (Gruzinov 2015), characterized by an outer gap, wider equatorial current sheet, and a thinner open flux tube from the polar cap due to the expansion of the volumetric return current (see Fig. 8). Consequently, such pulsars must be observed with smaller spindown luminosities, ⟨L⟩≲0.8L_{0}, and thinner radio beams, due to the compression of the polar cap, making them distinguishable from regular forcefree pulsars.
Although the volumetric plasma injection allows for a spatially regulated prescription of plasma, we point out that this plasma prescription method may lead to unphysical overinjection of plasma close to the rotation axis, as seen in Figs. 7 and 8, enabled by parallel electric field fluctuations. This effect is particularly noticeable for the k_{lim} = 0.005 simulation, which possesses the minimum value of E_{∥}c/r_{*}ΩB_{*}. Also, the efficiency of pair production must not be spatially uniform in the magnetosphere, as this model assumes. In particular, pair production at radii comparable to the lightcylinder radius is only possible for pulsars with sufficient optical depth for the activation of the γ − γ pairproduction channel (Cheng et al. 1986). It is then important to assess if more realistic injection and/or pairproduction schemes can provide the plasma supply required for the magnetosphere to be in the forcefree regime. In the next sections, we address this question by considering plasma supply schemes limited to regions close to the stellar surface.
3.2. Surface injection
In this section, we limit injection to occur only at the stellar surface. In doing so, we phenomenologically introduce the important role of the magnetic field amplitude in our treatment of the magnetospheric plasma supply. As in Sect. 3.1, we do not allow particles to emit photons and/or pairs. We adopt two different criteria for the injection and vary the density and velocity of the surfaceinjected plasma. The parametrization of the plasma flow injected from the stellar surface is similar to that presented in Cerutti et al. (2015). However, our criteria for injection differ slightly from that work, which also assumes a minimum threshold for the local plasma magnetization. In all simulations presented in this section, we use B_{*}er_{*}/m_{e}c^{2} = 8 × 10^{3}, N_{r} × N_{θ} = 500^{2} and Δtc/r_{*} = 3 × 10^{−3}.
The first injection criterion is based on the local value of E_{∥}. We inject a macroelectronpositron pair in each cell just above the stellar surface (r = r_{*}) that satisfies E_{∥}c/r_{*}ΩB_{*} > k_{lim}. In this case, we consider a fixed k_{lim} = 0.002 and vary the properties of the injected pairs, namely their density n_{s} = w_{s}/V_{p} = k_{s}n_{GJ} and poloidal velocity v_{s}. These pairs are also injected with a toroidal velocity that matches the local linear velocity of the stellar surface, v_{ϕ} = Ωrsin θ.
Despite the large range of injection parameters considered, k_{s} = n_{s}/n_{GJ} = {0.2, 0.5, 1} and v_{s}/c = {0, 0.1, 0.5, 0.99}, the solutions obtained for long times, t/T ≳ 2, always converge to the diskdome solution identified in Sect. 3.1. Figure 11 shows the charge density ρ and E_{∥} of two runs with k_{s} = n_{s}/n_{GJ} = 1 and v_{s} = {0, 0.99} after a time t/T ≃ 4. After an initial transient, the system settles to a chargeseparated solution and effectively screens E_{∥} at the stellar surface, precluding further injection.
Fig. 11. Magnetospheric solutions obtained with surface injection proportional to E_{∥}. Panels (a12) show the total charge density and E_{∥}, respectively for a simulation with v_{s}/c = 0 and panels (b12) show the same but for a simulation with v_{s}/c = 0.99. Solid lines represent magnetic field lines, and vertical dashed lines show the location of the light cylinder. 
The second injection criterion does not depend on the local surface field conditions. Instead, injection is allowed in all cells above the stellar surface in which the combined local number density of positrons and electrons satisfies n_{+} + n_{−} < 5 n_{GJ}, to ensure that enough plasma exists everywhere to screen the local electric field parallel to the magnetic field. We emphasize that n_{GJ} = ΩB_{*}/2πec is the pole GJ density and not its local value. This criterion allows injection to occur even if E_{∥} ∼ 0, and is thus harder to motivate from firstprinciples arguments. Here, we shall interpret it as a means of producing a set plasma density over a layer near the stellar surface of width smaller than the local resolution of the simulation grid. In pulsars, such layer can be as small as ∼100 m (Ruderman & Sutherland 1975). We consider that the injected electronpositron pairs carry a number density n_{s} = k_{s}n_{GJ} and poloidal velocity v_{s}.
In Fig. 12, we show the charge density distribution of the solutions obtained for a fixed k_{s} = n_{s}/n_{GJ} = 0.2 and varying v_{s} for a time t/T = 1. With v_{s} = 0, the system converges to the electrosphere solution. Particles injected at early times develop a spacecharge limited flow, driving E_{∥} to zero near the stellar surface and thus inhibiting freshly injected particles to be pulled away from or towards the star. For v_{s} > 0, we observe that the system develops a positively charged outflow along the equatorial plane. This outflow occurs in a narrower current sheet for larger values of v_{s}, which can be understood as a mechanism to support the stronger toroidal magnetic field driven by the stronger poloidal currents of these regimes. However, we do not observe a current sheet as thin as that characteristic of the forcefree regime. Instead, the current sheet remains wide even for v_{s}/c = 0.99. This may indicate that the plasma launched into this region is not dense enough, a question that we address below in this section.
Fig. 12. Magnetospheric solutions obtained with surface injection proportional to n_{GJ} with fixed k_{s} = n_{s}/n_{GJ} = 0.2 and varying v_{s}/c. 
Figure 13 shows the timeaveraged Poynting flux produced by the simulations described above with surface injection as a function of the radial coordinate r and its surface value as a function of time. We see once again that an electrosphere solution (v_{s}/c = 0) produces no spindown luminosity, and that it increases overall with increasing v_{s}. The same decrease for r > R_{LC} observed in Sect. 3.1 is observed here. We note that the v_{s}/c = 0.99 run shows a surface Poynting flux larger than L_{0}, which is a consequence of the smaller size of the corotation region (and thus a smaller effective lightcylinder radius and larger effective L_{0}). This particular simulation demonstrates that the overinjection drives the system to expand the effective polar cap area, supporting an artificially stronger than forcefree magnetospheric current, which agrees with the discussion on the polar cap adaptation to the global magnetospheric current in Sect. 3.1.
Fig. 13. Radial and temporal dependencies of Poynting flux in simulations with surface injection proportional to n_{GJ} with fixed k_{s} = n_{s}/n_{GJ} = 0.2 and varying v_{s}/c. a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in a) and b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 
We also performed a set of simulations with fixed v_{s}/c = 0.5 and varying k_{s} = n_{s}/n_{GJ} = {0.1, 0.2, 0.5}. The charge density obtained in the steadystate (or quasisteadystate) of these simulations is shown in Fig. 14. These results confirm that the denser the injected plasma is, the more the solution approaches the forcefree regime (see in particular the solution obtained for k_{s} = 0.5). This injection density requirement seems to be critical in the launching of highdensity plasma to large radii, in particular along the return current layers that connect the surface to the equatorial current sheet. This result highlights the critical role of populating/screening the outer gap in promoting the steadystate magnetospheric solution from chargestarved to forcefree, having implications on the observed spindown luminosity and radio beam width.
Fig. 14. Magnetospheric solutions obtained with surface injection proportional to n_{GJ} with fixed v_{s} and varying k_{s} = n_{s}/n_{GJ}. 
Consistently with simulations presented in Sect. 3.1, intermediate supply of plasma leads to a subforcefree magnetospheric solution with a thick current sheet configuration (see Figs. 8b, 12bd and 14ab). At the location of the lightcylinder radius, we estimate the most energetic particles to possess a gyroradius of about 2.8R_{*}, which is larger than the measured width of the current sheet. This difference may indicate that the most energetic particles are lost and the current sheet is mainly supported electrostatically by a colder population, similarly to the solution detailed in Contopoulos et al. (2014). The inclusion of radiation reaction cooling significantly reduces the gyration of these particles, which leads to a thinner current sheet configuration (Uzdensky & McKinney 2011; Philippov et al. 2015a; Hakobyan et al. 2019, 2023; Schoeffler et al. 2023). This effect on the global magnetospheric structure and current sheet configuration is still overlooked, particularly for weak pulsar magnetospheres.
In summary, some of the parameters used in simulations presented in this section yield active magnetospheric solutions, with L_{*}/L_{0} ∼ 1 and a global configuration similar to the forcefree regime. Weak partially chargestarved pulsar magnetospheric solutions with persistent outer gaps and lower timeaveraged luminosities were also presented. Both regimes are consistent with the results presented in Cerutti et al. (2015). However, it is hard to motivate the injection criteria and the choice of numerical parameters required to observe such regimes.
3.3. Pair production
The results presented in Sects. 3.1 and 3.2 are in good agreement with similar previous works. In particular, both Philippov & Spitkovsky (2014) and Cerutti et al. (2015) observe a transition from electrosphere to active solutions with more abundant plasma supply. While in Philippov & Spitkovsky (2014) pairs are injected up to large radii, in Cerutti et al. (2015) only surface injection is considered, showing trends with k_{s} and v_{s} very similar to our results.
The convergence to a forcefree regime in the asymptotic limit of large plasma supply with both volume and surface injection is reassuring. However, an important question remains open when translating global simulations with volume and surface injection schemes to realistic systems: how is this plasma supplied, if strong field pair production operates efficiently only near the stellar surface? Is this pairproduction channel enough to supply the plasma to fill the whole magnetosphere?
In young and rapidly rotating pulsars (e.g., the Crab pulsar and other gammaray pulsars), pairs can also be created via the γγ channel (Cheng et al. 1986). In this process, for which the crosssection peaks at around a center of mass energy ∼2 m_{e}c^{2}, gammarays produced via synchrotron emission and/or inverse Compton scattering in the equatorial current sheet collide with photons from a low energy bath (e.g., supplied by synchrotron radiation from leptons in the weaker magnetic fields of the outer magnetosphere), producing pairs (Cheng et al. 1986). However, slower pulsars are not expected to have a sufficiently dense lowenergy photon bath for this process to be relevant, and strong field pair production remains the main plasma supply channel (Sturrock 1971).
In this section, we use global simulations that include pair production only near the stellar surface to understand whether it can provide enough plasma to maintain an active magnetospheric solution. We use the heuristic pairproduction model described in Cruz et al. (2021b, 2022), in which a lepton emits a pair of combined energy γ_{pair}m_{e}c^{2} whenever it achieves a threshold Lorentz factor γ_{thr}. We keep the ratio γ_{thr}/γ_{pair} constant, and vary the ratio η ≡ γ_{max}/γ_{thr}, where γ_{max} = eΦ_{pc}/m_{e}c^{2} is the maximum energy achievable by the particles in the voltage Φ_{pc} = B_{*}r_{*}^{3}Ω^{2}/c^{2} induced by the rotating star across the polar cap. In general, γ_{pair} ≪ γ_{thr} ≪ γ_{max} in real systems; however, it is very hard to achieve a large separation between these scales in global PIC simulations. For instance, previous works, considering a similar pairproduction model (Chen 2017; Philippov et al. 2015a), have used η ∼ 10 and γ_{thr}/γ_{pair} ∼ 2, which severely limits the efficiency of the pair cascades and the plasma multiplicity. In this Section, we present simulations with fixed γ_{pair} = 16 and γ_{thr} = 25 and a range of large values of η. We achieve this by controlling the surface magnetic field amplitude B_{*}. In doing this, besides increasing the scale separation between pair production and the dynamical scales, we also decrease the plasma kinetic scales. For this reason, we adopt a varying number of grid cells and time steps in our simulations to be able to resolve these scales. For η = 5 we use N_{r} × N_{θ} = 500^{2} and Δtc/r_{*} = 3 × 10^{−3}, for η = {25, 50} we use N_{r} × N_{θ} = 1000^{2} and Δtc/r_{*} = 10^{−3} and for η = {100, 150} we use N_{r} × N_{θ} = 2000^{2} and Δtc/r_{*} = 5 × 10^{−4}.
In order to mimic the relevance of the large magnetic field required for pair production to occur, we limit pair production to only occur at radii r/r_{*} < 3. We also forbid pair production for θ < 0.01, to reproduce the suppression of the corresponding QED crosssection in this region (Cruz et al. 2021a). Seed electronpositron pairs are provided at the stellar surface whenever E_{∥}c/r_{*}ΩB_{*} > k_{lim}, with k_{lim} = 0.1. Each pair is injected at rest and carrying a density n_{s} = k_{s}E_{∥}/er_{*}, with k_{s} = n_{s}/n_{GJ} = 0.2. We stress that in these conditions, we obtained an electrosphere configuration in simulations without pair production (see section 3.2).
In Figure 15, we show the charge density obtained at a time t/T ≃ 2 for a relevant subset of the simulations performed. We observe a transition from electrosphere to forcefreelike configurations by increasing η. Physically, this corresponds to allowing the creation of more pairs per particle, increasing the plasma supply of the system. For η = 5, pair production is not efficient enough, and after an initial transient with some pair production, the accelerating electric field is screened and the system settles to an inactive solution. For η ∼ 10 − 50, the system is able to launch plasma towards the lightcylinder and produce a positively charged equatorial outflow. This plasma is launched along the return current layers due to pair production at r/r_{*} < 3; however, because of the limited effectiveness of the pair production in this range of η, the plasma produced is not dense enough to screen the outer gap and to confine the equatorial current sheet to a thin region, and it becomes wide for large distances from the stellar surface. For η ≳ 100, the system converges to a near forcefree regime, with magnetic field lines open to infinity and a thin equatorial current sheet. In these simulations, pair production is very effective, and launches a large density (n∼ few n_{GJ}), quasineutral plasma to the lightcylinder. In this region, part of the plasma escapes along the equatorial field lines; however, a fraction of the particles flows back to the star. The majority of these particles are electrons, such that the return current layers are negatively charged.
Fig. 15. Magnetospheric solutions obtained with pair production. Panels ad) show the total charge density for simulations with η = {5, 25, 50, 150}. Solid lines represent magnetic field lines, and vertical dashed lines show the location of the lightcylinder. 
The timeaveraged radial dependence of the Poynting flux and its surface value as a function of time for the simulations described above are presented in Figure 16. The observed radial dependence is similar to the regimes previously observed, with the η ≳ 100 simulations approaching the forcefree spindown luminosity L_{0} within the lightcylinder. In the equatorial current sheet, a fraction of 0.3 − 0.4 L_{*} is dissipated between r ∼ R_{LC} and r ∼ 2 R_{LC} and converted into particle kinetic energy. For all η < 100 runs, the surface luminosity decreases over time, and we expect them to eventually converge to the electrosphere solution for t/T ≫ 1. However, for η ≳ 100, the surface Poynting flux remains stable over time.
Fig. 16. Radial and temporal dependencies of Poynting flux in simulations with pair production with varying η. a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in a) and b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 
All simulations present some temporal variability. We see smallscale fluctuations on the charge and current densities in the open field line outflows, due to the E_{∥} screening process resulting from pair cascades. These fluctuations occur on a temporal scale ∼r_{*}/ηc. We also observe a quasiperiodic launch of plasma towards the lightcylinder region along the return current layers with a temporal scale ∼0.3 − 0.5 T. We show one of these events in Figure 17 for a simulation with η = 100. As plasma is injected along the last closed field lines, most of it escapes along the equatorial current sheet. As this happens, the return current density drops close to r ∼ R_{LC}, allowing E_{∥} to grow. Electrons flowing back to the star are thus accelerated along these field lines and produce a large number of pairs when they enter the pair producing region r/r_{*} < 3 – see for example Figure 17a1) and b1). The secondary particles then advect to large radii along the return current layers, reestablishing j_{r} and effectively screening the E_{∥} responsible for triggering the process – see Figure 17d13). This process produces a larger fraction of the total pair production events for 10 ≲ η ≲ 50. As in Sect. 3.1, the solutions obtained in this range resemble that of weak pulsars (Gruzinov 2015), with screened surface E_{∥} but unscreened outer gaps and wide equatorial current sheets as a result of inefficient pair production. The process presented here is similar to that described in Chen et al. (2020), Bransgrove et al. (2023).
Fig. 17. Cyclic pair production along the return current layers. Columns labeled 13 correspond to different times and rows labeled ad show the electron and positron densities, radial current and E_{∥}, respectively. Results obtained for η = 100. Solid lines represent magnetic field lines. 
The periodicity of the cyclic behavior driven by pair production along the return current layers is ∼0.3 − 0.5 T. We believe that this periodicity can depend on the multiplicity from the pair cascade near r/r_{*} ∼ 3, since if more pairs outflow during the active phase, more electrons can be stored in the Ypoint charge cloud, which takes longer to deplete. If this is true, a larger multiplicity should translate to a longer duty cycle. A detailed study of the importance of the cascade multiplicity on the cyclic behavior is deferred to a future work.
Finally, we note that apart from the effective pair discharges along the return current layers, we also observe abundant pair production within the polar cap region for all simulations with η > 5 – see Figure 18 for an illustrative example. This occurs because the density supplied from the stellar surface is insufficient to screen E_{∥} in this region. With stronger surface injection, we expect this pair production to be less significant. However, we do not expect the overall structure of the magnetosphere to be meaningfully modified. Interestingly, the polar cap pair production observed in this regime resembles that expected when general relativity effects are taken into account. When corrections due to the strong gravitational field of the neutron star are considered, we expect pair creation activity within the polar cap even if the surface can supply a charge density ±en_{GJ} (Philippov et al. 2015a; Chen et al. 2020; Bransgrove et al. 2023), since general relativity requires a current in this region j_{r}> en_{GJ} (Beloborodov 2008; Belyaev & Parfrey 2016; Gralla et al. 2016; Torres et al. 2024). This current starvation state was shown to promote efficient particle acceleration in the polar cap (Philippov et al. 2015a; Chen et al. 2020; Bransgrove et al. 2023; Torres et al. 2024), reducing the minimum η for which the temporal evolution of the surface Poynting flux in Fig. 16b is stable. Besides driving this difference in the timedependent nature of the polar cap, general relativity also reduces the polar cap poloidal extension (Belyaev & Parfrey 2016; Gralla et al. 2016), affecting the radio emission profile of more compact or rapidly rotating lowobliquity pulsars (Torres et al. 2024).
Fig. 18. Pairproduction sites for the same simulation and times shown in Figure 17. The color indicates the number of pairproduction events between data dumps (≃0.3 r_{*}/c) in each grid cell. 
4. Conclusions
In this work, we presented a systematic study of the different global regimes of pulsar magnetospheres. Namely, we performed simulations with three distinct plasma sources: in volume, from the stellar surface, and via pair production. Our results, presented in Sect. 3, show that all plasma sources produce near forcefree solutions in the regime of large plasma supply. In the opposite regime, we obtain inactive electrosphere solutions with all sources. These results are in overall good agreement with other works that independently consider volume (Philippov & Spitkovsky 2014; Belyaev 2015a; Kalapotharakos et al. 2018; Brambilla et al. 2018) or surfaceinjection schemes (Cerutti et al. 2015; Hakobyan et al. 2023), and those that use heuristic pairproduction models (Chen & Beloborodov 2014; Philippov et al. 2015b,a; Chen et al. 2020; Guépin et al. 2020; Bransgrove et al. 2023).
While volume and surfaceplasma injection serve as a means to efficiently fill the pulsar magnetosphere and produce a near forcefree configuration, as shown in Sects. 3.1 and 3.2, respectively, these are hard to motivate from firstprinciple arguments. On one hand, the pair cascades that these injection schemes aim to mimic develop only when the local magnetic field is close to the Schwinger field, and as such they should only operate near the stellar surface. On the other hand, these cascades produce plasma with a complex energy distribution that depends on the local electric and magnetic field geometry for example. Thus, any volume or surfaceinjection scheme is a substantial simplification of the highly nonlinear plasma supply from pair cascades in pulsars. Understanding whether or not and how pair production alone can fill the whole pulsar magnetosphere is thus crucial, in particular in order to reliably determine observational signatures.
The simulations including pair production presented in Sect. 3.3 show that pair discharges operating close to the stellar surface produce a range of solutions of the pulsar magnetosphere. The character of the solution depends critically on the ratio between the maximum attainable particle energy and the energy at which leptons emit pairproducing photons, η = γ_{max}/γ_{thr}, which quantifies the efficiency of the pair discharges. Our results show that when η ≳ 100, a sufficient number of pairs are created to fill the magnetosphere and reach a near forcefree surface Poynting flux, with dissipation occurring in an equatorial current sheet beyond the light cylinder. In the opposite limit, η ≲ 10, the magnetosphere settles to a fully chargeseparated, static solution, with E_{∥} = 0 near the surface, which produces a negligible Poynting flux. For η ∼ 10 − 50, we observe an intermediate solution (Gruzinov 2015), with a wide equatorial current sheet and with a surface Poynting flux that is 50 − 80% below that expected in the forcefree regime.
Our simulations show that the pair production along the return current layers is essential in order to feed plasma to the lightcylinder region and beyond in near forcefree regimes, in line with the results reported in other works, such as Chen & Beloborodov (2014). We also identify a timedependent mechanism similar to that presented in Chen et al. (2020), Bransgrove et al. (2023), which results from periodic openings of an outer gap in which particles flowing back to the star are able to accelerate, producing pairs when they get close to the stellar surface.
Although all simulations presented in this study use small values of R_{LC}/r_{*}, corresponding to millisecondtype pulsars, we believe that the presented results are applicable to a wider range of stellar rotation periods. In particular, we expect our simulations to model pulsars approaching their death line; that is, when the magnetosphere becomes fully chargestarved and pair production is no longer possible (disk–dome solution). From an observational standpoint, we may distinguish these solutions from the regular forcefree magnetospheres by exploring the compression of the openflux tube caused by the presence of a persistent outer gap and volumetric return currents, as seen for all plasmainjection methods. For more compact magnetospheres, capable of γ − γ pair production near the Ypoint (Cheng et al. 1986), these may lead to episodic switches of the whole magnetosphere between weak and near forcefree solutions, which could be observed as intermittent pulsars (Li et al. 2012). Following the discussion in Sect. 3.1, we hypothesize that the continuous adaptation of the polar cap poloidal extension to the fluctuating global magnetospheric current should also imprint similar spatiotemporal features in the observed radio beam properties (e.g., beam width).
The simulations presented here use a very simple heuristic model to describe pair production in strong magnetic fields. In this work, we only explored the influence of the parameter η on the magnetospheric structure and left the ratio γ_{thr}/γ_{pair} unchanged. This ratio plays an important role in the multiplicity of pair cascades, and was kept low to make simulations feasible. Larger values of γ_{thr}/γ_{pair} will likely provide even more abundant pairs to large radii, such that smaller values of η may be enough to set the magnetosphere in a forcefree regime. We defer further exploration of this topic to future work.
The pairproduction model considered here provides an adequate description of pair cascades when the curvature photon mean free path is negligible; that is, when pair production is local. In global models, however, it is easy to conceive that photons emitted in some regions of the magnetosphere may decay into pairs in others. For instance, photons emitted by electrons traveling toward the star along the return current layer may decay in the polar cap region. It would thus be interesting to include more sophisticated pairproduction models in these simulations to assess whether or not nonlocal pair production plays a significant role in, for example, coherent emission processes.
In this work, we also describe a spherical grid that can be used to perform global PIC simulations of pulsar magnetospheres. We provide details of (a) an electromagnetic field solver –based on the Yee solver– that uses an integral form of Maxwell’s equations (Sect. 2.2), (b) particle pushers that solve the particle equations of motion in Cartesian coordinates (Sect. 2.3), and (c) a chargeconserving current deposition scheme (Sect. 2.4) for a nonuniform, curvilinear spherical grid. While the field solver and particle pusher techniques are also implemented in other similar codes, the current deposition scheme presented here is a novel development. By ensuring that the continuity equation (and, consequently, Gauss’ law) is satisfied in the current deposition, this method does not require that other numerical algorithms be used to correct for artificial charges in the grid. For each of the numerical schemes presented here, we provide comprehensive benchmarks for a variety of test scenarios. All numerical schemes presented here have been implemented in the PIC code OSIRIS.
Acknowledgments
FC, TG, RT, RAF and LOS acknowledge support by the European Research Council (ERC2015AdG Grant 695088) and FCT (Portugal)  Foundation for Science and Technology (grants PD/BD/114307/2016 and PD/BD/142971/2018, in the framework of the Advanced Program in Plasma Science and Engineering APPLAuSE, grants PD/00505/2012 and PD/00505/2018, and the XMASER project no. 2022.02230.PTDC). AC acknowledges support from NSF grants DMS2235457 and AST2308111. AS is supported by the NSF grant PHY2206607 and the Simons Foundation grant MPSCMPS00001470. We acknowledge PRACE (PULSAR project) for granting access to MareNostrum, Barcelona Supercomputing Center (Spain), where the simulations presented in this work were performed.
References
 Beloborodov, A. M. 2008, ApJ, 683, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A. 2015a, MNRAS, 449, 2759 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A. 2015b, New Astron., 36, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A., & Parfrey, K. 2016, ApJ, 830, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (IoP) [CrossRef] [Google Scholar]
 Boris, J. 1970, Proceedings of the Fourth Conference on Numerical Simulation of Plasmas (Washington DC: Naval Research Laboratory), 3 [Google Scholar]
 Brambilla, G., Kalapotharakos, C., Timokhin, A. N., Harding, A. K., & Kazanas, D. 2018, ApJ, 858, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., Philippov, A. A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, A. Y. 2017, Ph.D. Thesis, Columbia University, USA [Google Scholar]
 Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, ApJ, 889, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, K. S., Ho, C., & Ruderman, M. 1986, ApJ, 300, 500 [CrossRef] [Google Scholar]
 Contopoulos, I., Kalapotharakos, C., & Kazanas, D. 2014, ApJ, 781, 46 [NASA ADS] [Google Scholar]
 Cruz, F., Grismayer, T., Chen, A. Y., Spitkovsky, A., & Silva, L. O. 2021a, ApJ, 919, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, F., Grismayer, T., & Silva, L. O. 2021b, ApJ, 908, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Cruz, F., Grismayer, T., Iteanu, S., Tortone, P., & Silva, L. O. 2022, Phys. Plasmas, 29, 052902 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, J. M. 1962, Phys. Fluids, 5, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, J. M. 1983, Rev. Mod. Phys., 55, 403 [CrossRef] [Google Scholar]
 Esirkepov, T. Z. 2001, Comput. Phys. Commun., 135, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Fonseca, R. A., Silva, L. O., Tsung, F. S., et al. 2002, in Computational Science  ICCS 2002, eds. P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, J. J. Dongarra, et al. (Berlin, Heidelberg: Springer, Berlin Heidelberg), 342 [CrossRef] [Google Scholar]
 Fonseca, R. A., Martins, S. F., Silva, L. O., et al. 2008, Plasma Phys. Control. Fusion, 50, 124034 [CrossRef] [Google Scholar]
 Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869 [Google Scholar]
 Gralla, S. E., Lupsasca, A., & Philippov, A. 2016, ApJ, 833, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Grismayer, T., Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Phys. Plasmas, 23, 056706 [NASA ADS] [CrossRef] [Google Scholar]
 Grismayer, T., Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2017, Phys. Rev. E, 95, 023210 [NASA ADS] [CrossRef] [Google Scholar]
 Gruzinov, A. 2015, arXiv eprints [arXiv:1503.05158] [Google Scholar]
 Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hakobyan, H., Philippov, A., & Spitkovsky, A. 2019, ApJ, 877, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, 943, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (CRC Press) [Google Scholar]
 Jackson, J. D. 1975, Classical Electrodynamics, 2nd edn. (New York, NY: Wiley) [Google Scholar]
 Jackson, E. A. 1976, ApJ, 206, 831 [NASA ADS] [CrossRef] [Google Scholar]
 Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [NASA ADS] [CrossRef] [Google Scholar]
 KrausePolstorff, J., & Michel, F. C. 1985, MNRAS, 213, 43P [Google Scholar]
 Li, J., Spitkovsky, A., & Tchekhovskoy, A. 2012, ApJ, 746, L24 [CrossRef] [Google Scholar]
 Mur, G. 1981, IEEE Trans. Electromagn. Compat., EMC23, 377 [CrossRef] [Google Scholar]
 Philippov, A. A., & Spitkovsky, A. 2014, ApJ, 785, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJ, 815, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015b, ApJ, 801, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51 [Google Scholar]
 Schoeffler, K. M., Grismayer, T., Uzdensky, D., & Silva, L. O. 2023, MNRAS, 523, 3812 [NASA ADS] [CrossRef] [Google Scholar]
 Sturrock, P. A. 1971, ApJ, 164, 529 [Google Scholar]
 Torres, R., Grismayer, T., Cruz, F., Fonseca, R., & Silva, L. 2024, New Astron., 112, 102261 [NASA ADS] [CrossRef] [Google Scholar]
 Touati, M., Codur, R., Tsung, F., et al. 2022, Plasma Phys. Control. Fusion, 64, 115014 [CrossRef] [Google Scholar]
 Uzdensky, D. A., & McKinney, J. C. 2011, Phys. Plasmas, 18, 042105 [NASA ADS] [CrossRef] [Google Scholar]
 Vay, J.L. 2008, Phys. Plasmas, 15, 056701 [NASA ADS] [CrossRef] [Google Scholar]
 Villasenor, J., & Buneman, O. 1992, Comput. Phys. Commun., 69, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Comput. Phys. Commun., 204, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302 [Google Scholar]
All Figures
Fig. 1. Schematic representation of a spherical PIC grid: Panel (a) shows the grid layout and identifies the coordinate system and boundary types, while panel (b) shows the edges of the grid cell where each field component is defined. 

In the text 
Fig. 2. Particle pusher benchmarks corresponding to particle motions in (a12) a uniform azimuthal magnetic field, (b12) crossed constant magnetic and electric fields, and (c12) the timevarying electric and magnetic field components of TM modes. 

In the text 
Fig. 3. Schematic representation of (a) the spherical particle shape and (b) the variation of its flattop density value with the radial coordinate. The blue shaded region in (a) represents the particle shape and identifies its widths in the radial and polar directions. 

In the text 
Fig. 4. Particle shape properties. Panel (a) shows the radial width and panel (b) the density and real particle number. 

In the text 
Fig. 5. Schematic representation of the currentdeposition algorithm in (a) Cartesian and (b) spherical coordinates (see text for details). 

In the text 
Fig. 6. Current deposition benchmarks, showing that both (a) the continuity equation and (b) Gauss’ law are satisfied to machine precision. 

In the text 
Fig. 7. Forcefree magnetosphere obtained with volume injection. Panels (af) show the electron and positron density, total charge density, radial current density, azimuthal magnetic field, and electric field component parallel to the local magnetic field, respectively. Quantities are multiplied by powers of r to enhance the features at large radii. Black and white solid lines represent magnetic field lines, and vertical dashed lines show the location of the light cylinder. 

In the text 
Fig. 8. Magnetospheric solutions obtained with volume injection. The panels show the total charge density after a stellar rotation period. 

In the text 
Fig. 9. Poynting flux in simulations with volume injection. Values are normalized to the theoretical value L_{0} = μ^{2}Ω^{4}/c^{3}. 

In the text 
Fig. 10. Radial and temporal dependencies of Poynting flux in simulations with volume injection. Panel (a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and panel (b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in (a) and (b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 

In the text 
Fig. 11. Magnetospheric solutions obtained with surface injection proportional to E_{∥}. Panels (a12) show the total charge density and E_{∥}, respectively for a simulation with v_{s}/c = 0 and panels (b12) show the same but for a simulation with v_{s}/c = 0.99. Solid lines represent magnetic field lines, and vertical dashed lines show the location of the light cylinder. 

In the text 
Fig. 12. Magnetospheric solutions obtained with surface injection proportional to n_{GJ} with fixed k_{s} = n_{s}/n_{GJ} = 0.2 and varying v_{s}/c. 

In the text 
Fig. 13. Radial and temporal dependencies of Poynting flux in simulations with surface injection proportional to n_{GJ} with fixed k_{s} = n_{s}/n_{GJ} = 0.2 and varying v_{s}/c. a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in a) and b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 

In the text 
Fig. 14. Magnetospheric solutions obtained with surface injection proportional to n_{GJ} with fixed v_{s} and varying k_{s} = n_{s}/n_{GJ}. 

In the text 
Fig. 15. Magnetospheric solutions obtained with pair production. Panels ad) show the total charge density for simulations with η = {5, 25, 50, 150}. Solid lines represent magnetic field lines, and vertical dashed lines show the location of the lightcylinder. 

In the text 
Fig. 16. Radial and temporal dependencies of Poynting flux in simulations with pair production with varying η. a) shows the timeaveraged luminosity ⟨L⟩ as a function of r after a stellar rotation period, and b) shows the temporal evolution of the surface Poynting flux L_{*}. The dashed lines in a) and b) identify the lightcylinder radius and the theoretical surface Poynting flux L_{0} = μ^{2}Ω^{4}/c^{3}, respectively. 

In the text 
Fig. 17. Cyclic pair production along the return current layers. Columns labeled 13 correspond to different times and rows labeled ad show the electron and positron densities, radial current and E_{∥}, respectively. Results obtained for η = 100. Solid lines represent magnetic field lines. 

In the text 
Fig. 18. Pairproduction sites for the same simulation and times shown in Figure 17. The color indicates the number of pairproduction events between data dumps (≃0.3 r_{*}/c) in each grid cell. 

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