<figure><img src="image_1.png" /></figure>

HAL
open science

# Radiative pulsar magnetospheres: oblique rotators

J Pétri

► To cite this version:

J Pétri. Radiative pulsar magnetospheres: oblique rotators. Monthly Notices of the Royal Astronomical Society, 2022, 10.1093/mnras/stac572. hal-03607664

HAL Id: hal-03607664

https://hal.science/hal-03607664

Submitted on 31 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_2.png" /></figure>

# Radiative pulsar magnetospheres: oblique rotators

J. Pétri

<figure><img src="image_3.png" /></figure>

Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

Accepted 2022 February 24. Received 2022 February 24; in original form 2021 May 5

## ABSTRACT

Pulsar magnetospheres are filled with relativistic pairs copiously emitting photons detected from the radio wavelengths up to high and very high energies in the GeV and sometimes in the TeV range. Efficient particle acceleration converts the stellar rotational kinetic energy into radio, X-ray, and gamma-ray photons. Force-free magnetospheres, being dissipationless, cannot operate this conversion. Some non-ideal plasma effects must set in within the magnetosphere. In this paper, we compute numerical solutions of pulsar radiative magnetospheres in the radiation reaction limit, where radiation fully balances single particle acceleration. Using an appropriate Ohm's law, the dissipation is only controlled by the pair multiplicity factor $κ$. Moreover, we allow for either a minimal radiative region where dissipation is added only where required or for a force-free inside radiative outside model. This approach naturally and self-consistently connects the particle dynamics to its radiation field in the ultra-relativistic regime. Our solutions tend to the force-free limit for moderately large multiplicities, $κ ≫ 1$, decreasing the spin-down energy conversion into radiation. Nevertheless, for sufficiently low multiplicity $κ ≲ 1$, a significant fraction of the spin-down energy flows into radiation via particle acceleration. The work done by the electromagnetic field on the plasma mainly occurs in the current sheet of the striped wind, right outside the light-cylinder. Nevertheless the impact on the magnetic topology is negligible whatever the model. Therefore, the associated sky maps and light curves are only weakly impacted as shown.

**Key words:** magnetic fields – radiation mechanisms: non-thermal – methods: numerical – stars: neutron – pulsars: general – stars: rotation.

## 1 INTRODUCTION

Pulsars are very efficient particle accelerators as witnessed by the broad-band electromagnetic spectrum from radio (Manchester et al. 2005) up to very high-energy, GeV (Abdo et al. 2013), and sometimes TeV emission like for the Crab (Ansoldi et al. 2016) and Vela (Djannati-Ataï & for the H.E.S.S. Collaboration 2017). Particle acceleration and therefore radiation is rooted to the fast rotation of a strongly magnetized neutron star. Rotational kinetic energy is converted into radiation by curvature, synchrotron, and inverse Compton emission leading to the stellar braking accounted by the spin-down rate derived from the period *P* and its derivative $\dot{P}$. However, where exactly within the magnetosphere those mechanisms occur is still unclear. Undoubtedly particles flow at very high Lorentz factor from the star to the interstellar medium, shaping the pulsar wind as a ballerina similar to the solar wind. The global magnetosphere electrodynamics is intimately related to the motion of these particles and their subsequent radiation. Some localized regions are prone to efficient conversion of the rotational kinetic energy into acceleration and radiation but where and how remains to be self-consistently determined from global magnetosphere simulations including dissipation.

The simplest approach to find such solutions starts with force-free regime (FFE) where an ideal plasma is considered, neglecting particle inertia and temperature, meaning that the electric field *E*

is orthogonal to the magnetic field **B**, $E \cdot \mathbf{B} = 0$ and the electric field is weaker than the magnetic field in normalized units, meaning $E < c \mathbf{B}$ where *c* is the speed of light. In this picture, the Poynting flux is conserved because the electromagnetic field does no work on the plasma via the electric current **j**, meaning $\mathbf{j} \cdot \mathbf{E} = 0$. Strictly speaking, such magnetospheres are invisible because no photons are produced.

Numerical simulations have been pioneered by Contopoulos, Kazanas & Fendt (1999) for the axisymmetric rotator that was extended to an oblique rotator by Spitkovsky (2006) and retrieved by other authors, whether only aligned (Komissarov 2006; Timokhin 2006; Parfrey, Beloborodov & Hui 2012; Chen & Beloborodov 2014; Cao, Zhang & Sun 2016) or oblique (Kalapotharakos & Contopoulos 2009; Pétri 2012; Tchekhovskoy, Philippov & Spitkovsky 2016). See also different solutions not requiring a current sheet like for instance in Lovelace, Turner & Romanova (2006).

The aforementioned fluid description offers a good starting point to understand the global electric circuit made of charge and current densities. However, it neglects some fundamental kinetic aspects required to self-consistently include single particle acceleration as well as radiation feedback. As kinetic simulations are much more demanding than fluid models, this approach was only scarcely investigated in the last century. Let us mention Krause-Polstorff & Michel (1985) who computed axisymmetric dead pulsar magnetospheres called electrospheres. Due to the axisymmetry of the problem they used rings of charges instead of point particles. Later with the advent on computational power, Smith, Michel & Thacker (2001) showed with slightly more sophisticated simulations that a fully

<footnote>* E-mail: jerome.petri@astro.unistra.fr</footnote>
field magnetosphere is unstable and collapse to an electrosphere. The first full three-dimensional electrosphere was constructed by McDonald & Shearer (2009), using an electromagnetic Particle in Cell (PIC) code. They neglect pair creation and therefore did not add any particle injection process. Eventually Philippov & Spitkovsky (2014) computed the first two-dimensional axisymmetric pulsar magnetosphere for an aligned rotator by permanently injecting particle supposed to be released from the surface, avoiding to end to an electrosphere configuration (Pétri, Heyvaerts & Bonazzola 2002). Depending on the volume injection rate, they were able to find any equilibrium between the force-free and the fully charge separated state. Chen & Beloborodov (2014) improved this model by adding a prescription for the pair creation, putting a threshold on the lepton Lorentz factor. Following the same lines, Cerutti et al. (2015) assumed particle injection only from the vicinity of the stellar surface. Belyaev (2015) injected particles from regions where a parallel electric field exists. The first full three-dimensional PIC simulations of a pulsar magnetosphere were performed by Philippov, Spitkovsky & Cerutti (2015b). For an aligned rotator Philippov et al. (2015a) also included general-relativistic corrections with frame-dragging. Soon after some observational signature predictions were added to compute light curves and spectra emanating from curvature and or synchrotron radiation like for instance Cerutti, Philippov & Spitkovsky (2016b) who then included polarization (Cerutti, Mortier & Philippov 2016a). This PIC simulations were then extended to the striped wind well outside the light-cylinder to study its dissipation (Cerutti & Philippov 2017; Cerutti, Philippov & Dubus 2020). The oblique magnetosphere with radiation and general-relativistic correction was eventually computed by Philippov & Spitkovsky (2018). Several other groups performed similar simulations like Brambilla et al. (2018) or Kalapotharakos et al. (2017, 2018) who tried to explicitly connect their simulation results to gamma-ray observations. Alternatively, more simply test particle trajectories can be explored within a fluid code (see for instance Brambilla et al. 2015).

Even if the PIC approach is now mature to include several ingredients like pair creation and its subsequent radiative signature, its main flaw resides in its inability to simulate neutron star magnetospheres with realistic stellar magnetic field strengths and rotation periods. For instance, the Larmor radius is 10–15 orders of magnitude smaller than the light-cylinder radius, putting stringent constraints on the time step than cannot be fulfilled with current computational technology. It is therefore difficult to connect straightforwardly the microphysics dynamics induced by the gyro motion to the dynamics of the global magnetosphere, although the time and spatial scale hierarchy is maintained in current PIC simulations. Consequently, a fluid description, as the one we employ in this paper, remains a valuable tool to explore the dynamics of pulsar magnetospheres. A hybrid approach using a particle kinetic description wherever necessary and a fluid model elsewhere would represent a good compromise. Recent developments indeed combine the PIC technique to the MHD evolution as for instance performed by Marle et al. (2018) for particle shock acceleration or by Bai et al. (2015) for investigation of cosmic rays interaction with a thermal plasma. Such hybrid modules are also implemented in available MHD codes like PLUTO (Mignone et al. 2018). This new trend highlights the need to pursue our ongoing effort on improving plasma fluid models jointly with particle and Vlasov techniques by adding more physics on macro and micro scales simultaneously.

Returning to a simpler fluid and not particle description of the magnetosphere, the next step requires a proper treatment of dissipation with radiation. Some resistive simulations have been

performed but with a resistivity not always based on pure physical grounds or with some arbitrariness leading to no unique prescription (Li, Spitkovsky & Tchekhovskoy 2012; Kalapotharakos et al. 2012; Gruzinov 2008). Although the best view would be a full kinetic description including, acceleration and radiation is conceptually possible, we believe it better at this stage to use a fluid description of the plasma by deriving an Ohm's law according to the radiation reaction limit and derived in for instance Mestel (1999). Contopoulos (2016), build on this work and found radiative magnetospheric solution for oblique rotators however in a simplified manner. Recently, a fully radiative solution has been computed for aligned rotators by Pétri (2020b) and extended to oblique rotators Pétri (2020a), assuming that force-free dynamics holds inside the light-cylinder. Cao & Yang (2020) found similar solutions but allowing also possible dissipation within the light-cylinder. Recently, Cao & Yang (2022) computed high resolution radiative magnetosphere solutions including some test particle dynamics in order to predict synchrotron spectra and light curves.

In the radiation reaction limit, the equation of motion is solved for a single particle in a stationary regime, where the Lorentz force is counterbalanced by a so-called radiative friction for ultra-relativistic speeds. Assuming that the particle moves at exactly the speed of light leads to a unique solution for the velocity vector given some decades ago by Mestel (1999). This expression is sometimes also called Aristotelian electrodynamics. Electrons and positrons, although possessing the same electric drift motion because being independent of the particle charge, will move in opposite direction with respect to the perpendicular plane. Their charge density will generate a current leading to a one parameter family of current prescription reminded in the next section. Therefore, it should be clear that radiative simulations as those we present in this paper include single particle dynamics concretized through the electric current density prescription.

The radiation reaction approximation solves the single particle equation of motion in the ultra-relativistic regime, where the speed is exactly equal to the speed of light. The solution to the Lorentz force leads to the Aristotelian dynamics, the velocity being only a function of the local electric and magnetic field. However in order to avoid complication due to finite masses, as in the force-free case, we neglect their inertia. In this picture, particles only contribute to the charge and current density required to solve Maxwell equations. The derived current density possesses a component along the electric field and leads naturally to a kind of resistivity.

A full description of the plasma in this limit requires knowledge of its lepton content, separating the contribution from the electrons and the positrons. This pair multiplicity factor $\kappa$ serves this goal and is the only free parameter in this radiative Ohm law. A fully self-consistent picture must however add pair-production but this small-scale physics is still difficult to reconcile with the global scale of the magnetosphere. No self-consistent simulations have been performed so far taking all the ingredients self-consistently into account. Nevertheless, non-thermal acceleration and radiative feedback is included in this work thanks to the Aristotelian dynamics.

Because the radiation reaction limit regime relies on assumptions note always met within the magnetosphere, it is worth keeping in mind several caveats of our approach. The velocity field in Aristotelian electrodynamics is derived in the limit of significant radiative friction in the Lorentz force, reaching on a short time-scale an asymptotic regime of exact balance between electric field acceleration and radiation damping. While this regime could be achieved in many places within the magnetosphere, there exist
localized regions where such intense radiation damping is not effective due to negligible radiation reaction. Indeed, in the vicinity of polar caps, pair production efficiently screens the electric field component parallel to the magnetic field and particles do not experience the radiation reaction force (Timokhin & Arons 2013). Moreover, radiation damping involves ultra-relativistic particles with very large Lorentz factors that fails to be produced at several places except maybe for a tiny population of highly energetic particles. The outcome is a complex particle distribution function resembling more to a power law than to a mono-energetic population we take in this paper. As will be shown, Aristotelian dynamics allows for large regions where the electric field $E$ dominates the magnetic field $B$, i.e. where $E > cB$. However, recent studies showed that this should only occur as a transition stage to a magnetically dominated regime where $E < cB$ (Li, Beloborodov & Sironi 2021, see also Beskin 2018). Moreover, our treatment neglects magnetic reconnection, especially within the current sheet of the striped wind outside the light-cylinder, although it has been observed in PIC and MHD simulations. Our scheme represents a simplified two stage process where dissipation is directly converted into radiation.

In this paper, we compute oblique pulsar magnetospheres in the radiation reaction limit taking into account the exact dissipative current containing the electric drift component as well as the components aligned with the electric and magnetic field. In Section 2, we describe the model of our radiative magnetosphere and the prescription for Ohm's law derived from the radiation reaction regime. Some examples of magnetic topologies for an aligned and an orthogonal rotator are presented in Section 3 for the ideal FFE field and for the radiative magnetospheres. Next, in Section 4 we compute the spin-down luminosity extracted from these models and compare it with previous works. The importance of dissipation is pointed out in Section 5. The importance of the parallel electric field component is stressed in Section 6 influences on the polar cap shape and size is explored in Section 7. Sky maps and light curves are presented in Section 8. Conclusions are drawn in Section 9.

## 2 MAGNETOSPHERIC MODEL

In this section, we present the underlying model to compute radiative pulsar magnetospheres starting from Maxwell equations and then explaining the electric current prescription.

### 2.1 Maxwell equations

In our models, the plasma only furnishes the required charge $\rho_e$ and current $j$ densities to evolve Maxwell equations written in standard MKSA units as

$$ \nabla \cdot \mathbf{B} = 0, \qquad (1a) $$

$$ \nabla \times \mathbf{E} = - \frac{\partial \mathbf{B}}{\partial t}, \qquad (1b) $$

$$ \nabla \cdot \mathbf{E} = \frac{\rho_e}{\varepsilon_0}, \qquad (1c) $$

$$ \nabla \times \mathbf{B} = \mu_0 j + \frac{1}{c^2} \frac{\partial \mathbf{E}}{\partial t}. \qquad (1d) $$

Apart from the obvious boundary conditions on the stellar surface, the current density $j$ is the only unknown of the problem. Once fixed according to a given plasma model, Maxwell equations can be solved numerically, leading to a magnetosphere solution. So let us describe the possibilities for this current.

### 2.2 Current prescription

#### 2.2.1 Force-free limit

The simplest model corresponds to an ideal plasma with infinite conductivity, leading to the force-free prescription as

$$ j = \rho_e \frac{\mathbf{E} \wedge \mathbf{B}}{B^2} + \frac{\mathbf{B} \cdot \nabla \times \mathbf{B}/\mu_0 - \varepsilon_0 \mathbf{E} \cdot \nabla \times \mathbf{E}}{B^2} \mathbf{B}. \qquad (2) $$

By construction, this current does not work on particles since $j \cdot \mathbf{E} = 0$. All the rotational kinetic energy goes into the Poynting flux of the low-frequency, large-amplitude electromagnetic wave. For completeness and comparison with other models shown in this paper, we compute again some force-free magnetospheres.

#### 2.2.2 Radiative solution

If some emission is taken into account, for instance like in the radiation reaction limit of ultra-relativistic particles, the velocity of particles is fully determined by the local electromagnetic field configuration. Indeed, the friction caused by a radiative term can be seen as an isotropic emission of photons in the particle rest frame and at a rate controlled by the radiated power $\mathcal{P} \ge 0$ such that the balance between Lorentz force and radiative friction is

$$ q(\mathbf{E} + \mathbf{v} \wedge \mathbf{B}) = \frac{\mathcal{P}}{c^2} \mathbf{v}. \qquad (3) $$

A justification and argumentation about the validity of this balance can be found in Mestel et al. (1985), starting from the Lorentz–Abraham–Dirac equation. This equation is solved explicitly with respect to the velocity vector $\mathbf{v}$ and given for positive charges as $\mathbf{v}_+$ and negative charges as $\mathbf{v}_-$ according to

$$ \mathbf{v}_{\pm} = \frac{\mathbf{E} \wedge \mathbf{B} \pm (E_0 \mathbf{E}/c + c B_0 \mathbf{B})}{E_0^2/c^2 + B^2}. \qquad (4) $$

This expression only assumes that particles move exactly at the speed of light. $E_0$ and $B_0$ are the strength of the electric and magnetic field deduced from the electromagnetic invariants and satisfying $I_1 = E^2 - c^2 \mathbf{B}^2 = E_0^2 - c^2 B_0^2$ and $I_2 = c \mathbf{E} \cdot \mathbf{B} = c E_0 B_0$. Explicitly solving for $E_0 \ge 0$ and $B_0$, we find

$$ E_0^2 = \frac{1}{2} (I_1 + \sqrt{I_1^2 + 4 I_2^2}) \qquad (5a) $$

$$ c B_0 = \operatorname{sign}(I_2) \sqrt{E_0^2 - I_1}. \qquad (5b) $$

$E_0$ and $B_0$ are interpreted as the electric and magnetic field strength in a frame, where electric and magnetic field are parallel to each other.

The radiated power is then simply

$$ \mathcal{P} = |q| E_0 c \ge 0. \qquad (6) $$

Therefore, within physical constants, $E_0$ is a direct measure of the radiated power. Single particles are therefore evolved according to equation (4). It corresponds to the exact solution of the Lorentz equation of motion for charges subject to a friction and moving at exactly the speed of light. The particle inertia has been neglected and all species with the same sign of charge possess the same velocity vector irrespective of their charge to mass ratio $q/m$ as long as their sign does not change. From the velocity expression in (4), we can derive the electric current associated with this particle flow. The detailed derivation is given by Pétri (2016), and the associated radiative current density $j$ with minimal assumption is explained in Pétri (2020b). The final expression reduces to

$$ j = \rho_e \frac{\mathbf{E} \wedge \mathbf{B}}{E_0^2/c^2 + B^2} + (|\rho_e| + 2 \kappa n_0 e) \frac{E_0 \mathbf{E}/c + c B_0 \mathbf{B}}{E_0^2/c^2 + B^2}, \qquad (7) $$
where $κ$ is the pair multiplicity. The background particle density number is depicted by $n_0$ and varies from point to point within the magnetosphere. For a self-consistent picture, this density should be constrained by the pair production rate. However, such task is out of the scope of this study. In order to go further, we replace $n_0$ by $|ρ_e/q|$ and therefore is directly connected to the electric field via Maxwell–Gauss equation. Consequently,

$$|ρ_e| + 2κn_0e = |ρ_e|(1 + 2κ), \quad (8)$$

but actually this factor could be any function in the most general situation.

We stress that in our simulations, individual particles follow the velocity given in equation (4). Each particle possesses its own velocity depending solely on the local electromagnetic field configuration. Therefore, particles are present indirectly in the simulations with an analytical expression for the velocity, the Aristotelian dynamics. Nevertheless, we do not follow individual particle trajectory because this would also require some assumption about the particle injection rate and its spatial dependence. To avoid such arbitrariness, we fixed the local charge density to the Gauss–Maxwell expectations (1c).

We observe a difference between positive and negative charges because they move in opposite direction with respect to the electric and magnetic field direction. Radiation feedback is taken into account by a friction term in the Lorentz force, opposite to the velocity, which is proportional to the radiative power $P$, see equation (3). We do not use any electron-ion plasma in thermal equilibrium, rather an electron–positron plasma, although the difference between ions and positrons is anecdotal because the particle mass does not intervene in the ultra-relativistic regime, like photons. We also do not have to worry about the fluid motion because particle fill the whole space with a charge density given by Gauss law. This procedure is very similar to its force-free avatar.

In the radiative regime, dissipation of electromagnetic field is controlled by the electric field strength $E_0$ as measured in the frame where $\vec{E}$ and $\vec{B}$ are parallel because

$$j \cdot E = |ρ_e|(1 + 2κ)cE_0 \geq 0. \quad (9)$$

Expressed in terms of the radiated power, we find

$$j \cdot E = n_0 (1 + 2κ) P. \quad (10)$$

where $n_0 = |ρ_e/q|$ represents the particle density number required for the minimalistic model of a totally charge separated plasma.

We do not expect reconnection to play any role in the dissipation of the magnetic energy. The radiative Ohm’s law behaves as a resistive term. FFE is therefore broken when switching to the radiative solution. All the losses funnel into the particle velocity component along the magnetic field, making particles moving approximately at the speed of light, therefore copiously radiating energy and momentum. The only requirement is that kinetic energy losses being compensated by the electric field work. The balance equation (3) connects the radiative losses to the Lorentz force in a stationary state. In this picture we neglect particle inertia compared to the electromagnetic energy and radiative losses. To make an analogy with FFE, particle not only produce the required charge and current density but now they also produce some emission. Their dissipation rate is completely controlled by the radiative term and no reconnection is observed.

In the force-free limit, $E_0$ vanishes and dissipation disappears. In the minimalistic view, the current (7) is imposed only where necessary that is in regions where the condition $E < cB$ is violated whether inside or outside the light-cylinder. We call this model the radiative solution (RAD). We could also allow for less dissipation, for instance outside the light-cylinder, in the spirit of the force-free

inside/dissipative outside approach of Kalapotharakos et al. (2017). We call it force-free inside/radiative outside (FIRO).

### 2.2.3 Force-free inside/radiative outside

Inside the light-cylinder, the corotating electric field $E$ remains less than $cB$. It is therefore always possible to set force-free conditions in this region. However, outside the light-cylinder, the $E$ field can easily surpass the $B$ field strength. In such cases, we can artificially decrease the $E$ strength as in the force-free model. However, as a less stringent method and more realistically, we let the system evolve by adding dissipation not requiring any condition on the $E$ field as allowed by the radiative current prescribed previously in equation (7). Therefore in a last regime, we enforce a force-free inside radiative outside (FIRO) model, allowing force-free conditions inside and radiative dissipation outside the light-cylinder.

### 2.2.4 Numerical set-up

We performed several sets of runs with the aforementioned three regimes leading to a priori different magnetosphere models. The neutron star radius is set to $R/r_L = 0.3$ and the outer boundary of the simulation sphere is located at $7r_L$ where the light-cylinder is defined by $r_L = c/Ω$. This allows us to clearly compute the base of the striped wind on almost one wavelength. The pair multiplicity is set by the user, we chose $κ = \{0, 1, 2\}$.

The pair multiplicity $κ$ must always be a positive integer. It quantifies the deviation from a purely charge separated plasma. Indeed, in the minimalistic regime, a fully charge separated plasma requires $κ = 0$. If some weak pair production occurs within the plasma we chose low values such as $κ = 1, 2, 5$ or any small integer. For large multiplicities $κ \gg 1$, we will show that the solution tends quickly to the force-free magnetosphere when $κ$ augments. With $κ = 2$ the magnetosphere configurations becomes already indistinguishable from the FFE case. The multiplicity $κ$ is intimately related to the pair production efficiency within the magnetosphere. PIC simulations have shown that the pair injection process, rate, and location, crucially determines the outcome, tending either to a charge separated plasma forming an electrosphere or to an almost neutral plasma leading to a force-free and completely filled magnetosphere. The pair multiplicity remains so far largely unconstrained by observations. However, detailed numerical simulations of pair cascades around the polar caps performed by Timokhin (2010), Timokhin & Arons (2013) showed that values up to $κ = 10^5$ can be expected. Such high multiplicity leads to an almost perfect force-free regime. Nevertheless, the current prescription containing $κ$ as the only free parameter in our model could be supplemented by the freedom in the background charge density $ρ_0$, possibly differing from the standard corotating prescription given by $ρ_e \neq ρ_0$. This change in the background dynamics dramatically impacts the magnetosphere electrodynamics. Having no way to constrain this density $ρ_0$ we kept minimal assumption by imposing $ρ_0 = ρ_e$.

The pulsar obliquity is denoted by the angle $χ$. We implemented absorbing outer boundary conditions, meaning that the solution becomes unrealistic at distances $r > 5r_L$. In the following sections, we derive important quantities related to the pulsar electrodynamics such as its electromagnetic field structure, its spin-down losses, the work done on the plasma and the observational outcome relying on the polar cap shape, their light curves and the slot gap/striped wind emission properties. A numerical grid of $N_r \times N_θ \times N_φ = 257 \times 32 \times 64$ was sufficient for obtaining accurate solutions in all cases.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Convergence study for the spin-down luminosity $L$ normalized to the vacuum luminosity $L_{vac}$ of an orthogonal rotator depending on the grid resolution with $N_r$ radial points and $N_\theta$ latitudinal points.</figcaption></figure>

The adopted resolution stems from a convergence study of the spin-down luminosity. Acceptable accuracy is reached whenever the luminosity at the light-cylinder has converged within 1 per cent. We ran simulations for an orthogonal rotator with several radial and latitudinal grid points from the lowest resolution of 65 × 16 × 32 to the highest resolution of 257 × 128 × 256. The spin-down luminosity is plotted in Fig. 1 and shows that a resolution of 257 × 32 × 64 or even 129 × 32 × 64 is already sufficient for acceptable accuracy. This seemingly low resolution is actually due to the very low dissipation of our fully pseudo-spectral method compared to finite volume or finite difference methods. The grid resolution can be much coarser for spectral methods, especially if the solution is continuous.

Eventually, as a check of our algorithm, including filtering, de-aliasing, absorbing boundary layers and resolution, we computed vacuum solutions comparing our results with expectations from the Deutsch (1955) solution.

## 3 MAGNETIC FIELD LINES

Electrodynamics of neutron stars relies heavily on its electromagnetic field. We therefore start by showing the magnetic field structure. A full 3D picture being difficult to visualize on a sheet of paper, we restrict ourself to the geometry of magnetic field lines in the meridional plane for an aligned rotator and in the equatorial plane for an orthogonal rotator. Such lines are shown in Fig. 2 for the FFE limit and different radiative regimes for an aligned rotator on the left-hand panel and an orthogonal rotator on the right-hand panel. Because the pair multiplicity factor only weakly impacts on the geometry, we only show the cases with $\kappa = 0$. Compared to the cases $\kappa \in \{1, 2\}$, we have not found any significant changes, therefore they are not shown in Fig. 2. Inside the light-cylinder, shown as a black dashed line on the left-hand panel and as a circle on the right-hand panel, the magnetic field of the FFE and radiative cases are very similar, whatever the pair multiplicity. As expected the radiative solutions close more field lines along the equator outside the light-cylinder.

Next, we diagnose quantitatively the effect of a radiative magnetosphere by computing relevant physical quantities such as the spin-down luminosity and the work done on the plasma.

## 4 POYNTING FLUX

A radiative magnetosphere has the interesting property to allow conversion of the Poynting flux into particle acceleration and radiation accounting for the feedback of this current on to the electromagnetic field in a self-consistent manner. In this section, we report the efficiency of Poynting flux decrease depending on the model and on the pair multiplicity. The electromagnetic flux must be compared to the reference situation of a force-free magnetosphere for which no dissipation is expected by construction. However, such solutions develop current sheets that are strongest for an aligned rotator. Numerically such discontinuities are tricky to handle especially for a spectral method where the Gibbs phenomenon easily arises. Some artificial dissipation must be introduced to avoid strong oscillations.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Magnetic field lines for an aligned rotator, left-hand panel, and for an orthogonal rotator, right-hand panel, in the force-free limit (FFE) in blue, a minimalist radiative magnetosphere (RAD) in red and a FIRO magnetosphere in green, both latter with $\kappa = 0$.</figcaption></figure>
However for radiative models, the current prescription naturally leads to some dissipation controlled by a physical parameter like resistivity or radiation damping.

For purely electromagnetic interactions, energy is shared between three quantities, the electromagnetic energy density u defined by

$$u = \frac{\varepsilon_0 E^2}{2} + \frac{B^2}{2 \mu_0} \qquad (11)$$

the Poynting flux defined by

$$S = \frac{E \wedge B}{\mu_0} \qquad (12)$$

and the work done on the plasma represented by current density interacting with the electric field at a rate

$$\mathcal{D} = j \cdot E. \qquad (13)$$

Particle are assumed to have zero inertia in these simulations. Their velocity is governed by Aristotelian dynamics according to the local electromagnetic field equation (4). In this respect, the lost energy is directly converted into radiation because of this zero lepton mass limit. The strength of the radiative feedback depends on the $E_0$ field which is proportional to the radiated power as shown decades ago by Mestel et al. (1985). Magnetic energy is dissipated not via reconnection but via radiation damping, impacting the particle velocity and leading to equation (4).

The dissipative term $\mathcal{D}$ vanishes for a force-free plasma and in the radiation reaction limit it reduces to expression (9). The energy conservation law then reads

$$\frac{\partial u}{\partial t} + \nabla \cdot S + j \cdot E = 0. \qquad (14)$$

In a stationary state, the electromagnetic energy density u remains unchanged. Without dissipation, the Poynting flux across a closed surface is conserved but with for instance radiative losses energy flows into the plasma. From the conservation law equation (14) integrated within a sphere $\Sigma$ of radius $r$, we get

$$\iint_{\Sigma} S \cdot e_r d\Sigma = - \iiint_{V} j \cdot E dV, \qquad (15)$$

where $S \cdot e_r = S_r$ is the radial component of the Poynting flux, $d\Sigma$ a surface element on the sphere and $dV$ a volume element inside the sphere $\Sigma$.

The radial evolution of the Poynting flux is shown in Fig. 3 for several models, the force-free (FFE), the force-free inside/radiative outside prescription (FIRO) and the minimalistic radiative approach (RAD). As a check, the vacuum solution (VAC) is also shown to estimate the numerical dissipation. The luminosity is normalized with respect to the vacuum point dipole orthogonal rotator

$$L_{\text{vac}} = \frac{8\pi}{3\mu_0 c^3} \Omega^4 B^2 R^6 \qquad (16)$$

such that the flux plotted is $\ell = L/L_{\text{vac}}$.

The angular dependence on the obliquity $\chi$ is the same for all regimes (see Fig. 4). All the plasma filled fits are well approximated by

$$L/L_{\text{vac}} \approx 1.3 + 1.5 \sin^2 \chi. \qquad (17)$$

As a check, for the vacuum case we get

$$L/L_{\text{vac}} \approx 0.96 \sin^2 \chi. \qquad (18)$$

We observe an important dissipation of the Poynting flux for the FIRO case. In order to better localize this dissipative effect within

<figure><img src="image_6.png" /><figcaption>Figure 3. Radial decrease of the Poynting flux depending on the model. FFE is shown in blue, RAD in red and FIRO in green, both latter for $\kappa = 0$. For reference, the vacuum case is also shown in yellow.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. The Poynting flux crossing the light-cylinder for oblique rotators in vacuum (VAC) in yellow, force-free (FFE) in blue, RAD in red and FIRO in green. The solid lines show the best fit.</figcaption></figure>

the magnetosphere, we show in the next section the work done on the plasma for an aligned and an orthogonal case.

The sensitivity to the pair multiplicity $\kappa$ is only weakly perceptible because the prescription for particle injection according to the local electric field via Maxwell-Gauss law already tends to the FFE limit for low multiplicities. What effectively controls the spin-down losses and the magnetosphere solution, either closer to an electrosphere or to the FFE regime is the particle density number. PIC simulations have also shown that the injection procedure is critical for the final outcome of the simulation.

Fig. 5 shows the different spin-down luminosities at the light-cylinder depending on the plasma model. The variation in luminosities remains very small except although the radiative solution seems slight more dissipative. The sensitivity to the pair multiplicity is also only weakly perceptible. Most of the dissipation occurs outside the light-cylinder as shown in the next section.

The spin-down luminosity is relatively insensitive to the dissipation term because the particle injection scheme follows the force-free scheme by dropping particles where the Maxwell-Gauss law imposes it. Therefore, the radiative model tends quickly to a nearly FFE state, even with a low to moderate multiplicity factor $\kappa$.
<figure><img src="image_8.png" /><figcaption>Figure 5. Comparison of spin-down luminosities at the light-cylinder depending on the plasma model. The integer within the bracket denotes the pair multiplicity factor (κ).</figcaption></figure>

## 5 DISSIPATION

Conservation of the total energy implies that some electromagnetic energy went into particle acceleration and radiation. The radial decrease in the Poynting flux L indicates a sink of electromagnetic energy imputed to the presence of a non-ideal plasma. The location where this conversion arises is important for the prediction of observational signatures such as radio and gamma-ray light curves and spectra. Within a spherical shell of radius r, this dissipation is given by the opposite of the Poynting flux radial derivative as

$$W_E = \iint_{\Sigma} j \cdot E \, d\Sigma = -\frac{dL}{dr}. \qquad (19)$$

Fig. 6 shows how fast dissipation occurs outside the light-cylinder depending on the radius. For FFE magnetospheres, dissipation vanishes but because of our numerical filtering procedure and of grid size effects, a small residual work is observed. In the FIRO model, dissipation starts at the light-cylinder, increasing to a maximum in the interval [1, 2] $r_L$ and then decreases slowly. In the RAD regime, the dissipation sets in with a delay, increasing significantly only beyond a radius $r \gtrsim 2 r_L$. This delayed dissipation impacts on the radio time

lag of gamma-rays photons with respect to radio photons. We expect to observe a decrease in this time lag for radio loud gamma-ray pulsar, helping to better jointly fit radio and gamma-ray light curves, as done in Pétri & Mitra (2021).

In order to better localize the radiative regions where the particle dynamics is the most important, we show maps of the work done locally on the plasma by computing the power defined in equation (9) for aligned and orthogonal rotators. Enlightening cases are shown for FIRO and RAD models in Fig. 7 for an aligned rotator and in Fig. 8 for an orthogonal rotator. The Poynting flux flows into the plasma outside the light-cylinder in the vicinity of the current sheet of the striped wind. With increasing distance, the power sharply decreases by two orders of magnitude at the outer boundary because of the decreasing electric field and current density. The thickness of this dissipative region is about $0.2 r_L$.

This dissipation layers are the privileged places, where high-energy radiation is produced and detected as pulsed gamma-ray emission. This fact is supported by the investigation of radio loud young gamma-ray pulsars for which the rotating vector model is consistent with gamma-ray light curve fitted for more than a dozen of pulsars (Pétri & Mitra 2021).

The acceleration and radiation processes are implicitly implemented by the Aristotelian velocity dynamics. Particles move at the speed of light because of zero inertia approximation and the radiated power is controlled by the $E_0$ field. Particles in the simulations are present but they only contribute to the charge and current density as in the FFE approximation.

In the closed field zone, within the light-cylinder, for the FIRO model, we found a dissipation rate $D$ which is less than $10^{-4}$ or even $10^{-5}$. Such small values are almost zero from a numerical point of view. The plasma really remains force-free as it should within the numerical error of the algorithm.

The connection between dissipation layers and radiation zones has been confirmed by kinetic simulations such as Philippov & Spitkovsky (2018) and Chen, Cruz & Spitkovsky (2020).

The current sheet thickness is governed by the local physics, the breakdown of the FFE conditions. It is controlled by the radiative term and not by numerical dissipation which has been checked to remain negligible compared to the dissipation introduced by the radiative Ohm law.

<figure><img src="image_9.png" /><figcaption>Figure 6. Efficiency of the dissipation according to the integral in equation (19) for the FIRO model on the left-hand panel and for the RAD model on the right-hand panel. The pair multiplicity is shown in the legend with different colours. For each multiplicity $κ$, the dissipation is shown for all obliquities in the same colour.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Work done on the plasma for $\kappa = 0$ as given by equation (9), for an aligned rotator with the FIRO model on the left-hand panel and the RAD model on the right-hand panel.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Work done on the plasma for $\kappa = 0$ as given by equation (9), for an orthogonal rotator with the FIRO model on the left-hand panel and the RAD model on the right-hand panel.</figcaption></figure>

## 6 PARALLEL ELECTRIC FIELD

In order to quantify the presence of a parallel electric field $E_\|$, we plot some maps of the strength of $E_\| = E \cdot B/B = E_0 B_0/B$ in the observer frame. This component of the electric field could be responsible for particle acceleration and therefore represents a good indicator of the deviation from force-free conditions. Actually, significant values of $E_\|$ are coincident with the dissipation maps shown in the previous section, in Figs 7 and 8.

Fig. 9 shows two $E_\|$ maps for an aligned rotator with $\kappa = 0$, for the FIRO model on the left-hand panel and for the RAD model on the right-hand panel. Fig. 10 shows the same quantities for an orthogonal rotator. Particles do not follow field lines any more due to dissipation. They are efficiently accelerated by the electric field *E* along the magnetic field *B* in regions where dissipation is maximal.

Particles are evolved implicitly according to equation (4). They follow trajectories imposed by the local electromagnetic field (*E*, *B*)

and possess a component along *E*. This electric field aligned acceleration captures non-ideal effects deviating from the pure force-free picture. *E* can be larger than *c* *B* and there exist a parallel electric field component $E_\|$ responsible for these non-ideal effects.

There are no strong FFE violation because the particle injection scheme resembles very much to the force-free scheme and therefore the radiative model is still able to tend to a nearly FFE state. The spin-down is not much affected but the electric field and the corresponding dissipation rate term equation (13) are very sensitive to the radiative mechanism however only in very localized areas where $E > cB$.

## 7 POLAR CAPS

As a preparation for the investigation of the radio and gamma-ray light curves of radio loud gamma-ray pulsars detected by Fermi/LAT and related to our radiative magnetosphere, we compute the shape of
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Absolute value of the parallel electric field $E_\|$ for an aligned rotator, right-hand panel, in the minimalist radiative magnetosphere (RAD) in red and a FIRO magnetosphere in green, both with $\kappa = 0$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Same as Fig. 9 but for an orthogonal rotator.</figcaption></figure>

the polar cap in the different plasma regimes, comparing them to the force-free limit. For better readability with different obliquities $\chi$, the origin of the plots corresponds to the location of the magnetic north pole, the axes are defined locally for each obliquity by performing a rotation from the rotation axis to the magnetic axis.

Illustrative examples are shown in Fig. 11 for the polar cap rim with $\chi = \{15^\circ, 45^\circ, 75^\circ\}$, in vacuum (VAC), force-free (FFE), FIRO, and RAD regimes. The vacuum polar cap shapes computed from our pseudo-spectral simulations are shown in orange solid line and checked against those polar caps found from the exact analytical Deutsch solution and shown in dashed blue lines. The agreement between both contours is excellent and gives us confidence about our results for plasma fields magnetospheres. The polar caps for the FFE, FIRO, and RAD regimes are also shown in Fig. 11, respectively, in blue, green and red solid line. We have not noticed any significant change in these caps and their rims almost overlap whatever the regime. Nevertheless, compared to vacuum, the area of these polar

caps is larger than in the vacuum case because magnetic field lines open up due to the magnetospheric current. The presence of the plasma inflates the caps.

We conclude that the impact of radiation on the polar cap shape remains rather weak when compared to the force-free solution. Even if the ideal plasma approximation does not produce any parallel electric fields that could accelerate particles and therefore does not produce any radiation, from a geometrical point of view, it nevertheless offers a reasonably faithful picture of the electromagnetic field of a dissipative magnetosphere, wherever the radiative dissipation occurs, everywhere outside the light-cylinder or only where required by local conditions imposed by the electromagnetic field. We stress that within a fluid model, as the one employed in our study, the acceleration mechanism cannot produce power law distribution functions. It only heats up particles by keeping them thermal. This acceleration is accounted for by the parallel electric fields $E_\|$ component for a resistive plasma for which magnetic
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Polar cap geometry for a rotator with obliquity $\chi = \{15^\circ, 45^\circ, 75^\circ\}$ in vacuum (VAC), force-free (FFE) or in the radiation reaction limit, respectively, FIRO and RAD. The blue dashed line corresponds to the polar cap computed from the exact Deutsch solution.</figcaption></figure>

energy is dissipated into particle kinetic energy. This requires an Ohm's law like equation (7) that allows the presence of a parallel electric field $E_\|$. Injecting test particles into this configuration would offer a good compromise between a fully kinetic simulation and a fluid approximation, leading to possible non-thermal acceleration and building power-law distribution functions for those particles.

The observation that force-free and dissipative magnetospheres look rather similar relies mainly on the fact that the particle injection scheme derives from Maxwell-Gauss law, sowing particles locally at the rate imposed by the electric field. In such an approach, a fully charge separated plasma present in the whole magnetosphere deviates only slightly from the force-free counterpart. This conclusion is not an artefact from the grid resolution that could impact on the physical dissipative term but a natural consequence of the lepton injection method. Replacing the charge density derived from Maxwell-Gauss by another explicit spatiotemporal injection dependence would lead to more effective dissipation. But this is at the expense of adding more arbitrariness into the model that we wanted to avoid in the present investigation.

Small changes on the stellar surface are amplified at the light-cylinder and beyond; therefore, we could expect a significant change in the multiwavelength light-curve predictions. This last point connecting our simulations to observations is touched in the next section through comparison of sky maps in the radio and the gamma-ray band.

## 8 PULSED EMISSION

Dissipative magnetospheres are necessary to produce some radiation as detected by a distant observer. The location and geometry of the emission regions strongly imprint on the multiwavelength light curves and the phase-resolved spectra. In this last section, we compute sky maps and light curves for the aforementioned models, highlighting the differences expected depending on the plasma regime, ideal, or radiative.

We consider the three main emission regions to be, first the polar cap for radio photons, second the slot gap for high-energy gamma-ray photons and third the striped wind model for high and very high-energy gamma-rays up to the TeV range. For radio photons, we assume a polar cap model with emissivity shaped by a Gaussian function centred on the magnetic axis. For gamma-ray photons, we assume a striped wind model with emissivity starting at the light-

cylinder and focused along the current sheet or a slot gap extending from the surface (meaning here $0.3 r_\perp$) up to the light-cylinder. Some more details about these emission models can be found in Pétri (2018).

In our model, the size of the radio cone emission is controlled by the polar cap rims computed in the previous section. We assume that photon are produced with an altitude in the interval $[0.3, 0.4] r_\perp$. In this region close to the stellar surface, the electric field remains weak compared to the magnetic field and particles follow almost magnetic field lines. Photons are therefore shoot in a direction tangential to the magnetic field lines, including retardation and aberration effects. Moreover, we implemented a Gaussian profile centred along the magnetic moment axis with a typical width equal to the cone supported by the last open field lines. We adopted this picture in order to mimic the true radio profiles observed in many pulsars.

Gamma-ray photons are produced in the equatorial current within the striped wind, outside, the light-cylinder. Because the wind is expanding almost radially, these photons are emitted in the radial direction. The emissivity is maximal at the centre of the current sheet and decreases following another Gaussian shape when deviating from this sheet. The emission zone extend from $1r_\perp$ to $3r_\perp$. Much more details can be found in Benli, Pétri & Mitra (2021) and Pétri & Mitra (2021) and references therein.

In order to probe the radio emission mechanism from the polar cap and the gamma-ray emission from higher altitudes, kinetic physics is required to capture the non-ideal electric field, the gap formation, and the pair production. However, such study is out of the present scope. Attempts to better understanding the radio emission generation have been pursued by for instance Philippov, Timokhin & Spitkovsky (2020) and for the gamma-ray emission by for instance Kalapotharakos et al. (2017). Here, we are only interested in the radio and high energy pulse profiles implied by the geometrical configuration of radiative magnetospheres. We occult the detailed energetics of individual particles, radiation, and electromagnetic interactions.

### 8.1 Sky maps

We start by reckoning a full set of light curves in radio and high-energy, following the three emission regions. The combined polar cap/striped sky maps are summarized in Fig. 12 for $\chi = \{15^\circ, 45^\circ, 75^\circ\}$ and the combined polar cap/slot gap is summarized in
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Sky maps of the ideal and radiative models compared to the vacuum case. From left to right column: Vacuum (VAC), force-free (FFE), FFE inside/radiative outside (FIRO), and radiative (RAD) on the right. Each line represents a different obliquity given from top to bottom by $\chi = \{15^\circ, 45^\circ, 75^\circ\}$. The gamma-ray emission assumes a striped wind model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Same as Fig. 12 but for the gamma-ray emission assumes a slot gap model.</figcaption></figure>

Fig. 13 for the same obliquities. It is also instructive to show the expectations from the Deutsch vacuum solution. The observer line-of-sight inclination $\zeta$ varies from 0° to 180°. Each column in this plot depicts a particular model. From left-to-right column, we have successively the vacuum (VAC), the force-free (FFE), the FFE inside/radiative outside (FIRO), and the radiative (RAD) regimes. Each line represents a different obliquity given from top to bottom by $\chi = \{15^\circ, 45^\circ, 75^\circ\}$.

Interestingly, the force-free maps are very similar to the radiative maps. The FFE case and radiative runs are different in their details. Radiative models clearly show some regions where the electric field exceeds the magnetic field $E > cB$ and some dissipative regions, where $D$ is non-negligible. Whereas the global magnetosphere is not drastically impacted by the local dissipative terms, the sky maps are.

Even the vacuum solution already produces maps resembling the plasma filled cases, especially when looking close to the equator ($\zeta \approx 90^\circ$). Inspecting more carefully Figs 12 and 13, we conclude that the

actual plasma regime only weakly impacts on the sky maps, whether in radio or in gamma-rays. From a geometrical point of view, when investigating light-curve shapes, the force-free limit allows for an accurate study of pulse profile without adding any free parameter into the game. The geometric dependence on $\chi$ and $\zeta$ is already faithfully reproduced in FFE. However, when energetic considerations come into play, the radiative models will generate very different phase-resolved spectra and multiwavelength light curves because of the varying parallel electric field acting on the particle dynamics. This requires deeper investigation of particle acceleration and radiation that we leave for future work.

## 8.2 Light curves

As a typical example of different light curves constructed from these models, we plot an atlas of gamma-ray light curves in Fig. 14 for $\chi = \{15^\circ, 45^\circ, 75^\circ\}$ and in steps of 10° for $\zeta \in [0^\circ, 90^\circ]$, according to the striped wind and the slot gap model in the force-free limit.

For the plasma filled models, we only notice a variation of several percent in phase lag between the light curves. To a large extent, the double peak gamma-ray separation remains insensitive to the model used. Discrepancies in phase are difficult to detect in real gamma-ray observations. However, more importantly are the variations in the peak maximal intensity between the models, notably the reversal of the dominant peak, leading or trailing with respect to radio, when switching from FFE/FIRO to RAD model (see Fig. 15) with $\{\chi, \zeta\} = \{45^\circ, 40^\circ\}$ on the left-hand panel and $\{\chi, \zeta\} = \{75^\circ, 60^\circ\}$ on the right-hand panel. The slot gap model produces much wider profiles, broader than what is observed by Fermi/LAT such that we favour the striped wind to explain GeV light curves although for some pulsar showing a kind of plateau emission the slot gap would better fit the observations. As a conclusion, we think that the majority of the gamma-ray pulsars fall into the striped wind emission model except for some outliers possibly also dominantly emitting in the slot gap sites.

## 9 CONCLUSIONS

Finding self-consistent dissipative pulsar magnetospheres is important to localize possible sites for their broadband pulsed emission. The underlying processes of particle acceleration and radiation remain to be properly identified. In this paper, we constructed self-consistent radiative pulsar magnetospheres in the radiation reaction regime, where ultra-relativistic particles flow around the neutron star. We introduced partially and minimalistic radiative models, demonstrating that the global electromagnetic topology remains mainly insensitive to the actual dissipation regime. In any case, most of the radiation occurs in the current sheet outside the light-cylinder but close to it.

Although different models produce marginal differences in the gamma-rays and radio light curves, we believed that current observations especially in high energy are not sensitive enough to disentangle for instance the minimalistic from the partially radiative magnetosphere. For better comparison with observations, the computation of phase-resolved spectra and multiwavelength light curves will undoubtedly help to segregate between competing dissipative models like the radiative or resistive magnetospheres introduced in the literature. The energetic of the magnetosphere model will leave the degeneracy contrary to a pure geometrical study of light curves. As phase-resolved spectra are available in radio and gamma-rays, we plan to compute multiwavelength light curves and spectra based on the above magnetosphere models.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Atlas of striped wind and slot gap light curves for $\chi = \{15^\circ, 45^\circ, 75^\circ\}$ and $\zeta \in [0^\circ, 90^\circ]$ in steps of $10^\circ$, see the inset in the format $\{\chi, \zeta\}$. The first, third, and fifth columns are for the striped wind, whereas the second, fourth, and sixth columns are for the slot gap.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Examples of gamma-ray light curves extracted from the sky maps shown in Figs 12 and 13 with the geometry given in the inset. On the top panel the results for the striped wind and on the bottom panel for the slot gap.</figcaption></figure>

The physics of radiative magnetospheres requires a more detailed investigation of the central role of particle injection and its crucial impacts on the magnetosphere energetics before to confront to the observations. Some more ingredients are required before investigating the gamma-ray light curves of individual pulsars.

A large amount of work has been carried out to identify these effects via in particle-in-cell simulations, taking the feedback between particle acceleration and radiation self-consistently into account. Relativistic magnetic reconnection is supposed to play a key role in the dissipation of the Poynting flux channelling into particle acceleration and emitting synchrotron photons as demonstrated by Cerutti et al. (2016b) and Philippov & Spitkovsky (2018). The particle production assumption, either from the surface or from the light-cylinder or from the whole magnetosphere volume decides what on the global solution found in numerical simulations. Chen et al. (2020) presented simulations not requiring efficient pair production in the vicinity of the light-cylinder and found quasi-periodic solution able to explain the cone versus and core emission characteristics of
radio pulses. Previously Chen & Beloborodov (2014) already found
significant curvature and synchrotron photon production around
the current sheet of an aligned rotator. These results are able to
reproduce the Fermi gamma-ray pulsar observations. In the same
vain, Kalapotharakos et al. (2017) computed test particle trajectories
and the associated radiation and particle efficiency in the so called
FIDO model assuming curvature emission in the radiation reaction,
in a similar way to the present study. Based on gamma-ray pulsar
luminosity, they found a positive correlation between the pair
multiplicity factor and the pulsar spin-down.

Nevertheless almost all these simulations rely on particle injection
prescriptions which are still not fully resolved from a physical point
of view. As the magnetospheric solution tend to depends crucially
on the injection rate, it is not yet clear how all these processes
operate. Moreover, the Lorentz factors obtained are still many orders
of magnitude below realistic values expected from observations. Our
approach represents a good alternative to tackle this important issue
of very high Lorentz factor within the magnetosphere.

ACKNOWLEDGEMENTS

I am grateful to the referee for helpful comments and suggestions.
This work has been supported by the CEFIPRA grant IFC/F5904-
B/2018 and ANR-20-CE31-0010. We acknowledge the High Per-
formance Computing centre of the University of Strasbourg for
supporting this work by providing scientific support and access to
computing resources.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request
to the corresponding author.

REFERENCES

Abdo A. A. et al., 2013, ApJS, 208, 17

Ansoldi S. et al., 2016, A&A, 585, A133

Bai X.-N., Caprioli D., Sironi L., Spitkovsky A., 2015, ApJ, 809, 55

Belyaev M. A., 2015, MNRAS, 449, 2759

Benli O., Pétri J., Mitra D., 2021, A&A, 647, A101,

Beskin V. S., 2018, Phys.-Usp., 61, 353,

Brambilla G., Kalapotharakos C., Harding A. K., Kazanas D., 2015, ApJ,
804, 84

Brambilla G., Kalapotharakos C., Timokhin A. N., Harding A. K., Kazanas
D., 2018, ApJ, 858, 81

Cao G., Yang X., 2020, ApJ, 889, 29

Cao G., Yang X., 2022, ApJ, 925, 130,

Cao G., Zhang L., Sun S., 2016, MNRAS, 455, 4267

Cerutti B., Philippov A. A., 2017, A&A, 607, A134

Cerutti B., Philippov A., Parfrey K., Spitkovsky A., 2015, MNRAS, 448, 606

Cerutti B., Mortier J., Philippov A. A., 2016a, MNRAS, 463, L89

Cerutti B., Philippov A. A., Spitkovsky A., 2016b, MNRAS, 457, 2401

Cerutti B., Philippov A. A., Dubus G., 2020, A&A, 642, A204

Chen A. Y., Beloborodov A. M., 2014, ApJ, 795, L22

Chen A. Y., Cruz F., Spitkovsky A., 2020, ApJ, 889, 69,

Contopoulos I., 2016, MNRAS, 463, L94

Contopoulos I., Kazanas D., Fendt C., 1999, ApJ, 511, 351

Deutsch A. J., 1955, Ann. d'Astrophys., 18, 1

Djannati-Ataï A. for the H. E. S. S. Collaboration, 2017, H.E.S.S. stereoscopic
observations of the Vela Pulsar above 100 GeV. Texas Symposium, Cape
Town, South Africa

Gruzinov A., 2008, preprint (arXiv:0802.1716)

Kalapotharakos C., Contopoulos I., 2009, A&A, 496, 495

Kalapotharakos C., Kazanas D., Harding A., Contopoulos I., 2012, ApJ, 749, 2

Kalapotharakos C., Harding A. K., Kazanas D., Brambilla G., 2017, ApJ,
842, 80

Kalapotharakos C., Brambilla G., Timokhin A., Harding A. K., Kazanas D.,
2018, ApJ, 857, 44

Komissarov S. S., 2006, MNRAS, 367, 19

Krause-Polstorff J., Michel F. C., 1985, MNRAS, 213, 43P

Li J., Spitkovsky A., Tchekhovskoy A., 2012, ApJ, 746, 60

Li X., Beloborodov A. M., Sironi L., 2021, ApJ, 915, 101

Lovelace R. V. E., Turner L., Romanova M. M., 2006, ApJ, 652, 1494

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Marle V., Jan A., Casse F., Marcowith A., 2018, MNRAS, 473, 3394

McDonald J., Shearer A., 2009, ApJ, 690, 13

Mestel L., 1999, Stellar Magnetism. Oxford, Clarendon

Mestel L., Robertson J. A., Wang Y. M., Westfold K. C., 1985, MNRAS, 217, 443,

Mignone A., Bodo G., Vaidya B., Mattia G., 2018, ApJ, 859, 13

Parfrey K., Beloborodov A. M., Hui L., 2012, MNRAS, 423, 1416

Pétri J., 2012, MNRAS, 424, 605

Pétri J., 2016, J. Plasma Phys., 82, 635820502

Pétri J., 2018, MNRAS, 477, 1035

Pétri J., 2020a, Universe, 6, 15,

Pétri J., 2020b, MNRAS, 491, L46

Pétri J., Mitra D., 2021, A&A, 654, A106

Pétri J., Heyvaerts J., Bonazzola S., 2002, A&A, 384, 414

Philippov A., Timokhin A., Spitkovsky A., 2020, Phys. Rev. Lett., 124, 245101,

Philippov A. A., Spitkovsky A., 2014, ApJ, 785, L33

Philippov A. A., Spitkovsky A., 2018, ApJ, 855, 94

Philippov A. A., Cerutti B., Tchekhovskoy A., Spitkovsky A., 2015a, ApJ,
815, L19

Philippov A. A., Spitkovsky A., Cerutti B., 2015b, ApJ, 801, L19

Smith I. A., Michel F. C., Thacker P. D., 2001, MNRAS, 322, 209

Spitkovsky A., 2006, ApJ, 648, L51

Tchekhovskoy A., Philippov A., Spitkovsky A., 2016, MNRAS, 457, 3384

Timokhin A. N., 2006, MNRAS, 368, 1055

Timokhin A. N., 2010, MNRAS, 408, 2092

Timokhin A. N., Arons J., 2013, MNRAS, 429, 20

This paper has been typeset from a *TEX/LATEX* file prepared by the author.