Hadronic models of the Fermi bubbles: Future perspectives

The origin of sub-TeV gamma rays detected by Fermi-LAT from the Fermi bubbles at the Galactic center is still unknown. In a hadronic model, acceleration of protons and/or nuclei and their subsequent interactions with gas in the bubble volume can produce observed gamma rays. Such interactions naturally produce high-energy neutrinos, and detection of those can discriminate between a hadronic and a leptonic origin of gamma rays. Additional constraints on the Fermi bubbles gamma-ray flux in the TeV range from recent HAWC observations restrict hadronic model parameters, which in turn disfavor Fermi bubbles as the origin of a large fraction of neutrino events detected by IceCube along the bubble directions. We revisit our hadronic model and discuss future constraints on parameters from observations in very high-energy gamma rays and neutrinos.


Introduction
The Fermi bubbles (FB), discovered in the gamma-ray data of the Fermi Large Area Telescope (LAT), are two giant structures extending up to 55 degrees (∼ 9 kpc) above and below the Galactic center (GC) [1][2][3][4]. The associated multi-wavelength observations in microwaves (WMAP haze) [5,6], X-ray [7,8] and polarized radio waves [9] provide comprehensive information for studying the physical origin of the bubbles. Several theoretical models have been proposed to explain the morphology and spectral properties of detected gamma rays, typically classified as hadronic and leptonic models. Both mechanisms can reproduce the observed hard spectrum, sharp edges and uniform emission at latitudes |b| > 10 • . In the hadronic models, gamma rays are generated by inelastic collisions of accelerated cosmic rays on thermal nuclei in the bubble gas [3,[10][11][12]. While in the leptonic models, inverse Compton scattering of relativistic electrons on optical and UV photons produce the gamma rays [2,3,13,14]. The true origin of gamma rays, whether hadronic or leptonic, has profound implications for the history of star-formation activities in the central region of the Milky Way and/or the activity of the super-massive black hole at the center, as well as particle acceleration.
In hadronic models, TeV to PeV neutrinos, resulting from charged pion/kaon decays, as counterparts of GeV gamma rays, resulting mostly from neutral pion decays, should exist [10,12]. Detection of these high energy neutrinos can serve as one of the major discriminator between the hadronic and leptonic models. It was pointed out in the past that a fraction of the high-energy astrophysical neutrinos detected by IceCube [15] could originate from the FB in hadronic models [16][17][18][19]21]. The corresponded neutrino flux from the bubbles is consistent with the hadronic flux model that reproduce FB gamma ray data. If that is the case, the bubbles could be the first multi-messenger study source in our milky way [20,21]. Recently, the High Altitude Water Cherenkov (HAWC) telescope reported their analysis of gamma-ray data from the high-latitude (b > 6 • ) Northern FB region collected over 290 days [24] and found no evident excess above ∼ 1 TeV. The HAWC upper limits are consistent with the gamma-ray flux from Fermi-LAT but constrained our previous "neutrino-inspired" hadronic models that predict high gamma-ray flux level in the ≥ 1 TeV range from the whole FB region [20].
The morphology of the FB gamma-ray emission, however, is not fully defined, especially toward the GC region, which has also not been clearly observed in other wavelength. In the latest Fermi-LAT analysis [23], it was found that the spectra of the low-latitude (|b| < 10 • ) and high-latitude (|b| > 10 • ) bubbles are similar in the 100 MeV to 100 GeV energy range, but different at energies above 100 GeV. In particular the spectrum of the low-latitude bubbles is hard and without any apparent fall-off as compared to the high-latitude bubble spectrum which shows a sharp cutoff above ∼ 100 GeV [23]. As mentioned, the HAWC upper limits [24] do not constrain the low-latitude FB spectrum. The HAWC upper limits derived for the North bubble also do not strictly apply to the South bubble in case the South bubble has different gamma-ray morphology and spectrum than the North bubble. The sample of the high-energy astrophysical neutrinos detected by IceCube has grown to 82 in the meantime (year 2010-2016) [26] and several more of these events are now from the FB region as compared to what we had reported in [18,19]. Given these developments, it is relevant now to critically review and update the FB hadronic models with relevant constraints.
In this article we present updated hadronic models for the FB, based on the recent gamma-ray spectral analysis by the Fermi-LAT Collaboration with 6.5 years data [23], using the high-and low-latitude bubble templates, latest sample of astrophysical neutrinos detected by IceCube and HAWC constraints at TeV energy gamma ray observations. A new feature in the latest Fermi-LAT data is that the gamma-rays from the FB at the high-and low-Galactic latitudes have different spectra. Therefore we present two hadronic models in the present work fitting those spectra. We take (in this work and in Refs. [18][19][20]) both the gamma-ray and neutrino data into account to constrain the primary proton spectra instead of calculating the proton distribution function fitting the gamma-ray data only (see, e.g., Ref. [22]). The neutrino data to is especially useful to constrain the high-energy cutoff in the primary proton spectrum for the low-latitude case. We discuss details of the latest gamma-ray and neutrino data in Sec. 2, develop updated hadronic models in Sec. 3 and discuss possible future constraints from gamma-ray observations in Sec. 4. We summarize and conclude in Sec. 5.

Recent gamma-ray and neutrino data
The Fermi bubbles as analyzed with Fermi-LAT data [3,23], were found to have uniform spectra above 10 • in Galactic latitude with sharp edges. In the most recent update with 6.5 years of Fermi-LAT data, the spectra and morphology of central bubbles have also been derived with the spectral components analysis procedure as in [3]. In this latest analysis, the spectra of low-and high-latitude bubbles are treated separately, but with the assumption that they behave similarly in the energy between 1 GeV and 10 GeV. In addition, the center of the bubbles is found to become brighter as they are close to the Galactic plane, which is different from previous models [25] with isotropic emission in the center region. In our study, we adopted the templates derived by Fermi Collaboration with the region of interest |b| < 60 • , |l| < 45 • , as shown in gray contours in Fig. 1, obtaining a solid angle of Ω FB =1.04 sr and fit our models with respect to the low-and high-latitude components separately as well.
Recently, the HAWC Collaboration reported their search for very high energy (VHE) gamma-ray emission from the Northern FB region with latitude above 6 • which corresponds to a FB solid angle of 0.42 sr [24]. No significant excess was found in the analysis and the 95% C.L. upper limits consistent with HAWC detection power on the differential flux in four energy bins were obtained as shown in Fig. 2. The limits agree with gamma ray flux above ∼ 2 TeV from the FB regions at high-latitudes.
Lately, the IceCube Neutrino Observatory has also updated its high-energy starting events (HESE) search with the neutrino interaction vertex inside the IceCube detector volume and energies above 30 TeV, using six-year dataset, which corresponds to a total livetime of 2078 days. In the full sample, there are 82 events detected, with expected atmospheric muon background of 25.2 ± 7.3 events [26]. Of these events, eight (Nos. 2,12,14,15,36,56,69,76) are spatially strongly correlated (the best-fit arrival direction is within the bubble geometry) and six (Nos. 17,22,24,25,49,68) are weakly correlated (the median positional error circles overlap with the bubble geometry) with the FB, as seen in Fig. 1. Except for the event No. 76, which is a track, all the other events are showers. In particular, within these eight strongly correlated neutrino events, three of them are from high-latitude northern bubble, four from high-latitude southern bubble, one with the highest energy ∼ 1 PeV is from the Galactic Center region. Fig. 1 also shows the median positional error for the neutrino arrival directions.
In Fig. 2 we show the neutrino spectrum (all three flavors) we have calculated from the 8 events in three energy bins from the FB region using the livetime and effective areas of IceCube for HESE analysis [26]. There are 3, 4 and 1 event(s), respectively, in the energy ranges (

Hadronic models
Prolonged star-formation activity near the Galactic center, which forms a strong bipolar wind above and below the Galactic plane, has been proposed as a plausible mechanism to inflate the FB [10]. Alternately, activity of the Galactic super-massive black hole in the past could also be responsible for the FB origin [2]. In the hadronic models, the gamma rays from the FB are produced primarily from neutral pion decays which are created by proton-proton (pp) interactions of energetic cosmic rays with the dilute gas inside the bubble volumes. Cosmic rays, accelerated to ∼ 10 15 eV and possibly beyond in supernova remnants, can be carried by the wind to fill the FB volume. Decays of charged pions and kaons, co-produced with neutral pions, by pp interactions guarantee a neutrino flux associated with the gamma-ray flux in the hadronic models.
In our hadronic model formalism [12,18,20] we assume a primary proton spectrum in the form of a power law with exponentially cutoff as N p (E) = N 0 E −k exp(−E/E 0 ). Here N 0 , k and E 0 are the normalization factor, spectral index and cutoff energy, respectively. The normalization factor N 0 is adjusted to fit the gamma-ray spectra by using a particle density of 10 −2 cm −3 for the gas inside the bubble volumes. The spectral index and the cutoff energy are crucial parameters which can be utilized to constrain the hadronic models from multi-wavelength and multi-messenger observations. In Fig. 2 we show our hadronic model fits for both high-latitude (FB |b| > 10 • ) and low-latitude (FB |b| < 10 • ) spectra of the FB. The high-latitude spectrum is rather well-fitted with k = 2.0 and E 0 = 1.6 TeV (blue dotted lines). This cutoff energy is consistent with the HAWC upper limits in the ∼ 2-200 TeV range, which apply dominantly to the high-latitude part (b > 6 • ) of the North bubble [24]. We show two model fits (orange dot-dashed lines) for the low-latitude spectra in Fig. 2, both fitting low-energy (≤ 10 GeV) data rather poorly. The low-energy part of the spectrum is less reliable due to uncertainties in the Galactic emission template at low latitudes and possible contamination with Galactic center emission [23], such as the point source detection, morphology study of bubbles and analysis of Galactic Interstellar Emission. The model with lower cutoff energy has parameters k = 2.15 and E 0 = 30 TeV and the one with higher cutoff energy has parameters k = 2.2 and E 0 = 3 PeV. All values of k and E 0 are consistent with shock-acceleration scenario of cosmic rays in supernova remnants.
The neutrino flux (all three flavors) arising from the hadronic model with k = 2.2 and E 0 = 3 PeV is shown in Fig. 2 with magenta dashed lines. Although such a high cutoff energy is not required by the Fermi-LAT gamma-ray data alone, an explanation of the neutrino data could be provided with this model. This is also motivated by the hard spectrum of the low-latitude FB without showing any cutoff, unlike the high-latitude spectrum, and that the HAWC upper limits do not strictly apply in this region. Of course the neutrino events that are strongly correlated with the FB should mostly be arriving from the |b| < 10 • regions in this scenario, which could be plausible given large uncertainties in reconstructing their directions [26]. Alternately, only a fraction of these 8 events could originate from the FB and the rest would be from a diffuse astrophysical background. We elaborate on this scenario in the following subsection.

Neutrino events from the FB
In the recent two-year dataset, IceCube Observatory detected 28 more HESE neutrinos, and notably all of them are with energy below 200 TeV [26]. Altogether there are now 82 events detected so far and as Nos. 20 and 55 are in coincident with background muons, they are excluded from the analysis. The likelihood fitting was performed in the deposited energy between 60 TeV and 10 PeV by IceCube [26], resulting a best fit single power law flux of dN/dE = 2.46 ×  (E/100 TeV) −2.92 GeV −1 cm −2 s −1 sr −1 per flavor, which is softer than the previous results [15,29] due to the additional events at low energies. Initially IceCube focused on the higher-energy search around PeV [15], later the selection threshold of deposited energy was reduced down to ∼ 1 TeV [27], where cosmic ray muon background is dominant. Applying the veto technique as in [27] and also the first-level online filters, resulted in an increase in the effective area [26] for both cascade-like and track-like events.
To estimate the neutrino signal excess from the FB in 2078 days operation of IceCube, the major backgrounds are atmospheric muons, neutrinos of astrophysical origin and atmospheric neutrinos in the energy range of 40 TeV to 70 PeV. The first case strongly relies on the efficiency of background rejection and reconstruction technique (see, [15,26] for more information). The predicted atmospheric ν µ +ν µ flux from [28], averaged over the declination angle of the FB as seen from IceCube, is adopted and extrapolated at high energy. Moreover, the atmospheric ν e +ν e fluxes are greatly suppressed compare to ν µ +ν µ by a factor of 14 [30], which has also been taken into account. The astrophysical neutrino fluxes from the IceCube best fit adopted here are shown in Fig. 3(a) as blue line with one-sigma uncertainties (shaded bands), and the hadronic fluxes from the FB region as in Fig. 2 is shown as well. All plotted fluxes in Fig. 3(a) are for a solid angle of FB Ω FB =1.04 sr for all three flavors (combined neutrino and antineutrino).
Since all the neutrinos spatially correlated with FB are cascade events, except the No. 76, we calculate the number of cascade events only with updated HESE effective areas as seen in Fig. 3.
The number of events in each energy bin i can be calculated as where A eff and T exp are ν e or ν τ effective area and exposure time (2078 days), respectively. These event distributions are plotted in Fig. 3. We expect a total of 8.8 and 7.9 cascade-like events from the backgrounds (atmospheric + astrophysical) and FB, respectively, in the first five energy bins. Given the harder spectrum of the FB flux than the diffuse astrophysical flux, the signature of the FB is more prominent in the ∼ 100 TeV-1 PeV range. Within one sigma uncertainty, the number of background events can be reduced to 5.5. The estimation agrees with the IceCube observation, where 8 (6) events are strongly (weakly) correlated with the FB. One note of caution, however, is that the diffuse flux fit is based on all events, including plausible contribution from the FB.
Detection of track-like neutrino events, which have ≤ 1 • angular resolution, by the upcoming KM3NeT Neutrino Telescope located in the northern hemisphere [33] will be helpful to separate the FB from other possible Galactic sources, as well as complimentary to the IceCube data dominated by shower-like events.

Future VHE gamma-ray constraints on hadronic models
The HAWC observatory found no significant excess towards the FB using the data between 2014 November 27th to 2016 February 11th, resulting in a 90% upper limits on high-latitude northern bubble with b > 6 • in the energy between 1 and 100 TeV [24]. These upper limits agree with our hadronic emission model of the high-latitude FB, which is strongly suppressed beyond ∼ 1 TeV (see Fig. 2). The spectrum of the low-latitude FB is harder and without any apparent cutoff in the energy range of the Fermi-LAT. Therefore both our nominal model with parameters k = 2.15 and E 0 = 30 TeV and the neutrino-inspired model with parameters k = 2.2 and E 0 = 3 PeV for the low-latitude spectrum are potential targets for future VHE observations. We expect that the HAWC collaboration will further analyze data from the FB in future, including the sensitivity at lower energies and larger field of view that includes the Galactic Center region. This will provide more stringent constraints on both the spectra at low-and high-latitude bubbles, the full image of non-uniform intensity and the shape of the central region.
On the other hand, there is good hope that the Cherenkov Telescope Array, the next generation of ground-based VHE observatory [31] with a northern (La Palma) and a southern (Chile) site, will be able to observe both bubbles. A detailed strategy needs to be developed about pointing different parts of the FB as well as better controlling the background. In Fig. 4 we show the differential sensitivity of the CTA-South for 50 hour observation of a point source [31], for illustration purpose. Moreover, the Large High Altitude Air Shower Observatory (LHAASO) being built at 4410 m altitude in Sichuan Province of China, detecting cosmic rays and gamma rays in the energy range of 100 GeV to 1 PeV [32], will be complementary to the HAWC and CTA. With a high duty cycle and wide filed of view, ≈ π sr, LHAASO will be able to observe the full northern bubble with a solid angle of 0.45 sr. We show the differential sensitivity of LHAASO in Fig. 4 for one year observation of a Crab-like point source [32], again for illustration purpose. A detailed spectral and morphological study with simulations will be required to calculate the sensitivities of the CTA and HAWC to the FB. The differential point source flux sensitivities of the CTA southern site [31] with 50 hour observation (gray dashed curve) and of the LHAASO [32] with 1 year observation (red dot-dashed curve) are shown with the FB flux models for illustration purpose.

Summary and Outlook
In this paper we have carried out updated hadronic modeling of the 0.1 GeV -1 TeV gamma-ray emission from the FB, based on new data from the Fermi-LAT, upper limits from the HAWC and expanded astrophysical neutrino sample from the IceCube. Latest analysis of Fermi-LAT data prefers a spectrum with cutoff above 1 TeV in the high Galactic latitude (|b| > 10 • ) region of the FB while a hard spectrum without any apparent cutoff in the low Galactic latitude (|b| < 10 • ) region of the FB [23]. Combined with the latest IceCube astrophysical neutrino data [26] our updated hadronic model can be extended to the 1 PeV range for the low-latitude FB emission and is not directly constrained by the HAWC upper limits derived from observations of the northern bubble at high-latitude (b > 6 • ) [24]. Our updated hadronic model for the high-latitude FB spectrum is consistent with the HAWC upper limits.
The observation of high-energy neutrinos from the FB would strongly support the hadronic origin of gamma rays. Detection of several more astrophysical neutrino events by IceCube from the direction of the FB in the latest dataset is therefore intriguing. Our estimated FB neutrino flux based on these events is harder than a diffuse astrophysical flux estimate by the IceCube collaboration, although it is not possible to distinguish them using the current data set. Future observation of neutrinos by KM3NeT [33] from the FB region can shed insights, especially if many track-like events are detected with good angular resolution. Observations of VHE gamma rays from the FB with upcoming CTA [31] and LHAASO [32], combined with deeper observations by HAWC, will be crucial to probe the spectra of FB in the 1-100 TeV range that ties with the lower energy range of the IceCube neutrino events. The detection of a cutoff in this energy range for the low-latitude FB spectrum could strongly constrain hadronic models that can explain the IceCube neutrino data and critically test the FB as a multi-messenger source.