Very-High-Energy Emission From Pulsars

Air-Cherenkov telescopes have detected pulsations at energies above 50 GeV from a growing number of Fermi pulsars. These include the Crab, Vela, PSR B1706-44 and Geminga, with the first two having pulsed detections above 1 TeV. In some cases, there appears to be very-high-energy (VHE) emission that is an extension of the Fermi spectra to high energies, while in other cases, additional higher-energy spectral components that require a separate emission mechanism may be present. We present results of broad-band spectral modeling using global magnetosphere fields and multiple emission mechanisms that include synchro-curvature (SC) and inverse Compton scattered (ICS) radiation from accelerated particles (primaries) and synchrotron-self Compton (SSC) emission from lower-energy pairs. Our models predict three distinct VHE components: SC from primaries whose high-energy tail can extend to 100 GeV, SSC from pairs that can extend to several TeV and ICS from primary particles accelerated in the current sheet, scattering pair synchrotron radiation, that appears beyond 10 TeV. Our models suggest that H.E.S.S.-II and MAGIC have detected the high-energy tail of the primary SC component that produces the Fermi spectrum in Vela, Geminga and PSR B1706-44. We argue that the ICS component peaking above 10 TeV from Vela has been seen by H.E.S.S. Detection of this emission component from the Crab and other pulsars is possible with HAWC and CTA, and directly measures the maximum particle energy in pulsars.


INTRODUCTION
The stunning increase in quality high-energy (HE) pulsar data facilitated by the operation of the Fermi Large Area Telescope (LAT; Atwood et al. 2009) provided ample opportunity to scrutinise our understanding of the emission mechanisms and magnetospheric geometry responsible for pulsed gamma-ray emission from these enigmatic stellar objects.A generic feature of the GeV spectra observed by the Fermi LAT is that they cut off, usually sub-exponentially, in a narrow band around a few GeV.Moreover, a stacking analysis involving 115 Fermi -detected pulsars (excluding the Crab pulsar) found no emission above 50 GeV (McCann 2015), indicating that one should probably in general not expect detectable levels of very-high-energy (VHE) photons from pulsars.
Earlier pulsar models indeed predicted detectable TeV fluxes from these sources, subject to some modelling uncertainties and approximations.For example, simulations invoking the standard outer gap (OG) model yielded spectral components beyond 100 GeV for some pulsars when estimating the inverse Compton scattering (ICS) flux of primary electrons on synchrotron radiation (SR) by secondary pairs (Cheng et al. 1986a,b;Romani 1996).This resulted in a natural bump around a few TeV that reached 5% of the GeV flux (Hirotani 2001;Takata et al. 2006), depending on model assumptions.However, it was not clear if such components would survive magnetic one-photon or two-photon absorption, so their detectability remained uncertain.Most recently Rudak & Dyks (2017) used an OG model to predict ICS emission for Vela above 10 TeV from primary gap-accelerated particles scattering optical/IR emission from secondary pairs.Some polar cap (PC) or pair-starved PC models predicted curvature radiation (CR) spectral cutoffs of up to 100 GeV (Bulik et al. 2000;Harding et al. 2002a).Later such models (Harding et al. 2005;Venter & De Jager 2005;Frackowiak & Rudak 2005), a slot gap (SG) model (Harding et al. 2008) and an OG model (Tang et al. 2008) predicted CR cutoffs around ∼ 10 − 50 GeV (still prior to the launch of the Fermi LAT, though).One could conclude that this CR component may just fall short of the reach of ground-based Cherenkov detectors, the energy threshold of which exceeded 100 GeV.
Given the uncertain and low expectations for pulsed TeV spectral components by some models and observations, it was thus surprising that MAGIC announced the detection of pulsations from the Crab pulsar at energies up to ∼ 25 GeV (Aliu et al. 2008;Aleksić et al. 2011Aleksić et al. , 2012)), followed by the VERITAS detection of pulsed photons up to ∼ 400 GeV (Aliu et al. 2011), and finally, MAGIC's detection of pulsations up to 1.5 TeV (Ansoldi et al. 2016).This opened a new window into pulsar science.H.E.S.S.-II next detected pulsed emission from the Vela pulsar in the sub-20 GeV to 100 GeV range (Abdalla et al. 2018).The latest observations by H.E.S.S. unveil pulsed emission from this pulsar at a few TeV (H.E.S.S. Collaboration, in preparation).Additionally, H.E.S.S.-II detected pulsed emission from PSR B1706−44 in the sub-100 GeV energy range (Spir-Jacob et al. 2019).MAGIC recently detected pulsed emission from the Geminga pulsar between 15 GeV and 75 GeV at a significance of 6.3σ (Acciari et al. 2020).However, only the second light curve peak is visible at these energies.The measured spectrum is smoothly connected to that measured by the Fermi LAT, and may perhaps indicate a transition from a CR to ICS spectral component.
Some trends are becoming apparent.Perhaps the most striking is the continuation of the decrease in the ratio of first to second light curve peaks (Abdo et al. 2010a(Abdo et al. , 2013) ) toward higher energies: as the photon energy E γ is increased beyond several GeV, the intensity ratio of the first (P1) to second (P2) light curve peak decreases for Vela and Geminga.In the GeV band, we attribute this to the systematically larger local curvature radius, ρ c , of trajectories of particles responsible for the second peak, compared to the first, leading to a relatively larger spectral cutoff for P2 (Barnard et al., submitted).The highest-energy particles may also be responsible for this effect in the TeV band via ICS emission.Furthermore, one notices that the γ-ray peaks of Crab, Vela and Geminga remain at the same phase positions, that the inter-peak (bridge) emission evolves in the case of Vela, and that the peak widths decrease for Crab (Aliu et al. 2011), Vela (Abdo et al. 2010b) and Geminga (Abdo et al. 2010c) at increased photon energy.
In this paper, we attempt to explain the new subTeV phenomenology and also make predictions of additional VHE emission components at higher energy to guide further observations by ground-based facilities.We model various spectral components, taking into account several recent developments, such as force-free (FF; Kalapotharakos & Contopoulos 2009) magnetic structures, and also indications from the FF inside and dissipative outside (FIDO) models (Brambilla et al. 2015;Kalapotharakos et al. 2014Kalapotharakos et al. , 2017) ) and particle-in-cell simulations (PIC; Brambilla et al. 2018;Cerutti et al. 2016aCerutti et al. ,b, 2020;;Kalapotharakos et al. 2018;Philippov & Spitkovsky 2018) that predominant dissipation is happening in the current sheet, rather than interior to the light cylinder.The added data from the TeV window may play an important role in constraining particle energetics and break degeneracies that exist when only considering lower-energy data.
The structure of the rest of the paper is as follows.In Section 2, we give a brief description of our model and calculation of the broadband spectra.In Section 3, we present phase-averaged and phase-resolved spectra for several pulsars: Vela, Crab, Geminga, PSR B1706-44 and PSR J0218+4232.Our discussion and conclusions follow in Section 4.

DESCRIPTION OF MODEL
Much of the magnetic field geometry, particle trajectory and radiation calculations are the same as in Harding et al. (2018), (hereafter H18), which is an updated version of the model used in Harding & Kalapotharakos (2015) (hereafter HK15).Since the details may be found in those papers, we will give only a summary of the main elements of the model and detailed descriptions of only what is new to the present calculations.
There has been recent rapid progress in understanding the structure, current flow and particle acceleration in global pulsar magnetospheres.Simulations of FF (where E • B = 0 and E < B is assumed) and dissipative magnetospheres (Kalapotharakos et al. 2012;Li et al. 2012), where a macroscopic conductivity is assumed, have shown that reproduction of Fermi γ-ray pulsar light curves requires near-FF magnetospheres (Kalapotharakos et al. 2014(Kalapotharakos et al. , 2017)).The region of highest accelerating electric field components are found near the current sheet outside the light cylinder at a radius R LC = c/Ω, where Ω is the pulsar angular rotation rate and c is the speed of light, indicating that most of the particle acceleration takes place there.More recent kinetic or PIC simulations of global magnetospheres (Kalapotharakos et al. 2018;Philippov & Spitkovsky 2018) confirm this finding by showing that the highest-energy particles are found in the current sheet.In our calculations, we adopt the 3D magnetic field geometry of an FF magnetosphere (Kalapotharakos et al. 2012) in which we inject both particles that enter the current sheet and are accelerated to high energy (primaries), and electron-positron pairs assumed to come from pair cascades near the NS surface (see §2.2).All particles are injected at foot points of open field lines at the PC (see §2.1) and radiate synchro-curvature (SC) and inverse-Compton scattering (ICS) emission (see §2.3) along their trajectories in the inertial observer frame (IOF).The radiated photons from each type of particle, pairs and primaries, are collected in arrays of energy, observer angle ζ with respect to the rotation axis and rotation phase φ, for a fixed magnetic inclination angle α.Sky maps of the intensity distribution in ζ vs. φ can then be constructed in different energy ranges and from these we can make energy-dependent light curves and phase-resolved spectra.

Magnetosphere geometry
The FF magnetosphere magnetic fields 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 on the neutron star (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 grows 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 §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 that use dynamic step sizes then interpolate from these pre-computed tables.This method of computing trajectories, detailed in Barnard et al., in prep., gives a more accurate determination of the trajectory radii of curvature.

PC pair spectra
Electron-positron pair cascades at the PCs are necessary to supply plasma to the magnetosphere, since extraction of charges of a single sign induce 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 5) 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 non-steady pair cascades, where pairs are created in bursts that completely screen the 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 Eqn (23) of Timokhin & Harding (2019) which is similar to the expression derived by Harding & Muslimov (1998) for steady gaps, where P is the pulsar period in s., B 12 is the surface magnetic field strength in units of 10 12 G, and ρ c,7 is the field line radius of curvature at the surface in units of 10 7 cm.The major characteristic is that the gap particle energy is very insensitive to surface magnetic field and period but is most sensitive to field line radius of curvature.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 above of the gap particle energy 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, 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 6 .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.MSPs such as J0218+4232 have pair spectra with much higher low and high-energy cutoffs, ranging from γ ∼ 5000 to nearly 10 7 .Since they have much lower surface field strengths compared with 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 assumes 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 Vela dashed-line spectrum is a product of a pair cascade model that assumes a non-dipolar field structure in which a toroidal perturbation has been added to a dipole (Harding & Muslimov 2011), parametrized 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 §3.2).
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 r ovc range, we further limit the pair injection to the regions where the global return current J/J GJ < 0 and connect to the current sheet (by limiting injection to a range of l ovc = 1.5 − 4.78 radians).For oblique rotators, this region is located on field lines orientated toward the spin equator.

Calculation of broadband photon spectra
Our modeling of the broadband radiation of the electron-positron pairs and accelerated particles (primaries) follows closely to 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 closed field line zone) 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 −3 and primaries ψ = 10 −5 .
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 7 .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) dP where is the effective curvature radius, is the gyroradius, ψ is the pitch angle, and 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 particle position 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 Eqns (39) and (40) of HK15, with the only difference being the replacement of separate SR and CR loss rates with the SC loss rate where 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 with with σ T the Thomson cross section.The VHE photons of energy ε scatter primarily with photons of energy since θ 1/γ x , where γ x is the Lorentz factor of the Optical/IR emitting pairs which we set to 10.In Eqn ( 13), n γ is the soft-photon density which we approximate as where L s is the luminosity at energy ε s and A 0.1 R 2 LC is the approximate area.Also in Eqn ( 13) we set ∆s ∼ (2 R 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 6 − 1.5 × 10 7 cm < R LC , while Vela has γγ ∼ 10 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: 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) The 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 factor ∼ 2 in the ICS radiation, where the energies in mc 2 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 SEDs and high-energy light curves are the magnetic inclination angle, α, the observer angle ζ, the low and high values of accelerating electric field, R low acc = eE low /mc 2 and R high acc = eE 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.

RESULTS
Here we discuss the results of our modeling, first presenting the calculated phase-averaged spectral energy distributions (SEDs) for all sources, followed by the energy-dependent sky maps and light curves.For all sources, the E low , E high and J/J GJ parameters are adjusted so that the SC component matches the hard X-ray to GeV data, particularly the high-energy cutoff in the Fermi SEDs which is especially sensitive to E high .The hard X-ray SED is most sensitive to the value of E low .The M + and r radio are adjusted so that the pair SR matches the observed IR to soft X-ray data.The viewing angle, ζ, is chosen to best match the observed Fermi and VHE light curves.After these parameters are adjusted, the ICS component fluxes and SEDs are determined without further parameter adjustments (the ICS intensity thus being anchored by the lower-energy components).

Vela pulsar
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 who used the same pair spectrum for a pure dipole and the same E low , E 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 primarily depends 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 which 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.
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 non-dipolar field but the same M + .Since the non-dipolar 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 low and E high values are the same.Viewing angles in the range 50 • to 70 • make minor differences to the primary ICS component.
Figure 4 shows a model SED for the same parameters and pair spectrum, but with a lower pair multiplicity, M + = 1 × 10 3 .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 the ζ ∼ 64 • measured from the pulsar wind nebula torus image fits (Ng & Romani 2008).
In Figure 5 we show the emission sky maps (intensity vs. ζ vs. φ) in three different energy ranges, as indicated in the figure.Also shown are light curves for the same energy ranges, that 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 which they emit parallel to their motion then arrive to the observer at the same phase, producing a sharp pulse each time our line of sight crosses the current sheet.Current sheet crossings can occur once, twice or not at all, depending on the viewing angle.In these models (just like in Two-Pole-Caustic/SG 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.
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.In the lowest energy range, the emission is more uniform along the caustic 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. (submitted).Emission primarily in the second peak at energies above 10 GeV is observed for Vela, as well as Geminga and B1706−44 (see §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. (submitted).This is confirmed in Figure 6, which shows a sky map of the maximum radius of curvature, ρ max c , of the trajectories of particles that emit photons in each direction ζ and φ.As in the intensity sky maps, 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 ρ max c closely reflects the caustic emission sky pattern.This indicates that the largest curvature radii occur in and near the current sheet at radii outside the light cylinder where the main caustics form.Also observe that the ρ max c 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 ρ max c 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 to P2 ratio and the narrowing of the peak width with increasing energy as a characteristic of CR in a FF field geometry, where the spectral cutoff of CR is predicted to be The above expression holds when the emitting particles are in 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, given the larger value of ρ c in this case.
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.
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 the ρ 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 §3.2).

Crab pulsar
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 HAWC.However, HAWC is primarily designed to operate in imaging mode and has presently not implemented searches for pulsed emission (B.Dingus, priv.comm.).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 the ζ = 72 • pair SR spectrum, but the ζ = 72 • pair SSC spectrum better matches the MAGIC points.
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 compared to that of Vela.

Geminga pulsar
In Figure 10 we show model SEDs for the Geminga pulsar for α = 75 • , ζ = 55 • and two different variations of R acc .The solid curve assumes two values inside and outside the light cylinder, R low acc = 0.04 and R high acc = 0.15, while the dashed spectra assumes a single value R low acc = R high acc = 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 low value is presently 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, both from 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 multi-valued 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 the R low acc = 0.04.As is apparent from Figure 1, the computed Geminga pair spectrum 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 compared with Vela.The light cylinder distance is also three times larger than for Vela, so that the 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, 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 sky map intensity distribution 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.

PSR B1706−44
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 current sheet caustics 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 to 100) than for Geminga, primarily because its pair spectrum is higher in multiplicity and extends to higher energy.Although the B1706−44 pair spectrum 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 leads 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 particle trajectory radii of curvature 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.
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 for this model ζ, the P1/P2 ratio decreases with increasing energy.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. 2019).In the highest energy band, > 70 GeV, we predict that only P2 will be present.

PSR J0218+4232
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 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.
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, as 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 hrs 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.

DISCUSSION
We have presented models of broadband emission, covering IR to VHE energies, from a number of rotation-powered pulsars assuming the global magnetic field structure of an FF magnetosphere.While much of the method of calculation is very similar to that of HK15, the expanded photon energy range that is modeled here (from 12.5 to 18 energy decades) fundamentally changes the results to reveal a number of new features.Prime among these is the appearance of a VHE emission component at photon energy of 0.1 -30 TeV from accelerated primary particles in the current sheet scattering very low energy photons produced by pair SR.There was a hint of this component in Figures 5 and 6 of HK15, but it was cut off above the maximum modeled energy of 5 TeV and suppressed in flux by a cutoff in the pair SR spectrum below the minimum modeled energy 1 eV, preventing the primary particles from scattering IR photons in the Thompson limit.Such a limitation was remedied for the Vela model in H18 and is now also for the Crab and several over sources in this paper.Another major improvement made possible by the expanded photon energy range is the ability to use the observed optical/UV spectrum to constrain the emission altitude and the multiplicity of the electron-positron pair SR component.This constrains the IR to UV photon density at the source that then determines the primary ICS flux above 10 TeV, as well as the pair SSC component.Since the pairs are injected from the surface, there is no free parameter for the IR/UV emission; only the radio emission height is free.
The results of our modeling predict that three main components contribute to pulsar VHE emission.First, we have shown that emission up to 100 GeV can be produced by primary SC and the tail of this spectrum can account for the pulsed emission observed by H.E.S.S.-II from Vela and B1706−44, by MAGIC from Geminga and by Fermi for J0218+4232.The extended sub-exponential, power-law-like high-energy tail of this component results from emission of particles in the CR regime of SC with different radiation-reaction limited energies.The particle trajectory radii of curvature vary at different altitudes along the current sheet, producing a range of CR cutoffs (see Eqn [16]) that all combine in the caustic peaks.The emission above 10 GeV is observed most strongly in the second light curve peak, extending the P1/P2 ratio decrease with increasing energy observed in many pulsars at lower energies by Fermi.This is caused by the different maximum radii of curvature of the particle trajectories that produce the two peaks (Figure 6).ICS may also account for this emission but its spectrum would need to join smoothly to the Fermi spectrum and the P1/P2 trend might not be smooth or may occur at higher energy.
The second predicted VHE emission component comes from pair SSC.The SSC from pairs produces a broad SED whose peak typically occurs between 1 − 10 GeV for young pulsars and around 100 GeV for MSPs.In all of the source models we have presented in this paper, the pair SSC component flux lies well below, and is thus obscured by, the primary SC emission below 100 GeV.Above 100 GeV, the pair SSC emission extends to around 1 TeV and may be detectable in some pulsars like the Crab and Crab-like pulsars.A pair SSC component high enough to be detectable requires 1) high pair multiplicity (the SSC flux is proportional to the square of M + ) 2) high magnetic field strength at the light cylinder, B LC (to boost the SR level), and 3) low pair energies so that the scattering is not in the extreme Klein-Nishina regime.
The third predicted VHE emission component comes from primary ICS, radiation from the most highly accelerated particles in the current sheet scattering low-energy pair SR.Our models predict primary ICS components for all pulsars but for many the flux is well below current detection thresholds.This component depends on all five free parameters of our model and also requires treating a very large range of photon energies, at least 18 decades.For some pulsars with large IR fluxes, it may be necessary to include one or two decades below what we have considered here.The most important constraint that can be derived from a measurement of this component is the maximum energy of the accelerated particles, which by conservation of energy is simply equal to the high-energy cutoff of the SED.Even without a precise measurement of the SED high-energy cutoff, the maximum photon energy detected puts a lower limit on the maximum particle energy.Particles with energies above 10 TeV are radiating in the CR regime of SC, thus ruling out models suggesting that the Fermi GeV emission is SR (Cerutti et al. 2016b;Philippov & Spitkovsky 2018).Most of the primary ICS emission should appear in the second light curve peak for most values of α and ζ.Since it has been shown that P2 is radiated by the highest-energy particles (Barnard et al., submitted), those same particles will produce the highest ICS energies, as shown in Figure 7. Pulsar characteristics needed for detectable primary ICS emission are 1) particle energies at least 10 TeV 2) high IR/optical flux (this is not necessarily correlated with pair multiplicity but with high radio luminosity) and 3) a relatively small distance between the low-altitude pair emission (radio altitude) and the current sheet, favoring shorter pulsar periods.
In our model, the primary ICS spectrum is predicted to sharply cut off at the maximum particle energy.The cutoff energy is constrained by fitting the Fermi high-energy cutoff, which is the cutoff in the CR spectrum of the same particles.Since the FF magnetic field structure determines the curvature radii of the particle trajectories, the Fermi spectral cutoff sets the maximum particle energy assuming the particles have reached radiation-reaction limit.Therefore, with the assumptions of a FF field, particles in CR-reaction limit and emission from the current sheet, the primary ICS cutoff is constrained to be 20-30 TeV.
In addition to a constraint on maximum particle energy, a detection of the primary ICS component combined with the cutoff in the SC/CR spectrum would provide separate constraints on maximum radius of curvature and E high .The primary ICS and SC cutoffs constrain different combinations of ρ max c and E high , the primary ICS cutoff through γ CRR by Eqn (17) and the primary SC cutoff through E c through Eqn (16).Measurement of the primary ICS cutoff for a variety of pulsars with different Fermi light curves could test whether the P1/P2 ratio decrease with energy is caused by spatial variation of ρ max c or E high , or a combination of both, in addition to probing the effect of α and ζ on this trend.Finally, the shape of the primary ICS spectrum is dependent on primary particle energies as a function of altitude, since it is reflected in the shape of the SC spectrum.Measurement of the ICS spectrum below the cutoff can therefore constrain the primary particle energy distribution.
Given the model requirements for production of the different VHE emission components, we can make predictions for the pulsars that offer the best prospects for detection by VHE telescopes.For the pulsars in the Second Fermi Pulsar Catalog (Abdo et al. 2013), there is a very weak correlation between the GeV cutoff E c (with a range less than a decade) and magnetic field at the light cylinder (with a range of about 3 decades).The sub-100 GeV primary SC component is thus not likely to be very sensitive to parameters of individual pulsars, so the VHE detectability will scale with the Fermi flux (as this indicates a higher flux of primary particles).Consequently, several of the pulsars with the largest Fermi GeV fluxes, Crab, Vela, B1706−44 and Geminga, have already been detected by VHE telescopes below 100 GeV.On the other hand J0218+4232, with γ-ray flux a factor of 20-50 times lower, has not been detected below 100 GeV by MAGIC.In this regard, Fermi pulsars with high fluxes that are good targets for VHE telescopes are PSR J0007+7303 (CTA1), PSR J2021+3651 and J2021+4026.Detection of the pair SSC component will be much more sensitive to intrinsic pulsar parameters, since high pair multiplicity and B LC are required.Consequently, Crab-like pulsars such as B0540−69 and B1509−58 are the best prospects for detection of the pair SSC.We argued that the pair SSC component has already been detected from the Crab by MAGIC.As we discussed in §3.5, MSPs will not be good targets for detection of this component since their pair energies are too high.Finally, the primary ICS component above 100 GeV is stronger for pulsars with pair spectra extending to low energies, producing high levels of pair SR, and smaller light cylinder radii.Based on these requirements, both Crab-like pulsars and Vela-like pulsars with smaller periods will offer the best prospects for detection of this component.We predict that the primary ICS emission of the Crab should be presently detectable by HAWC and MAGIC.We have noted that there is a predicted range of spectral fluxes, depending on selected model parameters.Thus, continued observations and modelling will allow us to constrain parameter degeneracies.
Our modeling of pulsar broadband spectra could certainly be improved in several respects.It would be more accurate to use the E distribution in the magnetosphere of dissipative models, such as those discussed in Kalapotharakos et al. (2018).These models show that E is not only a function of altitude but also of azimuth around the current sheet and of distance above the current sheet.The pair spectra could be computed more accurately by performing (at least) 2D simulations of the time-dependent acceleration gaps to determine both the time-averaged energy distribution and multiplicity of particles escaping the gap that seed the full cascade.This would eliminate the pair multiplicity as a free parameter and produce more accurate pair spectra.The photon-photon pair attenuation could be treated more accurately by using the stored SR emissivities, used now for the SSC and ICS emission, for the soft-photon densities.Ultimately, the field and particle distributions from PIC simulations could be used to model the whole broadband spectrum.At present, this option is not feasible for several reasons.First, scaling the low part of the particle energy spectrum of the PIC particle energies (that are already artificially low) up to those of real pulsars is problematic.Second, the current dynamic range of particle energies that can be treated in PIC models is very small (only a few decades), so that most low-energy particles are in a thermal distribution.Third, the PIC simulations do not have the resolution to treat the PC pair cascades self-consistently so the pair spectrum from local simulations would still need to be injected separately.Nevertheless, we plan in the future to add the improvements that are most feasible.Figure 1.Electron-positron pair spectra for the sources modeled in this paper, from a polar cap cascade simulation.The solid curves for Vela, B1706-44 and Geminga, and the dashed curve for the Crab assume centered dipole magnetic fields.The dashed curve for Vela assumes a non-dipolar surface field with offset parameter = 0.4 (Harding & Muslimov 2011), while the J0218+4232 spectrum assumes a non-dipolar field with = 0.6.The solid curve for the Crab is the dipole pair spectrum with a power-law extension.10 9 10 7 10 5 10 3 10 1 10 1 10 3 10 5 10 7 Energy (MeV)    (dashed lines) for the pair spectra with power-law extension shown in Figure 1, assuming M+ = 3×10 5 .Data points are from Kuiper et al. (2001), Abdo et al. (2013) (http://fermi.gsfc.nasa.gov/ssc/data/access/lat/2ndPSR catalog/), Ansoldi et al. (2016) and Sollerman et al. (2019).The HAWC sensitivity curve for 50 hr of observation is also shown (Abeysekara et al. 2017) .

Figure 6 .
Figure 6.Top: Sky map of particle trajectory maximum radius of curvature, ρ max c , in units of light cylinder radius, RLC.Bottom: Slice through the sky map (solid white line) showing the variation of ρ max c with pulse phase at ζ = 70 • (normalized to 1.0 at maximum).

Figure 7 .Figure 8 .
Figure 7. Model SED for the Vela model of Figure 3 for the phase ranges of the first and second peaks, P1 and P2.One can clearly see the relatively larger SC and ICS spectral cutoffs for P2.

Figure 9 .Figure 10 .
Figure 9. Emission sky maps and light curves for the Crab extended power-law pair model shown in Figure 8, for four different energy ranges as labeled.The light curves assume a viewing angle ζ = 72 • , as indicated by the horizontal white lines in the sky maps.

Figure 11 .Figure 12 .
Figure 11.Emission sky maps and light curves for the Geminga model shown in Figure 10, for three different energy ranges as labeled.The light curves assume a viewing angle ζ = 55 • , as indicated by the horizontal white lines in the sky maps.

Figure 13 .Figure 14 .
Figure 13.Emission sky maps and light curves for the B1706-44 model shown in Figure 12, for three different energy ranges as labeled.The light curves assume inclination angle α = 45 • and viewing angle ζ = 53 • , as indicated by the horizontal white lines in the sky maps.

Table 1 .
Source and Model Parameters