HAWC Observations Strongly Favor Pulsar Interpretations of the Cosmic-Ray Positron Excess
Dan Hooper Ilias Cholis Tim Linden Ke Fang


FERMILAB-PUB-17-061-A




0.2in























ORCID: http://orcid.org/0000-0001-8837-4127
dhooper@fnal.gov

icholis1@jhu.eduORCID: http://orcid.org/0000-0002-3805-6478
ORCID: http://orcid.org/0000-0001-9888-0971
linden.70@osu.edu
ORCID: http://orcid.org/0000-0002-5387-8138
kefang@umd.edu

[a]Fermi National Accelerator Laboratory, Center for Particle Astrophysics, Batavia, IL 60510
[b]University of Chicago, Department of Astronomy and Astrophysics, Chicago, IL 60637
[c]University of Chicago, Kavli Institute for Cosmological Physics, Chicago, IL 60637
[d]Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, Maryland, 21218
[e]Ohio State University, Center for Cosmology and AstroParticle Physics (CCAPP), Columbus, OH  43210
[f]University of Maryland, Department of Astronomy, College Park, MD, 20742
[g]Joint Space-Science Institute, College Park, MD, 20742


Recent measurements of the Geminga and B0656+14 pulsars by the gamma-ray telescope HAWC (along with earlier measurements by Milagro) indicate that these objects generate significant fluxes of very high-energy electrons. In this paper, we use the very high-energy gamma-ray intensity and spectrum of these pulsars to calculate and constrain their expected contributions to the local cosmic-ray positron spectrum.  Among models that are capable of reproducing the observed characteristics of the gamma-ray emission, we find that pulsars invariably produce a flux of high-energy positrons that is similar in spectrum and magnitude to the positron fraction measured by PAMELA and AMS-02. In light of this result, we conclude that it is very likely that pulsars provide the dominant contribution to the long perplexing cosmic-ray positron excess.

HAWC Observations Strongly Favor Pulsar Interpretations of the Cosmic-Ray Positron Excess
    and Ke Fang^f,g
    
=========================================================================================






§ INTRODUCTION



Measurements of the cosmic-ray positron fraction by the PAMELA <cit.> and AMS-02 <cit.> experiments (as well as HEAT <cit.>, AMS-01 <cit.> and Fermi <cit.>) have identified an excess of high-energy positrons relative to the standard predictions for secondary production in the interstellar medium (ISM). This result indicates that significant quantities of ∼0.01-1 TeV positrons are produced within the local volume (the surrounding ∼ kpc) of the Milky Way. The origin of these particles has been debated, and possible sources include nearby and young pulsars <cit.>, as well as annihilating dark matter particles <cit.>. This excess could also plausibly arise in nearby supernova remnants, which may be able to produce and accelerate secondary positrons before they escape into the ISM <cit.>. 

From the catalog of known pulsars, Geminga (PSR J0633+1746) and B0656+14 (PSR J0659+1414, thought to be associated with the Monogem supernova remnant) are particularly attractive candidates to account for the observed positron excess. These pulsars are each relatively young (370,000 and 110,000 years, respectively) and are located within a few hundred parsecs of the solar system (250^+230_-80 and 280^+30_-30 pc, respectively <cit.>). Much of what we know about these and other pulsars is based on gamma-ray observations. Geminga, in particular, is an extremely bright source of GeV-scale emission, thought to originate from the pulsar itself <cit.>. In contrast, observations of Geminga at very high energies reveal a significantly extended morphology. In particular, the Milagro Collaboration has reported the detection of gamma-ray emission at ∼ 35 TeV with a flux of (3.77 ± 1.07) × 10^-16 TeV^-1 cm^-2 s^-1 and an extended profile with a full-width at half-maximum of 2.6^+0.7_-0.9 degrees <cit.>.



Very recently, the HAWC Collaboration has presented their measurements of the Geminga pulsar, confirming its angular extension (∼ 2^∘ radius) and reporting a flux of (4.87 ± 0.68) × 10^-14 TeV^-1 cm^-2 s^-1 at 7 TeV, with a local spectral index of -2.23 ± 0.08 <cit.> (see also Refs. <cit.>). HAWC has also reported the detection of very high-energy emission from B0656+14, with a similar degree of spatial extension and with a flux and spectral index at 7 TeV of (2.30 ± 0.73) × 10^-14 TeV^-1 cm^-2 s^-1 and -2.03 ± 0.14, respectively <cit.>. 








These observations by HAWC and Milagro allow us to conclude that these pulsars, in fact, deposit a significant fraction of their total spindown power into high-energy leptons. Furthermore, we will show in this paper that the flux of leptons required to explain these observations is roughly equal to that required for the Geminga and B0656+14 pulsars to produce an order one fraction of the positron excess observed by PAMELA and AMS-02. The spectrum, morphology and intensity of the gamma-ray emission measured by HAWC and Milagro leave us with little choice but to conclude that nearby pulsars are likely to be the dominant source of the observed cosmic-ray positrons. 




§ INVERSE COMPTON SCATTERING OF VERY HIGH-ENERGY ELECTRONS AND POSITRONS NEAR PULSARS


The gamma-ray emission observed from Geminga and B0656+14 by HAWC and Milagro is almost certainly generated through the inverse Compton scattering of very high-energy leptons. The angular extension of this signal rules out other scenarios, with the possible exception of pion production. A pion production origin, however, would require an unrealistically large quantity ( 10^46 erg) of 𝒪(10^2) TeV protons to be confined to the region surrounding Geminga for  10^5 years. Such a scenario could also be constrained to some degree by the lack of TeV neutrinos detected by the IceCube experiment <cit.>. 



To study the diffusion and energy losses of electrons and positrons produced in nearby pulsars, we utilize the standard propagation equation:


    ∂/∂tdn_e/dE_e(E_e,r,t)    =   ▽⃗·[D(E_e) ▽⃗dn_e/dE_e(E_e,r,t) - v⃗_c  dn_e/dE_e(E_e,r,t) ] 
       +   ∂/∂ E_e[dE_e/dt(r)  dn_e/dE_e(E_e,r,t)    ]  + δ(r) Q(E_e,t),


where dn_e/dE_e is the differential number density of electrons/positrons at a distance r from the pulsar, D is the diffusion coefficient, v⃗_c is the convection velocity, and the source term Q describes the spectrum and time profile of electrons/positrons injected into the ISM. Energy losses are dominated by a combination of inverse Compton and synchrotron losses, and are given by <cit.>:














    -dE_e/dt(r)    =   ∑_i 4/3σ_T ρ_i(r) S_i(E_e) (E_e/m_e)^2 + 4/3σ_T ρ_ mag(r) (E_e/m_e)^2  
       ≡    b(E_e,r)  (E_e/ GeV)^2,

where σ_T is the Thomson cross section and


    b(r) ≈  1.02 × 10^-16  GeV/ s ×[ ∑_i ρ_i(r)/ eV/ cm^3  S_i(E_e) + 0.224  (B(r)/3 μG)^2 ].


The sum in this expression is carried out over the various components of the radiation backgrounds, consisting of the cosmic microwave background (CMB), infrared emission (IR), starlight (star), and ultraviolet emission (UV). Throughout our analysis, we adopt the following parameters: ρ_ CMB=0.260 eV/cm^3, ρ_ IR=0.60 eV/cm^3, ρ_ star=0.60 eV/cm^3, ρ_ UV=0.10 eV/cm^3, ρ_ mag=0.224 eV/cm^3 (corresponding to B=3 μG), and T_ CMB =2.7 K, T_ IR =20 K, T_ star =5000 K and T_ UV =20,000 K. For low to moderate electron energies, these parameters correspond to a value of b ≃ 1.8 × 10^-16 GeV/s. At very high energies (E_e  m^2_e/2T), however, inverse Compton scattering is further suppressed by the following factor:


    S_i (E_e) ≈45   m^2_e/64 π^2 T^2_i/(45   m^2_e/64 π^2 T^2_i)+(E^2_e/m^2_e).


To solve Eq. <ref>, we calculate the distribution of the electrons and positrons that were emitted a time t ago, and then sum the contributions produced over different periods of time. Considering an injected spectrum of the form Q(E_e,t) =  δ(t) Q_0 E^-αexp(-E_e/E_c), the solution to Eq. <ref> (neglecting convection) is given by:












    dn_e/dE_e(E_e,r,t) =Q_0   E^2-α_0/8π^3/2E_e^2  L^3_ dif(E_e,t) exp[-E_0/E_c]  exp[-r^2/4L^2_ dif(E_e,t)],


where E_0 ≡ E_e/(1-E_e b t) is the initial energy of an electron that has an energy of E_e after a time t, and the diffusion length scale is given by:


    L_ dif(E_e,t)    ≡   [   ∫^E_e_E_0D(E')/-dE_e/dt(E') dE'   ]^1/2,
    
       =   [D_0/b (E_e/ GeV)^1-δ (1-δ)(1-(1-E_e b t)^1-δ)]^1/2.


In the last step we had adopted a parameterization of D(E_e) = D_0 E_e^δ for the diffusion coefficient. Note that for E_e b t > 1, there are no electrons/positrons of energy E_e and the contribution to dn_e/dE_e is set to zero. 

To account for the time profile of the electrons and positrons injected from a given pulsar, we adopt a function proportional to the spin-down power (the rate at which the pulsar loses rotational kinetic energy through magnetic dipole braking) <cit.>:


    Ė   =    -8π^4 B^2 R^6/3 c^3 P(t)^4
       ≈    1.0 × 10^35  erg/ s×(B/1.6 × 10^12  G)^2  (R/15   km)^6  (0.23   s/P(t))^4,


where B is the strength of the magnetic field at the surface of the neutron star, R is the radius of the neutron star, and the rotational period evolves as follows:


    P(t) = P__0 (1+t/τ)^1/2,


where P__0 is the initial period, and τ is the spindown timescale:


    τ   =   3c^3 I P__0^2/4π^2B^2 R^6
       ≈    9.1 × 10^3   years (1.6× 10^12  G/B)^2 (M/1.4   M_⊙)  (15   km/R)^4  (P__0/0.040   sec)^2.


In Eqs. <ref> and <ref>, we had adopted benchmark values for the neutron star's magnetic field, radius and mass, chosen to match Geminga's observed period and its rate of change.


At energies within the range measured by MILAGRO and HAWC, inverse Compton scattering yields photons with energies not very far below that of the incident electrons and positrons, E_γ∼ E_e. Adopting this approximation, the angular profile of gamma rays generated through Inverse Compton scattering is given by:


    Φ_γ(E_γ=E_e,ψ)    ∝   ∫Ė  dt  ∫_ losdn_e/dE_e(E_e,r,t)   ρ_ rad(r)   dl 
       ∝   ∫τ^2/(t+τ)^2  dt ∫_ los  e^-r^2/4L^2_ dif(E_e,t) ρ_ rad(r)   dl,


where ψ is the angle observed away from the pulsar, and r^2=l^2+d^2-2ld cosψ, where d is the distance between the pulsar and the observer. If we adopt a uniform distribution of radiation in the vicinity of the pulsar, this reduces to a profile of the form Φ_γ(ψ) ∝exp[-d^2 sin^2 ψ/4L_ dif^2(E_e,t)]. Observations of Geminga by both Miligro and HAWC indicate that the very high energy gamma-ray emission from this source is extended over a region of a few degrees across the sky. This in turn requires a diffusion length given by L_ dif(E_e) ≃ (250   pc)  sin (0.5× 2.6^∘)/2(ln(2))^1/2≃ 2.6   pc. In contrast, adopting parameters appropriate for the ISM (D_0≃ 2× 10^28 cm^2/s, δ≃ 0.4, b=1.8× 10^-16 GeV/s), we find


    L_ dif(E_e,t) ≃ 200   pc  (35  TeV/E_e)^0.3 (1-(1-E_e b t)^0.6)^1/2.


Assuming conditions for cosmic ray transport that are similar to those found in the ISM, this calculation shows that we should have expected the Inverse Compton emission observed at very high-energies to be extended over a scale of ∼ 60^∘, dramatically more than the ∼ 2^∘ extension reported by both Milagro and HAWC (see Fig. <ref>). 






To resolve this puzzle, one might be tempted to consider the possibility that the pulsar is surrounded by a dense radiation field, which intensifies the resulting inverse Compton emission from the surrounding parsecs. The problem with this scenario, however, is that there is not nearly enough power available to generate the required density of radiation. More quantitatively, in order for a r^-2 profile of radiation to exceed the energy density of the CMB at a distance of 1 parsec from the pulsar would require an amount of energy equivalent to more than ten times the total spin-down power of Geminga. Adopting an extreme benchmark in which 100% of Geminga's energy budget is transferred into radiation, the profile of Inverse Compton emission is altered only very modestly; by less than 10% at ψ = 1^∘ (see Fig. <ref>). Based on these considerations, it does not appear that local concentrations of radiation play a significant role in explaining the angular extent of the gamma-ray emission observed from Geminga or B0656+14.




A more likely solution is that the conditions that dictate cosmic-ray diffusion around Geminga and B0656+14 are very different from those found elsewhere in the ISM, leading energetic leptons to escape from the surrounding regions much slower. In order to accommodate the observed extension in this way, we require that these particles only diffuse a distance of a few parsecs before losing most of their energy. For an energy loss time of 5000 years (corresponding to E_e=35 TeV and b=1.8× 10^-16 GeV/s), this requires a diffusion coefficient of D ∼10^27 cm^2/s (see Fig. <ref>). Although this is significantly smaller than that found in the bulk of the ISM, we note that it is similar to that predicted for standard Bohmian diffusion, D_ Bohm = r_L c/3 ≈ 1.2 × 10^27 cm^2/s ×  (E_e/35   TeV)(μ G/B).








If pulsars such as Geminga are typically surrounded by a region with inefficient diffusion (D ≪ D_ ISM), the volume of such regions must be fairly small to avoid conflicting with secondary-to-primary ratios in the cosmic-ray spectrum as measured at Earth. In particular, if such regions have a typical radius of r_ region, then such regions will occupy roughly the following fraction of the volume of the Milky Way's disk:


    f    ∼   N_ region×4π/3 r^3_ region/π R^2_ MW× 2 z_ MW
       ∼    0.25 ×(r_ region/100   pc)^3  (Ṅ_ SN/0.03   yr^-1)  (τ_ region/10^6   yr)  (20   kpc/R_ MW)^2  (200   pc/z_ MW),


where Ṅ_ SN is the rate at which new pulsars appear in the Galaxy, τ_ region is the length of time that such regions persist, and N_ region = Ṅ_ SN×τ_ region is the number of such regions present at a given time. The quantities R_ MW and z_ MW denote the radius and half-width of the Galaxy's cylindrical disk. Combined with Milagro and HAWC observations of Geminga and B0656+14, these considerations suggest 5  pc r_ region 50  pc, for which there will be little impact on the observed secondary-to-primary ratios (other than the positron fraction).





In Fig. <ref>, we plot the gamma-ray spectrum from the region surrounding Geminga, compared to the measurements by HAWC and Milagro (at 7 and 35 TeV, respectively). In performing this calculation, we utilized the full differential cross section for inverse Compton scattering <cit.>. In the left (right) frame, we have adopted a spectrum of injected electrons with an index of α=1.5 (1.9), and in each case selected a value of E_c that provides the best-fit to the combination of these two measurements. We have also allowed for the possibility that convective winds play a role in cosmic-ray transport <cit.>, moving particles away from the pulsar at a velocity given by either v_c = 55.4 km/s × (r_ region/10   pc) (blue dashed) or v_c = 554 km/s × (r_ region/10   pc) (green solid). In these four cases, the best-fits were found for E_c=44 TeV (α=1.5, low convection), E_c=35 TeV (α=1.5, high convection), E_c=67 TeV (α=1.9, low convection), and E_c=49 TeV (α=1.9, high convection). In each case, convection dominates over diffusion in transporting cosmic rays out of the region surrounding the pulsar.

In addition to their flux measurement, the HAWC Collaboration has also reported a value of -2.23 ± 0.08 for Geminga's spectral slope at 7 TeV. Among the models shown in Fig <ref>, those with a low convection velocity (v_c = 55.4 km/s × (r_ region/10   pc)) predict spectral slopes at 7 TeV of -2.47 (α=1.5) or -2.59 (α=1.9). Such values are highly inconsistent with that reported by HAWC. In contrast, for those models with a higher degree of convection (v_c = 554 km/s × (r_ region/10   pc)), we instead predict a spectral slope of -2.23 (α=1.5) or -2.32 (α=1.9), in excellent agreement with HAWC's measurement. This favors scenarios in which convection plays a very important role in the transport of high-energy leptons, especially among those with energies below ∼10 TeV. We lastly note that among this range of models, between 7.2% and 29% of Geminga's total current spindown power is being deposited into electron-positron pairs with E_e > 10 GeV. 





















§ IMPLICATIONS FOR THE COSMIC-RAY POSITRON EXCESS


Although the angular extent of the emission observed from Geminga and B0656+14 by HAWC and Milagro indicates that very high-energy (E_e  10 TeV) leptons are effectively trapped within a few parsecs of these sources, the same fate need not be experienced by lower energy electrons and positrons. In particular, even a modest degree of convection (i.e. the streaming of particles away from the source at a constant velocity) can remove sub-TeV leptons from the region before they lose a substantial fraction of their energy, while higher energy leptons lose the vast majority of their energy to inverse Compton scattering before escaping the same region <cit.>.






To address this more quantitatively, we can compare the timescale for energy losses, τ_ loss = 1/E_e b, to that for the escape from the region, τ_ escape = r_ region/v_c. For convective winds with a velocity of v_c, particles with E_e ≪ v_c /b r_ region are left largely unaffected, while those with E_e ≫ v_c /b r_ region lose the majority of their energy before escaping. In the left frame of Fig. <ref>, we plot the fraction of energy that an electron loses before escaping a region of radius 10 parsecs for several values of v_c. It is expected that future observations by the Cherenkov Telescope Array (CTA) will provide an important test of this transition by measuring the intensity and angular extent of the ∼ 0.1-10 TeV emission from Geminga and B0656+14.


Although we have focused here on a scenario in which convection is responsible for expelling sub-TeV leptons from the regions surrounding Geminga and B0656+14, other means may also be possible. For example, if we simply introduce a diffusion coefficient with no energy dependance (δ=0), we effectively mimic the effects of convection (see the right frame of Fig. <ref>). 





In Fig. <ref>, we plot the cosmic-ray positron fraction as measured by AMS-02, compared to the predictions from the Geminga pulsar, using the same choices of parameters as adopted in Fig. <ref>. In each frame, the solid black curve denotes the contribution from standard secondary production in the ISM, while the other curves include both this contribution and that from Geminga. We remind the reader that those models with only weak convection (dotted blue curves) do not lead to a spectral index compatible with the measurement of HAWC, and thus should be viewed as a poor fit to the data.

The positron fraction presented in Fig. <ref> includes a distinctive feature at 400-500 GeV, which is a result of energy losses. More specifically, a positron with an infinite energy will be reduced over a time t to an energy of E_e=(bt)^-1, which for t=370,000 years (the age of Geminga) yields a final energy of 476 GeV. Any positrons from Geminga above this energy were injected at later times and thus have not cooled to the same extent. 

The main lesson from the results shown in Fig. <ref> is that when the spectral shape and overall normalization of Geminga are fixed to reproduce the results of HAWC and Milagro, this pulsar is found to generate a non-negligible portion of the observed positron fraction. That being said, the overall size of this contribution to the cosmic-ray positron flux can vary by a factor of order unity depending on the precise values of the convection velocity, v_c, and spindown timescale, τ (see Eq. <ref>) that are adopted. The impact of the convection velocity is clearly evident in Fig. <ref>. The resulting positron flux scales approximately as τ^-1 (we have adopted a value of τ=9.1× 10^3 years). Furthermore, the time profile of a pulsar's emission could plausibly depart to some extent from that predicted from standard magnetic dipole braking <cit.>, potentially altering the normalization of the positron flux predicted here, as well as the inferred age of the pulsar.




Finally, in Fig. <ref>, we plot the contributions to the positron fraction from the Geminga and B0656+14 pulsars, as well as the average contribution from those pulsars located more than 500 parsecs away from the Solar System. For each source, we have adopted α=1.9, E_c=49 TeV, v_c=554 km/s × (r_ region/10   pc), and normalized their contributions with τ=4.3× 10^3 years, and adopting a total birth rate of two pulsars per century in the Milky Way.[We produce nearly identical results if we instead adopt our default value for τ≃ 9.1× 10^3 and a somewhat higher value for the convection velocity, v_c ≃1160 km/s × (r_ region/10   pc).] For other details regarding the calculation of the contribution from distant pulsars, we direct the reader to Ref. <cit.>. In reality, we expect many of these parameters to vary from pulsar-to-pulsar, making a detailed prediction of this kind difficult and possibly unreliable. That being said, this exercise clearly demonstrates that in light of the measurements by HAWC and Milagro, it appears very likely that a sizable fraction of the observed positron excess originates from this class of sources. In addition, we note that it was recently shown that the stochastic acceleration of cosmic-ray secondaries in supernova remnants is also likely to contribute to the local positron flux <cit.>. 













 
 
 

   
    
  






















§ SUMMARY AND CONCLUSIONS


In this paper, we have made use of measurements by the very high-energy gamma-ray telescopes HAWC and Milagro to better understand and constrain the injection of high energy electrons and positrons from the nearby pulsars Geminga and B0656+14. The angular extension of the  TeV gamma-ray emission observed from these pulsars indicates that very high-energy leptons are effectively trapped within the surrounding several parsecs around these sources. Furthermore, their very high-energy gamma-ray spectra indicate that lower energy leptons are able to escape more easily, suggesting the presence of strong convective winds, with velocities of several hundred kilometers per second. 

In models that are able to reproduce the characteristics of the gamma-ray emission reported by HAWC and Milagro, these pulsars invariably provide a significant contribution to the local cosmic-ray positron spectrum, and thus to the measured positron fraction. Although it is not yet possible to precisely predict the normalization of the positron flux from these sources, these results show that Geminga and B0656+14 are expected to generate a significant fraction of the observed high-energy positron flux. If we make the entirely reasonable assumption that other pulsars in the Milky Way behave in a similar fashion to Geminga and B0656+14, we find that it is very likely that pulsars are responsible for much, if not the entirety, of the reported positron excess. 

An important test of this conclusion will come from future Imaging Atmospheric Cherenkov Telescopes (IACTs), such as the Cherenkov Telescope Array (CTA). Although existing IACTs have not yet reported any significant detection of TeV-scale emission from Geminga or B0656+14 <cit.>, next generation telescopes will be far better suited to detect emission that is extended over the angular scales reported by Milagro and HAWC. Such a measurement is expected to be able to confirm the transition to convection-dominated transport at energies below several TeV, and enable us to produce a more detailed determination of the spectrum of electrons and positrons that are injected from this class of sources.
 


Acknowledgments. DH is supported by the US Department of Energy under contract DE-FG02-13ER41958. Fermilab is operated by Fermi Research Alliance, LLC, under contract DE- AC02-07CH11359 with the US Department of Energy. IC acknowledges support from NASA Grant NNX15AB18G and from the Simons Foundation. TL acknowledges support from NSF Grant PHY-1404311.



JHEP

