Pulsar gamma-ray emission in the radiation reaction regime

TL;DR

This paper models pulsar gamma-ray emission by computing particle trajectories in strong electromagnetic fields, reproducing observed spectral features and light-curve morphologies with minimal assumptions.

Contribution

It introduces a novel approach to simulate pulsar gamma-ray emission using radiation reaction limit trajectories in vacuum dipole fields, matching observations.

Findings

Spectral break at several GeV consistent with Fermi data

Emission sites localized within the pulsar magnetosphere

Single or double-peaked light-curves depending on geometry

Abstract

Since the era of the Fermi/LAT and atmospheric Cerenkov telescopes, pulsars are known to emit high and very high-energy photons, in the MeV-GeV range and sometimes up to TeV. To date, it is still unclear where and how these photons are produced. Nevertheless gamma-ray photons require particle acceleration to ultra-relativistic speeds. In this paper, we compute single particle trajectories for leptons in an arbitrary strong electromagnetic field in the so-called radiation reaction limit. In this picture, particle velocity only depends on the local electromagnetic field which we assume to follow the vacuum dipole rotator. From this velocity field, we compute the curvature radiation spectrum and light-curves. Sky maps and phase-resolved spectra are then deduced accounting for realistic pulsar periods and magnetic field strengths. Emission sites within the pulsar magnetosphere where most of radiation emanates are then localized. For standard parameters of millisecond and normal pulsars, we show that a break in the spectrum occurs at several GeV in agreement with the Fermi/LAT second pulsar catalogue. A sample of representative phase-resolved spectra and sky-maps are shown. A pair multiplicity of several tenths to several thousands is required to account for the total gamma-ray luminosity. Moreover depending on the geometry, single or double-peaked light-curves are found. Our model shows that minimalist assumptions are already able to reproduce salient features of pulsar emission.

Equations

$$\frac{{\omega }_B}{\Omega }=\frac{eB}{m_{\mathrm{e}}\Omega }=\text{numprint}2.8e18\left(\frac{P}{1\text{SIunits}\mathrm{s}}\right)\left(\frac{B}{\text{numprint}e8\text{SIunits}\mathrm{T}}\right).$$

$$r_{\mathrm{B}}=\frac{\gamma mc}{eB}=\text{numprint}1.7e-5\text{SIunits}\mathrm{m}\left(\frac{\gamma }{\text{numprint}e6}\right){\left(\frac{B}{\text{numprint}e8\text{SIunits}\mathrm{T}}\right)}^{-1}$$

$$\frac{{\omega }_B}{\Omega }=\text{numprint}2.8e9.$$

$$\text{mn@boldsymbol}{v}_{\pm }=\frac{\text{mn@boldsymbol}E\wedge \text{mn@boldsymbol}B\pm (E_0\text{mn@boldsymbol}E/c+c{B}_0\text{mn@boldsymbol}B)}{E_0^2/c^2+B^2}$$

$$\text{mn@boldsymbol}{E}^2-c^2\text{mn@boldsymbol}{B}^2$$

$$\text{mn@boldsymbol}E\cdot \text{mn@boldsymbol}B$$

$$\text{mn@boldsymbol}{v}_{\pm }=\pm c\frac{(\text{mn@boldsymbol}E\cdot \text{mn@boldsymbol}B)\text{mn@boldsymbol}B}{E_0(E_0^2/c^2+B^2)}.$$

$$\text{mn@boldsymbol}{v}_{\pm }=\pm c\text{sign}(B_0)\frac{\text{mn@boldsymbol}B}{B}$$

$$\text{mn@boldsymbol}{a}_{\pm }=\frac{d\text{mn@boldsymbol}{v}_{\pm }}{dt}=\frac{c^2}{{\rho }_c}.$$

$$\frac{dI}{d\omega }=\frac{\sqrt{3}{e}^2}{4\pi {\epsilon }_{0}c}\gamma \frac{\omega }{{\omega }_c}{\int }_{\omega /{\omega }_c}^{+\infty }{K}_{5/3}(x)dx$$

$${\omega }_c=\frac{3}{2}{\gamma }^3\frac{c}{{\rho }_c}$$

$$P_c=\frac{e^2}{6\pi {\epsilon }_0}{\gamma }^4\frac{c}{{\rho }_c^2}.$$

$$j_{\mathrm{cur}}(\mathbf{r},\omega )=\frac{\sqrt{3}}{2\pi }{\alpha }_{\mathrm{sf}}\frac{\hslash c}{{\rho }_c(\mathbf{r})}\gamma F\left(\frac{\omega }{{\omega }_c(\mathbf{r})}\right)$$

$${\alpha }_{\mathrm{sf}}=\frac{e^2}{4\pi {\epsilon }_0\hslash c}$$

$$F(x)=x{\int }_x^{+\infty }{K}_{5/3}(t)dt.$$

$${\gamma }^4=\frac{6\pi {\epsilon }_0}{e}{E}_0{\rho }_c^2.$$

$$\ell =\frac{\gamma m_{\mathrm{e}}{c}^2}{e{E}_0}.$$

$$\frac{\ell }{{\rho }_c}$$

$$=\text{numprint}4.7e-6{\left(\frac{E_0}{\text{numprint}e12\text{SIunits}V/m}\right)}^{-3/4}{\left(\frac{{\rho }_c}{12\text{SIunits}\mathrm{km}}\right)}^{-1/2}.$$

$$\frac{d{I}_{\mathrm{tot}}}{d\omega dt}={\iiint }_V{j}_{\mathrm{cur}}(\mathbf{r},\omega )n(\mathbf{r})d^3\mathbf{r}$$

$$n(\mathbf{r})=n_0{\left(\frac{R}{r}\right)}^q$$

$${\text{lambdabar}}_c=\frac{\hslash }{m_{\mathrm{e}}c}.$$

$${\chi }_{\mathrm{sync}}=\frac{3}{2}\gamma \frac{B}{B_{\mathrm{qed}}}\approx 1.$$

$${\chi }_{\mathrm{curv}}=\frac{3}{2}{\gamma }^2\frac{{\text{lambdabar}}_c}{{\rho }_c}$$

$$E_{\mathrm{Schw}}=\frac{m_{\mathrm{e}}^2{c}^3}{e\hslash }=10^{18}\text{SIunits}V/m$$

$$k_c=\frac{\hslash {\omega }_c}{m_{\mathrm{e}}{c}^2}=\frac{3}{2}{\gamma }^3\frac{{\text{lambdabar}}_c}{{\rho }_c}=\gamma {\chi }_{\mathrm{curv}}.$$

$${\gamma }^4=\frac{3}{2}\frac{\tilde{E}}{{\alpha }_{\mathrm{sf}}}\frac{{\rho }_c^2}{{\text{lambdabar}}_c^2}.$$

$$P_c=\frac{m_{\mathrm{e}}^2{c}^4}{\hslash }\tilde{E}\approx \text{numprint}6.36e7\text{SIunits}\mathrm{W}\tilde{E}.$$

$$\tilde{E}=\tilde{B}\frac{R}{r_{\mathrm{L}}}.$$

$$n_0=\frac{2{\epsilon }_0\Omega B}{e}=\frac{1}{2\pi {\alpha }_{\mathrm{sf}}{\text{lambdabar}}_c^2}\frac{\tilde{B}}{r_{\mathrm{L}}}\approx \text{numprint}1.46e21\text{SIunits}{\mathrm{m}}^{-3}{\left(\frac{r_{\mathrm{L}}}{10^5\text{SIunits}\mathrm{m}}\right)}^{-1}\tilde{B}.$$

$$\frac{d\tilde{I}}{d\tilde{\omega }d\tilde{t}}=\frac{\sqrt{3}}{2\pi }{\alpha }_{\mathrm{sf}}\frac{{\text{lambdabar}}_c}{{\rho }_c}\gamma F\left(\frac{\omega }{{\omega }_c}\right).$$

Full text

††pubyear: 2019††pagerange: Pulsar gamma-ray emission in the radiation reaction regime–Pulsar gamma-ray emission in the radiation reaction regime

1 Introduction

Within the last decade, important progresses have been made towards a better and deeper understanding of pulsar magnetospheric physics, particle acceleration and radiation properties. We have witnessed breakthroughs in numerical simulations of the full non-linear problem of force-free (Contopoulos et al., 1999; Spitkovsky, 2006; Parfrey et al., 2012; Pétri, 2012) and magnetohydrodynamic (MHD) magnetospheres (Komissarov, 2006; Tchekhovskoy et al., 2013). Accurate and detailed multi-wavelength observations of pulsar light-curves especially in the gamma-ray band are now available (Abdo et al., 2013). Gamma-ray pulsars are believed to furnish a faithful indirect view of the pulsar engine because most of their spindown luminosity goes into pulsed gamma-ray emission seen by a distant observer.

On the theoretical side, numerical simulations are now able to include dissipative effects in an heuristic way (Li et al., 2012; Kalapotharakos et al., 2012b) in order to incorporate self-consistently the acceleration zones and therefore localizing the emission sites within the magnetosphere. Nevertheless these simulations are inherently unable to follow single particle acceleration, preventing the building of power law distribution functions as required to fit observations. Recently, Particle In Cell (PIC) codes emerged to fully account for this single particle acceleration and its feedback onto the electromagnetic topology (Cerutti et al., 2015). Unfortunately, in these simulations, the neutron star size is unrealistically large with $r_{\rm L}=3\,R$, with $r_{\rm L}$ the light-cylinder radius and $R$ the neutron star radius. These parameters thus correspond to a sub-millisecond pulsar rotation period. However, the plasma magnetization is more realistic with $\sigma\approx\numprint{e3}$. Such simulations showed that particles reach Lorentz factors up to 103-104. However PIC simulations are still unable to catch neutron star electrodynamics for magnetic fields as high as those present in normal pulsars, that is about $10^{8}$ T. Even millisecond pulsars, believed to harbour fields of only $10^{5}$ T, are difficult to follow faithfully with current simulation techniques because of the large gap between the cyclotron frequency and the pulsar rotation frequency among others. Indeed, at the surface of the star, the ratio between gyro frequency $\omega_{B}$ and stellar rotation $\Omega$ is about

[TABLE]

$e$ is the electron charge, $m_{\rm e}$ its mass, $P=2\,\pi/\Omega$ the pulsar period and $B$ its magnetic field strength. Moreover, the Larmor radius associated to these fields for a particle with Lorentz factor $\gamma$ is

[TABLE]

where $c$ is the speed of light. This lengthscale remains much smaller than the typical size of a neutron star estimated to be about $R=12$ \SIunits k m . Thus the ratio between Larmor radius and neutron star radius is about $\epsilon=r_{\rm B}/R=\numprint{e-10}$, allowing to separate both scale. This clearly show that PIC codes will have tough time to solve the full span of dynamical ranges. Farther away from the star, the situation could get better. Actually, the magnetic field strength decreases quickly with radius like $B\propto r^{-3}$ because of its dipolar nature but for a millisecond pulsar, the ratio at the light cylinder is still

[TABLE]

The same ratio applies for a normal pulsar. Thus even at the light-cylinder, the microscopic and macroscopic timescales are too disparate to be caught by standard numerical techniques.

In order to circumvent these severe limitations, the magnetic field intensity is usually artificially decreased by several orders of magnitude to alleviate the stringent requirement about the time step for integrating the equation of motion. Unfortunately, such artefacts drastically minder the electric field strength too, electric field induced by the rotating dipole, thus disabling particle acceleration to ultra-relativistic speeds with $\gamma\gg 10^{4}$. The highest Lorentz factors obtained so far are usually around $\numprint{e3}-\numprint{e4}$. Consequently, photons are produced with artificially low energies. Kalapotharakos et al. (2018) were able to approach realistic values, using a pulsar period $P=0.1$ s but a still too low magnetic field of $B=\numprint{e2}$ T. In order to approach realistic $B$ fields of $B=\numprint{e8}$ T, they followed the high-energy tail of their particle distribution functions, assuming that their orbits are geometrically correct independently of $B$, integrating the energy conservation equation where electric acceleration is counterbalanced by radiation reaction. This should help to track the Lorentz factor evolution in time for real particles with $\gamma\approx\numprint{e8}$. Nevertheless some correcting factors are introduced to renormalize all energy scales including photon energies, a technique that is at least questionable. This represents the major flaw of direct numerical simulations of neutron star magnetospheres intended to compute realistic spectra from first principle particle in cell simulations.

From a geometrical point of view, several popular emission sites like the outer gaps, the slot gaps and the polar caps were often hypothesized to be privileged regions for producing photons. It is then possible to compare the merit of each zone and test their ability to reproduce the observed light-curves (Dyks & Rudak, 2003; Dyks et al., 2004). The presence of a plasma partially or completely screening the electric field shows up in distortions of the light-curves from a vacuum rotator (see for instance sky maps in Bai & Spitkovsky (2010b) but who surprisingly also assumed some force-free prescription for vacuum fields!) compared to a force-free model (Bai & Spitkovsky, 2010a). Often in the vacuum field investigations, widely used in the literature, the associated accelerating electric fields are not taken into account self-consistently. Let us however mention the work of Kalapotharakos et al. (2012a) and Kalapotharakos et al. (2014), who indeed used resistive plasma models with low conductivity to mimic almost vacuum electromagnetic fields. They also produced sky-maps and light-curves taking into account the accelerating electric field.

Watters et al. (2009) compiled an atlas of geometric light curves for young pulsars showing the essential characteristics of gamma-ray profiles depending on viewing angle and obliquity. Romani & Watters (2010) then designed a tool to constrain the magnetospheric structure from these gamma-ray light curves. Venter et al. (2009) investigated the special population of millisecond gamma-ray pulsars showing that two-pole caustics and outer gap models are favored. See also Pierbattista et al. (2015, 2016) for a large sample of pulsars fitted with several emission models and Johnson et al. (2014) for a similar study about millisecond pulsars. Obviously, more constraints can be obtained from simultaneous radio and gamma-ray fitting. Harding et al. (2011) produced atlases of two-pole caustics and outer gap emission models in force-free and vacuum retarded dipole field geometry to compare light curve features in symmetric and asymmetric slot gap cavities.

PIC codes are now able to follow particle trajectories including radiation reaction correction and therefore producing sky maps and light-curves assuming synchrotron emission. Unfortunately, the magnetic field strength as already pointed out is artificially decreased to too low values. Cerutti et al. (2016) got acceleration only up to $\gamma\approx\numprint{e2}$. This is clearly not enough to explain MeV or GeV photons produced by synchrotron radiation. The current sheet and the Y-point become the preferred site to produce high-energy gamma-ray photons (Philippov & Spitkovsky, 2018), showing light-curve features in agreement with Fermi/LAT observations (Abdo et al., 2013). However the maximum polar cap potential drop they used was set to get Lorentz factor at most of $\gamma=500$. This threshold is many orders of magnitude below any realistic pulsar acceleration efficiency. The neutron star period is also slighty to high with $R=4\,r_{\rm L}$. PIC simulations in their current development stage are unable to deal with real pulsar parameters. The derived spectra and light-curves are therefore also unrealistic as long as the down scaling operates to extrapole dangerously to $\numprint{e5}-\numprint{e9}$ \SIunits T . In all PIC simulations, the hierarchy of time scales is obviously respected but unfortunately not their ratio. We believe that such strong extrapolations must at least be verified on simple problems before dealing with the full complexity of a pulsar magnetosphere.

The second Fermi gamma-ray pulsar catalogue (Abdo et al., 2013) contains plenty of information about gamma-ray pulsars spectra and light-curves. The gamma-ray peak separation clusters around $\Delta\approx 0.5$ and the radio peak usually leads the first gamma-ray peak but with some outliers. Force-free or ideal MHD computations are unable to self-consistently accelerate particles and localize the emission site. Some kind of dissipation of the electromagnetic field is required in order to produce a signal detectable on earth. So dissipation within the magnetosphere and/or wind must occur, but the precise mechanism and its efficiency are difficult to predict from first principles. Nevertheless, some dissipative magnetospheres, called FIDO and introduced by Kalapotharakos et al. (2012b), were used by Brambilla et al. (2015) for computing the phase-averaged and phase-resolved $\gamma$-ray spectra of eight of the brightest Fermi pulsars. They used billions of test particles trajectories to compute curvature radiation spectra in realistic fields of $\numprint{e7}-\numprint{e9}$ T. Based on this work Kalapotharakos et al. (2017) constrained the dissipation mechanism by looking at curvature radiation in the equatorial current sheet outside the light-cylinder, using Fermi/LAT spectral data. This could put limits on the strength of the accelerating electric field. They used test particle integration in the radiation reaction limit regime in the global force-free dissipative magnetosphere which is basically a fluid description avoiding the stringent strong field constrain faced by PIC codes. In such a way, they were able to deduce realistic spectra for realistic pulsar field strength and period. Nevertheless, starting from PIC simulations, Kalapotharakos et al. (2018) recently found a relation between the particle injection rate and the spindown luminosity. This work shows the fruitful feedback between simulations and observations to extract useful information about the nature of particle acceleration and dissipation of the relativistic magnetized flow. Harding (2016) and Venter et al. (2018) gave recent reviews of the successful interplay between magnetospheric modeling and gamma-ray observations.

Other attempts to fit particular pulsars were carried out by other groups. For instance Takata et al. (2007) and Hirotani (2008) used the vacuum retarded dipole to model the outer gap of the Crab pulsar. Du et al. (2011) performed computation in the annular gap context for the Vela pulsar whereas Du et al. (2012) did it for the Crab pulsar. Several millisecond pulsars were also fitted by Du et al. (2013) using a static dipole.

Global magnetospheric simulations converge to a stationary picture of a corotating electromagnetic field and particle distribution function. Nevertheless the paradigm of pulsars being stable and constant broadband emitter in time has been invalidated in radio since their discovery fifty years ago. However, gamma-ray pulsars were though to still remain steady emitter. But this picture has recently also been challenged by some gamma-ray variability reported for instance in PSR J2021+4026 by Allafort et al. (2013).

All the above investigations started from more or less sophisticated numerical simulations of neutron star magnetospheres according to force-free, MHD, resistive/dissipative or PIC approximations. Observational signatures are then post-processed or self-consistently included for comparison with existing data in radio and gamma-rays. This is always the starting point to support any model of pulsar magnetosphere. Force-free magnetosphere simulations give quick and accurate answers to the global electromagnetic field produced by ideal presssureless and massless plasmas. It corresponds to the ultra-strong field limit where particles move at the speed of light. Unfortunately, these simulations cannot resolve for individual particle acceleration. PIC codes are therefore intended to catch all the physics, from macroscopic scales to microscopic scales, self-consistently. This formidable and laudable task is however hampered by the span in time and length scales. This forces simulations to run with unrealistically low values of the electromagnetic fields, which is the major drawback of full particle approaches. There is so far no way out to satisfactorily treat single particle acceleration with radiation reaction self-consistently in ultra-strong electromagnetic fields. This represents a major task towards a deeper and closer investigation of realistic pulsar electrodynamics but so far no numerical technique is able to deal with such regimes.

In the present paper, we decided to start the study of pulsar high-energy emission from a different perspective, trying as much as possible to shortcut any large and time consuming plasma simulations in realistic field strengths. Instead, we require a minimal amount of assumptions putting special emphasize to computing light-curves and spectra with realistic values of the electromagnetic field of about 108 \SIunits T for normal radio pulsars and 105 \SIunits T for millisecond pulsars. In the radiation reaction limit, particles follow a velocity field solely prescribed by the local electromagnetic field itself. It is sometimes called Aristotelian dynamics or zero mass dynamics (Gruzinov, 2013) but it is a simple consequence of particle motion with friction in the ultra-relativistic limit (Mestel, 2012). In section 2, we remind the velocity prescription derived from the radiation reaction limit and the method to compute light-curves and spectra. The limit of applicability of the classical curvature radiation formula is briefly discussed. In section 3 we show detailed results about spectra and sky maps for millisecond pulsars and normal pulsars. Possible future detection in the sub-TeV range from CTA is also briefly investigated. The limit of our approach is discussed in details in Sec. 4. Conclusions are then drawn in Sec. 5.

2 Magnetospheric emission model

We start with a description of the minimalist model used to compute light-curves, sky maps and spectra. Our primary target is to refrain from adding excessive a priori unconstrained parameters into the model in order to catch the essential physics required to fit gamma-ray pulsar data compiled in the second Fermi gamma-ray pulsar catalogue. The master physical quantities are the neutron star period $P$, its period derivative $\dot{P}$ (from which we deduce a fiducial magnetic field strength $B$ at the equator) and the inclination angle between its rotation and magnetic dipole axis depicted by the obliquity $\chi$. Apart from this obliquity $\chi$ which is not constrained by observations, $P$ and $\dot{P}$ are well quantified by pulsar timing campaigns. However, some other inputs are required like the particle distribution function and the extent of the emitting volume. We recall these inputs in the following paragraphs.

2.1 Radiation reaction

Pulsar magnetospheres are filled with ultra-relativistic electron/positron pairs copiously radiating while accelerated by the electric field. It is safe to assume in a first stage that they reach an equilibrium state between acceleration and braking, called radiation reaction limit regime. The photon back reaction onto the particle motion is therefore important. The radiative friction brakes the particle in a direction opposite to its motion such that in the ultra-relativistic limit its velocity depends only on the local value of the electromagnetic field $\mathbf{B}$ and $\mathbf{E}$ (also sometimes called Aristotilean electrodynamics). Electrons and positrons will not react the same way to the electric field $\mathbf{E}$ thus two expressions for the velocity are required. It can be shown, assuming that both particle species speeds are equal to $c$, that the velocity is given by (Gruzinov, 2013; Pétri, 2018)

[TABLE]

where the plus sign corresponds to positrons and the minus sign to electrons. Moreover, we introduced the two electromagnetic invariants $E_{0}$ and $B_{0}$ such that

[TABLE]

with the subsidiary condition $E_{0}\geqslant 0$ ensuring that the radiation reaction force is always directed oppositely to the velocity direction. As explained in Pétri (2018) these invariants are related to the electromagnetic field strength in a frame where $\mathbf{E}$ and $\mathbf{B}$ are parallel. The lepton motion can be decomposed into an electric drift part $\mn@boldsymbol{E}\wedge\mn@boldsymbol{B}$, a motion along magnetic field lines $\mathbf{B}$ and a motion along electric field lines $\mathbf{E}$. This last part of the motion is responsible for dissipation because the power of the Lorentz force is $q\,(\mn@boldsymbol{E}+\mn@boldsymbol{v}_{\pm}\wedge\mn@boldsymbol{B})\cdot\mn@boldsymbol{v}_{\pm}=q\,\mn@boldsymbol{v}_{\pm}\cdot\mn@boldsymbol{E}\geq 0$ where $q=\pm e$ depending on the charge, positron or electron.

In the near field zone, i.e. close to the neutron star surface, where $E\ll c\,B$, the particle velocity simplifies into a motion solely along $\mathbf{B}$ such that

[TABLE]

This expression can be reduced to

[TABLE]

by noting that in this weak electric field limit the magnitude of $\mn@boldsymbol{B}$ is almost equal to the invariant $B_{0}$, namely $B^{2}\approx B_{0}^{2}$. Particles are accelerated mostly by the electric component parallel to the magnetic field. The surface $\mn@boldsymbol{E}\cdot\mn@boldsymbol{B}=0$ are of particular interest because the velocity change sign when the particle cross this region. It is called a force-free surface and represents trapping regions for those particles (Finkbeiner et al., 1989). We return to this important point in Sec. 4.

2.2 Curvature radiation

Particle trajectories can be computed from the velocity field prescription given in eq. (4). These trajectories are obviously bent, leading to curved paths and therefore curvature radiation. The curvature radius $\rho_{c}$ is computed according to the acceleration following the expression

[TABLE]

This acceleration is evaluated by a simple second order finite difference scheme. The associated curvature radiation spectrum for a particle with Lorentz factor $\gamma$ is given in Jackson (2001) by

[TABLE]

where $K_{5/3}$ is the modified Bessel function of order $5/3$, $\varepsilon_{0}$ the vacuum permittivity, $I$ the intensity and $\omega$ the angular frequency. The fundamental frequency is $\omega_{0}=c/\rho_{c}$ and the characteristic curvature photon frequency therefore reads

[TABLE]

from which we deduce the curvature power as

[TABLE]

The curvature emissivity depends on the observation frequency $\omega$ as well as on the location in the magnetosphere $\mathbf{r}$. This emissivity is given by

[TABLE]

showing explicitly the spatial dependence of this emissivity. $\alpha_{\rm sf}$ is the fine structure constant defined by

[TABLE]

with $\hbar$ the reduced Planck constant. Curvature radiation is very similar to synchrotron radiation for which the function F is usually defined by

[TABLE]

The spectra and cut off frequency in both cases are described in Jackson (2001). In the radiation reaction limit, the power exerted by the electric field is simply $\pm e\,\mathbf{v}_{\pm}\cdot\mathbf{E}=e\,c\,E_{0}\geqslant 0$. The work done is always positive as it should be for a dissipative force. According to curvature radiation losses, the maximum Lorentz factor an electron or a positron can reach is

[TABLE]

This equilibrium Lorentz factor $\gamma$ weakly depends on the electric field and curvature radius. Aristotelian electrodynamics implies a particle speed exactly equal to the speed of light, by definition and construction of the velocity given in expression (4). Thus, technically, the Lorentz factor $\gamma$ is computed from the knowledge of the curvature radius deduced from the acceleration of ultra-relativistic particles, eq. (8). We will show that the actual Lorentz factors are $\gamma\gtrsim 10^{8}$ thus widely justifying the approximation of taking $v=c$. In Aristotelian electrodynamics, particles do not have memory about their past trajectory because the velocity is computed according to only the local current electromagnetic field at their position. This locality remains true as long as the particles are able to accelerate due to the electric field or decelerate due to radiation reaction on a length scale $\ell$ much smaller than electromagnetic field gradient and curvature radius. The distance required to gain energy up to $\gamma\,m_{\rm e}\,c^{2}$ is

[TABLE]

With the radiation reaction limiting Lorentz factor, we find

[TABLE]

This ratio is always much less than one for realistic pulsar parameters. In other words, particle emission at some location is not affected by the electric field the particle encountered at another position. It loses its memory within a short distance much smaller than any macroscopic length scale.

The total luminosity radiated by the magnetosphere is therefore

[TABLE]

where $n$ is the particle density number and integration goes along the emitting volume $V$. The inner and outer boundary of the integration is not specified. A natural choice for the minimum radius is the neutron star size and a possible maximum radius is the light-cylinder although other prescription are conceivable. For instance, photon production outside the light-cyliner is another interesting possibility.

The density of leptons is another important unknown parameters. As we want to stay minimalist in our model, we assume a spherically symmetric profile with a decrease in radius according to

[TABLE]

where $n_{0}$ is a normalisation factor and $q$ the exponent of the power law decrease in radius. In the same spirit of simplicity, we do not introduce any power law particle distribution function but straightforwardly choose the local Lorentz factor according to the radiation reaction limit regime prescribed in eq. (15).

2.3 Quantum corrections

Radiation processes are usually derived in a non QED framework where quantum corrections to emission are neglected. Such expressions remain valid as long as the magnetic field strength stays well below the quantum critical field of $B_{\rm qed}=\numprint{4.4e9}$ T. Quantum corrections arises because of the particle recoil and when the photon wavelength becomes comparable to the particle Compton wavelength

[TABLE]

Let us quantify when QED sets in to modify the photon spectra. Curvature radiation is very similar to synchrotron radiation. Both processes originate from the radiation of a charged particle subject to acceleration. The associated photon spectra are therefore similar if the cyclotron gyro-frequency is replaced by the instantaneous rotation frequency of the particle along its curved path. Due to the general law of conservation of energy, no charge can radiate more than its kinetic energy. It is well known that quantum synchrotron sets in whenever the following parameter reaches s close to unity (Erber, 1966; Aharonian et al., 2013)

[TABLE]

We stress that because of the emitting particle Lorentz factor intervening in the above expression, quantum effects manifest already at field strengths much less that $B_{\rm qed}$. This is of primary importance in pulsar magnetospheres because as will be shown later, $\gamma$ can go up to 108-109. Consequently, quantum synchrotron radiation is at work up to very large distances compared to the neutron star radius. In the same vain, looking for the curvature radiation, quantum effects become perceptible whenever the parameter

[TABLE]

approaches unity. Now the $\chi$ parameter is even more sensitive to the Lorentz factor. We will check a posteriori that $\chi_{\rm curv}$ remains weak or at least $\chi_{\rm curv}\lesssim 1$ in all our computations.

2.4 Normalisation

In order to simulate realistic value of electromagnetic field strengths, electron/positron energies and photon energies, we normalise the fundamental quantities of the problem. The magnetospheric distances are normalised to the light-cylinder radius $r_{\rm L}=c/\Omega$. Velocities are normalised to the speed of light $c$. The magnetic field normalisation is performed according to the critical field $\tilde{B}=B/B_{\rm qed}\equiv b$. The electric field typical value is given by Schwinger value of

[TABLE]

such that the normalised value of the electric field becomes $\tilde{E}=E_{0}/E_{\rm Schw}$. The characteristic curvature photon energy in normalised units is conveniently written in units of the electron rest mass energy such that

[TABLE]

The Lorentz factor balancing exactly acceleration against radiation therefore becomes

[TABLE]

The curvature power emitted by a single particle accelerated in the electric field is

[TABLE]

In orders of magnitude, the normalised electric field strength is

[TABLE]

As a characteristic particle number density, we use the expression deduced from the force-free condition, the Goldreich-Julian density, given by

[TABLE]

For the radiative properties, normalizing energies also to the electron rest mass energy, the curvature emissivity is given by

[TABLE]

We introduced normalised frequency and time such that $\tilde{\omega}=\hbar\,\omega/m_{\rm e}\,c^{2}$ and $c\,t=\lambdabar_{c}\,\tilde{t}$. In normalised units, the luminosity becomes by introducing the multiplicity factor $\kappa$

[TABLE]

or explicitly with the spatial dependence of curvature radius and Lorentz factor

[TABLE]

With the normalisation of the density $n_{0}$ we finally get

[TABLE]

This last expression is used to compute the whole information about emission in the magnetosphere. In particular, spectra and light-curves shown in the following section are derived from eq. (32). The flux restored with SI units therefore becomes

[TABLE]

For the neutron star radius, we take a fiducial value of 12 \SIunits k m (Özel & Freire, 2016). The pair multiplicity is fixed to $\kappa=1$ if not otherwise specified. The radial boundary radii are normalized to the light-cylinder, $r_{\rm in}=R_{\rm in}/r_{\rm L}$ and $r_{\rm out}=R_{\rm out}/r_{\rm L}$. Because in the following section there is no confusion possible between electric field $E$ and photon spectra $I_{\rm tot}$, we restore the usual notation, replacing $I_{\rm tot}$ by $E=\hbar\,\omega$ when discussing spectra. We therefore use the conventional notation again like $E^{2}\,dN/dE\,dt$ in eq. (33). Next we show a detailed analysis of the pulsed emission characteristics.

3 Simulations

High energy emission emanates from regions close to the neutron star surface because the electromagnetic field is largest there and therefore the invariant field quantity $E_{0}$ as well as the Lorentz factor required in the radiation reaction limit regime too. On one hand, TeV photons are produced in the innermost part of the magnetosphere. As a general comment, for normal pulsars, the field is strong enough to disintegrate these photons into electron/positron pairs, rendering the medium opaque to this light. Therefore the effective TeV photon flux, if any, is much weaker than in the case of a magnetically optically thin magnetosphere. On the other hand, sub-GeV and MeV photons are produced close to the light-cylinder and freely escape the magnetosphere with a low probability interaction with the magnetic field.

In this section, we show some typical mean and phase-resolved spectra, sky maps and light curves for realistic magnetic field strengths, rotation periods and geometries when particles radiate in the radiation reaction limit regime. We also discuss the cut-off energy and the gamma-ray luminosity dependence on these fundamental parameters. Results are shown for two archetypal classes of pulsars: millisecond pulsars with typical period of $P=5$ ms and normal pulsars with typical period of $P=100$ ms. The magnetic field strength is given in units of the critical field $B_{\rm qed}$ such that we used the normalized field given by the parameter $b=B/B_{\rm qed}$. The Lorentz factor used for beaming in the direction of motion of particle as imposed by eq. (4) is set to $\Gamma=10$. It should actually be beamed into a cone of opening angle $\propto 1/\gamma\ll 1$ but this would require a fantastic angular resolution in the volume integration of eq. (32). In any case, for $\gamma\gg 1$ light-curves and spectra become insensitive to the precise value of $\gamma$. They are shaped by the electromagnetic field topology that is a macroscopic scale.

3.1 High-energy spectra

High-energy spectra are easily compiled by computing the energy flux $E^{2}\,d^{2}N/dt\,dE$ for different energy bands. A typical example of spectra for a normal pulsar is shown in Fig. 1.

The pulsar obliquity is set to $\chi=60\degr$ and the magnetic field strength to $b=10^{-3}$. The particle density profile $n(r)$ is spherically symmetric and decreases with radius according to $n(r)\propto r^{-q}$ with $q$ an arbitrary constant taken for concreteness within the set $q\in\{1,2,3\}$. The spherical symmetry is clearly a crude approximation of the spatial distribution of particles. A better description would require a deeper understanding of pair creation within the magnetosphere. In our minimalistic approach, we bypass such refinements. Moreover, we assume a Goldreich-Julian corotation density normalisation at the surface such that $e\,n(R)=2\varepsilon_{0}\,\Omega\,B(R)$. The energy flux is measured at a distance of 1 kpc. Emissivity occurs within a spherical shell of inner radius $R_{\rm in}$ and outer radius $R_{\rm out}$ not necessarily equal to the light-cylinder radius $r_{\rm L}$. Emissivity is integrated within the volume located between $R_{\rm in}$ and $R_{\rm out}$. Because of magnetic photo-absorption efficiency close to the surface, very high energy photons preferentially come from outer regions such that $r\gtrsim 0.5\,r_{\rm L}$ rather than from regions $r\lesssim 0.1\,r_{\rm L}$. Emission from high altitude magnetospheric sites is also preferred following current wisdom. Although the particle distribution is not monoenergetic (the monoenergetic distribution we enforce is spatially variable due to the variability of the electromagnetic field and the corresponding radiation reaction rate), the volume integrated spectra resemble a monoenergetic distribution, typically a power law with an exponential cut-off. Moreover, the cut-off energy lies around several GeV as seen in the Fermi/LAT pulsar catalog (Abdo et al., 2013) for our special choice of $b=10^{-3}$. The slope of the power law below the cut-off agrees with the $1/3$ exponent of a single particle curvature radiation spectrum. The maximum energy flux at its peak depends on the size of the emission volume as expected, proportional to the total number of particles. Indeed, by inspection of Fig. 1, we deduce that increasing $q$ reduced the total energy flux because the particle density decreases also faster for $q=3$ compared to $q=1$. The location of the inner radius $R_{\rm in}$ impacts only on the shape of the exponential cut-off. This is because the most energetic photons are produced in the strongest accelerating field that is close to the surface. Cutting the emission volume at higher altitude removes these photons from the spectra as expected (compare the red, green and blue lines).

In the wave zone, outside the light-cylinder, when radiation emanates from distances $r>r_{\rm L}$, the average spectra remain very similar to those produced inside the light-cylinder, Fig. 2. The cut-off energy is slightly less outside but still very close to several GeV. Consequently, spectral features are insensitive to the precise extent of the emission regions. However, as will be shown later, light-curves shapes and pulse profiles are sensitive to the location of the photon production sites.

A second example of mean spectra is shown in Fig. 3 for a 5 ms period pulsar with a lower magnetic field of $b=10^{-6}$. Again a power-law with exponential cut-off is observed but with a sharp extinction above the cut-off energy of several GeV. The density profile mainly impacts on the maximum intensity level whereas the extension of the emitting region controlled by $R_{\rm in}$ slightly shapes the cut-off behaviour as for normal pulsars.

Mean spectra highlight the general trend of magnetospheric emission. For a peculiar pulsar, phase-resolved spectra offer valuable insight into the emission region, its shape and geometry within the magnetosphere. Therefore, Fig. 4 shows a phase-resolved spectrum for a 100 ms pulsar with $\chi=60$°, $b=10^{-3}$ and $r_{\rm in}=0.2$.

The total flux variation between the off-pulse and on-pulse peak intensity is about one order of magnitude. The spectral shape remains substantially the same for all phases during its rotation. Very similar results are found for a 5 ms pulsar for which the phase-resolved spectra are given in Fig. 5. We essentially observe the same trend as for the 100 ms pulsar that is similar spectra for all phases but with a shift in magnitude of at most one decade.

After this brief survey on the pulsed spectral features, we dig into the geometrical properties of the light-curves as depicted in sky map diagrams.

3.2 Sky maps

Sky maps are useful graphical representations of the light-curve profiles depending on the obliquity $\chi$ and inclination of the line of sight $\zeta$. Several pertinent subset of sky maps are shown in the following figures. First, Fig. 6 shows a sample of sky maps for a 100 ms pulsar with appropriate magnetic field strength as given in the previous paragraph. The density profile sharply decreases with $q=3$. In each plot of the panel, the $x$ axis depicts the phase of the pulsar (one period normalized to unity) and the $y$ axis depicts the line of sight inclination angle $\zeta$ going from $0\degr$ to $180\degr$.

Each plot is given for a specified energy, increasing from top left to bottom right. Precise values are shown in the top left labels of each map. The intensity level of each map is shown in the colorbox legend on the right. In the power law regime of curvature radiation scaling as $\omega^{1/3}$, sky maps are more or less the same at all energies below the cut-off energy of several GeV, but around and above this cut-off, the light-curves appreciably change their profile due to the exponential tail until a very low imperceptible flux and eventually an extinction at very high energies. A notable difference with respect to all other magnetospheric or wind emission models is that the position in phase of the two peaks is insensitive to the inclination angle $\zeta$. Moreover, the peak separation in phase remains invariably equal to 0.5. This is not typical for Fermi/LAT pulsars. However, such patterns reflects the symmetry of the electromagnetic field dictated by the rotating star. Shifting the location of the magnetic dipole with respect to the geometrical centre of the star would introduce an asymmetry in the field topology around the north and south poles, alleviating the phase separation exactly equal to half a period observed in the present work. Such extension of our model is left for future work but these asymmetries have already been reported in our previous works, showing asymmetric polar cap shapes (Kundu & Pétri, 2017), asymmetric wind structures (Pétri, 2016) as well as asymmetries in the polarization pattern (Pétri, 2017). In our particular set-up, we fixed the origin of phase $\phi_{0}$ in order to locate the peaks around phase $\phi=0.25$ and $\phi=0.75$. This value of the phase origin $\phi_{0}$ is chosen such that both peaks stay visible well within the phase interval $\phi\in[0,1]$. This artificial lag prevents an unintended cut at phase zero or one, that is, right in the middle of a pulse. It is performed just for graphical purposes, avoiding to plot light-curves on two periods as sometimes done for better readability. Obviously, this phase $\phi_{0}$ is arbitrary but with absolutely no impact on light-curves and spectra. Two strong spots are visible around $\zeta\approx\chi$ and $\zeta=\pi-\chi$ when observing below the cut-off. They reflect the location of both polar caps. The situation reverses at the highest energies above the cut-off. The two spots become invisible letting emerge a more diffuse emission away from the polar caps. The emission sites for high energy spread around the outer part of the light cylinder. Note also the drastic decrease in flux of several decades with respect to the low energy part. A second example of sky maps is shown in Fig. 7 for a slowly decreasing density profile with $q=1$. In such a scenario, the high energy flux remains significant and the emission appears less diffuse than for the case $q=3$. But for $q=2$ as shown in Fig. 8 the change in sky maps above several GeV is already apparent and resembles the $q=3$ case. The impact of $R_{\rm in}$ on these same sky maps is also investigated by inspection of Fig. 9 for which $r_{\rm in}=0.1$ and Fig. 10 for which $r_{\rm in}=0.5$. When emission is shifted to the outer parts of the light-cylinder, like in the case $r_{\rm in}=0.5$, the peaks broaden and show a shift with respect to the polar cap location, especially at highest energies above several GeV.

A last example of sky maps is shown in Fig. 11 for a millisecond pulsar. However there is no noticeable discrepancy to discuss between millisecond and normal pulsars.

We finish our discussion about sky maps by exploring their dependence on the inner and outer boundaries $R_{\rm in}$ and $R_{\rm out}$. Indeed, photons of given energies are produced in different radial shells. The sky-maps shown above are the sum of the radiation from all these shells. To better understand the physiognomy of this radiation, we separate the contribution from each spherical shell, assuming a thickness of $\Delta r/r_{\rm L}=0.1$ for each shell. A sample is shown in Fig. 12 for a 100 ms pulsar and a range of $r_{\rm in}$ from 0.1 to 0.9 and $R_{\rm out}=R_{\rm in}+\Delta r$ with $E=511$ MeV. The same plot around the cut-off energy $E=5.1$ GeV is shown in Fig. 13 and well , at $E=51$ GeV in Fig. 14.

Below the cut-off energy, sky maps are very similar, whatever the location of the emitting shell. The overall flux decreases but the peak intensity remains close to the position of the polar caps. However a small spread in the peak profiles is observed when approaching the light-cylinder. In the vicinity of this light-cylinder, electrons and positrons do not contribute symmetrically to the light-curve because the electric field becomes comparable in intensity to the magnetic field. Therefore their velocity field differ significantly, leading to different individual light-curves, see bottom right panel of Fig. 12. The situation is even more prominent at $E=5.1$ GeV, Fig. 13. Above a height of $r/r_{\rm L}>0.7$, significant emission is produced outside the polar cap window, leading to S-shape intensity maps, see the case $r_{\rm in}=0.8$. This leads to a possible phase lag between radio and gamma-ray peaks. At the highest energies, $E=51$ GeV, Fig. 14, a new pulsed component different from the polar cap region appears, especially close to the light cylinder, for $r/r_{\rm L}>0.7$.

Because light-curves are almost energy insensitive across the spectrum, except around the cut-off frequency, it is worth to compute light-curves at a given typical energy of the power law band. Moreover, the contribution from electrons and positrons to the total intensity are usually not symmetrical. Therefore we also show their respective light-curves for an energy in the $\omega^{1/3}$ regime with $r_{\rm in}=0.2$ and obliquities $\chi=\{30\degr,60\degr,90\degr\}$. A normal pulsar is shown in Fig. 15 and a millisecond pulsar is shown in Fig. 16 with respectively $\chi=30\degr$ in the first row, $\chi=60\degr$ in the second row and $\chi=90\degr$ in the third row. The total intensity is shown in red, the electron contribution in green and the positron contribution in blue. Electrons and positrons contribute similarly to the total flux. Three main light-curve profiles are observed. A first class of profiles showing an almost constant intensity where both electrons and positrons contribute in a symmetric manner. Pulsation is therefore very difficult to detect. A second class of profiles showing a prominent single pulse is observed mainly for obliquity much less than $\chi=90\degr$. A third class of double peaked structure is always seen when $\zeta\approx 90\degr$.

Millisecond and normal pulsars show similar trends. Moreover, electrons and positrons almost give the same contribution to the light curves, although some small discrepancies are seen in these particular cases.

If the inner boundary is shifted around the light-cylinder, then the pulses are no more aligned with the location of the magnetic poles. For instance, in Fig. 17, the dominant gamma-ray peak can lead or trail the radio peaks located at phase 0.25 and 0.75 for $E=5.1$ GeV. Electron and positron contributions also differ drastically, leading to highly asymmetric pulse profiles. In almost all cases, emission deviates from zero only around the peaks. Below the cut-off, at $511$ MeV, the light-curves look even more complex, Fig. 18, reflecting the complicated velocity field of the leptons.

From the point of view of possible light-curve profiles, our model can reproduce the same shapes as the one obtained by the competing models mentioned previously. Nevertheless, our new model naturally computes the evolution of a light curve with respect to energy and with realistic magnetic field strengths and rotation periods. Next we explore multi-wavelength light-curves in the following section.

3.3 Multi-wavelength light-curves

Fermi/LAT has shown that double peaked gamma-ray pulsar light-curves evolve with increasing energy towards a dominance of one pulse over the other and a possible shrinking of the pulse width. In the context of our model, we investigate the light-curve evolution with energies from MeV up to sub-TeV, picking out a subset of simulation parameters similar to the previous ones such that $\chi=60\degr$, $r_{\rm in}=0.2$, $r_{\rm out}=1$, $q=\{1,3\}$ and $\zeta=40\degr$. For a normal pulsar, the evolution with photon energy $E$ is shown in Fig. 19. In this particular case, the second peak is dominant at low energy $E\lesssim 10$ GeV, both peaks become equal in intensity at $E\approx 10$ GeV and the first peak dominates above $E=50$ GeV. Note also that the first pulse, almost undetectable and wide at lowest energy becomes intense and sharper at the highest energies. These conclusions do not depend on the density profile, the $q=1$ case on the right column shows similar trend as the $q=3$ case on the left column.

For a millisecond pulsar, the first peak intensity increases with photon energy too, as for the normal pulsar, but at the highest energies, the second pulse remains dominant, the first one almost disappearing, Fig. 20. Here also, not much differences are reported between the $q=1$ and the $q=3$ cases. Consequently, several scenarios are possible for main pulse/interpulse dominance depending on period and magnetic field strength.

The variation in the peak intensity ratio is a consequence of the particular behaviour of the phase-resolved spectra presented in section 3.1. In Fig. 4 for the normal pulsar, the second peak overtakes the first peak in the last visible point in $\log(E/MeV)$ between [4,4.5]. Such overtaking is not seen in the millisecond pulsar case plotted in Fig. 5.

The values of magnetic field strength and rotation period are adjusted to fit the second Fermi/LAT pulsar gamma-ray catalogue, especially the cut-off energy. But how are these cut-off energies related to the fundamental parameters of our model? This is the question we investigate in the next section.

3.4 Cut-off energy

The cut-off energy $E_{\rm cut}$ is usually obtained by fitting the Fermi/LAT spectra by a power-law with exponential or sometimes sub-exponential cut-off. As such fits seem not very robust against the sub-exponential coefficient and are subject to errors related to the power law index used below the cut-off, we prefer to use a more robust value simply by looking for the maximum in the spectral energy distribution. This definition is independent of any assumption about the fit. This maximum in the spectral energy distribution is similar to the apex energy defined by Renault-Tinacci et al. (2015). Results for the energy at the maximum flux $E_{\rm cut}=k_{\rm cut}\,m_{\rm e}\,c^{2}$ are shown in Fig. 21 for a normal pulsar and in Fig. 22 for a millisecond pulsar.

Both cut-off energies follow the law

[TABLE]

derived from eq. (10), (15) and (27). The zig-zag curve is an artefact due to the finite number of frequency points $\epsilon$ used in regular intervals of $1/2$ in a logarithmic scale of $\log(\epsilon/m_{\rm e}\,c^{2})$. Both cut-off energies follow the $\tilde{B}^{3/4}$ law but with different numerical constant in front of it related to the $R/r_{\rm L}$ ratio and to the curvature $\rho_{c}$. Note also that the curvature $\rho_{c}$ is insensitive to the pulsar period or in other words insensitive to the ratio $R/r_{\rm L}$ because it is related to the electromagnetic field topology that is almost the same whatever the pulsar period constrained by the ratio $R/r_{\rm L}$. At large distances, it reduces to a plane electromagnetic wave as shown in Pétri (2015) with a relative amplitude between electric part and magnetic part independent of $\Omega$. Indeed, inspecting Fig. 23 where two maps of curvature radius $\rho_{\rm c}$ are shown, one for a 100 ms pulsar on the upper panel and one for a 5 ms pulsar on the lower panel, we do not observe any significant difference, whether close to the surface nor around the light-cylinder or beyond it. Consequently, the cut-off scale as $k_{\rm cut}\propto(B/r_{\rm L})^{3/4}\propto(\Omega\,B)^{3/4}$ for any pulsar and because this product $\Omega\,B$ is very similar for all pulsars, we do not expect a large spread in cut off energies. Noting that the product $\Omega\,B$ is proportional to $\sqrt{\dot{P}/P}$ (assuming that $B\propto\sqrt{P\dot{P}}$), the cut off scales with the characteristics $\tau_{\rm c}$ according to $k_{\rm cut}\propto\tau_{\rm c}^{-3/8}$.

To finish our discussion, we need to adjust the particle density number to arrive at the correct energy flux and at the total gamma-ray luminosity. This is done in the next section, following the data from the second gamma-ray pulsar catalogue of Fermi/LAT (2PC).

3.5 Luminosity

Pulsar gamma-ray luminosities in the 2PC are estimated by integrating the flux in the energy range 100 MeV-100 GeV. Because of the sharp spectral cut-off around 1-5 GeV, the upper bound of this range, when well above this cut-off energy, does not impact on the total luminosity. A crude guess of the gamma-ray luminosity is simply given by the total flux radiated at this cut-off energy. In order to compare with Fermi/LAT data, we compute the total gamma-ray luminosity between 161 MeV and 161 GeV. The upper bound difference has no impact on the real gamma-ray luminosity as the cut-off is well below 100 GeV. The lower bound taken to be 161 MeV instead of 100 MeV has also little impact on the luminosity estimate as the spectrum peaks around 1-5 GeV. The value of 161 comes from our energy discretization grid which is uniform on a $\log$ scale thus getting sky-maps computed at energies in the form $10^{a/2}\times m_{\rm e}\,c^{2}$ where $a$ is an integer. Thus energies are decades in $\sqrt{10}\,m_{\rm e}\,c^{2}=1.61$ MeV or decades in $m_{\rm e}\,c^{2}=511$ keV.

In order to compute the total kinetic rotational energy losses, we remember that the spindown luminosity depends on the period $P$, its derivative $\dot{P}$ and the stellar moment of inertia $I$ according to

[TABLE]

From the magnetodipole losses, the magnetic field strength is related to pulsar timing by

[TABLE]

thus the spindown is evaluated to

[TABLE]

Fig. 24 shows the gamma-ray luminosity for normal and millisecond pulsars depending on the magnetic field strength. In all points, we used an obliquity $\chi=60\degr$ and a density profile with $q=3$. The plus sign $+$ depicts normal pulsars and the cross sign $\times$ depicts millisecond pulsars. The luminosity depends on $B^{2}$ for both kind of pulsars. For a fixed magnetic field strength, millisecond pulsars are more luminous than normal pulsars. This is because the cut-off energy in millisecond pulsar is unrealistically high compared to normal pulsars for a same magnetic field strength $B$. Moreover, according to the Fermi/LAT gamma-ray pulsar catalogue (Abdo et al., 2013), millisecond pulsars luminosities are mostly in the range $10^{25}-10^{27}$ W, requiring magnetic field strengths of $B=10^{-5}-10^{-6}\,B_{\rm qed}$, in accordance with observations. On the other side, normal pulsars possess luminosities mostly in the range $10^{26}-10^{30}$ W, requiring magnetic field strengths $B=10^{-2}-10^{-3}\,B_{\rm qed}$, also in accordance with current wisdom. The inflection of the curve at low and high magnetic field strength $b$ is an artefact due to the cut-off energy being respectively well below or well above the range used to computed the gamma-ray luminosity which is [161 MeV,161 GeV]. Correcting for this effect, the gamma-ray luminosity is proportional to the square of the magnetic field strength and shown by the line $L_{\gamma}=10^{\eta}\,b^{2}~{}\SIunits{\mathrm{W}}$ in the plot, with $\eta\gtrsim 30$.

Fig. 25 shows the gamma-ray luminosity for normal and millisecond pulsars depending on the spindown luminosity $\dot{E}$. Lines $L_{\gamma}=\eta\,\dot{E}$ are also shown for reference. For emission starting at an altitude $h>0.1\,r_{\rm L}$ the gamma-ray luminosity never exceeds the spindown power as required from the basic principle of energy conservation. It is actually 2 to 3 orders of magnitude less than $\dot{E}$. To increase $L_{\gamma}$, we can invoke the pair multiplicity factor $\kappa$ constraining it to $\kappa=\numprint{e2}-\numprint{e4}$ to reconcile this plot with Fermi/LAT second catalogue, depending on the particular pulsar fitted. Consequently, our simple model reproduce the main spectral and timing properties of gamma-ray pulsars without any violation of basic physical principles. Here again, the inflection at high spindown luminosities $\dot{E}$ is an artefact. Correcting for this effect, the gamma-ray luminosity is proportional to the spindown power for both millisecond and normal pulsars, shown by the line $L_{\gamma}=\eta\,\dot{E}$ in the plot, with $\eta\leq 1$.

3.6 Detection of the VHE component

Fermi/LAT already reported more than 200 gamma-ray pulsars. It offers a good sample to extrapolate emission at sub-TeV energies in the window of atmospheric Cerenkov telescopes such as HESS-II who detected Vela above 20 GeV (Collaboration et al., 2018) and the upcoming Cerenkov Telescope Array (CTA) observatory (Vercellone, 2014). It is therefore opportune to show expectation of gamma-ray fluxes in the sub-TeV range for emission in the radiation reaction limit. A representative sample of spectra for normal and millisecond pulsars is shown in Fig. 26 with comparison to Fermi and CTA sensitivities. The pair multiplicity is $\kappa=1$ by default, the obliquity is $\chi=60\degr$. Emission within the magnetosphere is reported as (m) meaning $r_{\rm in}=0.5$ and $r_{\rm out}=1$ whereas wind emission is reported as (w) meaning $r_{\rm in}=1$ and $r_{\rm out}=5$. Normal and millisecond pulsars with respectively $\log(b)=-2$ and $\log(b)=-5$ are marginally detectable above 100 GeV with CTA South in 50h whatever the emission location, within the magnetosphere or within the wind. Lower magnetic fields drastically reduce the flux as well as the cut-off energy. A more realistic multiplicity $\kappa\gg 1$ would increase the flux but not the cut-off. We do not expect any emission above several hundreds of GeV, irrespective of the leptons distribution function in space and momentum. Nevertheless our crude analysis still requires a precise and careful analysis for individual pulsars to make clear and sensible predictions in the sub-TeV range.

4 Discussion

The light-curves and spectra exposed in our simplistic view of pulsar magnetospheric acceleration and radiation mechanisms gives already interesting results worthwhile to extend although we used very little inputs and fitting parameters except for those well constrained by observations such as the period $P$, its derivative $\dot{P}$, and the neutron star radius $R$. There is no doubt that our approach could benefit from some improvement to better fit particular pulsars. Nevertheless, in the next section we focus on some problems and non elucidated electrodynamics about pulsar magnetospheres related to the complex geometry of trajectories and particle behaviour within this magnetosphere, pointing out the limitation of our work.

Our model deals with realistic pulsar periods and magnetic fields, producing spectral high-energy features in agreement with gamma-ray observations performed by the Fermi/LAT (Abdo et al., 2013). We are also able to produce double-peaked light-curves with variable peak intensity ratio, variable shapes and widths. Nevertheless, these light-curve profiles rely heavily on the underlying spatial lepton distribution within the magnetosphere. We certainly took a too simplistic view of spherically symmetric repartition decreasing in radius with a simple power law in radius depicted by $n\propto r^{-q}$. Assuming emission emanating from any point is admittedly too a crude approximation. However, the energy balance between electric acceleration and radiation reaction is not impacted by the spatial distribution of particles except for some corrections due to the back reaction of the plasma onto the field. However, the electric current generated by the plasma flow close to the surface remains too weak to appreciably perturb the electromagnetic field. We conclude that the results obtained about spectral shape and cut-off is robust and insensitive to the geometry except for small changes imprinted by the curvature radius.

4.1 Particle flow and trapping

A proper account of particle flow within this magnetosphere according to the radiation reaction limit prescription requires to solve for the particle number conservation law supplemented with an appropriate source term of electron/positron pair formation. The arbitrariness of our spatial distribution would then be transposed to the arbitrariness of pair creation efficiency. We could get various spatial distribution by changing this source function. We could impose pair cascading only around the polar caps, or spread out over the whole neutron star surface or even within some special regions in the magnetosphere up to the light-cylinder. Clearly, such conclusions would not be as robust as on the energy budget. Solving for the pair formation is a difficult task about the microphysics which inevitably translates into a geometrical problem of localising the source of leptons. It immediately reflects into the light-curve profiles as an observable. We will not go into such refinement but stress that the damped motion implied by radiation reaction produces three kind of particle flows:

outflowing particles escaping the neutron star and its magnetosphere, forming the base of the pulsar (striped) wind. 2. 2.

trapped particles, staying in a defined region close to the surface for a long time with respect to the pulsar period. These regions are identified as $\mathbf{E}\cdot\mathbf{B}=0$ surfaces or volumes and called force-free surface. 3. 3.

returning particles that hit the neutron star surface.

These classes of trajectories are found by a direct numerical integration of the velocity eq. (4) assuming a background Deutsch field.

The outflowing, trapped or returning motion depends on the initial position where the particle has been launched. The overall maps obtained by integration show large inner volumes where particles return to the star and sometimes are trapped. As a general trend however, for the outer most regions far from the stellar surface or well beyond the light cylinder, particles always escape to infinity because they essentially feel a plane vacuum electromagnetic wave with $\mathbf{E}\cdot\mathbf{B}\approx 0$. A typical cross section of these regions in the meridional plane $xOz$ is shown in Fig. 27 for several inclination angles with $\chi=\{0\degr,60\degr,90\degr\}$. The initial particle positions for escaping are shown with red dots, for trapping particles with green dots and for returning particles with blue dots. Obviously, electrons (left column) and positrons (right column) do not share the same returning and outflowing regions. Trapped regions can be large too and being trapped is the privilege of only one species when the geometry is close to an aligned or to a counter-aligned rotator. We emphasize that these regions do not correspond to places where particles are actually trapped but to the starting point of the trajectory for which particles move to trapping regions. The same interpretation holds for escaping and returning trajectories.

Trapped regions are spread around the so-called force-free surface defined by $\mathbf{E}\cdot\mathbf{B}=0$. The three different regions are not the same for electrons and positrons. We note however that the situation is symmetric with respect to the obliquity $\chi$ and $\pi-\chi$, meaning that the cavities for electrons with obliquity $\chi$ are the same as the cavities for positrons with obliquity $\pi-\chi$. If the content in electrons and positrons within the magnetosphere is not exactly the same, it is in principle possible to differentiate the obliquity $\chi$ from $\pi-\chi$ in the sky-maps and light-curves, departing therefore from the traditional symmetry with respect to the equatorial plane.

4.2 Particle trajectories

A detailed study of particle motion in the Deutsch electromagnetic field is out of the scope of this paper. However several interesting works about trajectories in this field or in the near quasi-static zone have been discussed in the literature. Some details can be found for instance in Laue & Thielheim (1986) for a perpendicular rotator with $\chi=90\degr$. The importance of trapping regions for arbitrary obliquity $\chi$ is emphasized by Finkbeiner et al. (1989). It is possible that a fraction of particles stay in the electrospheric configuration given by a dome+torus geometry Jackson (1976); Krause-Polstorff & Michel (1985); Pétri et al. (2002). Zachariades & Jackson (1989) even found bounded trajectories outside the light-cylinder showing the complexity of defining particle motion and density in the surrounding of neutron stars when damping is included.

To better grasp the complexity of trajectories available for particles already in our simple vacuum model, we computed the long term motion of particles launched from the surface of the star. This should mimic the effect of pair creation in the vicinity of the star. Already at this surface, three kind of trajectories emerged as explained above: escaping, trapping and returning particles. This is shown in Fig. 28 where the initial position of particles (on the stellar surface) is represented in spherical coordinates $(\phi,\theta)$. Escaping particles are depicted by red dots, trapping particles by green dots and returning particles by blue dots. Two symmetrical cases are shown, one for electrons $e^{-}$ with $\chi=60\degr$, left column, and the other for positrons $e^{+}$ with $\chi=120\degr$. The species of the opposite charge immediately returns to the star. We therefore obtained a dynamical structure that is not symmetric with respect to the angle $\chi$ and $\pi-\chi$. If the content in electrons and positrons differs to each other, we expect to get different magnetospheric configurations and therefore different emission properties (light-curves and spectra) with respect to the two geometries $\chi$ and $\pi-\chi$. The initial position of escaping, trapping and returning particles builds complicated shapes far from the standard polar cap geometry. It is actually related to the full electromagnetic field and not only to the magnetic part.

A sample of trajectories for escaping electrons with $\chi=60\degr$ is shown in Fig.29 where the anisotropic character of the filling is clearly visible. The left plot shows the projection onto the equatorial plane $xOy$ whereas the right plot shows the projection onto the meridional plane $yOz$. Close to the star, their trajectories can be complicated but outside the light-cylinder, they become almost radial as the electromagnetic wave tends to a plane wave propagating in the direction $\mn@boldsymbol{e}_{\rm r}\approx\mn@boldsymbol{n}\propto\mn@boldsymbol{E}\wedge\mn@boldsymbol{B}$.

At large distances, well outside the light cylinder, the distribution of leptons is concentrated in specific sky directions as shown in Fig. 30. Some regions are devoid of electrons whereas other regions are devoid of positrons. Many electron trajectories tend to preferred directions in the sky as shown by these maps.

Studying single particle trajectories in a background electromagnetic field of a rotating magnetized neutron star is a full topic by itself. Such refinement must be included in a comprehensive description of pulsar electrodynamics but the scope of this paper was to focus mainly on high-energy emission from a simplistic model without resorting to large scale particle simulations.

4.3 Invariants and particle dynamics

In order to better understand the radiative properties of the magnetosphere, we plot several important geometrical and dynamical properties in the meridional plane $xOz$. Normalized units are used as explained in the previous section. The electromagnetic invariants, $E_{0}$ and $B_{0}$ intervening in the velocity field are shown in the upper panel of Fig. 31 on a log scale, for a pulsar with $R=0.1\,r_{\rm L}$, $b=10^{-3}$ and $\chi=60\degr$. Because $B_{0}$ can be of either sign, we plot $\log(B_{0})$ for $B_{0}>0$ and $-\log(-B_{0})$ for $B_{0}<0$. This helps to identify the location where $B_{0}$ abruptly changes sign. In the plot, a negative invariant $B_{0}$ corresponds to $\log(-B_{0})>0$. The highest values of $E_{0}$ and $B_{0}$ are observed close to the neutron star. $B_{0}$ changes sign in the region around $x=z$. This implies a discontinuity in the velocity field eq. (4) along the direction of $\mathbf{B}$. The middle panel shows the curvature radius $\rho_{c}$ normalized to the light-cylinder radius $r_{\rm L}$, for electrons on the left and for positrons on the right. The curvature radius goes from $\rho_{c}/r_{\rm L}=10^{-2}$ to $\rho_{c}/r_{\rm L}=10^{3}$ for both species. It is shortest close to the star and along the rotation axis. The lower panel shows the Lorentz factor for electrons on the left and for positrons on the right. They attain similar speeds from $\gamma=10^{6.5}$ to $\gamma=10^{9.5}$, the highest values being obtained very close to the star along an axis inclined with respect to the rotation axis. From the knowledge of the curvature radius and the Lorentz factor, we compute $\chi_{\rm curv}$ and check that quantum corrections to curvature radiation remain negligible.

4.4 Influence of the electric charge

The total charge of the neutron star is an important parameter to determine the electric field at long distance. It produces a monopolar component decreasing very slowly with radius, therefore producing sensitive effects even around and outside the light-cylinder. Here we do not report on the full impact of this charge on the neutron star electrodynamics and radiation. We stress that assuming a dipolar field inside the star, an electric charge given by

[TABLE]

is located at the centre. The characteristic electric charge scale is

[TABLE]

A net charge shifts the mean spectra to higher photon energies making the cut-off less sharp. A typical example is shown in Fig. 32 where such shift is clearly seen by comparison with Fig. 1. The net charge also influences light-curves and spectra. To remain brief, we show a small sample in Fig. 33 to be compared with Fig. 10. If high energy emission emanates from regions around or beyond the light-cylinder, following current wisdom, this charge must be included for a self-consistent picture of spectra and light-curves. Our findings urge us to better take care of this electric charge breathing when modelling pulsar radiation. It opens up another road towards a better understanding of neutron star electrodynamics and on its intrinsic temporal variability as already demonstrated by numerical simulations of time-dependent pair creation as found by Timokhin (2010) and by Timokhin & Arons (2013).

4.5 Magnetic field strength estimates

In order to obtain realistic spectra with cut-off energies around several GeV for gamma-ray pulsars, we need to fix the magnetic field strength at the surface. We showed that the cut-off scales as $B^{3/4}$ thus it is possible to retrieve any cut-off value by simply adjusting the magnetic field strength $B$. In this work we used field strengths that seem slightly underestimated compared to what is usually assumed from magnetodipole losses. Actually, the cut-off also depends on the location of the inner boundary $r_{\rm in}$ where gamma-ray photons start to escape the magnetosphere without being magnetically absorbed. This effect is shown in Fig. 34 for a 100 ms pulsar with a density profile $n(r)\propto r^{-3}$, a variable inner boundary $r_{\rm in}$ and an outer boundary $r_{\rm out}=2$. The magnetic field strength is $b=10^{-3}$ and the obliquity $\chi=60$°. It is seen that an increase in $r_{\rm in}$ implies a decrease in cut-off energy. Switching from $r_{\rm in}=0.1$ to $r_{\rm in}=1.0$ the cut-off decreases by one order of magnitude. Therefore if emission starts only around the light-cylinder, the magnetic field strength must be augmented by at least a factor ten. Thus our estimates become closer to traditional field estimates from vacuum magneto-dipole losses (although that such estimates are not necessarily realistic when plasma, wind and multipolar components are taken into account). Precise values of $B$ would require fitting spectra on a case-by-case basis for each pulsar. This is however left for future work.

5 Conclusions

We have shown that a simple magnetospheric emission model can account for the gross features of gamma-ray pulsars light-curves and spectra. The spectral shape as reported in the second Fermi/LAT catalogue, cut-off energies and fluxes, are retrieved with realistic neutron star parameters such as its period and magnetic field strength consistent with millisecond and normal pulsars. Moreover light-curves with single or double peaked profiles are obtained depending on the viewing angle and obliquity. Although the gamma-ray luminosity falls well below the $L_{\gamma}=\dot{E}$ line in some cases, introducing a pair multiplicity factor $\kappa$ much larger than unity increases the total power radiated by the magnetosphere. We found estimates of the order $\kappa=\numprint{e2}-\numprint{e4}$.

Single particle trajectories have been computed, showing the complexity of escaping, trapping and returning motion allowed within the magnetosphere. Such trajectories induce complicated geometries for possible vacuum gaps and filled regions that require further and deeper investigation to fully understand their impact on real pulsar electrodynamics. Because in the radiation reaction limit the velocity field depends also on the electric field, we expect variation in light-curves and spectra due to fluctuating electric charge within the magnetosphere. The particle outflow need not be stationary neither exactly compensating one charge escape by the other charge escape. Thus a kind of magnetospheric breathing is induced, impacting also on the pair formation rate. The total charge of the neutron star indeed affects the spectra and light-curve as shown in depth in Pétri (2018).

Moreover, this study suffers from several flaws that need to be fixed in forthcoming works. Firstly, a mono-energetic particle distribution function is only able to reproduce precisely a small samples of the Fermi/LAT gamma-ray spectra below the cut-off energy. This is because the measured power-law spectra require a power-law distribution of emitting particles not restricting the spectra to a simple $\omega^{-1/3}$ law. But this introduces one more free parameter to our model, thus opposite to the philosophy we followed here for our minimalist model. Nevertheless such studies are planed in the future to fit several samples of millisecond and normal pulsars that do not belong to fluxes depicted by mono-energetic lepton distribution functions. Secondly, back reaction of the plasma onto the intially vacuum rotating electromagnetic field must be added, especially close and behind the light-cylinder where this retroaction is preponderant. Kinetic simulations are therefore unavoidable but, unfortunately, such codes are still unable to catch the full span of length and time scales, going from the Larmor frequency to the neutron star rotation frequency. However, this is absolutely compulsory to reach Lorentz factors as high as $\gamma=\numprint{e9}$ and thus realistic photon energies.

Another interesting class of gamma-ray emitting pulsars are the soft gamma-ray pulsar population discussed in (Kuiper & Hermsen, 2015). These pulsars must also be fitted by the same model, looking for spectra and light-curves. Moreover performing some predictions about phase-resolved polarization in high-energy, notably in X-rays in view of the coming IXPE mission (Weisskopf et al., 2016) will better constrain the location an topology of the photon production sites.

Acknowledgements

I am grateful to the anonymous referee for its useful insight and comments that improved the paper quality. This work has been published under the framework of the IdEx Unistra and benefits from a funding from the state managed by the French National Research Agency as part of the investments for the future program. It also benefited from grant No. ANR-13-JS05-0003-01 (project EMPERE) and from the computational facilities available at Equip@Meso (Université de Strasbourg). It also benefited from a CEFIPRA grant IFC/F5904-B/2018.