A CENSUS OF BARYONS AND DARK MATTER IN AN ISOLATED, MILKY WAY SIZED ELLIPTICAL GALAXY

2.1.1. Data Reduction

2.1.2. Chandra Image

We examined the image for evidence of morphological disturbances that may complicate deprojection, or possibly suggest perturbation from hydrostatic equilibrium. In Figure 1, we show the smoothed, flat-fielded Chandra image, having removed the point sources with the algorithm outlined in Fang et al. (2009). The data were then flat fielded with the exposure map and smoothed with a Gaussian kernel. Since the S/N of the image varies strongly with off-axis angle, the width of the Gaussian kernel was varied with distance according to an arbitrary power law, ranging from ∼1'' in the center of the image to ∼1' at the edge of the field.

Zoom In Zoom Out Reset image size

The image is smooth and slightly elliptical, showing no clear evidence of disturbances or asymmetries. This is consistent with our findings based on a shallower Chandra observation (Buote et al. 2002). To search for more subtle structure in the image, we used dedicated software, built around the MINUIT library,⁷ to fit an elliptical beta model (with constant ellipticity) to the central ∼25 wide portion of the unsmoothed image, (correcting for exposure variations with the exposure map).⁸ We show in Figure 1 an indication of deviations from this simple model by plotting (data − model)²/model, which corresponds to the χ² residual in each pixel. To bring out structure, we smoothed this image with a Gaussian kernel of width 15 (3 pixels). As is clear from this "residuals significance" image, there is no evidence of any coherent residual features. This is unsurprising given the lack of currently active central black hole, as suggested by the lack of significant radio emission (Fabbiano et al. 1987), nor a central, bright X-ray source.

2.1.3. Spectral Analysis

We extracted spectra in a series of concentric, contiguous annuli, placed at the X-ray centroid (which was computed iteratively and verified by visual inspection). The widths of the annuli were chosen so as to contain approximately the same number of background-subtracted photons and ensure that there were sufficient photons to perform useful spectral fitting. The resulting annuli all had widths larger than ∼8'', which is sufficient to prevent the instrumental spatial resolution from producing strong mixing between the spectra in adjacent annuli. The data in the vicinity of any detected point source were excluded, as were the data from the vicinity of chip gaps, where the instrumental response may be uncertain. We extracted products from all the active chips (excluding the S4 CCD, which suffers from considerable noise), making appropriate count-weighted spectral response matrices for each annulus with the standard CIAO tasks mkwarf and mkacisrmf. We extracted identical products individually for each data set, taking into account the astrometric offsets between them (determined from the image registration discussed above). The spectra were added using the standard Heasoft task mathpha, and the response matrices were averaged using the Heasoft tasks addrmf and addarf, with weights based on the number of photons in each spectrum. Representative spectra, without background subtraction, are shown in Figure 2, along with the best-fitting source and background models.

Zoom In Zoom Out Reset image size

Spectral fitting was carried out in the energy band 0.5–7.0 keV, using Xspec vers. 12.5.1n. We have shown previously that fits to Poisson distributed data which minimize χ² can yield significantly biased results, even if we rebin the data to more than the canonical ∼20 counts per bin (Humphrey et al. 2009b). In contrast, we found that the C-statistic of Cash (1979) typically gives relatively unbiased results. We therefore performed the fits by minimizing C, but nonetheless took care to verify that there is no significant bias, using the Monte Carlo procedure outlined in Humphrey et al. (2009b). Although not strictly necessary for a fit using the C-statistic, we rebinned the data to ensure a minimum of 20 photons per bin. This aids convergence by emphasizing differences between the model and data. The data in all annuli were fitted simultaneously, to allow us to constrain the source and background components at the same time.

We modeled the source emission in each annulus as coming from a single APEC plasma model with variable abundances, which was modified by Galactic foreground absorption (Dickey & Lockman 1990). Unlike several of our recent studies of galaxies (Humphrey et al. 2006, 2008, 2009a; Humphrey & Buote 2010) we did not deproject the data, but instead fitted directly the projected spectra, similar to Gastaldello et al. (2007b). We adopted this procedure since we were interested in accurately tracing the gas density to the largest available radii; deprojection requires the accurate subtraction of the emission from beyond the outermost radius, which is challenging given the very flat surface-brightness profile of this galaxy. In our previous studies, our approximate method for achieving this introduced significant uncertainties into the outermost density data point, which we have found to be prohibitive here. Since the temperature profile is relatively isothermal, the emission-weighted, projected gas density and temperature can be computed fairly accurately with the "response weighting" algorithm outlined in Appendix B of Gastaldello et al. We discuss the impact on our results of fitting the deprojected data in Section 5.4.

To model the source spectrum, we adopted a slightly modified version of the Xspec vapec model, for which the abundance ratios of each species with respect to Fe are directly constrained. We allowed the abundance of Fe (Z_Fe), and the abundance ratios of O, Ne, Mg, Si, and Ni with respect to Fe to vary. The other abundances (He) or abundance ratios with respect to Fe (for the other species) were fixed at their solar values (Asplund et al. 2006). To improve the constraints, Z_Fe was tied between multiple annuli, where required, and the abundance ratios were tied between all annuli. The best-fitting abundances (Figure 3) are consistent (in the central part of the system) with previous Chandra results from the shallow (17 ks) observation (Humphrey & Buote 2006) and are competitive with those obtained from the very deep Suzaku observation (Section 2.2; Tawara et al. 2008) as well as the XMM-Newton and Chandra measurements of the gas in much more massive (and X-ray bright) objects (e.g., Humphrey & Buote 2006; Kim & Fabbiano 2004; Buote et al. 2003). Intriguingly, we note evidence of a significant abundance gradient in this system (discussed in more detail in Section 6.5).

Zoom In Zoom Out Reset image size

To account for emission from undetected point sources, we included an additional 7.3 keV bremsstrahlung component. Since the number of X-ray point sources is approximately proportional to the stellar light (e.g., Humphrey & Buote 2008), the relative normalization of this component between each annulus was fixed to match the relative K-band luminosity in the matching regions, which we measured from the Two Micron All Sky Survey (2MASS) image. To determine the total X-ray binary luminosity, we added the total luminosity of this component, summed over all annuli, to the total L_X of the detected sources within the D₂₅ ellipse (de Vaucouleurs et al. 1991).⁹ The resulting L_X, (2.95 ± 0.12) × 10⁴⁰ erg s^-1, is in excellent agreement with that inferred from extrapolating the point-source luminosity function ((2.8 ± 0.8) × 10⁴⁰ erg s^-1; Humphrey & Buote 2008).

To account for the background, we included additional spectral components in our fits, in a variant of the modeling procedure outlined in H06, which were fitted simultaneously to each spectrum, along with the source. To account for the sky background, we included two (unabsorbed) APEC components (kT = 0.07 keV and 0.2 keV) and an (absorbed) power-law component (Γ = 1.41). The normalization of each component within each annulus was assumed to scale with the extraction area, but the total normalizations were fitted freely. (We discuss the likely impact of the "solar wind charge exchange" component in Section 5.2.1.) To account for the instrumental background, we included a number of Gaussian lines and a broken power-law model, which were not folded through the ancillary response file (ARF). We included separate instrumental components for the front- and back-illuminated chips and assumed that the normalization of each component scaled with the area of the extraction annulus which overlapped the appropriate chips. The normalization of each component and the shape of the instrumental components were allowed to fit freely. We included two Gaussian lines (at 1.77 and ∼2.2 keV), the intrinsic widths of which were fixed to zero. The energies of the ∼2.2 keV lines were allowed to fit freely, as were the normalizations of all the components.

The best-fitting models are shown in Figure 2 for a representative selection of spectra. To verify that the fit had not become trapped in a local minimum, we explored the local parameter space by stepping individual parameters over a range centered around the best-fitting value (analogous to computing error bars with the algorithm of Cash 1976). Error bars were computed via the Monte Carlo technique outlined in Humphrey et al. (2006), and we carried out 150 error simulations. In our previous work these simulations were used to generate the standard error on each measured temperature and density data point, which were assumed to be uncorrelated. In a low-temperature (≲1 keV) system, however, there is a strong degeneracy between the abundance and the gas density. Tying the abundance between multiple annuli, as done in the present work, will consequently result in correlated density errors. Therefore, rather than computing only the error on each data point from the error simulations, we also computed the covariance between each density data point, which we fold into our mass modeling (see Section 3). (In Section 5.8, we show that neglecting this effect, or adopting a full treatment that incorporates covariance between all the temperature and density data points, results in slightly different error bars and a slight bias, but does not strongly influence our conclusions.) In Figure 3, we show the derived gas temperature, abundance and projected density profiles, and the 1σ error bars on each data point (assuming that the data points are uncorrelated). The "projected density" is the mean gas density if all the emission measured in the annulus originates from a region defined by the intersection of a cylindrical shell and a spherical shell, both of which have the same inner and outer radii as the annulus. Defining it in this way is convenient as it has the units of density, but its precise definition is unimportant, so long as care is taken to compute the models appropriately.

As is clear from Figure 2, the instrumental component dominates the background. This reflects both the strong non-X-ray background of the telescope, but also the ability of Chandra to mitigate the cosmic component by resolving a significant fraction of it into individual point sources, which were excluded from subsequent analysis. Nevertheless, below ∼1 keV, the emission from the ∼0.5 keV gas in NGC 720 dominates at least out to ∼60 kpc, allowing the gas temperature and density to be constrained. Emission from unresolved point sources associated with NGC 720 is important only in the innermost few bins, and (as is clear from Figure 2) it can be clearly disentangled from the hot gas component, as it has a very different spectral shape. We note that we have not explicitly included a component to account for the combined emission from cataclysmic variables and hot stars (Revnivtsev et al. 2008). This component, however, is expected to be an order of magnitude fainter than the X-ray binary component, and so this omission will not affect our conclusions (Humphrey et al. 2009a).

The observation of NGC 720 with Suzaku was made early in the life of the satellite, while all four XIS units were operating. Data reduction was performed using the Heasoft 6.8 software suite, in conjunction with the XIS calibration database (Caldb) version 20090925. To ensure up-to-date calibration, the unscreened data were re-pipelined with the aepipeline task and analyzed following the standard data-reduction guidelines.¹⁰ Since the data for each instrument were divided into differently telemetered event file formats, we converted the "5 × 5" formatted data into "3 × 3" format and merged them with the "3 × 3" event files. The light curve of each instrument was examined for periods of anomalously high background, but no significant amount of data was found to be contaminated in this way. In Figure 1, we show the 0.5–7.0 keV image for the XIS0 image, excluding data in the vicinity of the calibration sources. A number of off-axis point sources, most likely background AGNs, were identified by visual inspection. To minimize their contamination, we excluded data in a 15 radius circular region centered on each of these sources. Given the large telescope point-spread function, such a cut will only eliminate ∼70% of the photons from these sources. However, we have chosen not to use a larger exclude region, since it will cut out a large portion of the field of view while only slightly reducing the contamination.

Spectra were extracted in three concentric annuli (0'–2', 2'–6', and 6'–10'), centered at the nominal position of NGC 720 in the field of view of each instrument. Data in the vicinity of the calibration sources and the identified point sources were excluded. For each spectrum, we generated an associated redistribution matrix file (RMF) using the xisrmfgen tool and an estimate of the instrumental background with the xisnxbgen task. Since the point-spread function of the telescope is very large (∼2' half-power diameter), even with such large apertures it is necessary to account for spectral mixing between each annulus. We did this by generating multiple ARFs for each spectrum with the xissimarfgen tool (Ishisaki et al. 2007), which models the telescope's optics through ray tracing. Formally, we wish to model the true spectrum of the hot gas, T^gas(α, δ, E), which is a function of right ascension (α), declination (δ), and energy (E). To make the problem tractable, it is a usual practice to adopt an approximation of the form T^gas(α, δ, E) ≃ ∑_ia^gas_i(α, δ)t^gas_i(E)ξ^gas(α, δ), where a^gas_i(α, δ) simply delineates discrete, non-overlapping regions of the sky (i.e., a^gas_i(α, δ) = 1 if (α, δ) is within region i, or 0 otherwise), and we have assumed that, within region i, the spectrum can be written as t^gas_i(E)ξ^gas(α, δ), where t^gas_i(E) is the normalized spectrum and ξ^gas is the normalization as a function of position on the sky. This is analogous to expanding T^gas(α, δ, E) on a set of basis functions that are separable in E and (α, δ). For our purposes, each i refers to an annular region on the sky (with outer radii 2', 6', and 10'), to correspond to the extraction regions. We approximated ξ^gas as a β-model (with parameters chosen to fit the Chandra image), correctly normalized to account for the gas emissivity. After passing through the telescope optics, a photon will be detected by the CCDs at some detector coordinate (x,y) and pulse height (h; corresponding to the energy). The pulse-height spectrum P(x, y, h) can be written as (e.g., Davis 2001)

$$\begin{eqnarray*} &&P(x,y,h) = \int d\alpha \ d\delta \ {\it dE} K(x,y,h,\alpha,\delta,E) T(\alpha,\delta,E), \end{eqnarray*}$$

where K encapsulates the effects of the telescope optics and the instrumental response. Now, we define D_j(h) as the pulse-height spectrum accumulated in an extraction region, j, such that D_j(h) = ∫dxdyb_j(x, y)P(x, y, h). (Typically we would choose i and j to correspond to the same set of regions on the sky, but that is not essential.) Here, b_j(x, y) is 1 if (x, y) lies within region j, or 0 otherwise. Thus, we can write

$$\begin{eqnarray*} D_j & = & \sum _i \int {\it dE}\ t^{\rm gas}_i \left[ \int dx\ dy\ d\alpha \ d\delta \ a^{\rm gas}_i b_j \xi ^{\rm gas} K \right]. \end{eqnarray*}$$

Now, we assume that the term in square brackets can be replaced by a product of two new terms, R_j(E, h) and A^gas_ij(E). With this definition, these new terms correspond, respectively, to the redistribution matrix and ancillary response matrix that can be created with the tasks xisrmfgen and xissimarfgen. Hence,

$$\begin{eqnarray} &&D_j (h) = \sum _i \int {\it dE}\ t^{\rm gas}_i(E) R_{j}(E,h) A^{\rm gas}_{ij}(E). \end{eqnarray} \tag{ 1 }$$

Spectral fitting involves adopting a parameterized form for the spectrum t_i(E) and comparing the corresponding model D_j(h) to the observed pulse-height spectrum. Unlike Chandra, for which the mixing between annuli was negligible (so that A^gas_ij(E) = 0 if i ≠ j), it was necessary to include all the ARFs A^gas_ij in the Suzaku fitting, which can be done with Xspec. The spectra from all extraction regions j were fitted simultaneously, to allow all three t_i(E) models to be constrained. We note that the Xspec suzpsf model can be used to perform a more approximate correction for the mixing between annuli. However, our approach more formally accounts for the spatial distribution of the source photons and is more practical for our present purposes. The spectra from all instruments were fitted simultaneously, allowing additional multiplicative constants in the fit to enable the relative normalization of the model for each instrument to vary with respect to the XIS 0.

It is necessary to add additional components to account for unresolved LMXBs and for the background components. The ARF terms (A_ij(E)) depend not only on the extraction region geometry but also the spatial distribution of the source photons, which clearly differs between the gas, LMXB, and background components. Therefore, it is necessary to generate additional ARFs for each component. For each of these new components, we assume that the spectral shape does not vary with α and δ, and so we can write T^lmxb(α, δ, E) = ξ^lmxb(α, δ)t^lmxb(E), and similar terms for the background. We adopt a 7.3 keV bremsstrahlung model for t^lmxb(E) and allow its normalization to be fitted freely. Since LMXBs should trace the stellar light, we used the K-band 2MASS image (appropriately normalized) to define the term ξ^lmxb. Folding all of the components into the fitting, Equation (1) is modified to

$$\begin{eqnarray} D_j(h) & = & \sum _i \int {\it dE}\ t^{\rm gas}_i(E) R_{j}(E,h) A^{\rm gas}_{ij}(E) \nonumber \\ & & + \int {\it dE}\ t^{\rm lmxb}(E) R_j(E,h) A^{\rm lmxb}_j(E) \nonumber \\ & & + \int {\it dE}\ t^{\rm sky}(E) R_j(E,h) A^{\rm sky}_j(E) \nonumber \\ & & + \int {\it dE}\ t^{\rm inst}(E) \overline{R}_j(E,h). \end{eqnarray} \tag{ 2 }$$

Here, we consider separately the sky and instrumental ("inst") background contributions. For the former, we adopted two APEC plasma models plus a power-law term, as described in Section 2.1.3 and assume that its surface brightness does not vary over the field of view when computing A^sky_j(E). For the latter term, we adopted two, complementary, approaches. One is to subtract off the spectrum generated for each annulus by xisnxbgen, in which case the last term in Equation (2) will be omitted. The other method, which we adopted by default as it provides greater generality in our fits, is to model the instrumental background, similarly to the treatment of this component with Chandra. We found that a single broken power law, plus Gaussian lines (at ∼1.5, 1.7, 2.2, 5.9, 6.5, and 7.5 keV), was sufficient for our purposes here, provided that the model was not multiplied by the effective area curve. Unfortunately, the quantum efficiency is folded into the RMF created by xisrmfgen, and so we first had to normalize the RMF through which this model was folded (the "renormalized" response matrix being indicated by $\overline{R}_j$ in Equation (2)). Still, we found that the precise shape of the instrumental term is not important since it is a relatively minor component in the spectral model (except in the vicinity of the calibration lines). We did not obtain significantly different results if we used the modeling method or the direct subtraction method for the instrumental contribution (Section 5.2).

We found that our model was able to capture the shape of the spectra well. The best-fitting normalization of the LMXB component corresponds to an LMXB luminosity of (3.08 ± 0.07) × 10⁴⁰  erg s^-1, in excellent agreement with the Chandra data, and with the extrapolation of the detected LMXB luminosity function. Representative spectra are shown in Figure 2, showing the various model components. It is immediately clear that the spectral mixing between annuli is significant. For example, ∼20% of the photons from the central bin are scattered into the second annulus. Failure to take account of this effect actually results in a significantly different temperature profile, which is almost isothermal at ∼0.5 keV and very inconsistent with the Chandra data. The background is dominated by the sky component, which is more significant than in Chandra since very few of the point sources that dominate it above ∼1 keV can be resolved and excluded.

The hot gas temperature, abundance, and projected density profiles are shown in Figure 3. In the two outer annuli, the data are in good agreement with the results from Chandra. In the innermost annulus, however, we found some discrepancies. Specifically, the abundance is significantly higher than observed with Chandra, while the temperature is slightly higher than an emission-weighted average of the Chandra temperature profile in this region. The Suzaku results in the central bin are, however, sensitive to systematic uncertainties in our modeling procedure (see Section 5). In particular, if we omitted data from the XIS1 unit from our fits, we obtained Z_Fe in much better agreement with Chandra (Figure 3). While this may point to problems with the calibration of the XIS1 unit, it may also be related to errors in the correction we applied to account for spectral mixing between the different annuli. If we omit this correction, for example, the discrepancy with Chandra is significantly reduced (Z_Fe ∼ 0.9). The real temperature profile is clearly more complex than we assume when performing this correction, which may be partially responsible, as there may be slight errors in the ray-tracing algorithm. Furthermore, the single-temperature model fitted to the large Suzaku regions, within which the Chandra data reveal temperature gradients, would not be expected to yield perfectly consistent results with the Chandra data. In any case, we omit the data from the central Suzaku annulus in our subsequent modeling, since the Chandra data are expected to be more reliable on these scales.

We first consider the case without a DM halo component. Although our previous work (e.g., H06) has found that the M/L ratio rises steeply with radius (implying that DM is needed at high significance), given the high-quality data in NGC 720, it is of particular interest to quantify how poorly the constant M/L ratio model fits these data. We find that the best-fitting model without DM is a very poor description of the temperature and density profiles (Figure 3; χ²/dof = 413/19, with a mean absolute fractional deviation between the data and models of ∼29%). As found by Buote et al. (2002), based on an analysis of the X-ray image in NGC 720, the stellar mass component is too centrally concentrated to produce the observed gas pressure profile.

As is clear in Figure 3, when a DM halo component was added, the hydrostatic model was able to capture the overall shape of the projected density and temperature profiles very well (χ²/dof = 12.7/17). The improvement to the fit when dark matter was included was therefore highly significant; the ratio of the Bayesian "evidence" returned by the fitting algorithm (see Feroz & Hobson 2008) for the two cases is 10⁻⁸⁵, implying that DM is required¹⁶ at ∼19.6σ.

In Figure 5, we show the best-fitting entropy profile (and its 1σ confidence range), which shows a flattening at large radius similar to two of the three galaxies studied in Humphrey et al. (2009a). As we will show later (Section 6.4), the entropy profile shape, even when extrapolated to large radius, is consistent with trends seen in galaxy clusters (Pratt et al. 2010). The derived radial distribution of the gravitating mass is shown in Figure 6 and the total mass-to-light ratio is shown in Figure 7. In Figure 6, we also show the contribution to the mass from each separate component. We find that the stars dominate the mass distribution within the central ∼5 kpc, but a significant DM halo is required to reproduce the data at larger radii, consistent with our previous study of this system, using data from a much shallower observation (Humphrey et al. 2006). The hot gas becomes the dominant baryonic component outside ∼200 kpc (∼R₅₀₀). Overlaid are a series of mass data points derived from fitting the deprojected Chandra data using a more traditional (but more uncertain) mass-modeling method (for more details see Section 5.4 and Humphrey et al. 2009a), which agree well with the inferred mass distribution, indicating that the resulting mass profile is not overly sensitive to the analysis method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A notable feature of the recovered mass distribution is its flattening outside ∼20 kpc, which allows us to constrain the scale radius, and hence the virial mass of the DM halo. In Figure 8, we show the relation between the concentration of the gravitating mass, c_vir, and the virial mass, M_vir. To be consistent with our past work (Buote et al. 2007), the M_vir and R_vir are derived from the distribution of the total gravitating mass, not just the DM, and the concentration, c_vir, is defined as the ratio of R_vir to the characteristic scale of the DM halo. These results are consistent with those obtained by H06 with Chandra, but the constraints are very much tighter, reflecting both the significantly deeper observation (∼5 times deeper Chandra observation, combined with the data from the very deep Suzaku observation) as well as improvements in our mass-modeling procedure. Overlaid on Figure 8 we show the theoretical concentration versus mass relation from DM only simulations (Macciò et al. 2008, adopting their "ΛCDM1" relation for relaxed halos) and the empirical relation from Buote et al. (2007). NGC 720 lies above both relations, but is within the 2σ scatter obtained by Buote et al. The best-fitting marginalized model parameter values and confidence regions are given in Tables 2 and 3.

Zoom In Zoom Out Reset image size

In contrast to our measurement in Humphrey et al. (2006), our adopted modeling procedure did not require us to adopt a strong prior on f_b in order to obtain interesting virial mass constraints. This allows us to measure f_b directly from our fit. In Figure 9, we show how the enclosed baryon fraction and gas fraction vary as a function of radius and tabulate the gas fraction (f_g) and f_b, measured at different overdensities, in Table 4. Outside ∼100 kpc (∼R₂₅₀₀), our model predicts that the f_b profile rises gently, from ∼10% to 15%. In Figure 10, we show the posterior probability density function for $f_{b,\rm vir}$ (f_b extrapolated to R_vir), which is tightly constrained to around ∼15%, very close to the Cosmological value. A potentially important source of uncertainty in this extrapolation is the shape of the entropy profile. In Figure 10, we show the equivalent posterior probability density functions for two different entropy extrapolations, which should bound all likely entropy profiles. We discuss these fits in detail in Section 5.3. Clearly, uncertainties in the entropy extrapolation do not produce a large systematic shift in our inferred $f_{b,\rm vir}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Baryon Fraction Results and Error Budget

Test	fg,2500	fg,500	$f_{g,\rm vir}$	fb,2500	fb,500	$f_{b,\rm vir}$
Marginalized	0.027+0.003−0.004	0.06+0.009−0.02	0.12+0.05−0.03	0.10 ± 0.01	0.11+0.02−0.02	0.16 ± 0.04
Best Fit	(0.024)	(0.053)	(0.119)	(0.096)	(0.102)	(0.154)
Systematic error estimates
ΔDM profile	−0.004 (±0.003)	−0.024 (±0.01)	−0.085 (+0.01−0.01)	−0.013 (±0.01)	−0.044 (±0.01)	−0.114 (±0.01)
ΔAC	−0.001 (±0.004)	−0.010 (+0.01−0.01)	−0.019 (+0.05−0.03)	−0.026 (±0.01)	−0.023 (±0.02)	−0.035 (±0.04)
ΔBackground	−0.003 (±0.003)	−0.001 (±0.01)	+0.03 (±0.04)	−0.005 (+0.01−0.01)	±0 (±0.02)	+0.02 (+0.04−0.03)
ΔSWCX	±0 (±0.002)	+0.001 (+0.01−0.01)	+0.03 (+0.03−0.04)	+0.003 (±0.01)	+0.005 (±0.01)	+0.03 (+0.03−0.04)
ΔEntropy	±0 (±0.004)	−0.018 (±0.01)	+0.02−0.04 (±0.05)	±0 (±0.01)	+0.01−0.01 (±0.02)	+0.01−0.05 (±0.05)
Δ3d	+0.005 (+0.004−0.005)	+0.010 (±0.02)	+0.04 (+0.05−0.08)	+0.01 (±0.01)	+0.01 (±0.03)	+0.03 (+0.06−0.08)
ΔSpectral	+0.001−0.008 (±0.003)	−0.026 (±0.01)	+0.02−0.07 (±0.05)	+0.00−0.03 (+0.01−0.01)	+0.01−0.04 (±0.02)	+0.03−0.09 (±0.05)
ΔWeighting	−0.001 (±0.002)	±0 (±0.01)	+0.01 (±0.05)	−0.006 (±0.01)	−0.003 (±0.01)	+0.008 (±0.05)
ΔFit priors	+0.007−0.004 (±0.004)	+0.01−0.02 (±0.02)	+0.02−0.04 (+0.06−0.05)	+0.01−0.01 (±0.01)	±0.02 (±0.02)	+0.01−0.06 (±0.06)
ΔInstrument	−0.008 (±0.005)	−0.025 (±0.02)	−0.067 (±0.06)	−0.021 (±0.01)	−0.032 (±0.02)	−0.081 (±0.05)
ΔStars	±0 (±0.004)	−0.006 (+0.02−0.01)	+0.008 (±0.05)	+0.02−0.01 (+0.03−0.05)	−0.008 (+0.03−0.02)	+0.01−0.01 (+0.05−0.06)
ΔDistance	−0.002 (+0.004−0.003)	−0.007 (+0.02−0.01)	+0.001 (+0.05−0.04)	+0.00−0.00 (±0.01)	−0.006 (+0.02−0.02)	−0.006 (±0.05)
ΔFit radius	−0.007 (+0.007−0.005)	−0.019 (±0.02)	−0.064 (+0.07−0.04)	−0.017 (±0.02)	−0.032 (+0.03−0.02)	−0.076 (+0.08−0.04)
ΔCovariance	−0.002 (+0.004−0.003)	−0.005 (±0.01)	+0.02 (+0.04−0.06)	+0.002 (±0.01)	±0.00 (±0.02)	+0.003 (±0.05)

Notes. Marginalized values and 1σ confidence regions for the gas fraction (f_g,Δ) and baryon fraction (f_b,Δ) measured at various overdensities (Δ). We also provide the best-fitting parameters in parentheses, and a breakdown of possible sources of systematic uncertainty, following Table 3. We find that f_b is reasonably robust to most sources of systematic uncertainty, especially within R₂₅₀₀.

Download table as:  ASCIITypeset image

One of the major sources of uncertainty on the recovered f_b is the choice of DM mass model. While our default model (NFW) is theoretically motivated, and supported by observations of galaxy clusters (e.g., Vikhlinin et al. 2006), we also experimented with a so-called cored logarithmic mass model, which has been widely used in studies of elliptical galaxies (Binney & Tremaine 2008), including some recent stellar dynamical studies (e.g., Thomas et al. 2007b; Gebhardt & Thomas 2009). This model tends to predict higher masses at large radii than NFW, resulting in a significantly larger M_vir, and correspondingly a smaller value of f_b when extrapolated within R_vir. Nevertheless, at radii corresponding to higher mass overdensities, the extrapolation range is smaller and the discrepancy is significantly reduced, so that the choice of DM model has little effect on f_b,2500 ("ΔDM profile" in Tables 3 and 4). In any case, we note that the cored logarithmic model is less well-motivated theoretically than the NFW profile. Furthermore, we find that the data slightly favor the NFW model; although the best-fitting χ² is comparable between the two cases, the ratio of the Bayesian evidence is ∼1.5 × 10⁻⁴, implying that the cored logarithmic model, with the adopted priors (a flat prior on the asymptotic circular velocity, between 10 and 2000 $\rm km\ s^{-1}$, and a flat prior on log₁₀r_c, where r_c is the core radius, over the range 0 ⩽ log₁₀r_c ⩽ 3), is a poorer description of the data at ∼3.8σ.

One modification to the (NFW) DM profile shape that is theoretically predicted is "adiabatic contraction" (Blumenthal et al. 1986; Gnedin et al. 2004), where the DM halo density profile reacts to the gravitational influence of baryons that are condensing into stars by becoming cuspier. There is little observational evidence of this effect (e.g., H06; Humphrey et al. 2009a; Gnedin et al. 2007; although, see Napolitano et al. 2010), and it may be overestimated by existing models (Abadi et al. 2010), so we did not employ it in our default analysis. Nevertheless, we experimented with modifying the NFW halo profile to account for this effect with the algorithm described by Gnedin et al. (2004).¹⁷ The principal effect of this modification was to lower slightly the stellar M/L ratio and, consequently, the inferred f_b ("ΔAC" in Tables 3 and 4).

Our treatment of the background is a potentially serious source of systematic error. To investigate this, we adopted an alternative approach in which we used, for the Chandra data, the standard blank-field event files distributed with the CALDB to extract a background spectrum for each annulus. Since the blank-field files for each CCD have different exposures, spectra were accumulated for each CCD individually, scaled to a common exposure time and then added. The spectra were renormalized to match the observed count rate in the 9–12 keV band. These "template" spectra were then used as a background in Xspec, and the background model components were omitted from our fit. We found that this approach gave a poorer fit to the data than the full modeling procedure used by default, and the overall density and temperature profiles were noticeably affected at large radii. For the Suzaku data, there are no blank-field files available, but we used the standard non-X-ray background spectra generated by xisnxbgen tool to remove the instrumental background, but maintained the sky background model components in our fit. Fitting the resulting Chandra and Suzaku temperature and density profiles with our mass-modeling apparatus, our results were not substantially affected ("ΔBackground" in Tables 3 and 4).

5.2.1. Solar Wind Charge Exchange

The extrapolation of the entropy profile is one of the prime sources of uncertainty in our determination of f_b at large scales. As we will show in Section 6.4, our default profile is, in fact, consistent with trends seen in massive galaxy clusters (Pratt et al. 2010), giving us confidence in its use. Nevertheless, it is important to explore plausible alternative shapes. Motivated by observed entropy profiles in galaxy groups and by the monotonically rising entropy profiles required by hydrostatic equilibrium, we investigated two pathological extrapolations that should bound any plausible model. Specifically, we adopted a power-law shape (entropy ∝R^γ) outside ∼80 kpc, with γ = 0 (i.e., a flat entropy profile), which is a secure lower boundary assuming stability against convection, or γ = 2.5. This latter value was chosen arbitrarily and is far larger than has been observed in real systems. This choice primarily affects the extrapolated temperature, rather than the density profile, and so f_b does not change substantially with these choices ("ΔEntropy" in Tables 3 and 4; Figure 10).

In the present work, we opted to fit the projected, rather than the deprojected data (as done, for example, in H06). In general, fitting the projected data leads to smaller statistical error bars, but potentially larger systematic uncertainties (e.g., Gastaldello et al. 2007b), and so it is important to investigate the likely magnitude of such errors. To do this, we examined the effect on our results of spherically deprojecting the data. We achieved this by using the Xspec projct model.¹⁹ To account for emission projected into the line of sight from regions outside the outermost annulus, we added an APEC plasma model to each annulus, with abundance 0.3 (consistent with the outermost annulus), and the temperature and normalization determined from projecting onto the line of sight the best-fitting gas temperature and density models described in Section 4.2, but considering the models only outside ∼70 kpc.

Although the error bars were increased when using the deprojected (rather than the projected) profiles ("ΔDeprojection" in Tables 3 and 4), the best-fitting results were not significantly changed. In Figure 6, we show a series of mass "data points" obtained from the deprojected density and temperature profiles, using the "traditional" mass analysis method described in Humphrey et al. (2009a, see also Humphrey & Buote 2010). These data agree very well with the best-fitting mass model found in our projected analysis. Similarly, the entropy profile derived directly from the deprojected data agrees well with the results of our projected analysis (Figure 5).

Since the choice of priors on the various parameters is arbitrary in our analysis, it is important to determine to what extent they could affect our conclusions. To do this, we replaced each arbitrary choice in turn with an alternative, reasonable prior. Specifically, for each parameter describing the entropy profile, we switched from a flat prior on that parameter to a flat prior on its logarithm. We replaced the flat prior on the stellar M/L ratio with a Gaussian prior, the mean and σ of which were determined from the SSP model results reported in H06. We used a flat prior on the DM halo mass, rather than on its logarithm, and, instead of the flat prior on log c_DM, we adopted the distribution of c around M found by either Buote et al. (2007) or Macciò et al. (2008) as a (Gaussian) prior. The effect of these choices is no larger than the statistical errors on each parameter, especially for the baryon fraction measured at R₂₀₀ or higher overdensities ("ΔFit priors" in Tables 3 and 4).

As discussed in H06, careful treatment of the stellar light is essential for obtaining an accurate measurement of the gravitating mass profile. Since the stellar light deprojection is not unique, in particular as the inclination angle is not a priori known for an elliptical galaxy, we investigated the associated systematic uncertainty by first varying it within reasonable limits (Section 3). Specifically, we lowered it as far as 70°, since the highly elliptical projected isophotes of NGC 720 make it unlikely to be observed at a much lower inclination. As a more extreme test of the influence of the stellar light on our fits, we also experimented with excluding all the data within the central ∼4 kpc, where the stars dominate the mass profile. We find that the stellar M/L ratio becomes (unsurprisingly) much more uncertain when we adopt this approach, but the best-fitting masses and f_b were not substantially affected ("ΔStars" in Tables 3 and 4).

In our default analysis, the projected temperature and density profile are weighted by the gas emissivity, folded through the instrumental responses (for details see Appendix B of Gastaldello et al. 2007b). Since the computation of the gas emissivity assumes that the three-dimensional gas abundance profile is identical to the projected profile (which is unlikely to be true), it is important to assess how sensitive our conclusions are to the emissivity correction. To do this, we adopted the extreme approach of ignoring the spatial variation of the gas emissivity altogether. We found that this had a very small effect on our results (ΔEmissivity in Tables 3 and 4).

We here outline the remaining tests we carried out, as summarized in Tables 3 and 4. First of all, since the inter-calibration of the Suzaku XIS units may not be perfect, we experimented with using only one of the units in the Suzaku analysis and cycled through each choice. We found that the XIS1 (back-illuminated) unit appeared most inconsistent with the Chandra data, and so also investigated using only units 0, 2, and 3. We also considered using only the Chandra data. While these choices had a significant effect on $f_{b,\rm vir}$, we found that f_b at higher overdensities was more resilient ("ΔInstrument").

To assess the impact of various spectral-fitting data analysis choices on our results, in turn we varied the neutral hydrogen column density by ±25%, performed the fit over a restricted energy range (0.7–7.0 keV or 0.5–4.0 keV), and replaced the APEC plasma model with a MEKAL model. The impact of these choices is comparable to the statistical errors on the parameters ("ΔSpectral").

To examine the error associated with distance uncertainties, we varied the distance to NGC 720 by the 1σ error given in Tonry et al. (2001). The effect on our results is minor ("ΔDistance"). To examine the sensitivity of our fits to the radial range studied, we excluded the data outside ∼40 kpc. This exclusion made the extrapolation more uncertain, especially for $f_{b,\rm vir}$ ("ΔFit radius"). Finally, to examine the possible errors associated with our treatment of the covariance between the density data points, we investigated adopting a more complete treatment that considers the covariance between all the temperature and density data points as well as adopting the more standard (but incorrect) approach of ignoring the covariance altogether. The effect is not large ("ΔCovariance" in Tables 3 and 4; Figure 8).

The ability of our hydrostatic model to reproduce the temperature and density profiles of the gas in NGC 720 strongly suggests that the gas is close to hydrostatic; despite highly nontrivial temperature and density profiles, a smooth, physical mass model and a monotonically rising entropy profile (required for stability against convection) were able to reproduce them well. If the hydrostatic approximation is seriously in error, this would require a remarkable conspiracy between the temperature and density data points. The closeness of the system to hydrostatic is unsurprising given its remarkably relaxed X-ray morphology (Figure 1; Buote et al. 2002; Buote & Canizares 1996, 1994) and lack of substantial radio emission (Fabbiano et al. 1987). Numerical simulations of structure formation suggest that deviations from hydrostatic equilibrium in morphologically relaxed-looking systems are not large, so that errors in the recovered mass distribution should be no larger than ∼25% (e.g., Tsai et al. 1994; Buote & Tsai 1995; Nagai et al. 2007; Piffaretti & Valdarnini 2008; Fang et al. 2009). Such a level of non-thermal support is comparable to the ∼10%–20% discrepancy found between X-ray and stellar dynamical mass measurements by Churazov et al. (2008) for two galaxies that are manifestly more disturbed than NGC 720.

In the case of NGC 720, there is additional, albeit circumstantial, evidence which supports the hydrostatic approximation. First, the best-fitting virial mass is in good agreement with the modest (although uncertain) velocity dispersion of the dwarf companions about NGC 720 (117 ± 54 km s^-1, corresponding to a virial mass of (4 ± 4) × 10¹² M_☉; Brough et al. 2006). In our studies of other systems, we have used the good correspondence between the measured stellar M/L ratio and that predicted by single-burst SSP models that assume a Kroupa (2001) IMF, to infer that there is little non-thermal pressure (Humphrey et al. 2009a). In the case of NGC 720, the measured M/L_K ratio (0.54 ± 0.05) is actually higher than the predictions of the SSP models (0.35 ± 0.07; H06). A natural way to bring these results into agreement is if the young (∼3 Gyr; Humphrey & Buote 2006; Rembold et al. 2005) stellar population that dominates the light only represents recent star formation superimposed on an older, underlying stellar population with a higher M/L_K ratio. Such a picture is supported by the multiple-aged stellar population inferred by Rembold et al. (2005). We note that the measured M/L_K value could also be reconciled with the predictions of SSP models if a Salpeter IMF is adopted (for which M/L_K∼ 0.54; taken in the context of our previous work this would indicate a non-universal IMF). Alternatively, if AC operates, making NGC 720 unique amongst the systems we have studied in the X-ray, the measured M/L_K would be 0.39 ± 0.04, also reconciling the fit result and the SSP predictions. While it is evidently difficult to draw strong conclusions based on the M/L_K ratio, none of these arguments indicates that the measured M/L_K is obviously underestimated. This is important since most models predict that non-hydrostatic effects lead to an underestimate of the mass (e.g., Nagai et al. 2007; Churazov et al. 2008; Fang et al. 2009).

Humphrey & Buote (2010) studied NGC 720, as part of a sample of galaxies, groups, and clusters, finding that the total gravitating mass profile in the central region (≲30 kpc) is close in shape to a singular isothermal sphere distribution. This shape is what is approximately inferred from lensing plus stellar dynamics studies of similar galaxies (e.g., Koopmans et al. 2009). Therefore, any deviations from hydrostatic equilibrium are unlikely to have substantially affected the overall shape of the gravitating mass profile, at least over this region. Nevertheless, it is plausible that a constant fraction of the pressure over this region comes from non-thermal effects, which would just reduce the overall normalization of the measured mass model. Such a scenario would imply that the true R_vir is underestimated, while the scale radius is approximately correct; hence c_vir would be underestimated, putting it in more tension with theory. Conversely, it is unlikely that deviations from hydrostatic equilibrium could lead to the scale radius also being underestimated, since that would require the non-thermal pressure support to increase dramatically outside ∼20 kpc. Since NGC 720 is an isolated system, it is difficult to imagine what processes could give rise to such an effect.

One effect which could produce deviations from hydrostatic equilibrium at the smallest scales is gas rotation induced by spin-up of subsonically inflowing gas. Considering the very relaxed galaxy NGC 4649, Brighenti et al. (2009) demonstrated that angular momentum conservation (if there is appreciable rotation in the stars) flattens the X-ray ellipticity significantly in the center of the system (where deviations from hydrostatic equilibrium are also most significant). NGC 720 exhibits sufficient major-axis rotation (Binney et al. 1990), but only a hint of a central ellipticity rise within the innermost ∼1 kpc (Buote et al. 2002). Since our results are not very sensitive to the exclusion of data within the central ∼4 kpc (Section 5.8), it seems unlikely that this effect could have an impact on our results, even if it is operating in NGC 720.

Our conclusions on the state of hydrostatic equilibrium in NGC 720 differ significantly from those of Diehl & Statler (2007). Those authors investigated the X-ray and optical isophotal ellipticities at a fixed, small scale in an heterogenous sample of early-type galaxies, concluding that the lack of a correlation between them (in conflict with would be expected if hydrostatic equilibrium holds exactly) indicates that hydrostatic equilibrium is not ubiquitous and, therefore, should never be assumed (even approximately) in mass analysis (and, in particular, for NGC 720). We consider their conclusions, however, to be extreme and misleading. First, their method does not allow them to say anything at all about individual objects, and only indicates that hydrostatic equilibrium is unlikely to hold perfectly in a galaxy randomly drawn from their sample, provided that it is chosen without any consideration of its morphology. Since their sample contained galaxies with such large-scale asymmetries that we would not recommend the routine application of hydrostatic mass methods (e.g., NGC 4636, Jones et al. 2002 and M 84, Finoguenov & Jones 2001), and since they focused on the smallest scales (where X-ray point-source removal is most challenging and where AGN-driven disturbances are most serious), in contrast to the large scales that are most important for DM analysis, the implications for morphologically relaxed systems such as NGC 720 or NGC 4649 (Humphrey et al. 2008) are unclear. Second, even if it does not hold perfectly, hydrostatic equilibrium can still be a useful approximation. With their hydrodynamical models, Brighenti et al. (2009) reconciled the X-ray and optical ellipticity profiles of NGC 4649 by requiring modest gas motions, but the gas remained very close to hydrostatic (especially outside the innermost ∼1 kpc). Substantial non-hydrostatic gas motions are likely to manifest themselves as sharp features in the surface brightness and temperature profiles of the hot gas, for example, those seen in the core of M 87 (Forman et al. 2007; Million et al. 2010). Still, even for M 87, Churazov et al. (2008) found that the gas was close to hydrostatic away from the regions of most significant disturbance. In the case of NGC 720, the profiles are significantly smoother than in M 87 (Figure 3), while the hydrostatic isophotal analysis undertaken by Buote et al. (2002) found evidence for a DM halo ellipticity in excellent agreement with the predictions of cosmological models (implying that the gas is not grossly out of equilibrium). In summary, in our detailed analysis of NGC 720 and other systems (Humphrey et al. 2009a; Brighenti et al. 2009), we find no evidence supporting Diehl & Statler's contention that hydrostatic equilibrium is a very poor approximation in morphologically relaxed galaxies (which seems little more than "guilt by association").

6.2.1. Robustness of f_b Measurement

6.2.2. Implications

Based on our hydrostatic mass models, the total stellar mass fraction within R_vir is small (∼0.035 ± 0.003), but in good agreement with that generally inferred from gravitational lensing studies of early-type galaxies (Hoekstra et al. 2005; Gavazzi et al. 2007). Most of the baryons within R_vir are, however, in the form of the hot gas. Including both gas and stars, the measured f_b profile is relatively flat at ∼0.10–0.15 between R₂₅₀₀ and R_vir (Figure 9), becoming completely consistent with the Cosmological value (0.17; Dunkley et al. 2009) when extrapolated to R_vir. In H06, we reported a lower value of f_b within R_vir for NGC 720 (∼0.04), based on simpler hydrostatic model fits to data from a much shallower Chandra observation. That constraint, however, was almost entirely driven by a strong prior that imposed a tight correlation between M_vir and f_b; when relaxing this prior, we found that f_b was very poorly constrained by those data.

A number of authors have investigated the gas (or baryon) fraction in galaxy groups and clusters, typically finding a shallow dependence on the virial mass of the halo. In the upper panel of Figure 11, we show the gas fraction measurements within R₂₅₀₀ (f_g,2500) versus M₂₅₀₀ for the group and cluster samples of Vikhlinin et al. (2006), Gastaldello et al. (2007b), and Sun et al. (2009).²⁰ Overlaid is the locus of NGC 720, which is in good agreement with an extrapolation of the power-law trend inferred for the more massive systems. Fitting a model of the form

$$\begin{equation} \log _{10} f_{g} = A + B \log _{10} \left(\frac{M}{2\times 10^{14}\,M_\odot }\right) \end{equation} \tag{ 3 }$$

to these data, using a method that takes into account errors on both axes and intrinsic scatter in the y-direction (specifically, that employed in Humphrey & Buote 2008), we found A = −1.08 ± 0.03 and B = 0.28 ± 0.04, with an intrinsic scatter of 0.11 dex. In the center panel of Figure 11, we show a similar comparison between f_g,500 and M₅₀₀, similarly finding that the locus of NGC 720 is consistent with an extrapolation of the trends for the massive systems. Fitting a log-linear regression line (Equation (3)), we obtained A = −1.02 ± 0.02 and B = 0.15 ± 0.03, with an intrinsic scatter of 0.06 dex, in good agreement with a fit by Giodini et al. (2009) to a similar group and cluster sample, but omitting the NGC 720 data point. Since the temperature varies only by ∼1 order of magnitude for the almost 3 orders of magnitude mass range between the galaxies and clusters, hydrostatic equilibrium generally requires that massive systems have a more rapidly declining gas density profile than less-massive objects, and hence a more centrally peaked gas mass distribution, as compared to the distribution of DM. Given this trend, it is not surprising that lower-mass objects have a smaller fraction of their total gas mass within R₂₅₀₀, as indicated by the dependence of B on the overdensity.

Zoom In Zoom Out Reset image size

Also shown in Figure 11 (center panel) are a set of data points taken from Dai et al. (2010), who fitted stacked ROSAT All-Sky Survey (RASS) data of optically selected galaxy groups and clusters. These points appear quite discrepant with this relation (and the other data) below ∼3 × 10¹³ M_☉, where they fall sharply. This may indicate real differences in the properties of optically selected groups from those selected from X-ray surveys, but it may also reflect the significant systematic errors that are involved with stacking the RASS data (see Dai et al. 2007).

While the gas is the dominant baryonic component for the most massive clusters, the stellar-to-gas mass ratio is actually a strong function of M₅₀₀. In NGC 720, approximately half of the baryons within R₅₀₀ are in the stars. Accounting for both gas and stars, Giodini et al. (2009) obtained average f_b constraints within R₅₀₀ for a galaxy group and cluster sample. In the lower panel of Figure 11, we show their data as a function of M₅₀₀, along with the best power-law fit to their data, and the locus of NGC 720. It is immediately clear that NGC 720 is in good agreement with the extrapolation of the (shallow) trend they observe to low masses. Dai et al. (2010) also reported total baryon fractions for their galaxy group and cluster sample, finding rather lower values of f_b,200 in their lowest-mass systems (f_b,200 ∼ 0.04 at M₂₀₀ ∼3 × 10¹³ M_☉). This is significantly below f_b,200 observed in NGC 720 (which has M₂₀₀ an order of magnitude lower), but is comparable to its stellar mass fraction. Still, as discussed above, the gas fraction estimates from their work appear discrepant with studies of X-ray-selected groups and this may explain the inconsistency.

Our study of NGC 720 clearly demonstrates that it is possible for a system with M_vir as low as a few times 10¹² M_☉ to retain a massive, quasi-hydrostatic hot halo. The total baryon fraction, when extrapolated to the virial radius, is consistent with the Cosmological value, indicating that massive DM halos around galaxies could be a significant reservoir for some fraction of the "missing baryons" (e.g., Fukugita & Peebles 2006). The biggest spiral galaxies are believed to reside in halos of comparable mass, making it particularly interesting to compare the hot halo around NGC 720 to those predicted around spiral galaxies. Intriguingly, the total gas content of the system (gas fraction within R₂₀₀ of 0.08 ± 0.03) does agree well with recent numerical disk galaxy simulations at this mass range (e.g., Crain et al. 2010). In these simulations, most of the gas is in the hot phase, suggesting that the cool baryon content of the most massive spiral galaxies (i.e., considering only stars and cool gas) should be systematically lower than f_b measured in NGC 720.

McGaugh et al. (2010, see also Dai et al. 2010) has assembled f_b,500 estimates for a sample of massive disk galaxies, omitting any putative hot halo from the estimate. Indeed, f_b,500 in these systems does lie slightly below that in NGC 720. Still, the stellar mass fraction in these spiral galaxies is substantially higher than in NGC 720 and the modest discrepancy in f_b,500 leaves little room for a massive hot halo around the spirals. This appears inconsistent with the simulations and, at face value, might suggest very different star formation efficiencies between the different types of galaxy. One difficulty with this interpretation, however, is the young (∼3 Gyr) dominant stellar population in NGC 720, suggesting a recent formation epoch, possibly from a merger event involving massive disk galaxies (Rembold et al. 2005). Since the gas fraction can only be reduced during a merger (as it may be ejected from the galaxy, funneled onto the central black hole or converted into stars), this implies that there must have been at least as much gas as stars within R₅₀₀ in the progenitors. Even allowing for stellar mass loss over the last ∼3 Gyr, this would require similarly inefficient star formation in the precursor (spiral?) galaxies. We caution that the M₅₀₀ values for the spiral galaxies in McGaugh et al. were obtained by scaling the rotation velocity, rather than formal modeling, which could plausibly introduce some uncertainties into both M₅₀₀ and f_b,500 reported in that work.

Finally, it is interesting to compare our results to measurements of f_b in the Milky Way, around which there is indirect evidence for an extended, hot halo. This system has M_vir in the range ∼(1–3) × 10¹² M_☉ (Klypin et al. 2002; Sakamoto et al. 2003), close to that of NGC 720. Based on the baryonic mass estimate of Flynn et al. (2006), which does not consider any putative hot halo, $f_{b,\rm vir}$ lies in the range ∼0.02–0.06. This is suggestively below our measurement for NGC 720, and the discrepancy is most acute toward the more massive end of the Milky Way mass range (i.e., almost the same mass as NGC 720).

We have obtained a high-quality measurement of the gravitating mass profile for the isolated elliptical galaxy NGC 720 between ∼0.002 R_vir and 0.2 R_vir (∼0.7 R₂₅₀₀). As found in past studies of the hot gas around galaxies (e.g., Buote et al. 2002; Humphrey et al. 2006, 2009a; Fukazawa et al. 2006; O'sullivan et al. 2007; Zhang et al. 2007; Nagino & Matsushita 2009), the data are inconsistent with a constant M/L ratio model. Given the high-quality constraints on the mass profile in NGC 720, this means that DM is required at 19.6σ.

With these data from deep Chandra and Suzaku observations, we have been able to obtain unprecedented constraints on the virial mass and DM halo concentration at this mass range (M_vir =(3.1 ± 0.4) × 10¹² M_☉). Such a low mass (only a few times more massive than the Milky Way) is consistent with our picture of this system as an isolated elliptical galaxy, rather than a galaxy group per se. These results are in agreement with, but very much tighter than, those found by Humphrey et al. (2006) for the same system. This improvement partially reflects the much better data used here (∼100 ks of Chandra data, plus ∼180 ks of Suzaku time, as compared to the ∼17 ks Chandra data set used by Humphrey et al.), but is also on account of various refinements made to our mass-modeling procedure in the interim. Of particular significance is the adoption of parameterized models for the entropy (rather than the temperature) distribution, allowing us to enforce a monotonically rising entropy profile (i.e., the Schwarzschild criterion for stability against convection). This efficiently eliminates unphysical portions of parameter space. Similarly, explicitly requiring the gas to remain approximately hydrostatic out to the virial radius also restricts parameter space to regions containing only physical models, whereas the fitting procedure of H06 was less restrictive.

In large part, the tight constraints on the virial mass and concentration are due to the ability of the data to constrain the flattening of the mass profile outside ∼20 kpc. This scale is much larger than the effective radius of the galaxy (3.1 kpc), and is thus unrelated to the stellar-DM conspiracy (Humphrey & Buote 2010), and instead must indicate the scale radius of the DM halo. As discussed in Gastaldello et al. (2007b), the ability to constrain the virial quantities is very sensitive to the accurate measurement of the scale radius. This flattening is not model dependent; when fitting the "cored logarithmic" DM model, we similarly found a flattening in the total gravitating mass profile at this scale.²¹

The small DM scale radius in NGC 720 indicates a high concentration; indeed, we find that it lies ∼3 times the intrinsic scatter above the predicted median relation found in DM only simulations for the "ΛCDM1" ("concordance") cosmology by Macciò et al. (2008). Buote et al. (2007) measured the c–M_vir relation empirically for a sample of early-type galaxies, groups, and clusters, finding higher concentrations in the galaxy regime (≲10¹³ M_☉) than Macciò et al. Still, we find that NGC 720 lies ∼2 times the intrinsic scatter above this relation.²² The extreme isolation of NGC 720 suggests an early epoch of formation (such that the brightest sub-halos have merged to form the central object), which should give rise to a higher DM halo concentration, and may provide a possible explanation for the measured c_vir (e.g., Zentner et al. 2005). The very young (∼3 Gyr) mean age for the central stellar population may simply reflect stars formed in a recent gas-rich merger rather than a realistic estimate of the age of the population as a whole (see discussion in Section 6.1).

One factor which could lead to an overestimate of c_vir is a failure to model correctly the stellar mass component (Mamon & Łokas 2005; Humphrey et al. 2006). If the young stellar population which dominates the light is, indeed, only a fraction of the total population in the galaxy, it is plausible that the assumption that mass follows light is not correct (although the modest color gradient of the galaxy is not consistent with a dramatic violation of this assumption; Peletier et al. 1990). However, excluding the central ∼4 kpc (∼1.3 R_e) does not significantly affect the concentration (Section 5.8). Alternatively, if adiabatic contraction operates, it should increase the cuspiness of the central DM halo profile, for a fixed stellar M/L ratio, which could also lead to c_vir being overestimated if it is not accounted for (e.g., Napolitano et al. 2010). Fitting an AC model (Section 5.1) also led to a slightly less concentrated halo, but the effect was not large enough to reconcile the best-fitting value with the Macciò et al. (2008) predictions. This comparative lack of sensitivity of the concentration to the details of the stellar modeling is not altogether surprising, since  R_e (3.1 kpc) is far smaller than the scale radius of the DM halo.

As expected for gas which is approximately hydrostatic, we obtain a good fit to the temperature and density profiles with a monotonically rising entropy profile. Although the entropy profile slope in the inner part of the system (β₁) is close to the canonical 1.1 (Table 2) predicted by gravitational structure formation simulations (considering only gravity without cooling or feedback; Voit et al. 2005), the profile flattens at smaller and larger radii. This is consistent with two of the three galaxies studied in Humphrey et al. (2009a), and similarly shaped profiles have also been observed in galaxy groups (e.g., Jetha et al. 2007; Finoguenov et al. 2007; Gastaldello et al. 2007a; Sun et al. 2009; H. M. L. G. Flohic et al. 2010, in preparation).

In Figure 12, we show the entropy profile as a function of R₅₀₀ and scaled by the "characteristic entropy," K₅₀₀ (e.g., Pratt et al. 2010). Clearly the profile is significantly flatter than the "baseline" model predicted by gravity-only simulations (Voit et al. 2005), and the normalization is significantly enhanced, reflecting the injection of entropy by non-gravitational processes, especially in the inner regions. In comparison with galaxy groups and clusters (e.g., Pratt et al. 2010; Sun et al. 2009), the entropy profile is more discrepant with the baseline model, reflecting the greater impact of non-gravitational heating on a lower mass halo. Nevertheless, the entropy injection still appears insufficient to evacuate a significant fraction of the baryons from the halo, as indicated by the approximate baryonic closure within R_vir. While this may reflect a finely balanced heating mechanism in the system, we note that most of the gas lies at large radii, where the entropy profiles are much less discrepant (for example, ∼65% of the gas lies outside R₅₀₀).

Zoom In Zoom Out Reset image size

If non-gravitational processes primarily redistribute the gas, rather than raising its temperature, we might expect the entropy profile to be simply related to the baseline model and the gas mass profile. Indeed Pratt et al. (2010) demonstrated with their cluster sample that the product S(R)(f_g(R)/f_b,U)^2/3, where S is the entropy profile, f_g(R) is the gas mass fraction profile (Figure 9), and f_b,U is the universal baryon fraction (0.17), is very close to the baseline model. It is interesting to investigate whether this also holds for lower-mass halos, where the impact of non-gravitational heating is systematically larger. We show the "f_g corrected" entropy profile for NGC 720 in Figure 12, finding that it also lies very close the gravity-only predictions, even in a ∼ Milky Way mass halo. We note that the corrected entropy profile agrees with this relation even when extrapolated outside the projected radii at which there are data, which provides support for our extrapolated entropy profile.

In order to interpret the impact of non-gravitational processes on the entropy distribution, it is helpful to consider the cooling time of the gas. In Figure 13, we show that, within ∼60 kpc, it is shorter than the Hubble time. Nevertheless, only ∼6% of the gas is currently within this region, implying that most of the gas in NGC 720 will not have had sufficient time to cool. In part this reflects the significant non-gravitational heating that has raised the entropy even in these outer regions significantly above the baseline model (Figure 12). The average stellar age in NGC 720 is young (∼3 Gyr), and if this reflects the time since NGC 720 formed, only the gas within ∼15 kpc would have had time to cool, within which there is presently only ∼0.5% (∼10⁹ M_☉) of the total gas mass.

Zoom In Zoom Out Reset image size

One of the intriguing results from our study is the detection, with both Chandra and Suzaku, of a strong negative abundance gradient in the hot gas of NGC 720 (Figure 3), making it arguably the lowest-mass system in which such a gradient has been definitively detected. Although the metallicity of the ISM in early-type galaxies with intermediate-to-high X-ray luminosities is now known to be ∼ solar, and comparable to the stellar metallicity (Humphrey & Buote 2006; Ji et al. 2009), simple chemical enrichment models typically predict substantially more metals in the central part of the halo, where enrichment from supernovae and stellar mass loss is concentrated (e.g., Pipino et al. 2005). Solutions to this problem may involve large-scale mixing due to AGN–ISM interaction or the assembly of structure through mergers (Mathews et al. 2004; Cora 2006; Cox et al. 2006). Clearly the spatial distribution of the metals in the hot gas is a useful constraint on any model predicting large-scale mixing.

While, to date, very little is known about the radial dependence of Z_Fe in the ISM of galaxy-scale halos, recent Chandra and XMM studies of relaxed galaxy groups have revealed similar, centrally peaked Z_Fe profiles (e.g., Humphrey & Buote 2006; Buote 2002; Buote et al. 2003; Rasmussen & Ponman 2007). Mathews et al. (2004) were able to explain the overall shape of these group-scale profiles by invoking large-scale flows driven by low-level AGN heating. It remains to be seen whether such a model could also explain the observed abundance profile of NGC 720 which, residing in a much shallower potential well, should be more susceptible to the effects of feedback.