© ESO 2019

1 Introduction

The Fermi bubbles (FBs) are one of the most spectacular and unexpected discoveries in the Fermi Large Area Telescope (LAT) data (Su et al. 2010). The FBs extend up to 55° in latitude above and below the Galactic center (GC), they have a well-defined edge and a relatively uniform intensity across the surface, apart from a “cocoon” in the southeastern part of the bubbles (Su & Finkbeiner 2012; Ackermann et al. 2014). The intensity spectrum is ~E⁻² at gigaelectronvolt energies with a cutoff or a softening around 100 GeV at latitudes |b| > 10° (Ackermann et al. 2014). The origin of the FBs is attributed either to an emission from the supermassive black hole (SMBH) at the center of our Galaxy or to a period of starburst activity that resulted in a combined wind from supernova (SN) explosions of massive stars (Su et al. 2010). The gamma-ray signal up to 100 GeV can be produced either by interactions of hadronic cosmic rays (CRs) with gas (hadronic model) or by inverse Compton (IC) scattering of high-energy electrons and positrons and interstellar radiation (leptonic model).

Although the FBs were detected about eight years ago, their origin is still unresolved. Important insights into their origin can be obtained from the study of the morphology and the spectrum of the FBs near the GC. The spectrum of gamma rays can provide information on the composition of the CRs that produce the gamma-ray signal (hadronic versus leptonic CRs), as well as the age of the CRs (through a cooling cutoff in the leptonic scenario or a break from the escape of high-energy CRs) and the spectrum of the CRs at the source. The morphology of the emission can point to the source of the bubbles: either the SMBH Sgr A* or a recent star-forming region. Previous analyses of the FBs at low latitudes indicated higher intensity of emission near the Galactic plane (GP) and a displacement to negative longitudes (Acero et al. 2016; Ackermann et al. 2017; Storm et al. 2017). The spectrum of the FBs for |b| < 10° is consistent with a power-law ∝ E⁻² without a cutoff up to 1 TeV (Ackermann et al. 2017).

The study of the FBs in the GP is complicated owing to bright Galactic diffuse emission components. The π⁰ component of the gamma-ray emission is well traced by the distribution of gas. But this component has large uncertainties toward the GC because of a lack of kinematic information from the motion of the gas in the Galaxy, which is used to reconstruct the gas distribution; the velocity of the gas in the direction of the GC is perpendicular to the line of sight. These uncertainties are also a consequence of the uncertainties in CO emission along the line of sight to the GC, which is used as a tracer of molecular hydrogen, and large dispersion of velocities of some molecular clouds near the GC. There are also large uncertainties in the distribution of the CR sources and the propagation model near the GC, which make it rather difficult to predict a priori the distribution of the propagated CRs in the Galaxy. For the IC component of the gamma-ray emission, there is a significant uncertainty in the interstellar radiation field (ISRF) density nearthe GC (Popescu et al. 2017; Porter et al. 2017; Niederwanger et al. 2019) in addition to the uncertainties in the CR distribution. There should also be undetected point-like and extended sources, which nevertheless contribute to the total flux. The distribution of CRs in the Galaxy can be computed with CR propagation tools, such as GALPROP (Strong et al. 2007). The maps of gamma-ray emission from interactions of CRs with gas in different Galactocentric rings and from interaction of CR electrons and positrons with the ISRF can be used as templates for the corresponding components of gamma-ray emission. The agreement of the gamma-ray data with modelsbased on templates derived with the CR propagation tools is rather poor in the GP (e.g., Porter et al. 2017).

In this paper, we analyze the gamma-ray emission at the base of the FBs and estimate the uncertainties related to modeling of the Galactic foreground and background components. We focus on the morphology and spectrum of the FBs at energies from 10 GeV to 1 TeV, for which the intensity of the gamma-ray emission from the FBs relative to the other Galactic components is higher than at low energies as a consequence of softer spectra of the Galactic components. The study of the FBs at high energies will be also important for future searches with neutrino and Cherenkov telescopes, such as the Cubic Kilometer Neutrino Telescope (KM3NeT) and the Cherenkov Telescope Array (CTA).

We use several methods to determine the Galactic foreground and background emission. Unfortunately the notion of the foreground and background emission may not be well defined if we search for an extended gamma-ray emitting source based on the gamma-ray data itself rather than using multiwavelength data, which are used, for example, to trace the distribution of gas. In this work, by foreground and background emission we mean a steady state (or average) diffuse emission of gamma-rays. The steady state distribution of CRs is obtained by averaging in time over many sources; it results in the local CR proton spectrum of ~ E^−2.7 for energies from ~5 GeV to ~1 PeV and ~ E^−3.3 for the local CR electrons from ~1 GeV to ~1 TeV energies. If the spectrum of the CRs at the source is ~E^−2.0−2.2, then the softening of the spectrum of CR protons by ~E^−0.3 − E^−0.6 is due to energy-dependent escape from the Galaxy, while the softening for the CR electrons by ~ E⁻¹ is due to cooling. A distinguishing property of a source of CRs is that the spectrum of CRs at or near the source is much harder than the average (propagated) spectrum of CRs. Thus, it is possible to look for the presence of a population of freshly accelerated CRs by searching for areas of gamma-ray emission with a spectrum harder than average. In particular, one can use an “on-off” analysis to subtract the steady state component of the gamma-ray emission. One of the caveats of this analysis is that even the steady state gamma-ray emission in the “on” region can be more intense than the emission in the “off” region, but in this case the difference of the fluxes would have a soft spectrum characteristic of the propagated CRs, unless the difference in flux is much smaller than the flux in both on and off regions. In this case the difference can have a harder spectrum than the two terms. We use the on-off technique as a preliminary check in a search for a population of freshly accelerated CRs: the hard spectrum of the difference is a necessary condition for the presence of a population of CRs with a spectrum harder than the steady state distribution of CRs. It is also a sufficient condition for the presence of a population of CRs with a hard spectrum if the difference has intensity comparable to the overall intensity in the on and off regions.

Consequently, considering that there is a tentative asymmetry in diffuse gamma-ray emission near the GC at high energies with a spectrum harder to the west of the GC relative to the emission to the east of the GC, as a first step, we estimate the amount of the asymmetry by taking the difference in the gamma-ray data to the west and east of the GC after masking bright point sources (PS). Second, since the spectrum of the FBs at high latitudes and the inferred spectrum at low latitudes are harder than the spectra of the other components of diffuse emission, the contribution of the FBs at energies ≲ 1 GeV is relatively small. We use the data below 1 GeV as a template of the Galactic foreground emission, which we fit to the data at energies above 1 GeV. The residual emission is used to estimate the contribution from the components with spectra harder than the typical spectra of the steady state diffuse gamma-ray components. Finally, we determine a model for the Galactic gamma-ray diffuse emission using templates for the emission components based on GALPROP calculation¹ (Moskalenko & Strong 1998; Strong et al. 2000, 2004, 2007; Ptuskin et al. 2006; Porter et al. 2008; Vladimirov et al. 2011). In this model, we allow many free parameters, such as rescaling of the π⁰ and bremsstrahlung emission in Galactocentric rings and refitting of bright PS near the GC. With the many free parameters, this model absorbs as much of the gamma-ray emission at the base of the FBs as it can. As a result, the residual emission at the base of the FBs in this model is smaller than the residuals in the other models considered in the paper.

The similarity of the energy spectrum below 100 GeV of the hard and bright emission at the base of the FBs and the FBs at high latitudes, as well as the spatial coincidence of the emission with a continuation of the FBs from high latitudes, make the physical correspondence of the two objects very plausible. Nevertheless, their alignment along the line of sight can be accidental and the distance to the two objects can be different, i.e., the FBs may be above and below the GC while the hard component may be at a closer distance, for example 1 kpc, from us so that the two objects would be physically unrelated. Thus, we refer to the hard and bright component as emission at the base of the FBs to keep both possibilities open: that the emission is the base of the FBs and that the emission is at a different location along the line of sight toward the base of the FBs. We discuss possible origins of the hard component of emission in Sect. 5. Section 6 contains conclusions.

2 Data selection

We used nine years of the Fermi-LAT Pass 8 Source class events between August 4, 2008 and August 3, 2017 (Fermi Mission Elapsed Time 239557418 s–523411376 s) with energies between 316 MeV = 10^2.5 MeV and 1 TeV separated in 21 logarithmic energy bins (six bins per decade). The selection of the events was performed with the standard quality cuts (DATA_QUAL >0)&&(LAT_CONFIG==1). In order to avoid contamination from gamma rays produced in interactions of CRs in the Earth atmosphere, we selected events with zenith angles θ < 100°, which is sufficient for energies above 316 MeV. We calculated the exposure and point-spread function (PSF) using the Fermi-LAT Science Tools package version 10-01-01 available from the Fermi Science Support Center² with the P8R2_SOURCE_V6 instrument response functions. Figure 1 shows the energy flux of the gamma-ray emission integrated in three energy ranges³. For all mapsshown in this paper, we used Galactic coordinates centered on the GC in Mollweide projection. For spatial binning we used the HEALPix ⁴ (Górski et al. 2005) scheme with a pixelization of order 7 (≈ 0. °46 pixel size).

Fig. 1 Energy flux of the Fermi-LAT 9 yr of Source class data integrated in energy ranges 1–10 GeV, 10–100 GeV, and 100 GeV–1 TeV. The maps are represented in logarithmic scale. The map in the energy range 10–100 GeV (100–1000 GeV) is smoothed with a Gaussian kernel with radius σ =0. °5 (σ =0. °7). The plots are created with the HEALPix package (Górski et al. 2005) in Mollweide projection with nside = 128 (≈ 0. °46 pixel size).

3 Modeling of the Fermi bubbles at low latitudes

3.1 East-west asymmetry in the data

In order to investigate the asymmetry of the emission at the base of the FBs relative to the GC, we compared the Fermi-LAT data east and west of the GC. We masked the 200 point sources (PS) in the Third Fermi-LAT source catalog (3FGL; Acero et al. 2015) with the largest flux above 1 GeV. Each PS is masked with a circle of radius $\frac{\delta }{\sqrt{2}}+$0. °5, where δ =0. °46 is the characteristic size of the pixels so that $\delta /\sqrt{2}$ is half of the pixel diagonal (if it were a square), while 0. °5 is comparable to the 68% containment radius of the PSF above ≈ 2 GeV. Thus the total radius is about 0. °82. In order to avoid possible bias by masking more pixels on one side of the GC in the GP, we symmetrized the PS mask relative to the GC, i.e., we set m_−b,−ℓ = 0 if m_b,ℓ = 0 and vice versa. We also use this PS mask in the following sections (apart from Sect. 3.4). The data are averaged over regions to the east (0° < ℓ < 10°) and to the west (− 10° < ℓ < 0°) of the GC at different latitudes. The regions have a latitude width of 4°. The fraction of masked pixels is about 50% within |b| < 2°, about 10% for 2° < |b| < 6°, and less than 5% for |b| > 6°.

The difference of the averaged intensity of emission west minus east of the GC is shown in Fig. 2. The error bars represent the statistical errors. The emission for latitudes b ∈ (−6°, −2°) and b ∈ (−2°, 2°) shows excess emission to the west of the GC, which remains significant at high energies. The Fermi-LAT exposure within |b| < 10° and |ℓ| < 10° is rather uniform; the maximal fractional variation of the exposure in this region for E > 1 GeV is less than 4%.

3.2 Low-energy data as a background model

In the previous section we have shown that the west-minus-east difference in the Fermi-LAT data in the GP relative to the GC has a hard spectrum ~E⁻² up to 1 TeV. One of the simplest ways to separate the soft component of emission from the hard component of emission is to use the low-energy data as a model of the soft component and subtract those data from the data at higher energies. Gamma rays produced in interactions of the steady state Galactic population of CRs with gas and ISRF, i.e., the π⁰, bremsstrahlung, and IC components, dominate the gamma-ray emission in the GP in the energy range E ≲ 1 GeV. Consequently, the low-energy Fermi-LAT data are a good tracer for the soft components of the diffuse gamma-ray emission in the GP and can be used to create a spatial template for the Galactic foreground. The 68% containment for the Fermi-LAT photons between 316 MeV and 1 GeV (averaged with an E⁻² spectrum) is about 1. °5, which is much larger than the sub-degree angular resolution above 1 GeV. In order to compensate for the difference in the angular resolution, we smoothed the data in each high-energy bin above 1 GeV with a Gaussian kernel of σ = 1°, which corresponds to 68% containmentangle of 1. °5 in two dimensions⁵.

We separate the whole sky in latitude stripes of 4° width where the HEALPix pixels are assigned to stripes that contain the centers of the pixels. We parametrized the model in each stripe and each energy bin independently, i.e., in the latitude stripe b for energy bin E and HEALPix pixel with index i the model consists of the following two terms: $$M_b(E,i)=k_b(E)\cdot {\tilde{N}}_b^{\text{low}}(E,i)+c_b(E)\cdot \tau (E,i),$$(1)

where k_b(E) and c_b (E) are coefficients to be determined from the fit to the data. The term c_b(E) ⋅ τ(E, i) is proportional to the exposure τ(E, i): it accounts for the isotropic extragalactic background and partially compensates for the latitude-dependent IC emission. The term $k_b(E)\cdot {\tilde{N}}_b^{\text{low}}(E,i)$ is proportional to the low-energy photon counts summed over n_low = 3 energy bins between 316 MeV and 1 GeV rescaled by the ratio of exposures at low and at high energies, i.e., $${\tilde{N}}_b^{\text{low}}(E,i)=\frac{1}{n_{\text{low}}}\left({\displaystyle \sum _{\alpha =1}^{n_{\text{low}}}}\frac{N(E_{\alpha },i)}{\tau (E_{\alpha },i)}\right)\tau (E,i),$$(2)

where N(E_α, i) is the number of photons in energy bin E_α and in pixel i.

We determine the parameters c_b(E) and k_b(E) by fitting the model to the Fermi-LAT data in energy bins E > 1 GeV. We maximize the log likelihood of the gamma distribution with respect to parameters c_b (E) and k_b (E) using the Python iminuit optimizer⁶. In order to avoid an overcompensation for the flux from the FBs, the region − 20°< ℓ < 20° is excluded from the fit, i.e., the fit is performed for the total remaining length of the stripe 20° < ℓ < 340°. In the following, the region − 20° < ℓ < 20° is referred to as the region of interest (ROI) of the FBs. Bright PSs are masked with the mask described in Sect. 3.1.

After we fit the model in each latitude stripe, we interpolated it inside the ROI of the bubbles. The intensity of the model energy flux integrated in three energy ranges is shown in Fig. 3, while the residuals after subtracting the model from the data are presented in Fig. 4. The FBs are clearly visible in the first two energy ranges, E = 1−10 GeV and E = 10−100 GeV. For E = 100 GeV−1 TeV the statistics are low, but we can still see an excess near the GP. There are other regions of large residuals in the GP and at high latitudes. These are diffuse and point-like sources with spectra harder than the average spectra of Galactic and extragalactic sources. In particular, there is a residual in diffuse emission along the GP within |ℓ| ≲ 60° at energies below 100 GeV. This residual can be due to IC emission, which has a harder spectrum than the average spectrum of the π⁰ component; to a harder spectrum of CRs in the inner Galaxy (Gaggero et al. 2015; Acero et al. 2016; Yang et al. 2016); or to a population of hard PSs in the GP in the inner parts of the Galaxy. At energies above 100 GeV the residuals are resolved into localized sources, some of which are extended. Apart from the Galactic foreground and background components, there are gamma rays emitted by the interactions of CRs with the Sun and its radiation field (Abdo et al. 2011b) and the Moon (Ackermann et al. 2016). The intensities of the emission from the Sun and the Moon are about an order of magnitude less than the intensity of the emission from the FBs. As a result, these components are not visible on the residual maps at energies above 10 GeV.

Fig. 2 Difference west minus east in the Fermi-LAT intensity relative to the GC after masking bright PSs. ThePS mask is symmetrized relative to the GC. The west (east) region is defined between − 10° < ℓ < 0° (0° < ℓ <10°). Gray lines show spectra for latitudes |b| > 10°.

Fig. 3 Energy flux of the diffuse emission model based on the low-energy data (Sect. 3.2) integrated over three energyranges.

Fig. 4 Energy flux of the residuals of the low-energy model derived in Sect. 3.2 integrated over three energy ranges.

3.3 Rectangles model of the bubbles

In the previous subsection, we excluded the ROI of the FBs. In this subsection, we perform the fit over the whole sky and model the emission from the FBs using rectangles. In order to explore the east-west asymmetry of the FBs, we introduce two rectangular templates in each latitude stripe b and energy bin E: one east, ℓ ∈ (0°, 20°), and one west, ℓ ∈ (−20°, 0°), of the GC. The width of the rectangles is 4°, i.e., the same as the width of the latitude stripes. We use the same foreground model as in Sect. 3.2 plus the isotropic template. Overall, the model has four terms in each energy bin as follows: $$\begin{array}{ll}{M}_b(E,i)=\hfill & k_b(E)\cdot {\tilde{N}}_b^{\text{low}}(E,i)+c_b(E)\cdot \tau (E,i)\hfill \\ \hfill & +R_b^{\text{east}}(E)S_b^{\text{east}}(i)+R_b^{\text{west}}(E)S_b^{\text{west}}(i),\hfill \end{array}$$(3)

where the scaling parameters k_b(E) and c_b(E) and the FBs model parameters $R_b^{\text{east}}(E)$ and $R_b^{\text{west}}(E)$ are determined independently in each 4° latitude stripe and in each energy bin E > 1 GeV by fitting the model to the smoothed Fermi-LAT data (see Sect. 3.2 for details on the smoothing of the data). The values $S_b^{\text{east}}(i)$ and $S_b^{\text{west}}(i)$ are step function templates, which are equal to 1 in each latitude stripe b for 0° < ℓ < 20° and 0° > ℓ > −20°, respectively, and 0 otherwise. The energy fluxes of the model and the residuals integrated in the range E = 10−100 GeV are shown in Fig. 5.

Fig. 5 Energy flux of the diffuse model excluding the rectangles model of the FBs (top left), the rectangles model of the FBs (top right), the residual after subtracting the total model (bottom left), and residuals plus the rectangles model of the FBs (bottom right) derived in Sect. 3.3 integrated between 10 and 100 GeV.

3.4 GALPROP model of the foreground and PS refitting

In this section, we fit the gamma-ray data using a model derived with the GALPROP Galactic CR propagation code v54.1. The model for the diffuse emission components is the same as the “sample” model in Ackermann et al. (2017). In particular, we used the GALPROP code to calculate templates for the Galactic diffuse emission components. We determine five templates in each energy bin for gamma rays produced in interactions of hadronic CRs with gas and bremsstrahlung emission corresponding to five Galactocentric rings: 0–1.5, 1.5–3.5, 3.5–8, 8–10, and 10–50 kpc. We assume propagation halo height of 10 kpc and radius of 20 kpc and spin temperature 150 K. We use three IC templates corresponding to the three ISRF components: cosmic microwave background (CMB), infrared (IR) emission of dust, and starlight (SL). In addition to GALPROP templates, we used a flat template for the FBs at latitudes |b| > 10° (Ackermann et al. 2014). We model the Loop I feature using a geometric template (e.g., Fig. 2 of Ackermann et al. 2014) based on a polarization survey at 1.4 GHz (Wolleben 2007). Templates for the gamma-ray emission from the Sun and the Moon (Orlando & Strong 2008; Abdo et al. 2011b; Johannesson et al. 2013; Ackermann et al. 2016) were obtained with the Fermi Science Tools⁷. We also include, as a GC excess template, dark matter annihilation with a generalized Navarro-Frenk-White (gNFW) profile with index γ = 1.25 (Goodenough & Hooper 2009; Abazajian et al. 2014; Calore et al. 2015) and scaling radius r_s = 20 kpc, as well as the isotropic template to model the non-resolved extragalactic sources.

For the PS, we added sources from the 3FGL catalog (Acero et al. 2015) to form a template in each energy bin. We created independent templates for the Large Magellanic Cloud and the Cygnus region. We also added the remaining extended sources in a separate template. The spatial templates for the extended sources are provided by the 3FGL catalog. The PS mask in this section is different from the mask that we used in the previous sections. In this case, we mask the 200 sources with largest flux above 1 GeV outside a 10° circle from the GC. Each source is masked with a circle of $\frac{\delta }{\sqrt{2}}+1^{°}$ radius (δ = 0. °46 is the characteristic size of the pixels), which is larger than the radius of $\frac{\delta }{\sqrt{2}}+$0. °5 used in theprevious sections. We do not mask PS within 10° from the GC. Instead, we include independent templates for 40 sources with largest flux between 10 and 100 GeV taken from the column “Flux10000_100000” in the 3FGL catalog. Each source is fit independently to each energy bin, i.e., we do not use the spectral information from the 3FGL catalog. In Fig. 6 we show the residual emission plus the FBs model and the gNFW model summed over energies between 10 and 100 GeV. Owing to the model having more free parameters in the GP relative to the models considered in the previous sections, there are fewer positive residuals in the GP. There are, however, negative residuals in several locations, which indicates errors in the model. In particular, there is a rather significant negative residual to the east of the GC. Nonetheless, the FBs are clearly visible along with the bright emission at the base of the FBs, which also displays a shift to negative longitudes relative to the GC.

Fig. 6 Energy flux of the background model based on rescaled GALPROP templates (left) and of the residuals plus the FBs and thegNFW component (right) derived in Sect. 3.4 integrated between 10 and 100 GeV.

Fig. 7 Latitude profiles of the energy flux between 10 and 100 GeV for the total data excluding the PS mask and for the FBs models for different foreground diffuse emission models to the east (left panel) and to the west (right panel) of the GC.

4 Morphology and spectrum of the gamma-ray emission at the base of the FBs

4.1 Latitude and longitude profiles

In order to quantitatively estimate the difference in the intensity of the emission of the FBs at high and at low latitudes, as well as the asymmetry of the gamma-ray emission at the base of the FBs, we plot the latitude profiles to the east and west of the GC. In Fig. 7 we show the profiles of the energy flux between 10–100 GeV as functions of the Galactic latitude for different models of Galactic foreground emission. For the low-energy model, the profiles show the residual gamma-ray emission; for the rectangles model the profiles show the intensity of the rectangles templates model. In case of the GALPROP model, we show the sum of the residual emission, the FBs template, and the GC excess template, as described in Sect. 3. For comparison we show the latitude profile of the total gamma-ray data (excluding pixels in the PS mask). For positive longitudes, some models overpredict the gamma-ray data near the GP, which leads to negative residuals. For negative longitudes, there is an increase of flux at the base of the FBs by at least a factor of 2 to 6 for all models relative to the gamma-ray emission from the FBs at high latitudes. The results for the GALPROP model differ from the low-energy and rectangles model in the GP by a factor of 2 to 3. This is due to additional freedom in the GALPROP model related to the usage of several templates in the GP.

In Fig. 8 we show the longitude profiles of the residual flux for the background model based on low-energy data. The excess at ℓ ≈−13° for |b| < 2° is probably due to residual emission around the SNR RX J1713 (Abdo et al. 2011a). Consequently, we do not considerthis source as a part of the emission at the base of the FBs. We associate the excess between − 8°≲ ℓ ≲−2° to the base of the FBs. In this paper, we use ℓ ∈ (0°,   − 10°) and |b| < 6° in the calculations of the spectrum of the gamma-ray emission at the base of the FBs.

Fig. 8 Longitude profiles of the energy flux between 10 and 100 GeV and between 100 GeV and 1 TeV for the residuals in the low-energy data background model (Sect. 3.2).

Fig. 9 Gamma-ray data SEDs of the residuals (LowE model), rectangles model of the FBs (rectangles model), residuals plus FBs model plus gNFW DM annihilation template (GALPROP model), and of the difference in the data west minus east of the GC (excluding pixels in the PS mask).

4.2 Comparison of the spectra at different latitudes

In this section we quantify the hardening of the gamma-ray spectrum at the base of the FBs. We first compare the spectral energy distribution (SED) of the emission at the base of the FBs for different foreground models in Fig. 9. The differential flux is averaged over regions to the east, ℓ ∈ (0°,  10°), and to the west, ℓ ∈ (−10°,  0°), of the GC in a thin stripe covering the GP b ∈ (−2°,  2°). For comparison, we show the total data (excluding pixels masked by the PS mask) and, on the plot for ℓ ∈ (−10°,  0°), the difference inthe data west minus east of the GC.

For negative longitudes, all models give similar results. The differential flux of the GALPROP model is smaller than the differential flux of the other models above 10 GeV, which is consistent with the profile plots in Fig. 7. The difference ofthe data west minus east of the GC is similar to the fluxes at the base of the FBs in the low-energy model and in the rectangles model. The spectra at positive longitudes show large oversubtractions and softer spectra.

To compare the behavior of the energy spectra at high energies for different latitudes, we fit a log-parabola plus the fixed foreground model counts to the total smoothed gamma-ray counts in each latitude stripe using the likelihood based on the gamma distribution. We use the following parametrization of the log-parabola function: $$f(E)=N_0{\left(\frac{E}{1\text{\hspace{0.17em}GeV}}\right)}^{-\alpha -\beta \ln (E/1\text{GeV})}.$$(4)

The local “index” of the spectrum at energy E is $$n\equiv -\frac{\text{d\hspace{0.17em}}\ln f}{\text{d\hspace{0.17em}}\ln E}=\alpha +2\beta \ln \left(\frac{E}{1\text{\hspace{0.17em}GeV}}\right).$$(5)

In Fig. 10 we show this log-parabola index n as a function of latitude at E = 500 GeV. We plot (2 − n), which corresponds to the SED index. For positive longitudes, the index is relatively soft (n > 2) for most of the latitudes except high latitudes where the gamma-ray statistics are small. For negative longitudes (i.e., west), the index near the GC is significantly harder (n ≈ 2) than the index at higher latitudes.

Fig. 10 Index of the log-parabola defined in Eq. (5) at E = 500 GeV as a function of latitude for the FBs in the different foreground diffuse emission models.

4.3 Parametric model of the gamma-ray spectrum at low latitudes

In this section we study the spectrum at the base of the FBs at latitudes |b| < 6°. As a baseline model we use the rectangles model of the FBs. We compare a power-law model of the energy spectrum with a power-law and an exponential cutoff model $$\frac{\text{d}N}{\text{d}E}=N_0{\left(\frac{E}{1\text{GeV}}\right)}^{-n}{e}^{-E/E_{\text{cut}}}.$$(6)

We fit the spectral model added to the fixed background to the total gamma-ray data counts. The best-fit parameters for the rectangles model are reported in Table 1. If the improvement in $2\Delta \text{log}\mathcal{L}$ is less than 1, then we report the power-law model parameters; it turns out that the improvement in $2\Delta \text{log}\mathcal{L}$ is less than 0.1 in the cases when the cutoff is not significant. For the model with a cutoff, we also find the 95% statistical confidence lower limit for the cutoff energy. In the last column we report the minimum among the 95% confidence lower limits for the three models considered in the paper (the low-energy data model, the rectangles model, and the GALPROP templates model). The spectrum of the residuals including the FBs at ℓ ∈ (−10°,  0°) is consistent with a power law with an index 2.1−2.2 up to E ≈ 1 TeV. The 95% confidence lower limit on the cutoff energy is ~500 GeV for latitudes |b| < 6° within − 10° < ℓ < 0°.

As seen from the spectra in Fig. 9, there are large (relative to statistical uncertainty) fluctuations in some of the models. These fluctuations are due to modeling uncertainty of the Galactic diffuse emission. As a result, a fit with a simple function may be sensitive to the starting point (there can be several local minima). In the GALPROP model of the foreground, we also had to constrain the range of the index of the spectrum of the FB between 1.5 and 2.5 in the derivation of the 95% confidence lower limit on the cutoff energy.

4.4 Inverse-Compton model of the gamma-ray emission

In this section we model the gamma-ray emission at the base of the FBs with an IC scattering model. The SED of the gamma-ray intensity is parametrized as $$E_{\gamma }^2\frac{\text{d}{F}_{\gamma }}{\text{d}{E}_{\gamma }}=\frac{E_{\gamma }c}{4\pi }{\displaystyle \int \frac{\text{d}{n}_{h\nu }}{\text{d}{E}_{h\nu }}}{\sigma }_{\text{IC}}\frac{\text{d}{\Sigma }_{\text{e}}}{\text{d}{E}_{\text{e}}}\text{d}{E}_{h\nu }\text{\hspace{0.17em}d}{E}_{\text{e}},$$(7)

where $\frac{\text{d}{\Sigma }_{\text{e}}}{\text{d}{E}_{\text{e}}}={\displaystyle \int \frac{\text{d}{n}_{\text{e}}}{\text{d}{E}_{\text{e}}}}\text{d}R$ is the column density of CR electrons integrated along the line-of-sight distance R, $\frac{\text{d}{n}_{h\nu }}{\text{d}{E}_{h\nu }}$ is the number density of ISRF photons of energy hν, and σ_IC (E_γ, E_hν, E_e) is the differential IC scattering cross sectionin units of $E_{\gamma }\frac{\text{d}\sigma }{\text{d}{E}_{\gamma }}$ (Blumenthal & Gould 1970). For details on the parametrization of σ_IC see Appendix B of Ackermann et al. (2014). The ISRF number density for SL and IR photons is taken from Porter et al. (2008), which is available with the distribution of GALPROP v54.1, assuming that the emission at the base of the FBs is near the GC. For the CMB, we use the thermal spectrum with a temperature of 2.73 K. There is a significant uncertainty in the Galactic ISRF energy density near the GC (e.g., Popescu et al. 2017; Porter et al. 2017; Niederwanger et al. 2019). We find (Appendix C) that this uncertainty can lead to a factor ~3 uncertainty in the overall normalization of the inferred CR electrons (CRe) energy density. Additional uncertainty on the CRe distribution and the associated IC emission comes from the distribution of the Galactic magnetic field (e.g., Orlando 2018, 2019). We model the column density of electrons as a power law with a cutoff $$\frac{\text{d}{\Sigma }_{\text{e}}}{\text{d}{E}_{\text{e}}}=n_{\text{e}}{\left(\frac{E_{\text{e}}}{1\text{\hspace{0.17em}GeV}}\right)}^{-{\gamma }_{\text{e}}}{e}^{-E_{\text{e}}/E_{\text{cut}}}.$$(8)

In order to determine the normalization n_e, the spectral index γ_e, and the cutoff E_cut, we fit the IC model of the FBs plus the (fixed) foreground model counts to the total smoothed gamma-ray counts in different latitude stripes using likelihood based on the gamma distribution. As a baseline case, we take the gamma-ray spectrum derived in the rectangles model of the FBs in Sect. 3.3. The best-fit spectra for the rectangles model of the bubbles are shown in Fig. 11 in latitude stripes b ∈ (2°, 6°), b ∈ (−2°, 2°) and b ∈ (−6°, −2°). If the improvement in $2\Delta \text{log}\mathcal{L}$ with and without the cutoff is less than 0.1, then we show only the parameters for the power-law model without a cutoff. For example, fornegative longitudes the cutoff is not significant.

The 95% statistical lower limit on the cutoff energy for the rectangles model and the minimum among the three models presented in Sects. 3.2–3.4 of the 95% confidence lower limits for the cutoff energy are presented in Table 2. For negative longitudes, the 95% confidence lower limit on the cutoff in the spectrum of electrons in the rectangles model is about 4 TeV, while the minimal value of the 95% confidence limit for the three models of the foreground emission is about 3 TeV.

4.5 Hadronic model of gamma-ray emission

In the hadronic model, the gamma rays are produced as a result of collisions of hadronic CRs with the interstellar gas. The gamma-ray intensity in the hadronic model is parametrized as $$E_{\gamma }^2\frac{\text{d}{F}_{\gamma }}{\text{d}{E}_{\gamma }}=\frac{E_{\gamma }}{4\pi }{\displaystyle \int n_{\text{H}}}{\sigma }_{\text{p}}{v}_{\text{p}}{\left(\frac{\text{d}\Sigma }{\text{d}T}\right)}_{\text{}\text{p}}\text{d}{T}_{\text{p}},$$(9)

where ${\left(\frac{\text{d}\Sigma }{\text{d}T}\right)}_{\text{}\text{p}}={\displaystyle \int {\left(\frac{\text{d}n}{\text{d}T}\right)}_{\text{}\text{p}}}\text{d}R$ is the column density of CR protons, v_p is the velocity of the protons, $T_{\text{p}}=\sqrt{{(qc)}^2+{(m_{\text{p}}{c}^2)}^2}-m_{\text{p}}{c}^2$ is the kinetic energy of the protons as a function of momentum q, σ_p (E_γ, T_p) is the differential cross section in units of $E_{\gamma }\frac{\text{d}\sigma }{\text{d}{E}_{\gamma }}$ for gamma rays in proton-proton collisions (Kamae et al. 2006; Karlsson & Kamae 2008), and n_H is the density of gas. We use n_H = 1 cm⁻³ as a reference density, which is comparable with the gas surface density of ~ 30 M_⊙ pc⁻² (Miville-Deschênes et al. 2017) averaged over the approximate height of the enhanced emission at the base of the bubbles, i.e., ~ 900 pc above and below the GC. Since most of the gas density is concentrated within 100–200 pc from the midplane (e.g., Sofue 2013; Miville-Deschênes et al. 2017), the gas density in the plane near the GC is 5–10 cm⁻³, which is consistent with the enhancement of the intensity of the emission at the base of the FBs at ℓ = 0°. We model the proton spectrum as a power law of the momentum $\frac{\text{d}\Sigma }{\text{d}qc}=n_{\text{p}}{q}^{-{\gamma }_{\text{p}}}$; we note that $v_{\text{p}}\frac{\text{d}\Sigma }{\text{d}{T}_{\text{p}}}=c\frac{\text{d}\Sigma }{\text{d}qc}$.

We fit the hadronic model of the gamma-ray emission at the base of the FBs analogously to the leptonic model in Sect. 4.4. The dash-dotted line in Fig. 11 represents the best-fit hadronic spectrum (labeled as π⁰). The index of the proton spectrum is relatively hard γ_p ≲ 2.3, especially to the west of the GC.

We calculate the significance of a cutoff in the CR proton spectrum by adding an exponential cutoff factor and refitting the model to the gamma-ray data. The improvement in the model and the 95% confidence lower bound on the cutoff values are presented in Table 3. Within ± 6° from the GP west of the GC, the 95% confidence level for the lower bound on the cutoff among the models of the foreground emission that we have considered in Sects. 3.2–3.4 is about 6 TeV.

4.6 Summary of the spectral analysis

In Fig. 12 we show the envelopes of the gamma-ray spectra at |b| < 2° of the emission at the base of the FBs for the models of foreground emission that we have considered, including the models in which we change the selection of the low energies in the definition of the foreground emission model (Appendix A). In order to determine the maximal and minimal models of the FBs in the GP, we fit the maximal and minimal points in the envelope above 3 GeV with a power law with a cutoff function; we fit above 3 GeV because some of the models of the foreground emission in Appendix A are determined for energies between 1 and 2.2 GeV. The corresponding parameters are reported in the first row of Table 4. We also fit the IC and hadronic models to the maximal and minimal points in the envelope and report the corresponding parameters in Table 4.

5 Intepretation

The bright and hard gamma-ray emission at the base of the bubbles can be at any position along the line of sight. In this section we consider two characteristic scenarios: that the emission is near the GC or that the emission is produced by one or a few supernova remnants (SNRs) closer to us.

5.1 Emission near the GC scenario

In the estimates in this subsection we use the following characteristic sizes: the distance to the GC is 8.5 kpc, the size of the region with the enhanced gamma-ray emission interpreted as an additional population of CRs is − 10° < ℓ < 0° and |b| < 6°. If we assume that the geometry of the region is a simple box, then the size of the box along the GP is 1.5 kpc, while the half-size in the vertical direction (distance from the GP to the boundary) is Δh = 0.9 kpc. The total gamma-ray luminosity in this case can be estimated by summing the data points in the rectangles model in Fig. 11 for |b| < 6°, which yields the luminosity of L = 1.1 × 10³⁷ erg s⁻¹ above 1 GeV.

One of the most intriguing features of the gamma-ray spectrum at the base of the FBs is the absence of a significant cutoff up to 1 TeV and an inferred hard spectrum of underlying electrons or protons. In particular, if the emission originates in hadronic interactions, the spectrum of the CR protons (CRp) is ~ E^−2.2, which is significantly harder than the propagated spectrum ~E^−2.5−2.7 observed both locally (e.g., Casandjian 2015; Remy et al. 2017; Neronov et al. 2017; Prokhorov & Colafrancesco 2018) and throughout the GP (e.g., Gaggero et al. 2015; Acero et al. 2016; Yang et al. 2016; Aharonian et al. 2018). Although there is a hardening of the spectrum of CRs in the inner rings of the Galaxy and near the GC, the inferred spectra of the CRs at the base of the FBs are significantly harder than the propagated spectra of the CRs in the inner rings and near the GC (Gaggero et al. 2015; Acero et al. 2016; Yang et al. 2016; Aharonian et al. 2018). We assume that the CRp spectrum at the base of the bubbles is equal to the injection spectrum unaffected by the propagation softening. This is possible if the CRp were injected relatively recently and had insufficient time to escape from the region of the enhanced emission, i.e., the age of the CRp is less than the propagation time to cross the vertical distance of 0.9 kpc. If we assume a spatially constant diffusion coefficient $D(E)=D_0{\left(\frac{E}{1\text{\hspace{0.17em}GeV}}\right)}^{\delta }$ with D₀ = 3e28 cm² s⁻¹ = 100 pc² kyr⁻¹ and δ = 0.4 (Strong et al. 2007), then the escape time for the protons at E > 6 TeV is $$T_{\text{esc}}<\frac{\Delta h^2}{2D(E)}\approx 100\text{\hspace{0.17em}kyr}.$$(10)

This gives us an approximate upper bound on the age of the proton CRs, assuming that the diffusion coefficient near the GC is similar to the diffusion coefficient near Sun. The escape time can be significantly longer if the diffusion length near the GC is smaller than the local diffusion length or if the CRs are confined by a particular configuration of magnetic fields. The CRp energy density within |b| < 6° in the rectangles model of the FBs (Sect. 4.5) normalized to the reference density of n_H = 1 cm⁻³ is d E_tot∕dV = 360 meV cm⁻³ above 1 GeV. In Fig. 13 (left), we compare the corresponding flux to the local CRp flux. The total energy in CRp within |b| < 6° is E_tot = 7e52 erg. This population of CRp can be obtained from ~700 SNe, assuming that on average SNe inject ~10⁵⁰ erg in CRp, which corresponds to about 10% efficiency of CR acceleration by a SN with kinetic energy ~ 10⁵¹ erg (e.g., Spurio 2015). There are several populations of young stellar objects near the GC (Yusef-Zadeh et al. 2009; Immer et al. 2012); in particular there is a population of stars close to the GC with an age of about 6 Myr (Paumard et al. 2006). If we assume thatthe rate of SNe in the Galaxy is about one per century and that the GC region contributes at the level of a few percent to the total rate, for example a few SNRs per 10 kyr, then the 700 SNRs can be produced near the GC in the past few million years.

The best-fit spectrum of CRe is ~E^−2.7, which is softer than the spectrum of the protons. Nonetheless, this spectrum is harder than the spectrum expected in the presence of cooling for a stationary population of CRe. Thus, unless the injection spectrum is harder than E⁻², the population of the CRe at the base of the bubbles is not affected by cooling up to 3 TeV, which is the 95% confidence minimal lower bound on the cutoff in the CRe (Table 2). The cooling time for the electrons at 3 TeV is T_cool ≈ 200 kyr, which is comparable to the diffusive escape time from the ± 900 pc volume around the GP at 3 TeV. The cooling time can have a factor ~3 uncertainty due to the uncertainties in the ISRF near the GC (see Appendix C).

The energy density in electrons with energies above E₀ = 1 GeV, which is necessary to produce the gamma-ray emission in the rectangles model of the FBs (Sect. 4.4), is 4.0 meV cm^-3. In Fig. 13 (right), we compare the corresponding flux to the local CRe flux. The total energy content of the ROI in electrons above 1 GeV is E_tot = 3e51 erg, which corresponds to the CR energy output of 3000 SNe, assuming a 0.1% efficiency in converting the SNR kinetic energy to CRe energy. The ratio of the electron and proton acceleration efficiencies K_e/p = 10⁻² assumed in this work is relatively optimistic, given that the ratio of efficiencies can be as low as ~ 10⁻³ (e.g., Park et al. 2015). Similar efficiencies of K_e/p ∝ 10⁻² can be expected in the acceleration of electrons and protons in jets from SMBHs (e.g., Ball et al. 2018). Since the escape time is comparable to the electron cooling time and the required number of SNe in the hadronic scenario is smaller than the number of SNe in the leptonic scenario (for the ratio of acceleration efficiencies K_e/p ≲ 10⁻²), we conclude that the majority of the gamma-ray emission near the base of the FBs is likely to be produced by the hadronic interactions of CRp with the gas. The presence of the hadronic production of gamma rays at the base of the FBs, without any sign of a cutoff up to gamma-ray energies of 1 TeV, is important for the detectability of the associated neutrino signal by neutrino telescopes; see Fig. 14 (see also Razzaque & Yang 2018).

Fig. 11 Rectangles model SEDs of the FBs in the latitude stripes (2°, 6°), (−2 °, 2°) and (−6°, −2°) for negative (blue) and positive (red) longitudes. We compare the parametric model of the SED with the leptonic (Sect. 4.4) and hadronic (Sect. 4.5) models.

Fig. 12 Spectral energy distribution of the FBs in the GP. The shaded areas show the envelope of the spectra of the FBs in all foreground models considered in the paper including the changes in the choice of the low-energy range to model the foreground emission (discussed in Appendix A). The lines show the fits to the maximal and minimal points in the envelopes above 3 GeV.

Fig. 13 Proton (left) and electron (right) spectra fitted to the gamma-ray emission at the base of the FBs. Solid lines show the best-fit spectra in the rectangles model of the FBs. The shaded bands represent the systematic uncertainties, estimated from the maximal and minimal best-fit FBs spectra among the foreground models. Uncertainties in the ISRF energy density near the GC result in an additional factor of ~2 uncertainties in the CRe energy density (see Appendix C for details). For comparison with the local CR spectrum, we show the local fluxes of CRp measured by AMS-02 (Aguilar et al. 2015) and CRe measured by H.E.S.S. (Aharonian et al. 2008), AMS-02 (Aguilar et al. 2014), Fermi-LAT (Abdollahi et al. 2017), and DAMPE (Ambrosi et al. 2017). The error bars represent the systematic and statistical uncertainties added in quadrature.

Fig. 14 Comparison of the hadronic model of the gamma-ray emission at the base of the FBs with CTA and KM3NeT for a source with 2° radius (Ambrogi et al. 2018). The minimal (min) and maximal (max) models correspond to the minimal and maximal models in the hadronic scenario in Table 4 (for the maximal model we assume a cutoff in the proton spectrum at 1 PeV) integrated over a 2° circle to transform intensity to flux. The KM3NeT sensitivity is for muon neutrinos only. In estimating the sensitivity, we take intoaccount that the GC is below the horizon for about two-thirds of the time at the KM3NeT location in the Mediterranean Sea.

5.2 Nearby SNR or a superbubble scenario

In this subsection we discuss models where the hard and bright gamma-ray emission at the base of the FBs is created either by a single SNR or several SNRs at a distance closer than the distance to the GC. We start with a single SNR scenario. For the leptonic model we use the localISRF and calculate the IC emission as described in Sect. 4.4. Assuming one SNR, we find a distance of 50 pc. The radius of an SNR corresponding to 6° circle at a distance of 50 pc is about 5 pc. Provided that the SNR would be rather close to Earth, it should be detected in radio and X-ray observations, however no suitable candidates are currently known (Green 2014, 2017). For the hadronic model of the gamma-ray emission, we assume the local gas density of 1 cm⁻³. At the GC, we need about 700 SNRs for this reference gas density; thus the distance at which a single SNR can explain the excess emission is $d=\frac{8.5\text{kpc}}{\sqrt{700}}\approx 300\text{pc}$. We can compare the flux from the low-latitude bubbles with the emission of a known SNR, such as Tycho’s SNR. The differential energy flux of Tycho’s SNR at 10 GeV is E² d F∕dE ≈ 10⁻⁶ MeV cm⁻² s⁻¹ (Archambault et al. 2017). The solid angle of the ROI of 10° × 12° at the base of the FBs is Ω ≈ 0.037 sr. If we assume that Tycho’s SNR is at a distance of 3 kpc, then the average intensity inside Ω at a distance of 300 pc is E² d N∕dE ≈ 2.7 × 10⁻⁶ GeV cm⁻² s⁻¹ sr⁻¹, which is comparable to the FBs emission at latitudes |b| < 6° (e.g., Fig. 12). The radius of Tycho’s SNR is about 3.5 pc (for the angular diameter of 8′ at 3 kpc), while the radius corresponding to 6° angular sizeat 300 pc is about 30 pc, which is almost ten times larger than the radius of Tycho’s SNR. If we assume that only the core of the emission at the base of the FBs within |b| < 2° corresponds to an SNR, the corresponding linear size is ≈10 pc (for the same distance of 300 pc). There are several SNRs in Green’s catalog (Green 2014, 2017) between − 10° < ℓ < 0°, but all of these have relatively small angular diameters ≲30′. Thus a single SNR is not a plausible scenario for the emission at the base of the FBs, unless this SNR is missing in Green’s catalog.

Another possibility is that the emission at the base of the FBs is created by several SNRs and/or wind from massive stars, similar to the Cygnus cocoon (Ackermann et al. 2011) or the Loop I and the local superbubbles (Smith & Cox 2001; Breitschwerdt & de Avillez 2006; Wolleben 2007; Breitschwerdt et al. 2016; Schulreich et al. 2018; Dickinson 2018; Shchekinov 2018). The Cygnus cocoon (Ackermann et al. 2011) has a radius of about 2°, which at a distance of 1.4 kpc corresponds to a radius of about 50 pc. If a Cygnus-like region is responsible for the emission within |b| < 6° at the base of the FBs, then such region should be at a distance of about 500 pc. If only the core of the high-intensity region within |b| < 2° is explained by such a cocoon, then the distance would be similar to the distance to the Cygnus cocoon, i.e., about 1.4 kpc. In Fig. 15 we compare the gamma-ray intensity at the base of the FB integrated over a 2° circle to the flux from the Cygnus cocoon (Abdo et al. 2007; Ackermann et al. 2011). The flux in the maximal model is comparable (within a factor of 2) to the flux from the Cygnus cocoon between 1 and 100 GeV.

Massivestars, which can inflate cocoons or superbubbles by stellar winds and SN explosions, are usually found in young open stellar clusters. There are several young open stellar clusters in the direction of the base of the FBs, i.e., within − 10° < ℓ < 0°. For example, clusters Trumpler 27 (Moffat et al. 1977) or NGC 6383 (Lloyd Evans 1978) contain more than ten massive OB stars, have an age of about 10 Myr, and an estimated distance of 1.2 and 1 kpc, respectively, in the direction of ℓ ≈ 356°, b ≈ 0°. There is also a cluster NGC 6405 (Rohlfs et al. 1959), which has a smaller estimated distance of about 500 pc but an older age ~ 100 Myr. The turbulence created by winds from the massive stars and the past SNRs in these clusters can (re)accelerate CRs and it would also lead to smaller diffusion length, which would prevent the escape of CRs on timescales shorter than a few Myr (Ackermann et al. 2011). There are relatively few SNRs detected in gamma rays in that area of the sky. In particular, there is only one SNR, 3FGL J1741.1−3053, in the 3FGL catalog (Acero et al. 2015)that is between − 10° < ℓ < 0°. This SNR is associated with the Tornado nebula at a distance of 12 kpc (Castro et al. 2013). Green’s catalog of SNRs contains several SNRs between − 10° < ℓ < 0°, but SNRs with measured distances are more than 4 kpc away. The absence of detected SNRs in the direction of the base of the FBs at distances ≲ 1.5 kpc does not necessarily mean that there were no SNRs in the open clusters mentioned above. A possible reason is that the wind from the massive stars as well as the first SNRs in the cluster create a hot superbubble with smaller gas density than the average density in the GP. As a result, the subsequent SNRs quickly expand in the less dense environment until they reach the common envelope of the superbubble. For instance, this is observed in simulations of the formation of the local hot superbubble (Breitschwerdt et al. 2016; Schulreich et al. 2017, 2018). Whether the stellar clusters can provide the necessary CR injection rate to power the gamma-ray emission at the base of the FBs deserves a separate detailed study.

An interesting question is whether the FBs themselves can be inflated by a Cygnus-like superbubble, which accidentally happened along the line of sight towards the GC. If the FBs are at a distance of about 1 kpc, then their vertical size is also about 1 kpc, which is ten times smaller than the size of the FBs at a distance of 8.5 kpc. In the leptonic scenario of the FBs, the 10⁵² erg in CRe at 8.5 kpc (Ackermann et al. 2014) would correspond to 10⁵⁰ erg at 1 kpc. This power can be provided by several hundred SNe, which is a relatively large number given the few young stellar clusters at distances ≲1 kpc in the direction of the base of the FBs. This estimate is based on an assumed 0.1% efficiency in acceleration of CRe by SN shells; if there were an additional acceleration of the electrons by the turbulence either at the base of the FBs or inside the high-latitude FB volume, then the required number of SNRs would be smaller.

The electrons can be delivered to the FBs volume by an outflow that is possibly created by the pressure of the CRs themselves (e.g., Jacob et al. 2018). The measured red- and blueshift velocities of the gas outflow in the direction of the FBs are about 200–300 km s⁻¹, which implies the gas outflow from the GP in a biconical model of the FBs with velocities ≳ 900 km s⁻¹ and an age of the outflow at the position above the GC about 6–9 Myr (Fox et al. 2015; Bordoloi et al. 2017). If the distance is ten times smaller, then the age would be ~ 600 − 900 kyr, which is comparable to the cooling time of 1 TeV electrons necessary to explain the gamma-ray emission from the FBs at high latitudes (Ackermann et al. 2014). Thus the gamma-ray emission at high latitudes could be explained by an IC model with about 100 SNRs (or fewer if there is re-acceleration of electrons) created ~1 Myr ago, while the bright and hard emission at the base of the bubbles can be explained by the hadronic emission from about 10 SNRs. The actual age of the CRs is determined by the confinement within the superbubble and can be also on the order of 1 Myr or more. In this scenario, the star-forming region, which can power the emission at the base of the FBs and, possibly, the FBs themselves can be situated in the Sagittarius arm at a distance of 1–1.5 kpc in the direction of the GC.

Fig. 15 Comparison of the minimal and maximal models of the residual flux at the base of the FBs with the Cygnus cocoon spectrum. The minimal and the maximal models are the same as in Fig. 14. For the 3FGL spectrum we take the power-law parameters from the catalog (Acero et al. 2015) and add a cutoff in the gamma-ray spectrum at 10 TeV so that the SED is consistent with the Milagro source MGRO J2031+41 at 20 TeV (Abdo et al. 2007). The purple squares represent the measurement of the flux from the Cygnus cocoon by the Fermi-LAT collaboration (Ackermann et al. 2011).

6 Conclusions

In this paper we use nine years of Fermi-LAT Pass 8 Source class data to study the gamma-ray emission at the base of the FBs. We use multiple methods to construct the foreground diffuse gamma-ray components and to determine the properties of the residual gamma-ray emission at the base of the FBs. We confirm the earlier findings that the emission at the base of the FB has a higher intensity than the FBs emission at high latitudes. The spectrum at the base of the FBs is consistent with a single power law without a cutoff up to about 1 TeV. The emission is slightly displaced from the GC to negative longitudes, which favors a starburst scenario of the bubbles formation, unless there is a mechanism to shift the gamma-ray emission away from Sgr A* in the SMBH scenario.

The gamma-ray emission at the base of the FBs can be explained by either a hadronic or leptonic model of gamma-ray production with CR spectra consistent with a power law without a cutoff. The index of the CRe (CRp) spectrum is 2.6–2.9 (2.0–2.3). We derive the 95% confidence lower bound on the cutoff in the CR electron and proton spectra by selecting the lowest 95% statistical confidence value among the models of the foreground emission considered in Sects. 3.2–3.4. Within |b| < 6° and − 10° < ℓ < 0°, the 95% confidence lower bound on the exponential cutoff in the CRe (CRp) spectrum is about 3 (6) TeV.

If the location of the residual emission is near the GC, then the total gamma-ray luminosity is L ≈10³⁷ erg s⁻¹ and the required energy of CRe (CRp) above 1 GeV is E_e = 3e51 erg (E_p= 7e52 erg). Although the total required energy in CRe is smaller than in CRp, the efficiency of acceleration of hadronic CRs is expected to be significantly higher than for leptonic CRs. In particular, with 10% CRp acceleration efficiency, about 700 SNRs are needed to explain the gamma-ray emission in the hadronic scenario near the GC, while with 0.1% CRe acceleration efficiency, the required number of SNRs is 3000. For the GC scenario, the estimated diffusive escape time and the electron cooling time is ~ 100 kyr. In this case a relatively high SN rate of one per ~100 yr at the location of the gamma-ray emission is required. If the average ISRF energy density is a factor of ~2 lower or higher than the ISRF in the reference model (Appendix C), then both the cooling time and the required energy density of CRe is a factor of ~2 higher or lower than in the reference model, which nevertheless requires the same SN rate of one per ~100 yr.

If the residual gamma-ray emission is produced at a closer distance, for instance at about 1 kpc by a Cygnus-like cocoon or a superbubble in the Sagittarius arm, then the required number of SNRs is about two orders of magnitude smaller. For the superbubble scenario, the characteristic linear size is 10 times smaller, which leads to 100 times shorter diffusion time of a few kyr. This is shorter than the lifetime of individual SNRs. Both in the GC and in the superbubble scenarios, the age of the CR population can be greater than the diffusion time if there is an additional confinement by the SN shells or the shell of a superbubble.

Apart from the question of the location of the emission at the base of the FBs, we can revisit the question of the location of the FBs themselves: the FBs could have been inflated by a superbubble about 1 Myr ago at a distance of 1 kpc from the Sun. The CRs in the past superbubble may create a wind with velocities up to ~1000 km s⁻¹, which is sufficient to inflate a 1 kpc large bubble within 1 Myr so that the CR electrons in the wind do not cool below 1 TeV and can explain the gamma-ray emission from the FBs at high latitudes. The gamma-ray emission at low latitudes is then explained (most naturally) by a hadronic scenario from a more recent starburst episode at about the same location.

In the future, the bright and hard gamma-ray emission at the base of the FB should be detectable with CTA, or a version of the High-Altitude Water Cherenkov (HAWC) Gamma-Ray Observatory in the southern hemisphere (Mostafa et al. 2017; Assis et al. 2017; Razzaque & Yang 2018). Observations with Cherenkov telescopes should detect or constrain the cutoff energy in the gamma-ray spectrum in multi-TeV regime. In the hadronic scenario, for a sufficiently high cutoff value of the CRp spectrum, for instance around 1 PeV, the corresponding neutrino emission could be detected with the future KM3NeT telescope.

Acknowledgements

The authors would like to thank Seth Digel, Anna Franckowiak, Stefan Funk, Gudlaugur Johannesson, Tsunefumi Mizuno, Troy Porter, Andrew Strong, and Luigi Tibaldo for valuable comments and suggestions. Also we would like to thank Cristina Popescu for providing us the ISRF data files in fits format. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Appendix A Modeling and systematic uncertainties

In this appendix we study the dependence of the spectrum of gamma-ray emission at the base of the FB on the selection of the low-energy range in the definition of the diffuse foreground model and on the choice of the class of the events (UltraCleanVeto versus Source). We use the rectangles model of the FB in Sect. 3.3 as our baseline. In this model, the foreground emission template consists of Source data integrated over the energy interval 0.3−1.0 GeV. To probe the dependence on the choice of the low-energy range in the model, we pick three non-overlapping energy ranges, 0.3−0.5, 0.5−1.0, and 1.0−2.2 GeV, and use the same analysis as for the baseline model. We also compare with the spectra obtained using UltraCleanVeto data both in the definition of the model at low energies and in the analysis at high energies. The residual spectra resulting from the different energy ranges and data classes are shown in Fig. A.1. The agreement among the models is reasonably good at E > 10 GeV, especially for negative longitudes.

Fig. A.1 Residual SED in the rectangles model using Source class (dotted lines with markers) and UltraCleanVeto class data (solid lines) in four different energy ranges used to determine the foreground diffuse emission template. The baseline model, presented in Sect. 3.3, has the low-energy range in the definition of the foreground model 0.3− 1.0 GeV (black squares).

As shown in Fig. A.1, the highest energy spectral point for negative longitudes may have an upward statistical fluctuation, which can influence the conclusions about the lower limit in the cutoff values. In order to test the dependence on the highest energy spectral point, we exclude the energy bin 680 GeV−1 TeV and repeat the derivation of the best-fit parameters in the rectangles model of the FBs. The results are presented in Table A.1. Without the last data point, the four areas with infinite cutoff energy now have a finite cutoff, but with a significance 2Δ log L < 1, i.e., the models are still consistent with a simple power law. The 95% confidence lower limit on the cutoff value is also smaller for the dataset without the last point, but in this case we are using the data up to 680 GeV rather than 1 TeV.

Appendix B Gamma distribution as a likelihood for smoothed data

In this appendix we show that smoothing the data and using gamma distribution instead of the Poisson distribution is a reasonable procedure, at least in the case when the scale of variations in the diffuse emission is larger than the smoothing radius. We split the log likelihood into a sum over areas with size approximately equal to the smoothing scale and neglect the correlation between the different areas after smoothing the data. We suppose that one of these areas has n pixels with random counts k_i drawn from a Poisson distribution with mean λ₀; in this case we assume that the underlying diffuse flux is approximately constant within the smoothing radius. The likelihood function for parameter λ is $$L(\lambda )={\displaystyle \prod _i\frac{{\lambda }^{k_i}}{k_i!}}{e}^{-\lambda }.$$(B.1)

The log likelihood is $$\text{log}L={\displaystyle \sum _{i=1}^n}(-\lambda +k_i\text{log}\lambda )+\text{const}.$$(B.2)

The maximum likelihood value λ_* is determined from $$0=\frac{\partial \text{log}L}{\partial \lambda }=-n+\frac{1}{\lambda }{\displaystyle \sum _{i=1}^n}{k}_i,$$(B.3)

which gives ${\lambda }_{*}=\frac{1}{n}{\sum }_{i=1}^n{k}_i$. The uncertainty is $$\frac{1}{{\sigma }^2}=-{\frac{{\partial }^2\text{log}L}{\partial {\lambda }^2}|}_{\lambda ={\lambda }_{*}}=\frac{n}{{\lambda }_{*}},$$(B.4)

which gives σ² = λ_*∕n.

If the smoothing radius is comparable to the size of the region, then we can approximate the values k_i with the average in the region ${\tilde{k}}_i=\overline{k}=\frac{1}{n}{\sum }_{i=1}^n{k}_i$; we note that ${\tilde{k}}_i$ are not integers. We take the gamma distribution for the likelihood function $$\tilde{L}={\displaystyle \prod _i\frac{{\lambda }^{{\tilde{k}}_i}}{\Gamma ({\tilde{k}}_i+1)}}{e}^{-\lambda }.$$(B.5)

The log likelihood is $$\text{log}\tilde{L}={\displaystyle \sum _{i=1}^n}(-\lambda +{\tilde{k}}_i\text{log}\lambda )+\text{const}.$$(B.6)

The maximum likelihood solution is the same as in the Poisson case${\tilde{\lambda }}_{*}=\overline{k}={\lambda }_{*}$. The uncertainty is also the same as in the Poisson case, i.e., $$\frac{1}{{\tilde{\sigma }}^2}=-{\frac{{\partial }^2\text{log}\tilde{L}}{\partial {\lambda }^2}|}_{\lambda ={\tilde{\lambda }}_{*}}=\frac{n}{{\tilde{\lambda }}_{*}}=\frac{1}{{\sigma }^2}.$$(B.7)

Thus, smoothing the data and using the gamma distribution (in this case) gives the same result as using the original Poisson distribution, which shows that the procedure is well defined from the statistical point of view and it gives reasonable results.

Appendix C Modeling uncertainty of ISRF near the GC

In this appendix we discuss the uncertainty of the IC model of the gamma-ray emission at the base of the FBs related to modeling of the ISRF near the GC. In order to estimate this uncertainty, we compare the ISRF model of Porter et al. (2008; available with the distribution of GALPROP v54.1), the ISRF model of Popescu et al. (2017), and the following two ISRF models of Porter et al. (2017): R12, based on Robitaille et al. (2012), and F98, based on Freudenreich (1998). In Fig. C.1 we show the corresponding densities of ISRFs averaged over the cylinder with the radius R = 1.5 kpc and the height z = ±0.3 kpc around the GC, which approximately corresponds to the latitude stripe |b| < 2° and |ℓ| < 10°. We notice that there can be up to a factor of 2 differences in the ISRF energy density at the peak of the SL emission (around 1 μm) and an order of magnitude differences in the IR wavelengths (around 100 μm), which can affect the inferred populations of CR electrons producing the IC gamma rays (see also Porter et al. 2017; Niederwanger et al. 2019).

Fig. C.1 Comparison of ISRF SEDs at the GC (see text for more details).

In Fig. C.2 we show the IC models of the gamma-ray emission at the base of the FBs in the rectangle b ∈ (−2°, 2°), ℓ ∈ (−10°, 0°) for the different ISRF models near the GC. We separate the SL, IR, and CMB contributions. The SL andIR ISRFs are separated by splitting the radiation field energy densities at 0.1 eV (≈ 12 μm). We model the CR electrons spectra by a power-law function with an exponential cutoff. We fix the cutoff at 1 PeV and fit the normalization and the spectral index by fitting the IC model to the Fermi-LAT spectral points. The corresponding parameters are presented inTable C.1. We use the Porter et al. (2008) ISRF as the reference model and show the normalizations of the CRe spectra relative to the reference model (the normalizations are determined at 1 GeV). The overall spread in the normalizations is about a factor of 3, while the overall variation of the index of the CRe spectra is less than about 0.2.

Fig. C.2 Contribution of various components of ISRFs to the IC model of the gamma-ray emission. The data points correspond to the emission at the base of the FBs in the rectangle b ∈ (−2°, 2°), ℓ ∈ (−10°, 0°) (middle panel in Fig. 11). The spectrum of CR electrons is modeled as a power-law function with an exponential cutoff at 1 PeV. We separate the SL and IR contributions to the ISRF by formally splitting the ISRF at 0.1 eV.

References

All Tables

All Figures

Fig. 1 Energy flux of the Fermi-LAT 9 yr of Source class data integrated in energy ranges 1–10 GeV, 10–100 GeV, and 100 GeV–1 TeV. The maps are represented in logarithmic scale. The map in the energy range 10–100 GeV (100–1000 GeV) is smoothed with a Gaussian kernel with radius σ =0. °5 (σ =0. °7). The plots are created with the HEALPix package (Górski et al. 2005) in Mollweide projection with nside = 128 (≈ 0. °46 pixel size).

In the text

	Fig. 2 Difference west minus east in the Fermi-LAT intensity relative to the GC after masking bright PSs. ThePS mask is symmetrized relative to the GC. The west (east) region is defined between − 10° < ℓ < 0° (0° < ℓ <10°). Gray lines show spectra for latitudes |b| > 10°.
In the text

	Fig. 3 Energy flux of the diffuse emission model based on the low-energy data (Sect. 3.2) integrated over three energyranges.
In the text

	Fig. 4 Energy flux of the residuals of the low-energy model derived in Sect. 3.2 integrated over three energy ranges.
In the text

	Fig. 5 Energy flux of the diffuse model excluding the rectangles model of the FBs (top left), the rectangles model of the FBs (top right), the residual after subtracting the total model (bottom left), and residuals plus the rectangles model of the FBs (bottom right) derived in Sect. 3.3 integrated between 10 and 100 GeV.
In the text

	Fig. 6 Energy flux of the background model based on rescaled GALPROP templates (left) and of the residuals plus the FBs and thegNFW component (right) derived in Sect. 3.4 integrated between 10 and 100 GeV.
In the text

	Fig. 7 Latitude profiles of the energy flux between 10 and 100 GeV for the total data excluding the PS mask and for the FBs models for different foreground diffuse emission models to the east (left panel) and to the west (right panel) of the GC.
In the text

	Fig. 8 Longitude profiles of the energy flux between 10 and 100 GeV and between 100 GeV and 1 TeV for the residuals in the low-energy data background model (Sect. 3.2).
In the text

	Fig. 9 Gamma-ray data SEDs of the residuals (LowE model), rectangles model of the FBs (rectangles model), residuals plus FBs model plus gNFW DM annihilation template (GALPROP model), and of the difference in the data west minus east of the GC (excluding pixels in the PS mask).
In the text

	Fig. 10 Index of the log-parabola defined in Eq. (5) at E = 500 GeV as a function of latitude for the FBs in the different foreground diffuse emission models.
In the text

	Fig. 11 Rectangles model SEDs of the FBs in the latitude stripes (2°, 6°), (−2 °, 2°) and (−6°, −2°) for negative (blue) and positive (red) longitudes. We compare the parametric model of the SED with the leptonic (Sect. 4.4) and hadronic (Sect. 4.5) models.
In the text

	Fig. 12 Spectral energy distribution of the FBs in the GP. The shaded areas show the envelope of the spectra of the FBs in all foreground models considered in the paper including the changes in the choice of the low-energy range to model the foreground emission (discussed in Appendix A). The lines show the fits to the maximal and minimal points in the envelopes above 3 GeV.
In the text

Fig. 13 Proton (left) and electron (right) spectra fitted to the gamma-ray emission at the base of the FBs. Solid lines show the best-fit spectra in the rectangles model of the FBs. The shaded bands represent the systematic uncertainties, estimated from the maximal and minimal best-fit FBs spectra among the foreground models. Uncertainties in the ISRF energy density near the GC result in an additional factor of ~2 uncertainties in the CRe energy density (see Appendix C for details). For comparison with the local CR spectrum, we show the local fluxes of CRp measured by AMS-02 (Aguilar et al. 2015) and CRe measured by H.E.S.S. (Aharonian et al. 2008), AMS-02 (Aguilar et al. 2014), Fermi-LAT (Abdollahi et al. 2017), and DAMPE (Ambrosi et al. 2017). The error bars represent the systematic and statistical uncertainties added in quadrature.

In the text

Fig. 14 Comparison of the hadronic model of the gamma-ray emission at the base of the FBs with CTA and KM3NeT for a source with 2° radius (Ambrogi et al. 2018). The minimal (min) and maximal (max) models correspond to the minimal and maximal models in the hadronic scenario in Table 4 (for the maximal model we assume a cutoff in the proton spectrum at 1 PeV) integrated over a 2° circle to transform intensity to flux. The KM3NeT sensitivity is for muon neutrinos only. In estimating the sensitivity, we take intoaccount that the GC is below the horizon for about two-thirds of the time at the KM3NeT location in the Mediterranean Sea.

In the text

Fig. 15 Comparison of the minimal and maximal models of the residual flux at the base of the FBs with the Cygnus cocoon spectrum. The minimal and the maximal models are the same as in Fig. 14. For the 3FGL spectrum we take the power-law parameters from the catalog (Acero et al. 2015) and add a cutoff in the gamma-ray spectrum at 10 TeV so that the SED is consistent with the Milagro source MGRO J2031+41 at 20 TeV (Abdo et al. 2007). The purple squares represent the measurement of the flux from the Cygnus cocoon by the Fermi-LAT collaboration (Ackermann et al. 2011).

In the text

	Fig. A.1 Residual SED in the rectangles model using Source class (dotted lines with markers) and UltraCleanVeto class data (solid lines) in four different energy ranges used to determine the foreground diffuse emission template. The baseline model, presented in Sect. 3.3, has the low-energy range in the definition of the foreground model 0.3− 1.0 GeV (black squares).
In the text

	Fig. C.1 Comparison of ISRF SEDs at the GC (see text for more details).
In the text

Fig. C.2 Contribution of various components of ISRFs to the IC model of the gamma-ray emission. The data points correspond to the emission at the base of the FBs in the rectangle b ∈ (−2°, 2°), ℓ ∈ (−10°, 0°) (middle panel in Fig. 11). The spectrum of CR electrons is modeled as a power-law function with an exponential cutoff at 1 PeV. We separate the SL and IR contributions to the ISRF by formally splitting the ISRF at 0.1 eV.

In the text