11institutetext: GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior
Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
11email<EMAIL_ADDRESS>22institutetext: Inductiva Research
Labs, Rua da Prata 80, 1100-420 Lisboa, Portugal 33institutetext: Physics
Department and McDonnell Center for the Space Sciences, Washington University
in St. Louis, MO 63130, USA 44institutetext: Department of Astrophysical
Sciences, Princeton University, Princeton, NJ 08544, USA 55institutetext:
DCTI/ISCTE Instituto Universitário de Lisboa, 1649-026 Lisboa, Portugal

# Particle-in-cell simulations of pulsar magnetospheres: transition between
electrosphere and force-free regimes

F. Cruz 1122 T. Grismayer 11 A. Y. Chen 33 A. Spitkovsky 44 R. A. Fonseca 1155
L. O. Silva 11

(Received )

###### Abstract

Aims. Global particle-in-cell (PIC) simulations of pulsar magnetospheres are
performed with a volume, surface and pair production-based plasma injection
schemes to systematically investigate the transition between electrosphere and
force-free pulsar magnetospheric regimes.

Methods. A new extension of the PIC code OSIRIS to model pulsar magnetospheres
using a two-dimensional axisymmetric spherical grid is presented. The sub-
algorithms of the code and thorough benchmarks are presented in detail,
including a new first-order current deposition scheme that conserves charge to
machine precision.

Results. It is shown 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 plasma injection rates,
and with pair production-based plasma injection for sufficiently large
separation between kinematic and pair production energy scales.

###### Key Words.:

pulsars – electrosphere – force-free magnetosphere – particle-in-cell –
simulations

##  1 Introduction

Over the last decade, global kinetic simulations have been essential tools to
understand 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., 2015; 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., 2015; 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 with high
fidelity the kinetic plasma phenomena relevant in pulsars, such as the
evolution of highly non-thermal particle distributions or kinetic-scale
fluctuations (Touati et al., 2022). Recent extensions of the PIC method have
also allowed the inclusion of Quantum Electrodynamics effects such as pair
production (Grismayer et al., 2016, 2017) or general relativity corrections
(Philippov et al., 2015) relevant in pulsars.

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. Such
description can be as simple as injecting plasma in a significant fraction the
simulation domain (Philippov & Spitkovsky, 2014; Belyaev, 2015a;
Kalapotharakos et al., 2018; Brambilla et al., 2018), limiting this injection
to occur 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., 2015, 2015; Chen et al., 2020; Guépin et
al., 2020; Bransgrove et al., 2022).

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 auto-regulates 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 perform two-dimensional axisymmetric global simulations of
pulsar magnetospheres with three different pair injection schemes: over large
volumes of the magnetosphere, from the stellar surface only and using a
prescription model for pair production. We use these simulations to 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 sub-force-
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) developed for
magnetospheric models of compact objects, presented also 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 (Sect.
3.1), surface (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 aims to model the global plasma
environment surrounding neutron stars, i.e., the spatial volume between the
stellar surface and a few light cylinder radii above it. We describe this
system in spherical coordinates, with the radial coordinate $r$ measured from
the center of the neutron star and the polar angle $\theta$ measured from the
star’s rotation axis $\boldsymbol{\Omega}$. We assume that
$\boldsymbol{\Omega}$ is either parallel or anti-parallel to the star’s
magnetic axis $\boldsymbol{\mu}$, such that we can assume axisymmetry about
$\boldsymbol{\Omega}$, i.e., derivatives with respect to the azimuthal angle
$\phi$ can be dropped, $\partial/\partial\phi=0$.

Similarly to Chen & Beloborodov (2014); Cerutti et al. (2015), we discretize
the simulation domain $r\in[r_{\mathrm{min}},r_{\mathrm{max}}]$,
$\theta\in[0,\pi]$ in a grid with $N_{r}\times N_{\theta}$ cells. We adopt a
regular grid spacing in $\theta$, $\Delta\theta=\pi/(N_{\theta}+1)$, 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

$\log r_{n}=\log r_{\mathrm{min}}+(n-1)\Delta\ ,\ \ n=1,2,...,N_{r}+1\ ,$ (1)

with $\Delta\equiv\log(r_{\mathrm{max}}/r_{\mathrm{min}})/N_{r}$. Equation (1)
can be manipulated to write the useful relation
$r_{n}=r_{\mathrm{min}}\delta^{n-1}$, where
$\delta\equiv(r_{\mathrm{max}}/r_{\mathrm{min}})^{1/N_{r}}$ is a parameter
that combines all properties of the radial axis.

Figure 1: Schematic representation of spherical PIC grid: a) shows the grid
layout and identifies the coordinate system and boundary types, b) shows the
grid cell’s edges where each field component is defined.

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_{\mathrm{min}}=r_{*}$,
whereas the upper radial boundary is at $r_{\mathrm{max}}\sim$ tens of
$r_{*}$, and acts as an open boundary. The $\theta={0,\pi}$ boundaries enforce
axisymmetry, effectively serving as reflecting boundaries. More details about
these boundaries are provided in Sects. 2.2 and 2.3.

In Fig. 1b, we show a schematic representation of a typical grid cell, that 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}\equiv(r_{i}+r_{i+1})/2$ and
$\theta_{j+1/2}\equiv(\theta_{j}+\theta_{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,

$\displaystyle\mathbf{B}^{n+1/2}=\mathbf{B}^{n-1/2}=-c\Delta
t(\nabla\times\mathbf{E})^{n}\ ,$ (2)
$\displaystyle\mathbf{E}^{n+1}=\mathbf{E}^{n}+c\Delta
t(\nabla\times\mathbf{B})^{n+1/2}-4\pi\Delta t\mathbf{j}^{n+1/2}\ ,$ (3)

where quantities with integer/half integer superscripts are defined in
integer/half integer times and $\Delta 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 by using Stokes’ theorem to
evaluate the curl of electric and magnetic fields in a given cell as

$(\nabla\times\mathbf{E})_{\mathrm{cell}}=\left(\oint_{\mathcal{C}_{\mathrm{cell}}}\mathbf{E}\cdot\mathrm{d}\mathcal{C}_{\mathrm{cell}}\right)/\mathcal{S}_{\mathrm{cell}}\
,$ (4)

where $\mathcal{C}_{\mathrm{cell}}$ is the contour defining the edge of that
cell, $\mathcal{S}_{\mathrm{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 $\nabla\times\mathbf{E}$ as

${(\nabla\times\mathbf{E})_{r}}_{(i,j+1/2)}=\frac{\sin\theta_{j+1}{E_{\phi}}_{(i,j+1)}-\sin\theta_{j}{E_{\phi}}_{(i,j)}}{r_{i}(\cos\theta_{j}-\cos\theta_{j+1})}\
.$ (5)

This expression is derived by noting that, according to Eq. (2),
$(\nabla\times\mathbf{E})_{r}$ should be defined in the same position as
$\mathbf{B}_{r}$, i.e., at cell indices $(i,j+1/2)$. This defines the
integration surface relevant to Stokes’ theorem as $r=r_{i}$,
$\theta\in[\theta_{j},\theta_{j+1}]$. The numerator and denominator in Eq. (4)
then read respectively
$2\pi(r_{i}\sin\theta_{j+1}{E_{\phi}}_{(i,j+1)}-r_{i}\sin\theta_{j}{E_{\phi}}_{(i,j)})$
and $2\pi r_{i}^{2}(\cos\theta_{j}-\cos\theta_{j+1})$, where the $2\pi$ factor
comes from the integration along $\phi$. 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_{\theta}+1\\}$), the integration regions are adapted to fit within
the domain. For example, the $\theta$ integration is changed to
$\theta\in[0,\theta_{1+1/2}]$ and $\theta\in[\theta_{N_{\theta}+1/2},\pi]$ at
the $\theta=0$ and $\theta=\pi$ boundaries, respectively. We also apply
special rules to the field components at the boundaries, e.g. in the polar
boundaries we enforce the axisymmetry conditions
${\mathbf{E}_{\phi}}_{(i,1)}={\mathbf{E}_{\phi}}_{(i,N_{\theta}+1)}=0$ and
${\mathbf{B}_{\theta}}_{(i+1/2,1)}={\mathbf{B}_{\theta}}_{(i+1/2,N_{\theta}+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), i.e., a perfect absorber of perturbations
propagating perpendicularly to the boundary. We have also implemented static
conductor boundary conditions for both inner and outer radial boundaries, that
enforce tangent (normal) electric (magnetic) field components to be null,
i.e., ${\mathbf{E}_{\phi}}_{(1,j)}={\mathbf{E}_{\phi}}_{(N_{r}+1,j)}=0$,
${\mathbf{E}_{\theta}}_{(1,j+1/2)}={\mathbf{E}_{\theta}}_{(N_{r}+1,j+1/2)}=0$
and ${\mathbf{B}_{r}}_{(1,j+1/2)}={\mathbf{B}_{r}}_{(N_{r}+1,j+1/2)}=0$.

We have benchmarked our field solver implementation by studying stationary
electromagnetic TM modes between two spherical static conductors (Jackson,
1975). We have 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 sub-algorithms 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 $\mathbf{x}_{p}^{n}\equiv(r_{p},\theta_{p})$, an operation that we
write schematically as
$(\mathbf{E}^{n}_{(i,j)},\mathbf{B}^{n}_{(i,j)})\to(\mathbf{E}^{n}_{p},\mathbf{B}^{n}_{p})$.
This interpolation is done using a area/volume weighting scheme. For example,
the toroidal component of the electric field can be written as

${\mathbf{E}_{\phi}}_{p}=\sum_{i^{\prime}=i,i+1}\sum_{j^{\prime}=j,j+1}{f_{r}}_{i^{\prime}}{f_{\theta}}_{j^{\prime}}{\mathbf{E}_{\phi}}_{(i^{\prime},j^{\prime})}\
,$ (6)

with

$\displaystyle{f_{r}}_{i}$
$\displaystyle=1-{f_{r}}_{i+1}=\frac{r_{p}^{3}-r_{i}^{3}}{r_{i+1}^{3}-r_{i}^{3}}\
,$ $\displaystyle{f_{\theta}}_{j}$
$\displaystyle=1-{f_{\theta}}_{j+1}=\frac{\cos\theta_{j}-\cos\theta_{p}}{\cos\theta_{j}-\cos\theta_{j+1}}\
.$

After the interpolation, the field components are converted from spherical to
Cartesian coordinates,
$(\mathbf{E}^{n}_{p},\mathbf{B}^{n}_{p})\to(\mathbf{E}^{n}_{p,\mathrm{C}},\mathbf{B}^{n}_{p,\mathrm{C}})$,
a calculation that depends on the particle position at time $t^{n}$,
$\mathbf{x}^{n}$. Finally, the particle momentum and position are updated in
time,
$\mathbf{u}^{n-1/2}\equiv\mathbf{p}^{n-1/2}/m_{e}c\to\mathbf{u}^{n+1/2}\equiv\mathbf{p}^{n+1/2}/m_{e}c$
and $\mathbf{x}^{n}\to\mathbf{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). We note that advancing the particle position
in $(x,y,z)$ does not introduce any asymmetry in the azimuthal direction
$\phi$; in fact, each macro-particle in our simulation represents a charged
ring with azimuthal symmetry and $\phi$ is never used throughout the rest of
the numerical scheme.

Figure 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.

We have 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 $\sim$ 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, i.e., they do not include the fields self-consistently created due to
particle motion via plasma currents.

###  2.4 Current deposition

A current deposition algorithm computes the current density $\mathbf{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),

$\frac{\partial\rho}{\partial t}+\nabla\cdot\mathbf{j}=0\ ,$ (7)

where $\rho$ is the total plasma density. Solving Eq. (7) ensures also that
Gauss’ law, written as

$\nabla\cdot\mathbf{E}=4\pi\rho\ ,$ (8)

is satisfied. Finding a current deposition algorithm that satisfies Eq. (7),
and consequently Eq. (8), 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 non-uniform spherical grids, this challenge is
more substantial, as grid cells (and particle shapes, that we shall define
below) change across the grid. Other codes adopting such grids (Chen &
Beloborodov, 2014; Cerutti et al., 2015; Belyaev, 2015b; Chen, 2017) usually
do not seem to include charge-conserving current deposition algorithms, and
adopt instead numerical schemes to enforce the validity of Eq. (8), e.g.
Poisson solvers.

Figure 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.

Here, we propose a new current deposition scheme that conserves charge to
machine precision in the non-uniform grid defined in section 2.1. We start by
defining the volume occupied by a macro-particle centered at
$(r_{p},\theta_{p})$. The function that defines this volume is usually called
the particle shape, $S(r,\theta,r_{p},\theta_{p})$. Before writing the exact
form of $S$, let us define some of its important properties, that 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},\theta_{p})=(r_{i+1/2},\theta_{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}$, i.e., $\Delta r\equiv\Delta 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.

Figure 4: Particle shape properties: a) radial width and b) density and real
particle number.

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

$\int_{V_{i}}n(r_{i+1/2})\
\mathrm{d}V_{i}=\int_{V_{i^{\prime}}}n(r_{i^{\prime}+1/2})\
\mathrm{d}V_{i^{\prime}}=N_{p}\ ,$ (9)

where $V_{i,i^{\prime}}$ are the volumes of cells with radial labels
$i,i^{\prime}$ (see Figure 3 b)). For simplicity, we assume that the particle
density is only a function of $r$, and generalize it later to include the
natural dependence in $\theta$ as well. Assuming that $n(r_{i+1/2})$ is
constant within cell $i$, we can solve Eq. (9) to obtain

$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{\left(\delta+1\right)^{3}}{\delta^{3}-1}\frac{1}{r_{i+1/2}^{3}}\
,$ (10)

where we have used the relation
$r_{i+1/2}=r_{i}(1+\delta)/2=r_{i+1}(1+\delta^{-1})/2$. We note that Eq. (10)
defines $n(r)$ for any $r_{i+1/2}$, but not for $r\neq 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}$, i.e.,

$n(r_{p})=\frac{3N_{p}}{32\pi}\frac{\left(\delta+1\right)^{3}}{\delta^{3}-1}\frac{1}{r_{p}^{3}}\
.$ (11)

Eq. (11) ensures that $n(r)$ satisfies exactly Eq. (9) when $r_{p}=r_{i+1/2}$
and that the particle shape is a smooth function of $r_{p}$. The particle
width $\Delta r(r_{p})$ is determined in a similar manner; first, we express
the grid spacing in terms of $r_{i+1/2}$, $\Delta
r_{i}=r_{i+1}-r_{i}=2r_{i+1/2}(\delta-1)/(\delta+1)$, and we extend this
definition to an arbitrary radius $r_{p}$,

$\Delta r(r_{p})=2r_{p}\frac{\delta-1}{\delta+1}\ .$ (12)

This quantity is represented for a typical grid in Fig. 4a, together with the
grid spacing $\Delta r_{i}$. As expected, both quantities match exactly when
$r=r_{i+1/2}$, and $\Delta r$ is a smooth function of $r$. Equations (11) and
(12) ensure that the conservation law expressed in Eq. (9) can be extended to
any radius, which is shown in Fig. 4b.

The general particle shape $S$ can be inferred from this discussion, and in
particular from Eq. (11). It reads

$\displaystyle S(r,\theta,r_{p},\theta_{p})$
$\displaystyle=\frac{3}{16\pi}\frac{\left(\delta+1\right)^{3}}{\delta^{3}-1}\frac{1}{r_{p}^{3}}b_{0}\left(\frac{r-r_{p}}{\Delta
r(r_{p})}\right)\times$
$\displaystyle\times\frac{1}{\cos(\theta_{p}-\Delta\theta/2)-\cos(\theta_{p}+\Delta\theta/2)}b_{0}\left(\frac{\theta-\theta_{p}}{\Delta\theta}\right)\
,$ (13)

where $b_{0}(x)$ is the zeroth order b-spline function, defined as
$b_{0}(x)=1$ if $|x|<0.5$ and $0$ otherwise. Note that Eq. (13) generalizes
the particle shape to a two-dimensional $(r,\theta)$ grid, hence the
$\cos(\theta_{p}\pm\Delta\theta)$ terms resulting from the integral in Eq.
(9). With the shape function in Eq. (13), we can compute the charge density at
any point $(r,\theta)$ due to the presence of a macro-particle with $N_{p}$
real particles of charge $q_{p}$ and coordinates $(r_{p},\theta_{p})$ as
$\rho_{p}(r,\theta,r_{p},\theta_{p})=q_{p}N_{p}S(r,\theta,r_{p},\theta_{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

$\displaystyle\rho_{(i,j)}(r_{p},\theta_{p})$
$\displaystyle=\frac{\int_{V_{i,j}}\rho_{s}(r,\theta,r_{p},\theta_{p})\
\mathrm{d}V_{i,j}}{V_{i,j}}=$
$\displaystyle=\dfrac{q_{p}N_{p}\dfrac{3}{16\pi}\dfrac{\left(\delta+1\right)^{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$
$\displaystyle\times\Bigg{[}\frac{r_{>}^{3}-r_{<}^{3}}{r_{p}^{3}}\Bigg{]}\Bigg{[}\frac{\cos(\theta_{p}-\Delta\theta/2)-\cos\theta_{j+1/2}}{\cos(\theta_{p}-\Delta\theta/2)-\cos(\theta_{p}+\Delta\theta/2)}\Bigg{]}\
.$ (14)

We note that the special integration limits $r_{>}=\min(r_{p}+\Delta
r(r_{p})/2,r_{i+1/2})$ and $r_{<}=\max(r_{p}-\Delta 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.

Figure 5: Schematic representation of the current deposition algorithm in a)
Cartesian and b) spherical coordinates (see text for details).

The particle shape in Eq. (13) and the deposition rule in Eq. (14) are the key
ingredients in our charge-conserving current deposition scheme. This scheme is
inspired by the seminal work of Villasenor & Buneman (1992) (hereafter VB),
that presented a scheme that predecessed the widely used method of Esirkepov
(2001) for PIC current deposition in Cartesian grids. The VB method is
schematically represented in Fig. 5a. VB proposed that the current density
$\mathbf{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 $\mathbf{x}^{n}$ to a position
$\mathbf{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 $\mathbf{j}$ are defined
are shown in Fig. 5a in green and red lines and arrows, respectively. VB
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_{\mathrm{green}}$ and $A_{\mathrm{red}}$ and the total area corresponding
to the particle shape, $A_{\mathrm{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.

A schematic representation of the method equivalent to VB 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_{\mathrm{?}}$.
Because the particle expands during its motion from $\mathbf{x}^{n}$ to
$\mathbf{x}^{n+1}$, it is not trivial to determine which fraction of
$A_{\mathrm{?}}$ should be combined with $A_{\mathrm{green}}$
($A_{\mathrm{red}}$) to compute ${j_{r}}_{(i+1/2,j)}$
(${j_{\theta}}_{(i,j+1/2)}$). We circumvent this issue by generalizing the
geometrical interpretation of $\nabla\cdot\mathbf{j}$ proposed by VB. They
suggested that the total current divergence can be split as
$\nabla\cdot\mathbf{j}=(\nabla\cdot\mathbf{j})_{x}+(\nabla\cdot\mathbf{j})_{y}$
in a Cartesian grid, with $(\nabla\cdot\mathbf{j})_{y}\propto
A_{\mathrm{green}}/A_{\mathrm{total}}$ and $(\nabla\cdot\mathbf{j})_{x}\propto
A_{\mathrm{red}}/A_{\mathrm{total}}$, and that these terms could be computed
directly by evaluating $-\partial\rho_{(i,j)}/\partial t$ assuming that the
particle moves purely along the corresponding direction at an average position
along the orthogonal direction. Formally, this is expressed as

$\displaystyle{(\nabla\cdot\mathbf{j})_{x}}_{(i,j)}$
$\displaystyle=-\left.\frac{\partial\rho_{(i,j)}}{\partial
t}\right|_{x^{n},\bar{y}}^{x^{n+1},\bar{y}}=\frac{\rho_{(i,j)}(x^{n+1},\bar{y})-\rho_{(i,j)}(x^{n},\bar{y})}{\Delta
t}\ ,$ (15) $\displaystyle{(\nabla\cdot\mathbf{j})_{y}}_{(i,j)}$
$\displaystyle=-\left.\frac{\partial\rho_{(i,j)}}{\partial
t}\right|_{\bar{x},y^{n}}^{\bar{x},y^{n+1}}=\frac{\rho_{(i,j)}(\bar{x},y^{n+1})-\rho_{(i,j)}(\bar{x},y^{n})}{\Delta
t}\ ,$ (16)

where $\bar{x}=(x^{n+1}+x^{n})/2$ and $\bar{y}=(y^{n+1}+y^{n})/2$. From Eqs.
(15) and (16), 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, i.e., we can write
$\nabla\cdot\mathbf{j}=(\nabla\cdot\mathbf{j})_{r}+(\nabla\cdot\mathbf{j})_{\theta}$.
However, because the particle shape changes continuously in the radial
direction, $(\nabla\cdot\mathbf{j})_{\theta}$ cannot be computed assuming that
the particle moves purely along the polar direction with
$\bar{r}=(r^{n+1}+r^{n})/2$. Instead, we proceed as follows: first, we compute
$\nabla\cdot\mathbf{j}$ and $(\nabla\cdot\mathbf{j})_{r}$ using

$\displaystyle{(\nabla\cdot\mathbf{j})}_{(i,j)}$
$\displaystyle=-\left.\frac{\partial\rho_{(i,j)}}{\partial
t}\right|_{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})}{\Delta
t}\ ,$ (17) $\displaystyle{(\nabla\cdot\mathbf{j})_{r}}_{(i,j)}$
$\displaystyle=-\left.\frac{\partial\rho_{(i,j)}}{\partial
t}\right|_{r^{n},\bar{\theta}}^{r^{n+1},\bar{\theta}}=\frac{\rho_{(i,j)}(r^{n+1},\bar{\theta})-\rho_{(i,j)}(r^{n},\bar{\theta})}{\Delta
t}\ ,$ (18)

where $\bar{\theta}=(\theta^{n+1}+\theta^{n})/2$. Then, we compute
$(\nabla\cdot\mathbf{j})_{\theta}=\nabla\cdot\mathbf{j}-(\nabla\cdot\mathbf{j})_{r}$.
Finally, we invert the nabla operators,

$\displaystyle{(\nabla\cdot\mathbf{j})_{r}}_{(i,j)}$
$\displaystyle=3\Bigg{[}\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}}\Bigg{]}\
,$ (19) $\displaystyle{(\nabla\cdot\mathbf{j})_{\theta}}_{(i,j)}$
$\displaystyle=\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$
$\displaystyle\times\Bigg{[}\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}}\Bigg{]}\
,$ (20)

to find the current components. The inversion of
${(\nabla\cdot\mathbf{j})_{\theta}}_{(i,j)}$ is simple, because the second
term in the square brackets of Eq. (20) is always zero given that the particle
motion is restricted to cell $(i,j)$. The same is applicable to the inversion
of ${(\nabla\cdot\mathbf{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
${(\nabla\cdot\mathbf{j})_{r}}_{(i-1,j)}$ using Eq. (18), 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. (19). When particles cross two cells from
$\mathbf{x}^{n}$ to $\mathbf{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_{\phi}}_{(i,j)}=\rho_{(i,j)}v_{\phi}$, where $v_{\phi}$ is the
macro-particle velocity in the azimuthal direction.

Finally, we note that Eqs. (14) and (19)-(20) 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.
However, the special radial integration rule described in this section to
account for particle shrinking/expansion should be included to ensure that
those algorithms conserve charge to machine precision.

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

We have benchmarked the current deposition method presented here by
initializing particles all over the simulation domain with a random velocity,
depositing their current over a time step $\Delta t$ and evaluating

$\displaystyle\Delta_{\mathrm{Continuity}}$ $\displaystyle=\frac{\Delta
t}{\rho_{(i,j)}}\left(\frac{\partial\rho_{(i,j)}}{\partial
t}+(\nabla\cdot\mathbf{j})_{(i,j)}\right)\ ,$ (21)
$\displaystyle\Delta_{\mathrm{Gauss}}$
$\displaystyle=\frac{1}{\rho_{(i,j)}}\left((\nabla\cdot\mathbf{E})_{(i,j)}-4\pi\rho_{(i,j)}\right)\
.$ (22)

Both $\Delta_{\mathrm{Continuity}}$ and $\Delta_{\mathrm{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 $10^{-15}-10^{-11}$, i.e., of
the order of machine precision. The value of both
$\Delta_{\mathrm{Continuity}}$ and $\Delta_{\mathrm{Gauss}}$ tends to be
larger closer to the star, due to the larger number of operations subject to
round-off errors in this region, caused by particles crossing more cell
boundaries and depositing their current in more than one cell. We have
verified that the accuracy of the method is maintained over multiple time
steps by ensuring that the evolution of the grid integrals of
$\Delta_{\mathrm{Continuity}}$ and $\Delta_{\mathrm{Gauss}}$ remain at machine
precision level.

This current deposition method thus accurately conserves charge, avoiding the
need for other correcting algorithms. It is also inexpensive, since most
factors in Eqs. 17-20 can be precomputed and reused throughout a simulation.

###  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_{\mathrm{GJ}}r_{*}$, where $n_{\mathrm{GJ}}=\Omega B_{*}/2\pi 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,\mathrm{GJ}}=c/\omega_{p,\mathrm{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_{*}\simeq 10^{12}$ G and $r_{*}\simeq 10$ km, we have
$B_{*}(er_{*}/m_{e}c^{2})\sim 10^{15}$. Global simulations are not feasible
with such values, since they would have to resolve scales of the order of
$\sim$ tens of $r_{*}$ down to $d_{e,\mathrm{GJ}}\sim
10^{-7}\leavevmode\nobreak\ r_{*}$. For this reason, we use more modest values
of $B_{*}(er_{*}/m_{e}c^{2})\sim 10^{3}-10^{6}$, such that we respect the
ordering in these objects, $\Omega\ll\omega_{p,\mathrm{GJ}}\ll\omega_{c}$,
where $\omega_{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 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_{*}$, i.e.,
$B_{r}(r,\theta)=B_{*}(r_{*}/r)^{3}\cos\theta$ and
$B_{\theta}(r,\theta)=(1/2)B_{*}(r_{*}/r)^{3}\sin\theta$. The inner radial
boundary is treated as a rotating conductor of angular velocity
$\boldsymbol{\Omega}=\Omega\mathbf{\hat{z}}$; at the surface of the neutron
star, we impose the co-rotation electric field
$\mathbf{E}=-(\mathbf{v}_{\mathrm{rot}}\times\mathbf{B})/c$, with
$\mathbf{v}_{\mathrm{rot}}=\boldsymbol{\Omega}\times(r_{*}\mathbf{\hat{r}})$.
In all simulations, we consider the stellar rotation frequency to be initially
zero and increase it linearly over a time $t_{\mathrm{rise}}c/r_{*}=1.5$ to
$\Omega r_{*}/c=0.125$. For times $t>t_{\mathrm{rise}}$, the stellar frequency
is kept constant. The stellar period is $T=2\pi/\Omega=50\leavevmode\nobreak\
r_{*}/c$ and the light-cylinder radius is $R_{\mathrm{LC}}/r_{*}=8$. All
simulations use also $r_{\mathrm{min}}/r_{*}=1$ and
$r_{\mathrm{max}}/r_{*}=20$, such that the plasma dynamics can be captured up
to $r/R_{\mathrm{LC}}>2$. The value of $B_{*}$ is chosen to satisfy the
ordering $\Omega\ll\omega_{p,\mathrm{GJ}}\ll\omega_{c}$ described in Sect. 2.5
while maintaining simulations numerically feasible. This choice and others
regarding e.g., 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_{\parallel}c/r_{*}\Omega B_{*}>k_{\mathrm{lim}}$, where
$k_{\mathrm{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_{\parallel}$ is sufficiently large, i.e., it neglects any
role of the local magnetic field magnitude or curvature. Since all fields (and
in particular $E_{\parallel}$) decay with $r$, the choice of
$k_{\mathrm{lim}}$ can also be interpreted as a spatial limitation to the
plasma supply: infinitely small values of $k_{\mathrm{lim}}$ allow plasma to
be injected up to $r\gg r_{*}$, whereas $k_{\mathrm{lim}}\sim 1$ restricts the
plasma supply to radii $r\sim r_{*}$. A macro-electron-positron pair carrying
a number density $n_{\mathrm{vol}}=k_{\mathrm{vol}}E_{\parallel}/er_{*}$, with
$k_{\mathrm{vol}}=0.2$, is injected at rest in each cell and time step in
which the injection condition is met. The choice of $k_{\mathrm{vol}}$ is such
that a few macro-particles are required to supply the charge density that
screens $E_{\parallel}$ and stops the injection. We can also interpret
$k_{\mathrm{vol}}$ as a parameter proportional to the local GJ density, since
$E_{\parallel}/er_{*}\sim n_{\mathrm{GJ}}$. In all the simulations presented
in this section, $B_{*}er_{*}/m_{e}c^{2}=8\times 10^{3}$, $N_{r}\times
N_{\theta}=1000^{2}$ and $\Delta tc/r_{*}=10^{-3}$. In these conditions,
$c/\omega_{p,\mathrm{GJ}}r_{*}\simeq 0.022$, whereas the minimum grid spacing
is $\min(\Delta r_{i})/r_{*}\simeq 0.003$.

Figure 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 large radii features. White/black solid lines represent magnetic
field lines, and vertical dashed lines show the location of the light-
cylinder.

In Fig. 7, we present an overview of the quasi-steady-state solution obtained
with $k_{\mathrm{lim}}=0.005$. This solution is achieved after a time $\sim
25\leavevmode\nobreak\ r_{*}/c\sim T/2$111This is not a universal result. In
fact, the time required by the system to achieve a steady-state (or quasi
steady-state) solution varies with the injection scheme, the stellar ramp-up
time $t_{\mathrm{rise}}$ and other initial and/or boundary conditions.. 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\simeq R_{\mathrm{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_{\phi}$, 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_{*}\simeq 12$ in Fig.
7a-b. The plasma supply in this simulation is large enough such that
$E_{\parallel}$ 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.

The quasi-steady-state shown in Fig. 7 is sustained via intermittent
injection, mainly along the return current layers. In these regions,
$E_{\parallel}$ 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 $\simeq
0.3-0.5\leavevmode\nobreak\ T$, and it may play a significant role in the
temporal evolution of the magnetospheric state. However, for
$k_{\mathrm{lim}}=0.005$ the solution never deviates significantly from the
force-free regime.

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

In order to demonstrate how the magnetospheric solution changes with
$k_{\mathrm{lim}}$, in Fig. 8 we compare the total charge density of the
solutions obtained with $k_{\mathrm{lim}}=\\{0.005,0.01,0.1\\}$. We recall
that $k_{\mathrm{lim}}$ is the minimum value of $E_{\parallel}c/r_{*}\Omega
B_{*}$ for which we inject plasma. It is clear that the force-free regime is
only observed for $k_{\mathrm{lim}}=0.005$. For $k_{\mathrm{lim}}=0.01$, the
equatorial current sheet (positively charged region at $r\gtrsim
R_{\mathrm{LC}}$) is wide and the return current layers are not positively
charged everywhere, and for $k_{\mathrm{lim}}=0.1$ the solution does not even
produce an outflow. In fact, by increasing $k_{\mathrm{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_{\parallel}$ with the reference value $\Omega B_{*}r_{*}/c$ (i.e., the
surface magnitude of the electric field in vacuum). Since the typical
magnitude of $E_{\parallel}$ decreases with $r$, decreasing $k_{\mathrm{lim}}$
limits plasma injection to smaller radii. In the $k_{\mathrm{lim}}=0.01$ run,
this supply occurs only up to radii $r/r_{*}\simeq 3$, and the solution shows
the same intermittency observed for $k_{\mathrm{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_{\mathrm{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.

Figure 9: Poynting flux in simulations with volume injection. Values are
normalized to the theoretical value $L_{0}=\mu^{2}\Omega^{4}/c^{3}$. Figure
10: Radial and temporal dependencies of Poynting flux in simulations with
volume injection. a) shows the time-averaged luminosity $\langle L\rangle$ 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
$L_{0}=\mu^{2}\Omega^{4}/c^{3}$, respectively.

An important property of the magnetospheric solution is the integrated
Poynting flux $L(r)$, defined as

$L(r)=\frac{c}{2}\int_{0}^{\pi}(\mathbf{E}\times\mathbf{B})_{r}\
r^{2}\sin\theta\mathrm{d}\theta\ .$ (23)

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}=\mu^{2}\Omega^{4}/c^{3}$, with
$\mu=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_{\mathrm{lim}}=0.1$ simulation converges to a surface
Poynting flux $L_{*}/L_{0}\ll 1$, which is a consequence of the inactivity of
disk-dome solution. On the contrary, the simulations with lower
$k_{\mathrm{lim}}$ have $L_{*}/L_{0}\sim 1$. The Poynting flux remains
approximately constant within the light-cylinder for these runs, and decays
with $r$ for $r>R_{\mathrm{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\leavevmode\nobreak\ T$, which are correlated with the intermittency
of the solution identified above in this section. The time-averaged radial
dependence of the luminosity $\langle L\rangle$ after a stellar period and the
temporal dependence of $L_{*}$ is shown in Fig. 10.

The simulations presented in this section show that the efficiency of the
plasma supply critically determines the global structure of the pulsar
magnetosphere. It is expected that pulsar magnetoshere is in a regime close to
the force-free configuration identified with $k_{\mathrm{lim}}=0.005$ or
lower. However, pair production cannot operate in all regions of the
magnetosphere, in particular at radii comparable to the light-cylinder radius.
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,
that 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\times 10^{3}$, $N_{r}\times N_{\theta}=500^{2}$ and
$\Delta tc/r_{*}=3\times 10^{-3}$.

Figure 11: Magnetospheric solutions obtained with surface injection
proportional to $E_{\parallel}$. a1-2) show the total charge density and
$E_{\parallel}$, respectively for a simulation with $v_{\mathrm{s}}/c=0$ and
b1-2) show the same for a simulation with $v_{\mathrm{s}}/c=0.99$. Solid lines
represent magnetic field lines, and vertical dashed lines show the location of
the light-cylinder.

The first injection criterion is based on the local value of $E_{\parallel}$.
We inject a macro-electron-positron pair in each cell just above the stellar
surface ($r=r_{*}$) that satisfies $E_{\parallel}c/r_{*}\Omega
B_{*}>k_{\mathrm{lim}}$. In this case, we consider a fixed
$k_{\mathrm{lim}}=0.002$ and vary the properties of the injected pairs, namely
their density $n_{\mathrm{s}}=k_{\mathrm{s}}n_{\mathrm{GJ}}$ and poloidal
velocity $v_{\mathrm{s}}$. These pairs are also injected with a toroidal
velocity that matches the local linear velocity of the stellar surface,
$v_{\phi}=\Omega r\sin\theta$.

Despite the large range of injection parameters considered,
$k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=\\{0.2,0.5,1\\}$ and
$v_{\mathrm{s}}/c=\\{0,0.1,0.5,0.99\\}$, the solutions obtained for long
times, $t/T\gtrsim 2$, always converge to the disk-dome solution identified in
Sect. 3.1. Figure 11 shows the charge density $\rho$ and $E_{\parallel}$ of
two runs with $k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=1$ and
$v_{\mathrm{s}}=\\{0,0.99\\}$ after a time $t/T\simeq 4$. After an initial
transient, the system settles to a charge-separated solution and effectively
screens $E_{\parallel}$ at the stellar surface, precluding further injection.

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\leavevmode\nobreak\ n_{\mathrm{GJ}}$, to ensure that
enough plasma exists everywhere to screen the local electric field parallel to
the magnetic field. We emphasize that $n_{\mathrm{GJ}}=\Omega B_{*}/2\pi ec$
is the pole GJ density and not its local value. This criterion allows
injection to occur even if $E_{\parallel}\sim 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 $\sim 100$ m (Ruderman & Sutherland, 1975). We
consider that the injected electron-positron pairs carry a number density
$n_{\mathrm{s}}=k_{\mathrm{s}}n_{\mathrm{GJ}}$ and poloidal velocity
$v_{\mathrm{s}}$.

Figure 12: Magnetospheric solutions obtained with surface injection
proportional to $n_{\mathrm{GJ}}$ with fixed
$k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=0.2$ and varying
$v_{\mathrm{s}}/c$. Figure 13: Radial and temporal dependencies of Poynting
flux in simulations with surface injection proportional to $n_{\mathrm{GJ}}$
with fixed $k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=0.2$ and varying
$v_{\mathrm{s}}/c$. a) shows the time-averaged luminosity $\langle L\rangle$
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 $L_{0}=\mu^{2}\Omega^{4}/c^{3}$, respectively.

In Fig. 12, we show the charge density distribution of the solutions obtained
for a fixed $k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=0.2$ and varying
$v_{\mathrm{s}}$ for a time $t/T=1$. With $v_{\mathrm{s}}=0$, the system
converges to the electrosphere solution. Particles injected at early times
develop a space-charge limited flow, driving $E_{\parallel}$ to zero near the
stellar surface and thus inhibiting freshly injected particles to be pulled
away from or towards the star. For $v_{\mathrm{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_{\mathrm{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_{\mathrm{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.

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_{\mathrm{s}}/c=0)$ produces no spindown luminosity,
and that it increases overall with increasing $v_{\mathrm{s}}$. The same
decrease for $r>R_{\mathrm{LC}}$ observed in Sect. 3.2 is observed here. We
note that the $v_{\mathrm{s}}/c=0.99$ run shows a surface Poynting flux larger
than $L_{0}$, 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_{0}$).

Figure 14: Magnetospheric solutions obtained with surface injection
proportional to $n_{\mathrm{GJ}}$ with fixed $v_{\mathrm{s}}$ and varying
$k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}$.

We have also performed a set of simulations with fixed $v_{\mathrm{s}}/c=0.5$
and varying $k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{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_{\mathrm{s}}=0.5$). This
injection density requirement seems to be critical in the launching of large
density plasma to large radii, in particular along the return current layers,
that connect the surface to the equatorial current sheet.

In summary, some of the parameters used in simulations presented in this
section yield active magnetospheric solutions, with $L_{*}/L_{0}\sim 1$ and a
global configuration similar to the force-free regime. This is 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 regime.

###  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_{\mathrm{s}}$ and $v_{\mathrm{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 $\gamma$-$\gamma$ channel. In
this process, for which the cross-section peaks at around a center of mass
energy $\sim 2\leavevmode\nobreak\ m_{e}c^{2}$, gamma-rays produced via
synchrotron emission and/or inverse Compton scattering in the equatorial
current sheet collide with photons from a low energy bath, producing pairs.
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.

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 $\gamma_{\mathrm{pair}}m_{e}c^{2}$ whenever it
achieves a threshold Lorentz factor $\gamma_{\mathrm{thr}}$. We keep the ratio
$\gamma_{\mathrm{thr}}/\gamma_{\mathrm{pair}}$ constant, and vary the ratio
$\eta\equiv\gamma_{\mathrm{max}}/\gamma_{\mathrm{thr}}$, where
$\gamma_{\mathrm{max}}=e\Phi_{\mathrm{pc}}/m_{e}c^{2}$ is the maximum energy
achievable by the particles in the voltage
$\Phi_{\mathrm{pc}}=B_{*}r_{*}^{3}\Omega^{2}/c^{2}$ induced by the rotating
star across the polar cap. In general,
$\gamma_{\mathrm{pair}}\ll\gamma_{\mathrm{thr}}\ll\gamma_{\mathrm{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.,
2015), have used $\eta\sim 10$ and
$\gamma_{\mathrm{thr}}/\gamma_{\mathrm{pair}}\sim 2$, which severely limits
the efficiency of the pair cascades and the plasma multiplicity. In this
Section, we present simulations with fixed $\gamma_{\mathrm{pair}}=16$ and
$\gamma_{\mathrm{thr}}=25$ and a range of large values of $\eta$. 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 $\eta=5$ we use $N_{r}\times
N_{\theta}=500^{2}$ and $\Delta tc/r_{*}=3\times 10^{-3}$, for
$\eta=\\{25,50\\}$ we use $N_{r}\times N_{\theta}=1000^{2}$ and $\Delta
tc/r_{*}=10^{-3}$ and for $\eta=\\{100,150\\}$ we use $N_{r}\times
N_{\theta}=2000^{2}$ and $\Delta tc/r_{*}=5\times 10^{-4}$.

Figure 15: Magnetospheric solutions obtained with pair production. Panels a-d)
show the total charge density for simulations with $\eta=\\{5,25,50,150\\}$.
Solid lines represent magnetic field lines, and vertical dashed lines show the
location of the light-cylinder.

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 $\theta<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_{\parallel}c/r_{*}\Omega B_{*}>k_{\mathrm{lim}}$, with
$k_{\mathrm{lim}}=0.1$. Each pair is injected at rest and carrying a density
$n_{\mathrm{s}}=k_{\mathrm{s}}E_{\parallel}/er_{*}$, with
$k_{\mathrm{s}}=n_{\mathrm{s}}/n_{\mathrm{GJ}}=0.2$. We stress that in these
conditions, we obtained an electrosphere configuration in simulations without
pair production (see section 3.2).

Figure 16: Radial and temporal dependencies of Poynting flux in simulations
with pair production with varying $\eta$. a) shows the time-averaged
luminosity $\langle L\rangle$ 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 $L_{0}=\mu^{2}\Omega^{4}/c^{3}$,
respectively.

In Figure 15, we show the charge density obtained at a time $t/T\simeq 2$ for
a relevant subset of the simulations performed. We observe a transition from
electrosphere to force-free-like configurations by increasing $\eta$.
Physically, this corresponds to allowing more pairs per particle to be
created, hence increasing the plasma supply of the system. For $\eta=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 $\eta\sim 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 $\eta$, the plasma
produced is not dense enough to confine the equatorial current sheet to a thin
region, and it becomes wide for large distances from the stellar surface. For
$\eta\gtrsim 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\sim$ few $n_{\mathrm{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.

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 $\eta\gtrsim 100$ simulations approaching the force-free
spindown luminosity $L_{0}$ within the light-cylinder. In the equatorial
current sheet, a fraction of $0.3-0.4\leavevmode\nobreak\ L_{*}$ is dissipated
between $r\sim R_{\mathrm{LC}}$ and $r\sim 2\leavevmode\nobreak\
R_{\mathrm{LC}}$ and converted into particle kinetic energy. For all
$\eta<100$ runs, the surface luminosity decreases over time, and we expect
them to eventually converge to the electrosphere solution for $t/T\gg 1$.
However, for $\eta\gtrsim 100$, the surface Poynting flux remains stable over
time.

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_{\parallel}$ screening process resulting from pair
cascades. These fluctuations occur on a temporal scale $\sim r_{*}/\eta c$. We
also observe a quasi-periodic launch of plasma towards the light-cylinder
region along the return current layers with a temporal scale $\sim
0.3-0.5\leavevmode\nobreak\ T$. We show one of these events in Figure 17 for a
simulation with $\eta=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\sim R_{\mathrm{LC}}$, allowing
$E_{\parallel}$ 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 e.g., Figure 17 a1) and
b1). The secondary particles then advect to large radii along the return
current layers, reestablishing $j_{r}$ and effectively screening the
$E_{\parallel}$ responsible for triggering the process — see Figure 17 d1-3).
This process produces a larger fraction of the total pair production events
for $10\lesssim\eta\lesssim 50$. The solutions obtained in this range resemble
that of weak pulsars (Gruzinov, 2015), with screened surface $E_{\parallel}$
but with 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. (2022).

Figure 17:  Cyclic pair production along the return current layers. Columns
labelled 1-3 correspond to different times and rows labelled a-d show the
electron and positron densities, radial current and $E_{\parallel}$,
respectively. Results obtained for $\eta=100$. Solid lines represent magnetic
field lines.

The periodicity of the cyclic behaviour driven by pair production along the
return current layers is $\sim 0.3-0.5\leavevmode\nobreak\ T$. We believe that
this periodicity can depend on the multiplicity from the pair cascade near
$r/r_{*}\sim 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 behaviour is deferred to a future work.

Figure 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 ($\simeq 0.3\leavevmode\nobreak\ r_{*}/c$) in each grid cell.

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 $\eta>5$ — see Figure 18 for an
illustrative example. This occurs because the density supplied from the
stellar surface is insufficient to screen $E_{\parallel}$ 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 $\pm en_{\mathrm{GJ}}$ (Philippov et al., 2015; Chen et al., 2020;
Bransgrove et al., 2022), since general relativity requires a current in this
region $|j_{r}|>en_{\mathrm{GJ}}$ (Beloborodov, 2008; Belyaev & Parfrey, 2016;
Gralla et al., 2016; Torres et al., 2023). Apart from driving this difference
in the time-dependent nature of the polar cap, general relativity is not
expected to play a significant role in the overall magnetospheric
organization.

##  4 Conclusions

In this work, we have presented a systematic study of the different global
regimes of pulsar magnetospheres. Namely, we have 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
considering independently 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) or with
heuristic pair production models (Chen & Beloborodov, 2014; Philippov et al.,
2015, 2015; Chen et al., 2020; Guépin et al., 2020; Bransgrove et al., 2022).

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 e.g., the local electric and magnetic
field geometry. Thus, any volume or surface injection scheme is a substantial
simplification of the highly nonlinear plasma supply from pair cascades in
pulsars. Understanding if and how pair production alone can fill the whole
pulsar magnetosphere is thus crucial, namely 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,
$\eta=\gamma_{\mathrm{max}}/\gamma_{\mathrm{thr}}$, that quantifies the
efficiency of the pair discharges. Our results show that when $\eta\gtrsim
100$, enough 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, $\eta\lesssim
10$, the magnetosphere settles to a fully charge-separated, static solution,
with $E_{\parallel}=0$ near the surface, that produces a negligible Poynting
flux. For $\eta\sim 10-50$, we observe an intermediate solution (Gruzinov,
2015), with a wide equatorial current sheet and with a surface Poynting flux
$50-80\%$ below that expected in the force-free regime.

Our simulations show that the pair production along the return current layers
is key to feed plasma to the light-cylinder region and beyond in near force-
free regimes, in line with the results reported in other works, e.g. Chen &
Beloborodov (2014). We have also identified a time-dependent mechanism similar
to that presented in Chen et al. (2020); Bransgrove et al. (2022), that
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.

The simulations presented here used a very simple heuristic model to describe
pair production in strong magnetic fields. In this work, we have only explored
the role of the parameter $\eta$ on the magnetospheric structure and left the
ratio $\gamma_{\mathrm{thr}}/\gamma_{\mathrm{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
$\gamma_{\mathrm{thr}}/\gamma_{\mathrm{pair}}$ will likely provide even more
abundant pairs to large radii, such that smaller values of $\eta$ may be
enough to set the magnetosphere in a force-free regime. Such study is left for
future work.

The pair production model considered here provides an adequate description of
pair cascades when the curvature photon mean free path is negligible, i.e.,
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 travelling
towards 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 if nonlocal pair production
may play a significant role in e.g., coherent emission processes.

In this work, we have also described a spherical grid suitable to perform
global PIC simulations of pulsar magnetospheres. We have detailed 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
particles’ equations of motion in Cartesian coordinates (Sect. 2.3) and c) a
charge-conserving current deposition scheme (Sect. 2.4) for a non-uniform,
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 are
used to correct for artificial charges in the grid. For each of the numerical
schemes presented here, we have provided comprehensive benchmarks for a
variety of test scenarios. All numerical schemes presented here have been
implemented in the PIC code OSIRIS.

##  5 Acknowledgments

FC, TG, RAF and LOS acknowledge supported by the European Research Council
(ERC-2015-AdG Grant 695088) and FCT (Portugal)— Foundation for Science and
Technology (grant PD/BD/114307/2016, in the framework of the Advanced Program
in Plasma Science and Engineering APPLAuSE, grant PD/00505/2012, and project
no. 2022.02230.PTDC). AC acknowledges support from NSF grants DMS-2235457 and
AST-2308111. AS is supported in part by NSF grant PHY-2206607. We acknowledge
PRACE for granting access to MareNostrum, Barcelona Supercomputing Center
(Spain), where the simulations presented in this work were performed.

## References

  * Beloborodov (2008) Beloborodov, A. M. 2008, Astrophys. J. Lett., 683, L41
  * Belyaev (2015a) Belyaev, M. A. 2015a, Mon. Not. R. Astron. Soc., 449, 2759
  * Belyaev (2015b) Belyaev, M. A. 2015b, New Astron., 36, 37
  * Belyaev & Parfrey (2016) Belyaev, M. A. & Parfrey, K. 2016, Astrophys. J., 830, 119
  * Birdsall & Langdon (1991) Birdsall, C. K. & Langdon, A. B. 1991, Plasma physics via computer simulation (IoP)
  * Boris (1970) Boris, J. 1970, in Proceedings of the Fourth Conference on Numerical Simulation of Plasmas, Naval Research Laboratory, Washington DC, 3
  * Brambilla et al. (2018) Brambilla, G., Kalapotharakos, C., Timokhin, A. N., Harding, A. K., & Kazanas, D. 2018, Astrophys. J., 858, 81
  * Bransgrove et al. (2022) Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2022, Radio Emission and Electric Gaps in Pulsar Magnetospheres
  * Cerutti et al. (2015) Cerutti, B., Philippov, A. A., Parfrey, K., & Spitkovsky, A. 2015, Mon. Not. R. Astron. Soc., 448, 606
  * Chen (2017) Chen, A. Y. 2017, PhD thesis, Columbia University
  * Chen & Beloborodov (2014) Chen, A. Y. & Beloborodov, A. M. 2014, Astrophys. J. Lett., 795, L22
  * Chen et al. (2020) Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, Astrophys. J., 889, 69
  * Cruz et al. (2021a) Cruz, F., Grismayer, T., Chen, A. Y., Spitkovsky, A., & Silva, L. O. 2021a, Astrophys. J. Lett., 919, L4
  * Cruz et al. (2022) Cruz, F., Grismayer, T., Iteanu, S., Tortone, P., & Silva, L. O. 2022, Phys. Plasmas, 29, 052902
  * Cruz et al. (2021b) Cruz, F., Grismayer, T., & Silva, L. O. 2021b, Astrophys. J., 908, 149
  * Dawson (1962) Dawson, J. M. 1962, Phys. Fluids, 5, 445
  * Dawson (1983) Dawson, J. M. 1983, Rev. Mod. Phys., 55, 403
  * Esirkepov (2001) Esirkepov, T. Z. 2001, Comput. Phys. Commun., 135, 144
  * Fonseca et al. (2008) Fonseca, R. A., Martins, S. F., Silva, L. O., et al. 2008, Plasma Phys. Control. Fusion, 50, 124034
  * Fonseca et al. (2002) Fonseca, R. A. et al. 2002, in Computational Science — ICCS 2002, ed. P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, & J. J. Dongarra (Berlin, Heidelberg: Springer Berlin Heidelberg), 342–351
  * Goldreich & Julian (1969) Goldreich, P. & Julian, W. H. 1969, Astrophys. J., 157, 869
  * Gralla et al. (2016) Gralla, S. E., Lupsasca, A., & Philippov, A. 2016, Astrophys. J., 833, 258
  * Grismayer et al. (2016) Grismayer, T., Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Phys. Plasmas, 23, 056706
  * Grismayer et al. (2017) Grismayer, T., Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2017, Phys. Rev. E, 95, 023210
  * Gruzinov (2015) Gruzinov, A. 2015, No-Hair Theorem for Weak Pulsar
  * Guépin et al. (2020) Guépin, C., Cerutti, B., & Kotera, K. 2020, Astron. Astrophys., 635, A138
  * Hakobyan et al. (2023) Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, Astrophys. J., 943, 105
  * Hockney & Eastwood (1988) Hockney, R. W. & Eastwood, J. W. 1988, Computer Simulation Using Particles (CRC Press)
  * Jackson (1976) Jackson, E. A. 1976, Astrophys. J., 206, 831
  * Jackson (1975) Jackson, J. D. 1975, Classical electrodynamics; 2nd ed. (New York, NY: Wiley)
  * Kalapotharakos et al. (2018) Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, Astrophys. J., 857, 44
  * Krause-Polstorff & Michel (1985) Krause-Polstorff, J. & Michel, F. C. 1985, Mon. Not. R. Astron. Soc., 213, 43P
  * Mur (1981) Mur, G. 1981, IEEE Trans. Electromagn. Compat., EMC-23, 377
  * Philippov et al. (2015) Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015, Astrophys. J. Lett., 815, L19
  * Philippov & Spitkovsky (2014) Philippov, A. A. & Spitkovsky, A. 2014, Astrophys. J. Lett., 785, L33
  * Philippov et al. (2015) Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, Astrophys. J. Lett., 801, L19
  * Ruderman & Sutherland (1975) Ruderman, M. A. & Sutherland, P. G. 1975, Astrophys. J., 196, 51
  * Torres et al. (2023) Torres, R., Grismayer, T., Cruz, F., & Silva, L. O. 2023, (in prep.)
  * Touati et al. (2022) Touati, M., Codur, R., Tsung, F., et al. 2022, Plasma Phys. Control. Fusion, 64, 115014
  * Vay (2008) Vay, J.-L. 2008, Phys. Plasmas, 15, 056701
  * Villasenor & Buneman (1992) Villasenor, J. & Buneman, O. 1992, Comput. Phys. Commun., 69, 306
  * Vranic et al. (2016) Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Comput. Phys. Commun., 204, 141
  * Yee (1966) Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302


Applying the Cauchy-$\delta$ inequality with $\delta=k_{1}(e_{0}+K_{M})/2$,
followed by integration, yields

$\mathcal{E}_{T}^{2}\leq\mathcal{E}_{T}^{2}(0)\displaystyle e^{\displaystyle-
k_{1}(K_{M}+e_{0})t}+\cfrac{1}{(k_{1}(e_{0}+K_{M}))^{2}}\bigg{[}k_{1}\lambda^{2}+\cfrac{k_{2}\lambda
e_{0}}{e_{0}+K_{M}}\bigg{]}^{2}.$ (179)

Finally, after scaling and algebraic simplification, it holds that:

$|\mathcal{E}_{T}|\leq|\mathcal{E}_{T}(0)|e^{\displaystyle-
k_{1}(K_{M}+e_{0})t/2}+\lambda\bigg{[}\cfrac{\lambda}{e_{0}+K_{M}}+\nu\bigg{(}\cfrac{e_{0}}{e_{0}+K_{M}}\bigg{)}\bigg{(}\cfrac{K_{M}}{e_{0}+K_{M}}\bigg{)}\bigg{]}.$
(180)

In terms of $T_{z}$, the inequality (180) is

$\displaystyle|\bar{\mathcal{E}}_{T}|$
$\displaystyle\leq|\bar{\mathcal{E}}_{T}(0)|e^{\displaystyle-
s_{0}T_{z}/2\lambda\tilde{\epsilon}}+\cfrac{\eta(1+\nu\tilde{\epsilon})}{1+\eta},\qquad\qquad\;\;\text{if}\;\;e_{0}\leq
s_{0},$ (181a) $\displaystyle|\bar{\mathcal{E}}_{T}|$
$\displaystyle\leq|\bar{\mathcal{E}}_{T}(0)|e^{\displaystyle-
s_{0}T_{z}/2\lambda\tilde{\epsilon}}+\cfrac{\lambda}{e_{0}+K_{M}}+\eta\nu\tilde{\epsilon},\qquad\text{if}\;\;s_{0}<e_{0}$
(181b)

where $\bar{\mathcal{E}}_{T}=\lambda^{-1}\mathcal{E}_{T}$ and
$\tilde{\epsilon}\equiv K_{M}/(e_{0}+K_{M})$. Thus, the tQSSA as used in
practice is valid at very low $e_{0}$ as well as very high $e_{0}$.

## References

  * Crampin et al. [2004] E. J. Crampin, S. Schnell, P. E. McSharry,  Mathematical and computational techniques to deduce complex biochemical reaction mechanisms,  Prog. Biophys. Mol. Biol. 86 (2004) 77–112.
  * Sauro and Kholodenko [2004] H. M. Sauro, B. N. Kholodenko,  Quantitative analysis of signaling networks,  Prog. Biophys. Mol. Biol. 86 (2004) 5–43.
  * Turner et al. [2004] T. Turner, S. Schnell, K. Burrage,  Stochastic approaches for modelling in vivo reactions,  Comput. Biol. Chem. 28 (2004) 165–178.
  * Grima and Leier [2017] R. Grima, A. Leier,  Exact product formation rates for stochastic enzyme kinetics,  J. Phys. Chem. B 121 (2017) 13–23.
  * Bertram and Rubin [2017] R. Bertram, J. E. Rubin,  Multi-timescale systems and fast-slow analysis,  Math. Biosci. 287 (2017) 105–121.
  * Rubin and Wechselberger [2008] J. Rubin, M. Wechselberger,  The selection of mixed-mode oscillations in a Hodgkin-Huxley model with multiple timescales,  Chaos 18 (2008) 015105, 12.
  * Booth and Diniz Behn [2014] V. Booth, C. G. Diniz Behn,  Physiologically-based modeling of sleep-wake regulatory networks,  Math. Biosci. 250 (2014) 54–68.
  * Coullet and Spiegel [1983] P. H. Coullet, E. A. Spiegel,  Amplitude equations for systems with competing instabilities,  SIAM J. Appl. Math. 43 (1983) 776–821.
  * Roberts [1997] A. J. Roberts,  Low-dimensional modelling of dynamics via computer algebra,  Comput. Phys. Commun. 100 (1997) 215–230.
  * Roberts [1989] A. J. Roberts,  The utility of an invariant manifold description of the evolution of a dynamical system,  SIAM J. Math. Anal. 20 (1989) 1447–1458.
  * Crawford and Knobloch [1991] J. D. Crawford, E. Knobloch,  Symmetry and symmetry-breaking bifurcations in fluid dynamics,  Annu. Rev. Fluid Mech. 23 (1991) 341–387.
  * O’Malley [1991] R. E. O’Malley, Jr., Singular perturbation methods for ordinary differential equations, volume 89 of Applied Mathematical Sciences, Springer-Verlag, New York, 1991.
  * Tikhonov [1952] A. Tikhonov,  Systems of differential equations containing small parameters in their derivatives,  Mat. Sb. (N.S.) 31 (1952) 575–586.
  * Gradšteĭn [1953] I. S. Gradšteĭn,  Application of A. M. Lyapunov’s theory of stability to the theory of differential equations with small coefficients in the derivatives,  Mat. Sbornik N. S. 32 (1953) 263–286.
  * Stroberg and Schnell [2016] W. Stroberg, S. Schnell,  On the estimation errors of ${K}_{m}$ and ${V}$ from time-course experiments using the michaelis–menten equation,  Biophys. Chem. 219 (2016) 17–27.
  * Choi et al. [2017] B. Choi, G. A. Rempala, J. K. Kim,  Beyond the Michaelis–Menten equation: Accurate and efficient estimation of enzyme kinetic parameters,  Sci. Rep. 7 (2017) 17018.
  * Roberts [2018] A. J. Roberts, Backwards theory supports modelling via invariant manifolds for non-autonomous dynamical systems, 2018\. arXiv:1804.06998.
  * Heineken et al. [1967] F. G. Heineken, H. M. Tsuchiya, R. Aris,  On the mathematical status of the pseudo-steady state hypothesis of biochemical kinetics,  Math. Biosci. 1 (1967) 95–113.
  * Laidler [1955] K. J. Laidler,  Theory of the transient phase in kinetics with special reference to enzyme systems,  Can. J. Chem. 33 (1955) 1614–1624.
  * Reich and Sel’kov [1974] J. G. Reich, E. E. Sel’kov,  Mathematical analysis of metabolic networks,  FEBS Lett. 40 (1974) S112–S118.
  * Schauer and Heinrich [1979] M. Schauer, R. Heinrich,  Analysis of the quasi-steady-state approximation for an enzymatic one-substrate reaction,  J. Theor. Biol. 79 (1979) 425–442.
  * Segel [1988] L. A. Segel,  On the validity of the steady state assumption of enzyme kinetics,  Bull. Math. Biol. 50 (1988) 579–593.
  * Burke et al. [1990] M. A. Burke, P. K. Maini, J. D. Murray,  On the kinetics of suicide substrates,  Biophys. Chem. 37 (1990) 81–90.
  * Burke et al. [1993] M. A. Burke, P. K. Maini, J. D. Murray,  Suicide substrate reaction-diffusion equations: Varying the source,  IMA J. Math. Appl. Med. Biol. 10 (1993) 97–114.
  * Schnell and Mendoza [2000a] S. Schnell, C. Mendoza,  Enzyme kinetics of multiple alternative substrates,  J. Math. Chem. 27 (2000a) 155–170.
  * Schnell and Mendoza [2000b] S. Schnell, C. Mendoza,  Time-dependent closed form solutions for fully competitive enzyme reactions,  Bull. Math. Biol. 62 (2000b) 321–336.
  * Eilertsen et al. [2018a] J. Eilertsen, W. Stroberg, S. Schnell,  A theory of reactant-stationary kinetics for a mechanism of zymogen activation,  Biophys. Chem. 242 (2018a) 34–44.
  * Eilertsen et al. [2018b] J. Eilertsen, W. Stroberg, S. Schnell,  Phase-plane geometries in coupled enzyme assays,  Math. Biosci. 306 (2018b) 126–135.
  * Eilertsen et al. [2019] J. Eilertsen, W. Stroberg, S. Schnell,  Characteristic, completion or matching timescales? An analysis of temporary boundaries in enzyme kinetics,  J. Theoret. Biol. 481 (2019) 28–43.
  * Segel and Slemrod [1989] L. A. Segel, M. Slemrod,  The quasi-steady-state assumption: A case study in perturbation,  SIAM Rev. 31 (1989) 446–477.
  * Schnell and Maini [2000] S. Schnell, P. K. Maini,  Enzyme kinetics at high enzyme concentration,  Bull. Math. Biol. 62 (2000) 483–499.
  * Borghans et al. [1996] J. A. Borghans, R. J. D. Boer, L. A. Segel,  Extending the quasi-steady state approximation by changing variables,  Bull. Math. Biol. 58 (1996) 43–63.
  * Schnell and Maini [2002] S. Schnell, P. K. Maini,  Enzyme kinetics far from the standard quasi-steady-state and equilibrium approximations,  Math. Comput. Modelling 35 (2002) 137–144.
  * Tzafriri [2003] A. R. Tzafriri,  Michaelis–Menten kinetics at high enzyme concentrations,  Bull. Math. Biol. 65 (2003) 1111–1129.
  * Goeke et al. [2015] A. Goeke, S. Walcher, E. Zerz,  Determining “small parameters” for quasi-steady state,  J. Differ. Equations. 259 (2015) 1149–1180.
  * Goeke and Walcher [2014] A. Goeke, S. Walcher,  A constructive approach to quasi-steady state reductions,  J. Math. Chem. 52 (2014) 2596–2626.
  * Noethen and Walcher [2011] L. Noethen, S. Walcher,  Tikhonov’s theorem and quasi-steady state,  Dis. and Cont. Dyn. Sys. B. 16 (2011) 945–961.
  * Fenichel [1979] N. Fenichel,  Geometric singular perturbation theory for ordinary differential equations,  J. Differ. Equations 31 (1979) 53–98.
  * Hanson and Schnell [2008] S. M. Hanson, S. Schnell,  Reactant stationary approximation in enzyme kinetics,  J. Phys. Chem. A 112 (2008) 8654–8658.
  * Nguyen and Fraser [1989] A. H. Nguyen, S. Fraser,  Geometrical picture of reaction in enzyme kinetics,  J. Chem. Phys. 91 (1989) 186–193.
  * Roussel and Fraser [1990] M. R. Roussel, S. J. Fraser,  Geometry of the steady-state approximation: Perturbation and accelerated convergence methods,  J. Chem. Phys. 93 (1990) 1072–1081.
  * Roussel and Fraser [1991] M. R. Roussel, S. J. Fraser,  Accurate steady-state approximations: Implications for kinetics experiments and mechanism,  J. Phys. Chem. 95 (1991) 8762–8770.
  * Goeke et al. [2012] A. Goeke, C. Schilli, S. Walcher, E. Zerz,  Computing quasi-steady state reductions,  J. Math. Chem. 50 (2012) 1495–1513.
  * Pedersen et al. [2006] M. G. Pedersen, A. M. Bersani, E. Bersani,  The total Quasi-Steady-State Approximation for fully competitive enzyme reactions,  Bull. Math. Biol. 69 (2006) 433â€“–457.
  * Tzafriri and Edelman [2007] A. R. Tzafriri, E. R. Edelman,  Quasi-steady-state kinetics at enzyme and substrate concentrations in excess of the Michaelis-Menten constant,  J. Theoret. Biol. 245 (2007) 737–748.
  * Dingee and Anton [2008] J. W. Dingee, A. B. Anton,  A new perturbation solution to the Michaelis-Menten problem,  AIChE J. 54 (2008) 1344–1357.
  * Dell’Acqua and Bersani [2012] G. Dell’Acqua, A. M. Bersani,  A perturbation solution of Michaelis-Menten kinetics in a “total” framework,  J. Math. Chem. 50 (2012) 1136–1148.
  * Dhatt and Bhattacharyya [2013] S. Dhatt, K. Bhattacharyya,  Single-substrate enzyme kinetics: the quasi-steady-state approximation and beyond,  J. Math. Chem. 51 (2013) 1467–1477.
  * Pedersen et al. [2008] M. G. Pedersen, A. M. Bersani, E. Bersani, G. Cortese,  The total quasi-steady-state approximation for complex enzyme reactions,  Math. Comput. Simulation 79 (2008) 1010–1019.
  * Bersani et al. [2015] A. M. Bersani, E. Bersani, G. Dell’Acqua, M. G. Pedersen,  New trends and perspectives in nonlinear intracellular dynamics: one century from Michaelis-Menten paper,  Contin. Mech. Thermodyn. 27 (2015) 659–684. URL: https://doi.org/10.1007/s00161-014-0367-4. doi:10.1007/s00161-014-0367-4.
  * Battelli and Lazzari [1986] F. Battelli, C. Lazzari,  Singular perturbation theory for open enzyme reaction networks,  IMA J. Math. Appl. Med. Biol. 3 (1986) 41–51.
  * Verhulst [2007] F. Verhulst,  Singular perturbation methods for slow-fast dynamics,  Nonlinear Dynam. 50 (2007) 747–753.
  * Kuehn [2015] C. Kuehn, Multiple time scale dynamics, Springer, Cham, 2015.
  * Hek [2010] G. Hek,  Geometric singular perturbation theory in biological practice,  J. Math. Biol. 60 (2010) 347–386.
  * Roussel [2019] M. R. Roussel,  Heineken, Tsuchiya and Aris on the mathematical status of the pseudo-steady state hypothesis: A classic from Volume 1 of Mathematical Biosciences,  Math. Biosci. 318 (2019) 108274, 8.
  * Boulier et al. [2009] F. Boulier, F. Lemaire, A. Sedoglavic, A. Ürgüplü,  Towards an automated reduction method for polynomial ODE models of biochemical reaction systems,  Math. Comput. Sci. 2 (2009) 443–464.
  * Goeke and Walcher [2013] A. Goeke, S. Walcher,  Quasi-steady state: searching for and utilizing small parameters,  in: Recent trends in dynamical systems, volume 35 of Springer Proc. Math. Stat., Springer, Basel, 2013, pp. 153–178.
  * Goeke et al. [2017] A. Goeke, S. Walcher, E. Zerz,  Classical quasi-steady state reduction—A mathematical characterization,  Phys. D 345 (2017) 11–26.
  * Roberts [1989] A. J. Roberts,  Appropriate initial conditions for asymptotic descriptions of the long term evolution of dynamical systems,  J. Austral. Math. Soc. Ser. B 31 (1989) 48–75.
  * Lax and Walcher [2020] C. Lax, S. Walcher,  Singular perturbations and scaling,  Discrete Contin. Dyn. Syst. Ser. B 25 (2020) 1–29.
  * Roussel [1994] M. R. Roussel, A rigorous approach to steady-state kinetics applied to simple enzyme mechanisms, ProQuest LLC, Ann Arbor, MI, 1994. Thesis (Ph.D.)–University of Toronto (Canada).
  * Fraser [1988] S. J. Fraser,  The steady state and equilibrium approximations: A geometrical picture,  J. Chem. Phys. 88 (1988) 4732–4738.
  * Pavliotis and Stuart [2008] G. A. Pavliotis, A. M. Stuart, Multiscale Methods: Averaging and homogenization, Springer, New York, 2008.