© The Authors 2024

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 gamma-ray emission consistent with observations.

Particle-in-cell (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 kinetic-scale 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 pair-production model, the global asymptotic magnetospheric topology varies quite significantly: in some cases, the system autoregulates to a fully charge-separated configuration (also called electrosphere) that does not produce a Poynting flux, whereas in other cases the magnetosphere converges to a force-free 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 pair-production-based injection schemes.

In this work, we performed two-dimensional axisymmetric global simulations of pulsar magnetospheres with three different pair-injection 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 pair-production model parameters. We show that all plasma sources produce near force-free 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 subforce-free 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 two-dimensional 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 charge-conserving current-deposition 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 volume-based (Sect. 3.1), surface-based (Sect. 3.2), and pair-production-based (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 co-rotation 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 anti-parallel 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}_{\min }+(n-1)\mathrm{\Delta },n=1,2,...,N_r+1,\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/N_r 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 pair-production 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 pair-production 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}}^{n-1/2}=-c\mathrm{\Delta }t{(\mathrm{\nabla }\times \mathbf{E})}^n,\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},\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}·\mathrm{d}{\mathcal{C}}_{\mathrm{cell}})/{\mathcal{S}}_{\mathrm{cell}},\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})}.\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_isin θ_j + 1E_{ϕ(i, j + 1)} − r_isin θ_jE_ϕ(i, j)) and 2πr_i²(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 first-order 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 quasi-stationary 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ⁿ ≡ (r_p, θ_p), an operation that we write schematically as (E_(i, j)ⁿ, B_(i, j)ⁿ)→(E_pⁿ, B_pⁿ). 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 })},\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},\hfill \\ \hfill {f_{\theta }}_j& =1-{f_{\theta }}_{j+1}=\frac{cos{\theta }_j-cos{\theta }_p}{cos{\theta }_j-cos{\theta }_{j+1}}.\hfill \end{array}$$(6)

After the interpolation, the field components are converted from spherical to Cartesian coordinates, (E_pⁿ, B_pⁿ)→(E_p, Cⁿ, B_p, Cⁿ), a calculation that depends on the particle position at time tⁿ, xⁿ. Finally, the particle momentum and position are updated in time, u^{n − 1/2} ≡ p^{n − 1/2}/m_ec → u^n + 1/2 ≡ p^n + 1/2/m_ec and xⁿ → 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 macro-particle 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 time-varying 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 self-consistently created due to particle motion via plasma currents.

Fig. 2. Particle pusher benchmarks corresponding to particle motions in (a1-2) a uniform azimuthal magnetic field, (b1-2) crossed constant magnetic and electric fields, and (c1-2) the time-varying 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 macro-particles 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 }·\mathbf{j}=0,\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 }·\mathbf{E}=4\pi \rho ,\end{array}$$(10)

is satisfied. Finding a current deposition algorithm that satisfies Eq. (9), and consequently Eq. (10) (i.e., a charge-conserving 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 charge-conserving 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 macro-particle 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 macro-particle 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 macro-particle, such that its total charge remains constant.

Fig. 3. Schematic representation of (a) the spherical particle shape and (b) the variation of its flat-top 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 macro-particle as N_p, we formally wish to find a waterbag-like particle number density n(r) such that

$$\begin{array}{c}\hfill {\displaystyle {\int }_{V_i}}n(r_{i+1/2})\mathrm{d}{V}_i={\displaystyle {\int }_{V_{i^{\prime }}}}n(r_{i^{\prime }+1/2})\mathrm{d}{V}_{i^{\prime }}=N_p,\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{3N_p}{4\pi }\frac{1}{r_{i+1}^3-r_i^3}=\frac{3N_p}{32\pi }\frac{{(\delta +1)}^3}{{\delta }^3-1}\frac{1}{r_{i+1/2}^3},\end{array}$$(12)

where we have used the relation r_i + 1/2 = r_i(1 + δ)/2 = r_i + 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{3N_p}{32\pi }\frac{{(\delta +1)}^3}{{\delta }^3-1}\frac{1}{r_p^3}.\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)=2r_p\frac{\delta -1}{\delta +1}.\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),\hfill \end{array}$$(15)

where b₀(x) is the zeroth order b-spline function, defined as b₀(x) = 1 if |x|< 0.5 and 0 otherwise. We note that Eq. (15) generalizes the particle shape to a two-dimensional (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 macro-particle with N_p real particles of charge q_p and coordinates (r_p, θ_p) as ρ_p(r, θ, r_p, θ_p) = q_pN_pS(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)\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_{i-1/2}^3)(cos{\theta }_{j-1/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)}].\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}\mathrm{and}{\rho }_{(i,j)}(r_p,{\theta }_p)\propto \frac{r_{>}^3-r_{<}^3}{r_{i+1/2}^3-r_{i-1/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 charge-conserving 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ⁿ 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 macro-particle 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 current-deposition 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ⁿ 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 }·\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},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill {{(\mathrm{\nabla }·\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},\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 }·\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},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill {{(\mathrm{\nabla }·\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},\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 }·\mathbf{j})}_r}_{(i,j)}& =3[\frac{r_{i+1/2}^2{j_r}_{(i+1/2,j)}-r_{i-1/2}^2{j_r}_{(i-1/2,j)}}{r_{i+1/2}^3-r_{i-1/2}^3}],\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill {{(\mathrm{\nabla }·\mathbf{j})}_{\theta }}_{(i,j)}& =\frac{3}{2}\frac{r_{i+1/2}^2-r_{i-1/2}^2}{r_{i+1/2}^3-r_{i-1/2}^3}\times \hfill \\ \hfill & \times [\frac{sin{\theta }_{j+1/2}{j_{\theta }}_{(i,j+1/2)}-sin{\theta }_{j-1/2}{j_{\theta }}_{(i,j-1/2)}}{cos{\theta }_{j-1/2}-cos{\theta }_{j+1/2}}],\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ⁿ 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 macro-particle 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 first-order) 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 Cartesian-like 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 }·\mathbf{j})}_{(i,j)}),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{Gauss}}& =\frac{1}{{\rho }_{(i,j)}}({(\mathrm{\nabla }·\mathbf{E})}_{(i,j)}-4\pi {\rho }_{(i,j)}).\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⁻¹⁵ − 10⁻¹¹. 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 round-off 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 multi-particle 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_ec²/er_*, however we typically represent them in units of en_GJr_*, where n_GJ = ΩB_*/2πec is the surface Goldreich-Julian (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_ec²). For realistic parameters, B_* ≃ 10¹² G and r_* ≃ 10 km, we have B_*(er_*/m_ec²)∼10¹⁵. 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⁻⁷r_*. For this reason, we use more modest values of B_*(er_*/m_ec²)∼10³ − 10⁶, such that we respect the ordering in these objects, Ω ≪ ω_p, GJ ≪ ω_c, where ω_c = eB_*/m_ec 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 electron-positron 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)³cosθ and B_θ(r, θ) = (1/2)B_*(r_*/r)³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 co-rotation 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_risec/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 light-cylinder 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 electron-positron 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 macro-electron-positron pair with particle weight w_vol = k_volV_pE_∥/er_* and carrying a number density n_vol = w_vol/V_p, with k_vol = 0.2 and V_p being the macro-particle’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 macro-particles 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_ec² = 8 × 10³, N_r × N_θ = 1000² and Δtc/r_* = 10⁻³. In these conditions, c/ω_p, GJr_* ≃ 0.022, whereas the minimum grid spacing is min(Δr_i)/r_* ≃ 0.003.

In Fig. 7, we present an overview of the quasi-steady-state solution obtained with k_lim = 0.005. This solution is achieved after a time ∼25 r_*/c ∼ T/2¹. In the first half stellar period, the simulation undergoes a transient stage in which the vacuum co-rotation fields are established and plasma is created. The solution presented in Fig. 7 resembles the canonical force-free 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 light-cylinder 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. 7a-c, 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. 7a-b. 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 force-free regime for pulsar magnetospheres.

Fig. 7. Force-free magnetosphere obtained with volume injection. Panels (a-f) 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 quasi-steady-state 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 force-free 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 force-free 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 charge-separated configuration, with only electrons (positrons) in the poles (equatorial region). This solution is often denominated as the disk-dome or electrosphere solution (Jackson 1976; Krause-Polstorff & 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{r}^{2}sin\theta \mathrm{d}\theta .\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₀ = μ²Ω⁴/c³, with μ = B_*r_*³. 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₀ ≪ 1, which is a consequence of the inactivity of disk-dome solution. On the contrary, the simulations with lower k_lim have L_*/L₀ ∼ 1. The Poynting flux remains approximately constant within the light-cylinder 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 time-averaged 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 L0 = μ2Ω4/c3.

Fig. 10. Radial and temporal dependencies of Poynting flux in simulations with volume injection. Panel (a) shows the time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, 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 force-free 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 quasi-force-free 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 spin-down power, reflected in the observed luminosity and also in the poloidal extent of the polar cap. Here, we show that pulsars with quasi-force-free magnetospheres (e.g., k_lim = 0.005 simulation) must be observed with ⟨L⟩∼L₀ and the radio beam poloidal extent, if resolved, must be close to the polar cap angle provided by force-free models, ${\theta }_{\mathrm{PC}}\sim \sqrt{r_{\ast }/R_{\mathrm{LC}}}$. This is in contrast with subforce-free magnetospheres (e.g., k_lim = 0.01 simulation), resembling the weak-pulsar 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 spin-down luminosities, ⟨L⟩≲0.8L₀, and thinner radio beams, due to the compression of the polar cap, making them distinguishable from regular force-free 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 over-injection 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 light-cylinder radius is only possible for pulsars with sufficient optical depth for the activation of the γ − γ pair-production channel (Cheng et al. 1986). It is then important to assess if more realistic injection and/or pair-production schemes can provide the plasma supply required for the magnetosphere to be in the force-free 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 surface-injected 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_ec² = 8 × 10³, N_r × N_θ = 500² and Δtc/r_* = 3 × 10⁻³.

The first injection criterion is based on the local value of E_∥. We inject a macro-electron-positron 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_sn_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 disk-dome 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 charge-separated solution and effectively screens E_∥ at the stellar surface, precluding further injection.

Fig. 11. Magnetospheric solutions obtained with surface injection proportional to E∥. Panels (a1-2) show the total charge density and E∥, respectively for a simulation with vs/c = 0 and panels (b1-2) show the same but for a simulation with vs/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 first-principles 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 electron-positron pairs carry a number density n_s = k_sn_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 space-charge 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 force-free 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 nGJ with fixed ks = ns/nGJ = 0.2 and varying vs/c.

Figure 13 shows the time-averaged 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₀, which is a consequence of the smaller size of the co-rotation region (and thus a smaller effective light-cylinder radius and larger effective L₀). This particular simulation demonstrates that the over-injection drives the system to expand the effective polar cap area, supporting an artificially stronger than force-free 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 nGJ with fixed ks = ns/nGJ = 0.2 and varying vs/c. a) shows the time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, 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 steady-state (or quasi-steady-state) of these simulations is shown in Fig. 14. These results confirm that the denser the injected plasma is, the more the solution approaches the force-free regime (see in particular the solution obtained for k_s = 0.5). This injection density requirement seems to be critical in the launching of high-density 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 steady-state magnetospheric solution from charge-starved to force-free, having implications on the observed spin-down luminosity and radio beam width.

Fig. 14. Magnetospheric solutions obtained with surface injection proportional to nGJ with fixed vs and varying ks = ns/nGJ.

Consistently with simulations presented in Sect. 3.1, intermediate supply of plasma leads to a subforce-free magnetospheric solution with a thick current sheet configuration (see Figs. 8b, 12b-d and 14a-b). At the location of the light-cylinder 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₀ ∼ 1 and a global configuration similar to the force-free regime. Weak partially charge-starved pulsar magnetospheric solutions with persistent outer gaps and lower time-averaged 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 force-free 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 pair-production channel enough to supply the plasma to fill the whole magnetosphere?

In young and rapidly rotating pulsars (e.g., the Crab pulsar and other gamma-ray pulsars), pairs can also be created via the γ-γ channel (Cheng et al. 1986). In this process, for which the cross-section peaks at around a center of mass energy ∼2 m_ec², gamma-rays 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 low-energy 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 pair-production model described in Cruz et al. (2021b, 2022), in which a lepton emits a pair of combined energy γ_pairm_ec² 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_ec² is the maximum energy achievable by the particles in the voltage Φ_pc = B_*r_*³Ω²/c² 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 pair-production 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² and Δtc/r_* = 3 × 10⁻³, for η = {25, 50} we use N_r × N_θ = 1000² and Δtc/r_* = 10⁻³ and for η = {100, 150} we use N_r × N_θ = 2000² and Δtc/r_* = 5 × 10⁻⁴.

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 cross-section in this region (Cruz et al. 2021a). Seed electron-positron 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_sE_∥/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 force-free-like 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 light-cylinder 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 force-free 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), quasi-neutral plasma to the light-cylinder. 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 a-d) 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 light-cylinder.

The time-averaged 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 force-free spindown luminosity L₀ within the light-cylinder. 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 time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, respectively.

All simulations present some temporal variability. We see small-scale 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 quasi-periodic launch of plasma towards the light-cylinder 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, re-establishing j_r and effectively screening the E_∥ responsible for triggering the process – see Figure 17d1-3). 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 1-3 correspond to different times and rows labeled a-d 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 Y-point 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 time-dependent 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 low-obliquity pulsars (Torres et al. 2024).

Fig. 18. Pair-production sites for the same simulation and times shown in Figure 17. The color indicates the number of pair-production 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 force-free 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 surface-injection schemes (Cerutti et al. 2015; Hakobyan et al. 2023), and those that use heuristic pair-production 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 surface-plasma injection serve as a means to efficiently fill the pulsar magnetosphere and produce a near force-free configuration, as shown in Sects. 3.1 and 3.2, respectively, these are hard to motivate from first-principle 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 surface-injection 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 pair-producing 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 force-free 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 charge-separated, 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 force-free regime.

Our simulations show that the pair production along the return current layers is essential in order to feed plasma to the light-cylinder region and beyond in near force-free regimes, in line with the results reported in other works, such as Chen & Beloborodov (2014). We also identify a time-dependent 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 millisecond-type 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 charge-starved and pair production is no longer possible (disk–dome solution). From an observational standpoint, we may distinguish these solutions from the regular force-free magnetospheres by exploring the compression of the open-flux tube caused by the presence of a persistent outer gap and volumetric return currents, as seen for all plasma-injection methods. For more compact magnetospheres, capable of γ − γ pair production near the Y-point (Cheng et al. 1986), these may lead to episodic switches of the whole magnetosphere between weak and near force-free 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 force-free regime. We defer further exploration of this topic to future work.

The pair-production 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 pair-production 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 charge-conserving 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 (ERC-2015-AdG 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 X-MASER project no. 2022.02230.PTDC). AC acknowledges support from NSF grants DMS-2235457 and AST-2308111. AS is supported by the NSF grant PHY-2206607 and the Simons Foundation grant MP-SCMPS-00001470. We acknowledge PRACE (PULSAR project) for granting access to MareNostrum, Barcelona Supercomputing Center (Spain), where the simulations presented in this work were performed.

References

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 (a1-2) a uniform azimuthal magnetic field, (b1-2) crossed constant magnetic and electric fields, and (c1-2) the time-varying 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 flat-top 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 current-deposition 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. Force-free magnetosphere obtained with volume injection. Panels (a-f) 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 L0 = μ2Ω4/c3.
In the text

Fig. 10. Radial and temporal dependencies of Poynting flux in simulations with volume injection. Panel (a) shows the time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, respectively.

In the text

	Fig. 11. Magnetospheric solutions obtained with surface injection proportional to E∥. Panels (a1-2) show the total charge density and E∥, respectively for a simulation with vs/c = 0 and panels (b1-2) show the same but for a simulation with vs/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 nGJ with fixed ks = ns/nGJ = 0.2 and varying vs/c.
In the text

Fig. 13. Radial and temporal dependencies of Poynting flux in simulations with surface injection proportional to nGJ with fixed ks = ns/nGJ = 0.2 and varying vs/c. a) shows the time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, respectively.

In the text

	Fig. 14. Magnetospheric solutions obtained with surface injection proportional to nGJ with fixed vs and varying ks = ns/nGJ.
In the text

	Fig. 15. Magnetospheric solutions obtained with pair production. Panels a-d) 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 light-cylinder.
In the text

Fig. 16. Radial and temporal dependencies of Poynting flux in simulations with pair production with varying η. a) shows the time-averaged 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 light-cylinder radius and the theoretical surface Poynting flux L0 = μ2Ω4/c3, respectively.

In the text

	Fig. 17. Cyclic pair production along the return current layers. Columns labeled 1-3 correspond to different times and rows labeled a-d 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. Pair-production sites for the same simulation and times shown in Figure 17. The color indicates the number of pair-production events between data dumps (≃0.3 r*/c) in each grid cell.
In the text