© ESO, 2014

1. Introduction

Two huge, bubble-like structures have been reported by Su et al. (2010), Dobler et al. (2010), and Su & Finkbeiner (2012) in Fermi-LAT gamma-ray data to extend ~50° above and below the Galactic center. The gamma-ray emission from these structures, dubbed the “Fermi bubbles” (FBs), exhibits an E^-2 type spectrum which is significantly harder than the spectrum of the diffuse gamma-ray emission from the Galactic disk. Remarkably, structures coincident or similar to the FBs can be seen at other wavelengths, including the (total intensity) microwave haze found in WMAP (Finkbeiner 2004) and, most recently, in Planck data (Planck Collaboration 2013); the polarized microwave structures reported by Jones et al. (2012); the large-scale, biconical structures found (Bland-Hawthorn & Cohen 2003) in ROSAT X-ray data (Snowden et al. 1997); and the giant polarized radio lobes recently found at 2.3 GHz (Carretti et al. 2013).

Several models have been proposed to explain both the morphology and spectral properties of the FBs (Cheng et al. 2011; Mertsch & Sarkar 2011; Crocker & Aharonian 2011; Crocker 2012, 2013; Zubovas et al. 2011; Zubovas & Nayakshin 2012; Guo & Mathews 2012; Guo et al. 2012; Yang et al. 2012). These models predict different energy-dependent morphologies. For example, in the simplest IC model (e.g., discussed by Su et al. 2010) with a spatially independent electron spectrum described by a power law electron population with a low energy cutoff at 500 GeV, one expects a tendency for softening of the gamma-ray spectra at high latitudes. This is due to the gradual reduction of the IC components produced via upscattering of the optical and UV photons – the main target fields contributing to the production of gamma-rays above 10 GeV that also decline toward high latitudes. On the other hand, in a simple one-zone hadronic model (e.g., Crocker & Aharonian 2011) the protons’ steady state distribution should produce a similar spectrum of gamma-rays at all latitudes. Thus, studies of the energy-dependent morphology of the FBs may shed light on the nature of the gamma-ray emission mechanism(s). This is the basic motivation for the current study.

Here we present an analysis based on 60 months of Fermi-LAT data. We find a similar overall morphology and spectrum for the FBs to those obtained by Su et al. (2010). The FBs exhibit a fairly homogeneous surface brightness. Important is that this implies non-homogeneous volumetric emissivities that consider projection effects (Su et al. 2010). The FBs, however, are not completely uniform because they exhibit some “hot spots” (referred to by Su et al. 2010; and Su & Finkbeiner 2012 as the “jet”, “donut”, and “cocoon”).

In this paper the significantly larger photon statistics and the availability of the recently updated Fermi science software tools (Ackermann et al. 2012b) allow us to investigate the gross spectral morphology of the FBs. A somewhat similar study has recently been conducted by Hooper & Slatyer (2013). However, they concentrate spectral features in low latitudes while we focus on the high latitude, especially the top of the bubble. We proceed by dividing each bubble into several slices to investigate possible gamma-ray spectral change with latitude. Employing different background models, we find, robustly, a spectral hardening or, more specifically, a deficit of low energy flux toward the top of the SFB. We are also able to extend the spectrum of the FBs to lower energy than previously attempted because of the improvement in the analysis software and the instruments response functions.

The paper is structured as follows. In Sect. 2 we present the results of our data analysis, in Sect. 3 we discuss the fitting of leptonic and hadronic models to the gamma-ray data, and in Sect. 4 we set out our conclusions.

2. Data analysis

We use the publicly-available data obtained by the Fermi-LAT in the gamma-ray energy interval from 100 MeV to 300 GeV over the period of 4 Aug. 2008 and 17 July 2013 (MET 239557417–MET 395797613). To avoid contamination from charged particles we use the ULTRACLEAN data set in the analysis. We adopt the instrument response function version P7V6 (Ackermann et al. 2012b). All events are binned in the all-sky map in HEALPIX format with NSIDE = 256. The sources in the 2nd Fermi catalog (Nolan et al. 2012) are subtracted from the counts map using the flux in the catalog. To take into account the energy-dependent point spread function (PSF) of the Fermi-LAT, we use the functional form described in the LAT homepage¹ when masking sources in the catalog.

We use a likelihood method for determining the FBs’ apparent morphology and spectral features for different diffuse emission templates. The likelihood function has the form log  (L) = ∑ _kilog (μ_i) − μ_i − log (k_i !), where k is the number of photons in the i-th bin in the counts map and μ_i is the predicted number of photons within a particular linear combination of the templates. The sum is over all the spatial bins in the map. The likelihood function is determined in different energy bins to derive the energy spectrum. For fitting different energy bins, rather than smoothing the map to an universal full width at half maximum (FWHM) as used in Su et al. (2010), we only smooth the diffuse templates with the Fermi PSF and then fit the counts map as a linear combination of these smoothed diffuse templates. The normalization of each diffuse template is left free in the likelihood fitting. We do not assume a priori the existence of a pair of bubble-like structures. We only use the spatial templates for π⁰ decay and IC gamma-rays generated by GALPROP (Vladimirov et al. 2011) in addition to the isotropic template related to the extragalactic gamma-ray background and cosmic ray contamination. Finally, we take into account the large diffuse structure high above the plane to Galactic north which may be connected with the nearby ISM feature Loop 1 (Su et al. 2010). In generating the diffuse templates with GALPROP, we adopt the default Galdef webrun setting, the 2D plane diffusion model that was tuned to fit the ACE data.

After subtracting the best-fit linear combination of the diffuse templates, we find residual maps. By summing over all energy bins above 1 GeV we obtain the image in Fig. 1 in which two bubble-like structures are clearly seen. Next we generate spatial templates for the FBs from the residual map. In the second step of our analysis these bubble templates are included. We employ the likelihood method mentioned above once again and obtain the spectrum of all the diffuse templates.

Fig. 1 Residual map above 2 GeV. Background subtraction is described in the text. Two bubble-like structures can be seen. To render the picture clearer we mask the bright Galactic plane | b | < 5°.

Fig. 2 SEDs of different slices of the north and south FBs. The numbers 1 to 4 run from low to high latitudes. Left panel: north bubble (NFB). Right panel: south bubble (SFB).

To derive the spectrum in different parts of the FBs, we divide the SFB (−55°<b < −25°) into four slices and the north Fermi bubble (NFB; 25°<b< 50°) into three. To avoid contamination from the Galactic plane, in the fitting we mask out the inner ±25° region. The position of each analysis slice can be found in Fig. 3. The SED of each slice is shown in Fig. 2 where the numbers 1 to 4 run from lowest to highest latitude. From the SEDs it is evident that the spectrum of the highest southern slice is in deficit at low energies relative to the other slices. It is important to note that this same slice suffers little from geometric projection effects and is thus a true reflection of the spectrum at the top of the SFB. The spectrum at the top of the SFB is, therefore, significantly different from that in the interior. Figure 3 reveals the different morphologies of the FBs in different energy bins; at high energies the SFB is clearly more extended than at low energies. This is highlighted in Fig. 4. It should also be noted that in the top slice of the SFB (South 4) there are no known point sources from the 2nd Fermi catalog, making our conclusion about the spectral variation quite robust.

Fig. 3 Residual maps for different energy bins. The left and right panels correspond to energy intervals 1–2 GeV and 10–30 GeV, respectively. The inner disk of the Galaxy (| b | < 2°) is masked. The boxes laid over the residual maps show the position of each slice.

Fig. 4 Zoom in of Fig. 3 for the SFB region. Left panel 1–2 GeV and right panel 10–30 GeV.

Regarding the energy-dependent morphology of the FBs, we also flag the following point to be dealt with in further work: as evident from Figs. 3 and 4, the SFB is also relatively more extended to Galactic west at high energies than low energies implying a spectral hardening going from east to west. As far as we are aware, no model for the FBs currently accounts for this effect. It is interesting that, because of this extension to the west, at high energies, the FBs come to more closely resemble the polarized radio lobes detected at 2.3 GHz (Carretti et al. 2013)².

Accurate determination of the morphology and spectral features of the bubbles strongly depends on reliable modeling of the diffuse background. Unfortunately, because of uncertainties in the distributions of cosmic rays and interstellar gas, our knowledge about the diffuse gamma-ray background is imperfect. In such circumstances, one can only introduce different templates that instantiate different estimates of the background and try to derive the true bubble flux with them using a likelihood method. Unfortunately, different choices of background templates alter the final result significantly and this may lead to significant systematic errors. To study such possible errors, we investigate the 128 GALPROP models listed by Ackermann et al. (2012a). The aim of Ackermann et al. (2012a) was to study the origin and propagation of cosmic-rays and the distribution of the interstellar medium by simultaneously fitting diffuse gamma-ray emission and cosmic ray data. Although all their models underestimate the GeV emission at low latitudes (Ackermann et al. 2012a), this has little influence on our results given our masking of the inner ±25° when deriving spectra. The 128 models are different in cosmic ray (CR) source distribution, CR halo size, HI spin temperature and E(B − V) magnitude cut. Following perusal of the online material for Ackermann et al. (2012a), we found that 64 models with z = 8 kpc and z = 10 kpc do not fit the ⁹Be/¹⁰Be data well (the derived curves fall outside of the error bars of nearly all the data points), where z is the height of the CR halo. Therefore we abandon these, and investigate the other 64 models as diffuse emission templates.

With the same procedure as described above, we find the SEDs in different slices with all chosen templates; the results are summarized in Fig. 5. Because of the large systematic errors, we may not claim a deficit at low energies and high latitude in the NFB. However, in the SFB the low energy deficit of the top slice is significant for every template. It should be noted that, due to the slight differences in overall flux normalization for the different templates choices, differences of spectral shape are smeared out if we plot them all together as in Fig. 5 which, therefore, suggests a smaller difference between the top and bottom slices than exists in reality. In Fig. 6 we show 8 examples of SEDs for the slices in the SFB which give the extreme cases of the gamma-ray flux. From these individual examples the low energy deficit in the top slice is much more evident. Results for the northern slices are also shown in Fig. 7 which reveal the low energy part of the top slice has very large uncertainties. The large systematics in the NFB are likely due to the fact that it is partially coincident with the extended Loop 1 structure whose physical origin and exact morphology and spectral features are still uncertain.

3. Discussion

The gamma-ray residual maps obtained in the current study – based on six years’ Fermi-LAT data – confirm that the surface brightness of each bubble is homogeneous at a gross level. This implies that the volumetric emissivity is not homogeneous, otherwise we would detect a higher brightness in the center of each bubble due to projection (Su et al. 2010). One possible explanation of this may be higher turbulence in the bubbles’ edges generating more efficient particle acceleration or stronger confinement of high energy particles there with a resulting increase in the local gamma-ray emissivity; alternatively, in a non-saturation hadronic model, higher gas densities toward the edges may have the same effect (Crocker et al. 2013).

Below we assume the low energy deficit of the gamma-ray spectrum at high latitudes found above is a real effect and discuss its possible implications.

To simplify our modeling, we assume each bubble is a sphere with a radius 3.4 kpc. The 2D projection of each bubble can be approximated as a half sphere in the region b > 30° for the NFB and b < −30° for the SFB. The center of each bubble is located at (0, 0, ±5 kpc) in Galactic coordinates. Assuming an intrinsic north-south symmetry, we only analyze the SFB and, to reveal most starkly the spectral change with latitude, we consider only slice South 1 (− 33°<b < − 25°) and South 4 (−55°<b < − 47°) in detail. Slice South 1 covers the center of the SFB. We use the ISRF value in Galprop at (0, 0, 5 kpc) for the average value of the slice South 1 in the calculation, w_IR = 0.18 eV/cm³ and w_opt = 0.9 eV/cm³. The height of South 4 is about 7 kpc, thus we adopt the ISRF energy density from Galprop at (0, 0, 7 kpc), 0.7 eV/cm³ for the optical component and 0.15 eV/cm³ for the IR component. The IR and optical photon fields are modeled as diluted blackbody (gray body) spectra with temperatures of 100 K and 5000 K, respectively.

Fig. 5 SEDs for the seven slices with all 64 templates we use in our analysis. The shaded areas span all derived SEDs. The left panel is for the NFB while the right is for the SFB.

In the hadronic scenario, motivated by the relative low energy deficit of the gamma-ray flux in the top slice, we fit the SED in South 1 with a pure power law proton spectrum while for South 4 we use a power law with a low energy break that we find needs to be at 20 GeV to fit the data. We assume the proton flux below the break is zero for simplicity (see Fig. 8). Our results are shown in Fig. 9. The required low energy break may arise naturally due to energy-dependent diffusion effects in a non-steady-state scenario: given South 4 is far from the (assumed) injection source in the plane, it may be that only high energy particles (which diffuse faster) have had time to travel there since a previous injection event or given the age of the structure. The formalism describing these energy-dependent diffusion effects can be found in Aharonian & Atoyan (1996) where both impulsive and continuous injection cases are described. The position of the low energy break can be estimated by equating the age of the bubble (or the time since a previous burst) with the diffusion time, i.e., the break energy is ${\mathit{E}}_{\mathrm{bk}}\mathrm{\simeq }\mathrm{(}{\mathit{d}}^{\mathrm{2}}\mathit{/}\mathrm{(}{\mathit{D}}_{\mathrm{0}}\mathit{t}\mathrm{)}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathit{\delta }}}$, where d is the distance of the top slice from the injection source, t is the age of the bubble (or burst) and D(E) = D₀E^δ is the energy-dependent diffusion coefficient. If we assume protons are injected at the Galactic center, d is about 7 kpc. Notice that in the kinetic equations which govern proton diffusion only the combination D₀t appears. Thus any timescale can be obtained if we tune the diffusion coefficient D₀. If D(E) takes a value similar to that in the Galactic plane, say 4 × 10²⁸ cm² s^-1 at 1 GeV and δ = 0.3, the age of the FBs should be t ~ 10⁸ yr (cf. Crocker et al. 2013). Alternatively, if there is continuous proton injection into the FBs over 10¹⁰ yr we need a rather small diffusion coefficient, D₀ ~ 10²⁶ cm² s^-1 at 1 GeV, which is about 1/100 of that in the Galactic plane.

Fig. 6 Six examples of SEDs for slices in the SFB.

Fig. 7 Six examples of SEDs for slices in the NFB.

If we assume the injection is impulsive and the ISM density is 0.005 cm^-3 the total energy in protons needed to light the bubbles is of the order 10⁵⁶ erg. Note that if we assume the simple geometry mentioned above and a magnetic field B ~ 10 μG as suggested by Carretti et al. (2013), then the energy densities of cosmic rays and magnetic fields are close to equipartition. A burst-like injection event requires that the duration of the injection is much smaller than the age of the structure which is t ~ 10⁸ yr (for D₀ = 4 × 10²⁸ cm² s^-1 at 1 GeV), i.e., the duration of the injection event would have to be of order 10⁷ yr or less in this case. Then the injection rate is of the order 10⁴¹ erg/s, two orders of magnitude larger than the X-ray luminosity of a X-ray reflection nebulae near the Galactic center (e.g., Gando Ryu et al. 2012), but still only one thousandth of the Eddington Luminosity of Sgr A^∗, ~5 × 10⁴⁴ erg/s.

It should be noted that we neglect the hadronic secondary electrons (and positrons; hereafter simply secondary “electrons”) in the discussion above. In p − p collision the total emissivities of secondary electrons and gamma-rays account for $\frac{\mathrm{1}}{\mathrm{6}}$ and $\frac{\mathrm{1}}{\mathrm{3}}$ of the total pion emissivity (e.g., Kelner et al. 2006), respectively. Meanwhile, the energy density of the magnetic fields is ${\mathit{U}}_{\mathrm{M}}\mathrm{~}\mathrm{2.5}\mathrm{(}\frac{\mathit{B}}{\mathrm{10}\mathit{\mu }\mathrm{G}}{\mathrm{)}}^{\mathrm{2}}$ eV/cm³ and the energy density of ambient photon fields is U_ph ~ 1 eV/cm³ as discussed above. Thus synchrotron cooling is the dominant electron cooling process and only about 30% of the secondary electrons energy is transferred to gamma-rays via ICs. As a result the gamma-ray flux from secondary electrons is only 15% of the gamma-ray flux from π⁰ decay and we neglect this contribution for simplicity.

A possible concern for an energy-dependent diffusion scenario is that, in the high energy range for which the diffusion length is much larger than the distance to the source, the predicted proton density should be the same at all latitudes for the impulsive case, or scale as 1 /R (where R is the distance to the injection source) for the continuous case. Thus, given the (expected) relatively smaller line of sight through the top of the SFB, one naively anticipates that the gamma-ray flux in South 4 is smaller than that in South 1. The observed data reveal a similar total flux in both slices at high energy, however. One possible resolution to this tension may be that the target gas density is higher at high latitudes. It is also likely that the geometry of the FBs is different from our simple assumption of sphericity thereby entailing non-trivial projection effects.

In a leptonic scenario the low energy deficit of gamma-ray flux in the top slice also requires a corresponding low energy break in the electron distribution (see Fig. 8 right). This can, again, be produced by several mechanisms. Firstly, the energy-dependent diffusion effects mentioned above also work for electrons in principle (cf. McQuinn & Zaldarriaga 2011). However the much faster cooling of electrons relative to protons tends to prevent such effects playing an important role. As discussed above in the circumstances of the FBs synchrotron cooling dominates ICs cooling, and bremmsstrahulung can be neglected due to the low density of the FB. To fit the high energy gamma-ray data we need TeV electrons, whose cooling time scale is about 10⁵ yr. Thus, even if we assume very fast diffusion with the diffusion coefficient of 4 × 10²⁹ cm²/s at 1 GeV (and an index of 0.6), the diffusion distance of TeV electrons in 10⁵ yr is only 2 kpc, significantly less than the distance of South 4 from the plane (which is more than 5 kpc). This implies high energy electrons accelerated at low latitudes can never reach South 4 and, as a result, if we want to explain the radiation of the FBs in a leptonic scenario, the electrons should be accelerated in situ.

The second problem related to elxectrons concerns the interpretation of the very hard energy spectrum of radiation below 5 GeV. Such a hard spectrum formally can be explained by assuming a low energy break in the electron spectrum. However, if the electrons are cooled, no matter how hard the injection spectrum is, because of the radiative (dE/dt ~ E² type) losses the resulted steady-state distribution below the break in the injection spectrum, will have have a standard, E^-2 power law-spectrum. The latter will result in IC gamma-ray spectrum with a photon index 1.5 which is still not sufficient to explain the observed gamma-ray spectrum of South 4; below 5 GeV the photon index of the latter is as small as 1. One can formally overcome this problem assuming that the electrons stay uncooled in FBs. However, in order to reproduce the break in the gamma-ray spectrum, one would need a corresponding break at 100 GeV. The corresponding synchrotron cooling time at this energy, $\mathit{t}\mathrm{~}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\mathrm{\left(}\frac{\mathrm{100}\mathrm{GeV}}{\mathit{E}}\mathrm{\right)}\mathrm{(}\frac{\mathit{B}}{\mathrm{10}\mathit{\mu }\mathrm{G}}{\mathrm{)}}^{-2}$ yr, appears significantly shorter than the age of FBs within any reasonable model of the latter. Thus we need very efficient continuous acceleration acceleration of electrons which would dominate over the rate of energy losses and in that way keep very hard the steady-state spectrum. This in principle can be realized with stochastic acceleration which can produce Maxwellian type steady-state energy distribution of electrons like N(E) ~ E²exp( − E/E₀). The value E₀ can be derived by equating the acceleration time of stochastic acceleration and synchrotron cooling time (e.g., Lefa et al. 2011), ${\mathit{E}}_{\mathrm{0}}\mathrm{~}\mathrm{1.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{\left(}\frac{{\mathit{v}}_{\mathit{A}}}{\mathrm{0.001}\mathit{c}}\mathrm{\right)}\mathrm{(}\frac{\mathit{B}}{\mathrm{10}\mathit{\mu }\mathrm{G}}{\mathrm{)}}^{-0.5}$ GeV for scattering off Alfven waves. Mertsch & Sarkar (2011) considered scattering off turbulence generated by plasma instabilities at a hypothesised shock front at the edge of the bubbles. They obtain an E₀ of several hundred GeV. Indeed adopting an Maxwellian type electron distribution with E₀ = 300 GeV, we can have a reasonable fit to the detected gamma-ray spectrum in South 4, as it is shown in Fig. 9. For comparison, we also show IC emission at the position of South 1 and South 4 by adopting a broken power law electron distribution with an index 2 below the break at 200 GeV and 3 above. The results fit the SED of South 1 well, but fail to fit that of South 4, which are also shown in Fig. 9.

Anisotropy effects may, in principle, increase the IC gamma-ray flux at high latitudes since most of seed photons come from regions in the Galactic plane; thus gamma-rays we detect at high latitudes are a result of (almost) head-on collisions. To probe this effect we use the formalism described in Moskalenko & Strong (2000). The enhancement factor is defined by the flux ratio of the anisotropic IC to the isotropic IC assuming the the same ISRF energy density. Results for the positions of South 1 and South 4 are shown in Fig. 9. We can see that although the anisotropic IC effect may cause a spectral a hardening, however the difference between South 1 and South 4 is rather small, less than 10%. It appears not sufficient to explain the very hard spectrum of South 4.

Fig. 8 Assumed CR proton (left) and electron (right) distributions in Slices 1 and 4 of the SFB.

Fig. 9 Left and middle panels: fitting the FBs’ SED with the hadronic model (left panel) and the leptonic model (middle panel) described in the text. The red curve shows the fitted spectrum and the red hatched region spans all spectra resulting from different templates for South 1; blue curve and hatched region are for South 4. For leptonic models, for South 1 we assume the broken power law described in the text and the ISRF value given by GALPROP. For South 4, the black curve is for the same broken power law with the ISRF floating while for the gray curve we use the ISRF in GALPROP at the position of South 4 and a Maxwellian type electron distribution with Ebk = 300 GeV. In our hadronic scenario energy-dependent diffusion effects are considered as described in the text. Right panel: the enhancement factor due to anisotropic IC. The IR and Opt/UV components of the ISRF have effective temperatures of 100 K and 5000 K, respectively. Results for z = 5 kpc and z = 7 kpc are shown (approximate heights of South 1 and South 4 above the Galactic plane).

4. Summary

In this paper we re-analyze the Fermi-LAT data covering the FB region. With an improved instrumental response function and more data, we are able to extend the previously obtained spectrum to lower energies. Furthermore we divide each bubble into slices to investigate possible variation of the spectrum with latitude. Given the improved data, for the first time, we can determined robustly that the spectrum of the top of south bubble drops appreciably at low energies relative to the spectrum determined at lower latitudes. We also show that the morphology of the south bubble is energy-dependent; at high energy the structure is relatively more extended to both Galactic south and west.

We have also investigated the influence of different choices of background model on the results of our analysis. We found that background models may significantly alter the apparent gamma-ray signal from the FBs. Nevertheless, the spectral hardening with latitude in the south bubble remains a robust result. In neither bubble do we find a spectral steepening with latitude.

The relative suppression of the low energy gamma-ray spectrum may be explained within a hadronic model for the FBs wherein energy-dependent diffusion leads to a relative deficit of low energy protons at high latitudes. Specifically, if protons are injected in the plane, the finite age of the FBs may only allow high energy protons to propagate to high latitude thereby predicting a gradual hardening of the gamma-ray spectrum with latitude.

In leptonic models the fast cooling of electrons means they cannot move too far from their accelerators and distributed acceleration inside the FBs seems to be favored. For the case of stochastic acceleration the maximal electron energy might be expected to be 1 TeV, implying the optical/UV and IR component of the ISRF may play an important role in gamma-ray production inside the FBs. However, the attenuation of these components of the ISRF at high latitudes is, at least naively, in conflict with the observed hardening of the gamma-ray spectrum in the highest slice of the south bubble unless a very specific electron spectrum evolution is realized.

Acknowledgments

R.M.C. was supported by grant FT110100108 from the Australian Research Council.

References

All Figures

	Fig. 1 Residual map above 2 GeV. Background subtraction is described in the text. Two bubble-like structures can be seen. To render the picture clearer we mask the bright Galactic plane | b | < 5°.
In the text

	Fig. 2 SEDs of different slices of the north and south FBs. The numbers 1 to 4 run from low to high latitudes. Left panel: north bubble (NFB). Right panel: south bubble (SFB).
In the text

	Fig. 3 Residual maps for different energy bins. The left and right panels correspond to energy intervals 1–2 GeV and 10–30 GeV, respectively. The inner disk of the Galaxy (| b | < 2°) is masked. The boxes laid over the residual maps show the position of each slice.
In the text

	Fig. 4 Zoom in of Fig. 3 for the SFB region. Left panel 1–2 GeV and right panel 10–30 GeV.
In the text

	Fig. 5 SEDs for the seven slices with all 64 templates we use in our analysis. The shaded areas span all derived SEDs. The left panel is for the NFB while the right is for the SFB.
In the text

	Fig. 6 Six examples of SEDs for slices in the SFB.
In the text

	Fig. 7 Six examples of SEDs for slices in the NFB.
In the text

	Fig. 8 Assumed CR proton (left) and electron (right) distributions in Slices 1 and 4 of the SFB.
In the text

Fig. 9 Left and middle panels: fitting the FBs’ SED with the hadronic model (left panel) and the leptonic model (middle panel) described in the text. The red curve shows the fitted spectrum and the red hatched region spans all spectra resulting from different templates for South 1; blue curve and hatched region are for South 4. For leptonic models, for South 1 we assume the broken power law described in the text and the ISRF value given by GALPROP. For South 4, the black curve is for the same broken power law with the ISRF floating while for the gray curve we use the ISRF in GALPROP at the position of South 4 and a Maxwellian type electron distribution with Ebk = 300 GeV. In our hadronic scenario energy-dependent diffusion effects are considered as described in the text. Right panel: the enhancement factor due to anisotropic IC. The IR and Opt/UV components of the ISRF have effective temperatures of 100 K and 5000 K, respectively. Results for z = 5 kpc and z = 7 kpc are shown (approximate heights of South 1 and South 4 above the Galactic plane).

In the text