THE ORIGIN OF THE HOT GAS IN THE GALACTIC HALO: CONFRONTING MODELS WITH XMM-NEWTON OBSERVATIONS

The model spectra discussed here were generated from one-dimensional Lagrangian hydrodynamical simulations of SNRs evolving at a variety of heights above the Galactic disk, and hence in a variety of ambient densities. The simulations are described fully in Shelton (1998, 1999, 2006), and include radiative cooling, thermal conduction, and an effective ambient magnetic field, B_eff, which exerts a non-thermal pressure, in addition to the ambient gas pressure. The ionization evolution in the shocked gas is modeled self-consistently.

We will first concentrate on the models from Shelton (2006) with SN explosion energy E₀ = 0.5 × 10⁵¹ erg and B_eff = 2.5 μG (corresponding to a non-thermal pressure P_nt = 1800 cm⁻³ K), and then discuss varying these parameters. We consider the models evolving in ambient densities n₀ = 0.2, 0.1, 0.05, 0.02, 0.01, and 0.005 cm⁻³, corresponding to heights z = 190, 310, 480, 850, 1300, and 1800 pc, using the interstellar density model from Ferrière (1998). Note that we are ignoring the model from Shelton (2006) at z = 76 pc, as we do not consider this to be in the halo. Density and temperature profiles from the model with n₀ = 0.01 cm⁻³ are shown in Figure 5. The SN explosion blows a hot rarefied bubble in the ambient medium; this hot bubble is the source of the X-rays. The model shown in Figure 5 corresponds to model A in Shelton (1999); see that paper for more details.

Zoom In Zoom Out Reset image size

For each epoch of each SNR model, we calculated spectra for a range of impact parameters through the model remnant, using the Raymond & Smith (1977) spectral code (updated by J. C. Raymond & B. W. Smith 1993, private communication with R. J. Edgar). The spacing between the impact parameters, Δb, was adjusted according to the size of the remnant. The spectral calculation takes into account the possible non-equilibrium ionization calculated during the hydrodynamical simulation.

We wish to compare the distributions of observed temperatures and EMs with those predicted by the model. In order to do this, we need to calculate the probability of a sightline passing through a remnant at a given height, of a given age, and at a given impact parameter. For example, a sightline is more likely to pass through a remnant at a larger impact parameter than a small impact parameter, and more likely to pass through a remnant closer to the disk, where the SN rate is larger. If R(z₁, z₂) is the rate per unit area of SNe at heights between z₁ and z₂, then the probability P(z₁, z₂, t, Δt, b, Δb) of intercepting a remnant between these heights with an age between t and t + Δt at an impact parameter between b and b + Δb is

$$\begin{equation} P (z_1,z_2,t,\Delta t,b,\Delta b) = 2 \pi b R(z_1,z_2) \Delta b \Delta t. \end{equation} \tag{ 5 }$$

In order to calculate the above probabilities, we followed Shelton (2006) and divided the halo into six plane-parallel slabs. We assumed that all the SNRs within a given slab experience a single ambient density, and so can be represented by one of our six SNR models. We calculated R(z₁, z₂) for each slab by integrating the volumetric SN rate at the solar circle, r(z), from Ferrière (1998):

$$\begin{eqnarray} r(z) &=& r_\mathrm{Ia}(z) + r_\mathrm{II}(z) \nonumber \\ &=&(4.0 e^{-|z|/325 \mathrm{pc}} + 14 e^{-|z|/266\,\mathrm{pc}}) \;\mbox{kpc}^{-3} \;\mbox{Myr}^{-1}, \end{eqnarray} \tag{ 6 }$$

where r_Ia(z) and r_II(z) are the rates of Type Ia SNe and isolated core-collapse SNe, respectively. The slab boundaries were placed at the midpoints between the nominal heights of the model remnants, except for the lower boundary of the lowest slab, which was placed at 130 pc (the scale height of the Galactic H i layer), and the upper boundary of the highest slab, which was placed at infinity. The boundaries, ambient densities, and SN rates for the six slabs are shown in Table 3. We calculated Δt from the timestamps of the output files from the hydrodynamical simulations. As noted above, Δb is the spacing between the impact parameters for which we calculated spectra from a given remnant.

We used a Monte Carlo method to construct model spectra for 2000 sightlines, using the above-calculated probabilities. In our Monte Carlo simulation, 64% of the model sightlines intercepted no remnants, and 9% of the sightlines intercepted more than one remnant. For model sightlines intercepting more than one remnant, we summed the spectra of the individual remnants.

In our analysis of the XMM-Newton observations, we modeled the halo X-ray emission with a single-temperature CIE plasma model, whereas the true halo emission is likely from plasma at a range of temperatures and in a range of ionization states. Similarly, the emission predicted by this extraplanar SNR model is from plasma at a range of temperatures and in a range of ionization states. Therefore, in order to compare the predictions of this model with our observational results, we first characterize the model SNR spectra calculated above with 1T CIE plasma models, by simulating XMM-Newton observations of the SXRB.

Our procedure for characterizing the spectrum for each model sightline is as follows. We first multiplied the spectrum by a renormalization factor, k_rn, in order to give a 0.4–2.0 keV surface brightness of 2.06 × 10⁻¹² erg cm⁻² s⁻¹ deg⁻². This value is the median intrinsic halo surface brightness inferred from the best-fit models in Table 2. We then subjected the model to an absorbing column N_H = 1.7 × 10²⁰ cm⁻²; again, this is the median value used in our XMM-Newton analysis. To the absorbed SNR spectrum, we added components representing the foreground emission, the extragalactic background, the soft protons, and the instrumental lines. The parameters for the foreground emission and the extragalactic background were taken from our observational analysis (see Section 3.1). The normalization of the foreground component was chosen to give an R12 count rate of 600 counts s⁻¹ arcmin⁻² (the median value from Table 2). The parameters for the soft-proton model and the instrumental lines were taken from the best-fit model for obs. 0301330401; this observation was chosen because it has close to the median level of soft-proton contamination, as judged by the ratio of the observed 2–5 keV flux to that expected from the extragalactic background (F^2–5_total/F^2–5_exgal; see Paper I).

We simulated an observation of the model spectrum by folding the above-described multicomponent model through the XMM-Newton response, assuming a typical field of view of 480 arcmin², and adding Poissonian noise to the spectrum, assuming a typical observing time of 15 ks. Our simulations also took into account the QPB spectrum. For each model sightline, we simulated a MOS1 and a MOS2 spectrum. We binned the resulting spectra such that there were at least 25 counts per bin. Note that our assumed field of view is not as large as the full XMM-Newton field of view, as we removed bright sources from the fields of our observations, and for some observations not all chips were usable.

We fitted the resulting simulated spectra with the same multicomponent model that we used in our observational analysis. In particular, the halo emission was modeled with an absorbed Raymond & Smith model, with N_H fixed at 1.7 × 10²⁰ cm⁻². For each model sightline, we fitted to the simulated MOS1 and MOS2 spectra simultaneously. Figure 6 shows a simulated XMM-Newton spectrum from one of our model sightlines, along with the best-fitting model. Also shown is the 1T CIE halo component of the best-fitting model, as well as the input SNR model. The input SNR spectrum is very different from that of a CIE plasma. The SNR plasma is relatively cool (as shown by the differential EM, plotted as the solid line in Figure 7) and overionized, and the spectrum exhibits strong recombination edges. If the SNR plasma were in CIE, instead of recombining, its emission would be four orders of magnitude fainter in the XMM-Newton band. However, when we use a 1T CIE plasma to model the halo component in our simulated spectra, the fits are generally good. This is because much of the spectral detail in the input SNR spectrum is lost when it is combined with the other components of the SXRB and folded through the XMM-Newton response. There is some excess emission at ∼0.75 keV in the simulated spectrum in Figure 6, possibly due to O⁺⁷ → O⁺⁶ recombinations (the recombination edge is at 0.74 keV), but it would not be easy to unambiguously identify such a feature in an observed spectrum as recombination emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used the resulting fit parameters to compare with the observations. The temperatures were taken directly from the fits, while the best-fit EMs were first divided by the relevant values of k_rn before comparing with the observations. In what follows, we refer to these values as "X-ray temperatures" and "X-ray EMs," to emphasize that they are derived from simulated X-ray observations, rather than being derived directly from the hydrodynamical data.

Our model X-ray EMs were derived for sightlines looking straight upward from the disk, i.e., toward |b| = 90°. Our observations, however, are toward a range of Galactic latitudes, and so will sample different amounts of halo material. We therefore multiply our observed EMs by sin |b| before comparing them with the model predictions. This transformation assumes that the halo is, in a statistical sense at least, uniform in directions parallel to the disk.

Figures 8(a) and (b) show histograms comparing the halo X-ray temperatures and EMs predicted by the extraplanar SNR model with the corresponding halo properties obtained from our XMM-Newton observations. (Note that obs. 0406630201, which yielded T = 11.2 × 10⁶ K, is not shown in Figure 8(a).)

Zoom In Zoom Out Reset image size

The X-ray temperatures predicted by the model are in reasonably good agreement with the observed temperatures. However, the emission predicted by this model is clearly too faint: the predicted X-ray EMs are typically an order of magnitude smaller than the observed values. Among the model sightlines that intercept extraplanar SNRs, the median X-ray EM is 1.37^+0.09_−0.12 dex smaller than the median observed value of EMsin |b|: 0.081 × 10⁻³ versus (1.90^+0.41_−0.46) × 10⁻³ cm⁻⁶ pc.⁹ As the predicted X-ray temperatures are in reasonable agreement with the observed temperatures, the EM result implies that the X-ray surface brightnesses are also typically underpredicted by an order of magnitude. Note also that our Monte Carlo simulation predicts that about two-thirds of the sightlines would intercept no extraplanar SNRs, whereas we observe hot halo gas on most, if not all, of our sightlines.

The model values in Figures 8(a) and (b) were calculated from the SNR simulations in Shelton (2006) with E₀ = 0.5 × 10⁵¹ erg and B_eff = 2.5 μG. We find that increasing E₀, B_eff, or the assumed SN rate all increase the predicted X-ray EMs. Table 4 shows the median X-ray EMs predicted by extraplanar SNR models with different explosion energies, ambient magnetic fields, and SN rates. As can be seen, increasing B_eff has the largest effect on the predicted X-ray EMs. However, it should be noted that the median model values in Table 4 are only for the subset of sightlines that intercept at least one SNR. Increasing E₀ or the SN rate increases the fraction of sightlines that intercept at least one SNR, while increasing B_eff decreases that fraction.

Figure 8(c) shows the histogram of X-ray EMs predicted by the Shelton (2006) simulations with E₀ = 1 × 10⁵¹ erg and B_eff = 5.0 μG, with a Galactic SN rate that is twice that given by Equation (6) (i.e., the fifth model in Table 4). The histogram of observed values is also plotted for comparison. The model still significantly underpredicts the observed EMs (it underpredicts the median value by 0.91^+0.09_−0.12 dex). Furthermore, our new Monte Carlo simulation predicts that ∼40% of the sightlines would intercept no extraplanar SNRs—we reiterate that we observe hot halo gas on most, if not all, of our sightlines. We will discuss the results presented here in Section 5.2.2.

The hydrodynamical simulation used here is described fully in Joung & Mac Low (2006), and the reader is referred to that paper for more details. The simulation was carried out using Flash,¹⁰ a parallelized Eulerian hydrodynamical code with adaptive mesh refinement. The simulation box extends from z = −5 kpc to z = +5 kpc, and has a size of (0.5 kpc)² in the xy plane. The maximum spatial resolution is 1.95 pc. The upper and lower boundaries have outflow boundary conditions, while the vertical sides of the simulation box have periodic boundary conditions.

The simulation box was initialized with 1.1 × 10⁴ K gas in hydrostatic equilibrium. This gas was then heated by discrete SN explosions, each of which injected 10⁵¹ erg of energy into a small region of the grid. These explosions generally occurred randomly in time and space, with a rate

$$\begin{eqnarray} r(z) &=& r_\mathrm{Ia}(z) + r_\mathrm{II}(z) \nonumber \\ &=&(6.2 e^{-|z|/325 \mathrm{pc}} + 167 e^{-|z|/90 \mathrm{pc}}) \;\mbox{kpc}^{-3} \;\mbox{Myr}^{-1}, \end{eqnarray} \tag{ 7 }$$

although 3/5 of the Type II SNe occurred in clusters of seven to ≈40 explosions. The above rates differ from those in Equation (6) in two ways. First, Joung & Mac Low (2006) assumed higher Galactic SN rates than Ferrière (1998): 1/330 yr⁻¹ versus 1/445 yr⁻¹ (Type I) and 1/44 yr⁻¹ versus 1/52 yr⁻¹ (Type II). Second, Equation (6) considers only isolated Type II SNe, whose average height is larger than the average height of all Type II progenitors (266 pc versus 90 pc; Ferrière 1995). The model also includes diffuse heating, due to photoelectric emission from UV-irradiated dust grains, and radiative cooling. Figure 9 shows vertical slices of the density and temperature at t = 120.0 Myr.

Zoom In Zoom Out Reset image size

We calculated X-ray spectra for 242 sightlines looking vertically upward and downward from the Galactic midplane. The viewpoints of these sightlines formed an 11 × 11 grid in the xy plane, with a grid spacing of ≈49 pc. The spectra were calculated using the Raymond & Smith (1977) spectral code (updated by J. C. Raymond & B. W. Smith 1993, private communication with R. J. Edgar), assuming that the material along the line of sight is in CIE. We assumed that the gas on the grid is optically thin. We ignored the emission from the first 100 pc of each sightline, as such material is not in the Galactic halo. In our observational analysis, emission from within ∼100 pc of the midplane is attributed to our foreground model component, derived from ROSAT shadowing data.

As with the previous model, we used the method described in Section 4.2.2 to characterize the model spectra with 1T models, and used the resulting X-ray temperatures and X-ray EMs to compare with our observations. As in Section 4.2.2, our simulated XMM-Newton spectra included foreground, extragalactic background, and instrumental background components, as well as the halo emission from the SN-driven ISM model. The emission from the halo and extragalactic components was subjected to absorption with N_H = 1.7 × 10²⁰ cm⁻². This column density was also used in the subsequent spectral fitting, from which we obtained the characteristic X-ray temperatures and X-ray EMs. Figure 10 shows a simulated XMM-Newton spectrum for one of our model sightlines, along with the best-fitting model. The differential EM for the input halo spectrum is shown in Figure 7 (dashed line). Because we assume that the plasma is in CIE, only plasma with T ≳ 10⁶ K will contribute in the XMM-Newton band. Although the input model predicts emission from a range of temperatures, a model with a 1T halo generally fits the simulated spectra well.

Zoom In Zoom Out Reset image size

It should be noted that our model spectra do not take into account absorption by cold gas on the grid—when we characterized our model spectra, we assumed that the absorbing column (N_H = 1.7 × 10²⁰ cm⁻²) was located entirely beneath the hot X-ray-emitting gas. For ∼2/3 of the model sightlines, the column density of cold gas that is mixed in with the bulk of the X-ray-emitting gas is <10¹⁹ cm⁻² (i.e., an order of magnitude smaller than the column densities used in our XMM-Newton analysis). We have investigated the effect of ignoring on-grid absorption by creating simulated XMM-Newton spectra from an XSPEC model of the form $\mathtt {phabs} \ast (\mathtt {raymond} + \mathtt {phabs} \ast \mathtt {raymond})$, and characterizing the resulting spectra with a model of the form $\mathtt {phabs} \ast \mathtt {raymond}$. As with our previous simulations, the column density of the first phabs component was fixed at 1.7 × 10²⁰ cm⁻². The two raymond components in the input model had the same EM, and temperatures of 1.5 × 10⁶ K and 2.5 × 10⁶ K (it does not matter which component is the hotter—our conclusion is the same either way). We found that the characteristic X-ray temperatures and X-ray EMs obtained were not strongly affected by the column density between the two raymond components in the input model, at least up to column densities of a few × 10²⁰ cm⁻². We therefore conclude that ignoring on-grid absorption will not adversely affect our results.

In the following, we also compare the predicted X-ray surface brightnesses with our observations. These values were extracted directly from the model spectra.

Figure 11 shows the time variation of the X-ray spectral properties predicted by the SN-driven ISM model. The model has not settled down to a steady state—the typical X-ray temperature and X-ray surface brightness rise steadily throughout the period shown. However, in this temperature regime, an increase in X-ray surface brightness can be brought about by an increase in temperature as well as by an increase in EM. The typical X-ray EM does not rise steadily throughout the plotted period.

Zoom In Zoom Out Reset image size

Figure 11 also compares the predicted X-ray properties with the observed values from our XMM-Newton analysis. Similarly to Section 4.2.3, we multiplied our observed EMs and surface brightnesses by sin |b| before comparing them with the model values. At the earlier times, the model halo is too cool and too faint in the 0.4–2.0 keV band. Around t = 110 Myr, the predicted X-ray temperatures are in reasonable agreement with the observed value, but the halo is typically a factor of ∼4 too faint. At later times the model halo is too hot; by t = 155 Myr it is also somewhat brighter than is observed, although the median predicted surface brightness is within 50% of the observed median.

From Figure 11 it is not clear what state the halo is tending toward. de Avillez & Breitschwerdt (2004) estimated that a steady state halo in a simulation such as this should be reached in ≲180 Myr. However, it is possible that, instead of settling down to a steady state that is hotter and brighter than the observed halo, the variation in X-ray temperature and surface brightness is part of a slow oscillation about a mean state, with a period of at least ∼130 Myr. Determining the final state predicted by this model would probably require running it for at least several more tens of megayears, which is not currently practical. We therefore try two different approaches in comparing the simulation data to our observations.

Our first approach is to assume that the data shown in Figure 11 represent roughly half a cycle in the oscillation of the halo about a mean state. The simulation models only a narrow column of the halo. As our observation directions are not all toward the Galactic poles, our sightlines sample different spatial locations in the halo, which, in a statistical sense, should correspond to different times throughout the cycle. We therefore proceed by averaging the spectra for each sightline from seven of the eight time steps shown in Figure 11 (we do not use the data from t = 118 Myr, so as not to oversample the times around t ∼ 120 Myr). We characterize the resulting averaged spectra with 1T models, using the procedure described previously. The results are compared with our observations in Figure 11 (the rightmost data point, labeled "Ave"), and in Figures 12(a) and (b).

Zoom In Zoom Out Reset image size

Our second approach is to assume that the final time in Figure 11 (t = 155 Myr) gives our best estimate of the steady state that the halo is approaching. We use the predictions from that final time in Figures 12(c) and (d).

Whether we use the averaged spectra or the final time step, the predicted halo temperature is too high. The median predicted temperatures are 2.73 × 10⁶ K (averaged spectra) and 3.55 × 10⁶ K (final time step), against an observed median of (2.08^+0.11_−0.07) × 10⁶ K. The averaged spectra underpredict the halo EM, although the discrepancy is not as large as with the extraplanar SNR model. The median predicted X-ray EM is 0.44^+0.09_−0.12 dex smaller than the median observed value of EMsin |b|: 0.69 × 10⁻³ cm⁻⁶ pc versus (1.90^+0.42_−0.46) × 10⁻³ cm⁻⁶ pc. For comparison, the discrepancy for the extraplanar SNR model is 1.37^+0.09_−0.12 dex. For the spectra from the final time step, the median predicted X-ray EM (1.17 × 10⁻³ cm⁻⁶ pc) is in better agreement with the observed median (the discrepancy is 0.21^+0.09_−0.12 dex). The predicted 0.4–2.0 keV surface brightnesses are generally in reasonably good agreement with the observed values, although there are a few sightlines in the final time step of the model that exhibit much greater surface brightnesses. The median predicted values are within 50% of the median observed value of S_0.4–2.0sin |b|: 1.04 × 10⁻¹² (averaged spectra) and 2.37 × 10⁻¹² (final time step), versus (1.62^+0.30_−0.05) × 10⁻¹² erg cm⁻² s⁻¹ deg⁻² (observed). We will discuss the results presented here in Section 5.2.3.

In Section 4.1, we examined a disk galaxy formation model (Toft et al. 2002; Rasmussen et al. 2009), which predicts the existence of a hot halo extended over tens of kiloparsecs, formed from material falling into the galaxy's potential well. Using results from Rasmussen et al. (2009), we find that the emission-weighted mean temperature predicted for a Milky-Way-sized galaxy (v_c ∼ 220 km s⁻¹) is in good agreement with the median observed halo temperature, but the model underpredicts the X-ray luminosity of the halo by at least an order of magnitude.

The above results suggest that the extended hot halo predicted by disk galaxy formation models is not a major contributor to the halo X-ray emission that we observe with XMM-Newton. However, this interpretation is complicated by the simulations of Crain et al. (2010). They argue that the stellar feedback in the Rasmussen et al. (2009) simulations is too weak, resulting in too much mass ending up in stars and too little mass in hot gas. Ideally, we would compare the surface brightness of the accreted extragalactic material in the Crain et al. (2010) simulations with the halo surface brightness obtained from our XMM-Newton observations, but this is not possible for two reasons. First, in addition to infalling extragalactic material, the Crain et al. model includes interstellar gas that was heated by stars and SNe and transferred upward by an outflow (similar to the galactic fountain in the model described in Section 4.3). Second, their model predictions are in terms of luminosity rather than surface brightness. Because gas heated by stars and SNe will tend to be concentrated nearer the disk, compared with the more extended halo of accreted extragalactic material, it may provide a much larger surface brightness per unit luminosity than the extended halo gas. Predictions of the separate contributions to the surface brightness due to the gas heated by stellar processes and the extended halo are not currently available from the Crain et al. model. Such predictions are needed to determine whether or not accreted extragalactic material makes a major contribution to the observed halo X-ray emission.

Although surface brightness predictions are not currently available, we can make rough comparisons between the total X-ray luminosity predicted by Crain et al. (2010) and the halo luminosity expected from the XMM-Newton data. An L_X–T plot derived from their simulations shows that galaxies with emission-weighted halo temperatures of ∼2 × 10⁶ K have 0.3–2.0 keV luminosities of ∼(3–40) × 10³⁹ erg s⁻¹ (R. Crain 2010, private communication), compared with <0.2 × 10³⁹ erg s⁻¹ from Rasmussen et al. (2009). In Section 4.1, we derived a Galactic halo luminosity of 1.9 × 10³⁹ erg s⁻¹ from our observations, which is slightly smaller than the range of luminosities predicted by Crain et al. (2010). However, this observed luminosity was calculated assuming that the observed emission came from a uniform sphere of radius R_sph = 15 kpc. Figure 2 in Crain et al. (2010) implies that the emission in their model is much more extended than this, with R_sph ∼ 50 kpc. Using this larger radius leads to an observed halo luminosity of 1.3 × 10⁴⁰ erg s⁻¹, which lies within the range predicted by Crain et al.

In Section 4.2, we examined a model in which the hot halo gas is contained in an ensemble of isolated extraplanar SNRs (Shelton 2006). The X-ray temperatures predicted by this model are in good agreement with the observed values. However, the predicted X-ray EMs, obtained from model spectra above 0.4 keV, are typically an order of magnitude smaller than the observed values. We can increase the predicted X-ray EMs by increasing the SN explosion energy, E₀, the effective ambient magnetic field, B_eff, or the SN rate. However, doubling all three of these parameters still resulted in X-ray EMs that were too small.

Of the above three parameters, B_eff has the largest effect on the predicted X-ray EMs (see Table 4). Increasing B_eff increases the non-thermal pressure, and therefore increases the compression of the hot X-ray-producing bubble. This increased compression increases the density and temperature of the bubble, increasing the X-ray EM inferred from the predicted emission in the XMM-Newton band. We have not carried out simulations with B_eff>5.0 μG, but if we extrapolate the values in Table 4, we find we would need a B_eff of a few tens of microgauss to match the observed EMs. Such a magnetic field implies a non-thermal pressure P_nt/k ≳ 10⁵ cm⁻³ K. This is an implausibly high non-thermal pressure for the halo, as it is several times larger than the midplane value (Boulares & Cox 1990; Ferrière 2001). In addition, the model with B_eff = 2.5 μG predicts that hot gas would be seen on only ∼1/3 of sightlines. As noted above, increasing B_eff results in smaller, hotter remnants. In addition, these smaller remnants are brighter and so shorter lived. As a result, increasing B_eff means that even fewer sightlines would intercept remnants. In reality, we find hot gas on most, if not all, of our XMM-Newton sightlines. Therefore, we cannot bring the extraplanar SNR model into agreement with our observations by increasing the assumed effective ambient magnetic field.

Increasing the SN rate increases the number of SNRs expected to lie along a given sightline, and so increases the predicted X-ray EMs. However, we would have to increase the SN rate given by Equation (6) by a factor of ∼6 in order to give the same increase in the median X-ray EM that we see when we increase B_eff from 2.5 to 5.0 μG. To match the observations, we would have to increase the SN rate by an even larger factor. The Galactic SN rate of ∼2 per century is constrained to within a factor of ∼2 (see the online Supplementary Information¹¹ for Diehl et al. 2006). If the SN scale heights in Equation (6) are well constrained, the halo SN rate is also constrained to within a factor of ∼2. Even if we double the SN scale heights in Equation (6), the integrated SN rate above |z| = 130 pc only increases by a factor of 2.5. Therefore, to bring the extraplanar SNR model into agreement with our observations we would have to increase the SN rate by an unrealistic amount. We also cannot bring the model into agreement with our observations by increasing E₀, as Table 4 shows that the median X-ray EM depends only weakly on E₀.

We therefore conclude that the hot halo gas that we observe with XMM-Newton is not primarily due to an ensemble of isolated extraplanar SNRs. Most of this population of remnants would be relatively old (age>1 Myr) and faint in the XMM-Newton band. It should be noted, however, that this result does not imply that extraplanar SNRs do not contribute to the hot halo gas at all. In the XMM-Newton band (above 0.4 keV), their contribution is masked by brighter emission from an SN-driven galactic fountain (see Section 5.2.3). At lower energies, extraplanar SNRs could still contribute significantly to the hot gas observed in the 1/4 keV band, as originally suggested by Shelton (2006). In addition, young, bright remnants, although rare, do produce emission that is detectable by Suzaku (Henley & Shelton 2009) and that should also be detectable by XMM-Newton.

In Section 4.3, we examined a hydrodynamical simulation of vertically stratified interstellar gas, driven by SN explosions (Joung & Mac Low 2006). Unlike the extraplanar SNR model, the SNRs do not evolve in isolation, and the model results in a galactic fountain of hot gas up into the halo (Shapiro & Field 1976; Bregman 1980). Also, for this model we assumed that the gas is in CIE, whereas the extraplanar SNR model includes self-consistent modeling of the ionization evolution.

As with the extraplanar SNR model, we folded the model spectra through the XMM-Newton response, added photon noise, and characterized the resulting spectra with 1T models. Before going on to discuss the comparison of the resulting X-ray temperatures and EMs with our XMM-Newton observations, we shall look at how well these X-ray spectral properties reflect the properties of the gas on the hydrodynamical grid.

Figure 14(a) compares the X-ray temperatures derived from the 1T models with the mean temperature along each model sightline through the hydrodynamical grid, weighted by the 0.4–2.0 keV emission. We find that the X-ray temperatures are in reasonable agreement with the emission-weighted mean temperatures (typically within 0.4 × 10⁶ K). Figure 14(b) compares the X-ray EM for each sightline with the EM of the gas with T>10⁶ K along that sightline (we use this quantity because there is no obvious analog to the emission-weighted mean temperature). The X-ray EM underestimates the EM of gas with T>10⁶ K, although the two quantities are well correlated and generally agree within a factor of 2.5. Overall, we conclude that the properties derived from the 1T models are good characterizations of the hot gas on the hydrodynamical grid

Zoom In Zoom Out Reset image size

The 0.4–2.0 keV halo surface brightness predicted by this model is in good agreement with the observations (within ∼50%). As the emission predicted by this model is much brighter, it will mask the contribution from extraplanar SNRs. Note that this model does include some in situ heating by extraplanar SNe. The integrated halo SN rate above |z| ∼ 100 pc is ∼2 times larger in this model than in the extraplanar SNR model (although the volumetric SN rate falls off more rapidly with |z| in this model). However, this larger integrated SN rate is not enough to explain the differences in X-ray surface brightness between the models: this model predicts emission an order of magnitude brighter than the extraplanar SNR model. The majority of the halo X-ray emission in this model comes from hot gas that is driven from the disk into the halo by a galactic fountain. As mentioned in Section 5.2.1, it is possible that the emission from a galactic outflow or fountain also masks the contribution from an extended hot halo of accreted material.

Although the median predicted X-ray surface brightness and X-ray EM are in good agreement with the observed medians, the halo X-ray temperature predicted by this model is ∼1 × 10⁶ K larger than our measured value, independent of whether the simulated halo is undergoing a long-period (P ≳ 100 Myr) oscillation about some mean state, or is settling down to a steady state. Given the uncertainties in the model predictions (due to the time variability shown in Figure 11), these discrepancies are probably insufficient to rule out this model. Nevertheless, we discuss below three possible causes for the temperature discrepancy: (1) the simulation overpredicts the temperature of the hot gas in the halo, (2) the temperature of the hot gas does not accurately predict the X-ray spectrum, because the gas is not in ionization equilibrium, or (3) a bias in our spectral analysis causes us to underestimate the observed halo X-ray temperature. Correcting for whichever of these effects turn out to be important may also affect the X-ray surface brightnesses and EMs. However, it is not easy to foresee whether these corrections will increase or decrease the discrepancies between the predicted and observed surface brightnesses and EMs.

The first possible cause of the temperature discrepancy is that the simulation overpredicts the temperature of the hot gas in the halo. The simulation does not include thermal conduction, which could in principle lower the temperature of the hot X-ray-emitting gas, although thermal conduction in the ISM can be suppressed by magnetic fields. de Avillez & Breitschwerdt (2004) argue that turbulent diffusion is more efficient at mixing hot and cold gas than thermal conduction. However, the efficiency of this mixing will be underestimated if the mixing is not fully resolved. In the simulation used here, the spatial resolution is between 1.95 pc and 15.6 pc, but the adaptive mesh refinement criterion is chosen to focus maximum numerical resolution on the 200 pc above and below the midplane, so the gas is less well resolved in the halo than in the disk. As a result, it is possible that turbulent mixing of hot and cold gas is underresolved in the halo, resulting in an overestimate of the temperature of the X-ray-emitting gas (Fujita et al. 2009).

A lack of spatial resolution may also suppress the radiative cooling of the hot gas—averaging the hot gas density over large cells eliminates local denser regions that would radiate more efficiently (de Avillez & Breitschwerdt 2004). de Avillez & Breitschwerdt (2004) found that the filling factor of hot (T>10^5.5 K) gas was not significantly affected by the spatial resolution of their simulations. However, there is insufficient information to determine how the emission-weighted mean temperature of the X-ray-emitting gas would be affected. Simulations with higher spatial resolution in the halo will help determine whether or not the mixing of hot and cold gas and the radiative cooling of hot gas are adequately resolved for predicting the X-ray emission.

Joung & Mac Low (2006) pointed out that the average gas density at several disk scale heights and beyond (0.2 kpc ≲ |z| ≲ 2.5 kpc) is somewhat underpredicted in their model compared to observations. They suggested that additional components of pressure from the magnetic field and cosmic rays may contribute significantly to the support in the vertical direction (see Section 3.1 of their paper). Additional vertical pressure support will lead to larger disk scale heights and hence larger gas masses and lower temperatures at 1–2 kpc heights, which will reduce the discrepancy between the observations and the SN-driven ISM model. Preliminary results from magnetohydrodynamics simulations are in agreement with this expectation (A. Hill et al. 2010, in preparation).

The second possible cause for the discrepancy between the predicted and observed X-ray temperatures is that the hot gas is out of equilibrium. This gas was shock-heated by SNe—this rapid heating causes the ionization temperature (which determines the X-ray spectrum) to lag behind the kinetic temperature (which is the quantity obtained from the hydrodynamical simulation). This ionization evolution was followed self-consistently in the one-dimensional SNR models described in Section 4.2, but not in this three-dimensional simulation. CIE is reached on a timescale of t_eq ∼ 10¹²n⁻¹_e s, where n_e is the electron density in cm⁻³ (Masai 1994). In the simulation, the density of gas with T>10⁶ K is ≲10⁻³ cm⁻³, implying t_eq ≳ 30 Myr. This is similar to the dynamical timescale (5 kpc/c_s ∼ 30 Myr, where c_s is the sound speed in 10⁶ K gas), and so we would expect the hot gas to be at least partially underionized.

Calculating the X-ray spectrum of an underionized plasma requires knowledge of the ionization balance, which in turn depends on the history of the plasma. However, we note that an underionized plasma will have fewer higher-stage ions (e.g., O⁺⁷) relative to lower-stage ions (e.g., O⁺⁶) than we would expect from the kinetic temperature. Therefore, if the halo is underionized, modeling the observed emission with a CIE model will underestimate the kinetic temperature of the halo gas. We have confirmed this by simulating XMM-Newton observations of an underionized plasma (using the XSPEC nei model), and characterizing the resulting simulated spectra with 1T CIE models (for simplicity, here we just use the nei model as an input for our simulations, instead of the full multicomponent model used in Section 4.2.2). For n_et ≲ 10⁹ cm⁻³ s, CIE models do not fit the simulated spectra well (these models underestimate the low-energy flux). However, for n_et ∼ 10⁹–10¹² cm⁻³ s, 1T CIE models give good fits to the simulated spectra, but underestimate the input kinetic temperature.

The above discussion considers only the hot gas being underionized. An additional possible source of X-rays is delayed recombination from overionized gas—gas that has undergone rapid adiabatic cooling as it expands into the halo (Breitschwerdt & Schmutzler 1994, 1999). With our current simulation data, it is difficult to estimate the contribution of delayed recombination to the halo emission, relative to the emission from the hot gas. Simulations that can track the ionization evolution of the plasma are needed to determine the extent to which the assumption of CIE affects the predicted spectra, and hence the derived X-ray temperature.

The third possible cause of the temperature discrepancy is that something (other than the assumption of CIE) is biasing the temperatures measured from our XMM-Newton observations. The discrepancy is not because we characterize the halo emission with a 1T model, as opposed to a 2T model (e.g., Kuntz & Snowden 2000; see Section 5.1), because we characterize the model spectra in the same way. If we are underestimating the halo temperature, the most likely cause is that our foreground (LB and/or SWCX) model is too faint; in particular, that it underestimates the foreground O vii emission (underestimating the foreground O vii emission means that the halo O vii emission is overestimated relative to the halo O viii emission, causing the measured halo temperature to shift to a lower value). However, our halo temperatures are in good agreement with those from other recent studies (see Section 5.1), which used different methods to estimate the foreground emission. It therefore seems unlikely that a bias in the halo temperature measurements is causing the discrepancy between the observed temperatures and those predicted by the SN-driven ISM model. Nevertheless, an accurate model of SWCX emission will help evaluate our method for estimating the foreground emission.