IceCube and HAWC constraints on very-high-energy emission from the Fermi bubbles
The nature of the -ray emission from the Fermi bubbles is unknown. Both hadronic and leptonic models have been formulated to explain the peculiar -ray signal observed by the Fermi-LAT between 0.1-500 GeV. If this emission continues above 30 TeV, hadronic models of the Fermi bubbles would provide a significant contribution to the high-energy neutrino flux detected by the IceCube observatory. Even in models where leptonic -rays produce the Fermi bubbles flux at GeV energies, a hadronic component may be observable at very high energies. The combination of IceCube and HAWC measurements have the ability to distinguish these scenarios through a comparison of the neutrino and -ray fluxes at a similar energy scale. We examine the most recent four-year dataset produced by the IceCube collaboration and find no evidence for neutrino emission originating from the Fermi bubbles. In particular, we find that previously suggested excesses are consistent with the diffuse astrophysical background with a p-value of 0.22 (0.05 in an extreme scenario that all the IceCube events that overlap with the bubbles come from them). Moreover, we show that existing and upcoming HAWC observations provide independent constraints on any neutrino emission from the Fermi bubbles, due to the close correlation between the -ray and neutrino fluxes in hadronic interactions. The combination of these results disfavors a significant contribution from the Fermi bubbles to the IceCube neutrino flux.
Observations by the Fermi Large Area Telescope (Fermi-LAT) have discovered two extended “bubbles” filling the regions above and below the Galactic center Su et al. (2010). Assuming that these bubbles originate in the Galactic center, their angular size translates to a physical extent exceeding 10 kpc. Intriguingly, the -ray spectrum of the Fermi bubbles does not show significant variation up to 50 away from the Galactic plane. This suggests that the particles producing the Fermi bubbles emission do not cool significantly during cosmic-ray transport. Additionally, the morphology of the bubbles shows sharp edges suggestive of a transient origin.
The Fermi bubbles are connected to previous Wilkinson Microwave Anisotropy Probe (WMAP) observations of a bright excess in the inner galaxy, known as the WMAP “haze” Finkbeiner (2004); Dobler and Finkbeiner (2008). Subsequent observations of the WMAP haze have verified this correlation, finding that the sharp edges of the Fermi bubbles are also seen in the haze emission Dobler (2012). The spectrum of the WMAP haze is softer than free-free emission, but is harder than synchrotron emission from astrophysical electrons produced through standard diffusive processes. Models typically produce the WMAP haze via synchrotron emission from an unknown, hard-spectrum cosmic-ray (CR) electron population (Su et al., 2010). The morphological similarity of the WMAP haze and Fermi bubbles suggests a similar leptonic origin for the bubbles. This CR lepton population could be formed through an episodic event that produced an intensive energy injection near the Galactic center, such as a past accretion event onto Sgr A*, or a nuclear starburst (Guo and Mathews, 2012; Guo et al., 2012; Yang et al., 2012; Bland-Hawthorn et al., 2013).
However, TeV leptons interact with both the magnetic field and interstellar radiation field (ISRF) of the Milky Way, losing energy to synchrotron radiation and inverse-Compton scattering on timescales of only 0.2 Myr. Unless the TeV electrons were injected very recently, leptonic models would generically include a significant softening in the Fermi bubbles -ray spectrum at high energies, which is not confirmed due to limited statistics at those energies. The process that produces a consistent -ray spectrum throughout the Fermi bubbles is unknown, and may hold the key to understanding the origin of the bubbles. Two classes of models have been proposed.
In the first class of models, the -ray signal is created by primary electrons accelerated near the Galactic center. In this scenario, the hard spectrum of the Fermi bubbles is produced by modifying cosmic-ray transport. If CRs are transported via diffusion, a diffusion coefficient of 2.510 cms at 1 TeV would be required. This value exceeds the average Milky Way diffusion coefficient by nearly two orders of magnitude. A strong convective wind in the Galactic center could be employed to more effectively move electrons out of the Galactic plane. However for 1 TeV electrons, this corresponds to an unphysical wind speed of 50,000 km/s (Su et al., 2010). One alternative possibility involves slower CR electron transport coupled with significant re-acceleration forces (e.g. due to Alfvén waves) which offset energy losses and keep the electron spectrum in steady state (Cheng et al., 2014; Mertsch and Sarkar, 2011; Sasaki et al., 2015). It has also been proposed that CRs could be injected by shocks near the bubble edge and diffuse away from it at high energies (Keshet and Gurwich, 2016).
In the second class of models, the primary CR population is energetically dominated by hadrons. The cooling of hadronic CRs is significantly slower than leptons, with a characteristic cooling time of , where is the number density of the cold gas within the bubble, and is the inelastic part of the pp cross section at PeV energies (Kelner et al., 2006). This allows hadronic CRs to easily propagate to the outer edges of the bubbles without significant energy losses. The hadronic interactions of these CRs produce an energetic CR lepton population in situ. This solves the CR transport problem at the cost of invoking an intermediate stage.
There are two important implications of models utilizing hadronic cosmic-ray transport to explain the Fermi bubbles. First, the flux of -rays and ee pairs produced in hadronic interactions are similar. This indicates that the majority of the observed -ray emission from the Fermi bubbles is, in fact, produced via the decay of neutral pions formed in the initial hadronic interaction. The contribution of inverse-Compton scattering from the lepton population is subdominant. Additionally, models find that the synchrotron emission in a hadronic scenario is typically 3-4 times lower than WMAP and Planck measurements (Ackermann et al., 2014). Hadronic models of the Fermi bubbles thus require a sub-dominant component of primary electrons, either accelerated in situ or transported by Galactic winds, in order to explain the WMAP emission (Crocker and Aharonian, 2011; Fujita et al., 2014).
Second, while energy losses constrain the maximum energy of leptons in the Fermi bubbles to be TeV-scale, protons can be accelerated to much higher energies. Although the CR acceleration site is uncertain, hadronic models often extend the Fermi bubbles spectrum to PeV energies. Supernova remnants are believed to accelerate CRs up to PeV energies, and CRs may be injected by the past starburst activities (Crocker and Aharonian, 2011). Alternatively, SgrA may be an efficient PeV accelerator, and is expected to have a maximum energy output sufficient to produce the bubbles during a transient event (Fujita et al., 2016). Another possibility is that protons may be accelerated at termination shocks around the edge of the Fermi bubbles, the energy of these accelerated protons may reach PeV energies in G magnetic fields, although the maximum energy can also be much smaller (Lacki, 2014; Fujita et al., 2014).
While the limited field-of-view of atmospheric Cherenkov Telescopes makes them incapable of cosntraining the TeV emission from the Fermi bubbles, the high sensitivity and wide field-of-view offered by the High Altitude Water Cherenkov (HAWC) telescope is capable of testing the extension of the Fermi bubbles to the TeV regime. Early observations by the HAWC collaboration have found no evidence for TeV -ray emission originating from the Fermi bubbles (Abeysekara et al., 2017).
Hybrid models are also possible. In particular, models with a significant primary electron component typically accelerate a population of very-high-energy (VHE) protons as well. In any hadronic or hybrid model a correlation is expected between the -ray and neutrino fluxes. In particular, any hadronic -rays must be accompanied by high-energy neutrinos emitted via the decays of charged pions produced in the hadronic interaction. On the other hand, leptonic models produce no associated neutrino emission. Thus, a discovery or non-detection of high-energy neutrinos offers the potential to differentiate hadronic and leptonic models of the Fermi bubbles (Lunardini and Razzaque, 2012). Meanwhile, the origin of the high-energy neutrinos observed by the IceCube Observatory remains a mystery (Halzen, 2016). Because the Fermi bubbles are bright at GeV energies, they have been suggested as potential contributors to the TeV neutrino sky (Crocker and Aharonian, 2011). An excess around the Galactic center was tentatively indicated in the two-year high-energy starting event (HESE) data (IceCube Collaboration, 2013), and a possible contribution of Fermi bubbles to the diffuse neutrino flux has been discussed (Ahlers and Murase, 2014; Lunardini et al., 2014a; Taylor et al., 2014; Aartsen et al., 2015; Lunardini et al., 2015; Ahlers et al., 2016). It has been suggested that both neutrino observations with IceCube as well as multi-TeV -ray observations with HAWC and other air-shower arrays are crucial to test hadronic models of the Fermi bubbles (Ahlers and Murase, 2014; Lunardini et al., 2015).
In this paper, we re-analyze the flux of neutrinos from the Fermi bubbles, utilizing the most recent four-year IceCube constraints as well as HAWC observations of the TeV -ray emission from the bubbles. By considering both atmospheric and astrophysical neutrino backgrounds in the Fermi bubbles region, we show that there is no statistically significant excess of high-energy neutrinos correlated with the Fermi bubbles. Additionally, we show that HAWC results independently rule out models where a significant IceCube neutrino flux is produced in a purely hadronic model of the Fermi bubbles. While hybrid models are still possible, we show that future HAWC observations could exclude at least half of the overlapping neutrino events as having a Fermi bubbles origin.
The paper is organized as follows. In Section II we calculate the neutrino flux within the Fermi bubbles region, and use models for the atmospheric and astrophysical backgrounds to constrain the neutrino excess associated with the Fermi bubbles. In Section III we consider constraints from HAWC null-observations of the Fermi bubbles in both hadronic and hybrid leptonic/hadronic models. Finally, in Section IV we discuss the prospects for constraining the neutrino emission associated with the Fermi bubbles with current and future experiments.
Ii IceCube Observations of the Fermi bubbles
In four years of data (1347 days), observations by the IceCube Observatory have detected 54 neutrinos with verticies inside the IceCube detector and contained energies exceeding 30 TeV. These events are known as the high-energy starting events (HESE) (IceCube Collaboration, 2013; Aartsen et al., 2014; The IceCube Collaboration et al., 2015). Of these events, 8 spatially overlap with the Fermi bubbles, as shown in Figure 1. For each event we show an ellipse that corresponds to the angular uncertainty in the event reconstruction, and an event color that corresponds to the reconstructed energy of the neutrino. All of these events are shower events, which have a good energy resolution (15%), but suffer from large uncertainties in the reconstruction of their arrival direction (10). This is comparable to the size of the Fermi bubbles, and implies that all events near the location of the Fermi bubbles must be studied in detail. In particular, we find that five events (2, 12, 14, 15, and 36) have a reconstructed direction that is centered within the bubbles. An additional three events (22, 24, 25) are in close proximity, but have reconstructed directions that are centered outside the bubbles. Event 14, which is located close to the Galactic center, is the highest-energy event in the region and is one of only three events in the full HESE dataset with reconstructed energies exceeding 1 PeV.
The poor angular resolution of these HESE events makes it difficult to evaluate the spatial correlation between each event and the Fermi bubbles. In what follows, we utilize the best-fit directions of each event, and consider only candidate events that are centered within the Fermi bubbles region. We note that weighting each event by the fraction of the point-spread function that lies within the Fermi bubbles provides similar results. We take the solid angle of the Fermi bubbles to be sr (Ackermann et al., 2014). The total number of events in the bubbles is .
In Figure 2 we show a distribution of the deposited energy for events centered within the Fermi bubbles region. The error bars are calculated using Feldman-Cousins confidence intervals Feldman and Cousins (1998). We compare these events against the expected number of events within the bubbles from atmospheric backgrounds (red), including atmospheric neutrinos from /K decays, charmed meson decays, and background muons (The IceCube Collaboration et al., 2015), as well as the expected number of HESE events in the ROI from observations of sideband regions in a similar declination range (; orange).
Additionally, we show the sum of the expected atmospheric background along with a best-fit astrophysical diffuse background, which is estimated by:
where is the expected background number in the Southern sky, scaled from the all-sky background using the average effective area of the Northern and Southern sky (IceCube Collaboration, 2013),
is the expected all-sky background number, obtained using the best-fit spectral slope of E obtained from all HESE neutrino events with a deposited energy between 60 TeV and 3 PeV (The IceCube Collaboration
et al., 2015). We also show a comparison to the background number obtained with a fixed spectrum (indicated by the grey steps).
Note that we have adopted the Southern sky (where the bubbles are located) as a reference, since the effective area of the IceCube detector is different for upgoing (Northern sky) and downgoing (Southern sky) events
By comparing our background models with the observed neutrino counts from the Fermi bubbles region, we find that there is no significant excess coincident with the bubbles. While the number of neutrinos observed in the Fermi bubbles region significantly exceeds the atmospheric background, the flux within the Fermi bubbles is consistent with the summed contribution from atmospheric and astrophysical backgrounds. In particular, the Fermi bubbles neutrino flux is well-fit by both background models that are derived from the sidebands regions at similar declinations, as well as the observed isotropic neutrino flux.
To quantify the significance of any excess from the Fermi bubbles region, we compare the number of events observed in each energy bin to the number of neutrinos expected based on the combination of atmospheric and astrophysical neutrino backgrounds. We define a test statistic (TS), calculated as the sum of the logarithms of the probability of finding the observed number of neutrinos in each bin E, assuming that the number of neutrinos at each deposited energy is independent and follows a Poisson distribution:
To determine the distribution of TS values expected from Poisson fluctuations in the background model, we utilize Monte Carlo techniques to produce a large population of mock observations, based on fluctuations of the background model. We find that in 22% of these observations, the resulting TS exceeds that computed for the observed data. Considering that most of the nearby events are shower events with large uncertainties, we also test an extreme scenario that all overlapped events are from the Fermi bubbles. Even with this assumption, we still find that 5% of realizations with background-only hypothesis resulted a TS higher than the observed value. We thus conclude that the neutrino flux from the Fermi bubbles region is consistent with background fluctuations.
We find that the best-fit flux from the background model accounts for events in the Fermi bubbles region. This implies that the best-fit flux of the Fermi bubbles (while not statistically significant), accounts for events. Utilizing Feldman-Cousins statistics, we can set a 90% confidence upper limit on the total contribution of the Fermi bubbles to the IceCube neutrino flux of 5.99 predicted events Feldman and Cousins (1998). Of the events expected from our background model, we find that events are produced by the atmospheric background, while events are produced by the isotropic astrophysical background.
The upper limit for the neutrino flux from the Fermi bubbles can then be estimated from the all-sky astrophysical flux through (Ahlers and Murase, 2014):
where is the Southern-sky HESE event number in the four-year data, and is the per-flavor astrophysical flux obtained for events between 60 TeV and 3 PeV assuming a fixed spectrum (The IceCube Collaboration et al., 2015).
Iii HAWC Observations of the Fermi bubbles
Recently, the HAWC collaboration analyzed the -ray emission from the Fermi bubbles region over 290 days of HAWC observation time. They found no statistically significant excess, and produced strong constraints on the -ray flux between energies of 1-100 TeV (Abeysekara et al., 2017). Utilizing these results, the HAWC collaboration placed stringent upper limits on the maximum neutrino flux in the Fermi bubbles, within the context of purely hadronic models. Interestingly, these observations ruled out a previous analysis of the neutrino flux from the Fermi bubbles, which did not take into account the possibility that the neutrino flux could be produced by atmospheric or astrophysical backgrounds (Lunardini et al., 2014b).
Here we generalize the results of the HAWC collaboration, and consider the potential for HAWC data to place model-independent constraints on the contribution of the Fermi bubbles to the IceCube neutrino flux. In particular, we consider two scenarios, which we call “pure hadronic” and “hybrid leptonic-hadronic”. In the pure hadronic scenario, we assume that all -rays produced in the Fermi bubbles at GeV energies are produced via hadronic processes. In this case, the detection of bright Fermi bubbles emission by the Fermi-LAT, coupled with the strong upper limits on Fermi bubbles emission by HAWC, combine to force the -ray spectrum (and by extension the neutrino spectrum) to be extremely soft. In the hybrid leptonic-hadronic scenario we instead assume that the bulk of the GeV -ray signal observed by the Fermi-LAT is produced by primary leptons, while -rays from hadronic processes are sub-dominant. This allows the spectrum of hadronic -rays and neutrinos to be relatively hard, allowing for a larger very-high-energy flux.
In each case, we fit the -ray spectrum and intensity to Fermi-LAT and HAWC observations, and then calculate the resulting neutrino spectrum under the assumption that neutrinos and -rays from hadronic interactions are correlated via the relationship (Ahlers and Murase, 2014):
where is the production rate of neutrinos and -rays. The number of neutrino events in the bubble region can then be calculated by:
where is the effective area of IceCube for contained neutrino searches averaged over the Southern sky (IceCube Collaboration, 2013). This serves as a reasonable approximation for the average effective area of IceCube in the Fermi bubbles region.
In Figure 3 we show the results of this analysis. We compare our models to four datasets: the 0.1-500 GeV measurements of the Fermi bubbles by the Fermi-LAT collaboration (black squares (Ackermann et al., 2014)), the 95% confidence upper limits from the non-detection of very-high-energy -rays in the northern bubble by HAWC (black solid bars), the 90% confidence upper limits from the non-detection of ultrahigh-energy -rays by the Chicago Air Shower Array - Michigan Muon Array experiment (CASA-MIA, olive upper-limits (Chantell et al., 1997)), and the neutrino flux calculated in this paper using 4-year HESE data (red upper-limit). Because the neutrino data is extremely sparse, we combine all energy bins in our neutrino analysis into one energy bin, depicted at a central energy of 230 TeV. We have included uncertainties in the neutrino flux stemming from the difference between the true neutrino energy and the deposited energy following Blum et al. (2014); Laha et al. (2013). We show the predicted -ray (blue dashed) and neutrino (orange solid) fluxes from the pure hadronic scenario (thick lines), as well as the hadronic portion of the emission in the hybrid hadronic-leptonic scenario (thin lines). In the next three subsections, we will investigate these models in detail, using the latest HAWC and IceCube constraints.
iii.1 Pure Hadronic Models of the Fermi bubbles
In pure hadronic models of the Fermi bubbles, the -ray spectrum is fit through a comparison of the bright 100 GeV emission observed by the Fermi-LAT, compared with the stringent upper limits on 10 TeV -ray emission produced by HAWC. We employ a proton spectrum modeled as a power-law with an exponential cutoff (i.e. ). We adopt the best-fit hadronic model with index as found byAckermann et al. (2014). We then saturate the HAWC upper limit at 10 TeV. We find that this requires an exponential cutoff in the bubbles spectrum at an energy TeV. We note that we could theoretically soften the injected -ray proton spectrum to allow higher values of E, which would consequently allow for a higher PeV neutrino flux. However, a softer -ray spectrum would provide a poor fit to the Fermi-LAT data.
Translating this -ray flux into a neutrino flux following Equation 5, we find that IceCube should detect no more than neutrinos from the Fermi bubbles over a four-year observation. In particular, this model cannot account for the observation of the PeV neutrino (Event 14) located within the Fermi bubbles region, which is the least likely event to be produced through fluctuations of the astrophysical background.
iii.2 Leptonic-Hadronic Models of the Fermi bubbles
In this scenario, we conversely assume that the bright GeV emission observed from the Fermi bubbles is produced predominantly by leptonic interactions ( below 100 GeV). Thus, the -ray spectrum from the hadronic portion of the Fermi bubbles emission can be made arbitrarily hard, so long as it does not exceed HAWC limits. We note that harder -ray spectral indices will inevitably produce larger contributions to the IceCube neutrino flux, because IceCube observes messengers at higher energies than either the Fermi-LAT or HAWC. In the case of arbitrarily hard hadronic CR injection spectra, a large IceCube neutrino flux could be produced while remaining consistent with Fermi-LAT or HAWC upper limits.
However, in this paper we adopt a reasonable lower limit on the CR proton energy spectrum of = 2.0, motivated by the spectrum obtained via diffuse shock acceleration (Axford et al., 1977; Krimsky, 1977; Bell, 1978; Blandford and Ostriker, 1978). We consider a pure power-law spectrum, and normalize the neutrino flux to the neutrino data point. Above 500 TeV, -ray spectrum is exponentially suppressed due to attenuation by the cosmic microwave background.We note that an additional exponential cutoff in the injection spectrum around 7 PeV or less is required to respect the CASA-MIA -ray upper limit around 10 PeV.
We find that this hybrid leptonic-hadronic scenario is consistent with current HAWC upper limits. Moreover, the injected cosmic-ray spectrum can be slightly softened compared to the theoretical upper limit, before the observed -ray signal would saturate the HAWC upper bound.
iii.3 Future Constraints from HAWC
So far, our results employ less than one year of HAWC data. If future HAWC observations do not find emission consistent with the Fermi bubbles, we expect these upper limits to fall as the square root of time. In Figure 4, we show the maximum number of neutrinos that would be consistent with HAWC upper limits assuming additional years of HAWC data, and no detection of -ray emission from the Fermi bubbles (red line). We calculate this upper limit using the maximum = 2.0 spectrum and normalize the -ray flux to the center of the last bin in the HAWC data, which is located at 69.5 TeV.
Additionally, we show the total number of IceCube events predicted in the Fermi bubbles region, given more years of IceCube data
We find that, at present, a -ray flux that saturates the HAWC upper limits in our hybrid leptonic-hadronic model produces neutrinos in the Fermi bubbles. This exceeds the 5 neutrinos that currently overlap with the Fermi bubbles region, and indicates that HAWC can not currently rule out a Fermi bubbles origin for the IceCube neutrinos in the hybrid leptonic-hadronic model. However, the number of allowable neutrino events falls approximately as the square-root of time, as additional HAWC constraints more stringently rule out a sizeable IceCube neutrino flux. Eventually, we note that this upper-limit becomes flat, as the additional flux sensitivity of HAWC is offset by the increase in the IceCube acceptance over time.
Figure 4 shows that if HAWC does not detect -ray emission from the Fermi bubbles within seven additional years of operation, approximately 50% of the neutrino events in the bubbles region can be excluded as having a bubbles origin. While we note that this constraint appears weak, it is based solely on the synergy between very-high-energy -rays and high-energy neutrinos, and thus has independent systematics from our models of the IceCube background. If the true neutrino flux from the Fermi bubbles is instead represented by the best fit event rate of events per four years of HESE data, HAWC will be unable to substantially constrain this model within 15 years of data. However, if the true neutrino event rate from the Fermi bubbles saturates the current IceCube upper limit, HAWC should begin to observe, or constrain a bubbles origin of this emission within 5 years. This indicates that HAWC observations are capable of testing currently viable models for the Fermi bubbles, even for extremely difficult hybrid leptonic-hadronic models. In particular, we note that throughout this section, we have assumed an E spectrum for both the hadronic -ray and neutrino flux in the hybrid leptonic-hadronic model. Softer spectral indices will produce more stringest limits from existing HAWC data.
We show that in four-years of IceCube data, eight HESE shower events overlap with the Fermi bubbles region. Five of those events are centered within the Fermi bubbles. The flux and spectrum of these events are consistent with expectations from a combination of atmospheric and astrophysical background neutrinos (p-value 0.22), allowing us to set strong upper limits on the maximum neutrino flux produced by the Fermi bubbles. These observations produce tension with models where the Fermi bubbles are produced by pure hadronic processes, indicating that the CR proton spectrum in these models must be exponentially suppressed above energies of 100 TeV.
While IceCube observations are likely to provide increasingly stringent constraints on hadronic models of the Fermi bubbles, our results indicate that future HAWC observations may be able to place even stronger constraints on the maximum neutrino flux from the bubbles. Already, the combination of Fermi-LAT and HAWC -ray observatinos rule out models where a significant IceCube neutrino flux is produced by the Fermi bubbles in a pure hadronic model. In the case of a hybrid leptonic-hadronic model, the constraining power of HAWC observations relies on the assumed spectrum of the hadronic emission component. However, even for CR spectral indices as hard as 2.0, null-observations with a few more years of HAWC data will strongly constrain models where the IceCube neutrino flux stems from true Fermi bubbles emission.
The total cosmic-ray energy in a hybrid leptonic-hadronic scenario can be estimated as follows. Fitting to the IceCube data point suggests a gamma-ray flux , corresponding to an integrated gamma-ray flux before attenuation, assuming an spectrum and PeV as discussed in Sec. III.2. Then the -ray luminosity is , where is the distance to the bubbles. This implies a total CR proton energy . The energy contained in the electron population needed to produce the bulk of the observed GeV -ray flux is (Ackermann et al., 2014). This leads to a CR electron-to-proton ratio , which is consistent with the value of expected from diffusive shock acceleration (Park et al., 2015). Future HAWC observations may be able to resolve this hadronic component based on its TeV emission.
Finally, we note that future instruments, such as the Large High Altitude Air Shower Observatory (LHAASO) (ZHA, 2012; Di Sciascio and on behalf of the LHAASO Collaboration, 2016), the Hundred Square km Cosmic ORigin Explorer (HiScore) (Tluczykont et al., 2012), and the Cherenkov Telescope Array (CTA) (Actis et al., 2011) will soon provide additional high-sensitivity observations of the Fermi bubbles, allowing the combined -ray constraints from these models to strongly constrain the spectrum and intensity of the -ray and neutrino emission from the bubbles.
We thank E. Blaufuss, D. Fiorino, C. Kopper, C. Rivière, D. Malyshev, and A. Smith for helpful comments. K. F. acknowledges the support of a Joint Space-Science Institute prize postdoctoral fellowship. TL acknowledges support from NSF Grant PHY-1404311. KM is supported by NSF Grant No. PHY-1620777.
- As a cross-check, we note that using the average effective area of the Southern sky, the best-fit isotropic per-flavor neutrino flux of (The IceCube Collaboration et al., 2015) predicts the observation of 35.5 events in the Southern sky. This is comparable to the actual value of 37.
- We note that IceCube has recorded more years of HESE data than have been publicly released, so the improvement in data can not be placed on a standard timeline for both HAWC and IceCube observations.
- M. Su, T. R. Slatyer, and D. P. Finkbeiner, Astrophys. J. 724, 1044 (2010), eprint 1005.5480.
- D. P. Finkbeiner, Astrophys. J. 614, 186 (2004), eprint astro-ph/0311547.
- G. Dobler and D. P. Finkbeiner, Astrophys. J. 680, 1222 (2008), eprint 0712.1038.
- G. Dobler, Astrophys. J. 750, 17 (2012), eprint 1109.4418.
- F. Guo and W. G. Mathews, Astrophys. J. 756, 181 (2012), eprint 1103.0055.
- F. Guo, W. G. Mathews, G. Dobler, and S. P. Oh, Astrophys. J. 756, 182 (2012), eprint 1110.0834.
- H.-Y. K. Yang, M. Ruszkowski, P. M. Ricker, E. Zweibel, and D. Lee, Astrophys. J. 761, 185 (2012), eprint 1207.4185.
- J. Bland-Hawthorn, P. R. Maloney, R. S. Sutherland, and G. J. Madsen, The Astrophysical Journal 778, 58 (2013), URL http://stacks.iop.org/0004-637X/778/i=1/a=58.
- K. S. Cheng, D. O. Chernyshov, V. A. Dogiel, and C. M. Ko, The Astrophysical Journal 790, 23 (2014), URL http://stacks.iop.org/0004-637X/790/i=1/a=23.
- P. Mertsch and S. Sarkar, Phys. Rev. Lett. 107, 091101 (2011), URL http://link.aps.org/doi/10.1103/PhysRevLett.107.091101.
- K. Sasaki, K. Asano, and T. Terasawa, Astrophys. J. 814, 93 (2015), eprint 1510.02869.
- U. Keshet and I. Gurwich, ArXiv e-prints (2016), eprint 1611.04190.
- S. R. Kelner, F. A. Aharonian, and V. V. Bugayov, Phys. Rev. D 74, 034018 (2006), URL http://link.aps.org/doi/10.1103/PhysRevD.74.034018.
- M. Ackermann, A. Albert, W. B. Atwood, L. Baldini, et al., ApJ 793, 64 (2014), eprint 1407.7905.
- R. M. Crocker and F. Aharonian, Physical Review Letters 106, 101102 (2011), eprint 1008.2658.
- Y. Fujita, Y. Ohira, and R. Yamazaki, Astrophys. J. 789, 67 (2014), eprint 1405.5214.
- Y. Fujita, K. Murase, and S. S. Kimura, ArXiv e-prints (2016), eprint 1604.00003.
- B. C. Lacki, MNRAS 444, L39 (2014), eprint 1304.6137.
- A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, J. D. Álvarez, R. Arceo, J. C. Arteaga-Velázquez, H. A. Ayala Solares, A. S. Barber, N. Bautista-Elivar, et al., ArXiv e-prints (2017), eprint 1703.01344.
- C. Lunardini and S. Razzaque, Physical Review Letters 108, 221102 (2012), eprint 1112.4799.
- F. Halzen, Nature Phys. 13, 232 (2016).
- IceCube Collaboration, 342 (2013), ISSN 0036-8075, URL http://science.sciencemag.org/content/342/6161/1242856.
- M. Ahlers and K. Murase, Phys. Rev. D 90, 023010 (2014), eprint 1309.4077.
- C. Lunardini, S. Razzaque, K. T. Theodoseau, and L. Yang, Phys. Rev. D 90, 023016 (2014a), eprint 1311.7188.
- A. M. Taylor, S. Gabici, and F. Aharonian, Phys. Rev. D 89, 103003 (2014), eprint 1403.3206.
- M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M. Ahrens, D. Altmann, T. Anderson, C. Arguelles, T. C. Arlen, et al., European Physical Journal C 75, 116 (2015), eprint 1409.4535.
- C. Lunardini, S. Razzaque, and L. Yang, Phys. Rev. D 92, 021301 (2015), eprint 1504.07033.
- M. Ahlers, Y. Bai, V. Barger, and R. Lu, Phys. Rev. D 93, 013009 (2016), eprint 1505.03156.
- M. G. Aartsen, M. Ackermann, J. Adams, et al. (IceCube Collaboration), Phys. Rev. Lett. 113, 101101 (2014), URL http://link.aps.org/doi/10.1103/PhysRevLett.113.101101.
- The IceCube Collaboration, M. G. Aartsen, K. Abraham, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M. Ahrens, D. Altmann, T. Anderson, et al., ArXiv e-prints (2015), eprint 1510.05223.
- G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998), eprint physics/9711021.
- M. C. Chantell, C. E. Covault, J. W. Cronin, B. E. Fick, L. F. Fortson, J. W. Fowler, K. D. Green, B. J. Newport, R. A. Ong, S. Oser, et al., Physical Review Letters 79, 1805 (1997), eprint astro-ph/9705246.
- M. Lemoine-Goumard, PoS ICRC2015, 012 (2016), eprint 1510.01373.
- M. Tluczykont, D. Horns, D. Hampf, R. Nachtigall, U. Einhaus, M. Kunnas, T. Kneiske, and G. P. Rowell, Nuclear Instruments and Methods in Physics Research A 692, 246 (2012).
- M. ZHA, International Journal of Modern Physics: Conference Series 10, 147 (2012), URL http://www.worldscientific.com/doi/abs/10.1142/S2010194512005867.
- C. Lunardini, S. Razzaque, K. T. Theodoseau, and L. Yang, Phys. Rev. D 90, 023016 (2014b), eprint 1311.7188.
- M. Ackermann, A. Albert, W. B. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, R. Bellazzini, E. Bissaldi, R. D. Blandford, et al., The Astrophysical Journal 793, 64 (2014), URL http://stacks.iop.org/0004-637X/793/i=1/a=64.
- K. Blum, A. Hook, and K. Murase, ArXiv e-prints (2014), eprint 1408.3799.
- R. Laha, J. F. Beacom, B. Dasgupta, S. Horiuchi, and K. Murase, Phys. Rev. D 88, 043009 (2013), eprint 1306.2309.
- W. I. Axford, E. Leer, and G. Skadron, International Cosmic Ray Conference 11, 132 (1977).
- G. F. Krimsky, Dokl. Akad. Nauk SSSR 234, 1306 (1977).
- A. R. Bell, MNRAS 182, 147 (1978).
- R. D. Blandford and J. P. Ostriker, ApJ 221, L29 (1978).
- J. Park, D. Caprioli, and A. Spitkovsky, Physical Review Letters 114, 085003 (2015), eprint 1412.0672.
- G. Di Sciascio and on behalf of the LHAASO Collaboration, ArXiv e-prints (2016), eprint 1602.07600.
- M. Actis, G. Agnetta, F. Aharonian, A. Akhperjanian, J. Aleksić, E. Aliu, D. Allan, I. Allekotte, F. Antico, L. A. Antonelli, et al., Experimental Astronomy 32, 193 (2011), eprint 1008.3703.