Fermi-bubble Bulk and Edge Analysis Reveals Dust, Cooling Breaks, and Cosmic-Ray Diffusion, Facilitating a Self-consistent Model

The full, radio to γ-ray spectrum of the Fermi bubbles is shown to be consistent with standard strong-shock electron acceleration at the bubble edge, without the unnatural energy cutoffs and unrealistic electron cooling of previous studies, if the ambient interstellar radiation is strong; the γ-ray cooling break should then have a microwave counterpart, undetected until now. Indeed, a broadband bubble-edge analysis uncovers a pronounced downstream dust component, which masked the anticipated ∼35 GHz spectral break, and the same overall radio softening consistent with Kraichnan diffusion previously reported in γ-rays. A self-consistent bulk and edge model implies a few Myr old bubbles, with fairly uniform radiation fields and enhanced magnetization near the edge.

Figure 1 shows the specific I ϵ γ-ray and I ν microwave intensities, where ϵ and ν are the photon energy and frequency.The above arguments indicate that the γ-rays arise from the inverse-Compton (IC) scattering of the ambient radiation field by the CREs.Starlight photons of mean ε ≃ 1 eV energy 1 dominate the interstellar radiation budget above the Galactic center, out to a height of at least 5 kpc (e.g., Ackermann et al. 2014;Popescu et al. 2017).The massive FB-driven outflows, not incorporated in typical radiation models (which may underestimate starlight even near the Galactic center; e.g., Ajello et al. 2016), have likely elevated the starlight and infrared energy densities considerably, out to the ∼ 10 kpc height of the bubbles, by removing much of the ambient absorbers.Hence, the ϵ ∼ 50 GeV break is naturally explained as the ϵ ≃ 3(m e c 2 ) 2 /16ε ≃ 50ε −1 1 GeV, Klein-Nishina (KN; e.g., Rybicki & Lightman 1979) suppression of IC-scattered starlight, whereas the ϵ ≳ few GeV softening is shown below to be the signature of a cooling break.Here, m e is the electron mass, c is the speed of light, and ε 1 ≡ ε/1 eV.In such a model, the ∼ 50 GeV break should become less pronounced at higher latitudes, where it is increasingly washed out by the spectrally-flat, IC-scattered cosmic microwave background (CMB) photons; indeed, such a trend is hinted by spatially resolved (Hooper & Slatyer 2013) and edge (Keshet & Gurwich 2017) FB spectra.

Leptonic model
We assume that CREs are injected at the shock with index p ≃ (r + 2)/(r − 1), thought to depend only on the shock compression ratio r (Axford et al. 1977;Krymskii 1977;Bell 1978;Blandford & Ostriker 1978) if scattering is sufficiently isotropic (Keshet et al. 2020), implying a shock velocity Here, m ≃ 0.59m p and m p are the mean and proton masses, T u ≡ 0.2T 0.2 keV/k B (Miller & Bregman 2013;Kataoka et al. 2013;Keshet & Gurwich 2018) is the upstream temperature, k B is Boltzmann's constant, and we adopted an adiabatic index Γ = 5/3 in the last, p ≃ 2.1 approximation where q ≡ (p−2)/0.1 approaches unity.Such a p ≃ 2.1 spectrum reproduces the hard synchrotron fit in Fig. 1 for non-cooled CREs, and the flat ∼ 3-30 GeV spectrum for cooled CREs.The FB age for a simple, R ∝ t 2/3 Primakoff scaling (Courant & Friedrichs 1948;Keller 1956;Bernstein & Book 1980) of the R ≡ 10R 10 kpc FB height is then (Keshet & Gurwich 2018) Myr . (2) FB CREs cool by synchrotron and IC emission, thus softening to an index p + 1 above the cooling-break energy implying a synchrotron break at frequency and an IC-scattered starlight break at energy where the energy density accounts respectively for the ambient ε s ≃ 1 eV starlight, ε i ≃ 10 −2 eV infrared (e.g., Ackermann et al. 2014), ε c ≃ 10 −3 eV CMB, and magnetic, components.
Here, t 3 ≡ t/3 Myr, B 3 ≡ B/3µG is the normalized magnetic field, e is the electron charge, and σ T is the Thomson cross-section.The ratio between synchrotron emissivity at frequency ν from a non-cooled, n ∝ E −p CRE distribution, and starlight-IC emissivity at energy ϵ from the cooled CREs, is where we defined ν 10 ≡ ν/10 GHz and ϵ 10 ≡ ϵ/10 GeV; the dimensionless function ϕ(p) is given in Appendix §A.
Adopting the u i + u c ≪ u s limit (henceforth), the ratio between radio and ≳ 10 GeV γ-ray intensities should approximately be given by Eq. ( 7).
Figure 1 demonstrates a few one-zone leptonic-model fits to the bulk γ-ray spectrum, constrained also by the ν ≃ 23 GHz WMAP data point.Here, we fix p = 2.1 and u i = 0.1 eV cm −3 , and vary the other parameters; additional model fits are provided in Appendix §A.Interestingly, for the best-fit parameters, we find that the cooling break and KN turnaround approximately merge to yield a stronger spectral break.While the u s = 0.5 eV cm −3 fit marginally agrees, within statistical and systematic uncertainties, with the γ-ray data, much better fits are obtained for 1 eV cm −3 ≲ u s ≲ 3 eV cm −3 .Such starlight energy densities are higher than usually assumed at R ≳ 4 kpc heights, but are plausible as the mean starlight inside the FB outflows, given conservative u s (R = 2 kpc) ≃ 1.4 eV cm −3 estimates (Acker-mann et al. 2014), as discussed above and in §6.For such u s values, the γ-ray fit is surprisingly good, considering the oversimplified, one-zone approximation, and can be arbitrarily improved if this assumption is relaxed.In contrast, the microwave fit is poor in all our successful γ-ray models, as they require a ν b < 40 GHz cooling break.These quantitative conclusions are in line with the above qualitative estimates, and are not sensitive to our assumptions; in particular, we obtain similar results for different 2.0 ≤ p ≲ 2.2 CRE spectra, different u i < 0.2u s prescriptions, and two-zone models.For p ≳ 2.3, the γ-ray spectral fit becomes unacceptable for plausible t and u s values, as shown in Appendix §A.

Bubble-edge analysis
Next, we analyze the FB edges, identified by gradient filters (Keshet & Gurwich 2017), and bin the data parallel to the edge in the same method used to measure the γ-ray (Keshet & Gurwich 2017) and X-ray (Keshet & Gurwich 2018) spectra in the near downstream.Here, we study lower frequencies, ranging from 45 MHz radio to 60 µm infrared, focusing in particular on cleaned microwave WMAP and Planck maps with point sources masked and the CMB subtracted, as described in Appendix §B. Figure 2 demonstrates (RGB colors) the γray, infrared, and microwave sky in and around the FBs, and outlines (white contours) the north and south sectors used in our analysis.These sectors, crossing the FB edge at the top of each bubble, were chosen to contain large upstream and downstream regions with minimal extended contamination.
Figure 3 presents the spectrum of the downstream regions of both north (upper panels) and south (bottom panels) FBs, obtained by binning the data parallel to the edge and subtracting the upstream-based foreground and background (henceforth foreground).The spatial distribution is demonstrated (left panels) for select channels, all found to sharply brighten as one crosses downstream, facilitating a measurement of the spectrum (right panels) at different outward-oriented angular distances ψ from the edge.While the signal strengthens downstream with increasing |ψ|, the spectral shape is found to be robust as a function of ψ, and quite similar in the north and south bubbles, substantiating the validity of the measurement.As expected from previous studies, the north bubble shows more surrounding substructure, both upstream (Loop-I) and downstream (mainly at low frequencies); it also presents a more gradual γ-ray jump, possibly attributed in part to edge misalignment.
The cleaner, south bubble provides a good measurement of the synchrotron signal, found for the first time to persist down to frequencies as low as 45 MHz.In the present method, the spectral index can be accurately determined to be α = −0.77± 0.03 downstream (see Table 1).This index is softer by an η = 0.20 ± 0.07 offset with respect to the α = −0.57± 0.06 bulk synchrotron index (Fig. 1), consistent with both Kraichnan and Kolmogorov diffusion of escaping CREs.A similar but slightly stronger softening was identified in the IC γ-ray emission from the edge, and was interpreted as arising from CRE diffusion (Keshet & Gurwich 2017); its detection in radio validates this diffusive model.
A pronounced hard signal is found to dominate the downstream spectrum above ∼ 100 GHz, extending up to an exponential turnaround in the infrared.The signal is well-fitted by an I ν ∝ ν 7/2 /(e hν/k B T −1) spectrum (solid curves in Fig. 3), and is thus identified as thermal dust emission, of temperature T ≃ 29 ± 1 K (33 ± 1 K) in the north (south) FB (see Table 1).Such temperatures, already hinted by coincident hot dust regions in Planck maps, may underestimate the true postshock T due to the inaccurate subtraction of ≲ 20 K (Planck Collaboration et al. 2016) foreground dust.As the left panels of Fig. 3 show, the radial profiles of dust emission are nearly identical to those of IC and especially synchrotron emission, in both bubbles.We conclude that this dust component is directly associated with the FB a The downstream regions are defined in Fig. 3; in north bubble regions, ν < GHz data were excluded from the fit due to interfering substructure.
shock, probably due to dust entrainment or formation in the postshock region, in resemblance of AGB stars, novae, and core-collapse supernovae.Indeed, Fig. 2 shows some of the dust emission (green and blue) extending to the top of the γ-ray (red) bubbles.A dust FB component should therefore be included in template analyses that aim at measuring the microwave spectrum.

Joint bulk and edge model
The unexpected, hard microwave dust-emission downstream explains in part why a cooling break was not identified until now.Indeed, Fig. 3 shows that if the dust model is subtracted from the data, the remaining synchrotron component (dashed curves) presents a spectral break (insets) around ∼ 35 GHz, as anticipated above.In accordance, fitting a superposition of synchrotron emission of index α and a frequency ν b cooling break, and thermal dust emission of temperature T , as summarized in Table 1, gives ν b = 34 ± 5 GHz in the southern bubble.Curiously, both WMAP (Dobler 2012b) and joint WMAP and Planck (Planck Collaboration 2013) template-based analyses found a similar bulk spectral break in the total synchrotron signal, but attributed it to combined hard+soft foregrounds rather than to the bubbles.
While the local spectral measurement near the edge does not require strong assumptions concerning extended foregrounds, and does not suffer from templateremoval systematics (Keshet & Gurwich 2017, 2018) which can be highly detrimental (e.g., Keshet et al. 2004), it does show larger statistical errors than the volume-integrated, bulk spectrum.Our results indicate that this bulk spectrum too must show a cool-ing break, which would fix the FB parameters even without using the hydrodynamic Eq. ( 2).Requiring similar FB ages inferred from edge and bulk spectra suggests a ν bulk ≃ 23 GHz cooling break in the latter, unidentified so far due to modelling systematics and unaccounted FB dust.We then obtain for the bulk t ≃ 5.2ε ν 1.55 b,23 Φ −1 0.1 eV cm −3 , u/u s ≃ 1 + 0.077ν 0.45 b,23 Φ 0.1 , and where we defined ν b,23 ≡ ν bulk /23 GHz.Note that while a ν bulk < 60 GHz cooling break would modify the previously reported α ≃ 0.55 synchrotron spectrum, a p ≃ 2.1 CRE spectrum is nevertheless needed for a self-consistent model, as well as for an acceptable γ-ray fit (see Appendix A).

Discussion
We conclude that the spectra of both north and south FB edges, and of the bulk, both with and without hydrodynamic considerations, produce consistent results, thus providing strong constraints on the FBs and their environment.The results indicate younger FBs than obtained in early spectral models, an enhanced magnetic field near the edges, and strong and fairly uniform radiation fields throughout the bubbles.Indeed, while we cannot directly measure the ambient starlight intensity, Fig. 3 shows an infrared νu ν ≃ 0.1 eV cm −3 peak at high-latitudes, very similar to its estimated value near the Galactic plane (e.g., νu ν ≃ 0.15 eV cm −3 at R = 2 kpc; Ackermann et al. 2014).
The success of our oversimplified, one-zone bulk model, which attributes the synchrotron and IC signals to the same pure power-law CRE injection, in-validates the justification for previous hadronic or leptonic models, which invoked ad-hoc cutoffs on the CR spectrum.Without ad-hoc assumptions, future better measurements of the broadband, radio to γ-ray signal would robustly constrain the CRE distribution and diffusion, coincident radiation fields, magnetic fields, and FB dynamics and energetics.The edge spectra presented in this work validate the bulk model, directly confirming the CRE softening expected from diffusion (Keshet & Gurwich 2017), the anticipated microwave cooling break, and, indirectly, the persistence of strong radiation fields throughout the FBs.  1, but for p = 2 and model parameters {us cm 3 /eV, t/Myr, B/µG} = {0.5, 10, 1.7} (solid blue), {1, 8, 2} (dashed green), {2, 3.6, 2.6} (dot-dashed red), and {3, 2.7, 3} (double-dot-dashed orange).
Binning along the FB edge removes most foregrounds, but to obtain better results, we first subtract the CMB contribution, estimated using the Planck HFI internal linear combination (PILC; Planck Collaboration 2013).This map is constructed by a superposition of 143-545 GHz maps, designed to spectrally isolate the CMB from the dust template of Finkbeiner et al. (1999).
We mask the point sources from WMAP 's nine-year point-source catalog and Planck 's second catalogue of point sources (PCCS2) in the relevant bands, and the following large-scale structures: the Small and Large Magellanic Clouds, M31, NGC5090, NGC5128, Orion-Barnard Loop, and ζ Oph.We also mask pixels where dust extinction at Hα frequencies is larger than 1 mag, and pixels where the Hα intensity is greater than 10 Rayleigh.All maps are masked at the highest, N side = 2048 HEALPix (Górski et al. 2005) resolution available.

C. Fermi-LAT sky map
We use the archival Pass-8 LAT data from the Fermi Science Support Center (FSSC) 5 , and the Fermi Science Tools (version v10r0p5).Weekly all-sky files are used, spanning weeks 9 through 789 for a total of 781 weeks (∼ 15 years), with ULTRACLEANVETO class photon events.A 90 • zenith angle cut is applied to avoid CR-generated γ-rays from the Earth's atmospheric limb, according to the appropriate FSSC Data Preparation recommendations.Good time intervals are selected using the recommended expression (DATA QUAL==1) and (LAT CONFIG==1).Point-source contamination is minimized by masking pixels within the 95% total-event containment area of each point source in the LAT fourth source catalog (4FGL; Abdollahi et al. 2020).

Figure 3 .
Figure 3. Spectral analysis of north (top panels) and south (bottom) FB edges.Left panels show (error bars) the profiles of select (see legend) channels (normalized by their maxima) and the 23 GHz to 10 GeV brightness ratio Φ d , as a function of the angular distance ψ from the edge, after subtracting the foreground based on the upstream (ψ > 0) gray-shaded region.Right panels show (symbols enclosing error bars, with dotted lines to guide the eye) the spectrum of excess emission in the near (red diamonds), mid (green squares), and far (blue circles) downstream regions designated by the respective color-shaded regions in the left panels.Also shown are the synchrotron profiles (dashed curves) obtained by subtracting the respective best-fitted dust models (see text; solid curves) from the downstream excess, and a synchrotron model (of arbitrary normalization; labeled dot-dashed black) of p = 2 CRE injection, Kraichnan diffusion, and a ν b = 35 GHz cooling break (black disk); the insets focus on the cooling-break region in the mid-downstream.Data are extracted from 45 MHz(Alvarez et al. 1997;Maeda et al. 1999), reprocessed 408 MHz(Haslam et al. 1982;Remazeilles et al. 2015), 1.4 GHz(Reich 1982;Reich & Reich 1986;Reich et al. 2001), WMAP and Planck microwave, 160 µm and 140 µm AKARI(Doi et al. 2015), 100 µm and 60 µm IRIS (Miville-Deschênes & Lagache 2005), and 3-30 GeV Fermi-LAT maps.