A UNIFIED MODEL OF THE FERMI BUBBLES, MICROWAVE HAZE, AND POLARIZED RADIO LOBES: REVERSE SHOCKS IN THE GALACTIC CENTER'S GIANT OUTFLOWS

The Fermi bubbles are enigmatic, giant γ-ray structures emanating north and south from the nucleus of the Milky Way (Dobler et al. 2010; Su et al. 2010; Su & Finkbeiner 2012; Ackermann et al. 2014; Yang et al. 2014) and extending to heights above the Galactic plane approaching 8 kpc. Features coincident with the bubbles or some component of the bubbles have been discovered in total-intensity radio and microwave emission, polarized radio and microwave emission, and X-ray emission (Bland-Hawthorn & Cohen 2003; Finkbeiner 2004; Dobler & Finkbeiner 2008; Dobler 2012a; Ade et al. 2013; Carretti et al. 2013; Kataoka et al. 2013; Tahara et al. 2015).

The energy source of the bubbles remains in contention, as does the γ-ray emission process that illuminates them. Broadly speaking, the bubbles might be the signatures of a fairly recent (within a ∼ Myr timescale) explosive outburst (Su et al. 2010; Mertsch & Sarkar 2011; Zubovas et al. 2011; Guo & Mathews 2012; Yang et al. 2012; Zubovas & Nayakshin 2012; Fujita et al. 2013; Barkov & Bosch-Ramon 2014; Mou et al. 2014; Wardle & Yusef-Zadeh 2014) from the central supermassive black hole with their γ-ray emission originating in inverse Compton (IC) collisions of a nonthermal electron population accelerated in association with this activity. Alternatively, the bubbles may result from the integrated effect of more or less secular processes associated with the inner Galaxy (Thoudam 2013) or Galactic nucleus, e.g., tidal disruption events that occur there every 10⁴–10⁵ yr (Cheng et al. 2011, 2015) or the concentrated and sustained star formation activity (e.g., Lacki 2014) that occurs throughout the nuclear region (inner 200–300 pc diameter region around the black hole).

Invoking sustained nuclear processes as the origin of the bubbles complements the possibility that their γ-ray emission originates not from IC emission but from hadronic collisions experienced by a population of cosmic-ray protons and heavier ions (Crocker & Aharonian 2011; Crocker 2012; Crocker et al. 2014). Such particles lose energy much more slowly than the high-energy electrons requisite for IC γ-ray production, thereby allowing the total current energy content of the bubbles to be accumulated more gradually from relatively less powerful processes such as star formation (the Eddington luminosity of the Galaxy's 4 million solar mass black hole is 5 × 10⁴⁴ erg s⁻¹; in contrast, the Milky Way's nuclear star formation releases a mechanical power of ∼3 × 10⁴⁰ erg s⁻¹; see below).

At present, there are strengths and weaknesses to either of these broad scenarios. IC scenarios share the attractive feature that the very same electron population that putatively supplies the IC γ-rays can, for magnetic field strengths approaching 10 μG, also supply the hard microwave haze emission via synchrotron radiation (Dobler et al. 2010; Su et al. 2010; Ackermann et al. 2014). On the other hand, if described by a single power law, synchrotron emission from this population of electrons would not explain the comparatively steep radio spectrum measured between 2.3 and 23 GHz for the bubbles' polarized, radio lobe counterparts (Carretti et al. 2013). Furthermore, the fact that the γ-ray spectrum of the bubbles does not steepen (Ackermann et al. 2014)—rather, it seems to harden (Selig et al. 2014; Yang et al. 2014)—with increasing Galactic latitude is difficult to explain naturally within IC scenarios given the attenuation of the optical and infrared background relative to the cosmic microwave background (CMB) as the distance from the Galaxy increases.

Another apparent weakness of many explosive, active galactic nucleus (AGN) related scenarios is that they require that we are viewing the bubbles at a privileged time at most a few million years since the outburst that created them, with the bubbles currently undergoing rapid expansion at speeds ≳2000 km s⁻¹. Rapid expansion is required in order either to transport the ≳TeV electrons (required to produce ≳GeV γ-rays via IC on the CMB) from the nucleus to the bubbles' peripheries within their ≲10⁶ yr loss times (in which case speeds ≳7000 km s⁻¹ ≃ 0.02 × c are actually required) or to have sufficiently energetic shocks at the expanding edges of the bubbles that cascading turbulence injected there may energize similarly energetic electrons in situ (Mertsch & Sarkar 2011). Yet X-ray observations of the bubbles have only indicated weak features associated with the bubble edges (Su et al. 2010; Kataoka et al. 2013; Fang & Jiang 2014; Tahara et al. 2015), implying rather mild expansion speeds, ≲300 km s⁻¹ (close to the sound speed). Closer to the base of the bubbles, data obtained in UV absorption observations by Fox et al. (2014) along the (l, b) = (104, +112) sightline (to quasar PDS 456) reveal cool (∼10⁴ K) gas entrained into a biconical, ∼900 km s⁻¹ nuclear outflow 2.5–4.0 Myr ago.

Moreover, there are other indications that the bubbles are rather long-lived structures.

• 1.  
  The fact that both bubbles lean away from the l = 0 meridian toward Galactic west and, moreover, there are extended plumes of polarized, 2.3 GHz emission leaving northwest (south west) from the high-latitude regions of the north (south) bubble (Carretti et al. 2013) suggests that the structures are distorted by a pressure gradient within the Galactic halo set up by the motion of the Milky Way through the local medium toward Andromeda. The ∼50 km s⁻¹ relative speed (van der Marel et al. 2012) of this motion implies a rather gentle effect that requires long timescales/low bubble gas velocities to take effect, consistent with the ∼100–150 km s⁻¹ flow velocities we determine below for gas in the bubbles (see Yang et al. 2012, where it was shown that an unreasonably large speed of the Galaxy through the local medium was required to produce the bubbles' distortion in the case that they are expanding rapidly owing to AGN activity).
• 2.  
  The ∼100 km s⁻¹ rotation of the bubbles' waist in the plane is apparently inscribed into giant radio ridges on the near surface of the bubbles (Carretti et al. 2013), again requiring that the structures not grow too quickly (here meaning ≲10⁷ yr).
• 3.  
  The steep spectrum of the bubbles' polarized radio counterparts between 2.3 and 23 GHz suggests that electrons emitting in this frequency range have had time to cool, mix, and accumulate; this requires ≳3 × 10⁷ yr for a reasonable magnetic field amplitude.
• 4.  
  Finally, observations of UV absorption features toward background quasar and AGB stars have established the existence of individual parcels of ionized gas at low Galactic longitude but varying latitude (Keeney et al. 2006; Zech et al. 2008). Kinematic analysis points to this gas participating in a nuclear fountain flow. Individual gas parcels seem to have been launched from the nucleus at different times in the past, ranging from 20–50 Myr to more than 800 Myr. The parcels are all now either rising to or falling from a maximum height of ∼12 kpc, similar to the vertical extent of the bubbles. Some of this material (Zech et al. 2008) appears significantly supersolar in metallicity. This phenomenology points to the existence of a long-standing, star-formation-driven, nuclear fountain; as we expand on later in the paper, we believe that this material is condensing in the shell (R3 in Figure 1) leading up to the contact discontinuity (CD) in each bubble.

Zoom In Zoom Out Reset image size

Turning now to the star-formation-driven hadronic model (Crocker & Aharonian 2011; Crocker 2012; Crocker et al. 2014), this has the appealing feature that the power currently being delivered into cosmic-ray ions by nuclear star formation is elegantly sufficient to explain the bubbles' (putatively hadronic) γ-ray luminosity in steady state. On the other hand, in the original instantiation of this scenario, the low density of target gas nuclei in the bubble plasma implied a very long proton–proton loss time of a ≳few Gyr, this being the timescale required to establish a steady state in the first place (i.e., approximately the timescale required to accumulate a sufficient total energy in cosmic-ray ions in the bubbles to power their instantaneous γ-ray luminosity given the low gas density). Some of us have recently shown, however, that accumulation of hadronic cosmic rays into gas condensations formed by local thermal instability in the bubbles can significantly reduce the timescale associated with the hadronic scenario to ∼few × 10⁸ yr (Crocker et al. 2014). Another point concerning the hadronic scenario is that although secondary electrons (and positrons) are a natural concomitant of hadronic γ-ray production, it is difficult to explain the spectrum of the microwave haze with such secondaries since IC and synchrotron losses would be expected, in a single-zone model, to steepen the steady-state secondary distribution such that it is inconsistent with the hard haze emission (e.g., Ackermann et al. 2014; Cheng et al. 2014b). It seems, then, that even the hadronic scenario requires primary electrons at some level to explain the totality of giant outflow phenomenology.

Indeed, the power in cosmic-ray primary electrons accelerated in concert with nuclear star formation is a good match to the bubbles' 2.3 GHz synchrotron luminosity (Carretti et al. 2013). The complement to these points is that the nucleus is underluminous in both diffuse γ-ray and radio continuum emission given how much star formation occurs there (Crocker et al. 2011a, 2011b; Crocker 2012), and this seems to be because both cosmic-ray electrons and ions are advected out of the region in the large-scale outflow before they can radiate. Thus, these particles radiate on the size scale of the bubbles into which the nuclear outflow delivers them. As we show below, however, the passage from nucleus to bubbles is not a smooth journey for these particles.

Here we briefly explain our model for the bubbles' structure and nonthermal emission phenomenology. A schematic showing the configuration of the northern bubble according to our model is shown in Figure 1 for reference.

• 1.  
  Co-entrained hot plasma, magnetic fields, and cosmic rays are energized within the central ∼100 pc (in radius) region of the Galactic nucleus. The pressure gradient supplied by these interstellar medium (ISM) phases (with the hot plasma dynamically dominant) accelerates a nuclear wind that quickly reaches ≳500 km s⁻¹. The dense and massive molecular gas torus circumscribing the injection region collimates this outflow and directs it perpendicular to the plane. Shocks driven into the hot ISM, which is more dense at low latitudes, may account for the X-ray edges observed by Bland-Hawthorn & Cohen (2003). Radiative shocks driven into atomic and/or molecular gas may also contribute.
• 2.  
  The bubbles are inflated by the injection of hot plasma in the nuclear outflow over timescales of at least a few × 10⁸ yr. Characteristic scales for the mean mechanical power and mass flux into the bases of the bubbles are ∼10⁴⁰ erg s⁻¹ and ∼0.1 M_⊙ yr⁻¹, respectively.
• 3.  
  This outflow, however, is not powerful enough—given the strong gravitational field, the thermal pressure of the overlying halo gas, and the occurrence of gas cooling—to break out of the Galaxy.
• 4.  
  Instead, the outflow generates the two Fermi bubbles and associated phenomena. The bubbles inflate against an atmosphere of non-negligible pressure. An important ingredient of our scenario is that recent work indicates that the gas density and pressure in the halo, especially in the vicinity of the bubbles, are larger than has often been assumed in the past (e.g., Stocke et al. 2013; Miller & Bregman 2014).⁷
• 5.  
  The bubbles develop the classical wind-blown bubble structure of (starting from the smallest scales and working outward) (i) a reverse shock, (ii) a CD, and (iii) a forward shock; see Figure 3.
• 6.  
  The CD separates ISM material swept up and shocked in the downstream of the forward shock from shocked bubble material in the downstream of the internal reverse shock. (In our model, the γ-ray edge of the bubbles is identified with the CD.)⁸
• 7.  
  The bubbles quickly (within ∼10⁷ yr) reach pressure equilibrium with their environment; at this point (neglecting, as we do in our current model, cooling in the region between forward shock and CD) the forward shock degenerates into a sound wave. The CD decelerates even faster and is subsonic after only ≲10⁶ yr; at present it is expanding at only ∼ few km s⁻¹.
• 8.  
  Thermal plasma and cosmic rays initially suffer adiabatic and radiative energy losses as they are blown out of the nucleus on the conically expanding outflow (R1 in Figure 1). However, they encounter giant reverse shocks (one at ∼1 kpc height in each bubble) where the ram pressure of this flow equilibrates with the thermal pressure of the shocked, downstream gas (contrast the model of Lacki 2014, where the ram pressure of the expanding outflow equilibrates with the thermal pressure of the background Galactic halo atmosphere). At these shocks, the plasma and cosmic rays are reheated/reaccelerated. (In the remainder of this paper, unless otherwise noted, the shocks we refer to generically are these reverse shocks interior to the CD.)
• 9.  
  In the immediate downstream of the shocks, there is a systematic subsonic flow of ${v}_{\mathrm{down}\ }\sim 1/4{v}_{{\rm{w}}}\;\sim$ 100–150 km s⁻¹ (with v_w the wind speed immediately upstream of the shock). There is energy-independent transport of the nonthermal particles into the bubbles mediated by this flow in combination with cosmic-ray streaming at approximately the Alfvén speed. High-energy cosmic-ray electrons produce IC and synchrotron emission on photon and magnetic fields in the shock downstream. Given that there is a systematic flow (such that, at any particular distance from the shock, all electrons are the same "age"), the spectrum of these electrons (and their resultant emission) is hard (mirroring the injection distribution) up to some time-dependent (and thus distance-dependent) maximal energy above which the spectrum is exponentially suppressed. Given the rather slow transport speeds and relatively short cooling times, microwave-synchrotron-emitting electrons reach only ∼1–2 kpc past the shocks and ≳GeV IC-emitting electrons only reach a few × 10² pc. (In contrast, cosmic-ray protons and 2.3 GHz synchrotron-emitting electrons survive to reach the full extent of the bubbles.) Altogether, we show below that the hard-spectrum synchrotron radiation from cosmic-ray electrons reaccelerated on the shocks and transported downstream provides an excellent match to the luminosity, spectrum, distribution, and morphology of the microwave haze.
• 10.  
  As the bubbles age and their volume slowly increases (more than compensating for their density decline), the radiative losses of the plasma gradually become more significant. In fact, there is a steady-state limit (at constant radius and mass) where the bubbles reach pressure equilibrium with the environment and the cooling radiation from their interior plasma completely saturates the injected mechanical energy. As this limit is approached, an increasingly large fraction of the hot plasma mass injected into the bubbles is lost by dropout as cooling leads to local thermal instability that causes individual parcels of gas to cool, collapse, lose buoyancy, and fall out of the bubbles under gravity.
• 11.  
  Because of the density-squared dependence of cooling, this mass dropout process dominantly occurs in the somewhat overdense shell of material that forms inside of the CD, where the slow flow of plasma away from the reverse shock is gradually arrested (R3 in Figure 1).
• 12.  
  While the plasma evolves radiatively in the shell (as it cools down into localized, overdense condensations), cosmic-ray hadrons and magnetic fields, in constrast, continue to evolve adiabatically. Their compression in the shell material leads to adiabatic energy gain; they are compressed until they equilibrate with the thermal pressure of the volume-filling plasma in R2 between the shell and the reverse shock. In other words, adiabatic evolution of these nonthermal phases leads them to become the dominant dynamic agent in the shell.
• 13.  
  The ultimate extent and mass of shell R3 are determined by the pressure balance condition explained above and the further condition that the shell's cooling radiation saturates the enthalpy flux represented by the slow flow of hot plasma into it from R2.
• 14.  
  Given that both the mass of shell target material and the hadronic cosmic-ray energy density in the shell can be determined from these considerations, we can calculate the hadronic emission from the shell; we find that this essentially saturates at the observed γ-ray luminosity of the bubbles for reasonable parameter values. The fact that the dominant hadronic emission occurs mostly in the shell explains why the volumetric emission from the bubbles peaks toward their edges (as required to reproduce their flat projected surface brightness; Su et al. 2010).
• 15.  
  As for the protons and heavier ions, the lower-energy cosmic-ray electrons that live long enough to reach the shell are also adiabatically compressed and reenergized there. At the same time, collisions between the adiabatically compressed hadron population and the shell gas inject secondary electrons (+positrons) into R3. Unlike for the region downstream of the reverse shocks, however, in the shells there is a mixture of electron populations with different ages and thus different characteristic energies. This leads to the integrated, steady-state shell electron spectrum being steepened by the canonical difference in spectral index of 1 associated with synchrotron/IC losses. This cooled primary + locally produced secondary electron population synchrotron-radiates on the adiabatically compressed shell magnetic field, explaining the steep-spectrum, polarized emission measured between 2.3 and 23 GHz.
• 16.  
  The shell/CD does not constitute a perfect boundary, however: observationally, in 2.3 GHz polarization data, there are plumes of radio emission that extend northwest and southwest beyond the γ-ray boundaries of the bubbles (R4). Recent X-ray evidence (Tahara et al. 2015) seems to indicate hot plasma venting from the north bubbles, roughly coincident with the northern radio plume. We speculate that these plumes are symptomatic of the breakout of the light fluid of co-mixed cosmic rays, plasma, and magnetic fields through the heavy shell, mediated by the Rayleigh–Taylor instability. However, the CD generates the (required) γ-ray edge because, even if some cosmic-ray protons also escape through this boundary, this delimits the region of dense target gas accumulated into the shell from the interior flow.

The remainder of this paper is structured as follows. We describe the input data that frame the problem in Section 2. In Section 3 we show that, incorporating the effects of gravity, the bubbles should contain giant, reverse shocks at heights of ∼1 kpc into the outflows north and south from the nucleus. In Section 4 we set out a semianalytic calculation describing the evolution of a bubble fed by an outflow of fixed mechanical power and mass flux expanding into an atmosphere of finite, non-negligible pressure. We show that, adopting values for these quantities apposite to the Fermi bubbles, such a bubble reaches the size scale of the Fermi bubbles after a few ×10⁸ yr and that the asymptotic radius of such a bubble for these parameters is not much larger than that of the observed Fermi bubbles. This finding removes any "Why now?" problem and justifies our subsequent steady-state treatment of the bubbles' nonthermal emission. In Section 5 we describe our procedure for delimiting the region of the parameter space defined by the energy and mass fluxes (${\dot{E}}_{0}$ and $\dot{M}$, respectively) injected at the nucleus that describe the present-day bubbles. In Section 6 we demonstrate explicitly that, in the small region of parameter space consistent with prior constraints, the bubbles' (i) predicted hadronic γ-ray emission, (ii) predicted microwave synchrotron emission, and (iii) radio continuum synchrotron emission all match the observed values. In other words, we show that our scenario consistently explains, in a unified model, the nonthermal phenomena described as the "Fermi bubbles," the "microwave haze," and the "S-PASS lobes." In Section 7 we review specific observations that seem to support (or, in one case, may challenge) our model. Finally, in Section 8 we set out implications of our model and future work.

Adopting a distance to the Galactic center of 8 kpc, geometrically the bubbles can be approximated as a pair of spheres with radii ∼3.2 kpc with a total volume of ${V}_{\mathrm{FB}}\simeq 8.4\times {10}^{66}$ cm³ and a combined surface area of ${A}_{\mathrm{FB}}\simeq 2.5\times {10}^{45}$ cm². A better characterization of the bubbles' real geometry seems to be that they are bi-cones emanating north and south from the nucleus to heights of 2–3 kpc above the plane, above which they collimate into cylinders (Bland-Hawthorn & Cohen 2003; Carretti et al. 2013). The Fermi bubbles reach ∼7 kpc maximal height above the Galactic center accounting for projection effects; the polarized radio lobes extend up to ∼8 kpc. A schematic showing the configuration of the northern bubble according to our model is shown in Figure 1.

Kataoka et al. (2013) have performed scans across the limbs of both the northern and southern bubbles at high latitudes with the Suzaku X-ray satellite. In neither scan do these authors find evidence for a change of temperature between the exterior and interior plasma; rather, the plasma temperature seems constant at ∼(3.0–3.5) × 10⁶ K. For the northern scan, starting in the bright North Polar Spur (NPS) feature (located east in projection from the northern bubble) and scanning inward across the bubble edge, these authors find a drop in emission measure (EM) at the edge from (5–6) × 10⁻² cm⁻⁶ pc outside to (3–4) × 10⁻² cm⁻⁶ pc inside. In the south, the emission measure is approximately constant at (1–2) × 10⁻² cm⁻⁶ pc across the scan. Adopting a total cross-sectional area of ${A}_{\mathrm{FB},{\rm{X}}}=7.0\times {10}^{44}$ cm², we find an average line-of-sight distance through the bubbles of ${l}_{\mathrm{LOS}}\equiv {V}_{\mathrm{FB}}/{A}_{\mathrm{FB},{\rm{X}}}\simeq 3.7$ kpc. With this distance scale and assuming for the moment a uniform density distribution for the bubble plasma, a characteristic ionized hydrogen number density in the bubbles can be estimated as $\langle {n}_{{{\rm{H}}}^{+}}^{\mathrm{FB}}\rangle \sim$ $\sqrt{\mathrm{EM}/{l}_{\mathrm{LOS}}}$ ≃ (1–3) × 10⁻³ cm⁻³. Using the volumetric density estimate, the total plasma mass of the bubbles is then estimated as

$$\begin{eqnarray}&&{M}_{\mathrm{FB}}\simeq 1.4\ {m}_{{\rm{p}}}\ {n}_{{\rm{H}}}\ {V}_{\mathrm{FB}}=2.1\times {10}^{7}\ {M}_{\odot }\left(\displaystyle \frac{\langle {n}_{{{\rm{H}}}^{+}}^{\mathrm{FB}}\rangle }{0.003\;{\mathrm{cm}}^{-3}}\right).\end{eqnarray} \tag{ 2.1 }$$

In reality, we show below that most of the emission measure through the bubbles (and therefore most of their observed thermal X-ray emission) is contributed by gas in a shell whose real density is larger than the volumetric average density estimated above by a factor $\langle {n}_{{{\rm{H}}}^{+}}^{\mathrm{shell}}\rangle /\langle {n}_{{{\rm{H}}}^{+}}^{\mathrm{FB}}\rangle$ ∼ ${(2{\rm{\Delta }}{r}_{\mathrm{shell}}/{l}_{\mathrm{LOS}})}^{-1/2}$, and thickness (conservatively) r_shell ≲ 2 kpc, inside the CD. The fitted X-ray temperature is thus likely to be characteristic of this shell; we expect the volume-filling, shocked plasma between the shell and the reverse shock to be hotter and more rarefied.

Whatever the exact configuration of plasma within the bubbles, consistent with our modeling, they are not strongly overpressured with respect to the surrounding halo (consistent with the absence of obvious X-ray limb brightening in their higher reaches; Kataoka et al. 2013; Tahara et al. 2015). Again working with the characteristic, volumetric-average number density for the moment, the pressure in the bubbles ${P}_{\mathrm{FB}}={n}_{\mathrm{tot}}^{\mathrm{FB}}\ {k}_{B}\ {T}_{\mathrm{FB}}\ \simeq 2.0$ eV cm⁻³ ${n}_{{\rm{H}}}^{\mathrm{FB}}/(0.003$ cm⁻³)T_FB/(3.5 × 10⁶ K) with the total plasma number density ${n}_{\mathrm{tot}}^{\mathrm{FB}}\simeq 2.3\ {n}_{{\rm{H}}}^{\mathrm{FB}}$ (adopting a universal n_He = 0. 1n_H in a totally ionized, electrically neutral plasma and ignoring metals for the moment).

This bubble pressure should be compared with the determination of halo pressures derived from the recent modeling of Miller & Bregman (2014) of O viii emission line data. Adopting halo metallicities of ∼0.3 Z_⊙ as suggested by independent evidence, these authors determine a halo density profile that gives n_e ≃ 2 × 10⁻³ at a halo height of z = 8 kpc and n_e ≃ 6 × 10⁻³ at 4 kpc. Together with their determination of a fairly constant halo gas temperature of ∼2 × 10⁶ K, these results suggest a total atmospheric pressure around the bubbles of ${P}_{\mathrm{ext}}={n}_{\mathrm{tot}}^{\mathrm{halo}}{T}^{\mathrm{halo}}$ ∼ 2.3/1.2 × (6 × 10⁻³ cm⁻³) × 2 × 10⁶ K = 3.2 × 10⁻¹² dyn cm⁻² = 2.0 eV cm⁻³ at their ∼4 kpc half-height.

In the limit that radiative losses can be neglected and the atmospheric density treated as constant, the total enthalpy of a slowly inflated bubble is

$$\begin{eqnarray}&&H\equiv \displaystyle \frac{\gamma }{\gamma -1}{P}_{\mathrm{ext}}V=\dot{E}\ t,\end{eqnarray} \tag{ 2.2 }$$

where $\dot{E}$ is the average mechanical power fed into the bubble over time t and γ is the effective adiabatic index of the bubble fluid. Employing the atmospheric pressure of P_ext = 2.0 eV cm⁻³ motivated above, this suggests H_FB ≃ 6.7 × 10⁵⁵ erg for the bubbles (probably a lower limit as it neglects the z-dependence of the halo density and pressure). Accounting only for core-collapse supernovae (SNe), the minimum mechanical power injected by nuclear star formation can be estimated from the Galactic center (GC) star formation rate SFR_GC as

$$\begin{eqnarray} {\dot{E}}_{\mathrm{GC}\;\mathrm{SF}} & \gtrsim & {\mathrm{SFR}}_{\mathrm{GC}}\times 1\ \mathrm{SN}/(90\;{M}_{\odot })\times {10}^{51}\;\mathrm{erg}/\mathrm{SN}\\ & = & 2.8\times {10}^{40}\;\mathrm{erg}\ {{\rm{s}}}^{-1}\left(\displaystyle \frac{{\mathrm{SFR}}_{\mathrm{GC}}}{0.08\;{M}_{\odot }\;\;{\mathrm{yr}}^{-1}}\right), \end{eqnarray} \tag{ 2.3 }$$

where a single core-collapse SN requires the formation of ∼90 M_⊙ in stars (for a Kroupa 2001 initial mass function if the least massive star to explode has zero-age main-sequence mass of 8 M_⊙), we assume 10⁵¹ erg mechanical energy release per SN, and we have normalized to a nuclear star formation rate of SFR_GC = 0.08 M_⊙ yr⁻¹ (Crocker 2012).

Note that there are other significant star-formation-related sources of mechanical energy, including stellar winds and protostellar outflows, that could, in a full accounting, double this estimate. However, we adopt Equation (2.3) as our fiducial value for the injected power (see Lacki 2014). With this power, a conservative minimum timescale to inflate the bubbles (derived assuming that an unrealistically high 100% of mechanical energy injected in the nucleus is transferred to the bubbles and neglecting cooling and gravitational losses, which we show below are both actually important) is $t={H}_{\mathrm{FB}}/{\dot{E}}_{\mathrm{GC}\ \mathrm{SF}}\simeq 8\times {10}^{7}$ yr. Given the current surface area of the bubbles, the expansion speed where ${\dot{H}}_{\mathrm{FB}}$ saturates the star formation mechanical power is ∼14 km s⁻¹; this is an upper limit on the current expansion velocity of the bubbles' CD.

The total plasma mass in the bubbles of M_FB ∼ 2.1 × 10⁷ M_⊙, if fed by the nuclear outflow at a mass inflow rate, $\dot{M}$, similar to the SFR, requires a minimum timescale to establish of ${M}_{\mathrm{FB}}/\dot{M}=2.6\times {10}^{8}$ yr × (0.08 M_⊙ yr⁻¹$)/\dot{M}$. Normalizing to volumetric average values for plasma temperature and density obtained in Section 2.2, the cooling time of the plasma in the bubbles is ${t}_{\mathrm{cool}}=3/2\ {k}_{B}T/({\rm{\Lambda }}[T]\ {n}_{\mathrm{tot}})\;\simeq 5\times {10}^{8}$ yr × T_FB/(3.5 × 10⁶ K) × (0.007 cm${}^{-3}/{n}_{\mathrm{tot}}^{\mathrm{FB}}$), where⁹  ${\rm{\Lambda }}[T]$ is the cooling function (Dopita et al. 2013; Sutherland et al. 2013) for plasma at temperature T. Thus, even the minimal inflation time is comparable to the radiative cooling time, indicating that an adiabatic treatment is inadequate. Note the confluence of timescales around the formation/maintenance of the bubbles at order ∼few × 10⁸ yr (see Crocker et al. 2014).

We show below (Section 4) that the bubbles approach a steady state such that their radii do not increase and that the freshly heated plasma mass flowing into them is matched by the dropout of cooling plasma, which condenses under a local thermal instability.

In such a steady state, material is advected from the cylindrical acceleration zone in the nucleus, of volume V_GC and radius r₀, by the wind, with speed v_w, in a timescale t_w, implying that a rate of mass loss of

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{GC}}\equiv \mu {\rm{m}}\;{n}_{\mathrm{tot}}^{\mathrm{GC}}\ {V}_{\mathrm{GC}}/{t}_{{\rm{w}}}=\mu {\rm{m}}\;{n}_{\mathrm{tot}}^{\ \mathrm{GC}}\ 2{v}_{{\rm{w}}}\ \pi {r}_{0}^{2}\end{eqnarray} \tag{ 2.4 }$$

(where m is the atomic mass unit and μm is the mean mass per particle in the plasma) is identical to the rate of mass injection into the bubbles, which is, in turn, identical to the rate at which mass is cooling out of the structures,

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{cool}}=\mu {\rm{m}}\;{n}_{\mathrm{tot}}^{\mathrm{FB}}{V}_{\mathrm{FB}}/{t}_{\mathrm{cool}}.\end{eqnarray} \tag{ 2.5 }$$

Hence, we estimate the radius of the nuclear region from which the outflow feeding the bubbles emerges as

$$\begin{eqnarray} {r}_{0} & \simeq & \sqrt{\displaystyle \frac{{V}_{\mathrm{FB}}\;{n}_{\mathrm{tot}}^{\mathrm{FB}}}{2\pi \;{n}_{\mathrm{tot}}^{\mathrm{GC}}\;{v}_{\mathrm{wind}}\;{t}_{\mathrm{cool}}}}\simeq 70{\left(\displaystyle \frac{{v}_{\mathrm{wind}}}{500\;\mathrm{km}\;{{\rm{s}}}^{-1}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{0.004\;{\mathrm{cm}}^{-3}}\right)}^{1/2}{\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{GC}}}{0.1\;{\mathrm{cm}}^{-3}}\right)}^{-1/2}\;\mathrm{pc}, \end{eqnarray} \tag{ 2.6 }$$

where we normalize ${n}_{\mathrm{tot}}^{\mathrm{GC}}$ to a value informed by GC X-ray observations (Muno et al. 2004). This is close to the radius of the central molecular zone (CMZ) molecular torus, 80–100 pc (Molinari et al. 2011). The CMZ is the region of enhanced molecular gas density and star formation activity surrounding the central supermassive black hole, and the good agreement of our estimate of the outflow radius and the radius of this region supports our model that the bubbles are inflated by star formation over this region.

Let ρ_w and v_w be the density and velocity, respectively, of the Galactic center wind. The location of the reverse shock is determined from the condition that the ram pressure of the wind and the internal bubble pressure, p_b, are related by the Rankine–Hugoniot condition for a strong shock¹⁰

$$\begin{eqnarray}&&{p}_{{\rm{b}}}=\displaystyle \frac{2}{\gamma +1}{\rho }_{{\rm{w}}}{v}_{{\rm{w}}}^{2},\end{eqnarray} \tag{ 3.1 }$$

where γ = 5/3. On each side of the nucleus, we approximate the wind as conically expanding with solid angle ${{\rm{\Omega }}}_{{\rm{w}}}$ so that the mass flux is

$$\begin{eqnarray}&&\dot{M}={{\rm{\Omega }}}_{{\rm{w}}}{R}^{2}{\rho }_{{\rm{w}}}{v}_{{\rm{w}}},\end{eqnarray} \tag{ 3.2 }$$

where R is the spherical radius. Hence, the location of the shock is given by

$$\begin{eqnarray}&&{R}_{\mathrm{shk}}={\left(\displaystyle \frac{2}{\gamma +1}\displaystyle \frac{\dot{M}{v}_{{\rm{w}}}}{{p}_{{\rm{b}}}{\rm{\Omega }}}\right)}^{1/2}.\end{eqnarray} \tag{ 3.3 }$$

Normally for a starburst wind, with initial power $\dot{E}$ and mass flux $\dot{M}$, we would assume that the wind velocity, ${v}_{{\rm{w}}}\simeq {(2\dot{E}/\dot{M})}^{1/2}$, is constant. In this case the shock location would be given by

$$\begin{eqnarray}&&{R}_{\mathrm{shk}}=\displaystyle \frac{{2}^{3/4}}{{(\gamma +1)}^{1/2}}\displaystyle \frac{{(\dot{E}\dot{M})}^{1/4}}{{({p}_{{\rm{b}}}{\rm{\Omega }})}^{1/2}}.\end{eqnarray} \tag{ 3.4 }$$

However, the Galactic center wind, with a kinetic power ∼10⁴⁰ ergs s⁻¹, is relatively weak, compared to, say, the ∼10⁴² erg s⁻¹ wind from the dwarf starburst galaxy M82 (Shopbell & Bland-Hawthorn 1998; Cooper et al. 2008), and in the case of the Galactic center, gravity reduces the velocity of the wind on scales ∼kpc.

In order to estimate the effect of gravity on the energy flux and the wind velocity, we begin with the equation for the conservation of energy. In the following _w is the internal energy density of the wind, ϕ is the gravitational potential, and h_w is the specific enthalpy. We have

$$\begin{eqnarray} &&\displaystyle \frac{\partial }{\partial t}\left[\displaystyle \frac{1}{2}{\rho }_{{\rm{w}}}{v}_{{\rm{w}}}^{2}+{\epsilon }_{{\rm{w}}}+{\rho }_{{\rm{w}}}\phi \right]\\ &&\quad +\displaystyle \frac{\partial }{\partial {x}_{i}}\left[\left(\displaystyle \frac{1}{2}{v}_{{\rm{w}}}^{2}+{h}_{{\rm{w}}}+\phi \right){\rho }_{{\rm{w}}}{v}_{{\rm{w}},i}\right]=0. \end{eqnarray} \tag{ 3.5 }$$

Integrating over a volume bounded by streamlines and capped at the top by the reverse shock (S) and below by a surface near the base of the flow (S₀), we have, for a stationary flow, a conserved energy flux,

$$\begin{eqnarray} {F}_{E} & = & {\displaystyle \int }_{S}\left(\displaystyle \frac{1}{2}{v}_{{\rm{w}}}^{2}+{h}_{{\rm{w}}}+\phi \right){\rho }_{{\rm{w}}}{v}_{{\rm{w}},i}{n}_{i}\;{dS}\\ & = & {\displaystyle \int }_{{S}_{0}}\left(\displaystyle \frac{1}{2}{v}_{{\rm{w}},0}^{2}+{h}_{{\rm{w}},0}+{\phi }_{0}\right){\rho }_{{\rm{w}},0}{v}_{{\rm{w}},i,0}{n}_{i,0}\;{{dS}}_{0}. \end{eqnarray} \tag{ 3.6 }$$

Denoting the hydrodynamic part of the energy flux by

$$\begin{eqnarray}&&\dot{E}={\displaystyle \int }_{S}\left(\displaystyle \frac{1}{2}{v}_{{\rm{w}}}^{2}+{h}_{{\rm{w}}}\right){\rho }_{{\rm{w}}}{v}_{{\rm{w}},i}{n}_{i}\;{dS},\end{eqnarray} \tag{ 3.7 }$$

we have that

$$\begin{eqnarray}&&\dot{E}={\dot{E}}_{0}-{\displaystyle \int }_{S}{\rho }_{{\rm{w}}}\phi {v}_{{\rm{w}},i}{n}_{i}\;{dS}+{\displaystyle \int }_{{S}_{0}}{\rho }_{{\rm{w}},0}{\phi }_{0}{v}_{{\rm{w}},0,i}{n}_{i,0}\;{{dS}}_{0}.\end{eqnarray} \tag{ 3.8 }$$

Assuming that the bounding surfaces are approximately equipotential, and taking ${\rm{\Delta }}\phi =\phi -{\phi }_{0}$ as the difference in potential between the two surfaces, we have

$$\begin{eqnarray}&&\dot{E}={\dot{E}}_{0}-\dot{M}{\rm{\Delta }}\phi ={\dot{E}}_{0}\left(1-\displaystyle \frac{\dot{M}}{{\dot{E}}_{0}}{\rm{\Delta }}\phi \right).\end{eqnarray} \tag{ 3.9 }$$

The factor $1-\dot{M}{\rm{\Delta }}\phi /{\dot{E}}_{0}$ determines the reduction in hydrodynamic power resulting from the gravitational field.

We now use Bernoulli's equation to relate the velocity, specific enthalpy, and potential on S with the corresponding

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{v}_{{\rm{w}}}^{2}+{h}_{{\rm{w}}}+\phi =\displaystyle \frac{1}{2}{v}_{{\rm{w}},0}^{2}+{h}_{{\rm{w}},0}+{\phi }_{0}\end{eqnarray} \tag{ 3.10 }$$

and approximate ${\dot{E}}_{0}$ by

$$\begin{eqnarray}&&{\dot{E}}_{0}\simeq \dot{M}\times \left(\displaystyle \frac{1}{2}{v}_{{\rm{w}},0}^{2}+{h}_{{\rm{w}},0}\right).\end{eqnarray} \tag{ 3.11 }$$

This implies that

$$\begin{eqnarray}&&{v}_{{\rm{w}}}^{2}+2{h}_{{\rm{w}}}\simeq 2\displaystyle \frac{{\dot{E}}_{0}}{\dot{M}}-2{\rm{\Delta }}\phi .\end{eqnarray} \tag{ 3.12 }$$

As the flow expands, the enthalpy term diminishes in importance compared to the kinetic term (i.e., ${h}_{{\rm{w}}}\ll {v}_{{\rm{w}}}^{2}/2$). Hence,

$$\begin{eqnarray}&&{v}_{{\rm{w}}}\simeq {\left(\displaystyle \frac{2{\dot{E}}_{0}}{\dot{M}}\right)}^{1/2}{\left(1-\displaystyle \frac{\dot{M}}{{\dot{E}}_{0}}{\rm{\Delta }}\phi \right)}^{1/2}.\end{eqnarray} \tag{ 3.13 }$$

Without the gravitational field, the asymptotic wind speed would be ${(2{\dot{E}}_{0}/\dot{M})}^{1/2}$ ≃ $630{({\dot{E}}_{0}/5\times {10}^{39}\;\mathrm{ergs}\;{{\rm{s}}}^{-1})}^{1/2}{(\dot{M}/0.04\;{M}_{\odot }\;{\mathrm{yr}}^{-1})}^{-1/2}\;\mathrm{km}\;{{\rm{s}}}^{-1}$. The factor ${\left(1-\dot{M}{\rm{\Delta }}\phi /{\dot{E}}_{0}\right)}^{1/2}$ determines the reduction in wind speed resulting from the gravitational potential.

The factor $\left(1-\dot{M}{\rm{\Delta }}\phi /{\dot{E}}_{0}\right)$, appearing in Equation (3.13), depends on the location of the reverse shock through the dependence of the potential difference ${\rm{\Delta }}\phi$ on cylindrical coordinates, based on the Galactic center. We now describe how we determine the approximate location of the reverse shock and the hydrodynamic energy flux there.

Using Equation (3.1), relating the bubble pressure (p_b) and the wind ram pressure, and (3.2) for the mass flux, together with Equation (3.13) for the wind velocity, we have for the radial location of a spherical shock in a conically expanding wind:

$$\begin{eqnarray}&&{r}_{\mathrm{shk}}={\left(\displaystyle \frac{2}{\gamma +1}\right)}^{1/2}\;\displaystyle \frac{{(2{\dot{E}}_{0}\dot{M})}^{1/4}}{{\left({p}_{{\rm{b}}}{\rm{\Omega }}\right)}^{1/2}}{\left(1-\displaystyle \frac{\dot{M}}{{\dot{E}}_{0}}{\rm{\Delta }}\phi \right)}^{1/4}.\end{eqnarray} \tag{ 3.14 }$$

For the potential of the Galaxy, we use the analytic forms and parameters provided by Breitschwerdt et al. (1991). For the combined disk and bulge potential difference from the center of the Galaxy, expressed in (km s⁻¹)², as a function of cylindrical coordinates r and z, both in kpc, this gives

$$\begin{eqnarray} {\rm{\Delta }}\phi [r,z] & \simeq & 1.60\times {10}^{5}\;-\displaystyle \frac{8.82\times {10}^{4}}{\sqrt{{r}^{2}+{z}^{2}+0.245}}\\ & & -\displaystyle \frac{1.10\times {10}^{6}}{\sqrt{{r}^{2}+{\left(\sqrt{{z}^{2}+0.270}+7.26\right)}^{2}}}\\ & & +\displaystyle \frac{5.81\times {10}^{5}}{13+\sqrt{{r}^{2}+{z}^{2}}}+4.47\\ & & \times {10}^{4}\mathrm{ln}\left[13+\sqrt{{r}^{2}+{z}^{2}}\right]\quad {(\mathrm{km}\;{{\rm{s}}}^{-1})}^{2}. \end{eqnarray} \tag{ 3.15 }$$

Notwithstanding the expression for ${\rm{\Delta }}\phi$ being in terms of cylindrical coordinates, the contours of ${\rm{\Delta }}\phi$ are approximately spherical in the region of interest (r_shk ≲ 1.5 kpc). Hence, in solving Equation (3.14) we have approximated ${\rm{\Delta }}\phi$ by its value on the z-axis, i.e., ${\rm{\Delta }}\phi [0,{r}_{\mathrm{shk}}]$.

In Figure 2 we show, as a function of the mechanical power injected at the nucleus into one hemisphere (for the parameters specified in the caption), our numerical solution for (i) the distance to the shock in kpc (with and without gravitational deceleration of the flow) and (ii) the relative mechanical power available at the shock (given gravitational deceleration). It is evident that for the chosen parameter values, which are apposite to the nuclear outflow that feeds the Fermi bubbles, the flow can lose 30%–50% of its mechanical power to gravity before reaching the shock.

Zoom In Zoom Out Reset image size

As a first approximation to modeling the Fermi bubbles, we treat them as two hemispheres of a single spherical bubble driven by a spherically symmetric wind propagating into a constant-density, constant-pressure medium. In the observed Fermi bubbles the morphology close to the Galactic plane is probably constrained by the larger density and pressure there. In comparing the radius of this model with the extent of the Fermi bubbles, we calculate an equivalent spherical radius ${r}_{\mathrm{FB}}={\left[3/(4\pi ){V}_{\mathrm{FB}}\right]}^{1/3}\simeq 4.1\;\mathrm{kpc}$, where V_FB = 285 kpc³ is the estimated volume. Our model generalizes the usual treatment of a wind-driven bubble (e.g., Weaver et al. 1977), in that the bubble is not assumed to be strongly overpressured with respect to the ambient medium and is also subject to significant optically thin radiative cooling. The model in this section does not incorporate a gravitational potential (so that there is no internal or external pressure stratification) and, in common with other bubble models, neglects the small but finite fraction of the total bubble volume in the region upstream of the internal reverse shock.

Figure 3 lays out the main features of such a bubble: there is a reverse shock interior to the bubble, where the decreasing ram pressure of the spherically expanding, supersonic nuclear wind is balanced by the internal pressure of the bubble as described in Section 3, and a region of shock-heated bubble material downstream of the reverse shock, in which the internal energy dominates the kinetic energy. This region is bounded by a CD, external to which is a shell of shocked, ambient ISM between the forward shock, which is expanding into the external atmosphere, and CD. We approximate this region as one of constant spatial pressure and velocity and neglect cooling there. We also in this section approximate the region between CD and reverse shock as constant in density. The bubble is fed by mechanical energy injected at a rate $\dot{E}$ and mass at a rate $\dot{M}$. We assume a constant γ-law equation of state, with the pressure, p, and internal energy density, , related by $p=(\gamma -1)\epsilon$ and γ = 5/3. It is important to note that the growth of the total mass of the bubble M_b does not satisfy ${\dot{M}}_{{\rm{b}}}=\dot{M}$ because of the phenomenon of mass dropout driven by cooling (Crocker et al. 2014). We estimate the rate of mass dropout by equating the post-reverse-shock enthalpy flux of that mass with the cooling inside the bubble.

Zoom In Zoom Out Reset image size

The equations governing the pressure, radius, mass, and density of the bubble are derived in the Appendix, using energy, mass, and momentum conservation integrated over the various regions identified above. We initiate the evolution of the bubble at t = 10⁵ yr with initial values defined by an adiabatic bubble (see the Appendix for details).

Cooling sets up local thermal instability of the bubble plasma, leading it to collapse into overdense condensations that, having lost buoyancy (in a real situation that includes a gravitational potential), fall out of the bubble (see Crocker et al. 2014 and references therein).

As the bubble's pressure decreases as a result of the combined effect of the work done by expansion and the radiant loss of energy, it can reach a steady state where its pressure equilibrates with the surrounding atmosphere, the energy input is balanced by thermal radiation, and the mass flux is balanced by mass dropout. We show below that, for parameters appropriate to the Fermi bubbles, the system is approaching this limit.

On the basis of the numbers presented in Section 2.3, we choose as representative values for the total mass and mechanical power injected at the nucleus $({\dot{E}}_{0},\dot{M})$ ∼ $({10}^{40}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$, 0.08 M_⊙ yr⁻¹) (so that ${\dot{E}}_{1/\mathrm{2,0}}\equiv 1/2\times {\dot{E}}_{0}$ and ${\dot{M}}_{1/2}\equiv 1/2\times \dot{M}$ is injected into the flow into each hemisphere); the exact parameter range allowable for these quantities is delimited in Section 5. One must remember, however, that in Section 3.3 we showed that, for parameter values in this range, the flow into each bubble can lose ∼30%–50% of its mechanical power to gravity before reaching the shocks. For purposes of comparing our current modeling—which neglects gravity—with our previous (and subsequent) results, we now therefore choose an injected mechanical power that is reduced by 40% with respect to the fiducial value nominated above, namely, 6 × 10³⁹ erg s⁻¹. Our purpose here is to show that, for parameters relevant to the Fermi bubbles, the expansion of a radiative bubble (into an atmosphere of pressure similar to that of the Galactic halo) quickly slows down (consistent with the lack of X-ray brightening at the edges of the Fermi bubbles) and, moreover, there is an asymptotic radius for such a bubble that is not much larger than the equivalent radius of the observed Fermi bubbles.

As justified above, another relevant parameter is an ambient pressure of 2.0 eV cm⁻³ (=3.2 × 10⁻¹² dyn cm⁻²). Note that we assume that the plasma has solar metallicity but verify that, given that the metallicity affects the cooling, assuming twice solar metallicity does not change our results significantly (see Figure 4).

Zoom In Zoom Out Reset image size

Figure 4 shows the evolution of the bubble radius (i.e., CD) for input parameters $\dot{E}=6.0\times {10}^{39}$ erg s⁻¹ and $\dot{M}=0.08\;{M}_{\odot }$ yr⁻¹ (i.e., 1.2 × 10¹⁵ erg g⁻¹ as the energy per mass content of the outflow), as well as the locus of the forward shock. The radius of a highly overpressured bubble and the equivalent spherical radius of the Fermi bubbles are also shown for comparison. Clearly, for the fiducial parameters the bubble radius asymptotes to approximately 5 kpc after about 1 Gyr, and from about 3 × 10⁸ yr onward, the bubble radius is not far from the spherical equivalent radius ≃4.1 kpc of the observed bubbles. For comparison, for fixed $\dot{M}=0.08\;{M}_{\odot }$ yr⁻¹, cases of $\dot{E}=1.2\times {10}^{40}$ erg s⁻¹ (2.4 × 10¹⁵ erg g⁻¹) and $\dot{E}=3.0\times {10}^{39}$ erg s⁻¹ (6.0 × 10¹⁴ erg g⁻¹) are also shown in Figure 4. The former does not asymptote to a final radius of ∼10 kpc until 10 Gyr (with volume ∼15 × V_FB), and it spends only a brief time in the vicinity of the currently observed size, tending to suggest a return to the "Why now?" problem encountered by explosive scenarios. The latter never reaches the size (or minimum temperature) of the observed Fermi bubbles.

The expansion of the bubble decelerates from approximately 2 × 10⁶ yr on, and this is clearly shown in Figure 5. From about 10⁸ yr on, the bubble expansion velocity decelerates to less than 10 km s⁻¹ as it approaches pressure equilibrium with the environment. The approach to pressure equilibrium is shown in Figure 6. At the same time, the forward shock becomes a sound wave propagating into the background medium at a sound speed of ≃280 km s⁻¹ (again see Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the bubble's growth saturates, its total mass and energy content are fixed, with the post-reverse-shock enthalpy flux balanced by the radiative luminosity. Figure 7 shows that the mass dropout saturates at the injection rate of 0.08 M_⊙ yr⁻¹ and that, consistent with this, the net mass injection into the bubbles saturates at zero, maintaining the total mass interior to the bubbles at a constant value.

Zoom In Zoom Out Reset image size

For comparison with other sections, we can estimate the final, asymptotic volume, V_f, in this circumstance, with the input hydrodynamic energy flux, $\dot{E}\equiv 5/2\dot{M}/(\mu {\rm{m}}){k}_{B}T$, arriving at and processed through the shocks, balanced by radiative cooling and with the pressure in the downstream bubble region in equilibrium with the halo atmosphere. Here the post-reverse-shock bubble temperature is $T=3/16\;\mu {\rm{m}}/{k}_{B}\;{v}_{{\rm{w}}}^{2}\simeq 3.5\times {10}^{6}{({v}_{{\rm{w}}}/500\;\mathrm{km}\ {{\rm{s}}}^{-1})}^{2}\;{\rm{K}}$; the total post-shock number density in the bubble is ${n}_{\mathrm{tot}}=4\dot{M}/(\mu {\rm{m}}{\rm{\Omega }}{R}_{\mathrm{sh}}^{2}{v}_{{\rm{w}}})$ (see Section 3 for the definitions of these symbols). Balancing the energy injected into the bubble with thermal cooling and normalizing to fiducial values for temperature and number density (such that pressure equilibrium with the halo is satisfied), we then have that the final volume of the bubble is

$$\begin{eqnarray} {V}_{{\rm{f}}} & \simeq & \displaystyle \frac{\dot{E}}{{n}_{\mathrm{tot}}^{2}{\rm{\Lambda }}[T]}=\displaystyle \frac{5/2\ \dot{M}/(\mu {\rm{m}})\ {k}_{B}{T}_{\mathrm{FB}}}{{n}_{\mathrm{tot}}^{2}{\rm{\Lambda }}[T]}\\ & = & 2.0\times {10}^{67}\;{\mathrm{cm}}^{3}\;\left(\displaystyle \frac{{T}_{\mathrm{FB}}}{3.5\times {10}^{6}\;{\rm{K}}}\right)\left(\displaystyle \frac{\dot{M}}{0.08\;{M}_{\odot }\;\;{\mathrm{yr}}^{-1}}\right)\;\\ & & \times {\left(\displaystyle \frac{{n}_{\mathrm{tot}}}{0.0066\ {\mathrm{cm}}^{-3}}\right)}^{-2}{\left(\displaystyle \frac{{\rm{\Lambda }}[T]}{{10}^{-23.15}\ \mathrm{erg}\ {\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}}\right)}^{-1}, \end{eqnarray} \tag{ 4.1 }$$

where we normalize to an $\dot{M}$ favored by the modeling presented in Section 5. This is ∼2.3× the observationally estimated volume of the Fermi bubbles, V_FB = 8.2 × 10⁶⁶ cm⁻³, consistent with our expectation that the bubbles are still slowly expanding. Note that the estimated volume would be smaller if the increase in density of cooling material in the shell region (R3) were taken into account.

We now explain our procedure for scanning over the steady mass and mechanical energy flux, $\dot{M}$ and ${\dot{E}}_{0}$, respectively, into the central acceleration region. As shown above, these parameters determine the temperature reached within the central acceleration region following thermalization of the mechanical energy. The high pressure in this region accelerates the plasma into a flow that is collimated into a biconical wind transverse to the Galactic disk by the surrounding dense molecular torus of the CMZ. In this expanding flow, the plasma cools adiabatically until it reaches the reverse shock in each bubble, where it is shock (re)heated. For every $\{{\dot{E}}_{0},\dot{M}\}$ in the parameter space, the asymptotic velocity of the wind, v_asym, is determined. Note that v_asym is reached soon after the flow escapes from the central acceleration zone; subsequent to reaching v_asym, the adiabatically cooling flow decelerates in response to the gravitational potential, ϕ (Breitschwerdt et al. 1991), as described above. Given the high (i.e., supersonic) upstream flow speeds in which we are interested, we treat the flow evolution in the ballistic limit, where we ignore small corrections due to the temperature evolution of the upstream gas.

In our modeling, for every sampled $\{{\dot{E}}_{0},\dot{M}\}$ point, we determine (i) the height of the reverse shocks R_sh and (ii) the temperature of the shocked material, T_sh. Given that, observationally, the bubbles do not expand latitudinally very much at heights greater than the expected position of each shock, we neglect further possible adiabatic plasma temperature losses with height above the shocks. We also assume a uniform plasma temperature in R2 given by that for a strong shock with upstream flow velocity v_w and γ = 5/3. Hence, ${T}_{\mathrm{FB}}={T}_{\mathrm{sh}}=3/16\ \mu {\rm{m}}\ {v}_{{\rm{w}}}^{2}/{k}_{B}$.

The shocks deliver a fraction, _CR, of their total mechanical power (accounting for gravitational losses), $\dot{E}$, into freshly accelerated cosmic rays: ${L}_{\mathrm{CR}}\equiv {\epsilon }_{\mathrm{CR}}\times \dot{E}$; the spectral index of the cosmic rays' power-law distribution is controlled by the Mach number of the shocks. We shall see subsequently that, a posteriori, we are warranted in assuming strong shocks: the hard, nonthermal spectra of both the γ-ray emission from the Fermi bubbles and the microwave haze require shock Mach numbers of M ≳ 4 (see Section 6). We also find (see Section 6.2.2) that given that nonthermal ISM phases dominate the dynamics of the compression shell R3, its size is dependent on the efficiency with which the shocks accelerate cosmic rays, _CR (see Equation (6.9)). We have explored a range of values for this parameter. Also important are the pressure of the atmosphere into which the bubbles are expanding and the effective temperature of the compression shell. For the former our fiducial value is P_halo = 2 eV cm⁻³ = 3.2 × 10⁻¹² erg cm⁻² adopted from Section 2.2, and, for the latter, following Section 2.2, we fix ${T}_{\mathrm{shell}}=3.5\times {10}^{6}$ K (with T_FB somewhat larger in general).

With the downstream temperature and shock height corresponding to every $\{{\dot{E}}_{0},\dot{M}\}$ point determined, we can determine a number of subsequent quantities, including the Mach number of the shock, the bolometric thermal luminosity of the bubbles, and the thickness of cooling material in the compression shell leading up to the CD from the inside (region R3 in Figure 1). By confronting these modeled values with physical or observational constraints, we can exclude much of the $\{{\dot{E}}_{0},\dot{M}\}$ parameter space; see Figure 8. In this figure we show the case that the downstream, shocked plasma in region R2 has reached pressure equilibrium with an exterior atmosphere of P_halo = 2 eV cm⁻³, the fraction of mechanical power available at the shocks going into cosmic rays is ${\epsilon }_{\mathrm{CR}}=0.15$, and the compression shell at R3 has an effective temperature of T_shell = 3.5 × 10⁶ K. Most of the potential parameter space for these parameters can be excluded by only three considerations.

• 1.  
  The compression shell should not be too thick; the requirement of a flat, projected surface brightness in γ-rays leads to a geometric restriction that the shell thickness, d_shell, satisfies d_shell < 0.456r_B ∼ 1.4 kpc (see Section 6.2.5). This excludes the green region in Figure 8. Note that the shell thickness is determined by the generalization of Equation (6.10) below.
• 2.  
  The time for the gas to flow from the shocked interior (R2) into the compression shell (R3) (at speed given by Equation (6.4) below) should be shorter than the radiative cooling time in R2 (otherwise, cooling takes place in R2, not in R3, at variance with the assumptions of the model).
• 3.  
  The total cooling luminosity, ${L}_{\mathrm{cool}}^{\mathrm{shell}}$, of the compression shell, R3, should not exceed the upper limit on the bolometric luminosity of the X-ray bulge, which, conservatively, is 10⁴⁰ erg s⁻¹ (Snowden et al. 1997; Almy et al. 2000).

Zoom In Zoom Out Reset image size

Note that while, in principle, the allowed parameter space seems to extend off the plot toward the lower left, there are two considerations that exclude ${\dot{E}}_{0}\lesssim 3\times {10}^{39}$ erg s⁻¹, at least in the context of the current work: (i) within a hadronic model, the minimum (thick target) cosmic-ray luminosity necessary to sustain the bubbles' γ-ray luminosity is ∼3 × 10³⁸ erg s⁻¹ (Crocker & Aharonian 2011); a mechanical power at least a factor of ∼3 (for high cosmic-ray acceleration efficiency) larger than this has to be be delivered at the shocks in order to sustain their emission, and one must then factor in gravitational losses (see Figure 2); (ii) we show below (see Figure 12) that, in order to reproduce the giant outflows' synchrotron luminosity, an unreasonably large fraction of the total cosmic-ray luminosity of the shocks has to go into cosmic-ray electrons for ${\dot{E}}_{0}\lesssim 3\times {10}^{39}$ erg s⁻¹.

We show below that, within the restricted parameter space we have mapped out, our model well reproduces the γ-ray, microwave, and radio phenomenology of the bubbles. First, however, we remark on three interesting findings.

First, the $\dot{M}/{\dot{E}}_{0}$ ratio indicated by this restricted parameter space is again entirely consistent with—and, indeed, indicative of—the bubbles being driven by nuclear star formation. In Figure 8 the overplotted black curves delineate expected star-formation-driven mass flux versus star-formation-driven mechanical energy input for mass loading rates, as marked, of β = 2 and 4 that bracket our allowable parameter space (scaling results of Strickland & Heckman 2009). Here mass loading is defined with respect to the mass injection rate from SN and stellar winds, ${\dot{M}}_{\mathrm{SN}+\mathrm{SW}}$, that accompanies a given rate of star formation, SFR, not with respect to this star formation rate itself: ${\dot{M}}_{\mathrm{SN}+\mathrm{SW}}={\dot{M}}_{\mathrm{SN}+\mathrm{SW}}[\mathrm{SFR}]$ with ${\dot{M}}_{\mathrm{tot}}\equiv \beta \times {\dot{M}}_{\mathrm{SN}+\mathrm{SW}}$ (Strickland & Heckman 2009). From modeling of the diffuse, hard X-ray emission emanating from the starburst exemplar M82, Strickland & Heckman (2009) find—for centralized mass loading only—mass loading rates in the range 1 < β < 2.8 with a practical upper limit at β ∼ 10 for a very hot, (30–70) × 10⁶ K, wind. The near coincidence of our allowable parameter space with a scale for the β parameter derived from modeling other star-formation-driven outflows (Strickland & Heckman 2009, and references therein) supports the notion that the outflow is star formation driven; a priori there is no reason for an AGN-driven outflow to be anywhere close to this β range.

Note that there is, in general, the possibility for an outflow to also be subjected to distributed (in addition to centralized) mass loading. The population of H i clouds entrained into the nuclear outflow detected by McClure-Griffiths et al. (2013) and the high-speed, warm ionized material found by Fox et al. (2014) together constitute plausible evidence that exactly such distributed mass loading is taking place in the nuclear outflow, presumably as it scrapes along the inside of the CMZ molecular gas torus. The fact that our allowable parameter space extends into regions corresponding to very slightly higher values of β than determined by Strickland & Heckman (2009) for M82 may reflect the occurrence of both centralized and distributed mass loading in generating the final effective β. Note that the pressure in the CMZ H₂ is very close to that in the central X-ray-emitting plasma; the molecular gas is thus able to collimate this plasma and cause it to escape down the steepest density gradient, into the Galactic halo.¹¹

A second interesting finding comes from consideration of the dashed red curves in Figure 8. These show where the final volumes (obtained from the conditions of pressure equilibrium and balance between cooling radiation and incoming enthalpy flux according to Equation (4.1)) of spherical, radiative bubbles fed with energy and mass at the nominated rates are, respectively, equal to the current estimated volume of the observed Fermi bubbles and 15× the current volume. Left and upward of the V_FB curve, radiative bubbles never reach the current size of the observed Fermi bubbles, therefore also ruling out this region of parameter space. With increasing distance right and downward of the 15 V_FB curve, there is an increasingly severe temporal fine-tuning problem as radiative bubbles spend only a brief time over their entire evolution close to the current observed size. We emphasize, then, that for the region of $\{{\dot{E}}_{0},\dot{M}\}$ parameter space delimited by the independent considerations set out above, the Fermi bubbles are naturally found to have a size close to that observed. In other words, there is no "Why now?" problem in this model.

A final interesting finding is that allowable values of P_halo and _CR = 0.15 are quite tightly constrained within our model: we are not able to reproduce the observed γ-ray luminosity of the bubbles (see Section 6 below) for P_halo ≲ 1.5 eV cm⁻³; for _CR ≳ 0.3 we find that the parameter space allowed by the three considerations listed above closes off completely, and for _CR ≲ 0.05 we cannot simultaneously reproduce both the observed γ-ray and radio+microwave luminosities of the giant outflows (see below). This allowable range for _CR is consistent with independent determinations of the fraction of shock mechanical power that goes into such particles in other systems (e.g., Zirakashvili & Völk 2006, on the nuclear outflow from local starburst NGC 253).

We lastly note that, in the allowed parameter space, cooling in the shell does not saturate the mechanical energy injected in the nucleus; this means that the CD must be expanding and thereby doing work on its surroundings at a rate $P\dot{V}$ that saturates the residual mechanical energy. The bubbles, therefore, have not yet reached their asymptotic, final volume. In fact, we find that the CD must be expanding at about 2 km s⁻¹ to saturate the injected power not lost to cooling radiation in the shell. Comparing against the results from Section 4, the spherical radiative bubble is expanding this quickly at an age of 5 × 10⁸ yr when, self-consistently, it has reached a volume equivalent to the observed Fermi bubbles (see Figure 5); thus, the model hangs together in this parameter region.

In the test particle limit (where the dynamics are completely dominated by thermal particles), the spectral index of the nonthermal particle population accelerated at a shock is given by (Bell 1978)

$$\begin{eqnarray}&&{\gamma }_{\mathrm{CR}}=\displaystyle \frac{R+2}{R-1}=\displaystyle \frac{(\gamma +1){M}^{2}}{(\gamma -1){M}^{2}+2},\end{eqnarray} \tag{ 6.1 }$$

where R is the shock compression ratio. Observations of γ-ray and microwave emission reviewed below indicate populations of nonthermal protons and electrons with power-law distributions ${dN}/{dE}\propto {E}^{-\gamma }$ with spectral indices of γ < 2.3 requiring Mach numbers M ≳ 4. Consistent with this, our modeling favors upstream flows at the giant shocks with Mach numbers in the range 6–9. This will generate nonthermal particle populations with spectral indices of ∼2.1, close to the central values required to explain both the bubbles' γ-ray and the haze's hard microwave emission.

Energy scales for individual, nonthermal particles of interest to us can be reckoned as follows. Adopting a fiducial magnetic field of 7 μG downstream of the shock (Section 6.3.3), the characteristic energy of electrons synchrotron-radiating at ν = 23 GHz is E_c = 18 GeV(B/7 μG)^−1/2(ν/23 GHz)^1/2. For IC emission, γ-rays at E_γ = 1 GeV are upscattered by primary electrons of energy ${E}_{e}\simeq 640$ GeV(E_γ/1 GeV)^1/2(E_γ,0/6.3 × 10⁻⁴ eV)^−1/2 adopting an energy for the photon background characteristic of the CMB. Finally, proton primaries of hadronic γ-rays in the observed range of 1–100 GeV have energies covering the approximate range 10–1000 GeV.

These characteristic energies should be compared with upper limits to the energies that the giant shocks might accelerate nonthermal particles. The Hillas criterion that a particle's acceleration time be at least equal to the maximal diffusive escape time from an acceleration region, of fiducial size r_acc ∼ 100 pc,¹² sets a maximum energy limit of E_max ≃ (10¹⁵ eV)(r_acc/100 pc)(B/7 μG)(v/500 km s⁻¹). For synchrotron-loss-limited acceleration of electrons over a timescale t_acc ≡ 4D/v², where D is the diffusion coefficient near the shock and v is a characteristic velocity, in the most optimistic case of Bohm diffusion D_Bohm = 1/3 r_gc, with r_g the particle gyroradius, we find a maximum upper limit to the accelerated electron energy of E_e^max ≃ (4 × 10¹³ eV)(B/7 μG)^−1/2(v/500 km s⁻¹). In this case, the shocks can easily accelerate electrons up to the requisite energies to produce the observed microwave emission and even the IC γ-rays by upscattering the CMB (though whether such high-energy electrons can reach the top of the bubbles is another matter).

In our modeling we have explored different values for the fraction, _CR, of shock mechanical power that is injected into freshly accelerated, nonthermal particles, L_CR, finding that ranges of ∼5%–20% are acceptable. For the magnetic field amplitudes explored below, we find that electrons are required to take up a fraction 5%–30% of L_CR in order to explain the radio and microwave emission from the bubbles, with the remainder going into ions, predominantly protons.

6.2.1. Observations

6.2.2. Expected Hadronic Gamma-Ray Emission

6.2.3. Hadronic Gamma-Ray Luminosity

Adapting the calculations presented by Crocker et al. (2014), we can calculate the hadronic luminosity of the bubbles. Our essential argument is that, given the expansion timescale, the evolution of the plasma in the bubbles cannot be treated as adiabatic, but rather the implication of plasma cooling must be considered. This cooling leads to local thermal instability and the formation of overdense, cool condensations that, having lost buoyancy, fall down through the bubbles under gravity. (In the steady-state limit the dropout of these condensations completely saturates the injection rate of fresh plasma mass.) Globally there is a slow flow in the bubbles of material downstream of the reverse shock and toward the CD. As the CD is approached, this flow must be decelerated with a consequent increase in plasma density. Given the n² dependence of the cooling, it is in this overdense shell region (R3 in Figure 1) that thermal instability will preferentially form the condensations. Thus, there will be a flow of gas at a rate (saturating $\dot{M}$ in the steady-state limit) of shocked plasma from the interior volume of the bubbles behind the reverse shocks, R2, into this shell, R3, with subsequent thermal collapse, and dropout under gravity. We review the observational evidence for formation of cool gas condensations near the edges of the bubbles in Section 7.2.

In our model, in the interior region, R2, the gas evolves adiabatically, whereas in the denser shell, R3, it evolves radiatively. However, a crucial point is that cosmic rays and magnetic fields evolve adiabatically in both shell and volume gas. Adiabatic compression of these nonthermal ISM components thus raises their relative pressure within the shell. In fact, compression of the shell by the plasma in R2 is arrested when the adiabatic energy gain of magnetic fields and cosmic rays results in their reaching pressure equilibrium with the interior plasma. The size of the shell is dually controlled by this and the requirement that the cooling radiation taking place within it (strictly, in the steady-state limit) balances the enthalpy flux represented by the flow of hot plasma from the hot interior.

As we now show, we can use this pressure equilibrium between the adiabatically compressed cosmic rays and magnetic fields in R3 and the thermal plasma in R2 and the size of the shell determined by the enthalpy flux argument to determine the hadronic γ-ray luminosity of the bubbles. This is dominated by collisions in the shell between the adiabatically compressed cosmic rays and the compressed, cooling gas, which can be calculated as shown below. The Fermi-LAT band hadronic ("pp") luminosity of the shell is

$$\begin{eqnarray} {L}_{\gamma }^{\mathrm{pp}} & \simeq & \displaystyle \frac{3}{2}\times \displaystyle \frac{1}{3}\times {f}_{\mathrm{bol}}\ {u}_{{\rm{p}}}^{\mathrm{shell}}\ {V}_{\mathrm{shell}}\ {n}_{{\rm{H}}}^{\mathrm{shell}}\ {\sigma }_{\mathrm{pp}}\ {\kappa }_{\mathrm{pp}}\ c\\ & = & \displaystyle \frac{{f}_{\mathrm{bol}}}{2}\ {u}_{{\rm{p}}}^{\mathrm{shell}}\ \displaystyle \frac{{M}_{\mathrm{shell}}}{\mu {\rm{m}}}\ \displaystyle \frac{{n}_{{\rm{H}}}^{\mathrm{shell}}}{{n}_{\mathrm{tot}}^{\mathrm{shell}}}\ {\sigma }_{\mathrm{pp}}\ {\kappa }_{\mathrm{pp}}\ c. \end{eqnarray} \tag{ 6.2 }$$

Here, in the first line, the first 3/2 prefactor (approximately) corrects for the presence of heavy ions among beam and target nuclei (Mori 1997), 1/3 comes from the relative multiplicity of π⁰ among all daughter pions, f_bol is the fraction of the bolometric luminosity emitted in the relevant γ-ray band, ${u}_{{\rm{p}}}^{\mathrm{shell}}$ is the energy density of relativistic protons in the shell, σ_pp ≃ 4 × 10⁻²⁶ cm² is the total hadronic cross section, and κ_pp ≃ 0.5 is the hadronic inelasticity. In the second line ${n}_{{\rm{H}}}^{\mathrm{shell}}/{n}_{\mathrm{tot}}^{\mathrm{shell}}=0.44578$ is the H⁺ to total plasma number density ratio in the shell at solar metallicity, m is the atomic mass unit, and μ = 0.6039 is the mean mass of the particles in a solar-metallicity plasma.

We now assume pressure balance between the relativistic shell components ("p" for protons and heavier ions and "B" for the magnetic field), ${P}_{{\rm{p}}}^{\mathrm{shell}}+{P}_{B}^{\mathrm{shell}}$, and the thermal plasma in the interior, P_FB, as justified above (and also adopt cosmic-ray/magnetic field equipartition in the shell as justified observationally by Carretti et al. 2013, ${P}_{{\rm{p}}}^{\mathrm{shell}}\simeq {P}_{B}^{\mathrm{shell}}$), so that the pressures represented by these various components are related by

$$\begin{eqnarray}&&{P}_{{\rm{p}}}^{\mathrm{shell}}\equiv \displaystyle \frac{1}{3}{u}_{{\rm{p}}}^{\mathrm{shell}}=\displaystyle \frac{1}{3}{u}_{B}^{\mathrm{shell}}=\displaystyle \frac{1}{2}{P}_{\mathrm{FB}}.\end{eqnarray} \tag{ 6.3 }$$

We now calculate M_shell as follows: mass conservation implies that the flux of mass from interior volume, of number density ${n}_{\mathrm{tot}}^{\mathrm{FB}}$, into the shell is equal (in the steady-state limit) to the mass flux into the bubbles, $\dot{M}$. Thus, the effective speed, v, with which plasma is carried across the surface, of area A_shell, separating the interior region R2 from the shell R3 is

$$\begin{eqnarray} v & \equiv & \displaystyle \frac{\dot{M}}{\mu {\rm{m}}\ {A}_{\mathrm{shell}}\ {n}_{\mathrm{tot}}^{\mathrm{FB}}}\simeq 13\;\;\mathrm{km}\;{{\rm{s}}}^{-1}\ \left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{0.002\;{\mathrm{cm}}^{-3}}\right)}^{-1}\ {\left(\displaystyle \frac{{A}_{\mathrm{shell}}}{2.5\times {10}^{45}\;{\mathrm{cm}}^{2}}\right)}^{-1}. \end{eqnarray} \tag{ 6.4 }$$

The enthalpy flux represented by this flow is dissipated in cooling radiation from the shell; thus,

$$\begin{eqnarray}&&\displaystyle \frac{5}{2}\ {n}_{\mathrm{tot}}^{\mathrm{FB}}\ {k}_{B}\ {T}_{\mathrm{FB}}\ v\ {A}_{\mathrm{shell}}={({n}_{\mathrm{tot}}^{\mathrm{shell}})}^{2}\ {\rm{\Lambda }}[{T}_{\mathrm{shell}}]\ {V}_{\mathrm{shell}},\end{eqnarray} \tag{ 6.5 }$$

where ${\rm{\Lambda }}[{T}_{\mathrm{shell}}]$ is the cooling function at the shell temperature and metallicity (the latter assumed to be solar as for the hot plasma) and, geometrically, we assume a thin shell. As for the shell temperature, as we have emphasized, the X-ray spectrum across the edge of the bubbles, encompassing the region of the shell, has been successfully fit with a single-temperature, collisional ionization equilibrium plasma model of temperature ∼3.5 × 10⁶ K, and we therefore adopt T_shell = 3.5 × 10⁶ K (so that the temperature in the interior volume is necessarily >3.5 × 10⁶ K). In reality, the gas distribution in the shell is more complicated as the gas cascades from the high injection temperature, where it starts to cool, down to optical emission line temperatures. Such a multitemperature distribution is inferred, e.g., from the coincident X-ray and Hα emission observed from filaments around the bright central galaxies in clusters (e.g., Fabian et al. 2011) and is common to cooling flows. In any case, here we ignore the complication of the multitemperature nature of the shell gas and assume that the effective temperature apposite to the calculation of shell cooling losses is precisely that found from fitting to the observed X-ray radiation for a single-temperature model.

Rearranging and combining Equations (6.4) and (6.5), we find

$$\begin{eqnarray} {V}_{\mathrm{shell}} & = & \displaystyle \frac{5}{2}\displaystyle \frac{{k}_{B}\ {T}_{\mathrm{FB}}\ \dot{M}}{{({n}_{\mathrm{tot}}^{\mathrm{shell}})}^{2}\ {\rm{\Lambda }}[{T}_{\mathrm{shell}}]\ \mu {\rm{m}}}\\ & = & \displaystyle \frac{5}{2}\displaystyle \frac{{k}_{B}\ {T}_{\mathrm{FB}}\dot{M}}{{({n}_{\mathrm{tot}}^{\mathrm{FB}})}^{2}\ {\rm{\Lambda }}[{T}_{\mathrm{shell}}]\ \mu {\rm{m}}}{\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{{n}_{\mathrm{tot}}^{\mathrm{shell}}}\right)}^{2}. \end{eqnarray} \tag{ 6.6 }$$

Following Equation (6.3), specifiying pressure balance between the adiabatically compressed shell magnetic field and cosmic rays and the thermal gas in the interior, and given adiabatic evolution of the cosmic rays,

$$\begin{eqnarray}&&{u}_{{\rm{p}}}^{\mathrm{shell}}={\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{shell}}}{{n}_{\mathrm{tot}}^{\mathrm{FB}}}\right)}^{4/3}{u}_{{\rm{p}}}^{\mathrm{FB}},\end{eqnarray} \tag{ 6.7 }$$

then for ${u}_{\mathrm{therm}}^{\mathrm{FB}}$ as the internal energy density of the plasma in the bubbles, and L_CR and ${\dot{E}}_{\mathrm{therm}}$ the power fed, respectively, into freshly accelerated cosmic rays and freshly heated plasma at the reverse shocks, we have that

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{{n}_{\mathrm{tot}}^{\mathrm{shell}}}={\left(\displaystyle \frac{{u}_{{\rm{p}}}^{\mathrm{FB}}}{{u}_{\mathrm{therm}}^{\mathrm{FB}}}\right)}^{3/4}\simeq {\left(\displaystyle \frac{{L}_{\mathrm{CR}}}{{\dot{E}}_{\mathrm{therm}}}\right)}^{3/4}\simeq {\epsilon }_{\mathrm{CR}}^{3/4},\end{eqnarray} \tag{ 6.8 }$$

where the last two near equalities rely on adiabatic evolution of both cosmic rays and thermal plasma downstream of the shocks in the interior volume R2. We thus find, implicitly defining the thickness of the shell, d_shell (assumed to be composed of two concentric shells, one around each spherical bubble)

$$\begin{eqnarray} {V}_{\mathrm{shell}} & \equiv & {V}_{\mathrm{FB}}-2\displaystyle \frac{4\pi }{3}{({r}_{\mathrm{FB}}-{d}_{\mathrm{shell}})}^{3}\\ & \simeq & \displaystyle \frac{5}{2}\displaystyle \frac{{k}_{B}\ {T}_{\mathrm{FB}}\dot{M}}{{({n}_{\mathrm{tot}}^{\mathrm{FB}})}^{2}\ {\rm{\Lambda }}[{T}_{\mathrm{shell}}]\ \mu {\rm{m}}}{\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}\right)}^{3/2}, \end{eqnarray} \tag{ 6.9 }$$

which implies a characteristic thickness of the shell in the geometrically thin limit

$$\begin{eqnarray} {d}_{\mathrm{shell}} & \simeq & 510\ \mathrm{pc}\ \left(\displaystyle \frac{{T}_{\mathrm{FB}}}{6.0\times {10}^{6}\ {\rm{K}}}\right)\left(\displaystyle \frac{\dot{M}}{0.1\;{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{0.004\;{\mathrm{cm}}^{-3}}\right)}^{-2}{\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}/0.1\right)}^{3/2}, \end{eqnarray} \tag{ 6.10 }$$

where we set A_shell = 2.5 × 10⁴⁵ cm². The total mass of the shell, in all generality, is

$$\begin{eqnarray} {M}_{\mathrm{shell}} & \simeq & \displaystyle \frac{5}{2}\displaystyle \frac{{k}_{B}{T}_{\mathrm{FB}}\dot{M}}{{n}_{\mathrm{tot}}^{\mathrm{FB}}\ {\rm{\Lambda }}[{T}_{\mathrm{shell}}]}{\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}\right)}^{3/4}\\ & \simeq & 2.5\times {10}^{7}\ {M}_{\odot }\left(\displaystyle \frac{{T}_{\mathrm{FB}}}{6.0\times {10}^{6}\ {\rm{K}}}\right)\left(\displaystyle \frac{\dot{M}}{0.1\;{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{n}_{\mathrm{tot}}^{\mathrm{FB}}}{0.004\;{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}/0.1\right)}^{3/4}, \end{eqnarray} \tag{ 6.11 }$$

and we may finally write the 1–100 GeV hadronic luminosity of the shell as

$$\begin{eqnarray} {L}_{\gamma }^{\mathrm{pp}} & \simeq & \displaystyle \frac{15\ {f}_{\mathrm{bol}}}{8}\displaystyle \frac{{({k}_{B}{T}_{\mathrm{FB}})}^{2}\dot{M}}{\mu {\rm{m}}{\rm{\Lambda }}[{T}_{\mathrm{shell}}]}\displaystyle \frac{{n}_{{\rm{H}}}}{{n}_{\mathrm{tot}}}\\ & & \times {\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}\right)}^{3/4}{\sigma }_{\mathrm{pp}}\ {\kappa }_{\mathrm{pp}}\ c\\ & \simeq & 2.4\times {10}^{37}\ \mathrm{erg}\;{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{T}_{\mathrm{FB}}}{6.0\times {10}^{6}\ {\rm{K}}}\right)}^{2}\left(\displaystyle \frac{\dot{M}}{0.1\;{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{L}_{\mathrm{CR}}^{\mathrm{FB}}}{{\dot{E}}_{\mathrm{therm}}^{\mathrm{FB}}}/0.1\right)}^{3/4}, \end{eqnarray} \tag{ 6.12 }$$

which compares favorably to the ∼2 × 10³⁷ erg s⁻¹ 1–100 GeV luminosity of the bubbles (Su & Finkbeiner 2012) if ≳10% of the mechanical power at the shock is delivered to cosmic rays. Note that in the interior region, we are normalizing to a somewhat hotter temperature than the 3.5 × 10⁶ K measured in X-rays, assuming that this is characteristic of the cooling shell and that here f_bol ≃ 0.4.

The results of a full calculation of the hadronic luminosity of the filaments in the downstream are plotted in Figure 9. Note that there is a subdominant component of the bubbles' hadronic γ-ray emission due to collisions between the cosmic-ray protons and the rarefied, volume-filling plasma, but, for fiducial numbers, this provides only 5%–20% of the integrated shell luminosity.

Zoom In Zoom Out Reset image size

6.2.4. IC Gamma-Ray Emission

As has been noted, the same population of nonthermal electrons synchrotron-radiating at microwave frequencies to explain the haze may emit IC γ-rays in the Fermi band with the correct bubble spectrum. Within our scenario, however, we find that this process is subdominant to the hadronic γ-ray emission explaining ∼10%–30% of the total observed 100 MeV–500 GeV γ-ray luminosity of the bubbles for the fiducial case that the power fed into freshly accelerated electrons is 10% of the power fed into all CRs (which is, itself, assumed to be 15% of the mechanical power at the shock in this figure): L_e = 0.1 L_CR = $0.1\times 0.15\times \dot{E}$; see Figure 9. Of course, this conclusion is contingent on the assumed L_e/L_CR, but significantly higher values of this ratio can be rejected on two considerations.

• 1.  
  Within our model, IC cannot reproduce the observed morphology of the bubbles' overall γ-ray emission alone as the high-energy electrons never reach more than ∼300 pc downstream of the shocks (see Figure 10). On the other hand, this expected concentration of IC emission does imply that, within a few degrees of the shocks, it can contribute easily more than half of the total 100 MeV–500 GeV γ-ray intensity even for ${L}_{e}/{L}_{\mathrm{CR}}\sim 0.1$, and we thus expect γ-ray substructure around the shocks; preliminarily, we find evidence for such (see Section 7.1.1).
• 2.  
  Values L_e/L_CR ≳ 0.1 are disfavored on theoretical grounds (Bell 1978), on the basis of analysis of the Galaxy's global cosmic-ray and nonthermal radiation budget (Strong et al. 2010), and from observations and analysis of nonthermal radiation associated with giant cluster shocks (Pinzke et al. 2013; Brunetti & Jones 2014, and references therein).

Zoom In Zoom Out Reset image size

6.2.5. Projected Surface Brightness

As has been remarked, the bubbles exhibit a rather flat surface brightness. We consider a toy, two-zone geometric model of a spherical shell of luminosity ${L}_{\mathrm{shell}}\equiv {\epsilon }_{\mathrm{shell}}\ {V}_{\mathrm{shell}}$ enveloping a spherical interior region of luminosity ${L}_{\mathrm{interior}}\equiv {\epsilon }_{\mathrm{interior}}\ {V}_{\mathrm{interior}}$. The outer radius of the shell is r_B, and the inner radius of the shell/outer radius of the interior is r_shell. In the far field limit, we find from elementary geometric considerations that a roughly flat projected surface brightness for such a two-zone model requires that the volumetric emissivities are related as

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\epsilon }_{\mathrm{shell}}}{{\epsilon }_{\mathrm{interior}}}\right|}_{\mathrm{flat}}\simeq {\left[1+{\left(\displaystyle \frac{{r}_{B}}{{r}_{\mathrm{shell}}}\right)}^{3}\left(\sqrt{1-{\left(\displaystyle \frac{{r}_{\mathrm{shell}}}{{r}_{B}}\right)}^{2}}-1\right)\right]}^{-1}.\end{eqnarray} \tag{ 6.13 }$$

This implies the geometric restriction that r_shell > 0.544r_B ∼ 1.7 kpc (or d_shell ≲ 1.4 kpc; see orange exclusion region in Figure 8). Over the range $0.544\ {r}_{B}\lt {r}_{\mathrm{shell}}\leqslant {r}_{B}$, one finds ${\epsilon }_{\mathrm{shell}}/{\epsilon }_{\mathrm{interior}}{| }_{\mathrm{flat}}\gt 4.3$. In Figure 9 the gray dashed curve shows where Equation (6.13) is satisfied exactly; left and upward of the curve, the overall γ-ray surface brightness morphology tends to center brightened, and right and downward of the curve, the morphology tends to limb brightened according to our simplistic analysis. A constant-volume emissivity model would produce a center-brightened morphology, and the original instantiation of the hadronic model (Crocker & Aharonian 2011) was criticized on this ground. The updated hadronic model presented here, if anything, tends toward limb brightening in the favored parameter space according to Equation (6.13). However, we do not believe that this is a problem for our model for three reasons: (i) the observed surface brightness is not exactly flat, indeed there is, e.g., a large region close to the east edge of the southern bubble that appears brightened ("the cocoon"; Su & Finkbeiner 2012; Ackermann et al. 2014); (ii) we have assumed an unrealistically simple geometry in deriving Equation (6.13); and (iii) we expect, in reality, some relative diminution of shell hadronic emission with a compensating increase in interior emission due to condensations falling down through the bubbles under gravity, still emitting γ-rays (see Figure 1). A more accurate treatment of this geometrical aspect of the problem requires numerical modeling.

6.3.1. Observations

6.3.2. Synchrotron Emission in Our Model

Before proceeding to detailed discussion, we explain briefly how our model accounts for the complicated synchrotron emission phenomenology of the GC's giant outflows. For clarity, note that in general we shall refer to the 23–44 GHz band as "microwave" and 2.3 GHz emission as "radio."

The nucleus directly accelerates a population of hard-spectrum primary electrons that is advected from the region before it can lose much energy radiatively in situ (Crocker et al. 2011a, 2011b; Crocker 2012; Carretti et al. 2013). These electrons cool adiabatically and radiatively in the expanding outflow (R1) until they encounter the giant reverse shock in either bubble where they are reaccelerated.

At latitudes in the bubbles higher than the shocks, there are then effectively three zones for synchrotron emission in our model.

• 1.  
  A zone downstream of each reverse shock, partially encompassing R2, where the electrons are being conveyed systematically downstream by the combination of the downstream flow and their streaming along magnetic field lines. In this zone, there is a one-to-one correspondence between distance from the shock and electron "age." The frequency-dependent extent of this region is given approximately by the requirement that the transport time from the shock be less than the total (synchrotron+IC) cosmic-ray electron loss time. The synchrotron spectrum within this zone is thus "uncooled," F_ν ∼ ν^−0.5, reflecting that the electron spectrum itself is simply the hard injection spectrum produced at the shock, ${{dN}}_{e}/{{dE}}_{e}\sim {E}_{e}^{-2}$, up to some (position- and frequency-dependent) maximum energy above which it is exponentially cut off (see Figures 10 and 11). The hard-spectrum microwave haze is largely attributable to the synchrotron radiation from this electron population that is radiating on the volume-filling magnetic field in the base of R2.
• 2.  
  A zone corresponding to the compression shell that covers, in projection, the entire solid angle of the γ-ray bubbles out to the CD, i.e., R3. In this zone, the adiabatically amplified magnetic field pervading the cooling gas is rather strong and regular. Synchrotron emission here is by (i) the adiabatically compressed, lower-energy primary electrons that have long enough loss time that they can travel all the way from the shock and (ii) secondary electrons that result from hadronic collisions occurring in the denser shell gas. Two important considerations for synchrotron emission from R3 are that (i) given relative transport and loss times microwave-emitting primary electrons accelerated at the shock never reach this zone; and (ii) in this zone, there is a mixture of electrons of different ages as all shell electrons simply cool in situ (with the qualification that a fraction of the electrons apparently leak from R3 into the radio plumes, R4; see below). The latter implies that, above some break energy, E_bk, the integrated, steady-state shell electron spectrum of R3 is ${{dN}}_{e}/{{dE}}_{e}\sim {E}_{e}^{-3}$ (Kardashev 1962), implying a synchrotron spectrum above ν_bk of F_ν ∼ ν⁻¹. This explains the steep, polarized, and rather uniform spectrum measured between 2.3 and 23 GHz over the solid angle of the radio lobes. The break energy is implicitly defined by ${t}_{{\rm{c}}}({E}_{\mathrm{bk}})={t}_{\mathrm{FB}}$, where t_c is the energy-dependent electron cooling time (with cooling processes dominantly synchrotron radiation and IC for the relevant electron energy range) and t_FB is the age of the bubbles over which the electrons have been accumulated in the shell. At even higher energies (for typical parameters) the spectrum suffers a further quasi-exponential cutoff owing to the fact that higher-energy primaries can never reach this zone; this cutoff is softened by the presence of secondary electrons produced in situ, which can supply microwaves even at these large distances from the shock (see Dobler 2012b).
• 3.  
  A zone identified with the polarized radio plumes that extend beyond the boundaries of the north and south bubbles at high latitudes. For this zone, a one-to-one correspondence between distance from the bubbles' edges (the CD) and electron age can be reestablished as the electrons are systematically advected into the plumes. We thus expect, in principle, that the position-dependent electron distribution takes the form of the injection spectrum (already steepened to ∼E_e ⁻³), with an exponential cutoff that moves to progressively lower electron energy with increasing distance from the edges. With only two spectral data points currently available (at 2.3 and 23 GHz), this cutoff translates into a steepening with distance into the plumes (matching observations).

Zoom In Zoom Out Reset image size

6.3.3. Magnetic Fields

6.3.4. Synchrotron Emission from the Shock Downstream

Cosmic-ray transport downstream of the shock is a combination of the ∼100–150 km s⁻¹ flow of the shocked, downstream plasma and cosmic-ray streaming along the field lines at the Alfvén velocity (Kulsrud & Pearce 1969; Skilling 1971) of ∼240 km s⁻¹. The range microwave-emitting electrons can reach in their loss time is then a function of magnetic field amplitude, which controls both synchrotron losses and the Alfvén velocity and reaches a maximum for B ≃ 10 μG; see Figure 10. Here the cooling time, t_c, is the time taken for an electron injected at the accelerator with the maximal possible energy, E_max, to be cooled to the energy at which it synchrotron-radiates into the frequency range of concern, ${E}_{{\rm{c}}}\equiv {m}_{e}\ {c}^{2}\ \sqrt{\nu /{\nu }_{{\rm{g}}}[B]}$, where ${\nu }_{{\rm{g}}}[B]$ is the electron gyroradius in magnetic field B. The energy of an electron injected at E₀ into a constant magnetic field B evolves according to

$$\begin{eqnarray}&&E[t]=\displaystyle \frac{{E}_{0}}{1+t/\tau [{E}_{0}]},\end{eqnarray} \tag{ 6.14 }$$

where $\tau [E]$ = $1.3\times {10}^{8}{(B/10\;\mu {\rm{G}})}^{-2}$ ${(E/\;\mathrm{GeV})}^{-1}\;\mathrm{yr}$, ignoring non-synchrotron losses (particularly IC) for illustrative purposes and assuming pitch-angle scattering on timescales short with respect to the loss times (e.g., Reynolds 2009).

With this description of the energy evolution of the electrons, one obtains that an electron spectrum injected as a power law by the accelerator, ${{dN}}_{e}[{E}_{0}]/{dE}\equiv {{KE}}_{0}^{-s}$, evolves according to

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{e}[E,t]}{{dE}}={{KE}}^{-s}{(1-t/\tau [E])}^{s-2}.\end{eqnarray} \tag{ 6.15 }$$

The evolution of a cosmic-ray electron spectrum, injected with s = 2.1, as the electrons are conveyed downstream of a shock and, simultaneously, suffer synchrotron+IC losses is illustrated in Figure 11.

Once the transport time exceeds the cooling time t_c, a cutoff feature is introduced into the electron spectrum. Applied to the outflow population, GeV IC-emitting electrons only reach a few × 10² pc past the shocks and microwave-emitting electrons reach a few kpc.

6.3.5. Synchrotron Emission from the Compression Shell

In contrast, lower-energy ∼GHz emitting electrons survive longer, reach the shell before the CD, and are accumulated inside it (see Figure 10). Compression of the cooling plasma in this shell leads to adiabatic amplification of the frozen-in magnetic field and cosmic rays, implying a second reenergization of those electrons that actually survive to this point. The spatially coincident, enhanced gas and cosmic-ray ion populations close to the bubbles' edges or in the filaments also lead, however, to significant hadronic collisions, which, in addition to producing hadronic γ-ray emission, also inject secondary electrons (and positrons) via the decay of collisionally produced charged mesons. The edges/filaments are thus populated by both enhanced primary and secondary electron populations in addition to adiabatically compressed magnetic fields. The electrons do not escape from these regions but lose their energy in situ to synchrotron radiation. Accumulation and trapping of electrons in this region lead to mixing of cosmic-ray electron populations of various ages. Cooling + mixing (in projection) of differently aged electron populations will be enhanced if there are downdrafts toward the edges of the bubbles or convection within it. Mixing steepens the overall steady-state population into an $\propto {E}_{e}^{-3}$-type distribution (Kardashev 1962). The steep-spectrum, polarized emission detected at 2.3 GHz (Carretti et al. 2013) is largely attributable to this population that radiates in the enhanced magnetic field region associated with the cooling shell toward the edges of the bubbles (and possibly also swept-up halo magnetic field outside the CD). We can thus evaluate the 2.3 GHz emission from the bubbles in a thick-target, single-zone model; the emission thus calculated is displayed in Figure 12 (green curves).

Zoom In Zoom Out Reset image size

A prediction of our model is the existence of a break in the shell radio spectrum corresponding to where the synchrotron+IC cooling of shell electrons becomes less than the age over which they have been accumulated, roughly the age of the bubbles t_FB. The break occurs at a frequency of

$$\begin{eqnarray}&&{\nu }_{\mathrm{br}}=1.96\times {10}^{8}\;\mathrm{Hz}{\left(\displaystyle \frac{B}{10\ \mu {\rm{G}}}\right)}^{-3}{\left(\displaystyle \frac{{t}_{\mathrm{FB}}}{{10}^{8}\;\mathrm{years}}\right)}^{-2}.\end{eqnarray} \tag{ 6.16 }$$

The fact that the 23–2.3 GHz spectrum of the radio lobes is ∝ν^−1.1 with the break necessarily below 2.3 GHz implies corroborative evidence that the outflow structures have an age of >(2.9 × 10⁷ yr)[B_shell/(10 μG)]^−3/2(ν_br/2.3 GHz)^−1/2.

6.3.6. Synchrotron Emission from the Radio Plumes

6.3.7. Energetics of Synchrotron Emission and Secondary Electrons

7.1.1. South Bubble

7.1.2. North Bubble

In the northern bubble there is only a faint feature possibly corresponding to the depolarization feature described above (see Figure 13). The north, however, is highly confused and suffers from depolarization at 2.3 GHz significantly off the plane owing to ρ Ophiucus and other local, high-latitude gas structures. There is, however, a feature of the northern bubble that is highly suggestive of a shock. The Galactic Center Spur (GCS; Sofue et al. 1989; Carretti et al. 2013; see Figure 13) is a highly collimated total- and polarized-intensity radio continuum feature that extends in a northwest direction from just east of the Galactic center to Galactic latitudes of b ∼ 25°. Carretti et al. (2013) have claimed that this feature wraps around the forward surface of the northern bubble, an interpretation we adhere to.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the basis of a geometric analysis of the overall curvature of this structure, Carretti et al. (2013) also claimed, however, that the spur provides evidence that gas is being driven aloft into the bubbles at a speed of ∼1000 km s⁻¹ up to heights of at least a few kpc. This analysis was supported by an independent argument that if the 2.3 GHz emitting electrons were accelerated (exclusively) in the nucleus, they would have to be transported at a similar speed in order for them to reach the top of the bubbles in their loss time. This high velocity over the entire bubbles is not consistent with our present model, which suggests that while gas may be injected into the base of the bubbles at ≳500 km s⁻¹, it is slowed to 100–150 km s⁻¹ upon encountering the giant shocks at heights of ∼1 kpc.

Moreover, the assumptions leading up to the estimated ∼1000 km s⁻¹ presented by Carretti et al. (2013) are no longer met in our revised scenario. In particular, we now expect injection of high-energy electrons at large scales in the bubbles either via direct acceleration at the large-scale shocks or as secondaries. The shock in the north bubble changes our expectation for the geometry of the GCS. In fact, we expect an inflexion point in the GCS at the shock as the vertical speed of the flow is slowed by ∼1/4 crossing the shock. Such an inflexion does seem to be present in the GCS, and it has the correct geometry: at b ∼ 17° the angle between the long axis through the GCS changes from ∼17° to ∼48° (see Figure 13), suggesting a change of speed $\sim \mathrm{tan}[17^\circ ]/\mathrm{tan}[48^\circ ]=0.27$, consistent with a strong shock.

Furthermore, an independent determination (Vidal et al. 2014) of the spectral index of the polarized emission between 23 and 33 GHz from the GCS reveals a hardening at higher latitudes for the radio spectral index (b ≳ 18°) from 1.52 ± 0.51 to 0.62 ± 0.44 (with a similar hardening in a region coincident with part of the eastern edge of the north bubble), consistent with a giant shock—at the correct height above the plane—reaccelerating electrons, as we have described above.

Another piece of corroborative evidence in support of this interpretation of the shock geometry in the north is the following: consideration of Figure 13 shows that the southeast region of the north bubble is comparatively dim in 2.3 GHz polarized emission. While, in general, depolarization due to magnetized plasma along the line of sight may affect emission up to latitudes of ∼10°, this region extends to ∼20° close to the eastern edge. Moreover, the very same depolarization region is evident also at 23 GHz (see Figure 14), a frequency high enough to be completely unaffected by Faraday depolarization in this region of the sky. We expect, therefore, that this lack of detected polarized emission betokens a lack of intrinsic emission along this line of sight. Significantly, this dim region is below the putative shock feature. This can be naturally explained if, on the one hand, the high-energy primary electrons requisite to produce the emission do not reach these heights from the nucleus and, on the other hand, are only reaccelerated to microwave-emission energies at the shock that is further from the nucleus.

7.1.3. Why Are the Putative Shock Features Not Parallel to the Plane?

In Section 1 we described UV absorption evidence (Keeney et al. 2006; Zech et al. 2008) for higher-altitude warm ionized material in the Galactic bulge that is participating in a nuclear fountain flow. We believe that this material is the end state of gas condensing in the bubbles' shell. This warm ionized material seems to be concentrated toward the edges of the bubbles, just as required if the volumetric γ-ray emissivity is to peak here and as expected for our cooling shell model. Keeney et al. (2006) determine that their clouds are likely launched into a conical outflow that is collimated into a cylinder of radius ∼1.6 kpc; with this (somewhat uncertain) radius determination they find that their clouds reach maximum heights z_max ≃ 12 ± 1 kpc. Given the qualitative resemblance between this geometry and that of the bubbles, we note that, were the radius used in the Keeney et al. (2006) cloud trajectory calculations rescaled to 3 kpc,¹⁴ the clouds would be inside the bubbles or the radio plumes.

Zech et al. (2008) have measured the metallicity of the fountain material found by them to be significantly supersolar metallicity (∼1.6 Z_⊙). It is, thus, a very good prospect for having been launched from the GC. Moreover, the illuminating background UV source is a post-AGB star in the Messier 5 globular cluster at known distance (7.5 kpc) and height (5.3 kpc) above the plane, which, in projection, is coincident with the edge of the northern bubble. It seems more than plausible that the cloud is located in the compression shell (R3) and proffers a lower limit to the radius of the northern bubble of ∼2.7 kpc (at z = 5.3 kpc), entirely consistent with the ∼3 kpc radius we have adopted.

Fox et al. (2014) have recently reported on further UV absorption observations along the sightline to the quasar PDS 456 (l, b) = (104, +112), which passes through the conically expanding lower part of the northern bubble at a height of ∼2.3 kpc. These observations find red- and blueshifted warm ionized material, which kinematic analysis suggests is moving at ∼900 km s⁻¹, likely accelerated by the ram pressure of the dynamically dominant "wind fluid," i.e., the volume-filling, hot plasma phase. Given this supersonic velocity, this sightline must pass somewhere below the giant shocks. Indeed, the northern shock position we infer from the radio polarization observations discussed in Section 7.1.2 is safely above the sightline. The ∼900 km s⁻¹ speed inferred by Fox et al. (2014) is somewhat higher than the flow speeds indicated by our model and the heights somewhat above where we predict the shocks to lie on the basis of modeling the global flow parameters. It is tantalizing, however, that this sightline is on the east side of the north bubble, where, observationally, the shock seems to have been pushed furthest from the plane. This could be a spatially and temporarily localized upward excursion in the outflow velocity connected to the recent concentration of star formation activity on the east side of the Galactic center, as already discussed. To obtain a shock ∼2.3 kpc height with an upstream flow requires ${\dot{E}}_{0}/\dot{M}\simeq 5\times {10}^{15}$ erg g⁻¹ or a mass loading factor β ≃ 1.7, only a little smaller than the range of ∼2–4 we find for this parameter in the allowed parameter space from our modeling of the global outflow.

One possible tension between our model and X-ray observations is the following: modeling of data taken by Suzaku in high-latitude scans across the limbs of both the north and south bubbles (Kataoka et al. 2013) suggests a significantly subsolar metallicity $Z\simeq 0.2\;{Z}_{\odot }$. This is contrary to the simple expectation were this material to be supplied by the nucleus (in which case it would be expected to be at least solar in metallicity). Interestingly, the warm plasma found by Zech et al. (2008) along the Messier 5 sightline, relatively close to the Suzaku northern scan, is supersolar in metallicity, as already remarked. One possible explanation for this apparent tension is that the Suzaku observations predominantly reveal emission from swept-up halo material on the other side of the CD, which various pieces of evidence (as summarized by Miller & Bregman 2014) suggest has a metallicity of ∼0.3 Z_⊙. The northern Suzaku scan does commence significantly outside the edges of the bubbles, within the NPS feature; with six of the eight northern pointings outside the γ-ray edge and with the metallicity measured only with confidence in the aggregate (and not for individual pointings), a systematic bias toward the metallicity of exterior gas would be expected in our picture that the γ-ray edge is identified with the CD.

We also note here that the ∼300 km s⁻¹ expansion speed of the bubbles estimated by Kataoka et al. (2013) is also clearly in tension with our claimed expansion speed of only a ∼few km s⁻¹. The 300 km s⁻¹ estimate, however, is not model independent: it relies on the interpretation that the increase of the temperature of the plasma in the vicinity of the bubbles' limbs with respect to the background halo plasma by a factor of ∼3/2 can be explained by heating at a weak shock sweeping up the halo gas (which is assumed to be initially in hydrostatic equilibrium). In our model, the plasma inside the CD is heated by the reverse shocks low down in the nuclear outflows into both Galactic hemispheres; we do not specifically deal with the swept-up halo material here (which the Suzaku observations are biased toward, as already mentioned).

Another issue with the X-ray data may be the following: these have been interpreted by Kataoka et al. (2013) using a three-component model incorporating (i) an absorbed power-law isotropic background; (ii) a local (unabsorbed) single-temperature, collisional-ionization-equilibrium (CIE) plasma foreground (representing the local plasma bubble); and (iii) an absorbed, single-temperature, CIE plasma that represents the signal region (the NPS and the limb regions in either the north or south bubbles). Such a model may be too simple to adequately describe the putative cooling shell. Indeed, Weaver et al. (2000) and Strickland et al. (2002) demonstrate that fitting to X-ray plasma emission from starburst outflows often leads to unphysically low metallicity estimates. Possible circumstances leading to a too-low metallicity estimate are (i) the existence of neglected multitemperature plasma distributions, (ii) having a large and poorly constrained absorbing column, and/or (iii) nonthermal contribution to the continuum (that boosts apparent free–free emission relative to lines, thus making plasma appear to have smaller metallicity). Any or all of these may be relevant effects for the bubble plasma, in particular, we expect the gas in the cooling shells to be multitemperature, as already discussed, and catastrophic cooling of the plasma into atomic phase may also generate absorbing material.

Mention of halo material swept up ahead of the CD leads us to consideration of the forward shock in our model. In Section 4 we do not consider cooling in the region between the forward shock and the CD. The implication of this is that, with the attainment of pressure balance with the external atmosphere, the forward shock degenerates into a sound wave. Even traveling at only the ∼280 km s⁻¹ sound speed, however, this would have traveled ∼30 kpc t/10⁸ yr and, therefore, be currently out of the picture given our estimated age for the bubbles of ≳ few × 10⁸ yr. Yet, observationally, there are, in the northern Galactic hemisphere, very large angular scale X-ray and nonthermal radio and γ-ray structures that, at least in projection, are external to the bubbles' edges and appear concentric to them. One is tempted to identify these as shock or sound waves associated with the bubbles (e.g., Figure 4 of Su et al. 2010; Figure 1 of Jones et al. 2012.)

In particular, the NPS (Bolton & Westfold 1950) emerges north from the Galactic plane at l ∼ 20°–30° and rises to heights above b ∼ 80°, merging with the Loop I radio structure (Large et al. 1962; Haslam et al. 1982). It parallels the east edge of the north bubble over this extent; indeed, the NPS and Loop I, together, parallel the edge of the northern polarized radio lobe (which subsumes the north bubble) to a remarkable degree. The NPS exhibits total-intensity and polarized, nonthermal radio, and microwave emission (Sofue & Reich 1979; Sun et al. 2014; Vidal et al. 2014), thermal X-ray emission (Snowden et al. 1997; Willingale et al. 2003), and ∼GeV γ-ray emission (Casandjian et al. 2009; Ackermann et al. 2014).

It has long been argued on the basis of stellar polarization data that the NPS is a structure of the local ISM (within a few × 10² pc; Bingham 1967; Berkhuijsen et al. 1971). However, a number of recent works have claimed, on the basis of various independent analyses of X-ray, radio polarization, and other data, that it is actually a Galactic-scale feature (Kataoka et al. 2013; Fang & Jiang 2014; Sun et al. 2014; Sofue 2015; Tahara et al. 2015) or perhaps a confusion of two structures, one local and one at Galactic scales. This finding complements the idea that (some part of) the NPS is a signature of Galactic nuclear activity (Sofue 1977, 2000; Bland-Hawthorn & Cohen 2003; Totani 2006) and perhaps even related directly to the Fermi bubbles (Guo & Mathews 2012; Kataoka et al. 2013; Mou et al. 2014). In addition to the NPS, there are linear γ-ray and polarized radio arcs that run concentrically between the NPS and the eastern edge of the northern bubble (Su et al. 2010; Jones et al. 2012); fainter radio and X-ray spurs have also been claimed as the northwestern, southeastern, and southwestern complements of the NPS (Sofue 2000, 2015; Tahara et al. 2015). We set aside further discussion of this interesting phenomenology for a future publication. We do, however, make the claim here that accounting, in a fuller model, for cooling in the shell of material exterior to the CD implies a different channel for dissipation of the flux of internal energy of swept-up halo gas. This retards the advance of the forward surface of this shell (which, in the absence of cooling, would be the forward shock, degenerating into a sound wave with the attainment of pressure balance), leaving open the possibility that this feature could be identified with the NPS + Loop I feature. Note that the question of why the northeastern feature is so much more conspicuous than its northwest, southeast, or southwest analogs is not addressed in the current speculation.

In subsequent publications we will present or examine in more detail the detailed spectral morphology of the synchrotron emission in the bubbles, the observational evidence for the giant reverse shocks we have inferred, a detailed treatment of the dynamical production of the dense condensations, and a full treatment of the dynamic importance of cooling in the zone forward of the CD and the possibility to identify the NPS with the forward surface of a cooling shell of swept-up halo gas. We will also examine the origin of the γ-ray spectral downturn at ∼100 GeV (Ackermann et al. 2014) and the hardening of the γ-ray spectrum at high latitudes (Selig et al. 2014; Yang et al. 2014). The former is not explicitly dealt with in our current model but may signify the energy where the bubbles themselves cease to trap cosmic-ray protons effectively; the latter may signify an important role for diffusive transport of cosmic rays in the upper reaches of the bubbles (Yang et al. 2014), which is, again, not accounted for in the current model. Also awaiting a detailed treatment is the dynamical importance of the strong and coherent magnetic field that appears to be wrapped over the surface of the bubble, particularly in stabilizing this surface.

For fiducial parameters, the giant shocks are capable of accelerating cosmic-ray protons to energies around the knee in the observed cosmic-ray spectrum. Given the very hard spectrum, a total power in these particles that approaches 10³⁹ erg s⁻¹, and that the particles are delivered into the halo but not too far above the disk, these shocks are interesting candidates (see Parizot 2014) for the main accelerators of Galactic cosmic rays in this energy range (see Cheng et al. 2012; Lacki 2014; Taylor et al. 2014).

Our model selects a region of the $\{{\dot{E}}_{0},\dot{M}\}$ parameter space implying that, if star formation driven, the outflow is mass loaded at a level β ∼ 2–4, broadly consistent with observations of external, star-formation-driven outflows (Strickland & Heckman 2009, and references threrein). Further, it is remarkable that our results indicate that the mass efflux in the outflow is $\dot{M}$ ∼ 0.03–0.1 M_⊙ yr⁻¹, rather close to the current rate of nuclear star formation (Crocker 2012; Longmore et al. 2013, and references therein), and in turn close to the long-time-averaged nuclear SFR.¹⁵ Moreover, on the basis of the semianalytic model presented in Section 4, we can say that the Fermi bubble system is rather delicately poised in terms of allowable values of the ratio ${\dot{E}}_{0}/\dot{M}$: systems with ${\dot{E}}_{0}/\dot{M}$ ratios half those favored by our analysis asymptote to final radii within ∼10⁸ yr that are smaller than those observed (and contain plasma that is too cold); systems with twice the favored value for this ratio (of $\sim 1.6\times {10}^{15}$ erg g⁻¹) do not asymptote to a final radius until ∼10¹⁰ yr (and this radius is 10 kpc, significantly larger than the observed bubbles) and spend only a brief time in the "vicinity" of the observed size (implying a return to a "Why now?" fine-tuning); see Figure 4.

These observations suggest the operation of a negative feedback process to maintain the required delicate balance. An obvious mechanism to achieve this—completing the cycle—is that material fountaining back to the plane in the vicinity of the nucleus will trigger further star formation in this region. Star formation that is too vigorous or that produces a too large value of ${\dot{E}}_{0}/\dot{M}$ in the outflow (which may currently be happening on the east side of the nucleus) will tend to switch itself off because the outflowing plasma takes too long to condense and fall back from the downstream of the bubbles. The matter cycle will not be 100% efficient, of course, but the system is not closed: some amount of accretion to the nucleus will continue through the plane, and, moreover, interaction and mixing between bubble matter and the surrounding halo gas around the CD may catalyze the condensation of halo gas so that, overall, more gas fountains back to the plane than is lofted upward (Marinacci et al. 2010), some of it low metallicity.