OSIRIS-GR: General relativistic activation of the polar cap of a compact neutron star

We present ab initio global general-relativistic Particle-in-cell (GR-PIC) simulations of compact millisecond neutron star magnetospheres in the axisymmetric aligned rotator configuration. We investigate the role of GR and plasma supply on the polar cap particle acceleration efficiency - the precursor of coherent radio emission - employing a new module for the PIC code OSIRIS, designed to model plasma dynamics around compact objects with fully self-consistent GR effects. We provide a detailed description of the main sub-algorithms of the novel PIC algorithm, including a charge-conserving current deposit scheme for curvilinear coordinates. We demonstrate efficient particle acceleration in the polar caps of compact neutron stars with denser magnetospheres, numerically validating the spacelike current extension provided by force-free models. We show that GR relaxes the minimum required poloidal magnetospheric current for the transition of the polar cap to the accelerator regime, thus justifying the observation of weak pulsars beyond the expected death line. We denote that spin-down luminosity intermittency and radio pulse nullings for older pulsars might arise from the interplay between the polar and outer gaps. Also, narrower radio beams are expected for weaker low-obliquity pulsars.


Introduction
Neutron stars are exquisite physics laboratories coupling extreme magnetic fields with strong gravitational fields.These compact astrophysical objects were discovered in 1967 as radio pulsars (Hewish et al., 1968).Since then pulsars have been observed as bright sources of electromagnetic radiation, from radio to gamma rays.Despite recent advances in understanding the intricate dynamics of their magnetospheres, propelled in part by the available computational power, we still lack a complete picture that explains the observed radiation features.Goldreich and Julian (1969) noted that an isolated neutron star modelled by a rotating spheric perfect conductor that is threaded by a dipolar magnetic field would induce an electric field strong enough to extract charged particles from its surface, populating its magnetosphere with a dense chargeseparated plasma.In an ultra-strong magnetic field, the electromagnetic force dominates, leading to the force-free condition, where plasma shorts out the electric field parallel to the magnetic field.Consequently, particles move along magnetic field lines, which corotate rigidly with the neutron star.A crucial element in the magnetosphere is the light cylinder radius, R LC = c/Ω * , the distance from the rotational axis after which corotation is not possible due to its tangential speed exceeding the light speed.This boundary separates the closed corotating field lines from the open bundle originating from the polar cap of the neutron star.The outflowing particles generate poloidal currents that sweep back the open field lines, inducing a toroidal component.Beyond the light cylinder radius, an equatorial current sheet forms to support the discontinuity of the toroidal magnetic field component.In this region, the poloidal electric field and the toroidal magnetic field result in an outflowing Poynting flux that we observe as a spin-down luminosity.Contopoulos et al. (1999) confirmed this model by numerically solving the "pulsar equation" (Scharlemann and Wagoner, 1973;Michel, 1973), the time-independent version of force-free electrodynamics.This work established the standard magnetospheric picture, providing the complete magnetospheric current closure with a y-shaped current layer flowing between the open and closed field lines.
Observations of pulsar wind nebulae (PWNe), i.e. rotationpowered structures formed from the ejection of plasma advected from the inner pulsar magnetosphere (Kirk et al., 2009;Mitchell and Gelfand, 2022), support the force-free solution given that there is a dense plasma to screen the nonideal electric field everywhere (Gruzinov, 2005;Spitkovsky, 2006;Kalapotharakos, C. and Contopoulos, I., 2009;Timokhin, 2007;Bai and Spitkovsky, 2010;Pétri, 2012;Philippov and Spitkovsky, 2014).However, these models do not constrain the mechanism that provides the plasma to drive the required magnetospheric currents.Without external plasma sources, there must exist dissipative regions -gaps -within the magnetosphere where particles accelerate and produce pairs that cancel the accelerating electric field.These gaps would operate at much shorter time scales (microseconds) than the magnetospheric time scales (milliseconds to seconds, thus quasi-stationary), consistent with the force-free solution.Several studies explored the observational implications of the polar cap gap (Daugherty and Harding, 1982;Usov and Melrose, 1995;Daugherty and Harding, 1996;Ruderman and Sutherland, 1975), slot gap (Arons and Scharlemann, 1979;Arons, 1983;Muslimov andHarding, 2003, 2004), and outer gap models (Cheng et al., 1986;Romani, 1996;Takata et al., 2004;Hirotani, 2008), highlighting the importance of these gaps not only to provide the necessary plasma to the magnetosphere but also to explain the observed electromagnetic spectrum (e.g., Jun-Tao et al. (2021)).
The radio beam in radio pulsars originates from a region above the pair formation front of the polar cap gap due to the nonlinear screening of the accelerating electric field.The misalignment of the magnetic and rotation axis leads to a pulsed behaviour as the beam crosses our line of sight periodically, similar to a lighthouse.These radio pulses contain microstructures, some possessing single (Drake and Craft, 1968;Deshpande and Rankin, 1999;Vivekanand and Joshi, 1999) or multiple drifting components (Qiao et al., 2004;Basu et al., 2018).Models to explain this behaviour assume a limitation in the extracted charged particles from the stellar surface such that the spark filaments no longer corotate with it.Hence, the sub-pulses undergo an E × B drift different from the field lines due to the parallel electric field at the base of those field lines.Ruderman and Sutherland (1975) (hereafter RS75) assumed the case where no particles could be lifted from the stellar surface, thus forming a vacuum gap above it.The intermediate model of Gil, J. et al. (2003) allows a thermionic outflow from the surface forming a partially screened gap (PSG).This model accommodates a necessary reduction of the drifting speed and provides a better estimate of the X-ray luminosity from the backflow bombardment of particles.The space charge limited flow model (SCLF) by Arons and Scharlemann (1979), where particles can freely escape from the stellar surface, does not provide a natural explanation for this phenomenon.The binding energies of either ions or electrons should not be too high, thus motivating the SCLF model (Medin and Lai, 2010).Alternatively, studies have highlighted the possibility of strongly bound particles in bare strange quark stars (Ren-xin and Guo-jun, 1998;Xu et al., 1999;Yu and Xu, 2011), favouring the RS75 or PSG models.
First principle global magnetospheric models that capture the kinetic plasma scales are the best way to introduce the effects of these dissipative regions.Several works have proposed tools to tackle this very computationally demanding task in one-dimension (Timokhin, 2010;Timokhin and Arons, 2012), two-dimensions (Chen and Beloborodov, 2014;Belyaev, 2015b,a;Cerutti et al., 2013Cerutti et al., , 2015;;Philippov et al., 2015a;Cerutti, B. and Philippov, A. A., 2017;Cruz et al., 2021;Hu and Beloborodov, 2022;Bransgrove et al., 2022;Cruz et al., 2023) or three-dimensions (Philippov et al., 2015b;Cerutti and Beloborodov, 2016;Philippov and Spitkovsky, 2018;Kalapotharakos et al., 2018;Brambilla et al., 2018;Cerutti, Benoît et al., 2020;Hakobyan et al., 2023).The work by Philippov et al. (2015b) concluded that the polar cap of low obliquity rotators has weak particle acceleration, suppressing the pair production and shutting down the radio pulsar mechanism.Several authors later introduced the general-relativistic frame-dragging effect (Beskin, 1990;Muslimov and Tsygan, 1992) into Faraday's induction equation, demonstrating that this effect could potentiate the ignition of the polar cap even for the aligned rotator (Philippov et al., 2015a;Chen et al., 2020;Philippov et al., 2020;Bransgrove et al., 2022) as predicted by general-relativistic force-free models (Belyaev and Parfrey, 2016;Gralla et al., 2016;Huang et al., 2018).Although such kinetic models are able to capture the polar cap ignition, qualitative corrections from neglecting the lapse function and generalrelativistic magnetic field topology (Torres et al., 2023) are not resolved.One of such parameters of utmost importance is the spacelike angle θ SL , expected to limit the poloidal extension of the observed pulsar radio emission generation volume.This angle was obtained analytically for force-free models (Belyaev and Parfrey, 2016;Gralla et al., 2016), but was not yet verified using general-relativistic particle-in-cell simulations.Thus, ab initio simulations of neutron star magnetospheres with fully self-consistent general-relativistic effects still lack in the community.
The purpose of this paper is to present the first simulations of compact millisecond neutron stars including self-consistent general-relativistic effects, using a charge-conserving GR-PIC code.Also, we summarize all the necessary techniques required for the implementation of a GR-PIC code.As far as we know, a systematic description of the GR-PIC algorithm does not exist in the literature.We detail the methodology and algorithms used to extend the particle-in-cell code OSIRIS to include GR effects in Sect. 2. In Sects. 3 and 4, we perform numerical experiments of such magnetospheres by varying the plasma injection methods to explore the full range of magnetospheric solutions, from quasi-force-free to full-charge-separated solutions, expanding on the consequent observational implications.We show that the correct dimensions of the polar cap and spacelike current region can only be obtained when considering the lapse function effects.Also, we demonstrate that charge-separated magnetospheric solutions possess smaller values of the spacelike angle with respect to force-free (FF) estimates.Therefore, we conjecture the existence of narrower radio beams from weak pulsars, some still observable beyond the expected death line.Conclusions and future prospects are presented in Sect. 5.

GR-PIC method
In this section, we detail the extension module of general relativity to the OSIRIS particle-in-cell framework (Fonseca et al., 2002).We describe the main equations and integration algorithms implemented in OSIRIS including the field solver, equations of motion, charge conservation, and boundary conditions.A general description of the techniques and methodologies (independent of the internal specifics/implementation of the PIC code) is provided to facilitate their implementation in other PIC codes.
Throughout the paper, we use units in which c = G = 1, the (−, +, +, +) metric signature, and stellar properties are characterized by an asterisk in subscript (e.g., R * for the stellar radius).In addition, we adopt the 3+1 formalism of general relativity (Thorne and MacDonald, 1982).

Metric
The astrophysical objects we are interested in constitute the class of compact objects where the background metric is no longer a flat (Minkowski) metric.We focus our study on the dynamics of the magnetospheres of black holes and neutron stars.The electromagnetic energy density in such environments is negligibly small compared to their corresponding total mass density.Consequently, we can assume a fixed background metric for these systems.Black holes and neutron stars are modeled by the Kerr and the Hartle-Thorne metric (Hartle and Thorne, 1968), respectively.However, as neutron stars rotate at significantly smaller angular velocities, one can invoke the slowrotation approximation and simplify the corresponding metric (R * Ω * ≪ 1) for neutron stars.In this approximation, and keeping only the linear terms of the angular velocity, the exterior metric takes the form (Hartle, 1967): where α = √ 1 − R s /r and β i = (0, 0, −ω) are the lapse function and shift vector, respectively, γ i j is the spatial part of the metric, R s is the Schwarzschild radius, and dΩ 2 = dθ 2 + sin 2 θdϕ 2 .Interestingly, this metric is equivalent to the slow-rotation approximation of the Kerr metric, which leads us to adopt the Boyer-Lindquist coordinate system (t,r,θ,ϕ).Currently, OSIRIS has the Minkowski, the Schwarzschild, and the slow-rotation approximation of the Kerr metric implemented.The Minkowski metric is recovered in the zero compactness limit (R s → 0), and the Schwarzschild metric by neglecting the frame-dragging effect (ω(r) → 0).The latter effect is a correction that captures the differential rotation ω(r) associated with a free-falling inertial frame (Ravenhall and Pethick, 1994): We define fiducial observers (FIDOs) with world lines normal to spatial hypersurfaces of constant "universal" time.The lapse function relates the rate at which the local fiducial observer clock ticks in relation to the universal time.Hence, the lapse function translates local time-measured quantities to the universal time.The shift vector measures the angular velocity at which the coordinate frame shears in relation to the FIDOs.We select FIDOs corotating with the absolute space, thus considering zero angular momentum observers (ZAMOs).In these hypersurfaces, we use the orthonormal basis vectors and vectorial components, i.e. e î ≡ e i / √ γ ii and A î ≡ √ γ ii A i such that A = A i e i = A îe î.This basis allows for a trivial extension of three-dimensional vectorial operations to curved geometries.Also, the definition of line, area, and volume elements is modified according to eq. ( 1) -see Appendix A.

Grid and coordinates
So far, we have addressed the metric that defines the system dynamics.In this section, we exploit the usage of different coordinate systems for the same metric (e.g., cartesian or spherical for the Minkowski metric).
Compact objects are naturally described in spherical coordinates.Using body-fitted coordinate grids simplifies the inclusion of system symmetries and boundary conditions.Another degree of freedom is the usage of linear or non-uniform grid spacings.A common feature in similar codes is to allow uniform grid spacings in log r (logarithmic coordinate) and − cos θ (equal area coordinate).The first allows for higher resolution close to the central object, while the second enforces constant charge density of each macro-particle in the meridional direction (Belyaev, 2015a).These two properties have shown a reduction of numerical noise on the axis and increased code stability.Therefore, the OSIRIS-GR module has these two options available.Linear grid spacing options are also available, although these are not explored in the current paper.Another possibility available is to adopt logical coordinates, i.e. (r, θ)=(log r,− cos θ), which possesses the two properties described above and retains uniform grid spacings.The background metric in the new coordinate system takes the form: It is important to highlight that vectors and tensors are identical if written on the orthonormal basis, independently of the chosen coordinates: We should also mention that this is only possible for orthogonal metric systems, i.e. diagonal spatial metrics.Also, the OSIRIS-GR module introduces a new possibility towards modelling compact object magnetospheres via the halfdomain setup.This simulation mode differs from the standard pole-to-pole meridional extension (0 ≤ θ ≤ π) by reducing the poloidal plane in half, thus possessing a pole-to-equator meridional extension (0 ≤ θ ≤ π/2).This novel numerical setup captures the global magnetospheric current closure while reducing the computational cost of emulating the complete poloidal plane.In this way, a higher number of macro-particles can be employed, increasing the statistical significance and reducing the associated numerical noise.This simulation mode is particularly interesting to study the pulsar polar cap dynamics, where the asymmetric mode development of the equatorial current sheet may be neglected to leading order.

Maxwell's equations
In the 3+1 formalism, Maxwell's equations take the form (Komissarov, 2011): Figure 1: Spherical Yee cell adopted in OSIRIS-GR with the location in the numerical cell where each electromagnetic component is defined.
where E, B, j and ρ are physical quantities measured by fiducial observers in co-rotation with the absolute space (ZAMOs).In spherical coordinates, Maxwell's equations are better described in the integral form using the Kelvin-Stokes theorem on a Yeelattice (Belyaev, 2015b;Cerutti et al., 2015).The OSIRIS-GR module adopts the spherical cell shown in Fig. 1.In this way, the curl terms can be written as with Ẽ ≡ αE + β × B and B ≡ αB − β × E, respectively.Expressions for the line and area elements are detailed in Appendix A.
Up to this point, we have assumed β φ as the only non-zero component, with a purely radial dependence in α, and neglected all azimuthal derivatives (axisymmetry).The evolution equations are then expressed as in the typical Yee (1966) scheme with J ≡ α j − ρβ and n being the temporal index.Equation ( 16) also shows that to correctly push in time the electric field we need to gather the particle's half-step charge density and current (see Sect. 2.5 and Appendix B).Due to the complex coupling introduced by equations ( 11) and ( 14), the poloidal components are pushed in time before the azimuthal ones in equations ( 16) and (17).By applying the Newtonian limit, we retrieve the Minkowski version of the evolution equations in Belyaev (2015b); Cerutti et al. (2015).Notice that ( 9)-( 14) can be applied to any axisymmetric system given a diagonal spatial metric (e.g.spherical, cylindrical, toroidal).

Particle pusher
In the 3+1 formalism, the equations of motion for a charged particle (Thorne and MacDonald, 1982;Philippov et al., 2015a;Parfrey et al., 2019) are given by where g is the gravitational acceleration, H is the gravitomagnetic tensor, p = mΓv = mu is the FIDO-measured momentum, q and m are the charge and mass of each macro-particle, Γ = √ ε + u • u is the particle's special relativistic Lorentz factor, and ε = 0 (= 1) for massless (massive) particles.Equation (18) shows that the particle velocity is composed of the FIDO-measured velocity, αu/Γ, and the FIDO's velocity, β.Equation ( 19) details which forces dictate the particle dynamics.From left to right, we identify the Lorentz, the gravitational, the gravitomagnetic, the classical radiation reaction (Vranic et al., 2016;Landau and Lifschitz, 1975), and the "coordinate" forces.The last term comprises the generalized fictitious forces associated with the use of the orthonormal vector basis: where Γ i k j are the spatial Christoffel symbols defined in Appendix A. We advance in time the contravariant positions and the orthonormal contravariant momenta, i.e. r ≡ r i e i and u ≡ u îe î.To do so, we employ a novel variant of the Strang splitting method to advance in time the equations of motion (Strang, 1968), which uses leapfrogged positions and velocities.Thus, defining the positions centered in integer time steps and the u velocities in half-time steps, the discretized equations of motion read: with where the terms with a bar highlight the approximated terms.
The choice for v defines if the simulation will be using the Boris (Boris and Shanny, 1972), the Vay (Vay, 2008), or the Higuera-Cary (HC) method (Higuera and Cary, 2017) for the electromagnetic force.For the systems we are interested in, the HC method is preferred as it is volume-preserving in the phase space and accurately captures the E × B drift even when underresolving the cyclotron frequency (Higuera and Cary, 2017;Ripperda et al., 2018).The algorithm used to evaluate equations ( 22) and ( 24) is the second-order Heun's (or Runge-Kutta) method.This generalization of the leapfrog method allows using very efficient methods, such as the ones employed in flat spacetime particle-in-cell codes, while keeping the time centering of the evolution equations ( 21)-( 24).In each time step, we only need one interpolation of the electromagnetic fields.We check if the relativistic cyclotronic frequency is resolved locally, for each macro-particle, with at least ten temporal steps.If this is not the case, we subtract the particle's perpendicular momenta and add the E × B drift similarly to Bacchini et al. (2020).We perform this verification step when evaluating the Lorentz force.This algorithm to push particles in time is general and can also be applied to model photon dynamics in a fixed background metric (geodesic motion).
In neutron stars, the general relativistic effects in the particle push may be seen as minor corrections to the particle motion very close to the star.In future works, we intend to detail the differences between simulations using a complete GR description and simulations that use a flat spacetime push with the frame-dragging effect in the field solver (Philippov and Spitkovsky, 2018;Chen et al., 2020;Philippov et al., 2020;Bransgrove et al., 2022).This study may reveal if a numerically more expensive GR algorithm is necessary for neutron star magnetosphere simulations.Nevertheless, a complete implementation provides versatility to study plasma dynamics around more compact sources such as black holes.

Current deposition scheme
The only missing piece to complete the particle-in-cell algorithm is the evaluation of the charged current density in equation ( 16).This current term is the source term in the electric field evolution equation and captures the coupling between the electromagnetic fields and the charged particles.In the absence of sources, the Yee algorithm (Yee, 1966) ensures that the time-independent Maxwell's equations are always satisfied if initially satisfied.However, in the presence of sources, the only way to satisfy Gauss's law ( 5) is through the charge conservation condition: An important challenge is then to design an algorithm capable of satisfying this condition to machine precision in an arbitrary curvilinear grid.Such algorithms are readily available for Cartesian grids for any interpolation order (Esirkepov, 2001).
The usual way to go around this challenge for curvilinear grids is to employ a divergence-cleaning algorithm; however, these algorithms are significantly less computationally efficient and, from a numerical point of view, their properties are less well known and explored (in a systematic way) .Here, we present the general-relativistic generalization of the current deposition scheme first described in Cruz et al. (2023) that conserves charge to machine precision in the non-uniform grid in Fig. 1, even for very high stellar compactness values.It is important to note that this general-relativistic charge conservation scheme is the first of its kind in the literature.
The charge deposition scheme is placed after the particle push and before the field solver.In this way, it has full access to the initial and final positions of each macro-particle, i.e. r n and r n+1 , as well as their half-step velocities, u n+1/2 .The first allows us to deposit the initial and final charge density in the grid nodes, i.e. ρ n i, j and ρ n+1 i, j are known, following Appendix B. Reconstructing the charged current density from the half-step velocities and the charge densities to satisfy equation ( 28) is the challenge we need to overcome.The solution to this problem was inspired in the seminal work by Villasenor and Buneman (1992) (hereafter VB) that suggested the computation of the charged current density directly from inverting equation ( 28), thus enforcing it by construction.As a particle moves from its initial position r n to its final position r n+1 , it may cross grid cell boundaries, hence continuously changing its macro-particle shape.To address this, we split the inter-cell motion into several intra-cell trajectories and tackle each split individually, as shown in Fig. 2. Without loss of generality, we will describe

splits
Figure 2: Schematic of two particle trajectories on a non-uniform grid.The particle to the left has an inter-cell motion of 3 splits, represented in 3 colors.The particle to the right is used as an example of the current deposition scheme for a single split.
how to retrieve the charged current density contribution from a single macro-particle given that its motion is bound to a single cell.The algorithm starts by splitting the divergence operator into its radial and meridional components: The VB method suggests that each component of ( 29) can be used to evaluate J r i+1/2, j and J θ i, j+1/2 by computing the temporal derivatives of the charge density when the particle moves purely along the radial (meridional) direction with an average meridional (radial) position.Here, we take a different approach due to the continuous change in the particle shape.We assume that the particle shape changes more significantly along the radial direction; thus, (∇ • J) θ cannot be evaluated at a fixed average radius.To overcome this issue, we compute it via with the right-hand side terms given by with θ = (θ n+1 + θ n )/2 being the average meridional position of the particle during the split motion.Notice that the poloidal components of the divergence could also be written as where V i, j is the volume of the cell evaluated at the cell node, see Appendix A. To justify the disappearance of the second terms in equations ( 33) and ( 34), we use the scheme presented in Fig. 2. Note that each split deposits the charged current density in the grid cell its trajectory is bound to.Finally, we obtain the poloidal components of the current by combining equations ( 32) with (33), and ( 30) with (34): Due to the axisymmetric condition, equation ( 28) does not restrict the azimuthal component of the current.Nevertheless, we can construct it via where the bared quantities correspond to a temporal average over the initial and final positions, e.g.
Also, recall that the total current density is the sum of the contributions of all macro-particles, which, in turn, is a summation over all the splits for each macro-particle.The current deposit algorithm presented above conserves charge to machine precision independently of the selected coordinates (polar spherical or logical spherical, see Appendix C).The main difference between both is, in fact, the choice of the cell center position, modifying the effective particle shape.The logical spherical coordinate system possesses a uniform Cartesian-like grid and asymmetric particle shape, thus solving the cell expansion issue characteristic of the symmetric particle shape and non-uniform grid case (i.e., polar spherical case).The latter case requires some particles to deposit charge and current densities in more than two consecutive cells (Cruz et al., 2023).

Boundary conditions
We now discuss the boundary conditions available in the OSIRIS-GR particle-in-cell framework.As the OSIRIS code allows for two possible domains (full or half meridional domain), the field and particle boundary conditions change accordingly.

Field boundary conditions
The meridional boundaries are of two kinds: axial or equatorial.The axial boundary condition ensures the axisymmetric condition.Due to the locations of each field component within a spherical Yee grid cell (see Fig. 1), only the E r, J r, B θ, E φ and J φ components lie in the boundary, hence need to be corrected: All field component values are mirrored and saved in the guard cells for usage in the field evolution equations.For example, B θ changes direction when crossing the axis, thus changing the sign in the guard cells and being set to zero on the axis, while B r does not invert its direction on the guard cells, hence being copied.
As for the equatorial boundary condition, the physical condition to be satisfied is the magnetospheric up-and-down symmetry.In this boundary, we mirror B r, B φ, E θ, and J θ, and copy the rest of the components.
The radial boundaries also take two options: star or open.The star radial boundary condition mimics the surface of a perfect rotating spherical conductor threaded by a dipolar magnetic field.We impose on the interior guard cells and on the surface components (i.e.E φ, E θ and B r) the interior electromagnetic solution of an isolated compact neutron star (equations ( 144)-( 147) in Torres et al. (2023)).We compute the guard cell values of the current components such that their values are progressively filtered (i.e. with an increasing filter order starting from 0 for the first radial cell to the desired order n for the (n + 1) th radial cell -see Appendix D for more details).
The open radial boundary, placed at the outer frontier of the domain, is a Mur-like absorbing boundary condition (Mur, 1981) designed to absorb outwardly propagating electromagnetic waves.This condition for an arbitrary field component Φ in flat spacetime reads: which is the Sommerfeld radiation condition in spherical coordinates (Novak and Bonazzola, 2004;Espinoza et al., 2014).
In Appendix E, we present the extension of equation ( 40) for curved spacetime and use it to correct the fields that lie in the exterior boundary (i.e.E φ, E θ, and B r).This open boundary condition was already used in Torres et al. (2023) and demonstrated efficient absorption of transient electromagnetic waves launched from a compact neutron star in vacuum.We take the same procedure for the charged current density as in the inner radial boundary.

Particle boundary conditions
Particle boundary conditions can be of two types: open or reflecting.In the radial direction, we employ open (or absorbing) boundary conditions on both domain frontiers.Particles that cross these boundaries are subtracted from the simulation.In the meridional direction, we use reflecting (or specular) boundary conditions.Particles are reflected back to the simulation domain with no energy loss.

Particle injection mechanism
OSIRIS allows for several plasma injection mechanisms.
Here we limit the discussion to the ones used in the simulations presented in upcoming sections.We distinguish two types of simulations, with or without heuristic pair production (Cruz et al., 2022).For the latter, we inject pair plasma from the stellar surface at a constant rate with number density n inj every timestep.Electrons and positrons are initialized in corotation and with a parallel velocity of 0.5 [c].This injection procedure is limited by the local value of the magnetization parameter σ * , given by: where B * is the local amplitude of the magnetic field, Γ is the Lorentz factor of the injected particles, and n ± is the number density of positrons or electrons.Setting a minimum value for σ * yields a maximum value of the injected plasma density.With the parallel component of the velocity, we can control the injected current provided to the neutron star's magnetosphere.Then, the magnetosphere extracts the sign of the charges required to drive the magnetospheric currents and sustain the global torsion of the magnetic field lines.The goal is to mimic the plasma supply provided by the vacuum work function and the pair-producing cascade close to the stellar surface, i.e. the polar cap gap, while still populating both the closed field line region and the return current.Previous works using global particle-in-cell simulations have already taken this approach as it is a computationally efficient way to populate the magnetosphere without requiring expensive quantum electrodynamics (QED) processes (Cerutti et al., 2015;Philippov et al., 2015b).Nevertheless, these processes are already implemented in the OSIRIS framework, to model QED cascades in future ultraintense laser experiments (Grismayer et al., 2016(Grismayer et al., , 2017) ) or in local neutron star polar caps (Cruz et al., 2021).The second approach is to sustain a cold corotating atmosphere close to the surface, similarly to Hu and Beloborodov (2022); Bransgrove et al. (2022).We characterize the atmosphere by its value at the surface, n 0 , and the scale height of the exponentially decaying profile, σ atm : We ensure the non-depletion of this atmosphere up to k atm standard deviations from its base by injecting n inj plasma number density every timestep when required.In the current study, we consider a plasma composed of electrons and positrons, the injection of ions is deferred to future studies.This approach differs from the previous one as it allows the magnetosphere to regulate the flux of plasma extracted from the corotating atmosphere.This extraction mechanism would then seed the cascading QED process.To reduce the computational cost of each simulation, we model the cascade with a heuristic pair production (Philippov et al., 2015a,b;Philippov and Spitkovsky, 2018;Chen et al., 2020;Guépin et al., 2020;Cruz et al., 2022Cruz et al., , 2023)).This model for the operation of the vacuum gap defines a threshold energy γ thr for a particle to pair produce.Particles that satisfy this condition can pair produce if their radial position corresponds to the range R * ≤ r ≤ R PP .Also, this threshold energy has a polar dependence to inhibit pair production close to the polar axis, where the magnetic conversion does not occur due to the infinite magnetic field line curvature.The secondary pair is emitted at the same place and with the same momentum direction as the parent with a Lorentz factor of γ pair /2.These parameters are user-defined and can be adjusted to study different magnetospheric solutions.This model has shown to be a good model for the polar cap and the outer gaps (Philippov et al., 2015a,b;Philippov and Spitkovsky, 2018;Chen et al., 2020;Guépin et al., 2020;Cruz et al., 2022Cruz et al., , 2023)).

Simulation initialization and transient phase
Simulations start in vacuum with a prescribed generalrelativistic dipolar magnetic field (Rezzolla et al., 2001;Torres et al., 2023) threading the neutron star surface.Equivalently, the angular velocity of the star can be either in full rotation or gradually spun to its desired value.For the first scenario, the general-relativistic electric field (Torres et al., 2023) must be prescribed at t = 0.In the second scenario, the rotation imparted by the stellar boundary condition (see Sect. 2.6.1)generates the self-consistent electric field by launching Alfvénic torsion waves (Cerutti et al., 2015).At t = 0, we initiate the plasma injection mechanism described in Sect.2.7.Due to the high magnetization, the extracted charges from the surface/atmosphere are accelerated along field lines, generating a poloidal current that opens up the field lines that cross the lightcylinder radius, located at R LC ≡ c/Ω * .Beyond this distance, particles can no longer co-rotate with the field lines anchored to the stellar surface as they would need to co-rotate at speeds faster than the speed of light.These particles then sustain the torsion of the field lines that acquire a toroidal component.In the aligned rotator case, when the magnetic and rotation axes are aligned, the extracted current at the poles is negative due to the predominant outflow of electrons.The magnetosphere selforganizes to close the electric current, which generates the return current.This positive electric current neutralizes the charge density of the neutron star and flows between the closed and open field lines in a Y-shaped structure.Beyond the Y-point, the current sheet sustains the toroidal and poloidal magnetic field reversal at the equator.The magnetospheric solution starts to form in the back of the outwardly propagating Alfvénic torsion wave.We have assumed that after two rotation periods of the star, the magnetospheric solution has converged to the quasi-steady state.

General relativity and the polar cap emission
Coherent radio emission from the polar caps of a neutron star enabled the discovery of pulsars.The pulsed emission profile originates from the cosmic lighthouse effect, which relies on the misalignment angle χ between the magnetic and rotation axis.Recent works have demonstrated the generation of pulsar radio emission from non-stationary pair plasma discharges (Philippov et al., 2020;Cruz et al., 2021).Rotation-induced vacuum gaps, regions of an unscreened parallel electric field, potentiate the acceleration of particles to very high energies and subsequent emission of gamma-ray photons via curvature radiation.These hard photons are absorbed in super-strong magnetic fields and create pairs, generating a QED cascade that populates the magnetosphere with a dense electron-positron plasma.The pair production bursts are generated at an angle to the local magnetic field, favouring the induction of electromagnetic modes (Cruz et al., 2021).However, these models rely on efficient particle acceleration along the polar cap field lines for QED processes to kick in.Beloborodov (2008) and Timokhin and Arons (2012) addressed the particle accelerator problem at the polar cap.These works identified the α ≡ j m / j GJ parameter determinant in defining under which conditions particles accelerate up to relativistic energies and can produce secondary pairs.The motivation behind this parameter, which is the ratio between the electric current extracted from the polar cap and the local co-rotation current j GJ = −Ω * • B/2π, is that when 0 < α < 1, the accelerator is inefficient because the extracted current is enough to sustain the twist of the magnetic field lines near the light cylinder, given by When α > 1 or α < 0, charges extracted from the surface/atmosphere are insufficient to carry the magnetospheric required current.This current-starved state generates an electric field that accelerates particles.This flow becomes unstable to QED cascading, thus supplying the additional carriers.Under such conditions, the polar cap becomes an efficient accelerator, and pulsar radio emission is viable.However, threedimensional particle-in-cell simulations of low obliquity rotators (i.e., χ ≤ 40 • ) showed weak particle acceleration in the bulk of the polar region (Philippov et al., 2015b), confirming estimates using the α parameter (Bai and Spitkovsky, 2010;Timokhin and Arons, 2012).Therefore, low obliquity rotators could not generate pulsar radio emission, which is unsupported by observations.The aligned rotator (i.e., χ = 0 • ) constitutes the least favourable configuration for the pulsar mechanism occurrence.We can use this configuration as a robustness diagnostic for any other model that could explain the radio emission.Here we discuss the general-relativistic effects and adopt the aligned rotator configuration.
In the 3+1 formalism, we reformulate the α parameter analysis for the efficiency of the polar cap accelerator as the 4-norm of the current j µ j µ .Equivalently, we have an efficient accelerator for j µ j µ > 0 (i.e., spacelike current) or inefficient accelerator for j µ j µ < 0 (i.e., timelike current) (Belyaev and Parfrey, 2016;Gralla et al., 2016;Huang et al., 2018).In this way, we can use the spacelike current as an observable for field lines that will sustain pair creation and are viable locations for the pulsar radio emission.Therefore, characterizing the spacelike current region may lead to the direct characterization of the emitted polar radio beam, bridging the gap between numerical models and observations.
The wind region of the magnetosphere (the volume of open field lines that extend beyond the light cylinder radius) is somewhat similar for both the dipolar and split-monopolar magnetic field configurations near the rotation axis.The main difference is the existence of a closed field line region at the equator which reduces the number of open field lines crossing the light cylinder.As we are interested in the polar cap of the neutron star magnetosphere, we will neglect this effect and assume, hereafter, that the field configuration resembles that of a split monopole.Taking small angles from the axis, we can approximate the electric current as solely pointing in the radial direction.Hence, the 4-norm of the electric current is given by where we assumed the magnetospheric solution is close to the force-free regime where the charge density closely follows the co-rotation value ρ ∼ ρ GJ and the radial current approximates to the ratio between the azimuthal magnetic field and the polar component of the electric field, as obtained by Lyutikov (2011).Although simple, equation ( 44) is powerful enough to provide us with intuition on the interplay between magnetospheric plasma injection and general-relativistic effects.The first is made explicit through B φ that measures the torsion of the field lines, which is more efficient in the force-free regime where a strong poloidal current gets extracted from the polar cap.The second point relates to the frame-dragging effect (Beskin, 1990;Muslimov and Tsygan, 1992;Philippov et al., 2015a;Torres et al., 2023), which reduces the induced electric field due to the mismatch between the stellar and the spacetime rotation.In the force-free regime, the wind region possesses B φ ∼ E θ, corresponding to the null electric current case.This case is not particularly interesting as this would mean that, in flat spacetime, the aligned rotator would not be able to emit in the radio band.However, the frame-dragging effect reduces E θ raising the ratio above unity and leading to the possibility of pulsar radio emission (Philippov et al., 2015a).
To validate this intuitive picture drawn from equation ( 44), we simulated the neutron star magnetosphere with no pair production and injected plasma from the surface with parallel velocity as described in Sect.2.7.This injection method has two advantages: (1) direct control over the extracted poloidal current limited by the cold magnetization parameter; (2) successfully reproduces the force-free magnetospheric solution.Simulations ran in the half-domain setup for R s /R * ∈ [0.0, 0.6] and σ * ∈ [1000, 2000], presented in Fig. 4. As the maximum meridional extension predicted by force-free models for the spacelike current region (Belyaev and Parfrey, 2016;Gralla et al., 2016), θ SL , is a function of the stellar compactness, the meridional resolution changed accordingly for each simulation to resolve this region with at least 20 numerical cells.This angle also changes with the stellar angular velocity, fixed at Ω * = 0.1 [c rad/R * ] for this set of simulations.Table 1 contains the detailed list of parameters used.It is important to reiterate that our study takes into account all general-relativistic effects, not just the frame-dragging, unlike previous works on neutron stars (e.g., Philippov and Spitkovsky (2018) grove et al. ( 2022)).In fact, as is shown in Fig. 3, this approximation leads to an overestimation of both the polar cap and spacelike angles of ∼ 21.8% for R s /R * = 0.5 (typical stellar compactness selected), following the general-relativistic forcefree prediction given in Gralla et al. (2016).This overestimation leads to macroscopic magnetospheric changes, including an increased effective open flux tube area, spin-down luminosity, and wider generated radio beam.This point is particularly important when simulating magnetospheric environments of compact neutron stars approaching the Tolman-Oppenheimer-Volkoff mass limit (Kalogera and Baym, 1996;Romani et al., 2022) or long-lived hypermassive neutron stars (Falcke and Rezzolla, 2014;Chirenti et al., 2019;Chirenti et al., 2023;Selvi et al., 2024).The first striking feature in Fig. 4 is the inexistence of the spacelike region for some non-zero compactness simulations.The explanation lies in the poloidal current and, consequently, the degree of field line torsion.When σ * = 1000, the extracted poloidal current is very close to the expected for the force-free split-monopole case, thus ensuring B φ ∼ E θ.In this way, any general-relativistic reduction of the electric field results in a spacelike current, following equation ( 44).Lower plasma supply cases represented by higher magnetizations yield a decrease in the poloidal current extracted from the polar cap.The B φ reduction no longer guarantees the existence of the spacelike region at any compactness, thus explaining the appearance of cut-off conditions as shown in Fig. 5. Figure 5 highlights the first direct comparison between global particle-in-cell simulations and force-free model predictions for the meridional extension of the spacelike current volume.Results show a remarkable agreement in the force-free limit but start to deviate for more charge-starved magnetospheric solutions.Nevertheless, the force-free estimates limit the predicted spacelike extension from above, thus it is still a good estimate for the maximum generated radio beam width.
From the observational point of view, this shows the robustness of the pulsar radio emission mechanism at any inclination.We have just shown that the spacelike region appears at any compactness for force-free magnetospheric solutions, even for the least favourable inclination configuration (the aligned rotator).Nevertheless, the dependence of the radio beam generation mechanism on the overall magnetospheric solution could explain why some neutron stars are radio-quiet or may present intermittent behaviour.Weak pulsars or pulsars near the death line, where pair production is no longer very efficient, could be strong candidates for this scenario.Also, our model predicts narrower radio beam widths for older low-obliquity pulsars, whose magnetospheres are more charge-starved.In some cases, this width reduction can reach half the force-free value, which may be measurable through observations.In particular, observations of millisecond pulsars report a systematic tendency for smaller inferred luminosities and narrower emission beams (Kramer et al., 1998), which support our findings.Thus, the compactness of millisecond pulsar magnetospheres and deviations from force-free solutions may explain the reduction of the angular radius of the generated radio beams.

Transition from weak to energetic pulsars
In the previous section, we validated the intuitive picture given by equation ( 44).However, we adopted a controlled injection methodology that drove the system to the desired solution.In this section, we improve the magnetospheric model by capturing the gap dynamics through heuristic "on-the-spot" pair production (Cruz et al., 2022), i.e. in the limit of vanishing photon mean free path.Under these circumstances, solutions become sensitive to the upstream gap plasma supply from the neutron star surface.The charges lifted from the surface form an atmosphere maintained as described in Sect.2.7, feeding selfconsistently the magnetosphere with cold plasma.We follow the PSG non-stationary gap model by sustaining a minimum atmospheric plasma density at the base n 0 ≡ n + + n − ∼ 2n GJ .A weak parallel component of the electric field persists at the base of some magnetic field lines, providing an initial kick to those particles extracted from the atmosphere.
Figure 6 shows the quasi-stationary state of a millisecond neutron star with compactness R s /R * = 0.3 (e.g., Bhattacharyya et al. (2017)) and angular velocity Ω * = 0.2 [c rad/R * ] for different pair production efficiencies η given by: where γ max ∼ B * Ω 2 * is the estimated maximum achievable particle Lorentz factor.The magnetic field amplitude controls the pair production efficiency by keeping constant the threshold and secondary pair energy at γ thr = 25 [m e c 2 ] and γ pair = 8 [m e c 2 ].Also, we restrict pair production to be active only up to three stellar radii (i.e., R PP = 3R * ).Table 2 summarizes the simulation parameters used in this numerical experiment.The low pair production efficiency case, with η ∼ 38, highlights the typical solution of a weak pulsar with a persistent outer gap (Shinya and Shinpei, 2012;Gruzinov, 2012Gruzinov, , 2013Gruzinov, , 2015;;Chen et al., 2020).This gap is a consequence of the inability of the surface and the pair production in the volumetric return current to supply plasma to the region behind the null surface, i.e. the region where the poloidal magnetic field lines curve towards the equator.As the efficiency increases, the amplitude of the gap's parallel electric field reduces until it becomes almost zero, when the magnetosphere transitions to a force-free solution.The volumetric return current that occupied a portion of the polar cap for low-efficiency values now occupies lower latitudes, increasing the effective polar cap area up to the maximum poloidal extension angle θ PC , as predicted by force-free models (Gralla et al., 2016;Belyaev and Parfrey, 2016).Consequently, the magnetosphere can increase the electric current driven through the open flux tube, leading to a jump in the Poynting flux luminosity.In this model, and for simplicity, we kept the pair production efficiency constant in the magnetosphere.However, different QED processes play a role in distinct magnetospheric regions (Chen and Beloborodov, 2014), which may lead to differential pair production efficiencies (deferred for future work).Magnetospheres that continuously oscillate between weak and quasiforce-free states, depending on the morphology of the outer gap, can explain the observed phenomenon of pulsar luminosity intermittency (similar to Li et al. (2012)).
Several crucial changes accompany the transition from weak to force-free solutions: 1) the Y-point moves closer to the lightcylinder radius due to the more efficient screening of the equatorial superrotation region (Shinya and Shinpei, 2012;Hu and Beloborodov, 2022); 2) the stronger poloidal current driven through the magnetosphere induces a higher torsion of the magnetic field lines in the toroidal direction.The latter is vital for the activation of the polar cap.At fixed compactness (in this case R s /R * = 0.3), Fig. 5 shows that the polar cap may present a spacelike volumetric region only if the magnetosphere drives a stronger current, characteristic of higher plasma supply scenarios.With heuristic pair production, the efficiency controls the plasma supplied from the polar cap gap to the magnetosphere, thus having the same effect as the surface plasma supply from the previous section.In Fig. 6, we observe the appearance of the spacelike volume in the polar cap for efficiencies η ≳ 155.However, many processes affect this minimum efficiency value η min : 1) the frame-dragging effect for higher compactness neutron stars reduces η min due to the reduction of the poloidal electric field; 2) evident interplay between the outer and inner polar cap gap, more efficient pair production at the Y-point (e.g., including photon-photon processes) pulls the return current away from the polar cap thus increasing the field line torsion of the magnetosphere, which reduces η min ; 3) a denser neutron star atmosphere also facilitates the launch of a stronger poloidal current, in a similar way as in point 2. In particular, point 2 suggests a correlation between the radio beam and the existence of the outer gap.Intermittent pulsars may also present intermittent radio beams, potential candidates for explaining pulse nullings.We have also performed simulations where we inject an atmospheric plasma composed of electrons and ions, where the ions have the same mass and charge as positrons but are not allowed to pair produce (Philippov and Spitkovsky, 2018).As     expected, the injection of ions leads to a lower pair plasma supply to the outer magnetosphere of weak pulsars, thus increasing the minimum efficiency value η min .Apart from a shift of the magnetospheric solution towards charge-starvation, we do not observe significant changes in the outer gap dynamics or magnetospheric structure.However, this point may change if ions possess realistic q/m, which we differ to a future study.Also, the impact of ion injection on the magnetosphere is more significant in charge-starved pulsar solutions due to the high degree of charge separation, thus this effect is diluted for denser solutions approaching the force-free regime.
Solutions with very low pair production efficiencies are sustained solely by pair plasma generated at the return current region and progressively tend to the "dead pulsar" solution (electrosphere), corresponding to the no pair production case (Cruz et al., 2023).Solutions with even higher pair production efficiencies than the ones presented in this study tend to the forcefree regime and possess an active polar cap, i.e. in the accelerator regime.In this sense, we retrieve the same conclusions as in the previous section: low pair production efficiencies (i.e., low plasma supply) require higher stellar compactnesses for the activation of the pulsar mechanism; higher pair production efficiencies (i.e., high plasma supply) relax the compactness restriction.Nevertheless, for low-obliquity rotators, compactness is crucial in the generation phase of the radio beam along the magnetic axis.Also, as observed in the previous section, the width of the radio beam deviates from the force-free estimate depending on the charge-starvation state.Thus, for older and weaker low-obliquity pulsars, the radio beam width should be reduced.
Another important feature on the structure of the magnetosphere is the location of the Y-point, which for all simulations presented is located at a significant fraction of the light-cylinder radius (i.e.r Y /R LC ∼ 0.6 − 0.8), consistent with previous ki-netic models of neutron star magnetospheres (Chen et al., 2020;Guépin et al., 2020;Cruz et al., 2023;Hakobyan et al., 2023).The spatial restriction imposed on the heuristic pair production (R PP = 3R * ) limits the supply of plasma to the outer magnetosphere, in particular, to regions close to the Y-point.This weak supply of plasma beyond R LC explains the formation of the persistent outer gap, volumetric return current and wide equatorial current sheet.This is in contrast to the typical ideal force-free magnetosphere with abundant plasma supply, clear and localized Y-shaped return current and thin equatorial current sheet.Wider current sheets dissipate energy at a much slower rate, thus we expect our novel reduced-domain model (i.e., imposing equatorial symmetry) to be accurate in modelling weak pulsar magnetospheres in the aligned rotator configuration.We have also verified the convergence of the presented results by doubling the number of particles per cell injected every timestep.
Particularly interesting is the systematic formation of core component-associated beamlets at lower pair production efficiencies (see positronic charge density for η = 38 − 77 simulations in Fig. 6).These beamlets are generated by the cyclic opening of a fragile gap near the magnetic pole of the star.As pair production is prohibited close to the axis due to the small curvature of the magnetic field lines, some electrons from the magnetosphere are reversed back to the star, thus intermittently shorting out the accelerating gap from above (Lyubarsky, 2009).This cyclic behaviour allows the magnetosphere to drive the required quasi-steady-state polar current that supports the magnetic field line torsion in the outer magnetosphere.We conjecture that these beamlets are responsible for the narrow dwarf pulses observed during ordinary pulse nullings of PSR B2111+46 in Chen et al. (2023), referred to as "particle raindrops" by the authors.

Conclusions
Understanding neutron star magnetosphere dynamics is of utmost importance to reveal the locations of efficient particle acceleration and radiation generation.The connection between the rich pulsar observations and theoretical models is only possible through global ab initio simulations.In particular, due to the extreme nature of such astrophysical environments, kinetic models must capture self-consistent QED processes and general relativistic effects.This paper introduces a new generalrelativistic module for the particle-in-cell code OSIRIS (Fonseca et al., 2002), intended to study plasma dynamics in intense electromagnetic and gravitational fields, capturing selfconsistently GR effects in all components of the PIC algorithm (field solver, particle pusher, charge and current deposit scheme).
We designed two numerical experiments to gain a deeper insight into the role of general relativity in activating the polar cap, focusing on the conditions for efficient particle acceleration.Therefore, we used the 4-norm of the current to probe under what circumstances the emission of a coherent radio beam along the magnetic axis is viable.In the first experiment, we injected plasma at 0.5 [c] along the magnetic field lines, controlled by the local magnetization, to reproduce different quasi-stationary magnetospheric states.This setup allowed for a structured parametric scan over the plasma supply and stellar compactness.Without the general-relativistic effects, a low obliquity rotator would not be able to accelerate particles efficiently enough to potentiate the radio beam generation independently of the plasma supply.The general relativistic framedragging effect reduces the poloidal electric field, thus increasing the current amplitude.We show that this GR effect activates the polar cap of young and energetic neutron stars at any nonzero compactness, thus providing a robust mechanism for radio beam generation.For older neutron stars, it depends on both the compactness and plasma supply.Nevertheless, we show that GR could enable radio emission even for weak pulsars that drive weaker poloidal currents.Therefore, weak pulsars may be observable in the radio spectrum beyond their expected death line, some of which can present themselves with narrower beams concerning force-free estimates.We note, however, that measuring such deviations is challenging due to other competing effects, such as the observer's line of sight limitation or the misalignment angle χ, requiring a full three-dimensional model extension.Previous kinetic studies of pulsar magnetospheres at arbitrary inclinations and low plasma injection rates demonstrated a similar behaviour of the magnetospheric solution (Kalapotharakos et al., 2018;Brambilla et al., 2018), with decreased Poynting flux and an outer gap beyond the null surface and around the separatrix region.Consequently, the volumetric field-aligned magnetospheric current and the torsion of the magnetic field lines decrease for these weak pulsar solutions.Following the conclusions from this work, we expect a reduction of the spacelike volume ( j µ j µ > 0 or α > 1) in the outflowing region of the polar cap (Bai and Spitkovsky, 2010;Philippov and Kramer, 2022), where pair production and the generation of the radio beam are expected.However, this effect should be more prominent in low-obliquity or very chargestarved pulsars, where this deviation from the force-free solution produces clear differences in the expected beam width.Our findings and conjecture seem to agree with observations of millisecond pulsars, which do not follow the scaling predicted from a canonical pulsar model for the inferred luminosity and angular radius of the generated radio beam (Kramer et al., 1998).The compactness of their magnetospheres and deviations from the force-free solution may explain such systematic tendencies.
In the second experiment (Sect.4), we improve the accuracy of the previous model by resolving the polar cap gap dynamics using a heuristic "on-the-spot" pair production model.Here, we allow pair production to populate only regions of efficient particle acceleration.The magnetosphere self-regulates to drive the required currents that sustain the magnetic field line torsion.By varying the pair production efficiency, i.e. through an increase in the magnetic field amplitude, we observe a clear transition between the weak and force-free pulsar solutions at fixed compactness (fiducial value of R s /R * = 0.3).Once again, young and energetic neutron stars support strong poloidal currents, characteristic of a denser force-free magnetospheric solution, and possess active polar caps at any non-zero stellar compactnesses.Older pulsars, i.e. pulsars approaching their death line, transition to weak pulsar solutions and must rely on both the -10 0 -10 -1 0 10 -1 10 0 10 1 , , -10 -1 -10 -2 0 10 -2 10 -1 0.5 -0.5  pair production efficiency (plasma supply) and stellar compactness for the presence of efficient particle acceleration and consequent radio beam.In Fig. 6, we demonstrate the existence of a minimum pair production efficiency for which the spacelike current develops in the bulk of the polar cap.The previous numerical experiment also captured this transient behaviour of the spacelike current via varying the plasma supply at fixed compactness (see Fig. 5).Also, simulations highlight the poloidal current-mediated interplay between the polar cap and outer gaps.Magnetospheres with such null surface gaps develop a volumetric return current, constricting the polar cap area and limiting the poloidal current in the open magnetic field line bundle.Consequently, equation ( 44) requires a lower poloidal electric field to activate the polar cap discharge, achieved by increasing the stellar compactness.Magnetospheres with pair production near the Y-point populate these exterior regions with plasma, relaxing the polar cap poloidal extension and pair production efficiency at the poles.This interconnection between outer and inner gaps may explain the intermittency of certain observed pulsars or even justify the existence of pulse nullings.
In particular, weak pulsar solutions generate low-multiplicity plasma beamlets near the magnetic axis, compatible with the observed narrow-width dwarf pulses in Chen et al. (2023).
In this work, we have employed a novel reduced-domain approach to model the global magnetosphere of a neutron star.This reduced model allows for simulations with higher spatiotemporal resolution and better macro-particle statistics, at the expense of introducing an up-down equatorial symmetry.Consequently, by construction, this model does not capture the development of the kink instability in the equatorial current sheet (Cerutti et al., 2015;Philippov et al., 2015a).In that sense, we do not expect our method to model correctly the dissipation properties of thin current sheets, characterized by frequent magnetic reconnection events and plasmoid formation (Hakobyan et al., 2019(Hakobyan et al., , 2023;;Schoeffler et al., 2023).However, we expect good agreement for weak pulsar solutions with no pair production near R LC , which possess a volumetric returncurrent and wide equatorial current sheet (with no significant dissipation/reconnection events occurring).Also, the location of the Y-point is consistent with similar full-domain particlein-cell simulations (Chen et al., 2020;Guépin et al., 2020;Cruz et al., 2023;Hakobyan et al., 2023), indicating that this reduced model preserves the main elements of a neutron star magnetosphere in the aligned rotator configuration.
Future works with more accurate pair production models, i.e. heuristic models including non-zero photon mean free path and geodesic propagation, or even modelling QED processes more accurately with Monte Carlo techniques (Cruz et al., 2021), should explore this global outer-to-inner gap dynamics.Also, the variation of the upstream plasma supply from the stellar surface/atmosphere may affect the dynamics of the polar cap gap and its dependency on the driven magnetospheric current, requiring further investigation.In particular, with the surface/atmospheric injection of electron-ion plasma instead of pair plasma, which affects the discharge dynamics on the return-current, modifying the plasma supply to the outer magnetosphere.Also, three-dimensional models of weak pul-sar magnetospheres for low-obliquity rotators may unveil the importance of GR in the pulsar observational appearance.
The volume element can also be evaluated via which is used in the charge conservation scheme of Sect.2.5.The radial part of the volume integral is given by The same procedure can be repeated for any diagonal spatial metric to retrieve the metric elements of other metric or coordinate systems.
In Sect.2.4, the spatial Christoffel symbols used in equation ( 20) for the logical spherical coordinate system are given by:

Appendix B. Charge deposition scheme
In the 3+1 formalism, we can retrieve the charge conservation condition via the divergence of the Ampére's law: using Gauss's law (5) and the vector calculus property that states that the divergence of a curl of a vector field is always zero.
For the sake of clarification, we start by defining that each macro-particle, located at (r p , θp ), has the same shape as the cell it is in, which is centered at (r i+1/2 , θ j+1/2 ).Each macroparticle occupies a volume given by the particle shape function S (r, θ, rp , θp ), which we will define later.The number of particles in each macro-particle, N p , is given by where n(r, θ) is the particle number density, which we will assume to be constant within each cell (i.e.waterbag-like).In this way, we can obtain the particle number density in each grid cell: The charge density at any position due to a macro-particle located at (r p , θp ) with N p particles of charge q p reads ρ p r, θ, rp , θp = q p N p S r, θ, rp , θp . (B.9) Consequently, the contribution of the macro-particle's charge density to the charge density evaluated at the grid nodes is defined through the volume weighting technique as where S r(r p ) and S θ( θp ) are spline functions given by: , (B.12) specifying the limits taken for each integral function between square brackets.The new integration limits r< = max(r min , ri−1/2 ), θ< = max( θmin , θi−1/2 ), r> = min(r max , ri+1/2 ), and θ> = max( θmax , θi+1/2 ) arise from the integral of the b 0 functions present in equation (B.8).The total charge density at the grid position (i, j) is the sum of (B.10) over all the particles that contribute, i.e.We demonstrate that the code conserves charge to machine precision independently of the selected coordinate system.To exemplify, we use the standard spherical coordinate system with a non-uniform grid spacing in both radial and meridional directions, i.e. uniform spacing in log r and − cos θ, and initialize particles with a thermal distribution of velocities in all directions, i.e. u thermal /c = (2.0,2.0, 2.0), in the vicinity of a rotating neutron star with dipolar magnetic field B * = 1000 [m e c 2 /eR * ] and Ω * = 0.1 [c rad/R * ], see upper-left panel in Fig. C.7.Also, we selected a compactness parameter R s /R * = 0.5 to show that the scheme is valid for significant curved spacetimes.Every iteration, we compute the deviation from the continuity equation given by

Appendix D. Progressive filter
We explain in more detail the implementation of the generic progressive filter applied to the first and last radial domain cells.In particular, this progressive filter is used to filter the current density grid components, in this appendix designated f k i with i and k being the cell index and the filtering order, respectively.Outside the simulation domain, the current density is unknown, thus invalidating high-order smoothing near the boundaries.Instead of not filtering the boundary cells or filling the guard cells with unphysical field values, we compute the exterior current component values such that only the first and last cells remain unfiltered.
Within the OSIRIS framework, smoothing any field component to order n requires a single-passage kernel with 2n + 1 en-tries.For example, the single-passage first-, second-, and thirdorder kernels read with K being the weight given to each neighbour.The thirdorder filtered value of the first four cells yields where the terms in red are the unknown field values inside the guard cells.Equations (D.4)-(D.7)highlight that we can filter the fourth cell using the third-order single passage filter in eq.(D.3), the third cell using the second-order single passage filter in eq.(D.2), the second cell using the first-order singlepassage filter in eq.(D.1), and the first cell remains unfiltered, i.e.
All other field values are filtered using the third-order single passage kernel.

Appendix E. Mur open boundary
We discuss the implementation of the exterior radial boundary condition that mimics an open radiation boundary.We define the generic wave equation in spherical coordinates □Φ (t, r, θ, ϕ) = σ (t, r, θ, ϕ) ,   where σ is the source term and □ is the d'Alembertian operator.Since we place this absorbing boundary at the extreme of the domain, i.e. r = R max with R max ≫ R * , we can neglect the effect of the source term.Also, electromagnetic perturbations that propagate radially outwards are dominantly radial due to the geometry and dynamics of the system, i.e. ∂ θ Φ ≈ 0. With axisymmetry (∂ ϕ Φ = 0), we can neglect the d'Alembertian operator for the angular derivatives: The radiation condition impinging on a spherical shell placed far from the compact point source derived in (Sommerfeld, 1949) reads lim r→∞ (∂ r + ∂ t ) (rΦ) = 0, (E.3) describes only the outward propagating component of equation (E.2) and is hereafter designated as the Sommerfeld condition.At a finite distance from the source, it approximates to (Novak and Bonazzola, 2004;Espinoza et al., 2014) which is exact for pure monopolar waves.
In General Relativity, using the 3+1 formalism in spherical coordinates, the modification of the differential operators leads to the generalization of equation (E.4) as α −1 ∂ t Φ + α∂ r Φ + Φ r r=R max = 0, (E.5) which we discretize adopting a spatiotemporal centering at (r, t) = (r N+1/2 , t n+1/2 ): + Φ n+1/2 r r N+1/2 = 0, (E.6)where N and N + 1 are the radial indices of the last cell of the domain and the first guard cell outside of it, respectively.The field components that lie in the outer grid frontier, i.e. at r = r N+1 , are E φ, E θ, and B r, which means that equation (E.6) needs to be adapted for each of these components.Recalling that we are using the leapfrog method to push the electromagnetic field components in time, we see that we have the electric field components at times n and n + 1, while the magnetic field components are available at times n and n+1/2, for the first half push, and n + 1/2 and n + 1, for the second half push.Consequently, the algorithm to compute the correct field value outside the domain, i.e.Φ n+1 N+1 , differs for the electric and magnetic field components.

Appendix E.1. Electric field
Here we discuss the algorithm with Φ representing E φ and E θ.The first term in equation (E.6) is already in the right temporal positions but requires the correct spatial centering: In opposition, the second term is in the right spatial position and requires the correct temporal centering: The last term in equation (E.6) is more complicated because we need to center it in time and space: (E.9) Inputting equations (E.7)-(E.9)into (E.6) and solving for the corrected field outside the domain, we get: Notice that to compute this, we require knowledge of the term Φ n+1 N that is available on the grid, and the knowledge of the terms Φ n N+1 and Φ n N that need to be saved in memory as, otherwise, they would be overwritten on the grid.

Appendix E.2. Magnetic field
The approach taken to correct the out-of-bound B r is very similar but it must be split into two steps, as in the magnetic evolution algorithm.Therefore, we need to compute Φ n+1/2 N+1 (Φ n+1 N+1 ) after the first (second) half-step.To obtain these expressions, we need to repeat the same procedure as for the electric field with the extra step of adding and subtracting the half-step value, i.e. equation (E.7) reads (E.11) which we repeat for equations (E.8) and (E.9).Grouping terms that depend on the temporal indices n and n+1/2 yields the first half-step correction: (E.12)

Figure 3 :
Figure3: Force-free estimates of the maximum polar cap, θ PC , and spacelike, θ SL , poloidal extension at the neutron star surface for varying stellar compactness values, R s /R * .Curves are obtained using the model developed inGralla et al. (2016) for a general-relativistic description of force-free electrodynamics.The complete description is presented in blue.Estimates neglecting the framedragging contribution and the lapse function are presented in green and orange, respectively.We conclude that not considering the lapse function leads to an overestimate of the open flux tube and effective particle acceleration volume, thus a macroscopic inaccuracy.

Figure 4 :
Figure 4: Temporally averaged 4-norm of the electric current for a combination of different stellar compactness and cold magnetization values.This diagnostic allows for the identification of field lines where efficient particle acceleration is expected.The red lines represent the magnetic field lines that start at the maximum predicted angle for the spacelike region (θ SL ) and the polar cap (θ PC ) (as predicted by Gralla et al. (2016)).The grey lines represent the magnetic field lines.The results show efficient particle acceleration at any compactness for σ * = 1000, i.e. high plasma supply case.Also, at lower plasma supply, particle acceleration occurs only after a certain value of compactness, e.g.R s /R * ≳ 0.4 for σ * = 1500.

Figure 5 :
Figure 5: Maximum poloidal extension of the spacelike region (θ SL ) as function of the stellar compactness R s /R * .The grey line corresponds to the force-free estimate in Gralla et al. (2016).The circles correspond to the maximum and minimum measured angles extracted from numerical simulations presented in Fig. 4.This figure captures the transition from inactive to active polar caps with increasing compactness, e.g. for the σ * = 1500 (2000) case, the polar cap activates for R s /R * ≳ 0.4 (0.5).

Figure 6 :
Figure6: Panels showing the positronic charge density (ρ + ), 4-norm of the current ( j µ j µ ), FIDO-measured radial current (J r ), and parallel component of the electric field (E ∥ ) for increasing pair production efficiency (η).Field components are multiplied by the radial distance squared (r 2 ) to highlight magnetospheric features further away from the star.Magnetic field lines are in grey and the magnetic field lines that start at the maximum spacelike angle (θ SL ) and polar cap extension angle (θ PC ) are shown in red.This figure captures the transition from a weak to force-free solution, obtained by increasing the stellar magnetic moment and, hence, the pair production efficiency.Consequently, the Y-point (at r/R * ∼ r Y ) approaches the R LC , the polar cap extends to θ PC , the parallel electric field component almost vanishes everywhere, and the spacelike current and positrons from pair production appear in the polar cap (see j µ j µ zoom insets).The location of the low-multiplicity beamlets is highlighted in the positronic charge density insets close to the poloidal axes.
(C.2) and evaluate the grid average value of that quantity, shown in the right panel of Fig. C.7.Under the same conditions, but now using the logical coordinate system, i.e. using a uniform grid in (r, θ) = (log r, − cos θ), yields the results in Fig. C.8. Apart from small numerical fluctuations, the grid averaged deviation of the continuity equation stays close to machine precision over the test duration for both cases.

Figure C. 7 :
Figure C.7: Charge conservation diagnostic for a randomized sample of electrons initialized with thermal velocities.This setup uses nonlinear grid spacings in both r (logarithmic) and θ (equal area) to highlight that the proposed scheme works even when not using logical coordinate systems.Left panels show the initial and final states, the right panel displays the temporal evolution of the spatial averaged diagnostic.

Figure
Figure C.8: Charge conservation diagnostic for a randomized sample of electrons initialized with thermal velocities.This setup uses uniform grid spacings in both r and θ, i.e. using the logical coordinate system.Left panels show the initial and final states, the right panel displays the temporal evolution of the spatial averaged diagnostic.

Table 1 :
Simulation parameters used for the first numerical experiment (Sect.shown in Figs.

Table 2 :
Simulation parameters used for the second numerical experiment (Sect.4),shown in Fig.6.These simulations take on average 500k CPU core hours.