# Pulsed high-energy $\gamma$-rays from thermal populations in the current sheets of pulsar winds

I. Arka, G. Dubus

► To cite this version:

I. Arka, G. Dubus. Pulsed high-energy $\gamma$-rays from thermal populations in the current sheets of pulsar winds. Astronomy and Astrophysics - A&A, 2013, 550, 10.1051/0004-6361/201220110. insu-03617519

HAL Id: insu-03617519

https://insu.hal.science/insu-03617519

Submitted on 24 Mar 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pulsed high-energy γ-rays from thermal populations in the current sheets of pulsar winds

I. Arka and G. Dubus

Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), UMR 5274, 38041 Grenoble, France
c-mail: arkai@obs.ujf-grenoble.fr

Received 26 July 2012 / Accepted 27 October 2012

## ABSTRACT

**Context.** More than one hundred pulsars have been detected up to now at GeV energies by the Large Area Telescope (LAT) on the Fermi gamma-ray observatory. Current modelling proposes that the high-energy emission comes from outer magnetospheric gaps, but radiation from the equatorial current sheet that separates the two magnetic hemispheres outside the light cylinder has also been investigated.

**Aims.** We discuss the region outside the light cylinder, the “near wind” zone. We investigate the possibility that synchrotron radiation emitted by thermal populations in the equatorial current sheet of the pulsar wind in this region can explain the lightcurves and spectra observed by Fermi/LAT.

**Methods.** We used analytical estimates as well as detailed numerical computation to calculate the γ-ray luminosities, lightcurves, and spectra of γ-ray pulsars.

**Results.** Many of the characteristics of the γ-ray pulsars observed by Fermi/LAT can be reproduced by our model, most notably the position of these objects in the *P* − *P̂* diagram, and the range of γ-ray luminosities. A testable result is a sub-exponential cutoff with an index *b* = 0.35. We also predict the existence of a population of pulsars with cutoff energies in the MeV range. These have systematically lower spindown luminosities than the Fermi/LAT-detected pulsars.

**Conclusions.** It is possible for relativistic populations of electrons and positrons in the current sheet of a pulsar's wind immediately outside the light cylinder to emit synchrotron radiation that peaks in the sub-GeV to GeV regime, with γ-ray efficiencies similar to those observed for the Fermi/LAT pulsars.

**Key words.** pulsars: general – radiation mechanisms: non-thermal – stars: winds, outflows – gamma rays: stars – relativistic processes

## 1. Introduction

Since the launch of the Fermi space telescope, the sample of gamma-ray pulsars has grown to include more than one hundred objects. All these pulsars, whether young or millisecond, single or in binaries, show spectra consistent with power laws with exponential cutoffs, whose cutoff energy lies in the range 1–10 GeV (Abdo et al. 2010). A prominent exception is the Crab pulsar, which has been detected by ground-based Cerenkov arrays in the TeV regime, and whose spectrum can be well fitted by a broken power law (Aliu et al. 2008; McCann 2011). These observations point to the outer gap/slot gap models as the most probable explanation for the pulsed gamma-ray emission, based on population prediction statistics and lightcurve modelling (Gonthier et al. 2010; Decesar et al. 2011; Pierbattista et al. 2011; Venter et al. 2009; Watters & Romani 2011; Venter et al. 2012), although for millisecond pulsars one has to invoke non-dipolar field geometries or displaced polar caps to reach the required energies (Harding & Muslimov 2011). In those models, the emission originates within the light cylinder, defined by the cylindrical radius $r_{LC} = c/ω$ from the pulsar's rotational axis, where a corotating particle would reach the speed of light. However, it was pointed out by Bai & Spitkovsky (2010) that the emission region must be extended slightly outside the light cylinder to reproduce the double-peaked light curves using the magnetic field configuration from self-consistent, force-free simulations of pulsar magnetospheres.

The idea that high-energy pulsations might come from current sheets near or outside the light cylinder is not new. Lyubarskii (1996) predicted that particles accelerated through reconnection close to the light cylinder might emit gamma-rays through the synchrotron process. Lyubarskii's emission site is the point where the warped equatorial current sheet that separates regions of opposite magnetic field polarity in the pulsar's wind meets the current flowing to/from the pulsar's polar caps, the so-called Y-point (Spitkovsky 2006). Later it was pointed out that pulsed emission naturally arises from the periodicity of the pulsar wind, which is modulated by the star's rotation in combination with the large bulk Lorentz factors of the outflow (Kirk et al. 2002). This model has been successfully applied to calculate the polarization of optical emission from the Crab pulsar (Pétri & Kirk 2005) and it has also shown promise in explaining the gamma-ray light curves observed by Fermi/LAT (Pétri 2011). In these models the emission is attributed to the inverse Compton process and starts in the wind region far from the light cylinder, $r \gg r_{LC}$. In the same context it was proposed that the gamma-ray radiation of the Fermi/LAT band could be synchrotron radiation from power-law electrons in the pulsar wind's equatorial sheet far from the light cylinder (Kirk et al. 2002; Pétri 2012).

There is, however, no obvious reason why emission should be truncated outside the light cylinder (for magnetospheric gap models) or start at a minimum radius in the far wind zone (for the striped wind model). The emission region should evolve continuously from the outer magnetospheric gaps through the
light cylinder and into the equatorial current sheet, possibly affected by the reconnection process that will inevitably occur in the sheet (Bai & Spitkovsky 2010). The region beyond the light cylinder but still in the near wind zone is the area that we are trying to explore. The purpose of this paper is to demonstrate how emission in the 100 MeV–100 GeV range explored by Fermi/LAT can naturally arise from thermal populations of particles in the striped wind, without going into the details of reconnection physics, but just using a few simple assumptions about the local description of the current sheet. In Sect. 2 we describe our model, including our assumptions, and give some analytical estimates for the emitted radiation. In Sect. 3 we present examples of emission maps and phase-averaged spectra that arise from our model and explore the parameter space to arrive at more general conclusions about the pulsar population observed by Fermi/LAT. Finally we discuss our results and the potential of a more detailed description of the current sheet outside the light cylinder as a source of pulsed energetic radiation.

## 2. Equatorial current sheet

### 2.1. Description of the particle population

To describe the wind outside the light cylinder, we will use the “slow-rotator” solution to the magnetohydrodynamics equations that govern the physics of the pulsar wind, which was found by Bogovalov (1999). This solution refers to an electron-positron wind launched by a rotating neutron star, the magnetic axis of which is at an angle $\chi$ to its rotational axis. In this description the wind is purely radial and super-fast magnetosonic at launch, and the field can be described by a radial and an azimuthal component of magnitude

$$B_r = \frac{B_{LC}}{R^2} \qquad (1)$$

$$B_{\varphi} = \frac{B_{LC}}{\beta R} \sin \vartheta, \qquad (2)$$

where $R = r/r_{LC}$ is the spherical radius normalized to the light cylinder, $\vartheta$ is the polar angle, $\beta = (1 - 1/\Gamma^2)^{1/2}$ is the bulk speed of the wind, normalized to the speed of light, and $B_{LC}$ is a fiducial magnetic field magnitude at the light cylinder. The adequacy of this solution in describing approximately the structure of the wind for radii as small as the light cylinder has been confirmed by force-free simulations of pulsar magnetospheres (Spitkovsky 2006; Bai & Spitkovsky 2010).

A prominent feature of Bogovalov’s solution is the warped equatorial current sheet, separating the two magnetic hemispheres, across which the field changes sign. In the mathematical solution the sheet is just a discontinuity, but in reality it should have a finite thickness. Such a current sheet is populated by hot particles, whose pressure balances the magnetic pressure of the cold, strongly magnetized plasma outside the sheet. The sheet oscillates in space and time with a wavelength of $\lambda = 2\pi\beta/\omega$, and with the pulsar’s frequency $\omega = 2\pi/P$, where P is the pulsar’s period measured in seconds.

We assumed that in the wind frame (the frame that propagates radially outwards with a Lorentz factor equal to the bulk Lorentz factor of the wind, $\Gamma$), the sheet can be locally described by the relativistic Harris equilibrium (Hoh 1966). For this description to be valid, the segment of the sheet under consideration should be approximately flat. For a relativistic outflow, the hydrodynamically causally connected region of the flow has an opening angle of roughly ~1/Γ, centred on the direction of

motion of the outflow. Since for a pulsar wind $\Gamma \gg 1$, the hydrodynamically connected sheet segment can be considered to be flat and the local Harris sheet description should be a good approximation.

We denote quantities measured in the wind frame with a prime. If we denote the sheet normal direction in the wind frame with a capital $X'$, and the sheet midplane is at $X' = 0$, then locally the field inside and outside the segment can be described by a tangent hyperbolic profile

$$B'(X') = \pm B'_{0} \tanh \left( \frac{X'}{\delta'} \right), \qquad (3)$$

where $B'_{0}$ is the magnitude of the field outside the sheet, $\delta'$ is a measure of the sheet thickness and the field is in the direction parallel to the sheet and perpendicular to the direction of current flow, which is locally considered to be $Z'$. The particle population in the sheet consists of two counter-drifting relativistic Maxwellian distributions whose density falls with $X'$, and whose drift provides the net current. From the pressure balance across the sheet follows the dependence of the density of each distribution on $X'$ (Kirk & Skjæraasen 2003):

$$N'_{\pm} = N'_{\pm 0} \cosh^{-2} \left( \frac{X'}{\delta'} \right) \qquad (4)$$

$$N'_{\pm 0} = \frac{B'_{0}^{2}}{16\pi mc^{2}\Theta}, \qquad (5)$$

where $\Theta = k_B T'/(mc^2)$ is a dimensionless temperature associated with the particle distribution. The condition connecting the particle density, temperature, and current sheet thickness is (Kirk & Skjæraasen 2003; Lyubarsky & Kirk 2001)

$$\frac{\Theta mc^2}{4\pi N'_{\pm 0} e^2 \gamma_{\pm}} = \delta'^2 \beta_{\pm}^2, \qquad (6)$$

where $\gamma_{\pm}$ and $\beta_{\pm}$ are the Lorentz factor and the corresponding speed, normalized to the speed of light, of the drift of the distributions in the wind frame. Finally, by assuming that the ideal gas law holds for the distribution, we obtain

$$p' = 2N'_{\pm}T' = (\hat{\gamma} - 1)(e' - 2N'_{\pm}mc^2), \qquad (7)$$

where $e'$ is the energy density associated with the two counter-streaming distributions. For the relativistic Harris sheet solution, this equation becomes

$$\hat{\gamma} = 1 + \frac{\Theta}{\gamma_{\pm}^2 (1 + \beta_{\pm}^2) \Theta + \gamma_{\pm} (K_1(1/\Theta)/K_2(1/\Theta) - 1)}, \qquad (8)$$

where $K_1$ and $K_2$ are modified Bessel functions of the second kind. In the inner part of the pulsar wind, which is of interest to us, the thermal particles in the distribution are predicted to be highly relativistic, which allows us to set $\hat{\gamma} = 4/3$ in the above equation and also has as a result that the thermal energy of the particles in the distributions will greatly exceed their rest-mass energy: $\Theta \gg 1$. Using the approximations

$$K_1(x \ll 1) \approx x^{-1} \qquad (9)$$

$$K_2(x \ll 1) \approx 2x^{-2}. \qquad (10)$$

Equation (8) can be simplified to

$$\beta_{\pm} \approx \frac{1}{\sqrt{\Theta}}. \qquad (11)$$
This means that the drift of the distributions in the current sheet is not relativistic, and beaming effects caused by it can be ignored when calculating the emitted radiation. Therefore the distributions can for this purpose be considered to be isotropic in the frame in which the sheet is at rest. With this approximation we can solve Eqs. (5) and (6) for $\Theta$ and $N'_{\pm 0}$. In doing this, we have to take into account the fact that the angle arccos($\hat{n} \cdot \hat{r}$) between the sheet normal and the radius unit vector changes when moving from the lab to the wind frame.

The magnetic field from the solution of Bogovalov (1999) is always parallel to the sheet, if one ignores reconnection effects. This means that close to the light cylinder and for a radial outflow, one cannot ignore the radial component of the field. The full field in the wind frame will be

$$B'_{0} = \frac{B_{LC}}{R^2} \sqrt{1 + \left(\frac{R \sin \vartheta}{\beta \Gamma}\right)^2} \quad (12)$$

and the toroidal component will prevail in this frame only for $R > \Gamma$ (we call this region the far wind region). Conversely, in the region $R < \Gamma$ (which we call the near wind region), it is a good approximation to ignore the second term under the square root and approximate $B'_{0} \approx B_{LC}/R^2$.

If one considers a perpendicular rotator (a pulsar with $\chi = \pi/2$), the temperature and density in a current sheet segment are given by the simple expressions

$$\Theta = \left( \frac{a_{LC} \Delta \sin \vartheta}{2R} \right)^{2/3} \quad (13)$$

$$N'_{\pm 0} = \frac{B'_{0}^{2}}{16\pi mc^{2}\Theta} \quad (14)$$

The dimensionless parameter $a_{LC}$ appearing above is called the strength parameter, and is defined at the pulsar's light cylinder as

$$a_{LC} = \frac{eB_{LC}P}{2\pi mc} \quad (15)$$

Taking a pulsar's moment of inertia to be $I = 10^{45}$ gcm$^2$, the spindown luminosity is calculated by the period and the period derivative as (Abdo et al. 2010)

$$\dot{E}_{33} = 4\pi^2 10^{12} \dot{P} P^{-3}, \quad (16)$$

where we have normalized the luminosity to the value $10^{33}$ erg/s. The magnetic field (in Gauss units) can then be expressed as

$$B_{LC} = 46.83 \dot{E}_{33}^{1/2} P^{-1}. \quad (17)$$

The strength parameter then can be expressed as a function of the spindown luminosity only:

$$a_{LC} = 1.3 \times 10^8 \dot{E}_{33}^{1/2} P^{-1}. \quad (18)$$

All GeV pulsars detected so far have strength parameters in the range $10^8-10^{11}$.

Finally, $\Delta < 1$ is the fraction of a half wavelength that a sheet occupies, measured in the radial direction (i.e. not perpendicularly to the local sheet plane), in the lab frame. We assume that $\Delta$ is constant with radius for the sake of simplicity, but it is far from clear how this parameter evolves with radius and obliquity $\chi$ close to the light cylinder. We also assume that $\Gamma$ is constant, although generally reconnection in the current sheet has been shown to increase $\Gamma$ with radius in the far wind zone (Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003), a result that might also

apply to the near wind zone. These assumptions are based on the fact that only a very limited radius interval contributes to the gamma-ray radiation, as we will see below, therefore if the change in $\Delta$ and $\Gamma$ happens with a scale larger than $r_{LC}$, it is not relevant to the present estimations.

Substituting for $a_{LC}$ in Eq. (13), we obtain

$$\Theta = 1.6 \times 10^5 \dot{E}_{33}^{1/3} \left( \frac{\Delta \sin \vartheta}{R} \right)^{2/3} \quad (19)$$

The requirement that $\Theta \gg 1$ in the near wind region $R < \Gamma$ translates into

$$\Delta \gg 10^{-8} \frac{\Gamma}{\dot{E}_{33}^{1/2} \sin \vartheta}, \quad (20)$$

which, as we will see in the following, is easily satisfied for all γ-ray pulsars, provided the emitting region is not very close to the polar axis. In any case, the validity of the $\Theta \gg 1$ assumption has to be checked a posteriori for all cases of studied pulsars.

In the following, all results refer to the near wind region, $R < \Gamma$, but still outside the pulsar's light cylinder, close to which we assume that the current sheet is formed. Therefore we consider a minimum value for the radius of the radiating sheet of $R_{min} = 1$. Our model is valid only for relativistic outflows, therefore we also set a lower limit on $\Gamma$ of $\Gamma_{min} = 10$.

In the next paragraphs, we assume a perpendicular rotator ($\chi = \pi/2$) to estimate analytically the emitted radiation and the characteristics of the current sheet and its populations. These estimates can be generalized to the oblique rotator (see Appendix B), but the conclusions remain essentially the same as long as the line of sight is not close to the edge of the current sheet, where $\zeta = \pi/2 - \chi$.

## 2.2. Larmor radius of hot particles

An important condition for the consistency of our model comes from the requirement that the hot particles' Larmor radius in the full field between the current sheets should be smaller than the sheet width. In the wind frame this requirement translates into

$$\frac{\langle \gamma \rangle mc^2}{eB'_{0}} \gg \Delta r_{LC} R \left( 1 + \frac{R^2}{\Gamma^2} \right)^{-1/2} \quad (21)$$

This inequality is equivalent to $\Theta^{1/2} \gg 1$, which should hold close to the light cylinder, by the assumptions of the model.

## 2.3. Energy of emitted photons

The mean Lorentz factor of the electron/positron distribution in the current sheet is $\langle \gamma \rangle \sim 3\Theta$. Electrons of this $\langle \gamma \rangle$ gyrating in the field $B'_{0}$, radiate photons of energy (in the observer's frame):

$$\varepsilon \approx \frac{3}{2} \Gamma (3\Theta)^2 \hbar \frac{eB'_{0}}{mc}. \quad (22)$$

After some manipulation, one finds the highest observable photon energy

$$\varepsilon_{\text{GeV}} \approx 2 \times 10^{-4} \frac{\dot{E}_{33}^{7/6}}{P} \frac{\Gamma (\Delta \sin \zeta)^{4/3}}{R_{\text{min}}^{10/3}}, \quad (23)$$

where the photon energy is measured in GeV. In Eq. (23) we have kept only the poloidal component of the magnetic field,
since it dominates close to the light cylinder, at $R < \Gamma$. All pulsars in the first Fermi/LAT catalogue have $\dot{E}_{33}^{7/6}/P \gg 1$. This is why it is possible for most objects in the catalogue to find an appropriate combination of $\Gamma$, $\Delta$ and $R_{min} \ge 1$ that will bring $\mathcal{E}_{GeV}$ to the regime observed by LAT, while at the same time satisfying all restrictions on $\Delta$, $\Gamma$, and $R$. This is especially easy for the highest luminosity pulsars ($\dot{E}_{33} \gg 1$), or for the very low period ones ($P \ll 1$ s, millisecond pulsars), as one can deduce by inspecting Eq. (23). In other words, for pulsars of the same spin-down luminosity, millisecond pulsars are more likely to reach GeV energies. Alternatively, if $\Gamma$ and $\Delta$ are the same, millisecond pulsars with lower spindown luminosities than young pulsars can emit in the Fermi/LAT band, a trend that might be observable as our statistics on gamma-ray pulsars increase.

From the constraints on $\Delta$ and $R$, we can calculate an absolute maximum on the peak energy of the emitted spectrum for $\zeta = \pi/2$:

$$ \mathcal{E}_{GeV,max} \approx 2.3 \times 10^{-5} \Gamma \frac{\dot{E}_{33}^{7/6}}{P}. \qquad (24) $$

As we will see below, $\Delta$ is often more severely constrained by the radiation reaction limit, so the energy $\mathcal{E}_{GeV,max}$ will not necessarily be reached for all objects.

The rapid fall of the peak photon energy with radius means that in this simple model the most energetic photons come from close to the light cylinder. In a realistic situation it is to be expected that either $\Gamma$ or $\Delta$ or both will rise with radius, making the fall of $\mathcal{E}_{GeV}$ with radius less steep, but even if $\Delta$ and $\Gamma$ were to rise linearly with $R$, the emitted photon energy would still fall as $\propto R^{-1}$.

## 2.4. Energy losses

For the relativistic Harris equilibrium to hold in a quasi-steady state, equilibrium has to be established in the wind frame within a timescale of the order of magnitude $t'_R \sim R(\Gamma\omega)^{-1}$, which is the timescale on which the magnetic field change is comparable to its magnitude $\delta B' \sim B'$. Equilibrium is established in the current sheet within a timescale

$$ t'_{eq} \le \frac{\Delta r_{LC}}{c}. $$

Therefore the requirement $t'_{eq} < t'_R$ holds close to the light cylinder as long as $\Delta \le \Gamma^{-1}$, while this condition is relaxed with radius to $\Delta \le R/\Gamma$.

For the particles not to suffer catastrophic energy losses, their synchrotron cooling timescale $t'_s$ should be longer than $t'_R$. Comparing the two one finds

$$ \frac{t'_s}{t'_R} \approx 4.57 \frac{P}{\dot{E}_{33}^{4/3}} \frac{\Gamma R^{11/3}}{(\Delta \sin \zeta)^{2/3}}. \qquad (25) $$

For the lower $\dot{E}_{33}$ pulsars the inequality $t'_s \ge t'_R$ generally holds for $R > R_{min}$, but this may not be the case for some millisecond pulsars, because of the linear dependency of the ratio on the pulsar period $P$. The steep radius dependence, however, relaxes this condition rapidly, so that when it breaks down, it does so only for a very limited $R$-range close to the light cylinder.

If $t'_s/t'_R \ll 1$, as is often the case for high spindown pulsars and millisecond pulsars, particles will lose energy rapidly and cool, disturbing the pressure equilibrium between the populations in the current sheet and the magnetic field outside it. The pressure in the sheet will fall, thus causing its compression

by the external magnetic field. This can initiate compression-driven magnetic reconnection, a phenomenon previously studied in the context of the interaction of current sheets in a pulsar wind with the wind's termination shock (Lyubarsky 2003; Pétri & Lyubarsky 2007; Lyubarsky & Liverts 2008; Sironi & Spitkovsky 2011). Reconnection has been shown to be able to accelerate particles to non-thermal distributions above the thermal peak of the particle spectrum (and also, as already mentioned, to accelerate the bulk flow and cause the sheet width to rise). A detailed study of this radiation-loss-induced process is beyond the scope of this work. But if $t'_s/t'_R < 1$, some acceleration processes ought to be at work to supply the sheet with energetic particles. This will in all likelihood result in power-law distributions in the current sheet, a signature that might be observable in the spectrum above the MeV-GeV peak for the highest luminosity pulsars and for those millisecond pulsars with $t'_s/t'_R \ll 1$ (something that has not been modelled in the present article). Uzdensov & Spitkovsky (2012) recently reported on their investigation of this regime.

## 2.5. Radiation reaction limit and maximum emitted energy

If the particles of mean energy $\sim 3\Theta mc^2$ in the current are gaining energy through acceleration in an electric field $E \sim \xi B$ with $\xi < 1$, their acceleration is limited by radiation losses. For radiation-reaction-limited synchrotron emission (in the full-field amplitude $B'_0$), energy losses compensate possible energy gain:

$$ \left[ \frac{dW}{dt} \right]_{\text{syn}} = e \xi B'_{0} c \le e B'_{0} c. \qquad (26) $$

From this expression one can calculate an upper limit for $\Delta$:

$$ \Delta \le \Delta_{\text{lim}} = 256 \frac{R^{5/2} P^{3/4}}{\dot{E}_{33}^{7/8} \sin \zeta}. \qquad (27) $$

The $\dot{E}_{33}$ dependence in the denominator implies that for the most powerful pulsars, the sheet is thinner close to the light cylinder. For the Crab we obtain $\Delta_{\text{lim}} = 3 \times 10^{-4} \sin^{-1} \zeta$ at $R_{min}$, whereas for most of the weaker pulsars with $\dot{E}_{33} \le 10^2$, $\Delta_{\text{lim}}$ is greater than unity, i.e. it does not present a constraint. All millisecond pulsars detected by Fermi/LAT fall into this category.

Inserting $\Delta_{\text{lim}}$ into Eq. (23), one recovers the classical synchrotron limit:

$$ \mathcal{E}_{GeV,max} < 0.36\Gamma, \qquad (28) $$

where the factor $\Gamma$ comes from the boosting of the photon energy to the lab frame. The minimum of Eqs. (24), (28), then, can give an estimate of the maximum value of the GeV cutoff for any individual object.

Inserting $\Delta_{\text{lim}}$ into Eq. (25), one obtains a lower limit for the ratio $t'_s/t'_R$ (applicable mainly to higher-spindown pulsars, for which $\Delta_{\text{lim}} < 1$):

$$ \frac{t'_s}{t'_R} > \left[ \frac{t'_s}{t'_R} \right]_{\text{min}} = 0.23 \Gamma R^2 P^{1/2} \dot{E}_{33}^{-3/4}. \qquad (29) $$

We see that, mainly for young pulsars with a high spindown luminosity, the right-hand side can be less than unity. These pulsars either have smaller $\Delta$ to keep $t'_s/t'_R > 1$, or rapid acceleration mechanisms in their current sheets, as argued above (or both). Therefore there should be a trend in the $P-\dot{P}$ diagram to observe more objects with acceleration signatures (i.e. power-law tails) in their spectrum as one moves to higher spindown luminosities, or equivalently as one moves to higher $\dot{P}$ and lower $P$ (upward left region in a $P-\dot{P}$ diagram).
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Sky maps of the pulsar's total luminosity shown for two different bulk Lorentz factors $\Gamma = 10$ and $\Gamma = 50$ and for different $\chi$. In the *upper panels* $\chi = \pi/2$, in the *middle ones* $\chi = \pi/3$ and in the *lower ones* $\chi = \pi/6$. All the plots have been made using the parameters $\Delta = 0.01$ and $R_{\min} = 1$. The luminosity is given by integrating $L_v$ (Eq. (A.4)) over all frequencies and multiplying with $d^2$, where $d$ is the distance of the pulsar from the observer. A beaming factor of $(4\pi)^{-1}$ has been used.</figcaption></figure>

# 3. Lightcurves, spectra, and γ-ray luminosity

## 3.1. Lightcurves

The expected lightcurves and spectra produced in the current sheets of a pulsar can be computed numerically by integrating the emission coefficient of the particles in the magnetic field of the current sheet along the line of sight to the observer. The procedure is explained in the appendix.

In the wind model of the pulsar radiation two pulses per period are expected, the separation and the width of which varies with the obliquity $\chi$ and the line of sight angle $\zeta$. The width of the pulses depends on the bulk speed of the outflow, with wider pulses for lower $\Gamma$ (Kirk et al. 2002; Pétri 2011). Examples of the variations of skymaps with $\chi, \zeta$ and $\Gamma$ are shown in Fig. 1, where we used a model pulsar of surface magnetic field equal to $B = 10^{10}$ G and period $P = 0.01$ s to plot the sky maps of the pulsed total luminosity.

In the left column of Fig. 1 sky maps are shown for $\Gamma = 10$, while in the right column the value of the bulk Lorentz factor is $\Gamma = 50$. The width of the pulse strongly depends on $\Gamma$, as expected. Along with the gamma-ray intensity, we have indicated the position of the radio pulses on the maps with white shading. The discs that represent the radio pulses are only indicative. They are centred on the spot where a sharp radio pulse coming directly from the corresponding pole would appear, and serve to demonstrate the phase lag between the radio and the gamma-ray pulse for different obliquities. As is seen from the sky maps, both pulses are observable only for a limited range of observer angle $\zeta$, and for inclinations higher than $\chi = \pi/4$. For smaller inclination angles, like in the case $\chi = \pi/6$, the beam from the pole does not travel through the pulsar wind, the radio pulses are above and below the region of high gamma-ray luminosity in the sky map, and only one of the two kinds can be observed. The exact shape and width of the radio pulse will depend on the
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Peak flux $vF_v$ (multiplied by the square of the distance to give a luminosity estimate) as a function of phase $\phi$ and peak energy $\mathcal{E}$. $\mathcal{E}$ ranges from 10 MeV to 1 TeV. The same pulsar parameters were used as in Fig. 1, with $\chi = \pi/3$ and $\zeta = \pi/2$.</figcaption></figure>

geometry of the radiating particle beam above the polar cap. The details of the position calculation of the radio pulse can be found in Pétri (2011).

The upper maps in both columns represent the case of the perpendicular rotator, where $\chi = \pi/2$, the middle ones are for $\chi = \pi/3$, and the lower ones for $\chi = \pi/6$. As the obliquity decreases, the width of the pulse increases while the peak luminosity falls, which causes the overall phase-averaged luminosity to remain essentially constant for different obliquities. The widest pulses are observed for the smallest inclination angles. However, for very small $\chi$ the line of sight to the observer has to lie very close to the equatorial plane of the pulsar for the pulse to be observed, since the luminosity is significant only for $\zeta \ge \pi/2 - \chi$. This is a likely geometry for the objects that are observed to have two wide pulses with a phase separation of $\delta\phi \sim 0.5$.

From the sky maps for $\Gamma = 10$ it can be discerned that there is a slight substructure in the pulses. There is a slight dip in luminosity at the peak of each pulse, resulting in two sub-peaks appearing for all $\zeta$. This is caused by the structure of the current sheet. In the middle of the sheet the magnetic field is exactly zero, rising towards the sides, while the density of the hot particles is highest and falls towards the sides. Therefore the bulk of the radiation of each current sheet segment comes from two regions away from the sheet midplane, where the product of particle density and magnetic field is largest. The double-peaked shape of the pulse reflects the two luminosity peaks within the sheet. Because of the shrinking pulse thickness with rising $\Gamma$ this phenomenon should be observable only for the lowest bulk Lorentz factors.

It is useful to note here that the highest luminosities in the sky maps are correlated with the highest peak energies $\mathcal{E}_{\text{GeV}}$, resulting in pulses that become sharper for higher energy photons. This is shown in Fig. 2, where the flux $vF_v$ has been plotted in a colour map as a function of emitted energy $\mathcal{E}$ (in eV) and phase $\phi$. A horizontal slice of the map shows the lightcurve at a specific energy while a vertical slice shows the spectrum at a given phase. The pulses become narrower with increasing energy, as can be seen from the wider flux variation in a horizontal slice as one moves higher in the map. Also, the peak energy varies by many orders of magnitude within one phase, and is lowest between pulses, something that can be seen when comparing spectra at different phases.

In Fig. 3 we show how the lightcurves change when moving the minimum radius $R_{\text{min}}$ from the light cylinder to a distance in

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Variation of the pulsar's lightcurves for different $R_{\text{min}}$. The bulk Lorentz factor in these examples is $\Gamma = 10$ and the obliquity is $\chi = \pi/3$. The rest of the parameters are the same as in Fig. 1.</figcaption></figure>

the far wind region, $R_{\text{min}} > \Gamma$. For low $R_{\text{min}}$, close to the light cylinder, the pulse shape is almost symmetrical with respect to its peak (also seen in Fig. 1). As one moves towards larger $R_{\text{min}}$, the pulses stop being symmetrical and instead present an abrupt rise in luminosity followed by a gradual decrease, which makes them asymmetrical with respect to the peak. Also, the total luminosity in the pulse is reduced by many orders of magnitude as one moves outwards in the wind. In the example of Fig. 3 a change from $R_{\text{min}} = 1$ to $R_{\text{min}} = 30$ results in a decrease in luminosity by seven orders of magnitude. These results agree with previous investigations of pulses from the far wind region (for example Kirk et al. 2002).

From the sky maps we showed we can conclude that the separation of the two pulses depends on $\chi$ and $\zeta$, with a separation of $\Delta\phi \sim 0.5$ when the line of sight lies at the equatorial plane. However, as we have pointed out above, the phase-averaged luminosity changes significantly only if the line of sight lies above
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Phase-averaged spectra for the same model pulsar that was used in Figs. 1 and 3. All the plots have been made using $\chi = \pi/2$ (perpendicular rotator) and $\zeta = \pi/2$. Here the variation of the phase-average flux with the parameters $\Gamma$, $\Delta$ and $R_{\min}$ is shown.</figcaption></figure>

the wind, i.e. if $\zeta < \pi/2 - \chi$. In the opposite case the overall luminosity is not sensitive to variations of $\chi$ and $\zeta$.

## 3.2. Phase-averaged spectra

In Fig. 4 we used the same model pulsar as in the previous section to calculate phase-averaged spectra and their dependence on the model parameters. In the upper plot we show the variation of the spectrum when changing the bulk Lorentz factor $\Gamma$. As is expected from Eq. (23), the peak of the spectrum rises linearly with $\Gamma$, but the overall luminosity falls with rising $\Gamma$. This can be attributed to the diminishing of the area of the sheet from which doppler-boosted radiation is received: the boosted area scales as $1/\Gamma^2$, something that is only partly compensated for by the

square of the Doppler factor in the calculation of the received flux (given in the appendix). The most strongly boosted radiation comes from the line of sight, but the luminosity of the region that directly intersects the line of sight is not the major contribution to the overall observed luminosity. This is because one observes only the effects of the perpendicular (i.e. azimuthal) field at $\theta = \zeta$, while the effects of the much stronger poloidal field $B_r$ come from a region at angle $\sim 1/\Gamma$ to the line of sight. Therefore when the line of sight is directly aligned to the edge of the sheet, where $\zeta = \pi/2 - \chi$, still two pulses of the same amplitude are observed, that come from the two parts of the folded current sheet that are at angle $\sim 1/\Gamma$ to the line of sight, whereas the edge itself contributes only a little to the overall flux. This is different to what has been predicted in the past (Pétri 2011). However, a single wide pulse will be observed when looking over the edge of the sheet. The disadvantage in this case is that the luminosity decreases quickly with decreasing $\zeta$, which makes such pulses difficult to detect.

In the middle panel the change in the phase-averaged spectrum is shown when varying $\Delta$. The dramatic rise of the cutoff energy and of the luminosity with $\Delta$ can be explained by the dependency of the peak of the spectrum on $\Delta$, given in Eq. (23), as well as the rise of the relativistic temperature in the current sheet with $\Delta$ as seen in Eq. (19).

In the lower panel we demonstrate the dependence of the cutoff and luminosity on the minimum radius $R_{\min}$. The peak frequency of the spectrum and the overall luminosity both rise very strongly with decreasing $R_{\min}$, which agrees with Eq. (23) and is also expected because of the dependence of the temperature of the sheet particles on radius, given in Eq. (19). Therefore the bulk of the received high-energy radiation comes from a region of limited range in $R$ close to $R_{\min}$.

We would like to note at this point that the spectra that come from thermal distributions have a characteristic slope in the lower frequencies that is the same as the single-particle synchrotron spectrum (with a flux of $F_v \propto v^{1/3}$). This slope is not clearly seen in the 100 MeV–10 GeV range for most pulsars, because there is a broad peak in that region. This peak is created mainly by the spread of the magnetic field across the current sheet, which causes different regions in the current sheet to emit photons of different peak energies. Because the emissivity decreases rapidly with radius, the phase-averaged spectrum mainly consists of the radiation emitted by the small part of the current sheet that moves towards the observer, for a very short radius interval, and is dominated by the spectrum at the peak of the pulse. Different spectral slopes could be obtained by assuming the injection of a power-law rather than a thermal distribution, albeit at the price of introducing additional parameters. We return to the question of the particle distribution in Sect. 4.

## 3.3. Examples

In Fig. 5 we present two examples of spectra and lightcurves of Fermi/LAT-detected pulsars calculated with our model. The first example is the millisecond pulsar PSR J1614-2230, which has a period of $P = 3.2$ ms and a spindown luminosity of $\dot{E}_{33} = 5$. Its spectrum calculated according to our model can be seen in the upper left panel of Fig. 5 along with the best fit power law with exponential cutoff, as given in Abdo et al. (2009). The parameters used were $\Delta = 0.25$, $\Gamma = 20$, $\chi = \pi/2$, and $\zeta = \pi/2$. The lightcurve corresponding to the same parameters is depicted in the lower left panel. The phase of the two peaks is correctly predicted, but the width of the simulated pulses is narrower than observed. This is because in this example $\Delta > 1/\Gamma$, and in this
<figure><img src="image_6.png" /><figcaption>Fig. 5. Phase-averaged spectrum and simulated lightcurve of the millisecond pulsar J1614-2230 as calculated using our model, along with the best-fit power-law plus exponential cutoff fit (Abdo et al. 2009). The parameters that were used are $\Delta = 0.25$, $\Gamma = 20$, $\chi = \zeta = \pi/2$ for the millisecond pulsar (left). For the young pulsar in the supernova remnant 3C 58 we used $\Delta = \Delta_{lim}$, $\Gamma = 80$, $\chi = \pi/5$ and $\zeta = \pi/2$.</figcaption></figure>

case the lightcurve cannot be accurately predicted by our model (as explained in Appendix B).

The second example, shown in the plots on the right side of Fig. 5, is the pulsar PSR J0205+6449, which was discovered by Fermi/LAT in the supernova remnant 3C 58. This pulsar has a high spindown luminosity of $\dot{E}_{33} = 2.7 \times 10^4$ and a period of $P = 65.7$ ms$^{-1}$. It is a young, energetic pulsar for which $\Delta_{lim} \approx 1.5 \times 10^{-3}$. The parameters that we used for these plots are $\Delta = \Delta_{lim}$, $\Gamma = 80$, $\chi = \pi/5$, and $\zeta = \pi/2$. The relatively high Lorentz factor needed to reach the cutoff of $\mathcal{E} \approx 3$ GeV causes the pulses to be rather sharp, which is also observed with Fermi/LAT.

From these two examples and from the phase-averaged spectra shown in Fig. 4 it can be seen that our model predicts a substantially higher energy flux at photon energies above 10 GeV than would be expected using a power-law plus exponential cutoff fit, which implies that for the more energetic pulsars, which have larger bulk Lorentz factors, one could detect radiation in the sub-TeV to TeV regime explored by Cerenkov telescope arrays. Specifically, a good fit to our spectra is given by a power-law with a sub-exponential cutoff:

$$ \frac{dN}{d\mathcal{E}} \propto \mathcal{E}^{-p} \exp \left[ - \left( \frac{\mathcal{E}}{\mathcal{E}_{\text{cutoff}}} \right)^b \right]. \quad (30) $$

Our model predicts a cutoff with $b = 0.35$. This index is not sensitive to the obliquity $\chi$ or the observer angle $\zeta$, as long as the line of sight intersects the equatorial current sheet, i.e. when $\zeta > \pi/2 - \chi$. The index $b$ can reach the value $\sim 0.4$ for $\zeta < \pi/2 - \chi$,

but the very low luminosities that are expected in that case make the majority of those objects unobservable.

### 3.4. Peak energies and $\gamma$-ray luminosities

In Fig. 6 we show a prediction of the region in the $P-\dot{P}$ where spectra that peak in the range 0.1–10 GeV are expected, according to our model. Two shaded regions are shown that correspond to $\Delta = 0.1$ and $\Gamma = 100$ (yellow) and $\Delta = 0.001$ and $\Gamma = 100$ (green). We also plotted the lower and upper borders of the region that corresponds to $\Delta = 0.1$ and $\Gamma = 10$ (red lines). Obviously almost all pulsars already detected by Fermi/LAT fall in at least one of these regions, from which we conclude that one can accommodate the cutoffs of almost all pulsars with parameters in the range $0.001 \le \Delta \le 0.1$ and $10 < \Gamma < 100$.

In principle, the peak energies emitted by the near wind region of a pulsar wind are not constrained to the GeV regime. Peaks at TeV energies are likely excluded by the synchrotron limit given in Eq. (28), but there should be a population of pulsars with spectra peaking in the MeV regime (or lower). These pulsars can be found in a region of the $P-\dot{P}$ diagram roughly below the yellow shaded region of the $P-\dot{P}$ diagram in Fig. 6. The population of MeV pulsars will include objects with spindown luminosities lower than those detected by Fermi/LAT. This becomes obvious by inspecting Eq. (23): for lower spindown luminosities the same $\Gamma$ and $\Delta$ lead to lower peak energies. For example, a pulsar that is a perpendicular rotator with $P = 0.5$ s and $\dot{P} = 10^{-15}$ with a bulk Lorentz factor $\Gamma = 10$ and $\Delta = 0.1$ will radiate a spectrum that peaks at energies $\mathcal{E} \approx 0.3$ MeV.
<figure><img src="image_7.png" /><figcaption>Fig. 6. P − ˙ P diagram including the pulsars observed by Fermi / LAT and plot of the predicted γ -ray luminosity in the range 0 . 1 − 10 GeV for a random sample of pulsars. Gamma-ray luminosity L$_{γ}$ versus the spindown luminosity L$_{sp}$ for two random samples of millisecond and young pulsars. The lines L$_{γ}$ = L$_{sp}$ and L$_{γ}$ ∝ L 1 / 2 sp are also shown.</figcaption></figure>

The overall luminosity emitted by the near region of a pulsar wind from a perpendicular rotator according to our model can be estimated to within a factor of two by the expression

$$
\mathcal{L}_{\text{tot}} \approx \frac{2.7 \times 10^{31}}{\Gamma} \left(\frac{B_{\text{LC}}}{10^3}\right)^{4.7} \left(\frac{P}{0.1}\right)^{3.77} \left(\frac{\Delta}{0.1}\right)^{1.7} \text{ erg s}^{-1}, \quad (31)
$$

where B$_{LC}$ is given in G and P in s. In terms of spindown lumi- nosity, this can be expressed as

$$
\mathcal{L}_{\text{tot}} = \frac{7.6 \times 10^{29}}{\Gamma} \left(\frac{P}{0.1}\right)^{-0.9} \dot{E}_{33}^{2.35} \left(\frac{\Delta}{0.1}\right)^{1.7} \text{ erg s}^{-1}. \quad (32)
$$

It is obvious from the above expression that, since no trend is observed for γ -ray pulsars to have a γ -ray luminosity that is anti- correlated to their spindown luminosity, the values of Γ and Δ should vary significantly from object to object.

If the peak of the spectrum falls into the γ -ray regime, then Eq. (31) gives a good estimate of the γ -ray luminosity of the pul- sar. The dependence of L$_{tot}$ on χ and ζ is very weak as long as the line of sight intersects the wind, ζ > π/ 2 − χ , therefore one can use the above formula to reach approximate conclusions about the γ -ray luminosity of a pulsar. Equations (31) and (23) (which yield an estimate of the peak energy of the phase-averaged spec- trum) show that once P and ˙ P are known, the emitted luminosity and peak energy can be fixed essentially by the two parameters Δ and Γ. These two equations can also be inverted to deduce values of Δ and Γ for a sample of pulsars for which P , ˙ P , L$_{γ}$ and E$_{GeV}$ are known. This was done for the pulsars of the First Fermi / LAT pulsar catalogue and the results are shown in Fig. 7, where Δ and Γ have been plotted for each object of the sample.

The blue squares represent the young pulsars in the catalogue while the red triangles are the millisecond pulsars. For most ob- jects Δ ranges between 10 − 3 and 1, and a trend is observed for the millisecond pulsars to have higher values of Δ, slightly lower than unity. The bulk Lorentz factor Γ ranges roughly between 10 and 100, and most objects have relatively low Lorentz factors in the range 10 − 30.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Parameters Δ and Γ estimated for the pulsars in the Fermi / LAT one-year catalogue. Blue squares are the young pulsars and red triangles represent the millisecond pulsars.</figcaption></figure>

For the most part, the deduced Γ and Δ are reasonable, ex- cept for five objects for which values Δ > 1 are predicted, which are unphysical since Δ has to be less than unity for the current sheet to be narrower than the wind’s half wavelength. These are: J0633+1746 (the Geminga pulsar), J1057-5226, J1741-2054, J1836+5925, and J2021+4026. These five pulsars are distinct from the rest of the population, in that their estimated γ -ray lu- minosity is a large part of (and for the last two pulsars even larger than) their spindown luminosity. From Eq. (31) we see that a higher Δ is needed to accommodate a higher L$_{γ}$ , which explains the values deduced for these pulsars. The unphysical values of Δ, therefore, may arise because the γ -ray luminosity is estimated
from the observational data using too simple an assumption (a beaming factor of value unity). Also, our model is possibly too simple to accommodate all pulsars in the catalogue, and a contribution from another physical mechanism or emission region might result in higher luminosities in the exceptional cases of the objects with $\Delta > 1$.

In the diagram on the right side of Fig. 6 we have plotted the $\gamma$-ray luminosity $L_{\gamma}$ of two random samples of pulsars (millisecond and young pulsars) with respect to their spin-down luminosity $L_{\text{sp}}$. For this purpose two random samples of pulsars were used, with periods and period derivatives in the ranges $-2 < \log P < 0$ and $-15 < \log \dot{P} < -11$ (for the young pulsars) and $-3 < \log P < -2$ and $-20 < \log \dot{P} < -17$ (for the millisecond pulsars). We chose $\zeta$ and $\chi$ randomly in the interval $[0, \pi/2]$, the bulk Lorentz factor in the interval $1 < \log \Gamma < 2$, and the sheet thickness either within $\log \Delta_{\text{lim}} - 2 < \log \Delta < \log \Delta_{\text{lim}}$ in the case $\Delta_{\text{lim}} < 0.1$ or within $-3 < \log \Delta < -1$ if $\Delta_{\text{lim}} > 0.1$. The $\gamma$-ray luminosity was calculated by integrating the phase-averaged flux between 100 MeV and 100 GeV and using a beam correction factor $F_{\Omega} = 1$, to obtain a direct comparison with the results of Abdo et al. (2010).

About 12% of all young pulsars in our sample and 30% of the millisecond pulsars have a cutoff within the interval $100 \text{ MeV} - 10 \text{ GeV}$ and are shown in Fig. 6. These percentages reflect our choice of the range in $\Delta$. We constrained $\Delta$ to be within two orders of magnitude of the value min($\Delta_{\text{min}}$, 0.1). If the range of $\Delta$ were narrower, resulting in systematically higher values, the percentages quoted above would rise. Lower choices of $\Delta$ lead to lower $E_{\text{GeV}}$ so there are more pulsars that have peaks below 100 MeV and therefore do not make the cut.

The large beam correction factor used for this plot has as a result the occurrence in some rare cases of objects with $L_{\gamma} > L_{\text{sp}}$, whereas in most cases the $\gamma$-ray luminosity is a few percent of $L_{\text{sp}}$. The young pulsars tend to lie higher in the diagram than millisecond pulsars, which reflects their higher spindown luminosity in combination with the fact that the gamma-ray efficiency is similar in both samples.

In Abdo et al. (2010) a trend was reported that the gamma-ray luminosity is proportional to the square root of the spindown luminosity for the highest luminosities: $L_{\gamma} \propto L_{\text{sp}}^{1/2}$. We have no clear indication of such an effect in our simulated sample, meaning that the gamma-ray efficiency of the highest spindown pulsars is predicted to be similar to that of the lower spindown pulsars. The best-fit line to the points in Fig. 6 gives $L_{\gamma} \propto L_{\text{sp}}^{0.95}$ for young and millisecond pulsars, and if the statistical error is taken into account, the relationship is compatible with the proportionality $L_{\gamma} \propto L_{\text{sp}}$.

## 4. Discussion and conclusions

We have shown how gamma-ray pulses can naturally arise within the framework of the pulsar wind's equatorial current sheet outside (but close to) the light cylinder. The advantage of this emission model is that it is an intrinsic mechanism that naturally produces peak energies in the MeV-GeV range. It can give meaningful results for almost all pulsars, irrespective of age or environment, employing only a few parameters: the bulk Lorentz factor of the outflow $\Gamma$, the sheet thickness $\Delta$, the obliquity $\chi$, and the angle of the rotational axis to the line of sight $\zeta$. To these the minimum radius $R_{\text{min}}$ can be added (which, in the calculated examples, was set equal to $R_{\text{min}} = 1$). However, it is important to note that $R_{\text{min}}$ has to be sufficiently close to the light cylinder, so that the peak energy of the spectrum reaches the GeV regime. The parameters of the model are constrained

by the characteristics of the pulsar and by the physics of the current sheet. The novelty of our model is that it makes an attempt to account for the previously ignored region of the wind between $R=1$ and $R \approx \Gamma$, and is therefore able to take advantage of the strong poloidal field close to the light cylinder, something that has not been discussed in previous wind models.

The predictions of our model are the following:

1. Two pulses per pulsar period are expected. The pulses have the same amplitude.

2. The shape of the pulses is symmetric with respect to their peak, as long as emission starts close to the light cylinder, and becomes increasingly asymmetric with rising $R_{\text{min}}$.

3. There is no significant interpulse emission in the GeV range.

4. The width of the pulse decreases with increasing photon energy.

5. The separation of the peaks varies with the obliquity $\chi$ and the angle to the line of sight $\zeta$, as in previous models of pulsed radiation from the pulsar wind.

The features of the pulsed emission that cannot be explained by our model, such as the interpulse or the peaks of different intensity, could be accomodated by a model in which radiation comes both from within the light cylinder, possibly in an outer gap, and from the beginning of the wind. The difficulty of outer gap models to produce double-peaked lightcurves, as noted in Bai & Spitkovsky (2010), in combination with the difficulty of our near wind model to produce single pulses, point to the need for a combined model in which the high-energy emission starts within the light cylinder and continues in the equatorial current sheet outside the light cylinder. Therefore the physics of the formation region of the current sheet is very important in the endeavour to understand $\gamma$-ray emission from pulsars.

For pulsars in binaries, it is possible that the favoured radiation process is not synchrotron but rather inverse Compton on the low-frequency photons of the massive star, as investigated in Pétri & Dubus (2011). In this case the peak energies and emitted luminosities by the two mechanisms can be comparable, depending on whether the comoving energy density of the low-frequency photons in the wind frame is comparable to the magnetic field energy density close to the light cylinder. For the pulsar B1259-63, which is a member of a binary with a B2e star, the comoving stellar radiation energy density seen by the pairs moving towards the star is $u'_{\star} \approx 5 \Gamma^2 \text{ erg s}^{-1}$, while the comoving magnetic field energy density is $u'_{\text{B}} \approx 3 \times 10^7 R^{-4} \text{ erg s}^{-1}$ (the radial component dominates when $R < \Gamma$). Inverse Compton scattering is dominant if $\Gamma R^2 / \gtrsim 2500$, which requires emission to start at $5-15$ light cylinder radii for values of $\Gamma \approx 10-100$. The situation is more favourable to Compton scattering in the other, more compact, gamma-ray binaries where the radiation density is up to 200 times higher (LS 5039).

A possible problem of our model is that the supersonic solution of Bogovalov (1999) might not apply for the objects of lower bulk Lorentz factors. Particularly, the wind can accelerate, as is predicted in various magnetohydrodynamic models of pulsar winds (Beskin et al. 1998; Kirk et al. 2009). In this case the shape of the current sheet can in reality be different, since there is generally also a polar magnetic field component, which in our model was zero. However, as was mentioned above, the very strong dependence of $E_{\text{GeV}}$ on $R$ implies that even if the parameters $\Delta$ and $\Gamma$ changed with $R$, the present model would still be able to explain the main features of the observed spectra. Adding to these considerations the fact that the particle density within the current sheet and the magnetic field decrease as the wind expands to larger $R$, the conclusion is reached that the
emitted luminosity, which depends on $N'_0$ and $B'_0$, will also decrease with radius. This in combination with the decreasing peak of the spectrum means that the main contribution to the spectrum near the cutoff will still come from a very limited region close to $R_{min}$. Therefore, possible changes in the dynamic evolution of the wind should not have a strong effect on the observed spectrum and luminosity, while the exact shape of the current sheet should influence mainly the shape of the pulse, but not the cutoff energy or the emitted luminosity.

Another possible shortcoming of our model lies in the assumption of a thermal distribution in the current sheet. The Harris distribution has been extensively used in the past in the context of the current sheets in pulsar winds (Kirk & Skjæraasen 2003; Lyubarsky & Kirk 2001; Sironi & Spitkovsky 2011), which makes it the first reasonable assumption when one treats these sheets. It is generally thought, however, that the particle distribution close to the light cylinder should be non-thermal, since the particles are accelerated by electric fields within the light cylinder and possibly by the reconnection process at the Y-point, where the current sheet originates. Nevertheless, using the relativistic Harris equilibrium entails several advantages: it is the simplest self-consistent solution to a relativistic current sheet, and it has only two free parameters: a normalization constant, and the particle temperature. Because of these restrictions there are no degeneracies in our model, as we showed. The observational data uniquely predict the properties of the current sheet for a given pulsar. This makes our model powerful and able to give unambiguous predictions, but also refutable, should it not accomodate the data. In the light of these advantages, we kept the thermal distribution for our current sheets, and dit not resort to a power-law distribution, which would pose fewer restrictions on our model because of the additional free parameters associated with it. In principle, one way to distinguish between power-laws and thermal inputs is the low-energy spectral slope. For thermal distributions the slope asymptotes to 1/3, whereas a wider range is expected for power-laws depending on the distribution index.

It is expected that reconnection will generally alter the physics of the current sheet beyond the pulsar's light cylinder, and might result in an evolution of the current sheet thickness and the wind Lorentz factor, especially in objects for which the radiative timescales of the thermal particles are very short in comparison to the evolution timescale of the magnetic field (Uzdensky & Spitkovsky 2012). Furthermore, particle acceleration takes place during the reconnection process, which might result in a non-thermal tail to the thermal distribution, the characteristics of which will change as the current sheet evolves with radius. These phenomena will leave their imprint on the pulsed high-energy spectrum, and need to be investigated in a self-consistent way, which is beyond the scope of this article. We can, however, give a simple order-of-magnitude argument about the reconnection-accelerated particles: deep in the current sheet particles can be accelerated by reconnection electric fields to energies higher than the thermal peak. These particles would have Lorentz factors extending to $\gamma \sim \xi a_{LC}/\Gamma$ (in the wind frame), thus giving rise to photons up to an energy

$$ \mathcal{E}_{\text{GeV}} \sim 8.8 \frac{\xi^2 \dot{E}_{33}^{3/2}}{\Gamma P}, \qquad (33) $$

where $\xi = E'/B'_0 < 1$. This could extend beyond the thermal peak for sufficiently large $\xi$ or $\dot{E}_{33}$, whereas millisecond pulsars are again favoured by the dependence on the inverse of the period. If these high-energy non-thermal particles escape the acceleration site and radiate in the field within the current sheet,

their emission could give rise to a power-law tail extending beyond the GeV cutoff (Zenitani & Hoshino 2008). This mechanism should be prominent for objects for which $t'_s/t'_R < 1$, which tend to cluster at the upper left part of the $P-\dot{P}$ diagram, as discussed. It is for these pulsars that acceleration by reconnection could become prominent and be observed in the form of power-law tails in the GeV–TeV regime.

Another possibility is that the particles in the current sheet are already accelerated to a power-law distribution when the sheet starts radiating, in which case a different distribution function would have to be used to calculate the current sheet parameters (Balikhin & Gedalin 2008). This might apply particularly to millisecond pulsars, which have lower surface magnetic fields that lead to lower pair production rates and therefore less dense plasmas in their magnetospheres. In this case a non-thermal particle distribution seems likely to describe the physics of the current sheet in a more consistent way. We defer the investigation of such particle distributions to a future article.

Acknowledgements. We acknowledge support from the EC via contract ERC-StG-200911 and CNES. I.A. would like to thank Jérôme Pétri, Gilles Henri, and Anatoly Spitkovsky for useful discussions.

## Appendix A: Calculation of a received synchrotron spectrum by a relativistically moving sheet

The flux that an observer receives from a relativistically moving source is calculated by the formula (Lind & Blandford 1985):

$$ F_v = \frac{1}{d^2} \int \mathcal{D}^2 j' (\nu') dV. \qquad (A.1) $$

- Primed quantities are in the wind frame, i.e. the frame where the outflow is at rest. Since the wind is assumed to be flowing radially, at each azimuth and polar angle $(\varphi, \theta)$ there is a different, local, radially moving wind frame. Quantities pertaining to a point of the wind with coordinates $(r, \vartheta, \varphi)$ are calculated in this local wind frame.

- When calculating the emission coefficient $j'$ we should take care to calculate it in the direction of the line of sight, taking into account the aberration of photons in the wind frame (we return to this matter below).

- The volume element is taken in the observer's frame (or the frame in which the pulsar is at rest) and is $dV = r^2 \sin \vartheta dr d\vartheta d\varphi$.

- $\mathcal{D}$ is the doppler factor, which depends on the angle between the line of sight and the direction of motion of the wind. It is calculated by the expression

$$ \mathcal{D} = \frac{1}{\Gamma(1 - \beta \hat{O} \cdot \hat{r})} \qquad (A.2) $$

The unit vector $\hat{r}$ is the radial unit vector in a spherical coordinate system, whose centre is located at the pulsar. $\hat{O}$ is a unit vector in the direction of the observer. The cosine of the angle between these two vectors depends on $\vartheta$ and $\varphi$, and therefore $\mathcal{D} = \mathcal{D}(\vartheta, \varphi)$. The observed frequency is $\nu = \mathcal{D}\nu'$.

- We assume that the direction to the line of sight is at polar angle $\zeta$ and at $\varphi = 0$. In a cartesian system of coordinates, then, we have

$$ \hat{O} = \sin \zeta \hat{x} + \cos \zeta \hat{z}, $$

whereas the radial unit vector can be expressed as:

$$ \hat{r} = \sin \vartheta \cos \varphi \hat{x} + \sin \vartheta \sin \varphi \hat{y} + \cos \vartheta \hat{z}, $$
which gives us

$$ \hat{O} \cdot \hat{r} = \sin \zeta \sin \vartheta \cos \varphi + \cos \zeta \cos \vartheta. \quad (A.3) $$

- Finally, $d$ is the distance of the object from the observer. We define a luminosity associated with the pulsar as

$$ \mathcal{L}_v = 4\pi d^2 \mathcal{F}_{\Omega} F_v = 4\pi \mathcal{F}_{\Omega} \int \mathcal{D}^2 j' (v') \, d^3 x. \quad (A.4) $$

The parameter $\mathcal{F}_{\Omega}$ is the beam correction factor, which is taken to be equal to $\mathcal{F}_{\Omega} = (4\pi)^{-1}$ in the plots of the pulsar lightcurves and spectra, Figs. 1, 3, and 4. It is, however, taken to be equal to unity $\mathcal{F}_{\Omega} = 1$ for the $L_{\gamma} - L_{\text{sp}}$ plot in Fig. 6.

## Appendix B: Calculation of the emission coefficient

The emission coefficient can be calculated by multiplying the energy distribution of the particles by the single-electron synchrotron spectrum. This implies that the electrons are relativistic, which is, strictly speaking, not true for the whole population, since the distribution starts from $\gamma' = 1$. However, since the temperature of the distribution is relativistic, this approximation will not introduce a significant error. The distribution gives the number of electrons per unit gamma factor, per solid angle, and per unit volume:

$$ \frac{dN'}{d^3x'd\gamma'd\Omega'} = \frac{N'_0}{4\pi\Theta K_2(1/\Theta)} \gamma' \sqrt{\gamma'^2 - 1} e^{-\gamma'/\Theta} \cosh^{-2}\left(\frac{X'}{\Delta_{X'}}\right), $$

where $\Theta \gg 1$ is the dimensionless temperature of the distribution in units of $mc^2$ as defined in the text. This distribution is isotropic and its density decreases towards the edge of the sheet. The direction perpendicular to the sheet midplane is denoted by $X'$ and $\Delta_{X'}$ is the width of the sheet at its local rest frame.

Using the small argument approximation of the modified Bessel function $K_2$, Eq. (10), the distribution function can be expressed as

$$ \frac{dN'}{dV'd\gamma'd\Omega'} = \frac{N'_0}{2\Theta^3} \gamma' \sqrt{\gamma'^2 - 1} e^{-\gamma'/\Theta} \cosh^2\left(\frac{X'}{\Delta_{X'}}\right). \quad (B.1) $$

This then has to be doubled to take account of both species in the sheet (electrons and positrons). The single-electron synchrotron spectrum is

$$ \frac{dE'}{dt'dv'} = \frac{\sqrt{3}e^3 B'_{\perp}}{mc^2} F\left(\frac{v/D}{v'_{cr}}\right) \tanh\left(\frac{X'}{\Delta_{X'}}\right). $$

where $F(x)$ is the well-known synchrotron function. In the above expression, we encounter the following:

- The field $B'_{\perp} = B' \sin \alpha$ perpendicular to the direction that the photons that reach the observer travel along in the wind frame. $\alpha$ is the angle between the field in the wind frame and the photon trajectory, which can be computed through the scalar product $\vec{B}' \cdot \hat{O}'$. To calculate this, we need to find the relativistic aberration of photons that follow the line of sight in the lab frame. This is equivalent to calculating the aberration of the vector $\hat{O}$. The field has only a radial and an azimuthal component, so only the aberration of these two components of the unit vector $\hat{O}$ are needed. These components are

$$ O_r = \sin \zeta \sin \vartheta \cos \varphi + \cos \zeta \cos \vartheta \quad (B.2) $$

$$ O_{\varphi} = -\sin \zeta \sin \varphi. \quad (B.3) $$

The corresponding primed components in the local wind frame are

$$ O'_r = \frac{O_r - \beta}{1 - \beta O_r} \quad (B.4) $$

$$ O'_{\varphi} = \frac{O_{\varphi}}{\Gamma(1 - \beta O_r)} \quad (B.5) $$

and the angle $\alpha$ is

$$ \alpha = \arccos \left( \hat{O}' \cdot \vec{B}' / B' \right) = \arccos \left( \frac{O'_r B'_r + O'_{\varphi} B'_{\varphi}}{\sqrt{(B'_r)^2 + (B'_{\varphi})^2}} \right), \quad (B.6) $$

while the local magnetic field in the wind frame is

$$ B'_{r} = \frac{B_{LC}}{R^2} \quad (B.7) $$

$$ B'_{\varphi} = \frac{B_{LC}}{\beta \Gamma R} \sin \vartheta, \quad (B.8) $$

where the superscript LC refers to quantities at the light cylinder.

- The ratio $X'/\Delta_{X'}$ can be approximately converted to quantities in the lab frame to give:

$$ \frac{X'}{\Delta_{X'}} \approx \frac{R - R_0(\chi, \vartheta, \varphi, t)}{\Delta} \quad (B.9) $$

where $R_0(\chi, \vartheta, \varphi, t)$ is the location of the current sheet point under consideration at time $t$.

- The critical frequency $v'_{cr}$ is

$$ v'_{cr} = \frac{3eB'_{\perp}}{4\pi mc} \gamma'^2 \tanh\left(\frac{X'}{\Delta_{X'}}\right). \quad (B.10) $$

- In the general case of the oblique rotator the temperature is given by the expression

$$ \Theta = \Theta_{\perp} \left( \frac{1 + \left( \frac{\beta \Gamma}{R_0 \sin \vartheta} \right)^2}{1 + \left( \frac{\beta \Gamma}{R_0 \sin \vartheta} \right)^2 C(\vartheta, \chi)} \right)^{1/3} \quad (B.11) $$

$$ C(\vartheta, \chi) = 1 + \left( \frac{\csc \vartheta \cot \chi}{\sqrt{1 - \cot^2 \chi \cot^2 \vartheta}} \right)^2, \quad (B.12) $$

where $\Theta_{\perp}$ is the value of the temperature for the perpendicular rotator, given in Eq. (19).

The function $C(\vartheta, \chi)$ is very close to unity for a wide range of $\vartheta$ given a $\chi$, and deviates significantly from that value only for values $\vartheta \approx \pi/2 - \chi$, or for low obliquities, as is seen in Fig. B.1. Therefore, the perpendicular rotator is a good approximation to the general problem, and this is why we have used $\Theta_{\perp}$ in all our estimates.

The emission coefficient $j'$ has to be calculated at the retarded time

$$ t_{\text{ret}} \approx t - \frac{d}{c} + \frac{\vec{r} \cdot \hat{O}}{c}, $$

where $\vec{r}$ is the position vector of the radiating point in the current sheet, $d$ is the distance of the pulsar from the observer, and $t$ is the time that the observer measures. Normalizing time to $1/\omega$ and distance to $r_{LC}$ the retarded time becomes

$$ t_{\text{ret}}^* \approx t^* - \frac{d}{r_{\text{LC}}} + R_0 \hat{r} \cdot \hat{O} \quad (B.13) $$
<figure><img src="image_9.png" /><figcaption>**Fig. B.1.** Function C(ϑ, χ) plotted for 0 ≤ χ ≤ π/2, 0 ≤ ϑ ≤ π/2. Away from the sheet’s edge, where ϑ = π/2 − χ, its value is close to unity.</figcaption></figure>

(normalized time is denoted with an asterisk). The retarded time enters the computation via the equation for the motion of a current sheet:

$$R_0 = \beta (\pm \arccos(-\cot\chi \cot\vartheta) + t_{\text{ret}}^* - \varphi).$$

We replace the retarded time in the last equation from (B.13), and move the initial condition for time at $t^* = d/r_{LC}$ so that the new observer’s time is $t_d^* = t^* - d/r_{LC}$. We obtain an expression for the radius of the current sheet as a function of observer’s time, polar angle and azimuthal angle:

$$R_0 = \frac{\beta (\pm \arccos(-\cot\chi \cot\vartheta) + t_d^* - \varphi + 2k\pi)}{1 - \beta \hat{O} \cdot \hat{r}}. \quad (B.14)$$

Depending on how large or small the dot product $\hat{O} \cdot \hat{r}$ is, the radius $R_0$ where the current sheet part under consideration is located, the radiation from which is received at time $t_d^*$, can vary strongly. This means that at each moment $t_d^*$ the observer receives radiation that comes from different parts of the current sheet that propagates in the wind, with different coordinates $(R_0, \vartheta, \varphi)$. We note here that $R_0$ is considered to be the radius at the midplane of the current sheet.

Finally, the flux can be computed, according to the above, as follows:

$$\mathcal{F}_{\nu} = A \int \mathcal{D}^2 \frac{N'_0}{\Theta^3} \frac{B' \sin \alpha}{B_{cr}} F\left(\frac{\nu}{D\nu'_{cr}}\right) \frac{\tanh\left(\frac{R-R_0}{\Delta}\right)}{\cosh^{-2}\left(\frac{R-R_0}{\Delta}\right)} dv d\gamma', \quad (B.15)$$

with the volume element, the function $n(\gamma')$, the constant A, and the magnetic field $B'_0$ given by

$$dv = R^2 \sin \vartheta dR d\vartheta d\varphi \quad (B.16)$$

$$n(\gamma') = \gamma' \sqrt{\gamma'^2 - 1} e^{-\gamma'/\Theta} \quad (B.17)$$

$$A = \frac{\sqrt{3}}{4\pi d^2} \alpha_f mc^2 r_{LC}^3 \quad (B.18)$$

$$B'_{0} = \frac{B_{LC}}{R_{0}^{2}} \sqrt{\left(\frac{R_{0} \sin \vartheta}{\Gamma}\right)^{2} + 1}. \quad (B.19)$$

The doppler factor has been given in Eq. (A.2) and is calculated with the help of Eq. (A.3). The strength parameter $a_{LC}$ and the density in the middle of the current sheet $N'_0$ are given in the text. The expression for the flux also includes the sine of angle $\alpha$, which is calculated by Eq. (B.6). The critical frequency $\nu'_{cr}$ is given in Eq. (B.10), with $B'_0 = B'_0 \sin \alpha$. Finally, the above parameters are calculated at radius $R_0$ with the help of $t_d^*$ from Eq. (B.14). The purpose of the integration on $R$ is to calculate the integral of the emission coefficient across the sheet, when the sheet midplane is at the point $(R_0, \vartheta, \varphi)$. The numerical computation of the above integral for each point in time $t_d^*$, or equivalently each phase in a pulsar period, and each frequency $\nu$, produces the spectra and lightcurves needed.

Finally, the above prescription for the calculation of the pulsar lightcurves is accurate as long as $\Delta < 1/\Gamma$. In the opposite case, the pulse width is not governed by $\Gamma$, but is comparable to $\Delta$, which is not taken into account in the calculation presented here. However, the phase-averaged fluxes calculated by our model are reliable even for $\Delta > 1/\Gamma$.

## References

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 848

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460

* Aliu, E., Anderhub, H., Antonelli, L. A., et al. 2008, Science, 322, 1221

* Bai, X.-N., & Spitkovsky, A. 2010, ApJ, 715, 1282

* Balikhin, M., & Gedalin, M. 2008, J. Plasma Phys., 74, 749

* Beskin, V. S., Kuznetsova, I. V., & Rafikov, R. R. 1998, MNRAS, 299, 341

* Bogovalov, S. V. 1999, A&A, 349, 1017

* Decesar, M. E., Harding, A. K., & Miller, M. C. 2011, in AIP Conf. Ser. 1357, eds. M. Burgay, N. D'Amico, P. Esposito, A. Pellizzoni, & A. Possenti, 285

* Gonthier, P. L., Roberts, J. J., Nagelkirk, E., Stam, M., & Harding, A. K. 2010, in AAS/High Energy Astrophysics Division #11, BAAS, 42, 680

* Harding, A. K., & Muslimov, A. G. 2011, ApJ, 743, 181

* Hoh, F. C. 1966, Phys. Fluids, 9, 277

* Kirk, J. G., & Skjæraasen, O. 2003, ApJ, 591, 366

* Kirk, J. G., Skjæraasen, O., & Gallant, Y. A. 2002, A&A, 388, L29

* Kirk, J. G., Lyubarsky, Y., & Petri, J. 2009, in Ap&SS Lib., 357, ed. W. Becker, 421

* Lind, K. R., & Blandford, R. D. 1985, ApJ, 295, 358

* Lyubarskii, Y. E. 1996, A&A, 311, 172

* Lyubarsky, Y. E. 2003, MNRAS, 345, 153

* Lyubarsky, Y., & Kirk, J. G. 2001, ApJ, 547, 437

* Lyubarsky, Y., & Liverts, M. 2008, ApJ, 682, 1436

* McCann, A. for the VERITAS Collaboration 2011, in Proc. of the 32nd International Cosmic Ray Conference (ICRC2011), held 11–18 August, in Beijing, China, OG2.1–2.2: Cosmic Ray Origin and Galactic Phenomena, 7, 207

* Pétri, J. 2011, MNRAS, 412, 1870

* Pétri, J. 2012, MNRAS, 424, 2023

* Pétri, J., & Dubus, G. 2011, MNRAS, 417, 532

* Pétri, J., & Kirk, J. G. 2005, ApJ, 627, L37

* Pétri, J., & Lyubarsky, Y. 2007, A&A, 473, 683

* Pierbattista, M., Grenier, I., Harding, A., & Gonthier, P. L. 2011, in AIP Conf. Ser. 1357, eds. M. Burgay, N. D'Amico, P. Esposito, A. Pellizzoni, & A. Possenti, 249

* Sironi, L., & Spitkovsky, A. 2011, ApJ, 741, 39

* Spitkovsky, A. 2006, ApJ, 648, L51

* Uzdensky, D. A., & Spitkovsky, A. 2012 [arXiv:1210.3346]

* Venter, C., Harding, A. K., & Guillemot, L. 2009 [arXiv:0911.4890]

* Venter, C., Johnson, T. J., & Harding, A. K. 2012, ApJ, 744, 34

* Watters, K. P., & Romani, R. W. 2011, ApJ, 727, 123

* Zenitani, S., & Hoshino, M. 2008, ApJ, 677, 530