© C. Guépin et al. 2020

1. Introduction

Pulsars are rapidly rotating and highly magnetized neutron stars, that have been detected across the entire electromagnetic spectrum, from radio to gamma rays (see e.g. Abdo et al. 2013 for a Fermi-LAT catalog of gamma-ray pulsars and Aliu et al. 2008; VERITAS Collaboration 2011; Aleksić et al. 2012 for the detection of the Crab pulsar above 100 GeV). Their high-energy emissions have been associated with the radiation of accelerated leptons (e.g Arons 1983; Cheng et al. 1986; Romani 1996; Muslimov & Harding 2003 for acceleration by unscreened electric fields close to the neutron star surface).

A description is emerging of the structure of the magnetosphere from first principles accounting for the feedback of particles on the electromagnetic fields (e.g. Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Philippov et al. 2015; Cerutti et al. 2015; Belyaev 2015; Brambilla et al. 2018). The various particle interactions occurring in the magnetosphere and in particular the production of pairs (e.g. Daugherty & Harding 1982; Gurevich & Istomin 1985; Zhang & Harding 2000; Medin & Lai 2010; Timokhin 2010; Timokhin & Arons 2013) are processes that remain to be fully understood and self-consistently implemented into large-scale systems. This could have a critical influence on our understanding of pair multiplicities and high-energy emissions of pulsars, and might help to improve the models for energy dissipation and spin down. Another fundamental question is related to the nature of the wind around pulsars. The location of the energy dissipation where the Poynting flux is dissipated into particle kinetic energy is still to be clearly identified (e.g. Coroniti 1990; Kirk & Skjæraasen 2003; Komissarov 2013; Porth et al. 2013; Cerutti & Philippov 2017). Finally, the mechanisms for cosmic-ray acceleration in pulsars, studied for instance in Venkatesan et al. (1997), Blasi et al. (2000), Arons (2003), Fang et al. (2012, 2013), Lemoine et al. (2015) and Kotera et al. (2015), should be modelled from first principles in order to more precisely infer the contribution of these sources to the observed cosmic-ray flux. Many source categories have already been considered as potential cosmic-ray accelerators, such as active galactic nuclei, gamma-ray bursts, tidal disruption events, and pulsars. Despite the advent of time-domain and multimessenger astronomy, together with significant progress in source emission modelling, it is still difficult to discriminate between the various source scenarios.

In this context, a kinetic approach is required. Interestingly, most of the previous studies focused on magnetospheres filled with a plasma of electrons and positrons, without ion injection, and the injection of ions has only recently been considered (Chen & Beloborodov 2014; Philippov & Spitkovsky 2018). It is therefore timely to study the fate of protons in pulsar magnetospheres. Chen & Beloborodov (2014) study an axisymmetric pulsar magnetosphere, where ions and electrons are injected from the neutron star surface and pairs can be produced. These latter authors notice the different trajectories of pairs and ions as well as the acceleration and escape of ions in different configurations. Philippov & Spitkovsky (2018) consider a similar setup with an oblique rotator. In their simulations, ions have the same mass as positrons but they do not suffer radiative energy losses. They notice the acceleration of ions in the current sheet, mostly at the Y-point.

In the present study, we focus on the two following fundamental questions: (i) what is the maximum energy achievable for the ions, especially for the ones escaping the magnetosphere, and (ii) to what level can they contribute to the observed high- and ultra-high-energy cosmic-ray fluxes? To this end, we perform particle-in-cell (PIC) simulations of the aligned pulsar magnetosphere. With this series of numerical experiments, we aim to explore the transition between a charge-separated magnetosphere, or “electrosphere” in the following, and a force-free magnetosphere, by changing the yield of pair production, and assessing its impact on particle acceleration and escape. We describe the theoretical and numerical setup in Sect. 2. The structure of the magnetosphere is described in Sect. 3. The questions of proton acceleration in the simulations and how our results scale up to realistic pulsar parameters are addressed in Sect. 4, followed by a discussion and our conclusions in Sect. 5.

2. Simulating a pulsar magnetosphere

We make use of the PIC code ZELTRON (Cerutti et al. 2013) in its 2D axisymmetric version, with a non-uniform spherical grid, which is well suited to the study of aligned rotators.

2.1. Electromagnetic fields

In the following, r, θ, and ϕ are the usual spherical coordinates. The initial setup of our simulations is a perfectly conducting neutron star in a vacuum, with a magnetic dipole anchored at its surface

$$\begin{array}{cc}\hfill B_{\mathrm{r}}(r,\theta )& =B_{\star }{R}_{\star }^{3}cos\theta /r^3,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill B_{\theta }(r,\theta )& =B_{\star }{R}_{\star }^{3}sin\theta /2r^3,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill B_{\varphi }(r,\theta )& =0,\hfill \end{array}$$(3)

where R_⋆ is the radius of the neutron star, θ is the angle from the rotation axis, and B_⋆ is the polar magnetic field. For a perfect conductor rotating at angular velocity Ω, E′ = E + (Ω × r)×B/c = 0 in the co-rotating frame, where E and B are the electric and magnetic fields in the observer frame, and Ω is along the rotation axis. This allows us to estimate the electric field inside the star $(E_{\mathrm{r}}^{\mathrm{int}},E_{\theta }^{\mathrm{int}},E_{\varphi }^{\mathrm{int}})=(rsin\theta /R_{\mathrm{LC}})(B_{\theta },-B_{\mathrm{r}},0)$, where R_LC = c/Ω is the light cylinder radius defined as the distance at which the corotating speed reaches the speed of light. At t = 0, the rotation of the neutron star is forced by imposing at its surface the poloidal electric field induced by the rotation of a perfect conductor. The radial electric field can be discontinuous for a non-zero surface charge density. The outer boundary condition is defined by an absorbing layer in order to mimic an open boundary with no information coming back inwards (Birdsall & Langdon 1991; Cerutti et al. 2015). Apart from these boundary conditions, there are no constraints on the external electric field, which evolves self-consistently during the simulation.

2.2. Particle extraction

We note that at the surface of the star, the electric field is of the order E_⋆ ∼ B_⋆R_⋆/R_LC ∼ 10⁸ statV cm⁻¹ for a millisecond pulsar with B_⋆ = 10⁹ G. Due to this high electric field, charged particles can be extracted from the neutron star surface. In our work, we neglect the molecular or gravitational attraction (Pétri 2016, see however Ruderman & Sutherland 1975). We consider three particle species: electrons, positrons, and protons. Electrons and protons are extracted from the surface and positrons are created through a pair production process. In order to avoid overinjection, particles can be extracted when the local charge density does not exceed the local Goldreich–Julian (Goldreich & Julian 1969, GJ) charge density ρ_GJ. At the neutron star surface, for a dipole magnetic field and a rotation around the vertical axis, the GJ charge density is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{GJ}}=\frac{-B_{\star }(3{cos}^2\theta -1)}{4\pi R_{\mathrm{LC}}[1-{(R_{\star }sin\theta /R_{\mathrm{LC}})}^2]}·\end{array}$$(4)

The denominator adds a relativistic correction due to the modification of the magnetic field structure by currents, and is as small as R_⋆/R_LC ≃ 0.2 for a millisecond pulsar. Therefore, for electrons, the GJ number density at the surface of the neutron star reads n_GJ = B_⋆(3cos²θ − 1)/4πR_LCe.

2.3. Energy losses by radiation

In our simulations, the motion of a particle is governed by the Abraham–Lorentz–Dirac equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}=q(\mathit{E}+\mathit{\beta }\times \mathit{B})+\mathit{g},\end{array}$$(5)

where p = γmv is the particle momentum, $\gamma =1/\sqrt{1-{\beta }^2}$ the particle Lorentz factor, v = βc the particle 3-velocity, m the particle mass, and q the particle electric charge. The first right-hand side term is the usual Lorentz force and g is the radiation reaction force due to the radiation of accelerated particles given by the Landau–Lifshitz formula in the framework of classical electrodynamics (Landau & Lifshitz 1975):

$$\begin{array}{cc}\hfill \mathit{g}=& \frac{2}{3}\frac{q^4}{m^2{c}^4}[(\mathit{E}+\mathit{\beta }\times \mathit{B})\times \mathit{B}+(\mathit{\beta }·\mathit{E})\mathit{E}]\hfill \\ \hfill & -\frac{2}{3}\frac{q^4{\gamma }^2}{m^2{c}^4}[{(\mathit{E}+\mathit{\beta }\times \mathit{B})}^2-{(\mathit{\beta }·\mathit{E})}^2]\mathit{\beta },\hfill \end{array}$$(6)

where the terms containing the time derivative of the fields are neglected (Tamburini et al. 2010). This simplified view suffices for the current exploration study, as it accounts for synchrotron, synchro-curvature, and curvature regimes; but we note that detailed models that have recently been developed can lead to deviations from the standard curvature and synchrotron radiation spectra in the strong field regime (e.g. Voisin et al. 2017).

2.4. Pair production

The configuration of the magnetosphere, especially the plasma density and the existence of gaps, relies primarily on the production of electron and positron pairs. The pairs are also thought to contribute to the high-energy radiation of Pulsar Wind Nebulae (PWNe) through synchrotron and inverse-Compton radiation. A precise understanding of the pair production process is therefore critical for the modelling of pulsars. However, the amount of pair production in pulsar magnetospheres is poorly constrained. The pairs are thought to be mainly produced in the polar cap regions (Ruderman & Sutherland 1975) by the conversion of high-energy gamma rays into pairs in strong magnetic fields (i.e. B ≳ 10¹¹ G) and the subsequent development of a pair cascade. In the classical model, gamma rays are initially produced through curvature radiation. In the outer gaps, the interaction of gamma-ray photons with X-ray photons from the neutron star surface could also make a significant contribution to the production of pairs in the pulsar magnetospheres (Cheng et al. 2000). The pair multiplicity κ = (n⁺ + n⁻)/2n_GJ, which describes the number of electron and positron pairs produced by each primary particle, is a poorly constrained parameter that could range between 1 and 10⁷. From observations and PWNe emission models, the multiplicity has been estimated to be about 10⁵ − 10⁷ for the Crab PWN and 10⁵ for the Vela PWN (e.g. de Jager 2007; Bucciantini et al. 2011). However, recent theoretical predictions limit the pair multiplicity to about a few hundreds of thousands (Timokhin & Harding 2019) achieved for magnetic fields of 4 × 10¹² ≲ B ≲ 10¹³ G and hot neutron star surfaces T ≳ 10⁶ K, which questions the existing models of PWNe emissions requiring very high pair multiplicities.

Electron–positron pair plasma generation is the subject of active research (Timokhin & Arons 2013; Chen & Beloborodov 2014). In this study, we do not aim to model the full pair cascades; a simplified treatment is adopted, as described in Philippov et al. (2015). Pairs are directly produced at the location of the parent lepton if its Lorentz factor exceeds the threshold γ >  γ_min, pp and the produced pairs have a Lorentz factor γ_f ∼ f_γγ_i, which is a fraction f_γ = 0.1 of the Lorentz factor of the parent particle γ_i. This fraction is chosen for numerical reasons, namely to conserve a reasonable separation of scales. The threshold γ_min, pp = f_ppγ_0, e is a fraction f_pp of the maximum Lorentz factor of pairs γ_0, e = eΦ₀/m_ec² obtained by the acceleration of a particle through the full vacuum potential drop from pole to equator:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_0=-{\displaystyle {\int }_0^{\pi /2}\mathrm{d}\theta R_{\star }{E}_{\theta }(R_{\star })=\frac{R_{\star }^2{B}_{\star }}{2R_{\mathrm{LC}}}·}\end{array}$$(7)

For simplicity, the threshold is constant throughout the magnetosphere, and does not depend on the curvature of the magnetic field lines. This corresponds to a maximization of pair production as all the regions where leptons can be accelerated to sufficiently high energies are active pair-producing regions. Thus, pair production can take place for instance in the equatorial current sheet (Lyubarskii 1996) and not only in the polar cap regions near the surface of the neutron star. This simplified approach allows us to explore various regimes of the magnetosphere by adjusting the parameter f_pp, without entering into a detailed modelling of the radiative backgrounds that could significantly contribute to the production of pairs, and would require additional parameters. The connection between the implemented parameter f_pp governing the pair production in the entire magnetosphere, and the pair multiplicity κ at the polar cap is intricate. A straightforward comparison can be made by computing κ at the poles directly in the simulation and comparing the results with theoretical and observed values. Moreover, the amount of pair production in realistic populations of millisecond or newborn pulsars remains to be explored. The modelling of gamma-ray emissions could also help to establish a clearer link between these quantities.

Several additional effects such as photohadronic interactions could impact the particle motion and contribute to energy losses and pair production. For instance, we do not account for Bethe–Heitler processes that could contribute to the production of pairs. Inverse-Compton scattering could also lead to significant energy losses. These effects could be included in a future study.

2.5. Simulation features

The parameters of our simulations and notations are summarised in Table 1. In order to maintain acceptable computation costs, the radius R_⋆, the magnetic field B_⋆, and the mass ratio m_r  =  m_p/m_e are scaled down with respect to realistic values R_⋆  ∼  10⁶ cm, B_⋆  ∼  10⁸ − 10¹⁵ G and m_p/m_e ≃ 1836. In our simulations, R_⋆  =  10² cm, B_⋆ = 1.1 × 10⁵ G and m_r  =  18.36. However, several typical scales are preserved. For millisecond pulsars, we have $R_{\mathrm{LC}}/R_{\star }=cP/2\pi R_{\star }\sim 5P_{-3}{R}_{\star ,6}^{-1}$. We adopt this typical value of R_LC/R_⋆ in the simulations. Considering $n_{\text{e}}=n_{\text{GJ}}^{\star }$, the plasma skin depth of electrons at the star surface is $d_{\mathrm{e}}^{\star }=c/{(4\pi n_{\mathrm{e}}{e}^2/m_{\mathrm{e}})}^{1/2}\simeq 2\mathrm{cm}{B}_9^{-1/2}{P}_{-3}^{1/2}$ for millisecond pulsars. This value is similar in our simulations, and is resolved by several computational cells. The gyroradius of electrons r_g = γmcv_⊥/eB, where v_⊥ is the speed perpendicular to the magnetic field direction, is not resolved in our simulations for γv_⊥/c ∼ 1, that is, before they are accelerated close to the neutron star surface. This should not have an impact on the simulations, as electrons efficiently lose their perpendicular energy due to strong synchrotron radiation close to the neutron star surface. The radiation reaction force is amplified (within the limits of time resolution) by considering an effective magnetic field B_eff ∼ 10⁹ G in order to reduce the synchrotron cooling time.

At the beginning of the simulation, the magnetosphere is empty. The system is then strongly perturbed due to the sudden induced electric field, and the subsequent injection of particles and reconfiguration of electromagnetic fields. In order to capture the system behavior when the stationary regime is established, we evolve the system during at least five rotation periods. We check that the simulations have reached a steady state by looking at the temporal evolution of several key quantities, such as the total energy, the radial Poynting flux, and the electric charge of the neutron star. A non-uniform spherical grid of 2464 × 2464 cells is used; these are logarithmically spaced in r between r = R_⋆ and r = 5R_LC, and uniformly in θ between θ = 0 and θ = π/2.

3. Structure of the magnetosphere

3.1. Plasma density

We study the impact of changing the amount of pair production on the structure of the magnetosphere. Two extreme regimes of the magnetosphere have been described in the literature: the electrosphere, or disc-dome configuration (Krause-Polstorff & Michel 1985; Smith et al. 2001; Pétri et al. 2002; Spitkovsky & Arons 2002), and the force-free configuration (Goldreich & Julian 1969; Contopoulos et al. 1999; Timokhin 2006). In the first configuration, the production of pairs is not sufficient to completely populate the magnetosphere with plasma. Charges are separated in different regions, with a dome of negative charges and a disc of positive charges (if Ω ⋅ B >  0). In the second configuration, interesting features appear such as the transition between closed field lines and open field lines, the direct volume currents at the poles and the return currents along the last closed field line (the equatorial current sheet). The guideline studies illustrated in Figs. 1 and 2 mimic the transition between these extreme regimes. These examples are obtained for a high production of pairs (f_pp = 0.01) which mimics a force-free case and for a low production of pairs (f_pp = 0.05 and f_pp = 0.1) tending towards an electrosphere configuration, when a steady state is reached.

Fig. 1. Density maps for electrons, positrons, and protons (left to right), for fpp = 0.01, fpp = 0.05 and fpp = 0.1 (top to bottom), as a function of r/RLC and θ, for t = 5P. The densities are normalised by ${(r/R_{\star })}^{-2}{n}_{\text{GJ}}^{\star }$, where $n_{\text{GJ}}^{\star }=B_{\star }/2\pi R_{\text{LC}}e$ is the polar GJ number density. Solid black lines are the magnetic field lines. The dashed black line indicates the distance from the rotation axis r sin θ = RLC and the grey semi-disc represents the neutron star. The densities are in log scale in order to enhance the contrast.

Fig. 2. Same as Fig. 1 but for the zoomed-in simulation domain r <  RLC.

Low values of f_pp lead to a strong production of pairs and thus allow to study magnetospheres close to the force-free regime. We note that small gaps separate the polar flows and the current sheet. The magnetic field, initially in a dipolar configuration, is strongly affected by the dense plasma outflow. A closed magnetic field line region is maintained at low latitudes below the light-cylinder radius, whereas magnetic field lines open up at high latitudes, which is similar to the configuration in the force-free regime (e.g. Contopoulos et al. 1999; Timokhin 2007). Theoretically, the magnetic field lines have been predicted to present a Y-shape at the point where the last closed field line intersects the equatorial plane (called the Y-null point or Y-point), which seems to be the case in these simulations. We note however that in our simulations, the Y-point is located slightly below the light cylinder radius. We note that in the perfect steady state force-free MHD view, the Y-point is located at the light cylinder. However, in general, it does not have to be so and its location depends on the plasma content and magnetisation, as mentioned for instance in Belyaev (2015). We also observe that the position of the Y-point is highly time-dependent; it has a breathing motion around the light-cylinder radius which reflects the intermittent nature of reconnection. This phenomenon is emphasised by the production of plasma overdensities trapped in magnetic loops, i.e. plasmoids, within the equatorial current sheet (see top panels in Fig. 1).

Our simulation with f_pp = 0.01 is characterised by number densities of electrons in the polar regions of the order of the polar GJ number density $n_{\text{GJ}}^{\star }=B_{\star }/2\pi R_{\text{LC}}e$ multiplied by the distance scaling factor (R_⋆/r)². Therefore, only a small number of pairs are produced at the poles in agreement with previous studies (Chen & Beloborodov 2014; Philippov et al. 2015). In the equatorial region (the current sheet), the number density of pairs can reach higher values locally, above $10n_{\text{GJ}}^{\star }$, which is a signature of strong pair production. Protons are propagating in the equatorial region, with number densities around 1 to 10% of ${(R_{\star }/r)}^2{n}_{\text{GJ}}^{\star }$. Below the light cylinder, escaping protons flow along the last closed field line (the separatrices). At larger radii, protons form an equatorial flow, which is not confined inside the current sheet formed by the pairs, because of their larger gyroradii. We note that at the Y-point, the pair and proton skin depths are respectively $d_{\mathrm{e}}/R_{\star }=\sqrt{{\gamma }_{\mathrm{pc},\mathrm{e}}{m}_{\mathrm{e}}{c}^2/4\pi n_{\mathrm{e}}{e}^2}/R_{\star }\simeq 0.1$ and $d_{\mathrm{p}}/R_{\star }=\sqrt{{\gamma }_{\mathrm{pc},\mathrm{p}}{m}_{\mathrm{p}}{c}^2/4\pi n_{\mathrm{p}}{e}^2}/R_{\star }\simeq 4$, considering the polar cap Lorentz factors of pairs and protons and the number densities at the Y-point from the simulations. Therefore, protons are sensitive to larger field structures.

For high values of f_pp, the production of pairs is almost absent and allows us to study magnetospheres close to the electrosphere configuration. The simulation with f_pp = 0.05 illustrates the transition between the force-free and electrosphere regimes, and the simulation with f_pp = 0.1 illustrates the electrosphere regime. The magnetic field lines that open up at the beginning of the simulation because of the transitory dense plasma outflow tend to return to a dipolar configuration because of the low plasma density in the stationary regime. High electric fields contribute to maintaining the equatorial flow of protons, but the subsequent current is not sufficiently large to modify the dipolar structure of the magnetic field. The few open field lines are anchored to the star near the poles, where the electrons are extracted. In our simulation with f_pp = 0.1, there are no positrons in the magnetosphere. Electrons are essentially confined in the polar regions and are characterised by number densities of around 10% of ${(R_{\star }/r)}^3{n}_{\text{GJ}}^{\star }$, with higher number densities close to the star surface (as expected $\sim n_{\text{GJ}}^{\star }$ at the neutron star surface) and in high-latitude elongated regions. In these regions, it appears that electrons are trapped and are going back and forth before escaping or falling back to the star surface. Large gaps of densities separate the bulk of electrons and protons. A high density of protons is confined near the neutron star surface, with $n~n_{\text{GJ}}^{\star }$. Low number densities of protons, that is, below ${10}^{-3}{(R_{\star }/r)}^2{n}_{\text{GJ}}^{\star }$, propagate in the equatorial region and along the separation region between the bulk of electrons and the gaps. Due to the structure of the magnetic and electric fields, these protons escape from the disc of proton and swirl around the neutron star, with a large radial velocity.

3.2. Pair multiplicity at the pole and Y-point

The plasma densities estimated at the pole and at the expected location of the Y-point are illustrated in Fig. 3. The density of pairs divided by the typical local GJ density n_GJ = |B ⋅ Ω|/2πec (where the correction due to the modification of the magnetic field structure by currents is not accounted for) gives a local estimate of the pair multiplicity κ ∼ (n⁺ + n⁻)/2n_GJ, where n⁺ and n⁻ are respectively the densities of positrons and electrons. As noted previously, pair production is sub-dominant by several orders at the poles and occurs predominantly in the current sheet, which is consistent with several recent studies (Chen & Beloborodov 2014; Philippov et al. 2015). At the pole, f_pp ≤ 0.01 leads to a pair multiplicity of κ ∼ 1, which decreases slightly with increasing f_pp. At the Y-point, a high production of pairs leads to high multiplicities, for instance κ ∼ 10³ for f_pp = 0.01. When the production of pairs decreases, the pair multiplicity drops significantly, below κ = 1 for f_pp ≥ 0.05.

Fig. 3. Plasma densities (n+ + n−)/nGJ as a function of fpp at the pole and at the Y-point (respectively blue and orange lines).

3.3. Charge and current densities

The charge densities and radial currents are illustrated in Fig. 5 and have been normalised by the typical GJ charge and current densities at the poles, ${\rho }_{\text{GJ}}^{\star }=B_{\star }/2\pi R_{\text{LC}}$ and $J_{\text{GJ}}^{\star }=c{B}_{\star }/2\pi R_{\text{LC}}$, respectively, and multiplied by (r/R_⋆)². They present interesting features, due to the mixing of particle species. For f_pp = 0.01, the poles are dominated by negative charge densities, which carry a negative radial current out of the polar caps. The equator is mostly dominated by positive charge densities, and a positive current density. Just above the last closed field line, a small region is dominated by negative charge densities. The corresponding radial current shows that these negative charges are the main contributors to the return current (they carry a positive radial current), closing on the polar caps. The closed field line region is dominated by positive charge densities, which do not seem to contribute much to the radial current. For f_pp = 0.05 the situation is simpler, with negative charge densities at the poles and positive charge densities at the equator. The charge densities are smaller than for f_pp = 0.01, except in elongated regions around the poles for negative charge densities, and in a disc close to the neutron star surface for positive charge densities. Interestingly, only high latitudes seem to contribute to radial currents, with a small return current directly next to the negative current. Protons and electrons in this region therefore seem to contribute more to the return current than protons in the equatorial region.

Moreover, we illustrate in Fig. 4 the impact of pair production on the radial current density at r = 2R_LC averaged over θ and normalised by ${(R_{\star }/r)}^2{J}_{\text{GJ}}^{\star }$. As expected, the total radial current density is close to zero, that is, the star does not charge up. Electrons and positrons carry respectively most of the negative and positive radial currents, and protons only contribute marginally to the radial current. The radial current densities associated with electrons and positrons decrease for decreasing pair production, with a maximum value of $J_{\text{r}}(r)~0.1{(R_{\star }/r)}^2{J}_{\text{GJ}}^{\star }$ at r = 2R_LC, for f_pp = 0.01.

Fig. 4. Radial current density for r = 2RLC, averaged over θ and normalised by the polar GJ current density ${(R_{\star }/r)}^2{J}_{\text{GJ}}^{\star }$, as a function of fpp. We show the total radial current density (circles) and the contributions of electrons (crosses), positrons (plus signs), and protons (triangles).

Fig. 5. Charge density maps (top) and radial current maps (bottom) for a high production of pairs fpp = 0.01 (left) and a low production of pairs fpp = 0.05 (right), for r ≤ RLC. ${\rho }_{\text{GJ}}^{\star }$ and $J_{\text{GJ}}^{\star }$ are the GJ charge and current densities at the poles.

4. Particle acceleration and energy dissipation

A magnetized rotating conductor develops a potential difference between the pole and the equator. Particles that experience all or a fraction of the voltage drop can get accelerated through unipolar induction. This is the case for rotating and magnetized neutron stars, that are considered as perfect conductors in our model. As shown in Sect. 2.4, particles can be accelerated up to γ₀ = eΦ₀/mc² where ${\mathrm{\Phi }}_0=B_{\star }{R}_{\star }^2/2R_{\text{LC}}$ (see Eq. (7)) if they experience the full vacuum potential drop. A typical fraction of this full vacuum potential drop is given by the potential drop across the polar cap, the surface of the neutron star on which open field lines are anchored. As the typical polar cap angle is sin² θ_pc ∼ R_⋆/R_LC, it yields ${\mathrm{\Phi }}_{\text{pc}}=B_{\star }{R}_{\star }^3/2R_{\text{LC}}^2$ and γ_pc = eΦ_pc/mc². The detail of particle trajectories and structure of the electromagnetic field is important to precisely characterise their acceleration, which is the aim of our simulations. However, these theoretical estimates are useful to better understand and rescale the simulation outputs.

4.1. Particle spectra

The spectra of electrons, positrons, and protons are illustrated in Fig. 6 for f_pp = 0.01, f_pp = 0.05, and f_pp = 0.1. We compare the spectra obtained for all particles in the magnetosphere and for particles escaping the magnetosphere. To compute the spectra of escaping particles, we calculate the total number of particles comprised in the spherical shell between 0.8 r_max and 0.9 r_max, such as u_r >  0. Thus, the normalisations of the total and escaped spectra are only indicative and should not be compared.

Fig. 6. Spectra γdN/dγ normalised by the total number of particles in the magnetosphere Ntot, for electrons (blue), positrons (red), and protons (green), for all particles (dashed lines) and for particles escaping the magnetosphere (solid lines). We compare high, medium, and low production of pairs, fpp = 0.01, fpp = 0.05, and fpp = 0.1, respectively (top to bottom). Dotted lines indicate the maximum Lorentz factors γpc, p, γ0, p and γpc, e.

Each particle species shows different spectra, as they experience different fates in the magnetosphere. Moreover, the spectra obtained for a high and low pair production show large discrepancies. We note that for f_pp = 0.01, close to the force-free regime, the highest energy electrons and positrons are accelerated to Lorentz factors close to the vacuum polar cap Lorentz factor or pairs γ_pc, e. The electrons accelerated to the highest energies, close to the Y-point, do not escape as they fall back onto the neutron star surface (Cerutti et al. 2015). The maximum Lorentz factors of protons are close to the vacuum polar-cap Lorentz factor of protons γ_pc, p (Philippov & Spitkovsky 2018). For f_pp = 0.05, in the transition regime between force-free and electrosphere, electrons, positrons, and protons tend to reach higher energies due to higher unscreened electric fields. Regarding particles escaping the magnetosphere, as electrons are more confined in the polar flows, the ones escaping have lower energies. Moreover, the current sheet is thin and thus only a small fraction of positrons escape at the highest energies. A significant fraction of the highest energy protons that are not confined in the current sheet escape the magnetosphere. The proton spectrum shows a peak associated with these escaping protons, which was absent in the force-free regime. For f_pp = 0.1, close to the disc-dome configuration, electrons are accelerated to lower energies due to their confinement in the polar flows. The electrons accelerated to the highest energies do not escape as well. Protons are accelerated to higher Lorentz factors than in the force-free regime, close to the vacuum maximum Lorentz factor of protons γ_0, p.

4.2. Trajectories and acceleration

In all simulations, a fraction of the injected protons are systematically accelerated and escape the magnetosphere. Most of the protons injected at every time step directly fall back on the neutron star surface. Only protons injected at high latitudes escape, typically at θ ∼ 0.6−0.7 rad (and π − θ) for f_pp = 0.01 which coincides with the footpoints of the separatrix current layers, and at θ ∼ 0.9 rad (and π − θ) for f_pp = 0.1. In comparison, the polar cap angle is θ_pc ∼ 0.46 rad. Protons injected at lower latitudes are trapped in the closed field line region and wrap around the neutron star. Protons injected at the highest latitudes are the ones accelerated to the highest energies and have thus the highest chance of escaping. Protons that escape have quasi-radial trajectories at large distances.

Trajectories of protons, positrons, and electrons projected in a poloidal plane are illustrated in Fig. 7 for f_pp = 0.01. We select two particles for each species, which propagate in the equatorial region where particles are accelerated at the highest energies. Considering individual trajectories of protons and positrons f_pp = 0.01, we note that we retrieve the Lorentz factors by calculating γ ≃ [1 + (∫ dt(eE_∥/mc)²)]^0.5, where E_∥ = E ⋅ v/v, with v being the particle velocity. As expected, particles are accelerated by the electric field component parallel to their trajectory. We see that protons and positrons are efficiently accelerated along the separatrices, below the light cylinder radius. Protons can be slightly accelerated at larger distances when they cross the current sheet and experience unscreened parallel electric fields. Their Lorentz factors reach local minima at the crossing points with the current sheet. The positrons are more efficiently confined and thus accelerated in the current sheet, where they acquire a larger fraction of their final Lorentz factors. However, due to the development of instabilities (such as kink and tearing instabilities) and magnetic islands in the current sheet (see Fig. 1), a wide variety of trajectories can be observed for positrons and electrons. For instance we show the trajectory of one electron that escapes in the equatorial region despite its negative charge (dashed green line). Protons that are less sensitive to the structure of the current sheet display more similar trajectories.

Fig. 7. Example of proton, positron and electron trajectories (from top to bottom) in the equatorial region for fpp = 0.01. Markers of different colours link Figs. 7 and 8.

Fig. 8. Lorentz factors corresponding to the proton, positron, and electron trajectories shown in Fig. 7.

Protons gain a large fraction of their final energy within the separatrix current sheets inside the light-cylinder radius. For a high pair production with f_pp = 0.01, most of the escaping protons gain 75% of their maximum Lorentz factor below R_LC. For a low pair production with f_pp = 0.1, protons gain 75% of their maximum Lorentz factor below 0.6R_LC. The fate of positrons is different, as they can be confined in the current sheet and thus accelerated at larger distances, with a significant contribution from magnetic reconnection, which is consistent with previous studies (e.g. Cerutti et al. 2015). To more accurately interpret this important result, it is useful to look at the parallel component of the electric field in both magnetospheric regimes, E ⋅ B, as shown in Fig. 9. In the force-free-like solution, E ⋅ B ≈ 0 almost everywhere as it should except within the separatrix and equatorial current sheets. Protons fully experience the separatrix electric gap from their injection at the star surface up to the light cylinder where they flow perfectly along the field lines. Beyond the Y-point, they do not experience significant additional acceleration because their trajectories present large oscillations and therefore do not probe the scale of reconnection electric field in the equatorial layer set by the pairs. In contrast, positrons are mostly created outside the light cylinder and are well confined within the equatorial current sheet where they are accelerated by reconnection. In the electrosphere-like configuration, large vacuum gaps fill the magnetosphere outside of the electronic dome and proton disc. As a result, the few protons leaving the dome are quickly accelerated and experience the full vacuum potential.

Fig. 9. Parallel electric field E ⋅ B/B2 for r ≤ 2RLC, for fpp = 0.01 (left) and fpp = 0.05 (right).

The maximum Lorentz factor of escaping protons as a function of pair production efficiency is shown in Fig. 10 (top panel). We see that protons experience a fraction of the full vacuum potential drop (higher than the polar cap potential drop for f_pp >  0.01). This fraction is small for high production of pairs, increases when the production of pairs is reduced, and saturates at a maximum value for low or no pair production ∼0.75γ_0, p. The densities of electrons and positrons in the closed field line region, where protons are mostly accelerated, are high for high production of pairs, and therefore the high plasma multiplicities screen the parallel electric field and prevent protons from experiencing a large fraction of the full vacuum potential drop. For low production of pairs, only protons are present in the equatorial plane and can experience a large fraction of the full vacuum potential drop. We note that the magnetic field dependence of the proton maximum Lorentz factor γ_0, p seems to be well reproduced by the simulations.

Fig. 10. Maximum Lorentz factor of escaping protons, normalised by γ0, p, as a function of fpp for B⋆ = 1.1 × 105 G (top) and as a function of B⋆ for fpp = 0.05 (bottom). The dots correspond to simulation results, whereas the lines represent the maximum Lorentz factor of protons experiencing the vacuum potential drop, from pole to equator (orange, dashed) or across the polar cap (green, dotted).

In the case of escaping positrons, their maximum Lorentz factor also increases with decreasing pair production (an increasing f_pp), between f_pp = 0.01 and f_pp = 0.04. For lower pair productions f_pp ≥ 0.05, the number of positrons produced strongly decreases, the current sheet does not form, and thus the maximum Lorentz factor of escaping positrons drops.

4.3. Proton maximum energy in real pulsars

The estimates of the proton Lorentz factors cannot be directly related with realistic cases as the magnetic field, neutron star radius, and mass ratio are downscaled in our numerical experiments. A rescaling procedure is therefore required. In the formula that we use for extrapolation, several quantities intervene such as the radius of the star, the rotation frequency, and the magnetic field. In the range accessible with our simulations, we have performed several series of tests, varying the magnetic field strength B_⋆ and radius of the neutron star R_⋆ as well as the mass ratio m_r in several sets of simulations in order to check the impact of these parameters on the maximum energy and luminosity. These tests validate the dependencies that intervene in our extrapolation. For instance, as illustrated in the bottom panel of Fig. 10, the maximum Lorentz factor of protons max(γ) appears to be a nearly constant fraction of γ_0, p, which suggests that max(γ) is proportional to B_⋆. We note that the maximum Lorentz factor obtained for the lowest magnetic field is a higher fraction of γ_0, p, which is certainly due to the low maximum Lorentz factor and the confusion with thermal protons, as γ_0, p ≃ 3.5 for B_⋆ ∼ 10⁴ G.

Despite these tests, we caution that this rescaling procedure is a delicate process owing to the large difference between numerical and realistic scales. The quantities that we derive should therefore be considered with care. We assume that a constant fraction of the full vacuum potential drop can be channelled into proton acceleration. In our simulations, we obtain maximum Lorentz factors of between 15 and 75% of γ_0, p, from a high to a low pair production, respectively. As ${\gamma }_{0,\mathrm{p}}=3.3\times {10}^7{m}_{\mathrm{r},1836}^{-1}{B}_{\star ,9}{R}_{\star ,6}^2{P}_{-3}^{-1}$, we see that protons can be accelerated up to $E_{\mathrm{p}}\simeq 5\times {10}^{15}\mathrm{eV}{B}_{\star ,9}{R}_{\star ,6}^2{P}_{-3}^{-1}$ for a high pair production and up to $E_{\mathrm{p}}\simeq 2\times {10}^{16}\mathrm{eV}{B}_{\star ,9}{R}_{\star ,6}^2{P}_{-3}^{-1}$ for a low pair production. These estimates have been derived for typical properties of millisecond pulsars, with B_⋆ = 10⁹ G and P = 10⁻³ s, and a correct electron to proton mass ratio. Thus millisecond pulsars could produce cosmic rays at PeV energies. This could have interesting observational consequences, such as the production of gamma rays in the Galactic centre region (Guépin et al. 2018). For newborn pulsars with millisecond periods, we obtain $E_{\mathrm{p}}\simeq 5\times {10}^{19}\mathrm{eV}{B}_{\star ,13}{R}_{\star ,6}^2{P}_{-3}^{-1}$ for a high pair production and up to $E_{\mathrm{p}}\simeq 2\times {10}^{20}\mathrm{eV}{B}_{\star ,13}{R}_{\star ,6}^2{P}_{-3}^{-1}$ for a low pair production. Therefore, we show that newborn pulsars with millisecond periods could produce cosmic rays up to ultra-high energies, as proposed in several studies (Blasi et al. 2000; Fang et al. 2012, 2013; Lemoine et al. 2015; Kotera et al. 2015). We caution that the effect of curvature radiation is underestimated in our simulations because of the downscaled magnetic field and radius of the neutron star and the subsequent low Lorentz factors of accelerated particles. As curvature radiation can strongly limit particle acceleration below the light cylinder radius (Arons 2003), a realistic treatment could therefore impact the acceleration regions and maximum energies of particles, and should therefore be studied in future work. The cases of normal pulsars and millisecond magnetars are difficult to explore with our simulations due to the large distance between the star and the light cylinder radius, or the high magnetic fields. In particular, extreme magnetic fields could have consequences on pair production processes. These configurations therefore require dedicated studies.

4.4. Energy dissipation and luminosity

One last important quantity to infer is the total energy dissipated and channelled into particles, which allows us to estimate the proton luminosity. As illustrated in Fig. 11, the production of pairs has a strong impact on the outgoing Poynting flux; it decreases strongly with a decrease of the yield of pair production. Furthermore, we note that it can be larger than the analytical spin-down power of an aligned pulsar $L_0=c{B}_{\star }^2{R}_{\star }^6/4R_{\text{LC}}^4$ (e.g. Contopoulos et al. 1999; Spitkovsky 2006) for high pair production, as the Y-point is located below r = R_LC and thus a larger fraction of field lines are open. Moreover, it is less than 20% of L₀ for low pair productions. Therefore, aligned pulsars with low pair production barely spin-down, as expected for the disc-dome solution (Cerutti et al. 2015).

Fig. 11. Energy dissipation as a function of radius for fpp = 0.01, fpp = 0.05, and fpp = 0.1 (top to bottom). We show the radial Poynting flux integrated over a sphere of radius r (blue line), the luminosity in electrons (orange dashed line), in positrons (green dot-dashed line), and in protons (red dotted line), and the sum of all these components (black line), normalised by L0. We caution that the vertical scales of the figures are different.

Energy dissipation is illustrated in Fig. 11, where we show the radial outgoing Poynting flux and luminosity in electrons, positrons, and protons for f_pp = 0.01, f_pp = 0.05 and f_pp = 0.1 as a function of r/R_LC. We caution that the scales of the figures are different. These quantities are smoothed over several time-steps and radial bins in order to display the results clearly. In our simulations, the energy dissipated is self-consistently transferred to particles that are accelerated. A main and irreducible source of magnetic dissipation is via magnetic reconnection which operates in the equatorial current sheet. The separatrices are also a source of dissipation, as they form cavities that allow particle acceleration, with the electric field accelerating particles along the magnetic field lines. The total power shows significant variations with radius, which demonstrates the occurrence of nonstationary phenomena. These irregularities reflect the strong time dependency of reconnection and particle acceleration via the formation of plasmoids (see strong peaks in Fig. 11 for f_pp = 0.01) and kinks in the current sheet, and the related shifts of the Y-point position.

For f_pp = 0.01, the dissipation of radial Poynting flux into particle kinetic energy occurs mostly around and beyond the Y-point, which is located at approximately r = 0.8R_LC. The energy is mostly dissipated into positron kinetic energy. Energy is also dissipated below the Y-point along the gaps where the parallel electric field is not completely screened (see Fig. 9). The fraction of the Poynting flux dissipated into electron and proton kinetic energy decreases with increasing f_pp. For f_pp = 0.1, a significant fraction of the Poynting flux is dissipated into proton kinetic energy.

These simulations allow us to evaluate the typical proton luminosity. As illustrated in Fig. 12, the maximum proton luminosity is obtained for a high production of pairs. For a decreasing pair production, the proton luminosity first decreases, and then increases again. Given the limited number of simulations that we can perform, it is difficult to determine with certainly the minimum proton luminosity. We obtain the minimum proton luminosity L_p ≃ 2 × 10⁻³L₀ for f_pp = 0.03 and the maximum proton luminosity L_p ≃ 4 × 10⁻²L₀ for f_pp = 0.01. Assuming that we can use these fractions for typical pulsar properties, and considering the value of the spin-down power of an aligned pulsar $L_0=1.4\times {10}^{37}{\mathrm{erg}\mathrm{s}}^{-1}{B}_{\star ,9}^2{R}_{\star ,6}^6{P}_{-3}^{-4}$ for millisecond pulsar properties, we obtain $L_{\mathrm{p}}\simeq 3\times {10}^{34}-5\times {10}^{35}{\mathrm{erg}\mathrm{s}}^{-1}{B}_{\star ,9}^2{R}_{\star ,6}^6{P}_{-3}^{-4}$. For newborn pulsars with millisecond periods, we obtain $L_{\mathrm{p}}\simeq 3\times {10}^{42}-5\times {10}^{43}{\mathrm{erg}\mathrm{s}}^{-1}{B}_{\star ,13}^2{R}_{\star ,6}^6{P}_{-3}^{-4}$.

Fig. 12. Maximum proton luminosity normalised by L0, as a function of fpp for B⋆ = 1.1 × 105 G.

5. Discussion and conclusions

Here, we use 2D PIC simulations to study the impact of pair production on the acceleration of protons in aligned pulsar magnetospheres. These simulations confirm that pulsar magnetospheres are good candidates for the acceleration of protons: regardless of the yield of pair production, protons can be accelerated and escape. Interestingly, due to the mass ratio and large density contrast between protons and pairs, protons do not experience the same trajectories, and thus acceleration, as pairs; they are mostly accelerated below the light cylinder radius within the separatrix current layers but they are not confined in the equatorial current sheet when it exists, whereas pairs are accelerated at their highest energies at the Y-point and beyond in the equatorial current sheet. We note that higher magnetic fields could enhance pair production below the light cylinder radius, as mentioned in Philippov & Spitkovsky (2018), and screen the parallel electric field that accelerates protons in this region. Thus, protons could be mostly accelerated at larger distances in the current sheet. This effect could be investigated in future studies.

In our numerical experiments, magnetic energy is dissipated into particle kinetic energy through magnetic reconnection in the current sheet and acceleration by the unscreened electric field along the separatrices. Similar dissipation has been observed in different studies with different codes, even using different methods such as resistive force-free MHD (see, e.g., Parfrey et al. 2012). Gamma-ray observations give us a lower limit of the dissipation rate which is given by the GeV power output compared with the pulsar spindown. This rate is comprised between 1 and 10% and probably even higher in some cases. To this, a bolometric correction factor would be required as well as the efficiency factor of the accelerated particles. The dissipation reported here and in other studies is at least compatible with observed lower limits. For a high production of pairs, about 2% of the theoretical pulsar spin-down power of an aligned pulsar L₀ is channelled into protons, which is of the same order as the fraction required to fit the ultra-high-energy cosmic-ray (UHECR) spectrum (Fang et al. 2013). In this case, the maximum particle Lorentz is limited by the potential drop across the pulsar polar cap (Cerutti et al. 2015; Philippov & Spitkovsky 2018). For a low production of pairs, less than 0.05% of L₀ is channelled into protons. However, despite these low luminosities, it appears also that protons are accelerated to higher energies for a low production of pairs. The maximum Lorentz factor approaches and is limited by the theoretical maximum given by the total vacuum potential drop.

Comparisons between our simulation results and analytical estimates of the maximum Lorentz factor and spin-down luminosity allow us to estimate the maximum energy and luminosity of accelerated protons for realistic parameters. We caution that these estimates should be considered with care owing to the difference between numerical and realistic scales. For typical properties of millisecond pulsars, protons could be accelerated at E_p ≃ 1 PeV with luminosities of L_p ≃ 5 × 10³⁵ erg s⁻¹, which might have interesting observational consequences for the production of gamma rays (e.g. Guépin et al. 2018). Moreover, for typical properties of newborn pulsars with millisecond periods, protons could be accelerated to ultra-high energies, up to E_p ≃ 10 EeV with luminosities L_p ≃ 5 × 10⁴³ erg s⁻¹.

Typically, a few percent of the spin-down power of the total population of newborn millisecond pulsars is required to reproduce the observed UHECR spectrum (e.g. Fang et al. 2013), which seems achievable given our results and therefore supports the idea that newborn pulsars with millisecond periods are good candidate sources for the production of UHECRs. We note that heavy nuclei are required to explain the UHECR spectrum at the highest energies. The extraction of heavy ions from the neutron star surface, as already discussed in Chen et al. (1974) and Kotera et al. (2015) for instance, as well as their propagation and interaction in the magnetosphere could be explored in a subsequent study. As the neutron star crust is mostly composed of Fe or Ni elements, heavy elements could be naturally extracted from these objects. Extracting ions from the crust could be allowed by the high induced electric fields. Particles bombarding the neutron star surface could also help the extraction. Moreover, as mentioned in Sect. 4.3, we caution that the energy losses of particles that are due to curvature radiation are not fully taken into account in our simulations. As calculated in Arons (2003), as a consequence of curvature radiation, the energy of accelerated protons below the light-cylinder radius is limited to E_p ≃ 10^16.5 eV for pulsars with millisecond periods. This effect should be tested in future work in order to assess its impact on numerical simulations and whether or not protons and heavy nuclei can be accelerated at larger distances, in the current sheet or in the wind.

We note that this work is restrictive as only a small fraction of the pulsar wind is included in our simulations, and we do not account for energy losses or re-acceleration of protons at larger distances, for instance at a shock front. Moreover, we consider the case of an aligned pulsar, and the structure of the magnetosphere should be modified in the misaligned case (Spitkovsky 2006; Pétri 2016). However, the structure of the magnetosphere below the light cylinder radius, where most of the acceleration of protons seems to take place, should be similar for the misaligned configuration. We highlight that reconnection taking place in the striped wind or Fermi-type acceleration taking place at the termination shock between the pulsar wind and its nebula (e.g. Lemoine et al. 2015) may further increase the proton maximum energy in the outer regions.

In addition to studying the impact of the pair production strength on proton acceleration in pulsar magnetospheres, we performed several additional tests in order to assess the impact of the size of the simulation domain, the resolution, and the particle injection rate. First, simulations performed by completely shutting down the pair production process lead to the same steady state as the simulations with a low pair production (i.e. f_pp >  0.15). The structure of magnetospheres with a low production of pairs appears to be more affected by simulation parameters: the size of the simulation domain influences the extent in latitude of the polar flows of electrons, and the resolution influences the density of the equatorial flow of protons. Moreover, the injection rate influences these two characteristics. In particular, we noticed that simulations performed with a lower resolution show a larger number of protons escaping in the equatorial flow. This difference could result from numerical effects, or the longer times required to reach the steady state for higher resolutions, and is currently under study. However, despite small morphological differences, the general structure of the magnetosphere is similar and our main conclusions remain unchanged. In particular, the maximum energy of particles and the acceleration regions are not affected by these minor structural changes.

Further work will be required to better characterise the escape of protons by a detailed modelling of their trajectories. Moreover, the link between the simulated amount of pair production and pair multiplicities in realistic environments should be explored. For this purpose, a self-consistent modelling of pair production, but also of other types of interactions, will be required. The present study focuses on the aligned magnetosphere. The anti-aligned and inclined cases could be studied in future work. The anti-aligned configuration would reverse the charge densities at the surface of the neutron star, allowing for proton extraction at the poles and electron extraction in the equatorial region. This configuration could lead to interesting magnetospheric configurations because of the mass separation between protons and electrons. Given the fate of the electrons extracted from the poles in the aligned case, the acceleration of protons in the anti-aligned case could be suppressed. This interesting question should be investigated in a dedicated study.

Acknowledgments

The authors thank the anonymous referee for useful feedback on this work. We thank A. Spitkovsky, A. Philippov, F. Mottez and J. Pétri as well as the organizers and participants to the “Entretiens pulsars” for useful discussions. This work is supported by the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche as well as CNES. CG is supported by a fellowship from the CFM Foundation for Research and by the Labex ILP (reference ANR-10-LABX-63, ANR-11-IDEX-0004-02). This work was granted access to the HPC resources of CINES on Occigen under the allocations A0030407669 and A0050407669 made by GENCI. This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us.

References

All Tables

All Figures

Fig. 1. Density maps for electrons, positrons, and protons (left to right), for fpp = 0.01, fpp = 0.05 and fpp = 0.1 (top to bottom), as a function of r/RLC and θ, for t = 5P. The densities are normalised by ${(r/R_{\star })}^{-2}{n}_{\text{GJ}}^{\star }$, where $n_{\text{GJ}}^{\star }=B_{\star }/2\pi R_{\text{LC}}e$ is the polar GJ number density. Solid black lines are the magnetic field lines. The dashed black line indicates the distance from the rotation axis r sin θ = RLC and the grey semi-disc represents the neutron star. The densities are in log scale in order to enhance the contrast.

In the text

	Fig. 2. Same as Fig. 1 but for the zoomed-in simulation domain r <  RLC.
In the text

	Fig. 3. Plasma densities (n+ + n−)/nGJ as a function of fpp at the pole and at the Y-point (respectively blue and orange lines).
In the text

	Fig. 4. Radial current density for r = 2RLC, averaged over θ and normalised by the polar GJ current density ${(R_{\star }/r)}^2{J}_{\text{GJ}}^{\star }$, as a function of fpp. We show the total radial current density (circles) and the contributions of electrons (crosses), positrons (plus signs), and protons (triangles).
In the text

	Fig. 5. Charge density maps (top) and radial current maps (bottom) for a high production of pairs fpp = 0.01 (left) and a low production of pairs fpp = 0.05 (right), for r ≤ RLC. ${\rho }_{\text{GJ}}^{\star }$ and $J_{\text{GJ}}^{\star }$ are the GJ charge and current densities at the poles.
In the text

Fig. 6. Spectra γdN/dγ normalised by the total number of particles in the magnetosphere Ntot, for electrons (blue), positrons (red), and protons (green), for all particles (dashed lines) and for particles escaping the magnetosphere (solid lines). We compare high, medium, and low production of pairs, fpp = 0.01, fpp = 0.05, and fpp = 0.1, respectively (top to bottom). Dotted lines indicate the maximum Lorentz factors γpc, p, γ0, p and γpc, e.

In the text

	Fig. 7. Example of proton, positron and electron trajectories (from top to bottom) in the equatorial region for fpp = 0.01. Markers of different colours link Figs. 7 and 8.
In the text

	Fig. 8. Lorentz factors corresponding to the proton, positron, and electron trajectories shown in Fig. 7.
In the text

	Fig. 9. Parallel electric field E ⋅ B/B2 for r ≤ 2RLC, for fpp = 0.01 (left) and fpp = 0.05 (right).
In the text

	Fig. 10. Maximum Lorentz factor of escaping protons, normalised by γ0, p, as a function of fpp for B⋆ = 1.1 × 105 G (top) and as a function of B⋆ for fpp = 0.05 (bottom). The dots correspond to simulation results, whereas the lines represent the maximum Lorentz factor of protons experiencing the vacuum potential drop, from pole to equator (orange, dashed) or across the polar cap (green, dotted).
In the text

Fig. 11. Energy dissipation as a function of radius for fpp = 0.01, fpp = 0.05, and fpp = 0.1 (top to bottom). We show the radial Poynting flux integrated over a sphere of radius r (blue line), the luminosity in electrons (orange dashed line), in positrons (green dot-dashed line), and in protons (red dotted line), and the sum of all these components (black line), normalised by L0. We caution that the vertical scales of the figures are different.

In the text

	Fig. 12. Maximum proton luminosity normalised by L0, as a function of fpp for B⋆ = 1.1 × 105 G.
In the text