© The Authors 2025

1. Introduction

Rotation-powered young and millisecond pulsars represent the dominant class of identified high-energy gamma-ray sources in the Galaxy (Bulgarelli et al. 2019; Abdollahi et al. 2020). They have been firmly established as such via the detection of short pulses of light modulated with the pulsar spin period. The remarkable coherence and stability of the gamma-ray pulse profile (Kerr 2022) indicates that the emitting zone must be compact and close to the star where the energy densities are highest, naturally pointing towards a magnetospheric origin. How and where pulses are formed in the magnetosphere is a matter of interpretation. In this regard, the gamma-ray emission plays a special role because most of the power radiated away by pulsars is in the gamma-ray band. Current data estimate that about 1% to 100% of the pulsar spindown power (in other words, that of the full energy reservoir of the system) is channeled above 100 MeV gamma rays alone (Abdo et al. 2010a, 2013; Smith et al. 2023). In contrast, the radio emission represents a negligible fraction of the spindown power (≲10⁻⁶, Lorimer & Kramer 2004; Philippov & Kramer 2022). All things considered, understanding how the high-energy emission is produced is crucial to map the main dissipative zones in the magnetosphere.

This information is encoded in the gamma-ray pulse profiles. The majority of observed gamma-ray pulsars present two well-separated pulses per spin period. They are generally not aligned in phase and they are significantly wider than their radio counterparts – if any, setting notable exceptions apart, such as the Crab pulsar, where the radio and the gamma-ray pulses are nearly aligned (Abdo et al. 2010b). Based solely on statistical arguments, observations show that the gamma-ray beam also appears wider than in radio, that is to say, more pulsars are observed in gamma rays than in radio at a fixed sensitivity. These clues suggest that the gamma-ray and radio emitting regions are probably distinct. The polar cap of the star is not a favorable location for the gamma-ray emission because it does not generically produce a double-peaked light curve unless the magnetic axis is nearly perpendicular to the spin axis of the star. In addition, the polar cap is opaque via magnetic conversion to gamma rays above a few GeV. The discoveries of pulsed emission above ≳1 GeV, and up to the TeV range in the Crab and Vela pulsars (Ansoldi et al. 2016; H.E.S.S. Collaboration 2023) irrevocably push the gamma-ray emitting region far away from the stellar surface regardless of the radiation mechanism and, thus, favoring the outskirts of the magnetosphere.

In recent years, the current sheet forming beyond the light-cylinder surface (Michel & Tucker 1969; Michel 1971; Coroniti 1990; Bogovalov 1999; Contopoulos et al. 1999; Spitkovsky 2006) has become one of the main suspect for the location of the gamma-ray emission and, more generally, of the incoherent non-thermal pulsed emission (i.e., from optical to gamma rays, Pacini 1971; Pacini & Salvati 1983; Lyubarskii 1996, as opposed to the coherent radio emission). This scenario offers multiple attractive features, to cite a few: (i) it naturally explains the generic double-peaked structure and the phase-shift with the putative polar-cap emitting radio pulse (Kirk et al. 2002; Bai & Spitkovsky 2010; Contopoulos & Kalapotharakos 2010; Pétri 2011; Arka & Dubus 2013; Kalapotharakos et al. 2014, 2023); (ii) it is optically thin to magnetic conversion process; (iii) it is an obvious location for particle acceleration via magnetic reconnection (Lyubarskii 1996; Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003; Pétri 2012; Uzdensky & Spitkovsky 2014); and (iv) it can explain abrupt swings of the polarization angle (Pétri & Kirk 2005; Słowikowska et al. 2009; Contopoulos & Kalapotharakos 2010; Pétri 2013; Cerutti et al. 2016a; Harding & Kalapotharakos 2017). This scenario is further supported by global particle-in-cell (PIC) plasma simulations of pulsar magnetospheres (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Cerutti et al. 2015; Belyaev 2015; Philippov et al. 2015b; Brambilla et al. 2018; Hu & Beloborodov 2022; Cruz et al. 2024; Soudais et al. 2024). In the limit where the plasma supply is abundant in the magnetosphere, simulations consistently show that reconnection in the current sheet consumes ∼10% of the pulsar spindown power within a few light-cylinder radii. The amount of dissipation is solely governed by the universal relativistic reconnection rate (Cerutti et al. 2020; Hakobyan et al. 2023). The dissipated energy is then efficiently channeled into non-thermal particle acceleration and bright synchrotron radiation. Synthetic pulse profiles reconstructed from the global PIC models present a robust double-peaked pattern in qualitative agreement with the current sheet scenario (Cerutti et al. 2016b; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018, 2023). A pulse of light is received each time the current sheet crosses the observer line of sight, which occurs for most viewing angles twice per spin period.

This work is aimed at putting the global PIC model to the next phase of tests, in light of the most recent gamma-ray observations in the GeV and TeV bands. To this end, we produced synthetic observables of high enough quality, so that they could be directly compared to observations and, in particular, to the third Fermi-LAT catalog (Smith et al. 2023). Here, we focus our attention on the pulse profiles and their energy dependence in both the synchrotron channel and the inverse Compton channel, considering a uniform target photon bath and the anisotropic emission from the neutron star surface. Therefore, our objectives are here twofold: (i) to provide a new set of reference simulations and synthetic observables (Sects. 2–3) and (ii) perform light-curve fitting of the Fermi-LAT gamma-ray pulsar catalog to explore the validity and the implications of the PIC model at the scale of a population (Sect. 4). We summarize our findings in Sect. 5.

2. PIC simulation setup

2.1. Initial configuration

This new series of global models of inclined pulsar magnetospheres was performed with the Zeltron PIC code (Cerutti et al. 2013, 2016b). The numerical setup is similar to the one presented in Cerutti et al. (2016b), which we repeated in this work for completeness. The most important difference is how secondary pairs are injected in the simulation, as explained below.

We used a spherical coordinate system (r, θ, ϕ) where grid cells are equally spaced in log r, θ and ϕ. The grid was composed of 1024 × 256 × 512 cells along the r-, θ-, and ϕ-directions, respectively. The inner radial boundary is fixed at the neutron star surface, r_min = r_⋆. The outer radial boundary is fixed at r_max = 28r_⋆ but a buffer zone encompassing the whole domain absorbs all particles and fields beyond r_abs = 25r_⋆ (Cerutti et al. 2015), which effectively plays the role of an open boundary. Axial symmetry is imposed to the fields at both θ-boundaries (Holland 1983) and periodic boundary conditions are applied along the ϕ-direction.

The initial magnetic field is a pure dipole inclined at an angle χ from the rotation axis, the latter is aligned with the grid axis,

$$\begin{array}{c}\hfill \mathit{B}=\frac{3(\mathit{r}·\mathit{\mu })\mathit{r}}{r^5}-\frac{\mathit{\mu }}{r^3},\end{array}$$(1)

where μ = B_⋆r_⋆³ is the magnetic moment, and B_⋆ is the surface polar magnetic field. The dipolar field is frozen-into the neutron star surface and spins at the neutron star angular velocity, Ω. The light-cylinder radius, where the co-rotation velocity equals the speed of light, is set at R_LC ≡ c/Ω = 5r_⋆. The rotation of the field lines induces an ideal electric field, E = −(Ω×r) × B/c, which is enforced on the inner boundary at each time step. General relativistic effects are neglected in this work (see, however, Philippov et al. 2015a; Philippov & Spitkovsky 2018; Chen et al. 2020; Torres et al. 2024).

The magnetosphere is initially set in vacuum. Plasma is gradually injected into the simulation box via two channels. The first source directly originates from the star surface where a fraction of the surface charge density is delivered at each time step. This density is given by the mismatch between the radial electric field right above the star surface and the ideal co-rotation value according to the jump condition, as in Cerutti et al. (2015). To avoid over-injection, the maximum injected density is limited by the Goldreich-Julian number density, n_GJ = |Ω⋅B|/2πec, where e is the electron charge (Goldreich & Julian 1969). This population was modeled with one particle per cell per species injected at every time step within the first row of cells above the surface. The purpose of this source was to simulate the formation of the primary beam of charges ripped off the stellar crust by the surface electric field (Sturrock 1971; Ruderman & Sutherland 1975).

Once injected, these particles were accelerated to high energies as they undergo a large fraction of the vacuum potential drop of the star. Thanks to the strong surface magnetic field and field line curvature, they trigger intense pair creation (Sturrock 1971; Timokhin & Arons 2013; Timokhin & Harding 2015). This mechanism represents the second and dominant channel of plasma supply in active magnetospheres, which is the main focus of this work. Pair creation is modeled as in Chen & Beloborodov (2014), Philippov et al. (2015b): a pair is produced if the parent particle Lorentz factor, γ, reaches a fraction of the polar-cap potential drop, γ_pc = eμΩ²/m_ec⁴, where m_e is the electron mass. Here, this threshold is fixed at γ_th = 0.05γ_pc and is constant throughout the whole domain. Secondary pairs are boosted along the parent particle direction of motion with a Lorentz factor γ_s = 0.1γ_th. The energy of the secondary pair is removed from the primary particle energy. This simplified procedure is effective at filling the magnetosphere with abundant pairs and it differs from Cerutti et al. (2016b), where a high-multiplicity plasma is injected from the surface without pair creation elsewhere in the magnetosphere. Ions have a dynamically negligible role in the magnetosphere because they are significantly outnumbered by the pairs. In addition, ions do not participate to pair production and to the radiative output due to their heavier mass. Therefore, ions are not included in the simulations presented here (see, however, Chen & Beloborodov 2014; Guépin et al. 2020; Soudais et al. 2024).

The magnetic field strength is limited by the numerical resolution. The fiducial plasma skin-depth scale at the surface of the star, d_e = (m_ec²/4πn_GJe²)^1/2, can be resolved by two cells for a particle with γ = 1. We note that this is a conservative estimate since most pairs will have a Lorentz factor γ ≳ γ_s. In our simulations, we have γ_pc = 2.6 × 10³ so that γ_th = 1.3 × 10² and γ_s = 13. For secondary pairs, the plasma skindepth scale is increased by $d_{\mathrm{e}}^{\mathrm{s}}=\sqrt{{\gamma }_{\mathrm{s}}}{d}_{\mathrm{e}}$, and hence is resolved by about 7 cells. This series of new simulations is complemented by a large split-monopole 3D PIC simulation with χ = 60° already presented in Cerutti et al. (2020). In the latter simulation, the box size extends up to 50 light-cylinder radii. Our purpose is to reanalyse that simulation to explore the evolution of the emitted radiation at large radii in the wind region.

2.2. Synchrotron-curvature emission and feedback

The field strength, the curvature of field lines, and the ultra-relativistic nature of the plasma (γ ≫ 1) make the curvature and synchrotron cooling efficient and dynamically significant in pulsar magnetospheres. The radiative feedback on the particle motion is captured in simulations by including the Landau-Lifshitz radiation-reaction force to the particle equation of motion, g, that depends only on local quantities of the electromagnetic field (Landau & Lifshitz 1971),

$$\begin{array}{c}\hfill \mathit{g}\approx \frac{2}{3}{r}_{\mathrm{e}}^2[(\mathit{E}+\mathit{\beta }\times \mathit{B})\times \mathit{B}+(\mathit{\beta }·\mathit{E})\mathit{E}]-\frac{2}{3}{r}_{\mathrm{e}}^2{\gamma }^2{\stackrel{\sim }{B}}_{\perp }^2\mathit{\beta },\end{array}$$(2)

where r_e = e²/m_ec² is the classical radius of the electron, β = v/c is the particle 3-velocity divided by the speed of light, and

$$\begin{array}{c}\hfill {\stackrel{\sim }{B}}_{\perp }=\sqrt{{(\mathit{E}+\mathit{\beta }\times \mathit{B})}^2-{(\mathit{\beta }·\mathit{E})}^2},\end{array}$$(3)

is the effective magnetic field strength perpendicular to the particle momentum measured in the lab (simulation) frame (Cerutti et al. 2016b). Assuming that the emission is perfectly optically thin, the radiation power spectrum emitted by each simulation particle is given by the classical synchrotron-curvature formula (e.g., Blumenthal & Gould 1970):

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathcal{P}}_{\mathrm{syn}}}{\mathrm{d}\nu }=\frac{\sqrt{3}{e}^3{\stackrel{\sim }{B}}_{\perp }}{m_{\mathrm{e}}{c}^2}\left(\frac{\nu }{{\nu }_{\mathrm{c}}}\right){\displaystyle {\int }_{\nu /{\nu }_{\mathrm{c}}}^{+\infty }}{K}_{5/3}(x)\mathrm{d}x,\end{array}$$(4)

where K_5/3 is the modified Bessel function of 5/3 order, ν is the radiation frequency, and

$$\begin{array}{c}\hfill {\nu }_{\mathrm{c}}=\frac{3e{\stackrel{\sim }{B}}_{\perp }{\gamma }^2}{4\pi m_{\mathrm{e}}c}.\end{array}$$(5)

Given the extreme relativistic boosting effect for γ ≫ 1, we assume that photons are produced along the emitting particle’s momentum. The frequency-integrated radiated power per lepton is

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{syn}}=\frac{2}{3}{r}_{\mathrm{e}}^{2}c{\gamma }^2{\stackrel{\sim }{B}}_{\perp }^2.\end{array}$$(6)

This classical approach is valid if $\gamma {\stackrel{\sim }{B}}_{\perp }/B_{\mathrm{QED}}\ll 1$, where B_QED = m_e²c³/ℏe is the critical magnetic field, which is a good approximation in gamma-ray pulsars, except perhaps close to the star surface where this ratio can be close to unity. Because of the unrealistic field strength accessible in simulations, radiative cooling is rescaled by a large constant factor to achieve the relevant radiation-reaction-limited regime in the magnetosphere (e.g., Uzdensky & Spitkovsky 2014). This regime corresponds to the strong cooling limit in which the acceleration rate is balanced by radiative losses. It sets another energy scale in the problem, γ_rad, given by (Cerutti et al. 2012a)

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{rad}}=\sqrt{\frac{3eE}{2r_{\mathrm{e}}^2{\stackrel{\sim }{B}}_{\perp }^2}},\end{array}$$(7)

such that γ_rad ≪ γ_pc in active pulsars. In the simulations presented here, the classical radius of the electron in Eq. (7) is amplified by a factor f_rad = 5 × 10⁸, giving γ_rad ≈ 4 at the star surface, and ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\approx 45\lesssim {\gamma }_{\mathrm{th}}\ll {\gamma }_{\mathrm{pc}}$ at the light cylinder (assuming $E={\stackrel{\sim }{B}}_{\perp }$), and hence preserving a relevant hierarchy of scales. It was shown by Soudais et al. (2024) that this rescaling procedure does not affect the results for a fixed ratio between the above energy scales. The shortest synchrotron cooling time is resolved by 4 simulation time steps. A modified Boris pusher is used in the simulation to incorporate the full radiation-reaction force along with the Lorentz force in the equation of motion (Tamburini et al. 2010; Cerutti et al. 2013). The total number of iteration per pulsar spin period is then fixed to P/Δt = 2π/ΩΔt ≈ 3.3 × 10⁵.

2.3. Inverse Compton emission

In contrast to synchrotron and curvature radiation, inverse Compton cooling and its feedback onto the dynamic are negligible in isolated systems. The cooling timescale of TeV electrons onto CMB photons or those emitted by the surrounding nebula is way longer than the pulsar spin period. If, however, the background photon field is concentrated on magnetospheric scales, the IR and optical fluxes in the Crab and Vela pulsars can lead to non-negligible losses (H.E.S.S. Collaboration 2023). In both cases though, synchrotron radiation remains the main cooling channel, and thus the radiation-reaction force induced by Compton scattering is neglected in this work. The computation of the emitted spectrum is done under the same assumptions as above: the medium is optically thin and the upscattered photons are beamed along the emitting particle direction of motion. The background radiation field is fixed and photons produced in the simulation do not contribute to it.

While the curvature and synchotron emissivities are fully determined by the simulation native quantities, the geometry and the energy density of the background radiation field are not well constrained and must be prescribed. In this work, we explore two simple geometries: (i) the simulation is immersed into a uniform isotropic photon bath, which could represent CMB photons or diffuse emission emitted at larger scales such as the nebula surrounding the magnetosphere; and (ii) a radial point-like source located at the neutron star center and decaying as 1/r² meant to simulate thermal photons or polar-cap hotspots present on the star surface. The point-like approximation is valid at light-cylinder scales and beyond in the wind current sheet (since R_LC ≫ r_⋆). For simplicity, we assume that the target radiation field has a single photon energy, ϵ₀ (normalized to m_ec²), and we explore both the Thomson and Klein-Nishina regimes.

In the laboratory frame, the inverse Compton spectrum emitted by a relativistic electron scattering off a soft photon with a pitch angle θ₀ between the electron and incoming photon momenta, is given by (Aharonian & Atoyan 1981; Dubus et al. 2008):

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{ic}}}{\mathrm{d}t\mathrm{d}{ϵ}_1}=\frac{2\pi r_{\mathrm{e}}^{2}c}{{\gamma }^2{ϵ}_0}[1+2q(q-1)+\frac{1}{2}\frac{{\left({\mathrm{\Gamma }}_{\mathrm{ic}}q\right)}^2}{1+{\mathrm{\Gamma }}_{\mathrm{ic}}q}],\end{array}$$(8)

where

$$\begin{array}{c}\hfill q\equiv \frac{{ϵ}_1}{{\mathrm{\Gamma }}_{\mathrm{ic}}(\gamma -{ϵ}_1)},\end{array}$$(9)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{ic}}\equiv 2{ϵ}_0^{\prime }=2\gamma (1-cos{\theta }_0){ϵ}_0,\end{array}$$(10)

where ϵ₁ is the upscattered photon energy (normalized to m_ec²), and ϵ′₀ is the (normalized) soft photon energy in the rest frame of the electron. Kinematic constraints yield

$$\begin{array}{c}\hfill 0\le {ϵ}_1\le {ϵ}_{+}=\gamma \frac{{\mathrm{\Gamma }}_{\mathrm{ic}}}{1+{\mathrm{\Gamma }}_{\mathrm{ic}}}.\end{array}$$(11)

These formula are valid in the ultra-relativistic limit (head-on approximation, γ ≫ 1) for both the Thomson (Γ_ic ≪ 1) and the Klein-Nishina (Γ_ic ≫ 1) regimes. Integrating over all frequencies and assuming that ϵ₁ ≫ ϵ₀, the radiated power per lepton is given by (Jones 1965; Ball & Kirk 2000)

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{ic}}={\displaystyle {\int }_0^{{ϵ}_{+}}}{ϵ}_1\frac{\mathrm{d}{N}_{\mathrm{ic}}}{\mathrm{d}t\mathrm{d}{ϵ}_1}\mathrm{d}{ϵ}_1=\frac{2\pi r_{\mathrm{e}}^{2}c}{{ϵ}_0}\mathcal{F}\left({ϵ}_0^{\prime }\right),\end{array}$$(12)

where

$$\begin{array}{c}\hfill \mathcal{F}\left(x\right)=\frac{-10x^4+93x^2+51x(1+x^2)+9}{3x{(1+2x)}^3}\\ \hfill +\frac{(x-3)(x+1)}{2x^2}ln(1+2x).\end{array}$$(13)

Asymptotic limits give

$$\begin{array}{c}\hfill \mathcal{F}(x)\approx \frac{4}{3}{x}^2,x\ll 1,\end{array}$$(14)

valid in the Thomson regime and

$$\begin{array}{c}\hfill \mathcal{F}(x)\approx \frac{1}{2}(ln2x-\frac{5}{6}),x\gg 1,\end{array}$$(15)

that is relevant to approximate the deep Klein-Nishina regime.

For an isotropic photon field, the emitted spectrum is obtained by averaging Eq. (8) over all pitch angles,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{iso}}}{\mathrm{d}t\mathrm{d}{ϵ}_1}=\frac{1}{2}{\displaystyle {\int }_{-1}^1}\frac{\mathrm{d}{N}_{\mathrm{ic}}}{\mathrm{d}t\mathrm{d}{ϵ}_1}\mathrm{d}(cos{\theta }_0).\end{array}$$(16)

Using the variable q defined in Eq. (9) instead of cos θ₀ offers a straightforward way to perform the integration and to recover the well-known Jones (1968) formula,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{N}_{\mathrm{iso}}}{\mathrm{d}t\mathrm{d}{ϵ}_1}& =\frac{2\pi r_{\mathrm{e}}^{2}c}{{\gamma }^2{ϵ}_0}[2q_{j}ln{q}_j+(1+2q_j)(1-q_j)\hfill \\ \hfill & +\frac{1}{2}\frac{{\left({\mathrm{\Gamma }}_j{q}_j\right)}^2}{1+{\mathrm{\Gamma }}_j{q}_j}(1-q_j)],\hfill \end{array}$$(17)

where

$$\begin{array}{c}\hfill q_j\equiv \frac{{ϵ}_1}{{\mathrm{\Gamma }}_j(\gamma -{ϵ}_1)},\end{array}$$(18)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_j\equiv 4\gamma {ϵ}_0,\end{array}$$(19)

and

$$\begin{array}{c}\hfill 0\le {ϵ}_1\le \gamma \frac{{\mathrm{\Gamma }}_j}{1+{\mathrm{\Gamma }}_j}.\end{array}$$(20)

The frequency-integrated power does not admit a simple analytical expression in the isotropic case. In practice, we use the anisotropic expression in Eq. (12) and draw cos θ₀ as a random number uniformly distributed between −1 and 1 for each electron in the simulation.

3. Synthetic pulse profiles

We briefly outline the salient features of the magnetosphere (Sect. 3.1), before proceeding to a detailed analysis of the synthetic synchrotron and inverse Compton emission patterns (Sects. 3.2–3.5), which are the main focus of this part of the work.

3.1. Magnetospheric structure

A nearly steady state was established after a few star spin periods. Figure 1 shows a snapshot of a mature magnetospheric solution with a magnetic obliquity χ = 30°, reached after about two periods once the initial transient has left the domain. The plasma concentrates in the wind current layer that forms beyond the light cylinder where magnetic reconnection, particle acceleration and pair creation take place. The shape of the layer is well captured by the geometry of the magnetic null line predicted by the inclined split-monopole model (Bogovalov 1999), whose location fulfills the condition:

$$\begin{array}{c}\hfill cos\theta cos\chi +sin\theta sin\chi cos[\varphi -\mathrm{\Omega }(t-\frac{r}{c})]=0.\end{array}$$(21)

Fig. 1. Global plasma and magnetic structures forming in a pulsar magnetosphere of magnetic obliquity χ = 30°. Left panel: Three-dimensional rendering of the plasma density (color-coded) normalized by the surface Goldreich-Julian density and multiplied by (r/r⋆)2 to compensate for the effect of expansion. Closed and open field lines encompassing the equatorial current layer are shown by blue lines. The star is the gray sphere in the center. Right panels: Poloidal (top, y = 0) and toroidal (bottom, z = 0 or θ = 90°) slice of the plasma density. The light cylinder is shown with the white dashed line. The dotted lines represent the split-monopole solution, and the magenta arrow shows the orientation of the magnetic moment of the star at this time.

It is worth noting that this solution is a very good approximation even close to the light cylinder, where the dipolar field is still significant. However, we observe a discrepancy in the poloidal plane (y = 0, Fig. 1) where there is an offset in the maximum amplitude of the meanders of the layer. This is especially visible at low inclinations, possibly caused by kink modes growing in the layer (see also in Cerutti et al. 2020), or the effect of a finite magnetization in the wind zone. Overall, the portion of the pulsar wind containing the current layer, or striped region, is well delimited by the expected values: π/2 − χ < θ < π/2 + χ. In the rϕ plane, the current layer splits into two perfect Archimedean spirals of wavelengths 2πR_LC moving in solid rotation with the star. In the equatorial plane, both arms are diametrically opposed from each other (see the z = 0 plane, Fig. 1). The current layer is highly clumpy due to its fragmentation under the tearing instability (Cerutti & Philippov 2017; Cerutti et al. 2020; Hakobyan et al. 2023), which in turn mediates fast reconnection. Underdense regions correspond to magnetic null lines where the field reconnects, while dense regions are flux ropes collecting energetic pairs that have been processed by reconnection.

3.2. Emission pattern and geometric origin of pulses

In agreement with previous studies, dissipation mostly occurs beyond the light cylinder where reconnection kicks in and powers efficient particle acceleration and non-thermal radiation (Cerutti et al. 2016b; Philippov & Spitkovsky 2018). Taking into account time delays due to the propagation of photons towards a distant observer, we compute the emission pattern as a function of the viewing angle α defined with respect to the star rotation axis, and the normalized pulsar phase Φ whose origin is fixed by the plane containing Ω and μ (see Cerutti et al. 2016b for more details). Particles energized by reconnection primarily emit synchrotron radiation near the light cylinder where the field is strongest (see the next section). Synchrotron largely dominates over curvature radiation because, although the particle (guiding center) trajectories are bent near the light cylinder, their radius of curvature (of order R_c ≳ R_LC) is much larger than their Larmor radius (of the order of the current layer thickness). Figure 2 shows the ratio of the effective perpendicular magnetic field to the total field strength. From this figure, we can see that ${\stackrel{\sim }{B}}_{\perp }\lesssim B$ within the current layer (see also Chernoglazov et al. 2023 in their local Harris sheet study), indicating that synchrotron losses are dominant there. Indeed, the ratio of the Larmor radius to the light-cylinder scale is $\xi ={\gamma }_{\mathrm{th}}{m}_{\mathrm{e}}{c}^2/e{\stackrel{\sim }{B}}_{\perp }{R}_{\mathrm{c}}\lesssim 3\times {10}^{-4}$ in the simulation, and therefore the ratio of curvature to synchrotron radiation power scale as P_curv/P_syn = ξ² ∼ 10⁻⁷. In contrast, in the wind zone, particles stream with negligible synchrotron losses because they are drifting at the local E × B-drift velocity, thus ${\stackrel{\sim }{B}}_{\perp }\approx 0$ (see Fig. 2). High-energy particles accelerated in the current layer can also upscatter low-energy photons via inverse Compton scattering.

Fig. 2. Ratio of the effective perpendicular field, ${\stackrel{\sim }{B}}_{\perp }$, to the total field strength, B, in the pulsar midplane (θ = 90°) for the χ = 60° simulation. Blue streamlines represent the magnetic field lines in the plane, and the red dashed circle shows the location of the light cylinder.

Figure 3 shows the emission maps (or skymaps) in the following, for all the inclinations simulated in this work (except the aligned case) in the synchrotron channel (left panels) and in the inverse Compton channel for an isotropic radiation field (right panels). These maps include the contribution from both electrons and positrons. Synchrotron fluxes are frequency-integrated above the fiducial critical frequency, ν₀ = 3eB_⋆/4πm_ec. Similarly, to isolate the high-energy inverse Compton component, we only consider particles with γ > 10 (i.e., including all high-energy pairs and their secondaries with γ ≳ γ_s), and integrate the upscattered photon spectrum ϵ₁ > 30. These emission maps were reconstructed over multiple timesteps to wash out the effect of plasma irregularities and rapid intermittencies inherent to a reconnecting current layer (Andaç et al. 2022), and to collect a large enough sample of photons to limit the shot noise. This averaging is also relevant in comparing our results to observations since measured gamma-ray pulse profiles are reconstructed over a very large number of spin periods due to the poor photon statistics, and therefore they should be considered as mean profiles. It explains the qualitative differences with the results presented in Cerutti et al. (2016b) that were not of good enough quality to perform the quantitative comparison with observations that we propose below in Sect. 4. A synthetic light curve can be generated for any observer by cutting through these emission maps at a fixed value of the viewing angle α. For illustrative purposes, Figure 4 presents a collection of synthetic light curves for synchrotron and inverse Compton radiation generated from the χ = 30° emission maps (considering only the positrons, see the justification later below).

Fig. 3. Synthetic synchrotron (left panels) and inverse Compton (isotropic radiation field, right panels) fluxes as a function of the observer’s viewing angle, α, and pulsar phase, Φ. Each map is normalized by its maximum value and includes the emission from both electrons and positrons. The white dashed line represents the location of the current sheet in the split monopole geometry.

Fig. 4. Synthetic phase-resolved synchrotron (top) and inverse Compton (isotropic radiation field, bottom) fluxes (in arbitrary units) for χ = 30° color-coded by the observer’s viewing angle (emission from positrons only).

The synchrotron and inverse Compton emission maps are very similar and aligned in phase, confirming that both components are produced by the same particle population. Interestingly, the inverse Compton features are sharper, meaning that pulses are systematically narrower than their synchrotron counterpart. The origin of this important feature is further investigated and explained in Sect. 3.5. We also confirm the strong concentration of the emission within the equatorial plane of the pulsar regardless of the magnetic inclination (see Kalapotharakos et al. 2014 in their “FIDO” model, and Cerutti et al. 2016b; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018 for PIC studies), for both the inverse Compton and the synchrotron channels. This effect is due to the higher concentration of the Poynting flux in the equatorial regions and the nearly independent dissipation rate with latitude within the striped wind region (Cerutti et al. 2020). Figure 5 shows the beam correction factor defined as (Watters et al. 2009; Venter et al. 2009),

Fig. 5. Beam correction factor, fΩ, inferred from the synchrotron (top) and inverse Compton (isotropic radiation field, bottom) skymaps for all the magnetic obliquities simulated in this work.

$$\begin{array}{c}\hfill f_{\mathrm{\Omega }}(\chi ,\alpha )=\frac{{\displaystyle \int \int F}(\chi ,{\alpha }^{\prime },\mathrm{\Phi })sin{\alpha }^{\prime }\mathrm{d}{\alpha }^{\prime }\mathrm{d}\mathrm{\Phi }}{2{\displaystyle \int F}(\chi ,\alpha ,\mathrm{\Phi })d\mathrm{\Phi }},\end{array}$$(22)

where F is the photon flux received by a given observer. This quantity corrects for the anisotropy of the emission to infer the true luminosity. It is of the order f_Ω ∼ 1 in the equatorial regions where most of the emission comes from, and it quickly rises for any observer looking at the polar regions where little emission is predicted (see also Kalapotharakos et al. 2023). The emission pattern presents a bright sinusoidal-like component (or “caustic”) at low obliquities which breaks up into two bright hotspots at higher obliquities (χ ≳ 30°). The skymaps present a central symmetry with respect to the center of the maps, located here in the equatorial plane at the pulsar phase Φ = 0.5. Slight asymmetries between the left and the right caustics can be seen in some solutions like in the χ = 45° solution; these are due to statistical fluctuations and we anticipate that they would vanish for a longer integration time. A pulse of light is received by a distant observer each time the current layer passes through the line of sight, which happens in general twice per pulsar period. The central symmetry in the emission maps is due to the periodic change in the direction of the current density direction carried by the layer (mostly oriented towards ±θ directions) between two consecutive crossings, leading to a preferential emission in the northern/southern hemispheres in each hotspot.

Overall, the emission pattern closely follows the geometry of the current layer projected onto the sky. The latter is rather well approximated by the split-monopole solution (Eq. (21)) given by

$$\begin{array}{cc}\hfill \alpha (\mathrm{\Phi },\chi )& =arctan\left(\frac{1}{tan\chi cos2\pi \mathrm{\Phi }}\right),0<\alpha <\pi /2,\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \alpha (\mathrm{\Phi },\chi )& =\pi +arctan\left(\frac{1}{tan\chi cos2\pi \mathrm{\Phi }}\right),\pi /2<\alpha <\pi .\hfill \end{array}$$(24)

It is shown by the white dashed line in Figure 3. We recover the typical double-peaked structure reminiscent of gamma-ray pulsars. For α = 90°, both peaks are perfectly symmetric and separated in phase by ΔΦ ≈ 0.5, regardless of magnetic obliquity. The peaks become unequal due to the central symmetry highlighted above, and their separation in phase decreases for an observer looking away from the equatorial plane. The latter behavior is also consistent with the split-monopole solution that predicts (Pétri 2011)

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Phi }(\alpha ,\chi )=\frac{1}{\pi }arccos|cot\alpha cot\chi |.\end{array}$$(25)

As expected, pulsations disappear outside the domain of the striped wind region (α ≲ π/2 − χ or α ≳ π/2 + χ).

However, what the split-monopole model does not predict is the presence of a significant interpulse, which sometimes may even appear as a third peak, which is particularly visible in the synchrotron light curves (see Figs. 4–8). This is also a common feature observed in gamma-ray pulsars. This interpulse should not be confused with the bridge emission that arises when the observer intercepts the outer boundaries of the striped wind region (π/2 − χ ≲ α ≲ π/2 + χ) that is as visible at low obliquities, but that disappears at high obliquities (Fig. 3, Cerutti et al. 2016b; Philippov & Spitkovsky 2018; this feature is also observed in many pulsars. The presence of a third peak is a clearly recognizable, yet mysterious, property of the Vela gamma-ray pulse profile (Grenier et al. 1988; Pellizzoni et al. 2009; Abdo et al. 2010c). We checked that this third component is emitted in the wind zone by cutting off all the emission coming from inside the light cylinder. It is mostly prominent at high-inclinations (χ ≳ 30°) in the synchrotron channel and seems to originate from the wind region contained within the meanders of the layer for χ = 45° −75°, but it is unclear how particles are accelerated and heated (i.e., ${\stackrel{\sim }{B}}_{\perp }\ne 0$) in this sector. It is possible that this particular region is under the strong influence of the current layer that surrounds it. However, the χ = 90° solution shows that the interpulse emission is distributed at intermediate latitudes far from the boundaries of the layer that extends all the way to the poles, approximately at the same location as for χ = 60° and 75°. Another possibility is that high-energy particles accelerated in the sheet are escaping in the upstream medium, as recently shown in local 3D reconnection studies (Zhang et al. 2021), but the extension in phase of that component suggest that high-energy particles are distributed all over the wind region, which is difficult to reconcile with the latter scenario. In addition, this effect would probably be limited to ions only due to strong radiative losses experienced by the leptons (Chernoglazov et al. 2023).

What appears to be the most convincing scenario is that the Vela-like interpulse is produced within the volumetric return current. Figure 6 shows the spatial distribution of the synchrotron power emitted near the light cylinder by positrons (this feature is mostly produced by positrons). Outside the current sheet, the two hotspots coincide with regions of strong return current connecting the polar cap to the current sheet. This special region appears at large obliquities and is known to accelerate and create pairs because the amount of charges extracted at the polar cap are insuffisant to sustain the current required by the magnetosphere (Timokhin & Arons 2013; Philippov et al. 2015b). Interestingly, this scenario also explains the origin of the high-viewing angle components in the skymaps (Φ ∼ 0.4,  0.9, α ≳ 50° , 140°). These regions are consistent with polar-cap field lines carrying super-Goldreich-Julian current densities, and they are mostly sustained by electrons. In the orthogonal rotator, these components become the volumetric return current by symmetry. It is unclear though how particle acceleration along these field lines remains active this far from the star surface. We leave this question to a future study.

Fig. 6. Spatial distribution of the synchrotron power per unit of volume, d𝒫syn/dV, emitted by the positrons only at the light cylinder for χ = 60°. This figure highlights the role of the volumetric return current in explaining the Vela-like interpulse.

What is clear however, is that this component is emitted in a special part of the inter-layer wind. In contrast to the (unstriped) polar regions, the wind between two consecutive layers is not an homogeneous medium, there is a strong plasma density and magnetic field strength gradient. The wind region leading the current layer is highly depleted in plasma and magnetic field strength, while the wind region trailing the layer contains denser, more energetic plasma, and stronger fields (see Cerutti & Philippov 2017; Cerutti et al. 2020). All in all, these clues suggest that this mysterious third peak may be a by-product of the wind dynamic, plasma reorganization induced by the presence of the reconnecting layers, and the overall magnetospheric current density.

Figure 7 graphically summarizes how pulses form in the simulations. As pointed out above, the shape of the current layer is remarkably well fitted by the split monopole solution, even close to the light cylinder. Therefore, mapping the emission onto this surface is sufficient to render the main features of the light curves. A fixed observer probes a surface of constant latitude θ = α, then the equation of each layer is given by

Fig. 7. Schematic representation illustrating how pulses form relying on simple geometric considerations based on the split-monopole configuration. Here, χ = 60° for illustrative purposes, but this picture applies for all inclinations (except the aligned and orthogonal rotators). In this diagram, the pulse P1 (P2) is emitted once the sheet 1 (respectively 2) rotating at the angular velocity of the star, Ω, overlaps with the fixed black dashed line, an iso-phase contour set by the observer. Three typical configurations are presented: (i) an equatorial view (α = π/2) leads to two symmetric pulses separated by ΔΦ = 0.5, (ii) an off-equator view within the striped wind region (α = π/2 − 0.8χ) where the pulses are asymmetric and separated by ΔΦ = 0.28 and contains significant emission in between (a bridge or third peak, green region), and (iii) the degenerate single pulse case (ΔΦ = 0) where the observer’s line of sight grazes the outer boundary of the stripes (α = π/2 − χ).

$$\begin{array}{c}\hfill \varphi (r)=\pm arccos(-cot\alpha cot\chi )+\mathrm{\Omega }t-r/R_{\mathrm{LC}},\end{array}$$(26)

which corresponds to an Archimedean spiral (for sheet 1 and 2 in Fig. 7). The spiral pattern is moving in solid rotation with the star – but not the plasma – in virtue of Ferraro’s isorotation law. Thus, the emitting zone of photons arriving at the same phase at the observer has an Archimedean spiral shape as well (dashed line in Fig. 7). A simple way to see this is to consider a photon emitted radially at a time t towards the observer located along the ϕ = ϕ₀ direction. After a time δt, the photon will propagate by a distance δr = cδt while the emitting zone will have shifted by δϕ = ϕ − ϕ₀ = Ωδt. At all radii where new photons can be created, the location of all the emitting zones probed by the observer at a given time draws a spiral of equation ϕ = ϕ₀ − r/R_LC. Contrary to the sheets, this pattern is fixed in time since the observer is assumed static in the simulation frame. In this geometric construction, a pulse of light will be observed when the sheet overlaps with the observer’s spiral. This happens twice per spin period within the striped wind region. At α = π/2 ± χ, the layer degenerates into a single sheet leading to a single pulse. In this framework, interpulse emission can be produced in between the layers.

3.3. Asymmetries and energy dependence

Low-inclination solutions present additional emission features visible in the synchrotron maps. They appear as two humps, one above (in the Φ < 0.5 sector), and one below (in the Φ > 0.5 sector) the main emission pattern from the layer (Fig. 3). This feature was already reported in Cerutti et al. (2016b). It corresponds to high-energy electrons accelerated near the base of the current layer precipitating towards the star. This asymmetry between both species is explained by the overall polarisation of the magnetosphere: the plasma is negatively charged along magnetic polar regions and is positively charged in the equatorial regions (and the other way around, if Ω ⋅ μ < 0). Although this global polarization effect should remain at realistic scales, the relative contribution from precipitating electrons to the total emission synchrotron output may be overestimated due to the modest plasma multiplicities achieved in the simulations, κ ≡ n/n_GJ ≳ 5 (Fig. 1), as opposed to observed gamma-ray pulsars where κ ≳ 10² − 10⁶. A higher plasma multiplicity should lead to a smaller fraction of precipitating particles.

Although our scale separation is modest, we observe a significant dependence of the pulse profile with the photon energy. Figure 8 represents a synthetic light curve that is reminiscent of the Vela gamma-ray pulse profile, which is composed of two main pulses (P1 and P2 separation by ΔΦ ≈ 0.43) and a weaker broad peak in between (P3) obtained from the χ = 60° solution for a viewing angle α = 110°. At low frequencies (ν/ν₀ ≲ 0.1), P1 slightly dominates over P2 in height. Both peaks are connected by a smooth, asymmetric bridge. P3 appears for ν/ν₀ ≳ 0.1 close to P1 and drifts away towards P2 at higher frequencies. The amplitude of P1 relative to P2 dramatically decreases for ν/ν₀ ≳ 1 (see also Brambilla et al. 2015; Barnard et al. 2022). P1 and then P3 nearly completely vanish at the highest frequencies where only P2 remains. These features are in very good agreement with Fermi-LAT observations of Vela (Abdo et al. 2010c) and more generally of gamma-ray pulsars (Smith et al. 2023). This clean cut property can be drawn in the simulations only if considering the positronic emission from the magnetosphere. Adding the contribution from electrons makes this conclusion less robust because of the precipitating particles, which we argue may be overproduced in the simulation.

Fig. 8. Energy dependence of the synthetic synchrotron pulse profile with the photon frequency, ν normalized to ν0, for a Vela-like solution (χ = 60°, α = 110°). Top panel: Light curves equally spaced in log ν/ν0 from the lowest frequencies (bottom) to the highest frequencies (top). The bottom panel shows the full phase-frequency emission map. In both panels, light curves are normalized to the maximum flux for each frequency band. This figure only shows the contribution from positrons.

3.4. Radiated power and efficiencies

Figure 9 presents how is the radiated power distributed in the wind region at all radii probed by the simulations. The large-scale split-monopole simulation is shown for comparison and highlights trends at large radii. As expected, the synchrotron power rapidly falls off with radius because of the decrease of the magnetic field strength with radius, since ${\mathcal{P}}_{\mathrm{syn}}\propto {\stackrel{\sim }{B}}_{\perp }^2$ (Eq. (6)). The radial evolution of the effective perpendicular field felt by the particles within the sheets is well modeled by the superposition of the initial vacuum poloidal field, B_pol, and a purely toroidal component, B_ϕ, where both components are equal at the light cylinder, such that

$$\begin{array}{c}\hfill {\stackrel{\sim }{B}}_{\perp }^2\sim B_{\mathrm{pol}}^2+B_{\varphi }^2\sim \frac{B_{\mathrm{LC}}^2}{2}(\frac{1}{{\widehat{r}}^{2a}}+\frac{1}{{\widehat{r}}^2}),\end{array}$$(27)

Fig. 9. Radial distribution of the frequency-integrated radiated power in the synchrotron channel (top panel), in the inverse Compton channel for both an isotropic (middle panel) and an anisotropic (from the star, bottom panel) target photon field in the Thomson (red lines) and deep Klein-Nishina (blue lines) regimes. Simulations with an initial dipolar field are shown by solid lines (χ = 30°, 45°, 60°), while the split-monopole simulation is shown by dashed lines. Simple analytical expressions are also overplotted for comparison (black dashed-dotted and dotted lines). All quantities are normalized by their value at the light cylinder.

where $\widehat{r}=r/R_{\mathrm{LC}}$, a = 3 for the dipolar simulations, and a = 2 for the split-monopole simulation, and B_LC = B_⋆(r_⋆/R_LC)^a is an estimate of the field strength at the light cylinder.

The synchrotron power also depends on the square of the mean particle Lorentz factor in the current layer, but it increases at best as $ln\widehat{r}$ (Cerutti et al. 2020). This weak dependence with radius is well visible in the inverse Compton power for an isotropic radiation bath. In this case, the emitted power depends only on the particle Lorentz factor, as 𝒫_ic ∝ γ² in the Thomson regime and 𝒫_ic ∝ lnγ in the Klein-Nishina regime. In both regimes, the inverse Compton power slowly increases with radius. In contrast, if the target photon field originates from the star, the Compton power steeply drops with radius. This is, for one part, due to the decrease of the target photon number density with radius as $1/{\widehat{r}}^2$ and, for the other part, due to nearly rear-end collisions between the stellar photons and the particles that are mostly moving radially outward in the wind. The latter effect is most pronounced in the Thomson regime. This is explained by the dependence of the (anisotropic) inverse Compton power with ϵ′₀ = γ(1−cos θ₀)ϵ₀. Assuming that the particle direction of motion is given by the V_D = cE × B/B² drift velocity of the monopole wind and the purely radial stellar photons, then the pitch angle is given by (e.g., Cerutti & Beloborodov 2017)

$$\begin{array}{c}\hfill cos{\theta }_0={\mathit{V}}_{\mathrm{D}}·{\mathit{e}}_{\mathit{r}}/V_{\mathrm{D}}={(1+\frac{1}{{\widehat{R}}^2})}^{-1/2}\approx 1-\frac{1}{2{\widehat{R}}^2},\end{array}$$(28)

where $\widehat{R}=\widehat{r}sin\theta \sim \widehat{r}\gg 1$ is the cylindrical radius, and hence

$$\begin{array}{c}\hfill {ϵ}_0^{\prime }\propto \frac{1+ln\widehat{r}}{{\widehat{r}}^2},\end{array}$$(29)

where we assume once again that γ increases with $ln\widehat{r}$. In the Thomson regime, 𝒫_ic ∝ ϵ′₀² (Eq. (14)), while 𝒫_ic ∝ ln(2ϵ′₀) (Eq. (15)) when Klein-Nishina corrections are important. Thus, anisotropic effects are reduced in the Klein-Nishina regime and result in a slower decay of the radiated power than in the Thomson regime. Note that there is no obvious dependence with the magnetic inclination angle.

We can deduce from the above scaling laws an approximate expression for the synchrotron radiative efficiency and the radial evolution of the total (frequency-integrated) radiated power, and compare it with the other energy channels. Figure 10 shows the radial evolution of the outward radial Poynting (L_EM) and pairs kinetic energy (L_part) fluxes passing through a spherical shell of radius r as reported in Cerutti et al. (2020). In this work, we compute the contribution of the radiated synchrotron power, L_syn, and neglect the power lost via inverse Compton scattering. Assuming that the radiated power mostly depends on the evolution of ${\stackrel{\sim }{B}}_{\perp }^2$, using Eq. (27), we can derive a crude analytical expression for L_syn normalized to the pulsar spindown power, L₀,

Fig. 10. Radial evolution of all energy channels in the split-monopole solution, including the synchrotron luminosity, normalized to the pulsar spindown power (L0). The inferred reconnection rate is βrec ≈ 0.12 and the radiative efficiency is ηrad ≈ 0.05. The inverse Compton power is not shown since radiative energy losses through this channel are neglected in this work.

$$\begin{array}{c}\hfill \frac{L_{\mathrm{syn}}}{L_0}\approx {\eta }_{\mathrm{rad}}\frac{{\displaystyle {\int }_1^{\widehat{r}}}{\stackrel{\sim }{B}}_{\perp }^2\mathrm{d}\widehat{r}}{{\displaystyle {\int }_1^{+\infty }}{\stackrel{\sim }{B}}_{\perp }^2\mathrm{d}\widehat{r}},\end{array}$$(30)

$$\begin{array}{c}\hfill \frac{L_{\mathrm{syn}}}{L_0}\approx {\eta }_{\mathrm{rad}}(1-\frac{1}{\widehat{r}}),\end{array}$$(31)

for a pure toroidal field, or

$$\begin{array}{c}\hfill \frac{L_{\mathrm{syn}}}{L_0}\approx {\eta }_{\mathrm{rad}}(1-\frac{5}{6\widehat{r}}-\frac{1}{6{\widehat{r}}^5}),\end{array}$$(32)

for a toroidal and dipolar field components. The radiative efficiency, η_rad, is defined as the ratio of the asymptotic synchrotron power to the pulsar spindown power, ${\eta }_{\mathrm{rad}}=L_{\mathrm{syn}}^{\infty }/L_0$. The spindown power is measured here as the Poynting flux at the light cylinder, and thus assumes that there is negligible dissipation within the light cylinder. For the split-monopole simulation, dissipation depends mildly on the magnetic inclination and the radiative efficiency is η_rad ≈ 5% (Fig. 10). In contrast, the dipole simulations present a significant variation of the radiative efficiency that is not translated in Eq. (32), it varies from about ∼6% in the aligned case down to ≲1% for the orthogonal rotator (see also, Cerutti et al. 2016b).

Following Cerutti et al. (2020), the Poynting flux evolution with radius is approximately given by

$$\begin{array}{c}\hfill \frac{L_{\mathrm{EM}}}{L_0}=1-{\beta }_{\mathrm{rec}}ln\widehat{r},\end{array}$$(33)

where β_rec is the dimensionless reconnection rate that quantifies how efficient reconnection is at dissipating the Poynting flux. By virtue of energy conservation, the amount of power dissipated and channeled to particle kinetic and synchrotron power is then

$$\begin{array}{c}\hfill \frac{L_{\mathrm{part}}+L_{\mathrm{syn}}}{L_0}={\beta }_{\mathrm{rec}}ln\widehat{r}.\end{array}$$(34)

In the extreme case where all the dissipated power goes directly into radiation, that is to say, if L_part = 0, then using Eq. (31), we see that η_rad ≤ β_rec, meaning that the maximum bolometric radiative efficiency is bounded by the reconnection rate. Here, η_rad ≈ 0.05 < β_rec ≈ 0.12 in the split-monopole simulation. However, gamma-ray pulsars are most likely in this extreme cooling regime, as also indicated by recent full-scale hybrid simulations of millisecond pulsars (Soudais et al. 2024), such that η_rad = β_rec. Then, assuming a purely toroidal field geometry yields

$$\begin{array}{c}\hfill \frac{L_{\mathrm{syn}}}{L_0}={\beta }_{\mathrm{rec}}(1-\frac{1}{\widehat{r}}),\end{array}$$(35)

and thus using Eq. (34), the particle kinetic energy flux evolution can be inferred as

$$\begin{array}{c}\hfill \frac{L_{\mathrm{part}}}{L_0}={\beta }_{\mathrm{rec}}(ln\widehat{r}+\frac{1}{\widehat{r}}-1).\end{array}$$(36)

This important result can be recovered with the following argument for an aligned rotator. In the most efficient synchrotron cooling regime, the bulk of particles accelerated at the base of the current layer at the light cylinder will be in the radiation-reaction-limited regime, which implies that the power radiated away per electron, 𝒫_syn, balances the acceleration rate estimated as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{syn}}\sim ec{E}_{\mathrm{rec}},\end{array}$$(37)

where E_rec = β_recB_LC is the reconnection electric field accelerating particles in the layer. Assuming that the emitting particles are confined within a flat layer of inner radius R_LC and outer radius 2R_LC of thickness δ, the total number of particles is

$$\begin{array}{c}\hfill N_{\mathrm{part}}\sim \kappa n_{\mathrm{GJ}}(4\pi R_{\mathrm{LC}}^2-\pi R_{\mathrm{LC}}^2)\delta ,\end{array}$$(38)

where δ is the layer thickness. According to Eq. (35), 50% of the synchrotron power is radiated away within a few light-cylinder radii, suggesting that the synchrotron emission is highly concentrated close to the base of the layer as assumed here. Following Cerutti & Philippov (2017), the layer thickness at the light cylinder is given by δ ∼ R_LC/κ. With n_GJ = ΩB_LC/2πec, and ignoring factors of order unity yields

$$\begin{array}{c}\hfill N_{\mathrm{part}}\sim \frac{\mathrm{\Omega }{B}_{\mathrm{LC}}{R}_{\mathrm{LC}}^3}{\mathit{ec}}.\end{array}$$(39)

Then, the total synchrotron power can be estimated as $L_{\mathrm{syn}}\sim N_{\mathrm{part}}{\mathcal{P}}_{\mathrm{syn}}^{\mathrm{LC}}$. Noticing that B_LC ∼ B_⋆r_⋆³/R_LC³ = μ/R_LC³, we obtain

$$\begin{array}{c}\hfill L_{\mathrm{syn}}\sim {\beta }_{\mathrm{rec}}{L}_0,\end{array}$$(40)

where L₀ = μ²Ω⁴/c³ is the force-free spindown of the aligned rotator. A similar argument is also exposed in Uzdensky & Spitkovsky (2014).

The fact that all pulsars should have the same bolometric radiative efficiency may appear in contradiction with the Fermi-LAT data that show a clear decrease of the gamma-ray efficiency above 100 MeV with increasing spindown power, ranging from ∼10 − 100% for millisecond pulsars down ∼1% in powerful Crab-like pulsars. However, the third Fermi-LAT catalog also shows that high-spindown pulsars present significantly softer spectra than millisecond pulsars (Smith et al. 2023), suggesting that a sizeable fraction of the radiative output is below the Fermi-LAT threshold and thus missed in the estimation of the bolometric luminosity. This trend may be explained by a higher yield of pair production in high spindown power pulsars. Secondary pairs then redistributes the radiative power at lower frequencies over a broader spectral range. In this regard, a telescope operating in the MeV range would bring valuable constraints for the models. Unfortunately, given the limited scale separation achieved in the simulations and the simplified pair production model adopted in this work, we cannot elaborate more on this scenario at this stage.

3.5. Pulse width

Another important feature of gamma-ray pulsars is the narrowing of the pulse width with increasing photon energy (Abdo et al. 2010c,b,d). This property is also reproduced by the simulations. Figure 11 shows the evolution of the full width at half maximum for a Vela-like configuration with the photon energy for the synchrotron and the inverse Compton (for an isotropic radiation bath) emission channels. Below ν₀, the synchrotron pulse width represents about ΔΦ ≈ 0.1 in phase, and drops down to ΔΦ ≈ 0.02 at the high-end of the synchrotron spectrum (ν/ν₀ ≳ 10). The inverse Compton pulse is about twice thinner than the synchrotron pulse at low energy, and decreases in a similar manner as the synchrotron pulse width with increasing energy (ϵ₁ ≳ 10). Since the emission originates from high-energy particles accelerated in the wind current layer, this result indicates that the highest-energy particles are concentrated deep into the layer. This is reminiscent of the kinetic beaming effect reported in local reconnection studies (Cerutti et al. 2012b; Mehlhaff et al. 2020), where the angular distribution of the particles decreases with energy. This effect is operating in the strong synchrotron cooling regime because particles do radiate high-energy photons before they isotropize in plasmoids (in the wind comoving frame). This regime is relevant near the base of the current layer.

Fig. 11. Narrowing of the pulse full width at half maximum with the energy of the emitted photons in the inverse Compton (isotropic radiation field, red line) and synchrotron channels (blue line) for the dipole simulation with χ = 60° and α = 110°.

The acceleration of the wind bulk Lorentz factor with radius also leads to the sharpening of the pulses due to relativistic beaming effects (Kirk et al. 2002; Pétri 2011). Figure 12 shows the emitted pulse profile as a function of the radius of emission in the wind within a Δr = R_LC thick radial bin. The synchrotron and inverse Compton peaks become thinner with radius. In a given radial bin, the pulses have a similar width in both channels. However, the observed synchrotron light curve presents broader peaks because most of the synchrotron emission is produced near the light cylinder where the wind is slowest. In contrast, the observed inverse Compton light curve for a uniform photon field has narrower peaks because the majority of the flux is produced at larger radii where the bulk Lorentz factor is higher. This is consistent with the increase of the inverse Compton power with radius reported in Sect. 3.4 and in Figure 9. This effect explains the narrower pulses in the inverse Compton channel compared to the synchrotron light curve (Fig. 11). Conversely, if the target photon field originates from the star, the inverse Compton pulse widths are comparable to the synchrotron ones, hence confirming the important role of the photon field geometry in shaping the pulse profile.

Fig. 12. Evolution of the observed light curve as a function of the radius of emission. Top panel: Synthetic synchrotron and inverse Compton (isotropic radiation field) light curves as a function of the emitting radius, from the light-cylinder radius up to 5LC for the dipole simulation (χ = 60°, α = 90°). Bottom panel: Full width at half maximum of the inverse Compton pulse width as a function of radius for the dipole (magenta crosses) and split-monopole (blue dots) simulations. The red dotted and dashed lines shows the two asymptotic behaviors discussed here. The black vertical dashed line marks the location of the critical radius rm where the emitted pulses are the narrowest.

The above conclusion is valid as long as the geometric thickness of the layer is very small compared to the separation between two successive current layers in the wind, meaning Δ ≡ δ/πR_LC ≪ 1. In this regime, the pulse width is shaped by relativistic beaming effects. Following Pétri (2016), the pulse width emitted by an infinitely thin Archimedean spiral shape current sheet is given by

$$\begin{array}{c}\hfill \delta \mathrm{\Phi }\sim \frac{1}{\pi \mathrm{\Gamma }},\end{array}$$(41)

where Γ is the bulk Lorentz factor of the emitting particles. In the split-monopole solution, the latter evolves as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{M}}=\sqrt{1+{\widehat{R}}^2},\end{array}$$(42)

inside the fast magnetosonic point and saturates beyond this point to ${\mathrm{\Gamma }}_{\infty }\sim {\mu }_{\mathrm{M}}^{1/3}={(B^2/4\pi n{m}_{\mathrm{e}}{c}^2)}^{1/3}$ (Tomimatsu 1994; Beskin et al. 1998; Cerutti et al. 2020). The law given in Eq. (41) explains well the synthetic pulse width and its radial evolution near the light cylinder (r ≲ 5R_LC, see Fig. 12) if the particles in the layer begins with a mildly relativistic boost,

$$\begin{array}{c}\hfill \mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{LC}}{\mathrm{\Gamma }}_{\mathrm{M}},\end{array}$$(43)

where Γ_LC ≈ 5.5 best reproduce the reported simulation data. At larger radii, as the wind expands, the layer thickness increases with radius, then the pulse width due to the finite thickness of the layer is (Coroniti 1990; Cerutti & Philippov 2017)

$$\begin{array}{c}\hfill \delta \mathrm{\Phi }\sim \frac{\mathrm{\Delta }}{2}=\frac{1}{2\pi {\mathrm{\Gamma }}_{\mathrm{LC}}\kappa }\widehat{r},\end{array}$$(44)

meaning that the pulse width widens linearly with radius. In the simulations, we have $\mathrm{\Delta }/2\approx 0.004\widehat{r}$ (Fig. 12). Thus, there is a radius, r_m, where the pulse width is minimum and beyond which the effect of expansion dominates over relativistic beaming. Matching Eq. (41) with Eq. (44) gives this critical radius (for θ = π/2)

$$\begin{array}{c}\hfill r_{\mathrm{m}}=\sqrt{\frac{\sqrt{1+16{\kappa }^2}-1}{2}}\approx \sqrt{2\kappa },\mathrm{for}\kappa \gg 1,\end{array}$$(45)

if r_m is inside the fast magnetosonic point, else

$$\begin{array}{c}\hfill r_{\mathrm{m}}=2\kappa \left(\frac{{\mathrm{\Gamma }}_{\mathrm{LC}}}{{\mathrm{\Gamma }}_{\infty }}\right).\end{array}$$(46)

Only the split-monopole simulation is large enough to observe the effect of the layer expansion on the pulse width. In this simulation, the bulk Lorentz factor saturates $\widehat{r}\gtrsim 10$, thus r_m is best estimated with Eq. (44). With κ ∼ 7 (Fig. 1) gives r_m ≈ 3.7 (Fig. 12). For more realistic multiplicities (κ ∼ 10² − 10⁶), r_m ∼ 10 − 10³. Pulsations fully disappear if Δ = 1 at

$$\begin{array}{c}\hfill r_{\mathrm{c}}=\pi {\mathrm{\Gamma }}_{\mathrm{LC}}\kappa ,\end{array}$$(47)

so that r_c ∼ 10² − 10⁶ in active pulsars.

3.6. Implications for the pulsed TeV emission

This work constrains the location of the pulsed TeV emission discovered in the Crab and Vela pulsars within the inner parts of the pulsar wind. The target photon field must therefore remain compact. A promising candidate for emitting this low-energy radiation are secondary pairs produced in the vicinity of the current sheet – like a sheath surrounding the layer – near the light cylinder and shining synchrotron radiation (Lyubarskii 1996; Hakobyan et al. 2019). This photon field would appear isotropic seen from the high-energy pairs accelerated deep inside the reconnection layer, which could then be upscattered into the TeV range. In contrast, any diffuse target photon field on much larger scales, from the surrounding nebula, the CMB, or from a nearby main-sequence companion star, would smear out any pulsation by adding up a dominant contribution from large radii.

The detection of 20 TeV pulsed emission from the Vela pulsar implies that pairs must be accelerated up to at least γ_max ≳ 4 × 10⁷ (H.E.S.S. Collaboration 2023). This is a challenging constraint for the models because this energy scale is significantly above the radiation-reaction-limited particle energy estimated at the light cylinder, ${\gamma }_{\mathrm{rad}}^{\mathrm{LC}}\approx 4\times {10}^5\ll {\gamma }_{\max }$ (assuming $E={\stackrel{\sim }{B}}_{\perp }=B_{\mathrm{LC}}\approx 5\times {10}^4$G, see Eq. (7)). One possible explanation for this discrepancy involves X-point acceleration, where radiative losses are substantially reduced, meaning that particle acceleration could continue nearly unimpeded up to the upstream magnetization limit, γ_max ≳ σ_LC (Cerutti et al. 2012a; Kagan et al. 2016; Chernoglazov et al. 2023). This would in turn imply that the magnetization must be much higher than usually assumed in Vela.

Another effect that could alleviate this issue is particle reacceleration at larger radii, as reported in Sect. 3.5, because reconnection pursues magnetic dissipation while radiative losses diminish thanks to the weakening of the magnetic field strength with distance. Note that the spatial offset between the peak of the synchrotron and the inverse Compton emission would not lead to any significant changes in pulse phase, thanks to the caustic effect highlighted in Sect. 3.2. However, it remains uncertain whether this effect alone can account for the required hundred times increase in particle energy, from γ_rad ≈ 4 × 10⁵ to ≈4 × 10⁷, on reasonably compact scales (i.e., $\widehat{r}\lesssim r_{\mathrm{m}}$). Despite these challenges, a combination of both mechanisms may render this scenario viable, but further investigation is necessary to explore these questions in greater detail.

4. Comparison with the third Fermi pulsar catalog

Here, we compare the synthetic light curves predicted from the synchrotron emission with the observed phase-folded light curves of the 294 pulsars listed in the third Fermi-LAT pulsar catalog (3PC, Smith et al. 2023). We keep only the positron emission since the emission from the precipitating electrons is likely to be strongly overestimated (Sect. 3.3). Skymaps are frequency-integrated above ν₀ defined in Sect. 3.2, which represents the typical synchrotron frequency emitted by pairs energized at the light cylinder. We use this frequency as a proxy for the Fermi-LAT low-energy threshold.

4.1. Model fitting procedure

The synthetic skymaps were rebinned by a factor of 2 in phase and inclination α to increase their signal-to-noise ratio (S/N). We also assumed that both hemispheres have the same emission properties and, hence, enforced the skymap symmetry by averaging the values at (Φ,α) and (Φ + 180°,180° −α), further increasing the S/N of the skymaps and reducing the scope of the parameter range to explore.

The observed folded light curve consists of n phase measurements m_i with an associated error σ_i. The best-fit synthetic light curve was found by minimizing the chi-square,

$$\begin{array}{c}\hfill C^2={\displaystyle \sum _{i=1}^n}\frac{{(m_i-K{S}_i-\mathcal{B})}^2}{{\sigma }_i^2},\end{array}$$(48)

where K is a constant scaling factor, and ℬ is a constant background. Some pulsars present a strong DC component that may be of magnetospheric origin, but the model presented here does not predict any significant contribution to the background emission, we are thus considering ℬ independently from the physical model. S_i is the synthetic light curve at phase i for a given obliquity χ, inclination α, dephased by Φ₀ to adjust the data. The origin of the phase in the 3PC is usually set by the main radio pulse (for radio-loud pulsars), while it is set by the magnetic pole in the model. While these two definitions may be compatible if the radio pulse originates from the polar cap of the star (assuming a dipole), they should differ otherwise. Given the large diversity of the radio phenomenology (radio quiet/loud, number of pulses, phase-(mis)alignment with the gamma rays), we chose to treat the origin of phases as a free parameter in the fitting procedure. C² can be analytically minimised for the nuisance parameters K and ℬ. However, a full search must be performed over the parameter space in Φ₀, χ and α to find the minimum. We found that fitting with equal-sized phase bins gave acceptable results so we did not pursue more elaborate fitting routines aiming to take better into account narrow, high-amplitude pulses such as bin weighting (Romani & Watters 2010) or fixed count binning (Pierbattista et al. 2015). Although this allows a first confrontation of the model to the observations, an important caveat is that the fit is usually not good statistically given the quality of the Fermi-LAT light curves and the limitations of the model simulation (e.g. Romani & Watters 2010; Johnson et al. 2014; Pierbattista et al. 2015).

4.2. Best fit light curves

Figure 13 shows examples of best fits to the Fermi-LAT pulse profiles integrated above 100 MeV (see the full atlas). The obliquity χ and inclination α of the best fit are shown on top of each panel. The top row shows double-peaked pulse profiles, obtained when the line-of-sight is well within the opening angle of the current sheet. The best fit for the Vela pulsar (PSR J0835-4510) is consistent with the values used to illustrate the energy dependence of the pulse profile (Fig. 8) since α ≈ 70° is equivalent to α ≈ 110°  = 180° −70° with the map symmetry. Overall, our fit to the Vela pulse profile compares well with, e.g., the elaborate slot gap, outer gap, and two-pole caustic model fits in Barnard et al. (2016). Our values of α and χ are also consistent with their values. The bridge emission and the peak positions are well reproduced by the model. However, the observed peaks are narrower and their amplitude is greater than predicted. This is also seen in the fits for the two other pulsars on the top row of Figure 13, PSR J1028–5819 and PSR J1418–6058, both of which are seen at high inclination (leading to two pulses). This shortcoming is likely due to the limited range in radius and to the simulation scaling, which impact the bulk Lorentz factor and energy of the particles (Sect. 3.5). The choice of energy range for integrating the emission can also play a role since pulses are narrower at higher energies (Fig. 8).

Fig. 13. Best fit (full line) based on the model skymaps to 6 pulse profiles (data points) taken from the Fermi-LAT 3rd pulsar catalog. The top row illustrates double-peaked pulse profile fits, including the Vela pulsar (top left). The bottom row illustrates single-peaked fits.

The bottom row illustrates some of the challenging variety seen in pulse profiles. In those cases, the best fit ended up being a single-pulse profile. The low-amplitude pulse of PSR J1513–5908 is consistent with a solution in which the observer line of sight grazes the emission region (α ≈ χ). PSR J1709–4429 illustrates a difficulty of the model when there are two peaks close in phase. The expectation is that the two peaks seen at α = 90° get closer in phase as the inclination α approaches χ. However, the skymaps in Figure 3 show that this does not happen because the amplitude of the peaks drops quickly away from α ≈ 90° for large obliquities and because the peaks are broad and smeared into a single pulse for low obliquities. Thus, the model tends to fit a single pulse. Last, PSR J2229+6114 illustrates a slow rise and fast decay profile seen in several pulsars that is difficult to understand. Indeed, the current sheet has a sharp density jump when it moves into the line-of-sight, with the density gradually falling off as it rotates away, i.e. a fast rise and slow decay would be more natural. The fits for PSR J1709–4429 and PSR J2229+6114 are both localised in a low-emission area of the skymap. The low radiative efficiency is unlikely to be compatible with the observed flux even though both are young pulsars with high Ė (3.5 × 10³⁶ and 2.2 × 10³⁷ erg s⁻¹, respectively). Notwithstanding the statistically poor fit, we therefore regard the solutions for PSR J1709–4429 and PSR J2229+6114 as unphysical.

4.3. Angle distributions

Fitting all the 3PC pulsar profiles, we found that the global distribution in α or χ differs significantly from a randomly oriented distribution of angles. We verified through a Monte Carlo simulation that this is not due to a bias in the population selection or in the fitting procedure. To do this, simulated light curves were constructed from the synthetic skymaps, with a random orientation for α and χ, and scaled so that their baseline flux and amplitude correspond to one of the 294 observed pulsars. We then randomly picked an observed flux for each phase by assuming Poisson statistics at the predicted count rate for each phase, before adding a random phase shift. Simulating 100 observations of the 294 3PC pulsars clearly showed that our fits to the 3PC population have an excess of solutions with low inclination α ≈ 20° and high obliquity χ ≈ 60°. We attribute these excesses to unphysical fits such as those described above, which bias the distribution.

To get a measure of the biases introduced by the model on the distribution of angles, we investigated fitting using model skymaps with different energy selections or using only electron or all particle emission. Interestingly, the variance of the best fit C² is lowest when using only emission from the positrons, all other combinations yielded poorer results. Although the fits are relatively consistent regardless of the skymap when an acceptable solution exists (e.g. the pulsars on the top row of Fig. 13), no robust trends emerge from the global distributions. Restricting to the best fits or to sub-populations of pulsars (such as radio loud vs. radio quiet or millisecond pulsars vs. young pulsars) did not yield more significant results. Directly fitting the model to the observed pulsar population is not yet accurate enough to provide useful results on the distribution of their inclination and obliquity. Previous works, which also systematically fitted pulsar profiles to various models, encountered the same difficulty as they showed that not a single model fares well on all types of pulsar profiles (Romani & Watters 2010; Venter et al. 2012; Johnson et al. 2014; Pierbattista et al. 2015)

4.4. Number of peaks

To circumvent this limitation in directly fitting data, Pierbattista et al. (2016) compared the ability of various pulsar models to reproduce major morphological characteristics of the pulse profile such as the number of peaks. Similarly, we estimated from the skymaps that ≳62% of the observed pulse profiles should be double-peaked, assuming that α and χ are randomly oriented and that all pulsars are equally detectable. The 3PC catalog lists 74% pulsars with two or more peaks. Double peaks are associated with stronger emission (Fig. 3), so the detection bias is consistent with observing a higher fraction than expected. Pierbattista et al. (2016) pointed out a higher fraction of single-pulsed profiles in radio loud pulsars (≈30%) compared to radio quiet pulsars (≈10%) based on the sample young and middle-aged pulsars in the second Fermi-LAT pulsar catalog. Radio quiet pulsars are more likely to be seen edge-on since radio emission is associated with the polar caps, hence more likely to produce double-peaked pulse profiles. Alternatively, a more aligned population (χ → 0°) would be less likely to produce double-peaked profiles. In the 3PC catalog there are 93 pulsars with a spindown timescale ≲250 kyr and with a listed number of peaks : 51 are detected in radio and 42 are radio-quiet. Of those, 20/51 = 39% of the radio-loud and 11/42 = 26% of the radio-quiet have a single-pulse profile, which a Z-test shows is not a significant difference. In general, we did not find significant differences (> 2σ) in the number of peaks when comparing sub-populations of pulsars in the 3PC catalog.

5. Summary

To summarize, the salient features reported in the modeling of synthetic pulsar light curves are:

• Particle acceleration powered by reconnection in the wind current layer, leading to synchrotron and inverse Compton radiation.

• The synchrotron and inverse Compton emission patterns are similar, but not identical, and in-phase. They present two bright caustics leading, in general, to two main pulses per spin period.

• The separation in phase of both pulses is explained by the geometry of the reconnection layer forming in the pulsar wind and is well approximated by the split-monopole solution.

• The current layer is not uniformly bright at all latitudes, it is mostly shining in the equatorial plane. As a direct consequence, pulsars with nearly symmetric pulses separated by 0.5 in phase should present the highest apparent luminosities.

• An interpulse emission, or third peak, generically appears away from the equatorial regions at intermediate latitudes and it is most visible in the synchrotron light curves. This feature most convincingly corresponds to the volumetric return current densities regions and reflects the strong asymmetric nature of the wind between two consecutive layers. Additional features at high viewing angles are also consistent with super-Goldreich-Julian current densities regions. More works dedicated to fully characterize the origin of these extra components are required.

• The emission pattern has a central symmetry with respect to the center of the maps, located here in the equatorial plane at phase Φ = 0.5. This symmetry is related to the direction of the electric current carried by the layer, which leads to asymmetric pulse profiles seen by an observer looking away from the equatorial plane.

• The asymmetry between both pulses increases with the photon energy. In particular, Vela-like light curves with dominant P2 at high energies have been reproduced.

• The widths of the light-curve peaks are set by two competing effects: (i) the acceleration of the expanding wind implies relativistic beaming that narrows the pulse width and (ii) the geometric expansion of the layer which increases the pulse width. The emitted pulses are thinnest at a radius fully determined by the plasma multiplicity in the current layer, $r_{\mathrm{m}}=\sqrt{2\kappa }$.

• The inverse Compton pulses are about half narrower than their synchrotron counterpart, provided that the target photon field is uniform within r_m. Conversely, if the photon field is uniform beyond r_c = πΓ_LCκ, pulsations in the inverse Compton channel vanish.

• The synchrotron and inverse Compton light-curve peaks sharpen with increasing photon energy. The mean particle energy increases with radius so that more distant parts of the wind contribute to the high-energy emission, where beaming is more important and synchrotron cooling is weaker. This may also be due to kinetic beaming effects in the reconnection layer (the higher the particle energy, the deeper and the more focused its trajectory ends up in the layer).

• There is no notable differences between the synchrotron and inverse Compton light curves if the target photon field originates from the neutron star surface. Anisotropic effects and the dilution of the target photon field diminish the radiated power. Anisotropic effects are strongest in the Thomson regime, leading to a steeper decay of the radiated power with the distance to the star.

• The reconnection rate sets an upper limit to the bolometric synchrotron radiative efficiency, ${\eta }_{\mathrm{rad}}^{\max }\lesssim {\beta }_{\mathrm{rec}}\sim 0.1$.

This study paves the way for comparisons with observations of synthetic pulse profiles constructed from ab initio PIC simulations. The synthetic pulse profiles reproduce some of the salient features of observed gamma-ray pulsars, including the generic double-peaked structure and the presence of a bridge or third peak in between the main pulses, with the pulse narrowing with increasing energy. Hence, fitting the pulse profiles of the Fermi-LAT 3rd pulsar catalog shows promising results despite the inherent limitations of a “one-size-fits-all” approach with a single (millisecond) pulsar model sampling a limited scale range and where the only parameters are the inclination and obliquity. Using full-scale hybrid simulations (Soudais et al. 2024) may be the way forward to match observations more closely.

Data availability

The atlas of the 3PC best fit light curves is available on Zenodo.

Acknowledgments

We thank Alexander Philippov, Alexander Chernoglazov, Christo Venter, and the referee Constantinos Kalapotharakos for useful comments regarding this work. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 863412). Computing resources were provided by TGCC under the allocation A0130407669 made by GENCI and PRACE (call 17).

References

All Figures

Fig. 1. Global plasma and magnetic structures forming in a pulsar magnetosphere of magnetic obliquity χ = 30°. Left panel: Three-dimensional rendering of the plasma density (color-coded) normalized by the surface Goldreich-Julian density and multiplied by (r/r⋆)2 to compensate for the effect of expansion. Closed and open field lines encompassing the equatorial current layer are shown by blue lines. The star is the gray sphere in the center. Right panels: Poloidal (top, y = 0) and toroidal (bottom, z = 0 or θ = 90°) slice of the plasma density. The light cylinder is shown with the white dashed line. The dotted lines represent the split-monopole solution, and the magenta arrow shows the orientation of the magnetic moment of the star at this time.

In the text

	Fig. 2. Ratio of the effective perpendicular field, ${\stackrel{\sim }{B}}_{\perp }$, to the total field strength, B, in the pulsar midplane (θ = 90°) for the χ = 60° simulation. Blue streamlines represent the magnetic field lines in the plane, and the red dashed circle shows the location of the light cylinder.
In the text

	Fig. 3. Synthetic synchrotron (left panels) and inverse Compton (isotropic radiation field, right panels) fluxes as a function of the observer’s viewing angle, α, and pulsar phase, Φ. Each map is normalized by its maximum value and includes the emission from both electrons and positrons. The white dashed line represents the location of the current sheet in the split monopole geometry.
In the text

	Fig. 4. Synthetic phase-resolved synchrotron (top) and inverse Compton (isotropic radiation field, bottom) fluxes (in arbitrary units) for χ = 30° color-coded by the observer’s viewing angle (emission from positrons only).
In the text

	Fig. 5. Beam correction factor, fΩ, inferred from the synchrotron (top) and inverse Compton (isotropic radiation field, bottom) skymaps for all the magnetic obliquities simulated in this work.
In the text

	Fig. 6. Spatial distribution of the synchrotron power per unit of volume, d𝒫syn/dV, emitted by the positrons only at the light cylinder for χ = 60°. This figure highlights the role of the volumetric return current in explaining the Vela-like interpulse.
In the text

Fig. 7. Schematic representation illustrating how pulses form relying on simple geometric considerations based on the split-monopole configuration. Here, χ = 60° for illustrative purposes, but this picture applies for all inclinations (except the aligned and orthogonal rotators). In this diagram, the pulse P1 (P2) is emitted once the sheet 1 (respectively 2) rotating at the angular velocity of the star, Ω, overlaps with the fixed black dashed line, an iso-phase contour set by the observer. Three typical configurations are presented: (i) an equatorial view (α = π/2) leads to two symmetric pulses separated by ΔΦ = 0.5, (ii) an off-equator view within the striped wind region (α = π/2 − 0.8χ) where the pulses are asymmetric and separated by ΔΦ = 0.28 and contains significant emission in between (a bridge or third peak, green region), and (iii) the degenerate single pulse case (ΔΦ = 0) where the observer’s line of sight grazes the outer boundary of the stripes (α = π/2 − χ).

In the text

Fig. 8. Energy dependence of the synthetic synchrotron pulse profile with the photon frequency, ν normalized to ν0, for a Vela-like solution (χ = 60°, α = 110°). Top panel: Light curves equally spaced in log ν/ν0 from the lowest frequencies (bottom) to the highest frequencies (top). The bottom panel shows the full phase-frequency emission map. In both panels, light curves are normalized to the maximum flux for each frequency band. This figure only shows the contribution from positrons.

In the text

Fig. 9. Radial distribution of the frequency-integrated radiated power in the synchrotron channel (top panel), in the inverse Compton channel for both an isotropic (middle panel) and an anisotropic (from the star, bottom panel) target photon field in the Thomson (red lines) and deep Klein-Nishina (blue lines) regimes. Simulations with an initial dipolar field are shown by solid lines (χ = 30°, 45°, 60°), while the split-monopole simulation is shown by dashed lines. Simple analytical expressions are also overplotted for comparison (black dashed-dotted and dotted lines). All quantities are normalized by their value at the light cylinder.

In the text

	Fig. 10. Radial evolution of all energy channels in the split-monopole solution, including the synchrotron luminosity, normalized to the pulsar spindown power (L0). The inferred reconnection rate is βrec ≈ 0.12 and the radiative efficiency is ηrad ≈ 0.05. The inverse Compton power is not shown since radiative energy losses through this channel are neglected in this work.
In the text

	Fig. 11. Narrowing of the pulse full width at half maximum with the energy of the emitted photons in the inverse Compton (isotropic radiation field, red line) and synchrotron channels (blue line) for the dipole simulation with χ = 60° and α = 110°.
In the text

Fig. 12. Evolution of the observed light curve as a function of the radius of emission. Top panel: Synthetic synchrotron and inverse Compton (isotropic radiation field) light curves as a function of the emitting radius, from the light-cylinder radius up to 5LC for the dipole simulation (χ = 60°, α = 90°). Bottom panel: Full width at half maximum of the inverse Compton pulse width as a function of radius for the dipole (magenta crosses) and split-monopole (blue dots) simulations. The red dotted and dashed lines shows the two asymptotic behaviors discussed here. The black vertical dashed line marks the location of the critical radius rm where the emitted pulses are the narrowest.

In the text

	Fig. 13. Best fit (full line) based on the model skymaps to 6 pulse profiles (data points) taken from the Fermi-LAT 3rd pulsar catalog. The top row illustrates double-peaked pulse profile fits, including the Vela pulsar (top left). The bottom row illustrates single-peaked fits.
In the text