Very-high-energy Emission from Pulsars

The magnetic fields in the FF magnetosphere for various magnetic inclination angles, α, are generated in Cartesian grids with resolution 0.02 R_LC, from the NS "surface" at 0.2 R_LC to 2.0 R_LC. To compute particle trajectories and radiation, the local magnetic field values are obtained by 3D interpolation of the Cartesian grid. Since the actual NS surface is much smaller than the lower limit of the grid, we extend the field values to the real pulsar surface by a ramp function, joining the numerical FF solution with an analytic retarded vacuum dipole (Deutsch 1955; Dyks et al. 2004).

To compute the particle trajectories, we first determine the PC rims by finding the last open field lines at each magnetic azimuth by fourth-order Runge–Kutta integration (Dyks et al. 2004) and the open volume coordinates (ovc) that are the radial, r_ovc, and azimuthal, l_ovc, divisions of the distorted PCs. Particles are then injected at the NS surface at particular (r_ovc, l_ovc) points, where r_ovc = 0 labels the magnetic axis and r_ovc = 1 labels the PC rim. For the calculations in this paper we use between 360 and 720 azimuthal divisions and inject pairs from r_ovc = 0.8 to r_ovc = 0.9 and primaries from r_ovc = 0.9 to r_ovc = 0.96. For the Crab, the pairs are injected between r_ovc = 0.88 and r_ovc = 0.91 and the primaries between r_ovc = 0.91 and r_ovc = 0.94. We have observed in our PIC models (Kalapotharakos et al. 2018) that the width of the acceleration region in the current sheet becomes smaller with increasing pair injection. Since very energetic pulsars like the Crab have larger pair multiplicity, narrower widths over which particle are injected and accelerated are expected. Primary particles are injected with low Lorentz factors γ_p = 200 and pairs are injected with a spectrum of Lorentz factors (see Section 2.2).

We compute the particle trajectories in the IOF as in HK15, with the velocity being a sum of a component parallel to the magnetic field and the local drift velocity. The dominance of the drift velocity near and beyond the light cylinder allows the particle velocity to remain subluminal while drifting backward (in phase) along field lines. However, in the present calculation the particle trajectories and their radii of curvature, assumed to be independent of particle energy, are first computed using interpolation and smoothing, with fine spatial resolution, and stored in tables. The subsequent calculation of the particle dynamics and radiation at different Lorentz factors, which use dynamic step sizes, then interpolates from these precomputed tables. This method of computing trajectories, detailed in Barnard et al. (2021), gives a more accurate determination of the radii of curvature of the trajectories.

Electron–positron pair cascades at the PCs are necessary to supply plasma to the magnetosphere, since extraction of charges of a single sign induces large accelerating electric fields, E_∥, and cannot produce the current required by the global solution (Timokhin 2010; Timokhin & Arons 2013) over the whole PC. Pairs are also necessary to screen E_∥ throughout most of the magnetosphere to produce the near-FF global conditions. The cascades are initiated by high-energy photons that are produced by particles accelerated in the induced E_∥, creating pairs by the one-photon process in the strong magnetic field. The pairs produce SR photons that create multiple (up to five) generations of pairs (Daugherty & Harding 1982). Timokhin (2010) and Timokhin & Arons (2013) found that the need to supply the global, near-FF current as a boundary condition produces time-dependent pair cascades, where pairs are created in bursts that completely screen E_∥, which reappears after the pairs exit the gap. However, the multiplicity of pairs created to screen the gap is small in comparison to the multiplicity of the full cascade, most of which occurs above the gap (Timokhin & Harding 2015). The region of particle acceleration and screening is also much smaller than the typical distance over which the full cascade extends. We can therefore calculate the spectrum of pairs from the full cascade using the energy of the gap-accelerated particles.

To generate the pair spectra for the calculations of this paper, we use the estimate of gap particle Lorentz factor from Equation (23) of Timokhin & Harding (2019):

$$\begin{eqnarray}&&{\gamma }_{{\rm{gap}}}\simeq 3.2\times {10}^{7}\,{P}^{-1/7}\,{B}_{12}^{-1/7}\,{\rho }_{c,7}^{4/7},\end{eqnarray} \tag{ 1 }$$

which is similar to the expression derived by Harding & Muslimov (1998) for steady gaps, where P is the pulsar period in seconds, B₁₂ is the surface magnetic field strength in units of 10¹² G, and ρ_c,7 is the radius of curvature of the field line at the surface in units of 10⁷ cm. The major characteristic is that the gap particle energy is very insensitive to surface magnetic field and period but is most sensitive to the radius of curvature of the field line. Although there is a distribution of particles moving upward from the gap, the highest-energy particles that produce most of the pairs are nearly monoenergetic. The results of Timokhin & Harding (2015) adopted a 1D gap that did not take into account the variation of curvature radius across the PC, which causes a variation in accelerated particle energy (Harding & Muslimov 1998). We use the estimate of the gap particle energy above as the initial energy for the full cascade and compute the pair spectra from the whole PC, as in Harding & Muslimov (2011), varying ρ_c , which determines initial energy. Although we have an estimate for the initial cascade particle energy, the multiplicity, i.e., the ratio of the density of high-energy particles from the gap to the Goldreich–Julian density, ρ_GJ = −BΩ/π ec, is not as well determined. We will therefore take the pair cascade multiplicity M₊ (number of pairs created via the seed emission of a single primary) as a variable parameter. The pair spectrum, however, is determined by the initial particle energy and the magnetic field geometry.

Figure 1 shows the pair spectra for the sources we model in this paper. The different pulsars have a variety of spectral shapes, with very energetic pulsars like the Crab displaying the largest range of pair energies, extending from Lorentz factors of 20 to above 10⁶. Middle-aged pulsars like Vela and B1706−44 have pair spectra that cut off at somewhat lower energies, and older pulsars like Geminga have spectra with an even more limited energy range. Millisecond pulsars (MSPs) such as J0218+4232 have pair spectra with much higher low- and high-energy cutoffs, ranging from γ ∼ 5000 to nearly 10⁷. Since they have much lower surface field strengths than young pulsars, the photons must have much higher energies to reach the threshold for pair production, leading to higher pair energies. The spectra with solid lines in Figure 1 for Vela, B1706−44, Geminga, and J0218+4232 and the dashed-line spectrum for the Crab are products of pair cascade models that assume a pure dipole magnetic field. The dashed-line spectrum for Vela and the solid-line spectrum for the Crab are alternative models we explore. The dashed-line spectrum for Vela is a product of a pair cascade model that assumes a nondipolar field structure in which a toroidal perturbation has been added to a dipole (Harding & Muslimov 2011), parameterized by the quantity . The Vela pair spectrum shown in Figure 1 assumes = 0.4. The alternative Crab pair spectrum assumes a power-law extension of the dipole pair spectrum, which would be necessary to account for the COMPTEL soft γ-ray points (see Section 3.2).

Zoom In Zoom Out Reset image size

It was found that the pair cascades that match the global currents can only operate in PC regions where the current J/J_GJ > 1 or J/J_GJ < 0 (Timokhin & Arons 2013), where J_GJ = ρ_GJ c. Therefore, in addition to the pair injection regions that are limited by the range of r_ovc, we further limit the pair injection to the regions where the global return current J/J_GJ < 0 connects to the current sheet (by limiting injection to a range of l_ovc = 1.5–4.78 rad). For oblique rotators, this region is located on field lines oriented toward the spin equator.

Our modeling of the broadband radiation of the electron–positron pairs and accelerated particles (primaries) follows closely that of HK15, and details of the calculations can be found there. There are a number of differences in the present calculation that we will detail below. The primaries and pairs are injected at the NS surface as described in previous sections. The primaries are accelerated by E_∥ fields along their trajectories with assumed values that are lower below the light cylinder and higher in the current sheet. This two-tier E_∥ is guided by PIC models that show the highest acceleration in the current sheet starting at the so-called Y-point (i.e., the tip of the zone of closed field lines) and a low level of acceleration below the light cylinder (Brambilla et al. 2018). The pairs are injected with a range of energies and relative fluxes equal to those of the pair spectra shown in Figure 1 and are not accelerated along their trajectories. As in HK15, the code loops through all particle trajectories twice, the first time to compute the SC radiation and 3D grid of emissivities and the second time to compute the ICS emission using the SC emissivities from the first loop. Also as in HK15, we assume that both pairs and primaries are injected with zero pitch angle but can acquire finite pitch angles through resonant absorption of radio photons that are emitted at a specified altitude, r_radio, above the NS surface. However, the observed phase alignment of the light curves from radio to TeV γ-rays for the Crab implies that the radio emission also comes from the current sheet (Philippov et al. 2019). The cone beam geometry used in the radio absorption calculation of HK15 cannot be applied outside the light cylinder. Therefore, in our model for the Crab, to simulate absorption of radio emission in the current sheet, we give the particles a constant pitch angle above a radius of 0.8 R_LC, with pairs maintaining a pitch angle ψ = 2 × 10⁻³ and primaries ψ = 10⁻⁵.

The primaries will therefore radiate a mixture of SR and CR that is captured by treating the radiation as SC, with the emission below the light cylinder being mostly SR, and becoming mostly CR in the current sheet as particle Lorentz factors exceed 10⁷. Since the pairs have relatively low Lorentz factors, their radiation in the first loop is in the SR limit of SC emission. We use the SC expression for the radiated power per photon energy interval presented in Torres (2018):

$$\begin{eqnarray}&&\displaystyle \frac{{{dP}}_{\mathrm{SC}}}{{dE}}=\displaystyle \frac{\sqrt{3}{e}^{2}\gamma y}{4\pi {\rm{\hslash }}{r}_{\mathrm{eff}}}\,\left[(1+z)\,F(y)-(1-z)\,{K}_{2/3}(y)\right],\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&y=\varepsilon /{\varepsilon }_{c},\,\,\,\,\,z={({Q}_{2}\,{r}_{\mathrm{eff}})}^{-2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\varepsilon }_{c}=\displaystyle \frac{3}{2}{\rm{\hslash }}{{cQ}}_{2}{\gamma }^{3}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&F(y)={\int }_{y}^{\infty }{K}_{5/3}(y){dy},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{r}_{\mathrm{eff}}=\displaystyle \frac{{\rho }_{c}}{{\cos }^{2}\psi }\,{\left[1+\xi +\displaystyle \frac{{r}_{g}}{{\rho }_{c}}\right]}^{-1}\end{eqnarray} \tag{ 6 }$$

is the effective curvature radius,

$$\begin{eqnarray}&&{r}_{g}=\displaystyle \frac{{{mc}}^{2}\gamma \,\sin \psi }{{eB}}s\end{eqnarray} \tag{ 7 }$$

is the gyroradius, ψ is the pitch angle,

$$\begin{eqnarray}&&{Q}_{2}=\displaystyle \frac{{\cos }^{2}\psi }{{\rho }_{c}}\,{\left(1+3\xi +{\xi }^{2}+\displaystyle \frac{{r}_{g}}{{\rho }_{c}}\right)}^{1/2},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{\rho }_{c}{\sin }^{2}\psi }{{r}_{g}{\cos }^{2}\psi }.\end{eqnarray} \tag{ 9 }$$

The parameter ξ determines the dominant SC regime, with SR dominating when ξ ≪ 1 and CR when ξ ≫ 1.

To compute the ICS emission of both pairs and primaries, the photon density at each position of a particle along its trajectory is calculated from the stored SC emissivities from the first loop, as detailed in HK15. The particle dynamics follows the Lorentz factor and perpendicular momentum, p_⊥, along each trajectory, integrating Equations (39) and (40) of HK15, with the only difference being the replacement of separate SR and CR loss rates with the SC loss rate

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma }{{dt}}=\displaystyle \frac{2{e}^{2}\,{\gamma }^{4}\,c}{3{\rho }_{c}^{2}{{mc}}^{2}}\,{g}_{r},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{g}_{r}={\left(\displaystyle \frac{{\rho }_{c}}{{r}_{\mathrm{eff}}}\right)}^{2}\,\displaystyle \frac{[1+7{({r}_{\mathrm{eff}}{Q}_{2})}^{-2}]}{8{({Q}_{2}{r}_{\mathrm{eff}})}^{-1}}\end{eqnarray} \tag{ 11 }$$

is a factor that determines the dominance of CR (g_r ∼ 1) or SR losses (g_r ≫ 1).

In the first trajectory loop, the ICS losses are zero so only SC losses and E_∥ acceleration (for primaries) are at play. In the second trajectory loop, SC losses are included in the particle dynamics since they dominate over ICS losses, justifying the neglect of the latter in the first loop. Also in the second loop, we do not need to compute the perpendicular particle momentum p_⊥ along the trajectory, since the ICS emission depends only on γ.

In this calculation we introduce attenuation of the VHE emission due to γ–γ pair production. Our approximate treatment is to multiply the ICS flux by an attenuation factor

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{IC}}={\dot{N}}_{\mathrm{IC},0}\,\exp [-{\tau }_{\gamma \gamma }(\varepsilon )]\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }={\int }_{0}^{s}\,{n}_{\gamma }\,{\sigma }_{\gamma \gamma }\,{ds}^{\prime} \sim {n}_{\gamma }\,\displaystyle \frac{{\sigma }_{{\rm{T}}}}{3}\,{\rm{\Delta }}s,\end{eqnarray} \tag{ 13 }$$

with σ_T the Thomson cross section. The VHE photons of energy ε scatter primarily with photons of energy

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{s}}}=\displaystyle \frac{2}{{\varepsilon }_{\gamma }(1-\cos \theta )}\simeq \displaystyle \frac{4{\gamma }_{{\rm{x}}}^{2}}{{\varepsilon }_{\gamma }},\end{eqnarray} \tag{ 14 }$$

since θ ≃ 1/γ_x, where γ_x is the Lorentz factor of the optical/IR emitting pairs, which we set to 10. In Equation (13), n_γ is the soft-photon density, which we approximate as

$$\begin{eqnarray}&&{n}_{\gamma }\simeq \displaystyle \frac{{L}_{s}}{{\varepsilon }_{{\rm{s}}}{Ac}}\end{eqnarray} \tag{ 15 }$$

where L_s is the luminosity at energy ε_s and $A\simeq 0.1\,{R}_{\mathrm{LC}}^{2}$ is the approximate area. Also in Equation (13) we set Δs ∼ 2R_LC − r_p, where r_p is the position of the ICS-emitting particle, so that Δs is the path length that the VHE photon will travel to exit the computation box.

By estimating the mean free path ℓ_{γ γ} ≃ 1/n_γ σ_{γ γ} , the Crab has ℓ_{γ γ} ∼ 4 × 10⁶–1.5 × 10⁷ cm < R_LC, while Vela has ℓ_{γ γ} ∼10¹⁰ cm ≫ R_LC. The other sources have mean free paths similar to or larger than Vela's. Therefore we find that γ–γ attenuation will only be important for the Crab.

The major differences, then, between the model in this paper and that of HK15 are as follows. (1) The SC radiation spectrum and energy losses are implemented for primary particles. (2) The SC emissivity does not cover the full computational grid but has spatial upper limits that are tailored to the extent of the soft emission for each pulsar to provide better resolution. For Vela, Geminga, and B1706–44 the emissivity grid extends to 0.8 R_LC, while for the Crab and J0218+4232 it extends to 2.0 R_LC. (3) The spectral energy range has been extended to cover 18 decades, from IR to 100 TeV. (4) The computation of particle trajectories and their radii of curvature is more accurate. (5) E_∥ now has different constant values inside and outside the light cylinder, whereas HK15 assumed the same value from the NS surface to 2 R_LC. (6) The pair injection occurs only over part of the PC, where J/J_GJ < 0.

All of the above changes were implemented in H18. The model in this paper has the additional changes that were not in HK18. (1) We corrected an error of a factor ∼2 in the ICS radiation, where the energies in units of mc² were not converted to MeV. (2) We updated the radio luminosity values for Vela. (3) γ-γ pair attenuation of the ICS radiation is now included.

Table 1 lists both the observed and the model parameters for the sources we study in this paper. All of the radio fluxes listed are those measured at 400 MHz, except for Geminga, which has no measured radio emission. For Geminga we therefore assume that it produces radio emission that is not beamed in our direction and assume a value of 1000 mJy at 400 MHz to model the pair SR. The parameters that are varied to match the observed spectral energy distributions (SEDs) and high-energy light curves are the magnetic inclination angle, α, the observer angle, ζ, the low and high values of accelerating electric field, ${R}_{\mathrm{acc}}^{\mathrm{low}}={{eE}}_{\parallel }^{\mathrm{low}}/{{mc}}^{2}$ and ${R}_{\mathrm{acc}}^{\mathrm{high}}={{eE}}_{\parallel }^{\mathrm{high}}/{{mc}}^{2}$, the primary current, J/J_GJ, the pair multiplicity, M₊ and the radio emission altitude, r_radio. We do not consider the viewing angle ζ as a free parameter since, for a given α, it is adjusted to match the Fermi light curves.

Table 1. Source and Model Parameters

Pulsar	P	d	Radio Flux	α	${R}_{\mathrm{acc}}^{\mathrm{low}}$	${R}_{\mathrm{acc}}^{\mathrm{high}}$	J/JGJ	M+	rradio
	(s)	(kpc)	(mJy)		(cm−1)	(cm−1)			(RLC)
Vela	0.089	0.25	5000	75°	0.04	0.2	18	6 × 103	0.1/0.2
Crab	0.033	2.0	700	45°	0.04	0.4	5.0	3 × 105	>0.8
B1706−44	0.102	2.3	25	45°/30°	0.04	0.2	20	6 × 104	0.08
Geminga	0.237	0.25	1000	75°	0.04	0.15	10	2 × 104	0.1
J0218+4232	0.0023	3.1	100	60°	0.5	5.0	150	3 × 105	0.38

Download table as:  ASCIITypeset image

Figure 2 shows a model SED for the Vela pulsar using the solid pair spectrum in Figure 1, for three values of the viewing angle. This result is very similar to the model in H18, which used the same pair spectrum for a pure dipole and the same ${E}_{\parallel }^{\mathrm{low}}$, ${E}_{\parallel }^{\mathrm{high}}$, and M₊ values. While the pair SR spectra are sensitive to viewing angle, the high-energy primary SC and ICS spectra are much less sensitive. The primary ICS flux depends primarily on the low-energy part of the pair SR spectrum, which the high-energy primaries can scatter in the Thompson limit, while most of the sensitivity to ζ in the pair SR spectra is at the highest energies. In the soft X-ray range, the observed Vela spectrum is dominated by a strong thermal component that we have not plotted. There is a mismatch of the primary SC component with both the RXTE hard X-ray spectrum and the low-energy part of the Fermi spectrum around 100 MeV. This is a consequence of our simplified two-valued E_∥ distribution that creates an artificial dip in the SC spectrum where the E_∥ transition occurs at the light cylinder. A continuously variable E_∥ distribution that started at lower values and had no sharp jump at the light cylinder would produce a higher level of hard X-ray emission and a smoother transition to the Fermi spectrum. Future modeling could use the E_∥ distribution from dissipative or PIC global models as a guide. The shape of the primary SC spectrum determines the shape of the primary ICS spectrum below the cutoff (where the photon energy equals the maximum particle energy). Consequently, filling in the dip in the SC spectrum by a variable E_∥ to match the low-energy Fermi spectrum would raise the flux of the primary ICS spectrum.

Zoom In Zoom Out Reset image size

Figure 3 shows a model SED using instead the dashed pair spectrum in Figure 1, for three values of the viewing angle. The result is similar to that shown in H18 but with notable differences. This pair spectrum assumes a nondipolar field but the same M₊. Since the nondipolar pair spectrum extends to lower pair energy, the pair SR spectrum also extends to lower photon energy and the IR/optical flux is higher. Consequentially, this model's primary ICS flux is higher by almost an order of magnitude and is closer to the H.E.S.S.-II threshold. The primary SC component is the same as shown in Figure 2 since the ${E}_{\parallel }^{\mathrm{low}}$ and ${E}_{\parallel }^{\mathrm{high}}$ values are the same. Viewing angles in the range 50°–70° make minor differences to the primary ICS component.

Zoom In Zoom Out Reset image size

Figure 4 shows a model SED for the same parameters and pair spectrum, but with a lower pair multiplicity, M₊ = 1 × 10³. In this case, the primary ICS flux is about a factor of ten lower and closer to that predicted by H18, even though the pair SR spectrum here has a higher optical/IR flux. Although the pair SR spectrum for the ζ = 46° case in Figure 3 matches the observed optical spectrum as well, the higher values of ζ agree better with ζ ∼ 64° measured from the fits to the torus image of the pulsar wind nebula (Ng & Romani 2008).

Zoom In Zoom Out Reset image size

In Figure 5 we show the emission sky maps (intensity versus ζ versus ϕ) in three different energy ranges, as indicated in the figure. Also shown are light curves for the same energy ranges, which are horizontal cuts through the sky maps at ζ = 70°. The sky maps display bright caustics in the emission pattern that are characteristic of radiation from the FF current sheet (Bai & Spitkovsky 2010; Contopoulos & Kalapotharakos 2010; Barnard et al. submitted). The caustic patterns result from the particle acceleration being confined near the current sheet and their trajectories becoming nearly radial outside the light cylinder in the IOF. The photons that they emit parallel to their motion then reach the observer at the same phase, producing a sharp pulse each time our line of sight crosses the current sheet. Crossings of the current sheet can occur once, twice, or not at all, depending on the viewing angle. In these models (just like in two-pole-caustic/slot-gap models), light curves with two peaks are the result of emission from opposite rotational hemispheres, in contrast to traditional OG models (Romani & Yadigaroglu 1995) where the emission from both peaks originates from the same hemisphere. The particular pattern that appears in Figure 5 results from our assumption that E_∥ is constant in both radius and azimuth in the current sheet. Somewhat different emission patterns can appear for cases where E_∥ varies with radius and azimuth in dissipative magnetosphere (Kalapotharakos et al. 2014) or PIC (Kalapotharakos et al. 2018) models, for example. The light curves in those models provide a better match to the observed phases of the γ-ray peaks.

Zoom In Zoom Out Reset image size

Even though we have assumed a uniform accelerating field in the current sheet, the intensity of emission varies along the caustics and the intensity pattern varies with energy. The emission is more uniform along the caustic in the lowest energy range than in the higher energy ranges, where the highest-energy emission becomes more confined to particular locations. These changing patterns of emission produce peaks of similar flux at lower energy but increasingly unequal flux with increasing energy, with the first peak disappearing at the highest (TeV) energies (i.e., a continuation of the P1/P2 effect seen in the GeV band). The caustics also narrow with increasing energy, producing a narrowing of the peaks. Both of these effects—a harder second peak and a narrowing of both peaks with energy—are observed in many Fermi pulsar light curves. Similar results are presented by Barnard et al. (2021). Emission primarily in the second peak at energies above 10 GeV is observed for Vela, as well as Geminga and B1706−44 (see Sections 3.3 and 3.4).

Since we are also assuming that the MeV–GeV emission is SC, which is primarily CR at the energies displayed in Figure 5, the variation in intensity is caused by variations in local radius of curvature, as found by Barnard et al. (2021). This is confirmed in Figure 6, which shows a sky map of the maximum radius of curvature, ${\rho }_{{\rm{c}}}^{\max }$, of the trajectories of particles that emit photons in each direction ζ and ϕ. As in the sky maps of intensity, particles on a range of different field lines and radii can emit photons in the same direction. First of all, we notice that the sky pattern of ${\rho }_{{\rm{c}}}^{\max }$ closely reflects the sky pattern of caustic emission. This indicates that the largest radii of curvature occur in and near the current sheet at radii outside the light cylinder where the main caustics form. Also observe that the ${\rho }_{{\rm{c}}}^{\max }$ pattern is narrower than the emission caustics, indicating that the largest ρ_c values are more confined to the current sheet. The bottom panel in Figure 6 shows a cut at 70° through the sky map above to reveal that ${\rho }_{{\rm{c}}}^{\max }$ in the second light-curve peak (P2) is about three times higher than that in the first peak (P1). This figure then explains both the decreasing P1/P2 ratio and the narrowing of the peak with increasing energy as a characteristic of CR in an FF field geometry, where the spectral cutoff of CR is predicted to be

$$\begin{eqnarray}&&{E}_{\mathrm{CR}}\propto {{E}_{\parallel }}^{3/4}\,{\rho }_{{\rm{c}}}^{1/2}.\end{eqnarray} \tag{ 16 }$$

The above expression holds when the emitting particles are at the radiation–reaction limit where energy losses balance acceleration gains. Figure 3 of HK15 shows that accelerating particles in the current sheet have reached the curvature radiation–reaction (CRR) regime. The higher cutoff energy for P2 is due to the higher energy of the emitting particles in this regime,

$$\begin{eqnarray}&&{\gamma }_{\mathrm{CRR}}={\left(\displaystyle \frac{3\,{E}_{\parallel }\,{\rho }_{{\rm{c}}}^{2}}{2e}\right)}^{1/4},\end{eqnarray} \tag{ 17 }$$

given the larger value of ρ_c in this case.

Zoom In Zoom Out Reset image size

Figure 7 shows phase-resolved spectra for P1 and P2, illustrating that both the SC and ICS spectra of P2 are significantly harder. It is apparent that the P2 spectrum has a cutoff at higher energy, and this is reflected in the primary ICS spectrum. Since the high-energy cutoff of the P2 SC spectrum is due to the higher energy of the particles, these same particles are radiating an ICS spectrum that extends to higher energy.

Zoom In Zoom Out Reset image size

The decreasing P1/P2 ratio with increasing energy is dependent on α and ζ. For example, in Figure 6 the value of ρ_c in P1 can be larger, smaller, or similar to ρ_c in P2, for different ζ values. By symmetry, the ratio of ρ_c for P1 and P2 approaches 1 for ζ closer to 90°. For smaller α, the energy dependence of P1/P2 will be smaller for a wider range of ζ (see Section 3.2).

The model SED for the Crab pulsar is shown in Figure 8 for the pair spectrum with a power-law extension (solid line in Figure 1) and two different viewing angles. The pair SR component matches the optical through hard X-ray data reasonably well and requires a large pair multiplicity that is both suggested by nebular models and is theoretically possible in PC cascade models (Timokhin & Harding 2015). We find that the pair energy spectrum with a power-law extension, combined with the constant pair pitch angle, can account for the COMPTEL soft γ-ray points plus the lower part of the Fermi spectrum and produces a pair SSC spectrum that matches the MAGIC VHE points. The primary SC spectrum, combined with the high-energy part of the pair SR spectrum, can match the Fermi points up to about 50 GeV. The flux of the pair SSC component exceeds that of the primary SC component above 50 GeV. In this model, the MAGIC spectrum is not an extension of the primary SC spectrum, predicting a possible dip or discontinuity between the Fermi and MAGIC spectra. Our model also predicts a primary γ–γ attenuated ICS component peaking around 10–20 TeV that is potentially detectable by the High Altitude Water Cherenkov Observatory (HAWC). However, HAWC is primarily designed to operate in imaging mode and has not yet implemented searches for pulsed emission (B. Dingus, private communication). There are several MAGIC points and upper limits at the level our primary ICS spectral prediction, implying that more observation time by MAGIC may be able to detect part of this component. The pair SR spectrum for ζ = 60° matches the optical and X-ray data slightly better than that for ζ = 72°, but the pair SSC spectrum for ζ = 72° matches the MAGIC points better.

Zoom In Zoom Out Reset image size

Figure 9 shows the emission sky maps and light curves for ζ = 72° and for different energy ranges—one in the optical range and three at γ-ray energies. The model optical light curve shows two narrow peaks that are in phase with those of the γ-ray light curves, a result of locating all the emission in the high-altitude separatrix and current sheet. Even though the pairs and primaries are injected on different but adjacent sets of magnetic field lines at the NS surface, the radiation they emit near the current sheet produces light-curve peaks at the same phase. This is because the trajectories of particles outside the light cylinder become nearly radial. Radiation from these same sets of field lines inside the light cylinder would produce light curves at slightly different phases. As was also seen for our Vela model, the emission intensity distribution along the caustics changes with energy. However, for this value of α = 45° and ζ = 72° the relative intensity variation of P1/P2 is weaker. The decrease in P1/P2 flux with energy up to 50 GeV is slower than for Vela although we see that P1/P2 is smaller in the 10–50 GeV light curve. Indeed, the Fermi Crab light curves show a slower energy evolution than Vela.

Zoom In Zoom Out Reset image size

In Figure 10 we show model SEDs for the Geminga pulsar for α = 75°, ζ = 55°, and two different values of R_acc. The solid curve assumes two values inside and outside the light cylinder, ${R}_{\mathrm{acc}}^{\mathrm{low}}=0.04$ and ${R}_{\mathrm{acc}}^{\mathrm{high}}=0.15$, while the dashed spectrum assumes a single value ${R}_{\mathrm{acc}}^{\mathrm{low}}={R}_{\mathrm{acc}}^{\mathrm{high}}=0.15$ at all radii. Our assumed radio flux of 1000 mJy can account for the observed optical/UV emission through the pair SR component. The pair SR is very sensitive to ζ for ζ ≲ 50° since our line of sight must cut through enough of the low-altitude emission from the pairs, which is within about 30° of the magnetic pole at 75°. A viewing angle of ζ ∼ 50°–55° produces a flux of pair SR and shape of primary SC emission that best matches the data. The two E_∥ models produce similar primary SC spectra above about 1 GeV, where the CR regime dominates, but differ markedly below 1 GeV, where SC becomes a mix of SR and CR. The model with the lower E_∥ below the light cylinder suppresses acceleration at lower altitudes so that the primaries can temporarily maintain larger pitch angles and radiate a larger amount of SR. Since the hard X-ray spectrum of Geminga is not well measured, the ${E}_{\parallel }^{\mathrm{low}}$ value is currently unconstrained. We see that the MAGIC spectrum that is an extension of the Fermi spectrum is well explained as primary SC emission with no need to invoke any ICS component. Indeed, our predicted ICS emission, from both pairs (whose flux is below the level of the plot) and primaries, is well below any present detection thresholds. However, the two different models for E_∥ predict quite different primary ICS spectra, with the multivalued E_∥ model ICS spectrum extending to lower energies and having a slightly larger flux. Such an effect was also apparent in the Vela models presented by H18, where the observed hard X-ray flux constrains ${R}_{\mathrm{acc}}^{\mathrm{low}}=0.04$.

Zoom In Zoom Out Reset image size

As is apparent from Figure 1, the computed pair spectrum of Geminga is much lower than that of Vela, in total multiplicity and also in both low- and high-energy extent. Thus the predicted pair SR spectrum for Geminga has a lower flux and SED energy peak, both by about two orders of magnitude, and a narrower energy range than Vela. The distance of the light cylinder is also three times larger than for Vela, so that the high-energy particles in the current sheet are farther from the soft-photon source, making their local density lower. So even though the primary SC flux is comparable to Vela's, the lower soft-photon density produces a lower primary ICS component. The Geminga spectrum at soft X-ray energies is also dominated by a strong thermal component, which is not shown.

Figure 11 shows the sky maps and light curves in three different energy ranges: the lowest one in the Fermi range, the middle one for the MAGIC range, and the third for TeV energies. As in the case of Vela, for which we assumed the same α value, there is a strong dependence of the intensity distribution in the sky map on energy. This produces a strong decrease in the P1/P2 ratio with increasing energy, such that P1 should disappear above ∼20 GeV. Indeed, this is the case for both Fermi and MAGIC light curves, where P1 is not detected above 10 GeV, although the observed narrowing of the peaks is not as pronounced as indicated by our prediction.

Zoom In Zoom Out Reset image size

Figure 12 shows two SED models for PSR B1706−44, a Vela-like pulsar: one with α = 45° and ζ = 53° chosen to match the Fermi light curve and another with α = 30° and ζ = 60° that better matches the Chandra spectrum. The Fermi light curve is very different from that of Vela, with a narrow spacing of the observed γ-ray peaks that is better modeled by smaller inclination angles where our line of sight can cut through the caustics in the current sheet at a grazing angle as shown in Figure 13. There is no pulsed optical emission detected for this pulsar, since its distance is much larger than Vela's or Geminga's but there is a soft power-law X-ray detection from Gotthelf et al. (2002), which we use to match the flux of the pair SR component. The predicted pair SR spectrum has a higher flux (by two orders of magnitude) and SED peak energy (by a factor of 10–100) than for Geminga, primarily because its spectrum of pairs is higher in multiplicity and extends to higher energy. Although the pair spectrum of B1706−44 is very similar to that of Vela, the pair SR spectral flux is somewhat lower since the distance is larger and the radio luminosity is lower. However, the intrinsic luminosity of both the pair SR and primary SC must be larger than Vela's given that the Fermi peak SED flux is less than a factor of ten lower. Therefore, the luminosity of the primary ICS is predicted to be higher than Vela's, and it is only its larger distance and relatively lower predicted optical/IR flux that lead to the flux level being below detection thresholds. However, we predict that the flux of the primary ICS component is a factor of ten above that for Geminga. The two different model geometries produce two very different primary SC spectra since changing α and ζ changes the radii of curvature of the particle trajectory, which changes the dynamics and spectrum of the SC radiation. In some cases, the spectrum is more SR-dominated at low energies while in others it is more CR-dominated. Therefore a change in geometry, as well as a difference in E_∥ distributions as in two Geminga models we presented, can produce different primary SC spectra. We note that the primary ICS spectrum for the α = 30° model (dashed line) reaches a bit higher in flux, reflecting the low-energy enhancements of the primary SC spectrum and particle energy distribution. H.E.S.S.-II has detected pulsed emission up to 70 GeV, which we can explain as primary SC emission, as for Geminga. The pair SSC component is very weak and is not predicted to produce an additional component at VHE energies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show the energy-dependent sky maps and light curves for B1706−44 in Figure 13. The first two energy ranges are chosen to match the Fermi and H.E.S.S.-II detection ranges. As is the case for the other sources, the caustics in the emission sky map narrow with increasing energy, and the P1/P2 ratio decreases with increasing energy for this model ζ. This behavior is in agreement with the H.E.S.S.-II 10–70 GeV light curves in which P2 dominates. However, the narrowing of the peaks is not so apparent, either in the Fermi (Abdo et al. 2013) or H.E.S.S.-II light curves (Spir-Jacob et al. 1908). In the highest energy band, >70 GeV, we predict that only P2 will be present.

Our final source model is for J0218+4232, an energetic MSP with a very hard X-ray spectrum that extends to at least the highest NuSTAR energy, 80 keV. Figure 14 shows the model SED for inclination angle α = 60° and viewing angle ζ = 65° using the pair spectrum shown in Figure 1. The pair SR is very different from that of the younger pulsars, reflecting the much more energetic pair spectrum. The SR SED peak is predicted to be near 10 MeV, which matches the hardness of the observed X-ray spectrum. As a result, the optical/IR flux is predicted to be much lower than that of the younger pulsars by at least several orders of magnitude. Consequentially, the pair SSC component peaks at a VHE energy around 50 GeV but at a flux level well below any current detection threshold. The primary SC component can account for the higher-energy Fermi points, but falls short of the lower energy ones. It is possible that a stronger SR regime of primary SC radiation, or a pair spectrum extending to higher energy and producing a pair SR component extending to around 1 GeV, could reach the Fermi spectrum. A telescope with sensitivity at 1–100 MeV energies could address this question.

Zoom In Zoom Out Reset image size

Our model does not predict a detectable primary ICS component. The SED peaks around 1 TeV, a decade lower than for the younger pulsars, due to suppression by Klein–Nishina effects, because none of the ICS occurs in the Thompson limit and the flux is also suppressed. This effect was seen by HK15, who modeled pair SSC for two other energetic MSPs—B1821−24 and B1937+21—whose pair spectra are very similar to that of J0218+4232. Our model therefore suggests that MSPs will not be good targets for VHE telescopes. As shown in the figure, MAGIC has data totaling 90 hr observing J0218+4232 at 20 GeV—20 TeV, with only upper limits, although Fermi detects pulsed emission up to 25 GeV (Acciari et al. 2020). The discrepancy around 20–25 GeV is due to the high cosmic-ray background subtraction affecting MAGIC and other ground-based air-Cherenkov detectors.