THE INTERACTION OF THE FERMI BUBBLES WITH THE MILKY WAY'S HOT GAS HALO

We summarize the XMM-Newton emission line sample compilation here, but we refer the reader to HS12 for a more detailed description of their data reduction methods. Their initial sample included 5698 observations that had any EPIC-MOS exposure time. They reduced the data using the XMM-Newton Extended Source Analysis Software¹ (XMM-ESAS; Kuntz & Snowden 2008; Snowden & Kuntz 2011), which includes screening the 2.5–12 keV band count rate for soft proton flares. They also removed point sources from the spectral extraction regions using data from the Second XMM-Newton Serendipitous Source Catalog² (Watson et al. 2009), as well as visual inspection for bright sources in the images. The authors attempted to reduce geocoronal SWCX emission by excluding observing periods with high solar wind proton flux measurements (see Section 2.3 for the details of this procedure), which they called their "flux-filtered" sample. These reduction methods resulted in 1868 total observations and 1003 flux-filtered observations with $\geqslant 5$ ks of good observing time. Each of these observations includes an O vii and O viii emission line intensity measurement, and they have been used to analyze the known emission sources (i.e., Local Bubble (LB), extended hot halo, SWCX).

This sample has been used before to successfully model the Milky Way's global hot gas structure, thus motivating its use to model the Fermi bubbles. MB15 compiled a subset of the HS12 flux-filtered sample with additional spatial screening criteria to reduce any residual emission from sources other than the LB and Galactic hot halo. To achieve this, they removed observations within 05 of potential bright X-ray sources (see Table 1 in MB15 for the types of sources) or within 10° of the Galactic plane, and sight lines near the Fermi bubbles ($| l| \leqslant 22$°, $| b| \leqslant 55$°). This resulted in a subsample of 649 observations from the HS12 flux-filtered sample.

The XMM-Newton data used in this study are the same as those used used in MB15, but including the sight lines near the Fermi bubbles. We start with the HS12 flux-filtered sample and still exclude sight lines near bright X-ray sources or within 10° of the Galactic plane. These screening criteria result in a total of 683 XMM-Newton measurements distributed across the sky, with 34 measurements passing near or through the bubbles' gamma-ray edge.

We compiled an initial Suzaku target list of all observations that were publicly available as of 2015 January and near the Fermi bubbles. This included any observations with Galactic coordinates $| l| \leqslant 25$° and 10° $\leqslant | b| \leqslant 55$°. There were 143 observations in this region of the sky, which we inspected for usable spectra.

Each observation includes data from the three active X-ray Imaging Spectrometer (XIS) detectors on board Suzaku (Koyama et al. 2007). The detectors each have an $18^{\prime} \times 18^{\prime}$ field of view and a point-spread function of $\approx 2^{\prime}$. We only considered data from the back-illuminated XIS1 detector since it has a larger collecting area below 1 keV than the other detectors.

We processed all XIS1 data using HEAsoft version 6.17 and calibration database files from 2015 January. We followed the standard data reduction procedure described in The Suzaku Data Reduction Guide.³ This includes recalibrating the raw data files, screening for flickering or bad pixels, energy-scale reprocessing, and building good time interval (GTI) files. Fortunately, the FTOOL script aepipeline performs all these tasks for standard Suzaku observations. We used aepipeline version 1.1.0 to generate reprocessed images, light curves, and spectra for our analysis.

We extracted all data products using xselect version 2.4. Each data set combined data from 3 × 3 and 5 × 5 editing modes where applicable. We did not impose any non-standard criteria on the data extraction steps with the exception of the cutoff rigidity (COR) of the Earth's magnetic field. This parameter varies throughout Suzaku's low Earth orbit, and larger COR constraints result in fewer particle background counts. A higher cutoff than the standard COR > 4 GV has been used in previous SXRB spectral fits, and typically results in higher signal-to-noise ratio between 0.5 and 2.0 keV (Smith et al. 2007; Kataoka et al. 2015). Here, we follow the suggested value from Smith et al. (2007) and use a constraint of COR > 8 GV.

Since the observational goal of this study is to measure O vii and O viii intensities from SXRB spectra, we removed point sources from each observation before we extracted spectra. We inspected each XIS1 image for point sources to remove with self-defined region files. Observations of exceptionally bright sources (i.e., where there was clear emission extending over more than $\approx 5^{\prime}$), extended sources (galaxy clusters, star clusters, etc.), or other anomalous features were rejected from the data set. We also rejected any observations not taken in the standard observing window mode due to the reduction in field-of-view collecting area. For any remaining visible point sources, we defined circular exclusion regions between 1' and 4' in radius centered on each source. We then re-extracted the data products for each observation with the point source region excluded.

The next step in our data cleaning process was light curve inspection and filtering. We extracted 0.4–10.0 keV light curves and constructed count-rate histograms for each observation. Our default screening criterion was to remove observing periods that were $\gt 2\sigma$ from the mean count rate. This led to a small reduction in observing time since most count-rate histograms followed an approximately Gaussian distribution. We flagged observations that did not have an approximately Gaussian count-rate distribution, which is indicative of additional soft proton flares or residual point-source emission that may have been variable throughout the observation. Example light curves with the various filtering tasks can be seen in Figure 1. We also expand on this analysis step in Section 2.3 where we exclude observing periods with high solar wind proton flux measurements. This light curve filtering created new GTI files, which we used to compile our initial data products.

Zoom In Zoom Out Reset image size

The procedures for point-source exclusion and light curve filtering outlined above were the primary stages in our initial sample catalog. To summarize the main screening criteria, we excluded observations that showed anomalous features in either their images or light curves or that had exceptionally bright point sources, and we filtered the images for removable point sources. After these screening criteria, we kept observations with $\geqslant 5$ ks of good observing time. This resulted in 112 observations out of the original 143.

SWCX emission can occur in any X-ray observation, but it is difficult to predict or attribute the amount of SWCX emission in individual SXRB observations (e.g., Carter & Sembay 2008; Carter et al. 2011; Galeazzi et al. 2014; Henley & Shelton 2015). The emission can occur either at the interface between the solar system and interstellar medium (ISM) (heliospheric emission) or as solar wind ions pass near the Earth's neutral atmosphere (geocoronal emission). Heliospheric emission is thought to vary with the overall solar cycle, the observed direction relative to the Sun's orbit and ecliptic plane, and the hydrogen and helium density in the neutral ISM (Cravens et al. 2001; Robertson & Cravens 2003a; Koutroumpa et al. 2006, 2007, 2011; Galeazzi et al. 2014; Henley & Shelton 2015). Geocoronal emission is thought to depend on the solar wind proton flux and the observed direction relative to the magnetosheath (Snowden et al. 2004; Wargelin et al. 2004; Fujimoto et al. 2007; Carter & Sembay 2008; Carter et al. 2010, 2011; Ezoe et al. 2010, 2011; Ishikawa et al. 2013; Henley & Shelton 2015). It is still unclear what the typical amount of SWCX emission is in a given X-ray observation, and models predict a wide range of O vii and O viii intensities depending on the parameters listed above (Robertson & Cravens 2003b; Koutroumpa et al. 2006, 2007, 2011; Robertson et al. 2006). For the purpose of this project, SWCX emission is considered to be contamination, and our goal is to reduce the amount of potential emission as much as possible.

Following the work of HS12, we filter the observations for periods of high solar wind proton flux in an effort to reduce geocoronal SWCX. For each observation, we gathered solar wind data from the OMNIWeb database,⁴ which includes data from the Advanced Composition Explorer and Wind satellites. The database includes solar wind densities and velocities, which we convert to fluxes. We cross-correlated each solar wind proton flux light curve with the X-ray light curves. Periods with solar wind flux values $\gt 2\times {10}^{8}$ cm⁻² s⁻¹ were flagged and considered to potentially include SWCX emission. We illustrate how this screening works for several example spectra in Figure 2. We made new GTI files incorporating these filtering periods and the $\gt 2\sigma$ count rate periods discussed above. These GTI files were used for the final spectral extraction used in the fitting analysis.

Zoom In Zoom Out Reset image size

We point out that this filtering procedure is designed to reduce geocoronal SWCX emission, not heliospheric SWCX emission. The models suggesting that heliospheric SWCX varies with ecliptic latitude imply that applying an ecliptic latitude cut to an observing sample may help reduce that emission. Indeed, Henley & Shelton (2013) discuss this effect and employ this screening criterion for their study on fitting SXRB spectra. However, the analysis in MB15 argues that there does not appear to be a significant enhancement of O vii or O viii line emission within 10° of the ecliptic plane, part of which passes through the Fermi bubbles. Therefore, heliospheric SWCX is likely present at softer X-ray energies, but it does not appear to be a significant emission source for the oxygen lines of interest.

This additional screening procedure can only reduce the good exposure time in a given observation. Some observations occurred entirely during a period of high solar wind proton flux, in which case the observation was removed from the sample. Other observations occurred entirely during a period of low solar wind proton flux, in which case the observation was unaffected. The rest of the observations were partially contaminated, leading to a reduction of observing time. We enforced the same good exposure time requirement noted above of >5 ks to keep observations in the sample. After the default screening outlined in Section 2.2 and this additional SWCX filtering, our sample includes 58 of the original 143 observations.

This section outlines our spectral fitting procedure, including the response files used, our spectral model, and resultant data products. We extracted spectra in the 0.4–5.0 keV band, which is broad enough to model the known SXRB emission sources. Each observation had its own particle spectrum, or non-X-ray background (NXB) spectrum, and response files. We used Xspec version 12.9.0 for all spectral fitting, where we used the Cash statistic as our fit statistic (Cash 1979). Figure 3 shows an example observed spectrum and best-fit multi-component model. Our final result includes O vii and O viii line intensities for each observation.

Zoom In Zoom Out Reset image size

We generated response matrices and auxiliary response files using standard Suzaku scripts. The script xisrmfgen was used to generate response matrices (RMF files) for each observation. We used the ray-tracing program xissimarfgen to generate ancillary response files (ARFs) assuming a uniformly emitting source and a 20' radius circle for a simulated source region. The size of the source region acts as a normalization in our spectral fit values.

Each observation has an NXB spectrum collected by Suzaku observations of the Earth at night. We generated NXB spectra using the script xisnxbgen (Tawa et al. 2008). This creates an NXB spectrum using a weighted sum of Suzaku observations of the Earth at night based on the exposure times. The only unique parameter we supplied to the script is the same constraint of COR > 8 GV that we applied to the initial data extraction. Different SXRB studies have multiple treatments for these NXB spectra, with groups either subtracting the NXB counts from the observed spectrum (e.g., K15) or including the NXB spectrum as additional data and simultaneously modeling its contribution to the full spectrum (e.g., Smith et al. 2007). We follow the latter methodology, meaning we fit both our observed spectrum and the NXB spectrum as one process.

Our spectral model includes the following components: NXB spectrum, an absorbed cosmic X-ray background (CXB) or extragalactic power law, an absorbed hot gas continuum component with no oxygen emission lines, and the O vii and O viii lines of interest. The absorbed components were attenuated using the phabs model in Xspec (Balucinska-Church & McCammon 1992; Yan et al. 1998) and had column densities fixed to values from the LAB survey (Kalberla et al. 2005). We also assume metal abundances from Asplund et al. (2009) unless otherwise stated. The rest of this section details the spectral model assumed for each source.

The NXB model includes a contribution from particles hitting the XIS detectors that are not focused by the telescope and three instrumental lines. For the particle spectrum, we include a power law where both the normalization and slope can vary. The instrumental lines include an Al K line at 1.49 keV, an Si K line at 1.74 keV, and an Au M line at 2.12 keV.⁵ We model each of these lines as Gaussians with widths and normalizations left to vary as free parameters and the centroids fixed. This model component is only folded through the RMF response (as opposed to both the RMF and ARF), and it contributes to both the observed and NXB spectra.

The CXB spectrum is typically modeled as an absorbed power law or multiple broken power laws, and is thought to be due to unresolved active galactic nuclei (AGNs). The differences between these spectral shapes have minimal effects on the measured oxygen line intensities, since the CXB contributes $\lesssim 10 \%$ to the total SXRB flux below $\approx 1\,{\rm{keV}}$. Therefore, we adopt the spectral shape used by HS12, which is a power law with fixed spectral index of 1.46 (Chen et al. 1997). These authors discuss the uncertainty in the CXB power law normalization, but argue for a nominal value of 7.9 photons cm⁻² s⁻¹ sr⁻¹ keV⁻¹ at 1 keV after accounting for CXB sources with ${F}_{X}\lt 5\times {10}^{-14}$ erg cm⁻² s⁻¹ in the 0.5–2.0 keV band (Moretti et al. 2003). We allow the CXB normalization to be a free parameter in our spectral model, but place ±30% hard boundaries on the parameter to allow for field-to-field variation.

Although the oxygen lines are the measurement of interest, we still account for hot gas continuum emission in our model. We model this component as an absorbed thermal APEC plasma (Smith et al. 2001) with fixed solar metallicty and without the oxygen lines. We achieve the latter by setting the oxygen line emissivities to zero in the standard Xspec APEC files (Lei et al. 2009; Henley & Shelton 2012). The normalization and temperature were left as free parameters in the model. We expected the fitted plasma temperatures to be between 0.1 and 0.3 keV, but these temperatures are typically most sensitive to the oxygen lines that we disabled. Therefore, we let the plasma temperature vary outside this range, but with hard boundaries between 0.05 and 5 keV.

Our final spectral model components are the O vii and O viii emission lines. We model each of these components as Gaussian features with the widths fixed to the instrumental resolution. We fixed the O viii centroid to its laboratory value of 0.6536 keV, but let the O vii centroid vary since the line is an unresolved blend of the resonance, forbidden, and intercombination lines. The Gaussian normalizations were also free parameters, because these represent the line strength measurements. We point out that this spectral fitting method measures the total emission line strengths from all emission sources (residual SWCX, LB, absorbed hot gas halo, or Fermi bubble) for each sight line. This is why our model line intensities in Section 3.5 include all emission sources when comparing to the total observed line intensities.

The North Polar Spur (NPS) is an extended region of enhanced X-ray emission that is cospatial with part of a larger region of enhanced radio continuum emission known as Loop I (Berkhuijsen et al. 1971; Borken & Iwan 1977; Snowden et al. 1997). The X-ray enhancement is strongest at lower Galactic latitudes near l, b ≈ 25°, 25°, and gradually decreases in intensity toward the Galactic pole. Spectral observations with XMM-Newton and Suzaku indicate that the plasma is hotter than the surrounding medium with ${kT}\approx 0.3\,{\rm{keV}}$, is depleted in C, O, Mg, Ne, Fe, and enhanced in N (Willingale et al. 2003; Miller et al. 2008; Ursino et al. 2016). These measurements alone do not constrain the NPS distance, which makes it unclear whether the enhancement is associated with the Fermi bubbles.

Several independent methods place constraints on the NPS distance, but there is tension between the results. Sofue (2015) analyzed how the NPS X-ray emission varies with extinction near the Aquila Rift, and concluded that the NPS must be behind the rift and >1 kpc from the Sun. Puspitarini et al. (2014) compared three-dimensional models of the Milky Way's ISM and found no evidence for a low-density, hot cavity capable of producing the NPS X-ray emission within $\approx 200$ pc of the Sun. Alternatively, Sun et al. (2015) mapped the NPS in polarized radio emission and constrained the distance using maps of Faraday rotation measure and depth. They argue that the NPS emission at b ≳ 50° is likely within a few hundred parsecs of the Sun, while emission at lower latitudes can be either local or distant depending on the sign of the Milky Way's large-scale magnetic field.

These distance discrepancies are problematic because the NPS could be a region of compressed circumgalactic medium (CGM) material due to the Fermi bubbles' expansion, and thus a probe of their kinematics. This scenario was suggested long before the bubbles' discovery, and models of "biconical hypershells" expanding away from the Galactic center can reproduce the X-ray and radio morphologies of the NPS (Sofue 1977, 2000; Sofue et al. 2016). Alternatively, models of stellar winds and supernova remnants from the Sco–Cen OB association can also reproduce the X-ray and radio morphologies, in which case the emission occurs $\approx 200$ pc away from the Sun (Egger & Aschenbach 1995; Wolleben 2007). In this scenario, the NPS has no affiliation from the Fermi bubbles, so analyzing its emission would lead to inaccurate measured bubble properties.

In light of these uncertainties, we choose to exclude observations that are projected through the NPS in our model fitting analysis. This is a conservative approach that prevents us from fitting hot gas emission that may be unrelated to our emission model. Furthermore, the NPS appears to be a unique feature in terms of its location on the sky. There should be similarly enhanced X-ray-emitting regions in the other quadrants of the sky if the NPS were due to the Fermi bubbles, but these are not seen in the ROSAT maps. Thus, modeling the NPS emission could lead to biased and inaccurate results.

We exclude observations near the NPS in two square regions. Any observations within l = 20°–25°, b = 26°–34°, or l = 4°–20°, b = 40°–57° were excluded. This led to the exclusion of 27 observations from our model fitting procedure.

The final Suzaku sample with spectral fitting results can be seen in Table 1. We include the observation ID, the sight line in Galactic coordinates, the good XIS exposure time, and the oxygen line strengths with their 1σ uncertainties in Line Units (L.U.). The emission line measurements presented here are designed to have as little SWCX as possible, while containing only emission from astrophysical sources of interest (i.e., Galactic hot gas halo and Fermi bubbles). We also outlined data reduction, extraction, and cleaning procedures as similar as possible to the work by HS12, such that this sample and the XMM-Newton data are processed in a uniform way.

Our total data sample used in our astrophysical modeling combines the XMM-Newton and Suzaku measurements. There are 683 XMM-Newton measurements in total distributed across the sky, with 34 projected near the Fermi bubbles. The Suzaku data are exclusively projected near the Fermi bubbles, and there are 58 measurements included here. Figure 4 shows all-sky maps of the oxygen emission line strengths.

Zoom In Zoom Out Reset image size

The "local" emission source has been argued to include emission from both the LB and SWCX based on numerous shadowing experiments (Galeazzi et al. 2007; Koutroumpa et al. 2007, 2011; Smith et al. 2007) and studies of the ROSAT 1/4 keV band (Kuntz & Snowden 2000; Galeazzi et al. 2014). As discussed in Section 2.3, SWCX emission is difficult to predict or quantify; however, the flux-filtering techniques are designed to reduce its contribution to our measured line strengths. The physical properties of the LB have also been debated, with some studies arguing that the LB is volume-filled with ∼10⁶ K gas (Smith et al. 2007) and others arguing that the emission comes primarily from the bubble edges about 100–300 pc away (Lallement et al. 2003; Welsh & Shelton 2009). Regardless of these differences, our goal is to choose a parameterization that characterizes the emission from this source.

We parameterize the LB as a volume-filled plasma with a constant density and temperature and a size varying between 100 and 300 pc. This follows interpretation from Smith et al. (2007), who conducted SXRB modeling with Suzaku on the nearby molecular cloud MBM12. Under the assumption of a volume-filled plasma, these authors concluded that the LB has a temperature of $1.2\times {10}^{6}$ K and a density of 1–4 × ${10}^{-3}$cm⁻³. In our model, we fix the plasma temperature to this value and let the density, n_LB, be a free parameter.

While we include an LB emission source in our model for completeness, it is unlikely to have a significant impact on our results. The "local" plasma source is known to contribute more to the ROSAT 1/4 keV band than to the ROSAT 3/4 keV band (e.g., Snowden et al. 1997; Kuntz & Snowden 2000). This implies that the "local" emission source should produce more O vii than O viii in a given observation. Shadowing spectroscopic observations verify this, and show there is minimal ($\lesssim 0.5\,{\rm{L}}$.U.) O viii due to "local" sources (Koutroumpa et al. 2007, 2011; Smith et al. 2007). Furthermore, MB15 showed that this LB model effectively has no contribution to the O viii emission lines from the XMM-Newton sample discussed in Section 2.1. Our modeling work below focuses on O viii emission lines, so we do not believe that this LB parameterization will affect our results.

We assume that the Milky Way's "extended" hot gas plasma structure is dominated by a spherical, volume-filling halo of material extending to the virial radius, as opposed to the alternative assumption of an exponential disk morphology with scale height between 5 and 10 kpc. The latter structure is believed to form from supernovae in the disk (e.g., Norman & Ikeuchi 1989; Joung & Mac Low 2006; Hill et al. 2012) and can reproduce X-ray absorption and emission line strengths in several individual sight lines (Yao & Wang 2005, 2007; Yao et al. 2009; Hagihara et al. 2010). However, numerous studies have shown that a spherical, extended morphology due to shock-heated gas from the Milky Way's formation reproduces a multitude of observations (e.g., White & Frenk 1991; Cen & Ostriker 2006; Fukugita & Peebles 2006). These include ram-pressure stripping of dwarf galaxies (Blitz & Robishaw 2000; Grcevich & Putman 2009; Gatto et al. 2013), the pulsar dispersion measure toward the Large Magellanic Cloud (Anderson & Bregman 2010; Fang et al. 2013), and the aggregate properties of oxygen absorption and emission lines distributed in multiple sight lines across the sky (Bregman & Lloyd-Davies 2007; Gupta et al. 2012; Miller & Bregman 2013, 2015; Faerman et al. 2016). This distribution has been proven to reproduce most of the O viii emission line intensities from the XMM-Newton portion of the sample, thus justifying its use in this modeling work.

Our parameterized density distribution follows a spherical β-model, which assumes that the hot gas is approximately in hydrostatic equilibrium with the Milky Way's dark-matter potential well. The β-model has also been used to fit X-ray surface brightness profiles around early-type galaxies (e.g., O'Sullivan et al. 2003) and massive late-type galaxies (Anderson & Bregman 2011; Dai et al. 2012; Bogdán et al. 2013a, 2013b; Anderson et al. 2016). The model is defined as

$$\begin{eqnarray}&&n(r)\,=\,{n}_{o}{(1+{(r/{r}_{{\rm{c}}})}^{2})}^{-3\beta /2},\end{eqnarray} \tag{ 1 }$$

where r is the galactocentric radius, ${n}_{\circ }$ is the central density, r_c is the core radius ($\lesssim 5$ kpc), and β defines the slope (typically between 0.4 and 1.0). The previous modeling by MB15 was limited to using an approximate form of this model in the limit where $r\gg {r}_{{\rm{c}}}$, since they specifically did not include observations near the expected r_c. This resulted in constraints on a power-law density distribution:

$$\begin{eqnarray}&&n(r)\approx \displaystyle \frac{{n}_{o}{r}_{{\rm{c}}}^{3\beta }}{{r}^{3\beta }}.\end{eqnarray} \tag{ 2 }$$

The emission line sample in this study includes 33 sight lines that pass within 20° of the Galactic center, so we present model results assuming both distributions. The net effect of this will be for the power-law model to produce more halo emission for sight lines near the Galactic center than the usual β-model since the density continues to increase at small r instead of approach n_o for $r\lesssim {r}_{{\rm{c}}}$.

We assume the halo gas is isothermal with a temperature of log(T_halo) = 6.30, or ${T}_{\mathrm{halo}}\,=\,2\times {10}^{6}$ K. This temperature is characteristic of the Milky Way's virial temperature, and thus consistent with the picture in which the "extended" plasma is spherical and extended to r_vir. This temperature is also constrained by observations. Henley & Shelton (2013) provide the strongest observational constraints on the plasma temperature, as they fit SXRB spectra for 110 high-latitude ($| b| \gt 30$°) sight lines from the HS12 sample. They fit the spectra with thermal APEC plasma models and found a narrow range of plasma temperatures (median and interquartile range of $(2.22\pm 0.63)\times {10}^{6}$ K). These results suggest that the plasma is nearly isothermal, thus validating our assumption.

Our Fermi bubble structure includes two components: a volume-filled ellipsoid and a shell surrounding the ellipsoid. These components are designed to parameterize the outer regions of galactic outflows, where the volume-filled structure includes hot, shocked wind material, and the shell includes shocked ISM/CGM material (see Section 5.2.1 for an overview of galactic wind morphology). Hard X-ray emission bounds the bubbles at low Galactic latitudes, verifying that there is a distinct shell structure of thermal gas surrounding the bubbles (Bland-Hawthorn & Cohen 2003; Su et al. 2010). The flat gamma-ray intensity distribution on the sky indicates that non-thermal particles fill the bubbles in a quasi-spherical volume, and we include a thermal gas component in this region.

We define the bubble volume as a three-dimensional ellipsoid designed to match the bubbles' projected gamma-ray edge on the sky. Each bubble (positive and negative Galactic latitudes) is centered at $| z|$ = 5 kpc away from the Galactic plane, has a semimajor axis of 5 kpc, and has both minor axes set to 3 kpc. We also tilt each bubble $5$° toward negative longitudes to match the slight asymmetry observed in the bubble shape. Figure 5 shows this bubble volume in physical and projected coordinates.

Zoom In Zoom Out Reset image size

The shell volume is defined in the same way as the bubble volume, but with a thickness of $\approx 1\,\mathrm{kpc}$ away from the bubble surface. This implies that the shell ellipsoids are also centered at $| z|$ = 5 kpc from the Galactic plane and have semimajor axes of 6 kpc and minor axes of 4 kpc. The region inside this surface but outside the bubble surface is considered to be the shell region. We note that this parameterization is different from the modeling work presented by K15, who considered Fermi bubble emission from only two angled shells (one for each bubble) with inner and outer radii of 3 and 5 kpc. However, Figure 5 indicates that our bubble volume is consistent with the projected bubble outline, and the expected galactic wind morphology suggests that there should be at least two distinct outflow regions we can observe (the shocked wind and shocked ISM/CGM). Therefore, we feel that our choice of bubble volume is reasonable given the observational constraints available at this time.

We assume that the bubble and shell components each have constant electron densities, defined as n_FB and n_shell. This parameterization is useful since it is simple, yet still allows us to constrain the average thermal gas densities. Simulations suggest that the bubbles have some thermal gas substructure (e.g., Yang et al. 2012), and we did explore more sophisticated models with density gradients away from the Galactic plane or from the bubble edges. However, the data did not provide statistically significant constraints with these profiles (any gradient parameter was consistent with a constant-density profile within the 1σ uncertainties). These constant-density models should be considered a valuable first step when analyzing the bubbles' thermal gas structure.

The bubble and shell are likely hotter than the surrounding medium, so we assume that each component has a characteristic temperature $\geqslant 2\times {10}^{6}$ K. Like the bubble and shell densities, each component has a constant temperature (T_FB and T_shell). However, these temperatures are each initially fixed to $3\times {10}^{6}$ K during the model fitting process. The temperatures are not free parameters in our models because the calculated line intensity scales with density and temperature as $I\propto {n}^{2}\epsilon (T)$, where has a temperature dependence. Since we explicitly model a sample of O viii emission line intensities, the density and temperature parameters would be degenerate with each other. Section 4 discusses how we constrain the bubble and shell temperatures by looking at the distribution of O viii/O vii line ratios for different assumed temperatures.

Our modeling also implicitly assumes that the Fermi bubble and shell plasmas have solar metallicities. This implies that our bubble density parameters are degenerate with the assumed metallicity since the plasma emission measure scales linearly with metallicity. We make this choice because there are no direct observational constraints on the plasma metallicity inside the bubbles. The SXRB modeling from Kataoka et al. (2013) advocates for a sub-solar bubble metallicity of $Z\approx 0.2\,{Z}_{\odot }$, but this value is weakly constrained due to photon statistics, and the spectral fits represent the emission measure-weighted spectrum from the bubbles and Galactic hot halo (see discussion in Section 5.4.1). Thus, it is difficult to interpret whether the bubbles' plasma is enriched or sub-solar. The bubbles' abundance ratios can be diagnostics for how they formed (Inoue et al. 2015), but a detailed abundance analysis is beyond the scope of this work.

Calculating model line intensities depends on the density and temperature profile along each line of sight. For any given Galactic coordinate ($l,b$), we divide the line of sight into cells extending to the virial radius (r_vir = 250 kpc). Each cell position along the line of sight (s) is converted to Galactic coordinates (R, z, r) by the standard equations:

$$\begin{eqnarray}&&{R}^{2}={R}_{\odot }^{2}+{s}^{2}\cos {(b)}^{2}-2{{sR}}_{\odot }\cos (b)\cos (l),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{z}^{2}={s}^{2}\sin {(b)}^{2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{r}^{2}={R}^{2}+{z}^{2},\end{eqnarray} \tag{ 5 }$$

where R_⊙ = 8.5 kpc is the Sun's distance from the Galactic center. We assign a density and temperature to each cell based on its set of Galactic coordinates and the assumed model parameters. The hot halo profile described in Section 3.2 sets the density and temperature for cells outside the shell volume. The parameters n_shell and T_shell set the density and temperature for cells within the shell volume, while n_FB and T_FB set the density and temperature for cells within the bubble volume. Therefore, sight lines not passing through the bubbles include only halo emission, while sight lines passing through the bubbles include emission from the hot gas halo, bubble, and shell.

We assume an optically thin plasma in collisional ionization equilibrium to calculate all line intensities. Given a line-of-sight density and temperature profile, the model line intensity is defined as

$$\begin{eqnarray}&&I(l,b)\,=\,\displaystyle \frac{1}{4\pi }\int {n}_{{\rm{e}}}{(s)}^{2}\epsilon (T(s)){ds},\end{eqnarray} \tag{ 6 }$$

where n_e(s) is the line-of-sight electron density, T(s) is the line-of-sight temperature profile, and $\epsilon \,(T)$ is the volumetric line emissivity for a thermal APEC plasma. We use AtomDB version 2.0.2 for all line emissivities (Foster et al. 2012), and characteristic values for O viii (in photons cm³ s⁻¹) are $\epsilon ({T}_{\mathrm{halo}})\,=\,1.45\,\times \,{10}^{-15}$ and $\epsilon (3\times {10}^{6}\ {\rm{K}})\,=\,3.84\,\times \,{10}^{-15}$. Although believed to be minimal, the intensity contribution from the LB is defined as

$$\begin{eqnarray}&&{I}_{\mathrm{LB}}(l,b)\,=\,\displaystyle \frac{{n}_{\mathrm{LB}}^{2}L(l,b)\epsilon ({T}_{\mathrm{LB}})}{4\pi },\end{eqnarray} \tag{ 7 }$$

where $L(l,b)$ defines the LB path length (≈100–300 pc; see Lallement et al. 2003 and MB15). The total line intensity is thus defined as

$$\begin{eqnarray}&&{I}_{\mathrm{total}}(l,b)\,=\,{I}_{\mathrm{LB}}+{e}^{-\sigma {N}_{H{\rm{I}}}}({I}_{\mathrm{halo}}+{I}_{\mathrm{FB}}+{I}_{\mathrm{shell}}),\end{eqnarray} \tag{ 8 }$$

where the exponential term accounts for attenuation due to neutral hydrogen in the disk, N_{H i} is the same neutral hydrogen column assumed for each sight line in the spectral fitting procedure, and σ is the H i absorption cross section (Balucinska-Church & McCammon 1992; Yan et al. 1998). Thus, our model line intensities are comparable to the total observed line intensities.

The observed O viii/O vii ratios in our sample can be used as crude temperature diagnostics. If the observed emission lines come from a single, cospatial plasma source, Equation (6) indicates that the O viii/O vii ratio is a direct temperature diagnostic because ${I}_{{\rm{O}}{\rm{VIII}}}/{I}_{{\rm{O}}{\rm{VII}}}\propto {n}^{2}{\epsilon }_{{\rm{O}}{\rm{VIII}}}(T)/{n}^{2}{\epsilon }_{{\rm{O}}{\rm{VII}}}(T)\,=\,{\epsilon }_{{\rm{O}}{\rm{VIII}}}(T)/{\epsilon }_{{\rm{O}}{\rm{VII}}}$(T). The observations from SXRB spectra are more complicated since we know that multiple plasma sources exist along each line of sight. This implies that the total observed O viii/O vii line ratio probes the emission measure-weighted temperature due to the various plasma sources. However, we discussed in Section 3 how the LB is believed to produce little O viii emission with a variable amount of O vii emission, and the hot halo plasma is believed to be nearly isothermal at $\approx 2\times {10}^{6}$ K. The expected O viii/O vii line ratio for a thermal plasma at this temperature is $\approx 0.25$, so the observed line ratios in our sample would be $\lesssim 0.25$ if they included only emission from the LB and hot halo.

We explore this idea by examining the O viii/O vii distribution on the sky from our total observation sample. Figure 6 shows our line intensity ratios on the sky. Inspecting the sky projection alone suggests that the line intensity ratios are systematically higher for sight lines that pass through or near the Fermi bubbles (≈0.5) than for those farther away from the Galactic center (≈0.2). To quantify this, we bin the sight lines on the sky and calculate the median and interquartile range for the line ratios in each bin. The bin edges are defined as ellipses in $l,b$ space, where the first bin includes sight lines that pass though the Fermi bubbles and subsequent bins include sight lines extending farther into the halo (see dotted lines in Figure 6). Figure 6 shows the line ratio median and interquartile range for observations in each bin. These results clearly show that the line ratios are systematically higher for sight lines in the first bin, and are also higher than the characteristic ratio expected if the observations included just LB and hot gas halo emission (gray shaded band in Figure 6).

Zoom In Zoom Out Reset image size

These systematically larger line ratios near the Galactic center indicate the presence of hotter gas than the ambient $2\times {10}^{6}$ K plasma. This interpretation is model-independent and builds upon the fact that we know the Fermi bubbles occupy a significant volume above and below the Galactic center. While this is a useful result that relies only on observations, the observed line ratios do not encode the bubbles' detailed temperature structure due to additional emission from the LB and hot gas halo. Nevertheless, this result motivates the modeling work below and validates the assumption that the bubbles contain gas hotter than $2\times {10}^{6}$ K.

As a preliminary test, we explore an emission model including only contributions from the LB and hot gas halo. This model assumes that the bubble and shell volumes contribute no line emission, or equivalently that ${n}_{\mathrm{FB}}\,=\,{n}_{\mathrm{shell}}\,=\,0$. For the other emission components, we assume a parametric model distribution from MB15. This includes an LB density of ${n}_{\mathrm{LB}}\,=\,4\times {10}^{-3}$ cm⁻³ and a hot gas density profile described by Equation (2) with ${n}_{o}{r}_{{\rm{c}}}^{3\beta }$ = $1.35\times {10}^{-2}$ cm⁻³ kpc${}^{3\beta }$ and β = 0.5. This model likely overestimates any halo emission since it assumes a power law all the way to the Galactic center, as opposed to having a flat core density. We calculate model O viii emission line intensities for this limiting case and compute the residual emission defined as $({I}_{\mathrm{observed}}-{I}_{\mathrm{model}})/{I}_{\mathrm{error}}$. Figure 7 shows how the residual emission varies on the sky, with a particular emphasis on the strong ($\gtrsim 3\sigma$) positive residuals near the Fermi bubbles. We interpret these residuals as missing emission due to the Fermi bubbles, which motivates the modeling procedure outlined below.

Zoom In Zoom Out Reset image size

The goal of our modeling procedure is to find the density model that is most consistent with our observed data set, including contributions from the Fermi bubble and shell components. We quantify this consistency with the model ${\chi }^{2}$ or likelihood (${ \mathcal L }\propto$ exp $(-{\chi }^{2}/2)$). We use the publicly available Markov chain Monte Carlo (MCMC) Python package emcee (Foreman-Mackey et al. 2013) to explore our model parameter space and find the parameters that minimize the model ${\chi }^{2}$, or maximize the model ln(${ \mathcal L }$). The output chains for each model parameter are treated as marginalized posterior probability distributions. We define "best-fit" parameters as the median values for each binned chain distribution, which yields identical results to the Gaussian-fitting procedure outlined by MB15, assuming the distributions are approximately Gaussian. Thus, these best-fit parameters maximize the model likelihood, given the data.

We considered several different model parameterizations in our model fitting process. These included hot gas halo density models described by either a power law (Equation (2)) with two free parameters (the normalization and β) or a full β-model (Equation (1)) with r_c and β left to vary. We did not let n_o vary independently in this model since the previous modeling work from MB15 effectively constrained the halo normalization parameter ${n}_{o}{r}_{{\rm{c}}}^{3\beta }$. Our modeling procedure keeps this quantity fixed to ${n}_{o}{r}_{{\rm{c}}}^{3\beta }\,=\,1.35\times {10}^{-2}$ cm⁻³ kpc${}^{3\beta }$, while letting the core radius vary as the free parameter. We also experimented with fixing the hot gas halo profile with the fit values from MB15 or with r_c = 3 kpc, but we found that this made little difference in the best-fit parameters for either the halo density profile or the bubble/shell densities.

Table 2 summarizes our best-fit model parameters, including 1σ uncertainties encompassing the 68% probability ranges from the posterior probability distributions. There are several trends to note from these results. The LB density parameter is consistent with zero, validating the assumption that the LB contributes little emission to the O viii data. The halo density profile results are consistent with those reported in MB15 when considering the same power-law density parameterization (${n}_{o}{r}_{{\rm{c}}}^{3\beta }\,=\,1.35\times {10}^{-2}$ cm⁻³ kpc${}^{3\beta }$, β = 0.5). This implies that the hot gas density profile constraints are not biased due to observations near the Fermi bubbles. We also find characteristic best-fit core radii of 2–3 kpc, which is expected. The parameters of interest, n_FB and n_shell, have characteristic best-fit densities of (5–8) × ${10}^{-4}$ cm⁻³, assuming a temperature of log(T) = 6.50. The inferred densities are lower if we assume a power-law model for the halo gas density than if we assume a β-model. We expect to see this trend since the power-law model produces more halo emission near the Galactic center than a β-model with a core radius/density, thus resulting in less Fermi bubble/shell emission being required to produce the observed emission. After weighing these effects, we define our fiducial model to be one with a β-model and r_c fixed to 3 kpc. This results in best-fit parameters of ${n}_{\mathrm{FB}}=(7.2\pm 0.2)\times {10}^{-4}$ cm⁻³, and ${n}_{\mathrm{shell}}\,=\,(7.7\pm 0.2)\times {10}^{-4}$ cm⁻³. Figure 8 shows the marginalized posterior probability distributions and contour plots from our MCMC analysis assuming this parametric model (generated using the Python code corner.py; Foreman-Mackey et al. 2016).

Zoom In Zoom Out Reset image size

Table 2.  MCMC Fitting Results

nLBa	${n}_{o}{r}_{{\rm{c}}}^{3\beta }$	no	rc	β	nFBb	nshellb	${\chi }^{2}$ (dof)
(10−3 cm−3)	(10−2 cm−3 kpc${}^{3\beta }$)	(10−3 cm−3)	(kpc)		(10−4 cm−3)	(10−4 cm−3)
<5.81	1.01 ± 0.06	⋯	⋯	0.45 ± 0.01	6.70 ± 0.19	6.25 ± 0.30	2683 (736)
<5.42	1.35 (fixed)	4.47	2.12 ± 0.22	0.49 ± 0.01	6.67 ± 0.19	6.20 ± 0.29	2669 (736)
3.83 (fixed)	1.35 (fixed)	⋯	⋯	0.50 (fixed)	6.61 ± 0.18	6.01 ± 0.27	2716 (739)
3.83 (fixed)	1.35 (fixed)	2.60 (fixed)	3 (fixed)	0.50 (fixed)	7.17 ± 0.17	7.48 ± 0.22	2783 (739)

Notes. Unless noted otherwise, all best-fit values are defined as the most likely parameter values from the MCMC marginalized posterior probability distributions. The uncertainties encompass the 68% probability region relative to the maximum likelihood value for each parameter.

^aThe Local Bubble (LB) produces minimal O viii emission in our model. With this in mind, we report either the fixed value for n_LB or the 2σ upper limit from our MCMC analysis. ^bThese best-fit densities assume that each component has a temperature of log(${T}_{\mathrm{FB},\mathrm{shell}}$) = 6.50 when calculating line intensities. See Table 3 for densities with higher assumed temperatures.

Download table as:  ASCIITypeset image

We also explore models with either the bubble or shell distributions to determine each component's significance in our model fitting procedure. A bubble-only model is equivalent to setting the shell thickness to 0 kpc and not including n_shell as a model parameter. A shell-only model is the same as our original emission model (r_c fixed to 3 kpc) but with n_FB fixed to zero. When we refit the data with these models, we find that each density component increases to compensate for the lack of emission from the other component. For example, n_FB increased to $7.7\times {10}^{-4}$ cm⁻³ in the bubble-only model, and n_shell increased to $10.0\times {10}^{-4}$ cm⁻³ in the shell-only model. These best-fit models lead to changes in the overall fit quality, where our initial best-fit ${\chi }_{r}^{2}$ (dof) = 2783 (739). The bubble-only model leads to a marginal improvement in the overall fit quality (${\chi }_{r,\mathrm{bubble} \mbox{-} \mathrm{only}}^{2}$ (dof) = 2741 (740)), while the shell-only model leads to a worse quality of fit (${\chi }_{r,\mathrm{shell} \mbox{-} \mathrm{only}}^{2}$ (dof) = 3241 (740)). This implies that the volume-filled structure is more important than the shell structure, although we point out that this exercise fixes every other component of the emission model (halo emission, bubble/shell geometry, etc.). Thus, we still assume that the bubble and shell structures are each present in our discussion and temperature analysis.

In order to constrain the bubble and shell temperatures, we compare best-fit model line ratios for different temperature distributions with the observed O viii/O vii line ratios near the Fermi bubbles. To do this, we change the bubble and shell temperatures while keeping the product ${n}^{2}{\epsilon }_{{\rm{O}}{\rm{VIII}}}(T)$ fixed from the best-fit model results. This fixes the O viii emission coming from the bubble and shell, but changes the model O vii emission because ${\epsilon }_{{\rm{O}}{\rm{VII}}}(T)$ decreases faster than ${\epsilon }_{{\rm{O}}{\rm{VIII}}}(T)$ with increasing temperature. Thus, increasing the assumed temperature leads to an increase in the model line ratios, an increase in the inferred best-fit densities (${\epsilon }_{{\rm{O}}{\rm{VIII}}}$decreases for $T\gt 3\times {10}^{6}$ K), and a constant contribution to the O viii emission.

A model temperature distribution with log(T_FB) = 6.60 and log(T_shell) = 6.70 leads to a model line ratio distribution most consistent with the observed line ratios near the Fermi bubbles. This changes the inferred best-fit densities to ${n}_{\mathrm{FB}}\,=\,8.2\times {10}^{-4}$ cm⁻³ and ${n}_{\mathrm{shell}}\,=\,1.0\times {10}^{-3}$ cm⁻³ in order to keep the product ${n}^{2}{\epsilon }_{{\rm{O}}{\rm{VIII}}}(T)$ fixed for each component. Figure 9 shows histograms of the observed and new best-fit model line ratios for sight lines that pass within $\approx 5$° of the projected bubble edge (∼100 sight lines). These densities and temperatures produce an O viii/O vii ratio median and interquartile range of 0.52 (0.41–0.60), consistent with the observed median and interquartile range of 0.49 (0.38–0.62). We treat these densities and temperatures as the characteristic physical properties for the bubble and shell components in our subsequent analysis.

Zoom In Zoom Out Reset image size

We discuss our inferred bubble densities and temperatures in this section, and compare them to the assumed ambient structure. Overall, our constraints indicate that the bubbles are hotter and overpressurized compared to the surrounding medium, consistent with previous observations of the bubbles (Su et al. 2010; Kataoka et al. 2015). Figure 10 shows our best-fit model as two-dimensional maps of density, temperature, and pressure projected at the Galactic center. This visualizes the comparison with the surrounding medium that we discuss in the rest of the section.

Zoom In Zoom Out Reset image size

The derived structure indicates that the best-fit densities and temperatures for the bubble and shell components are nearly identical to each other. This suggests that there may not be two distinct outflow regions, which contradicts the global morphology predicted from Galactic outflow simulations (e.g., Yang et al. 2012; Sarkar et al. 2015; Sofue et al. 2016). This inference could be due to our assumed geometry, in terms of both the volume-filled shape and the shell thickness. Our emission model effectively constrains each component's emission measure, so changes in the bubble/shell path lengths along different sight lines could change the inferred densities. Additional substructure could also be present inside the volume-filled component, which might explain why we infer a relatively high density permeating the entire bubble volume. Modeling this substructure is beyond the scope of this work, but our modeling results still provide valuable constraints on the average bubble and shell properties.

The bubble and shell densities have characteristic values of $\sim {10}^{-3}$ cm⁻³, which are comparable to the surrounding medium at low z. Including a core radius for the hot gas halo of 3 kpc in our fiducial model implies a core density of $2.6\times {10}^{-3}$ cm⁻³, assuming a fixed power-law normalization of ${n}_{o}{r}_{{\rm{c}}}^{3\beta }\,=\,1.35\times {10}^{-2}$ cm⁻³ kpc${}^{3\beta }$. This suggests ${n}_{\mathrm{shell}}\sim {n}_{\mathrm{halo}}$ within $| z| \lesssim 5\,{\rm{kpc}}$. The hot gas halo density decreases by about a factor of 6 between r = 1 and 10 kpc, meaning that our bubble and shell densities are larger than the surrounding CGM density farther away from the Galactic plane. We also note that ${n}_{\mathrm{FB}}\approx {n}_{\mathrm{shell}}$ from our model fitting results, making it difficult to distinguish between volume-filling emission and limb-brightened emission. This might be due to our choice to parameterize the structures with constant densities and temperatures, but our constraints still probe the average densities associated with the bubbles.

Our inferred bubble and shell temperatures of log(${T}_{\mathrm{FB},\mathrm{shell}}$) = 6.60–6.70 are hotter than the surrounding medium (log(T_halo) = 6.30). This is broadly consistent with the bubbles injecting enough energy to shock-heat the surrounding medium, although simulations predict a wide range of shock strengths and bubble temperatures as high as $\sim {10}^{8}$ K (e.g., Guo & Mathews 2012; Yang et al. 2012). It is possible that the bubbles contain gas at this high temperature, but plasma at this temperature would not produce observable O viii emission. Our modeling results indicate that on average the bubble and shell are hotter than the surrounding medium, but still at low enough temperatures to produce observable signatures in the data.

Combining the density and temperature constraints indicates that the bubbles' are overpressurized compared to the surrounding medium. Our best-fit bubble and shell parameters indicate thermal gas pressures of ${P}_{\mathrm{FB}}\,=\,4.5\times {10}^{-13}$ dyn cm⁻², ${P}_{\mathrm{shell}}\,=\,$ $6.9\times {10}^{-13}$ dyn cm⁻², or P/k ≈ 3000–5000 cm⁻³ K. The surrounding thermal gas pressure varies with r and z due to the decreasing density profile, with characteristic values of $\approx 5000$ cm⁻³ K near the Galactic center and $\approx 1000$ cm⁻³ K at r = 10 kpc. In this picture, the bubbles are in approximate pressure equilibrium at lower z, but become overpressurized with increasing height away from the Galactic plane. These estimates are also a lower limit to how overpressurized the bubbles actually are, because they do not account for non-thermal or magnetic pressure contributions. Nevertheless, these constraints indicate that the bubbles are generally overpressurized, and thus expanding into the surrounding medium.

We use these quantities to infer a characteristic shock strength and instantaneous expansion velocity associated with the bubbles. The classic treatment of shocks propagating through the ISM yields specific pre- and post-shock jump conditions for the gas density, temperature, and pressure given a shock expansion velocity (e.g., Shull & Draine 1987). The Fermi bubbles' expansion is more complex than this traditional treatment since they do not appear to be spherical, and they are presumably expanding into a medium with varying density. For example, the ratio between n_shell (treated as post-shocked material) and the ambient halo gas density along the shell edge (treated as pre-shock material) ranges between $\approx 0.5$ near the Galactic center and $\approx 3$ at the maximum bubble height. On the other hand, our choice to parameterize the hot gas halo and shell with constant temperatures allows us to use the temperature ratio as a shock strength diagnostic similar to K15. Assuming a monotonic gas (γ = 5/3), an ambient gas sound speed of c_s = 212 km s⁻¹, ${T}_{\mathrm{shell}}\,=\,5\times {10}^{6}$ K, and ${T}_{\mathrm{halo}}\,=\,2\times {10}^{6}$ K, our temperature ratio implies a fiducial Mach number and corresponding expansion velocity of ${ \mathcal M }\,=\,2.3$ and v_exp = 490 km s⁻¹. The uncertainty in T_shell expands these constraints to ${ \mathcal M }\,=\,$ 1.9–3.4 and v_exp = 413–720 km s⁻¹. These shock parameters are broadly consistent with the range of density and pressure ratios we estimate, indicating that these constraints probe the bubbles' current expansion rate into the surrounding medium.

5.2.1. The Bubbles as a Confined Galactic Wind

We treat the bubbles in the framework of a continuous galactic outflow/superbubble with self-similar Sedov–Taylor solutions (e.g., Castor et al. 1975; Weaver et al. 1977; Mac Low & McCray 1988; Veilleux et al. 2005). The outflow morphology consists of five zones (from closer to farther from the outflow origin): the energy injection zone, a free-flowing outflow, shocked wind material, a shell of shocked ISM/CGM material, and the ambient ISM/CGM. Our model constraints probe the last three zones since we do not model observations in the inner $\approx 1\,{\rm{kpc}}$ from the Galactic center. The Sedov–Taylor solutions for this type of outflow relate the ambient density, bubble age, bubble size, expansion velocity, and average energy injection to each other. Assuming that the outflow is still in the energy-conserving phase (cooling time greater than the bubble age), the relations between these quantities are as follows:

$$\begin{eqnarray}&&{t}_{\mathrm{age}}=11.8\,\mathrm{Myr}\ \left(\displaystyle \frac{r}{10\,\mathrm{kpc}}\right){\left(\displaystyle \frac{v}{500\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-1},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\xi \dot{E} & = & 3.7\times {10}^{42}\,\mathrm{erg}\ {{\rm{s}}}^{-1}\\ & & \times \left(\displaystyle \frac{{n}_{o}}{{10}^{-3}\ {\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{r}{10\mathrm{kpc}}\right)}^{2}{\left(\displaystyle \frac{{v}_{\exp }}{500\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{3},\end{eqnarray} \tag{ 10 }$$

where n_o is the ambient density, r is the bubble radius, v_exp is the expansion velocity, t_age is the bubble age, $\dot{E}$ is the energy injection rate, and ξ is the thermalization efficiency of the mechanical energy. This thermalization efficiency is believed to vary with environment, but is estimated to be ≳10% in average galaxies with a typical assumed value of 0.3 (e.g., Larson 1974; Wada & Norman 2001; Melioli & de Gouveia Dal Pino 2004).

Our modeling results constrain three of the five variables in these equations. As discussed above, a hot gas halo model with r_c = 3 kpc and fixed power-law normalization of $1.35\times {10}^{-2}$ cm⁻³ kpc${}^{3\beta }$ results in an ambient core density of ${n}_{o}=2.6\times {10}^{-3}$ cm⁻³. The constraints on the bubble and shell temperatures suggest an outflow velocity of 490 km s⁻¹. The bubble size is not trivial to estimate in this framework, where the outflow is typically treated as a spherical shell with radius r centered on the injection source. The bubbles' shape is more complex than this since two lobes exist on each side of the Galactic plane. We use a characteristic bubble size defined as the geometric mean of the three ellipsoidal axes, resulting in r = 3.6 kpc.

Given these constrained values, we estimate the bubbles' age and average mechanical energy injection rate. The bubbles' dynamical timescale, ${t}_{\mathrm{dyn}}\,=\,d/{v}_{\exp }$, is a crude age estimate that does not incorporate the bubble environment or energy source. For v_exp = 490 km s⁻¹, the dynamical timescale for the full bubble height is ${t}_{\mathrm{dyn},h}$ = 10 kpc/490 km s⁻¹ = 20.0 Myr, and the dynamical timescale for half the bubble width is ${t}_{\mathrm{dyn},w}$ = 3 kpc/490 km s⁻¹ = 6.0 Myr. The superbubble model calculation (Equation (9)) is a refined age estimate, where we find t_age = 4.3 Myr for r defined as the geometric mean above and v_exp = 490 km s⁻¹. We also infer a combined energy injection for both bubbles ($2\times \xi \times \dot{E}$) of $2.3\times {10}^{42}$ erg s⁻¹ using Equation (10). Accounting for the uncertainty in v_exp leads to a characteristic age range of ≈3–5 Myr and energy injection rate of ≈(1–7) × 10⁴² erg s⁻¹. We compare this characteristic age and energy injection rate with possible bubble formation mechanisms.

5.2.2. Origin from Sgr A* Accretion

One suggested bubble formation mechanism has been a past accretion event onto Sgr A*, resulting in an AGN episode in the Milky Way. Sgr A* has an estimated mass of $4\times {10}^{6}\,{M}_{\odot }$ (Schödel et al. 2002; Ghez et al. 2003, 2008; Gillessen et al. 2009a, 2009b; Meyer et al. 2012), which is capable of producing significant amounts of energy during an accretion episode. We also know that accretion onto supermassive black holes can produce galactic outflows with significant energy injection rates and morphologies similar to the observed Fermi bubbles (e.g., McNamara & Nulsen 2007; Yuan & Narayan 2014). Here, we consider observations of Sgr A* and its possible accretion history, and compare the expected energetics with our modeling constraints.

Sgr A* is currently in a quiescent state with a bolometric luminosity of ${L}_{\mathrm{bol}}\sim {10}^{36}$ erg ${{\rm{s}}}^{-1}\sim 2\times {10}^{-9}\ {L}_{\mathrm{Edd}}$ (e.g., Yuan & Narayan 2014). Our proximity to Sgr A* allows for a combination of techniques to estimate current mass accretion rates. Chandra's resolution is comparable to the Sgr A* Bondi radius, and has constrained the Bondi accretion rate to be $\sim {10}^{-5}\,{M}_{\odot }$ yr⁻¹ (Baganoff et al. 2003). Polarized radio emission constrains the accretion rate near the event horizon, with limits being between $\gt 2\times {10}^{-9}\,{M}_{\odot }$ yr⁻¹ and $\lt 2\times {10}^{-7}\,{M}_{\odot }$ yr⁻¹ depending on the magnetic field orientation (e.g., Marrone et al. 2007). While this is a significant uncertainty in the current mass accretion rate, the consensus is that Sgr A* is accreting well below its Eddington rate, and has been a well-modeled source for radiatively inefficient accretion flows (RIAFs).

There are a number of observational indications that Sgr A* has been more active in the past (Totani 2006). Mou et al. (2014) summarizes these lines of evidence, which include: a higher Sgr A* luminosity is required to produce fluorescent iron emission and reflection nebulae seen in several nearby molecular clouds (Koyama et al. 1996; Murakami et al. 2000, 2001a, 2001b), there exists an ionized halo of material surrounding Sgr A* (Maeda et al. 2002), there are dynamic features indicating an outflow near the Galactic center in the form of the Galactic Center Lobe (Bland-Hawthorn & Cohen 2003) and the Expanding Molecular Ring (Kaifu et al. 1972; Scoville 1972), excess Hα emission seen in the Magellanic Stream (Bland-Hawthorn et al. 2013), and possibly the Fermi bubbles themselves. The RIAF modeling from Totani (2006) argues that Sgr A* should have had an accretion rate $\sim {10}^{3}$–10⁴ times larger than its current accretion rate over the past ∼10 Myr to reproduce these observations. This introduces additional scatter in the inferred past Sgr A* accretion rate, but motivates the assumption that Sgr A* has injected energy into the surrounding medium through an accretion event.

We estimate an energy injection rate due to past Sgr A* accretion and compare with our constrained energy input rate. The mechanical energy injection rate from black hole accretion is tied to the accretion power by the following relation:

$$\begin{eqnarray}{\dot{E}}_{\mathrm{BH}} & = & \varepsilon {\dot{M}}_{\mathrm{acc}}{c}^{2}\\ & = & 5.7\,\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,\left(\displaystyle \frac{\varepsilon }{0.1}\right)\left(\frac{{\dot{M}}_{\mathrm{acc}}}{{M}_{\odot }\ {\mathrm{yr}}^{-1}}\right),\end{eqnarray} \tag{ 11 }$$

where ${\dot{E}}_{\mathrm{BH}}$ is the mechanical energy injection rate, ${\dot{M}}_{\mathrm{acc}}$ is the accretion rate near the event horizon, and is the mechanical energy injection rate efficiency. If we assume a past accretion rate of ${10}^{-3}\,{M}_{\odot }$ yr⁻¹ (near the high end of the values discussed above), we find that ${\dot{E}}_{\mathrm{BH}}$ equals our inferred vale of $2.3\times {10}^{42}$ erg s⁻¹ for $\epsilon \approx 0.05$. This efficiency is larger than the typical values inferred from simulations (10⁻⁴–${10}^{-3};$Yuan et al. 2015), but this mechanical efficiency is often treated as a free parameter in simulations. We also point out that the required efficiency is less than one, indicating that this analysis does not violate energy conservation constraints. Thus, it is plausible that a past accretion episode onto Sgr A* could have produced enough energy to match our energy injection rate constraints.

The bubble age indicates that this Sgr A* accretion episode had a shorter active period than the typical AGN duty cycle. Studies constrain the AGN duty cycle by either comparing black hole mass functions (inferred from the M_BH–σ relation) to AGN luminosity functions at different redshifts (e.g., Yu & Tremaine 2002; Shankar et al. 2004; Hao et al. 2005; Schawinski et al. 2010) or through analytic models of black hole accretion (e.g., Hopkins & Hernquist 2006; Shankar et al. 2009). These techniques suggest that black holes with masses of $\sim {10}^{6}$ M_⊙ should have active periods of $\sim {10}^{8}\,{\rm{yr}}$ at z = 0, or ∼1% of a Hubble time. Our Fermi bubble age estimate is an upper limit to the active Sgr A* accretion time, and our constraint of 4.3 Myr is much smaller than the inferred duty cycle from AGN populations. One possible explanation for this discrepancy is that the Fermi bubble outburst was one of several past accretion events in the Galactic center. Our constraints imply that the Fermi bubbles are due to a relatively weak AGN event, and it is possible that multiple Sgr A* accretion events of comparable or lower energy have occurred over the past $\sim {10}^{8}\,{\rm{years}}$. Our results are also consistent with the overall decrease in AGN activity since $z\sim 2$ (e.g., Hopkins et al. 2007), as opposed to a prolonged Sgr A* accretion/growth phase.

5.2.3. Origin from Nuclear Star Formation

Numerous studies also suggest that the Fermi bubbles formed from a period of enhanced star formation activity near the Galactic center. The Galactic center hosts several young stellar clusters with ages ranging between 5 and 20 Myr and accounting for $\sim 5\times {10}^{5}\,{M}_{\odot }$ of material. The massive stars in these clusters could have generated a galactic-scale outflow due to stellar winds and type-II supernova explosions (e.g., Leitherer et al. 1999). Here, we compare the expected energy output from past star formation near the Galactic center to our energy injection rate constraints.

The Galactic center star formation history is complex and difficult to measure, but several studies argue for an average star formation rate (SFR) of $\approx 0.05\,{M}_{\odot }$ yr⁻¹ over the past ∼10 Myr. Crocker (2012) reviews these studies, most of which utilize Spitzer observations of young stellar objects within the inner ∼500 pc from the Galactic center. For example, Yusef-Zadeh et al. (2009) conducted a census of these objects using the Infrared Array Camera and Multiband Imaging Photometer on board Spitzer, and concluded the average SFR has been between 0.04 and 0.08 M_⊙ yr⁻¹ over longer timescales (∼10 Gyr). Immer et al. (2012) performed a similar analysis using data from the Spitzer Infrared Spectrogaph, and argue for an average SFR of $\approx 0.08\,{M}_{\odot }$ yr⁻¹ over the past ∼1 Myr. Others estimate the SFR to be ≈0.01–0.02 M_⊙ yr⁻¹ by counting the mass in young star clusters and dividing that by estimates for the period of star formation (Figer et al. 2004; Mauerhan et al. 2010). Thus, it appears that a characteristic SFR of $\approx 0.05\,{M}_{\odot }$ yr⁻¹ over the past ∼10 Myr is a reasonable assumption.

Similar to our argument concerning the black hole accretion energy above, we estimate the energy injection rate due to star formation in the Galactic center to compare with our constrained energy input rate. Assuming a Kroupa initial mass function (Kroupa 2001), and 10⁵¹ erg of mechanical energy input from a type-II supernova, the mechanical energy from type-II supernovae is related to the SFR as

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{nsf}}\,=\,1.1\times {10}^{40}\,\mathrm{erg}\ {{\rm{s}}}^{-1}\ \left(\displaystyle \frac{\epsilon }{0.3}\right)\left(\displaystyle \frac{\mathrm{SFR}}{0.1\,{M}_{\odot }}\right),\end{eqnarray} \tag{ 12 }$$

where ${\dot{E}}_{\mathrm{nsf}}$ is the mechanical energy input rate due to nuclear star formation and is an efficiency factor typically assumed to be $\approx 0.3$ (see Crocker et al. 2015 or Sarkar et al. 2015 for equivalent relations). This implies that an average SFR of 0.05 M_⊙ yr⁻¹ over the past ∼10 Myr produces an energy injection rate of $\approx 6\times {10}^{39}$ erg s⁻¹. This estimate falls $\approx 400$ times lower than our estimated input rate of $2.3\times {10}^{42}$ erg s⁻¹. It is possible that the SFR has been more variable over the past ∼10 Myr, but the upper limits are only $\approx 3$ times higher than the average value (Yusef-Zadeh et al. 2009). Thus, star formation in the Galactic center does not produce enough energy to inflate the bubbles based on our energy injection rate constraints.

We use our density constraints to estimate the thermal gas mass within the bubble and shell structures. This is a straightforward calculation since we assume that each component has a constant density and fixed volume. Thus, the mass in each component is defined as $M=\mu {m}_{{\rm{H}}}\times n\times V$, where $\mu \,=\,0.61$ is the average weight per particle, m_H is the mass of hydrogen, n is the inferred density, and V is the volume. Our geometric models imply a bubble volume of ${V}_{\mathrm{FB}}\,=\,2\times (4/3)\times \pi \times 5\times {3}^{2}=377$ kpc³ (the factor of two is for two ellipsoidal bubbles), and a combined shell volume of ${V}_{\mathrm{shell}}\,=\,411$ kpc³. The densities in Table 3 imply thermal gas masses of ${M}_{\mathrm{FB}}\,=\,4.6\times {10}^{6}\,{M}_{\odot }$ and ${M}_{\mathrm{shell}}\,=\,6.1\times {10}^{6}\,{M}_{\odot }$ for the bubble and the shell with a characteristic range between 5 and 10 × 10⁶ M_⊙ given the density uncertainties. These masses represent material that has been shock-heated by the bubbles or injected into the bubbles by the energy source.

We first explore whether the bubble and shell plasmas are predominantly shocked/mixed hot halo material by comparing the masses derived above to the inferred hot gas halo mass that would exist within the bubble and shell volumes. If the density inside the bubble+shell volume was defined by our hot gas halo density model with r_c = 3 kpc (instead of the Fermi bubble/shell densities), the halo mass in the volume would be $1.11\times {10}^{7}\,{M}_{\odot }$. The calculations above indicate that the combined thermal gas mass within the bubble+shell volumes is ${M}_{\mathrm{FB}}+{M}_{\mathrm{shell}}\,=\,1.07\times {10}^{7}\,{M}_{\odot }$. The consistency between these values suggests that most of the thermal gas associated with the bubbles is shock-heated ambient material, as opposed to material injected by the energy source.

As a consistency check, we estimate the amount of material injected by the energy source, either AGN or star formation, to compare with the above masses. The mass-loss rate due to nuclear star formation activity (or mass injection rate) is believed to be ${\dot{M}}_{\mathrm{inj},\mathrm{nsf}}\approx 0.3(\mathrm{SFR}/\,{M}_{\odot }\,{\mathrm{yr}}^{-1})$ (Leitherer et al. 1999). Mass-loss rates from black hole accretion events (from either jets or winds) are more uncertain, but simulations of RIAF accretion winds suggest values ranging between 2% and 20% of ${\dot{M}}_{\mathrm{Edd}}\,=\,10{L}_{\mathrm{Edd}}/{c}^{2}$ (Yuan et al. 2012, 2015). Assuming a nuclear SFR of 0.05 M_⊙ yr⁻¹, ${\dot{M}}_{\mathrm{Edd}}\sim {10}^{-1}\,{M}_{\odot }$ yr⁻¹ for Sgr A*, and an active period for the bubble of 4.3 Myr, we estimate that the injected mass is ${\dot{M}}_{{\rm{inj}}}\lesssim {10}^{5}\,{M}_{\odot }$ for both origin scenarios. This is significantly less than our constrained mass estimate of $\sim {10}^{7}\,{M}_{\odot }$, thus validating our claim that the bubbles contain predominantly shocked halo gas.

5.4.1. Analyses at Soft X-Ray Energies

The most direct comparisons to our analysis are the previous soft X-ray spectral analyses (Kataoka et al. 2013, 2015; Tahara et al. 2015). These studies follow a similar methodology and find similar results, with K15 being the most current and comprehensive work of the three. These authors compiled a sample of 29 Suzaku observations and 68 Swift observations distributed across the Fermi bubbles. They fit the Suzaku XIS data and Swift X-ray Telescope spectra with a multi-component thermal plasma model, where one component is typically fixed at ${kT}=0.1\,{\rm{keV}}$ to represent emission from the LB and residual SWCX, and the other represents the combined emission from the halo and Fermi bubbles. They systematically find a hot gas halo/Fermi bubble plasma temperature of ${kT}=0.3\,{\rm{keV}}$, and this is hotter than the characteristic value found for sight lines away from the bubbles (${kT}=0.2\,{\rm{keV}};$ Henley & Shelton 2013). From this temperature, they infer a relatively low Mach number of ${ \mathcal M }$ = 0.3 keV/0.2 keV = 1.5 and an expansion velocity of ${v}_{\exp }\,=\,300$ km s⁻¹. They also find emission measures that vary by over an order of magnitude, which they note is due to a combination of emission from the Galactic halo and from the Fermi bubbles. The authors model the emission measures in the northern Galactic hemisphere with a hot gas halo density model from Miller & Bregman (2013) and a Fermi bubble shell distribution with an inner radius of 3 kpc, an outer radius of 5 kpc, and a density of $3.4\times {10}^{-3}$ cm⁻³.

While our approach is similar to these studies, there are several differences that can explain why we derive a higher bubble temperature and lower bubble density. The main observational difference is that these studies fit the SXRB spectra for emission measures and temperatures, while we chose to measure oxygen emission line intensities. Emission measures and temperatures are more useful plasma properties to measure, but fitting a full 0.5–2.0 keV SXRB spectrum with a thermal plasma model requires more counts than fitting only the oxygen emission lines. This is why our sample is larger and has better sky coverage than the K15 sample. However, these observables should give consistent results for the inferred bubble temperature and density. The measured plasma temperature is most sensitive to the O viii/O vii ratio for temperatures between ≈0.1 and 0.3 keV, because the lines are strong and rapidly changing in strength in this regime (Figure 11). Thus, any temperature derived from fitting an SXRB spectrum with an APEC model should be consistent with the temperature inferred from fitting the oxygen lines separately (Yoshino et al. 2009). The analysis becomes more complicated when there are multiple emission components, and the interpretation of multiple SXRB sources is likely the bigger difference between approaches.

Zoom In Zoom Out Reset image size

The primary difference between our work and those discussed above is the treatment of combined X-ray emission from the hot gas halo and Fermi bubbles. The observed emission includes contributions from the hot halo and Fermi bubbles. Thus, the fitted plasma temperature of 0.3 keV is the emission measure-weighted temperature of the hot halo at 0.2 keV and a Fermi bubble plasma that is likely >0.3 keV. The K15 analysis assumes that the Fermi bubbles dominate the observed emission, while our analysis includes the combined emission from the hot halo and Fermi bubbles. This extension leads to a similar result to these previous works, but also explains why we infer a higher temperature for the Fermi bubble plasma than these studies (kT_FB,shell ≈ 0.4–0.5 keV).

A similar interpretation likely explains why the K15 analysis infers bubble densities 3–4 times higher than our constraints. Their Fermi bubble geometric distribution includes only a shell component that is 2 kpc thick, while ours includes both a volume-filled component and a shell component. This implies that our bubble+shell emission model has a longer path length along most sight lines near the bubbles than their model. The emission measure scales with density and path length as ${EM}\propto {n}^{2}L$, so a longer inferred path length would lead to a lower inferred density. We also point out that their hot gas halo density model extends only to r = 20 kpc. While the hot halo emission is likely dominated by gas within $r\lesssim 25\,{\rm{kpc}}$, failing to account for emission at greater radii can decrease the amount of modeled halo emission. The combined effect is that the K15 analysis assumes shorter bubble path lengths and less halo emission than our emission model, and this results in a higher inferred bubble density required to match the total observed emission.

5.4.2. Kinematic Estimates from UV Absorption Lines

5.4.3. Comparing with Simulations