The Extremely High Dark Matter Halo Concentration of the Relic Compact Elliptical Galaxy Mrk 1216

To make a quantitative analysis of the X-ray image morphology, we begin by computing the ellipticity (), position angle (PA), and centroid evaluated within elliptical apertures as a function of semimajor axis a. We apply an iterative scheme equivalent to diagonalizing the moment of inertia tensor of the image region (Carter & Metcalfe 1980; see Buote & Canizares 1994 for application to X-ray images of elliptical galaxies). Before applying this technique, we replaced detected point sources (Section 2) with local background using the CIAO dmfilth tool.

In Table 3 we list the ellipticity and PA as a function of a within 10 kpc (≈22''). Not listed in the table is the center position, which is quite steady; e.g., the center shifts by only 1.1 ± 0.4 pixels (05 ± 02) when the centroids of the a = 123, 2239 apertures are compared. Within the relatively large statistical errors, the PA within ≈10 kpc is consistent with the H-band value of 7015 reported by Y17, with some weak evidence that it increases near a = 10 kpc (also see below). The ellipticity, however, displays a clear radial variation. Within a ≈ 4'', the image is modestly flattened with  ≈ 0.20. The ellipticity then drops quickly for larger a to a small value ≈0.05 that is not inconsistent with  = 0. The X-ray morphology is thus broadly similar to that observed for the relaxed, fossil-like elliptical galaxy NGC 720 (e.g., Buote et al. 2002); i.e., within ≈1 R_e, the X-ray image is moderately flattened (though rounder than the stellar isophotes,  = 0.42—Y17) and consistent with being aligned with the stellar image before giving way to a much rounder X-ray image at larger radius. Hence, the moment analysis of the centroid, , and PA within a ≈ 10 kpc does not indicate the presence of irregular surface-brightness features.

Table 3.  Ellipticity Profile of the Central 10 kpc

a	a		PA
(arcsec)	(kpc)		(deg N-E)
1.23	0.55	0.20 ± 0.08	66 ± 45
3.69	1.66	0.22 ± 0.06	67 ± 30
6.15	2.77	0.06 ± 0.03	93 ± 26
8.61	3.88	0.09 ± 0.03	87 ± 12
11.32	5.09	0.07 ± 0.03	83 ± 20
14.27	6.42	0.05 ± 0.03	98 ± 28
17.96	8.08	0.04 ± 0.03	110 ± 28
22.39	10.07	0.05 ± 0.02	106 ± 22

Note. The ellipticity and position angle of the 0.5–2.0 keV surface brightness as a function of semimajor axis a obtained using a moment analysis (Section 3).

Download table as:  ASCIITypeset image

To search for more subtle features in the X-ray image, we construct a smooth two-dimensional model, subtract it from the image, and inspect the residual image using the sherpa fitting package⁵ within ciao. We initially defined a model consisting of an isothermal β model (Cavaliere & Fusco-Femiano 1978) for the hot gas and a constant background. Each model component was folded through the exposure map and fitted (using the C-statistic) to the full-resolution image within a radius of ≈100'' from the center of the galaxy (close to the edge of the S3). For our fiducial analysis here and throughout the paper, we defined the center of the X-ray image to be the centroid computed within a circular aperture of radius 3'' initially located at the emission peak. This gives (R.A., decl.) of (8:28:47.1410, −6:56:24.368), which is very consistent with the stellar center determined by Springob et al. (2014), although we note that there are differences at the 1'' level in the various position references collected by NED. (We examine the effect of choosing a slightly different center in Section 8.1.)

We found that the initial model fit produced significant residuals within the central region. We noted a substantial reduction in these residuals upon adding two Gaussian components with different widths; i.e., a crude multi-Gauss expansion (MGE). (Adding more Gaussian components produced comparatively minor changes.) The best-fitting parameters and 1σ errors are listed in Table 4. In Figure 2 we plot the best-fitting two-dimensional model (and data/model ratio) binned as a radial profile where each bin contains ∼200 counts. Note in particular the negligible residuals within the central region. The Gaussian components display and PA values similar to the moment analysis (Table 3); i.e., PA values broadly consistent with the H-band value with a moderately flattened  ≈ 0.30 that drops to a much smaller value ( ≈ 0.10), indicating nearly round isophotes at larger radius. (We emphasize that the individual components of our surface-brightness model—β model and two Gaussians—should not be thought of as distinct physically meaningful mass components. We describe the physical model(s) later; i.e., the fiducial HE model in Table 8.)

Zoom In Zoom Out Reset image size

Table 4.  Surface-brightness Model Parameters

	A	FWHM		PA	rc
Model	(cts s−1 arcmin−2)	(arcsec)		(deg N-E)	(arcsec)	β
Gauss 1	${6.01}_{-0.41}^{+0.44}$	${2.07}_{-0.12}^{+0.13}$	${0.29}_{-0.05}^{+0.05}$	${50}_{-6}^{+6}$	⋯	⋯
Gauss 2	${1.09}_{-0.12}^{+0.11}$	${6.75}_{-0.29}^{+0.34}$	${0.09}_{-0.04}^{+0.04}$	${92}_{-14}^{+14}$	⋯	⋯
Beta	${0.34}_{-0.10}^{+0.13}$	⋯	0	⋯	${6.2}_{-1.3}^{+2.0}$	${0.52}_{-0.02}^{+0.03}$
Const bkg	${0.0019}_{-0.0002}^{+0.0002}$	⋯	⋯	⋯	⋯	⋯

Note. Best-fitting parameters and 1σ errors of the two-dimensional, multicomponent surface-brightness (0.5–2.0 keV) model (Section 3).

Download table as:  ASCIITypeset image

When is allowed to vary for the β model, we obtain values  ≈ 0.13 and PA ≈ 135°, also fully consistent with the moment analysis, indicating that the PA begins to deviate significantly from the stellar value near a ≈ 30''. Proper assessment of potential systematic errors (e.g., from the treatment of embedded sources and the accuracy of the exposure map) on the values of and PA at these and larger radii is beyond the scope of our paper, and we defer such an analysis to a future investigation. Consequently, we fixed  = 0 for the β model for our present study, which we found had negligible impact on the fit residuals within the central ≈10 kpc region that is our focus here.

In Figure 3 we show the raw image in the central 20'' × 20'' region and the corresponding residual image constructed by subtracting the best-fitting two-dimensional model from the image. As is readily apparent, the residual image is noisy and lacks obvious bubbles or cavities or spiral features that are indicative of "sloshing." To guide the eye, we have overplotted smoothed contours (blue, linearly spaced) to trace subtle peaks and valleys. These regions are located within a radius of ∼5'' from the galaxy center without any obvious pattern in their locations.

Zoom In Zoom Out Reset image size

To study further the properties of these regions, we approximated the contour regions with simple boxes or circles, in some cases enlarging the regions to obtain more counts. We also added a few regions adjacent to the contours for comparison. Our region definitions are meant simply to provide a reasonable sampling of the contour regions and their surroundings and necessarily do not consider the location(s) of any extended radio jet emission, for which there is currently no evidence. Without having the radio jet emission as a guide, the statistical significance we quote below for the regions are overestimates because we do not account for the "look elsewhere effect." Therefore, while the absolute values of the quoted significances should be treated with caution, the relative significance values of the regions should still be useful for guiding future studies of the central image structure.

In all, we constructed 10 regions within a radius ∼8'' from the center. When compared to the model, 7 of 10 regions possess counts within 2σ of the model. We denote the three most significant deviations from the model by the dashed black regions in Figure 3 and list some of their properties in Table 5. Region 1, indicated by the rectangular region to the SW of the center, possesses by far the most significant (4.1σ) difference from the model, and its effect is even readily apparent in the raw image as a distortion in the third contour from the center. Regions 2 and 3 have significances just below 3σ, and their manifestations are not obvious in the raw image. Region 2 located to the SE is an excess of similar size (∼30% surface-brightness deviation from the model) to Region 1, but is less significant.

Table 5.  Residual Map Region Properties

	Center		Counts			Ratio	sign(Ratio)$\sqrt{| \mathrm{Ratio}| }$
Region	R.A.	Decl.	Image	Model	Residual	(%)	(%)
1	8:28:46.950	−6:56:26.092	239	175.5	+63.5	+36.2 ± 8.8	+16.7 ± 4.3
2	8:28:47.282	−6:56:27.420	143	108.1	+34.9	+32.2 ± 11.1	+15.0 ± 5.4
3	8:28:46.799	−6:56:24.885	36	57.1	−21.1	$-{36.9}_{-10.5}^{+12.3}$	$-{20.6}_{-5.1}^{+6.0}$
4	8:28:47.002	−6:56:22.486	147	169.3	−22.3	−13.2 ± 7.1	−6.8 ± 3.5

Note. Properties of notable regions of the 0.5–2.0 keV residual map (Section 3 and Figure 3). Regions 1–3 represent those with the largest residuals we studied; i.e., dashed regions in Figure 3. Region 4 displays the most interesting spectral deviations (Section 3.3); i.e., solid circle in Figure 3. The "Image" column gives the counts in the raw image, "Model" gives the counts predicted by the best-fitting model (Table 4 and Figure 2), and "Residual" gives the counts resulting when the model is subtracted from the image. The ratios refer to the data divided by the model expressed as percentages above (positive) or below (negative) the model. The final column takes the square-root of the data/model ratio and converts it into a percentage while preserving the sign. The result will indicate the hot gas density ratio if the differences between the plasma emissivities (especially the temperatures and iron abundances) of the data and model are negligible. The error bars on the ratios derive from Gaussian noise, except for Region 3, where we used the Poisson error bars tabulated by Gehrels (1986). The definitions of the regions are as follows: Region 1 is a rectangle with sides of lengths 35 and 17 rotated by 1218 NE; Region 2 is a rectangle also with sides of lengths 35 and 17, but rotated by only 09 NE; Region 3 is a circle with radius 11; and Region 4 is a circle with radius 13.

Download table as:  ASCIITypeset image

Region 3 has a deficit of ≈37% and a size that is well consistent with those seen in well-studied cavity systems in clusters (e.g., McNamara & Nulsen 2007). The fact that Region 1, an excess, is adjacent to Region 3 is intriguing. The configuration might be a rim bordering a cavity, although the relative placement (cavity at larger radius than the rim) would not obviously favor this interpretation. Below in Section 3.3 we examine the spectra of these regions and find gas parameters consistent with annular averages within the sizable error bars due to the relatively few counts in these regions. The most significant spectral difference we found is located in the region denoted by the solid black circle in Figure 3 and Region 4 in Table 5. This region, however, corresponds only to a ≈13% deficit (1.8σ), and we discuss it further in Section 3.3.

We conclude that at present, the Chandra X-ray image of Mrk 1216 does not reveal obvious features of AGN feedback in the form of bubbles, cavities, weak shocks, or other irregularities in the central surface brightness. Nevertheless, in this section we have identified regions of the most prominent surface-brightness deviations from a smooth two-dimensional model as leading candidates for such feedback signatures to be studied with future high-sensitivity X-ray and radio observations.

We model the interstellar plasma ("hot gas") using the vapec optically thin coronal plasma model with version 3.0.9 of the atomic database AtomDB⁶ and the solar abundance standard of Asplund et al. (2006). In our implementation of the vapec model, for elements heavier than He (which is kept fixed at solar abundance), we fit the ratios of the metal abundances with respect to iron; e.g., for Si, we fit Z_Si/Z_Fe rather than Z_Si itself. The free parameters we considered for the hot gas component in each spectrum are k_BT, Z_Fe, Z_Mg/Z_Fe, Z_Si/Z_Fe, and the normalization. All other elements heavier than He are fixed in their solar ratios with iron. We do not fit the other elements individually because they are too strongly blended with iron (e.g., Ne), are affected by background (e.g., S), or are simply too poorly constrained.

For several reasons, throughout most of this paper, we fit the plasma models directly to the observed spectra without performing any onion-peeling–type deprojection. While spectral deprojection to some extent can help to mitigate possible biases associated with fitting a single-temperature model to a multitemperature spectrum, this advantage is outweighed by some key disadvantages. Deprojection generally, and onion-peeling in particular, amplifies noise especially in the outer regions where the background dominates. Standard deprojection procedures also do not easily self-consistently account for the gas emission outside the bounding annulus, which can be a sizable source of systematic error (e.g., Nulsen & Bohringer 1995; McLaughlin 1999; Buote 2000a). They also typically assume that the gas properties are constant within what are often wide spherical shells (especially for systems like Mrk 1216), introducing additional systematic error for the radially varying gas properties. (The assumption of constant properties per circular annulus on the sky also applies to our default approach, but the errors associated with this assumption do not propagate between annuli in the same way as the spectral deprojection in which the model spectrum in any given shell depends on all of those exterior to it.) Consequently, we relegate deprojection analysis using the projct mixing model in xspec to a systematic error check (Section 4.2.2).

Although the emission from unresolved LMXBs and other stellar sources is a small fraction of the X-ray emission of Mrk 1216, we still included a 7.3 keV thermal bremsstrahlung component (e.g., Matsumoto et al. 1997; Irwin et al. 2003) to account for this emission. We restricted the normalization of this component to lie within a factor of 2 of the L_x–L_K scaling relation for discrete sources of Humphrey & Buote (2008) using the K-band luminosity (L_K = 1.7 × 10¹¹ L_⊙) from the Two Micron All-sky survey (2MASS) as listed in the Extended Source Catalog (Jarrett et al. 2000). Using the global L_x from the scaling relation, we assigned the expected range of L_x for each annulus according to the fraction of the total 2MASS K-band emission falling into that annulus. Hence, for each annulus, the flux of unresolved discrete sources (with range restricted as noted) is a free parameter.

As described in Section 3.1.2 of B17, we model the cosmic X-ray background (CXB) emission with multiple thermal plasma components for the "soft" CXB and a single power law for the "hard" CXB. By default, we fixed the soft CXB normalizations to those obtained from fitting ROSAT data using the HEASARC X-ray Background Tool.⁷ We examine the systematic error associated with this choice by allowing the normalizations of the soft CXB components to vary within a factor of 2 of the ROSAT values (Section 8). Hence, in our default model the normalization of the power law of the hard CXB contribution for all annuli is the only free parameter of the CXB.

For the particle background, we adopted a multicomponent model consisting of a power law with two break radii along with three Gaussians. Unlike the other models described above, we do not fold the particle background model through the ARF; i.e., it is "unvignetted." However, because the particle backgrounds of the BI and FI chips are not identical, we fitted separate versions of the model to the data on the BI (S3) and FI chips.

The Cycle 16 and Cycle 19 data are fitted separately. For each data set we fitted all annuli simultaneously, including large apertures at the largest radii (not listed in Table 6) dominated by background.

4.2.1. Analysis of Projected Spectra

The spectral model described above describes the Cycle 19 data well, with a minimum C-statistic value of 2023.5 (in 2002 pha bins) with 1932 degrees of freedom (dof). The success of the model is shown in Figure 4, where we plot the spectra and best-fitting models for five representative annuli. The most noticeable spectral features are the broad bump near 1 keV dominated by a forest of emission lines from the Fe L shell and the strong Si Kα line near 1.85 keV. Less noticeable, but still prominent in the inner annuli, is the Mg Kα line complex near 1.4 keV. (Note that for the Cycle 16 data we achieve a fit consistent with that obtained in Paper I.)

Zoom In Zoom Out Reset image size

The innermost and outermost annuli deserve special mention. Annulus 1 displays the most significant residuals from the best-fitting model, resulting in a fit there that is of formally marginal quality. In Appendix C we study the central spectrum in detail and examine several possibilities to improve the fit, although at present we cannot confidently recommend a specific modification to the fiducial model. Because the background dominates in the outermost annulus (Annulus 10), the properties of the gas component there cannot be robustly determined. We found it necessary in the spectral fitting to restrict the gas parameter ranges there more strongly: i.e., k_BT = 0.4–0.9 keV. Guided by the average radial profiles of Z_Fe for groups and clusters obtained by Mernier et al. (2017), by default we fixed Z_Fe = 0.36 Z_⊙ in Annulus 10. Mernier et al. (2017) quote a scatter of ∼±0.09 for Z_Fe over the radial range corresponding to Annulus 10, and we therefore use the range ${Z}_{\mathrm{Fe}}=0.27-0.45$ Z_⊙ as a systematic error in our HE models (Section 8.3).

We list the gas parameters measured for each annulus in Table 6 for the data sets. (We express the emission measure of the gas as a surface-brightness unit—see the notes to the table.) In Figure 5 we plot the radial profiles of the surface brightness (Σ_x) and temperature (k_BT). As expected, Σ_x and k_BT are consistent in their overlap region (as is Z_Fe). The lack of a great temperature jump in Annulus 1 of the Cycle 19 observation has implications for the mass of the SMBH (Section 7.3). The Cycle 19 data confirm and strengthen the similarity of the temperature profile of Mrk 1216 to that of the fossil group NGC 6482 (B17).

Zoom In Zoom Out Reset image size

The profile of Z_Fe decreases with radius but is nearly constant over large stretches; i.e., Z_Fe ≈ 1 Z_⊙ for R ≲ 2 kpc and Z_Fe ≈ 0.7 Z_⊙ for R ≈ 4–60 kpc. This negative gradient in Z_Fe is very similar to gradients seen in several X-ray bright, massive elliptical galaxies and small groups like NGC 6482 (B17), NGC 5044 (Buote et al. 2003b), and others (e.g., Buote 2000a; Humphrey & Buote 2006; Rasmussen & Ponman 2009; Mernier et al. 2017).

Although the Mg and Si abundances are not as well constrained as Fe, we find that the Cycle 19 observation does place interesting constraints on the radial variation of either Z_Mg/Z_Fe or Z_Si/Z_Fe. The Z_Mg/Z_Fe ratio appears to peak in the central R ≈ 2 kpc with a value that is at least solar and is consistent with a constant ratio of ≈0.7 solar at larger radius. The Z_Si/Z_Fe profile is broadly similar to Z_Mg/Z_Fe out to R ≈ 10 kpc, after which it increases significantly to Z_Si/Z_Fe ≈ 1.3 solar. Because at large radius the Si abundance measurement becomes especially more challenging due to the increasing background level (both the continuum and the presence of an instrumental line), we regard our measurement there as provisional. Nevertheless, the Z_Mg/Z_Fe and Z_Si/Z_Fe profiles are broadly similar to the mean profiles obtained from XMM-Newton for groups by Mernier et al. (2017).

For comparison, if we do not allow for a radial variation in either Mg or Si, we obtain Z_Mg/Z_Fe = 0.83 ± 0.06 solar and Z_Si/Z_Fe = 0.97 ± 0.07 solar for the Cycle 19 data (which also gives Z_Fe ≈ 0.80 Z_⊙ in Annulus 1). We use these results to perform a systematic error check on our fiducial models in Section 8. Note also that these constant values of Z_Mg/Z_Fe and Z_Si/Z_Fe for the Cycle 19 data are consistent within the ≈1.5σ errors with the values obtained for the Cycle 16 data (Table 6), for which interesting constraints on the radial variation were not obtained.

Finally, the Cycle 19 results we have described in this section for the combined data are fully consistent with those obtained when performing a joint fit of the individual exposures (see Appendix A and Table 14).

4.2.2. Spectral Deprojection in Spherical Shells

4.2.3. Central Region: Quadrants

Because the main objective of our paper is to infer the gravitating mass distribution using an HE analysis, ideally we would like kinematic information for the hot plasma to inform and correct our analysis, especially for the central region where AGN feedback is expected to periodically inject energy into the hot gas. High spectral resolution observations of the Perseus cluster with Hitomi (Hitomi Collaboration et al. 2016) and in massive elliptical galaxies/small groups with the XMM-Newton RGS (e.g., Ogorzalek et al. 2017) indicate low amounts of turbulent pressure at the centers of these systems.

Another manifestation of non-hydrostatic behavior is via spatial fluctuations in the gas properties, in particular, through the azimuthal scatter of properties within subregions of circular annuli (e.g., Vazza et al. 2011). The high spatial resolution combined with the moderate energy resolution of Chandra ACIS is well suited for studying such azimuthal fluctuations in the central, high S/N regions of Mrk 1216. We were able to obtain useful constraints on the gas properties when dividing up Annuli 2, 3, and 4 into four quadrants, where for each annulus we fixed the metal abundances and background levels to the best-fitting results obtained for the whole annulus in Section 4.2.1.

Because HE is a balance between pressure and gravity at any point, we use the inferred projected gas properties to construct three-dimensional proxies for the gas density, entropy, and pressure as follows. The normalization of the vapec model in xspec,

$$\begin{eqnarray}&&{\rm{}}\mathrm{norm}\equiv {10}^{-14}\int {n}_{e}{n}_{{\rm{H}}}{dV}/\left(4\pi \left[{D}_{A}{\left(1+z\right)}^{2}\right]\right),\end{eqnarray} \tag{ 1 }$$

is proportional to the emission measure, $\int {n}_{e}{n}_{{\rm{H}}}{dV}$. We define a pseudo-electron number density (${\tilde{n}}_{e}$) by dividing this emission measure by ${V}_{{ij}}^{\mathrm{int}}={({r}_{j}^{2}-{r}_{i}^{2})}^{3/2}$, taking the square root, and converting n_H into n_e. Here r_i and r_j are the inner and outer radii of the annulus on the sky representing the volume of the spherical shell intersected by the cylindrical annulus along the line of sight of the same radii (e.g., Kriss et al. 1983). We then use this pseudo-electron number density and the projected k_BT to define correspondingly a pseudo-entropy and pseudo-pressure.

In Table 7 we list the results for these gas properties obtained for each quadrant for Annuli 2–4; see the caption to the table for the definitions of the quadrants. (We obtain consistent results for the scatter if we rotate the quadrants by 45°.) To quantify the scatter between the quadrants of each sector, we follow Vazza et al. (2011) and define the quantity,

$$\begin{eqnarray}&&\mathrm{Scatter}\equiv \sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{y}_{i}-\bar{y}}{\bar{y}}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where y_i is the quantity in quadrant i, $\bar{y}$ is the average value of the quantity over the whole annulus, and N = 4. The errors quoted for the scatter assume normal error propagation.

Table 7.  Hot Gas Properties in Quadrants

		kBT	norm	${\tilde{n}}_{e}$	$\tilde{S}$	$\tilde{P}$
Annulus	Quadrant	(keV)	(10−5 cm−5)	(cm−3)	(keV cm2)	(10−10 erg cm−3)
2	1	0.951 ± 0.038	1.07 ± 0.07	0.156 ± 0.005	3.29 ± 0.15	2.375 ± 0.120
2	2	0.842 ± 0.040	1.26 ± 0.07	0.169 ± 0.005	2.75 ± 0.14	2.279 ± 0.127
2	3	0.898 ± 0.045	1.08 ± 0.07	0.156 ± 0.005	3.10 ± 0.17	2.250 ± 0.132
2	4	0.892 ± 0.037	1.28 ± 0.07	0.171 ± 0.005	2.90 ± 0.13	2.440 ± 0.122
2	Scatter	0.043 ± 0.022	0.085 ± 0.030	0.043 ± 0.015	0.067 ± 0.024	0.032 ± 0.027
3	1	0.839 ± 0.033	1.04 ± 0.07	0.091 ± 0.003	4.13 ± 0.19	1.230 ± 0.063
3	2	0.828 ± 0.038	1.10 ± 0.07	0.094 ± 0.003	4.00 ± 0.20	1.249 ± 0.069
3	3	0.812 ± 0.039	0.97 ± 0.06	0.088 ± 0.003	4.10 ± 0.22	1.147 ± 0.067
3	4	0.913 ± 0.048	0.83 ± 0.06	0.081 ± 0.003	4.86 ± 0.28	1.191 ± 0.076
3	Scatter	0.045 ± 0.026	0.105 ± 0.032	0.054 ± 0.017	0.080 ± 0.030	0.033 ± 0.028
4	1	0.782 ± 0.032	1.58 ± 0.09	0.044 ± 0.001	6.27 ± 0.28	0.551 ± 0.027
4	2	0.827 ± 0.035	1.37 ± 0.08	0.041 ± 0.001	6.96 ± 0.32	0.543 ± 0.028
4	3	0.807 ± 0.036	1.19 ± 0.08	0.038 ± 0.001	7.10 ± 0.35	0.494 ± 0.027
4	4	0.808 ± 0.036	1.24 ± 0.08	0.039 ± 0.001	7.01 ± 0.34	0.506 ± 0.028
4	Scatter	0.020 ± 0.020	0.110 ± 0.032	0.054 ± 0.015	0.048 ± 0.022	0.046 ± 0.026

Note. See Section 4.2.3 for definitions of the xspec norm parameter, pseudo-electron number density ${\tilde{n}}_{e}$, pseudo-entropy $\tilde{S}$ = k_BT ${\tilde{n}}_{e}^{-2/3}$, and pseudo-pressure $\tilde{P}={\tilde{n}}_{e}$k_BT. The annuli numbers correspond to those in Table 6. The quadrants are defined with respect to the H-band stellar position angle, θ_⋆ = 7015 (NE) from Y17: (θ_⋆, θ_⋆ + 90°) for quadrant 1, (θ_⋆ + 90°, θ_⋆ + 180°) for quadrant 2, (θ_⋆ + 180°, θ_⋆ + 270°) for quadrant 3, (θ_⋆ + 270°, θ_⋆ + 360°) for quadrant 4. The "Scatter" is the mean fractional scatter of the quadrants for a particular annulus defined by Equation (2).

Download table as:  ASCIITypeset image

The largest scatter (≈0.10) occurs for norm, while the other parameters have smaller scatters (≲0.05) in most cases. The HE equation may be expressed in terms of the entropy and pressure having a dependence S^3/5 and P^2/5, respectively (e.g., Humphrey et al. 2008). The scatter in the pseudo-entropy and pseudo-pressure listed in Table 7 suggest ≲5% azimuthal fluctuations in the HE equation in the central region of Mrk 1216. Errors of this magnitude are smaller than the statistical errors on the inferred mass properties (e.g., Table 10).

4.2.4. Central Region: Image Residuals

The primary method we adopt to fit the HE models to the Chandra data employs a Bayesian nested sampling procedure based on the MultiNest code v2.18 (Feroz et al. 2009; see B17 for details). In Table 8 we list the priors adopted for each free parameter. We use flat priors in most cases and flat priors on the decimal logarithm for a few parameters with ranges spanning multiple orders of magnitude. For the black hole mass, our fiducial prior is based on the M_BH–σ relation of van den Bosch (2016); i.e., a median value M_BH = 2.1 × 10⁹ M_⊙ using the stellar velocity dispersion for Mrk 1216 listed in Table 1. Because van den Bosch (2016) quote an intrinsic scatter of 0.49 in log₁₀ M_BH, we adopt a lognormal prior with mean ${\mathrm{log}}_{10}$ M_BH = 9.32 and standard deviation equal to the intrinsic scatter.

For the Bayesian analysis, we quote two "best" values for each free parameter: (1) the mean parameter value of the posterior, which we call the "Best Fit", and (2) the parameter value that maximizes the likelihood, which we call "Max Like." Unless otherwise stated, all errors quoted for the parameters are the standard deviation (1σ) of the posterior.

We also perform frequentist χ² fits of the HE models for comparison to the results obtained with the Bayesian fits. Despite its shortcomings (e.g., Andrae et al. 2010), we also prefer the frequentist χ² approach for model selection, which, unlike Bayes factors, does not depend on the priors. For the frequentist fits, we use the minuit fitting package that is part of the root v6.10 software suite.⁸

In Paper I we found that the Cycle 16 data were described well by an entropy profile consisting of only a constant and a power law with no breaks in radius. Fitting such a model jointly to the Cycle 16 and Cycle 19 data still results in a formally acceptable fit: χ² = 28.3 for 26 dof ("0 Brk Entropy" in Table 9). However, when we add a break radius, the fit improves significantly (χ² = 17.5 for 24 dof, "1 Brk Entropy" in Table 9). The one-break entropy profile has a poorly constrained break radius r_b,1 = 7 ± 5 kpc, and its best-constrained parameter is the power-law exponent exterior to the break radius (α₂ = 0.93 ± 0.05). Adding more breaks reduces χ², but does not lead to a statistically significant improvement in the fit.

Nevertheless, we add more break radii to provide greater flexibility in the entropy model for the following reason. When we compare models with different assumptions in the mass components (e.g., fixed overmassive SMBH, Einasto DM halo, and AC), we require that potential differences in the fits are determined by differences in the form of the mass components, not the precise form of the assumed entropy profile. HE and convective stability only require that the entropy profile increase monotonically with radius. Consequently, we added break radii until the value of χ² changed by less than 1.

We followed this procedure and arrived at an entropy profile with four breaks having break radii ≈0.7, 16, 22, and 58 kpc. (We obtained this final result after initially using larger prior ranges than indicated in Table 8.) The inner break adds flexibility when different SMBH priors are tested, while the others allow for flexibility primarily for the DM halo models. The two break radii in the middle are rather close together, indicating a fairly abrupt jump in entropy near 20 kpc. The slope parameter at large radius is not well constrained and is consistent with the r^1.1 profile from gravity-only simulations (Voit et al. 2005). Although this four-break entropy profile is the fiducial model we employ (Section 5, Table 8), in Section 8 we compare it to results obtained with the one-break entropy profile as a systematics error check. (The best-fitting one-break model is plotted in Figure 6 as the red dotted line.)

In Figure 6 we plot the entropy profile of the Bayesian fiducial model along with the ∼r^1.1 theoretical entropy profile produced by gravity-only cosmological simulations (Voit et al. 2005). The entropy profile of Mrk 1216 is remarkably similar to the profiles we have obtained previously for other massive, fossil-like elliptical galaxies, NGC 720 (Humphrey et al. 2011), NGC 1521 (Humphrey et al. 2012), and NGC 6482 (B17). That is, within a radius ∼r₂₅₀₀, the entropy profile greatly exceeds the theoretical gravity-only profile. In addition, when the observed entropy profile is rescaled by $\sim {f}_{\mathrm{gas}}^{2/3}$ (Pratt et al. 2010), the result broadly matches the gravity-only profile. The success of this rescaling indicates that for r ≲ r₂₅₀₀, the feedback energy responsible for raising the entropy in that region is consistent with having spatially rearranged the gas rather than raising its temperature, as also is observed for galaxy clusters (Pratt et al. 2010).

We see, however, that as the entropy approaches r₂₅₀₀, this rescaling becomes increasingly less successful. Very similar behavior is observed for NGC 1521 (Humphrey et al. 2012) and NGC 6482 (B17). Perhaps for Mrk 1216 and these galaxies at larger radius, the feedback energy heated the gas by raising its temperature, indicating that a different feedback mechanism prevails for r ≳ r₂₅₀₀. It is also noteworthy that interior to ≈0.6r₂₅₀₀ (77 kpc), the cooling time is shorter than the age of the universe (see Figure 9), and thus the details of gas cooling may be responsible for the apparent transition in the rescaling behavior as the radius approaches r₂₅₀₀.

The pressure profile of the fiducial Bayesian HE model (not shown) is extremely similar to that of NGC 6482 (see Section 6.3 and Figure 5 of B17). Between radii ≈0.1–0.7 r₅₀₀, the pressure profile broadly matches the "universal" pressure profile that was determined for galaxy clusters by Arnaud et al. (2010), while at smaller and larger radii, the observed profile is clearly a different shape. While the results for r ≳ 0.7r₅₀₀ are provisional owing to the limited data range, the different shape at smaller radius is robust, implying a breakdown in the mass scaling of the universal pressure profile, and presumably, the greater importance of nongravitational energy in shaping the thermodynamical properties of the gas in galaxy/group-scale halos.

In Paper I we analyzed the Cycle 16 observation and obtained a Bayesian constraint on the black hole mass, M_BH = (5 ± 4) × 10⁹ M_⊙, for a flat prior very consistent with the stellar dynamical measurement, M_BH = (4.9 ± 1.7) ×10⁹ M_⊙, of Walsh et al. (2017). When we instead use a flat prior on log₁₀ M_BH ("flat logspace prior"), we obtained a smaller value, M_BH = (1.4 ± 1.7) × 10⁹ M_⊙, which is still consistent with the stellar dynamical measurement within the large errors. Our best-fitting model in Paper I for the flat prior on M_BH predicted a projected emission-weighted temperature, k_BT ∼ 1.18 keV, within the central R = 1''; i.e., the region corresponding to Annulus 1 of the new Cycle 19 observation.

As is readily apparent from the temperature profile displayed in Figure 5 and the k_BT values listed in Table 6, we measure a much lower temperature in Annulus 1 with the Cycle 19 observation: k_BT = 0.969 ± 0.027 keV, implying a smaller value of M_BH than indicated in Paper I. (Note that only for Annulus 1 was the predicted k_BT from Paper I inaccurate; e.g., for Annulus 2, the predicted value was 0.910 keV compared to the value measured of k_BT = 0.905 ± 0.018 keV.) For a flat prior, our joint fit of the Cycle 16 and Cycle 19 data gives M_BH = (1.0 ± 0.6) × 10⁹ M_⊙, with a 99.9% upper limit, M_BH ≤ 3.9 × 10⁹ M_⊙. The flat logspace prior fits to even lower values approaching the lower limit of the adopted prior range, M_BH = (0.3 ± 0.2) × 10⁹ M_⊙ with a 99.9% upper limit, M_BH ≤ 1.1 × 10⁹ M_⊙.

As mentioned in the notes to Table 9, the frequentist fit also gives M_BH consistent with the adopted lower fit boundary. However, when M_BH is fixed to the stellar dynamical value, χ² only increases by 2.3 ("fixed overmassive BH" in Table 9). According to the F-Test, the fiducial model is ≈92% more probable; i.e., the preference for a small value of M_BH is lower than a 2σ effect for the frequentist fit.

Fixing M_BH to the stellar dynamical value leads to a smaller inferred value for M_⋆/L_H ("fixed overmassive BH" in Table 10). The Bayesian fit gives M_⋆/L_H = 0.73 ± 0.13 solar with a 99.9% upper limit, M_⋆/L_H ≤ 0.99 solar, significantly lower than M_⋆/L_H = 1.2 solar expected for a Kroupa IMF (Section 7.4). The frequentist fit gives M_⋆/L_H = 0.96 ± 0.13 solar, which is ∼2σ discrepant.

Table 10.  Stellar and Total Mass

	M⋆/LH	MBH	c2500	M2500	c500	M500	c200	M200
	(${M}_{\odot }\,{L}_{\odot }^{-1}$)	(109 M⊙)		(1012 M⊙)		(1012 M⊙)		(1012 M⊙)
Best Fit	1.19 ± 0.11	0.8 ± 0.4	11.1 ± 1.6	3.2 ± 0.3	21.1 ± 3.0	4.5 ± 0.5	30.4 ± 4.3	5.3 ± 0.6
(Max Like)	(1.22)	(0.6)	(10.9)	(3.2)	(20.7)	(4.4)	(29.9)	(5.3)
1 Brk Entropy	−0.04	−0.1	0.8	−0.1	1.5	−0.2	2.1	−0.2
BH Flat Prior	−0.03	0.2	0.4	−0.1	0.7	−0.1	1.0	−0.1
BH Flat Logspace Prior	0.03	−0.5	−0.1	0.0	−0.3	0.0	−0.4	0.0
Fixed Overmassive BH	−0.46	4.1	3.7	−0.4	6.8	−0.6	9.8	−0.8
Fixed MBH–σ BH	−0.13	1.3	1.0	−0.1	1.8	−0.2	2.6	−0.3
Sérsic Stars 2MASS	0.04	−0.0	−1.5	0.3	−2.8	0.5	−4.0	0.6
Einasto	−0.07	0.0	−1.0	0.4	−2.1	0.5	−3.2	0.5
Strong AC	−0.40	−0.0	−3.3	0.6	−6.1	1.1	−8.8	1.4
Weak AC	−0.06	−0.0	−2.5	0.3	−4.6	0.6	−6.5	0.8
Weak AC Einasto	−0.13	0.0	−3.6	0.9	−6.6	1.6	−9.6	1.9
Joint Fit of Cycle 19 Obs.	−0.04	−0.0	1.6	−0.2	2.9	−0.4	4.2	−0.5
Constant Zα/ZFe	0.13	−0.0	−0.6	0.1	−1.1	0.2	−1.6	0.2
Annulus 10 ΔZFe	${}_{-0.00}^{+0.02}$	−0.0	${}_{-0.4}^{+0.4}$	${}_{-0.1}^{+0.2}$	${}_{-0.7}^{+0.7}$	${}_{-0.2}^{+0.2}$	${}_{-1.1}^{+1.0}$	${}_{-0.2}^{+0.3}$
Deproj	−0.16	0.0	1.3	−0.2	2.5	−0.3	3.9	−0.2
Distance	−0.16	0.0	−1.1	0.2	−2.1	0.4	−2.9	0.6

Note. Best-fit values and 1σ error estimates for the free parameters of the mass components of the fiducial Bayesian hydrostatic equilibrium model; i.e., stellar mass-to-light ratio (M_⋆/L_H), black hole mass (M_BH), concentration, and enclosed total mass (BH+stars+gas+DM). We show the concentration and mass results obtained within radii r_Δ for overdensities Δ = 200, 500, and 2500. In addition, we provide a budget of systematic errors where only the most significant/interesting systematics are shown (Section 8). For each column we quote values with the same precision. In several cases, an error has a value smaller than the quoted precision, and thus it is listed as a zero; e.g., "0.0" or "−0.0", where the sign indicates the direction of the shift. Briefly, the various systematic errors are as follows: ("1 Brk Entropy"): Entropy profile has only one break radius. ("BH Flat Prior"): Flat prior on M_BH ranging from 10^8–10 M_⊙. ("BH Flat Logspace Prior"): Flat prior on log₁₀M_BH ranging from 8–10 ("Fixed Overmassive BH"): M_BH fixed to 4.9 × 10⁹ M_⊙ (Walsh et al. 2017). ("Fixed M_BH–σ BH"): M_BH fixed to 2.1 × 10⁹ M_⊙, the median of the M–σ relation (van den Bosch 2016). ("Sérsic Stars 2MASS"): The stellar mass is represented by a Sérsic model with R_e from 2MASS (see Section 7.4). ("Einasto"): Einasto DM profile. ("Strong AC"): Adiabatically contracted DM profile following Blumenthal et al. (1986). ("Weak AC"): Weak adiabatically contracted DM profile following Dutton et al. (2015); i.e., their "Forced Quenched" model implemented as described in B17. ("Weak AC Einasto"): Same as "Weak AC", but here the Einasto DM profile is used. ("Joint Fit of Cycle 19 Obs."): Spectral fitting of Chandra Cycle 19 observation performed simultaneously on each observation segment (Section 2). ("Constant Z_α/Z_Fe"): Z_Mg/Z_Fe and Z_Si/Z_Fe not allowed to vary with radius in the spectral analysis (Section 4). ("Annulus 10 ΔZ_Fe"): Spectral fitting of Chandra Cycle 19 observation performed for different choices of Z_Fe in Annulus 10 (Section 8.3). ("Deproj"): Spectral analysis performed with deprojection assuming constant spectral properties within each annulus (Sections 4.2.2 and 8.2). ("Distance"): Use luminosity distance of 113 Mpc (NED; Springob et al. 2014).

Download table as:  ASCIITypeset image

Therefore, the Bayesian fits clearly disfavor the overmassive SMBH, as indicated by the values inferred for M_BH when it is fitted freely and for M_⋆/L_H when M_BH is fixed to the stellar dynamical value. (The significantly larger c₂₀₀ value provides more evidence against the overmassive SMBH model—see Sections 7.6 and 9.1.) The frequentist analysis also disfavors the overmassive SMBH, but at a much lower significance level (∼2σ).

Given the dependence of the constraints on the priors, by default, we adopt a prior on M_BH consistent with the M_BH–σ relation of van den Bosch (2016)—see Section 6.1. In Section 9.3 we discuss the implications of the smaller values inferred for M_BH compared to the stellar dynamics measurement.

The Chandra data clearly require the stellar mass component. If the stellar mass is omitted, the frequentist fit gives χ² = 40.5 for 19 dof (Table 9), which is an increase of 28.6 over the fiducial model with one extra dof. The strong need for the stellar component translates to a good constraint on stellar mass-to-light ratio; i.e., M_⋆/L_H = 1.19 ± 0.11 for the fiducial Bayesian model (Table 10). This value agrees very well with that expected from single-burst stellar population synthesis (SPS) models with a Kroupa IMF (M_⋆/L_H = 1.2), but is significantly lower than predicted for a Salpeter IMF (M_⋆/L_H = 1.7 solar); see Walsh et al. (2017) for more information about these SPS estimates.

Our fiducial measurement also agrees very well with the value 1.3 ± 0.3 solar (statistical error) obtained by the stellar dynamical study of Walsh et al. (2017), with a factor of ≈2.5 higher precision. Y17's stellar dynamical measurement gives a value about 50% larger, M_⋆/L_H = ${1.9}_{-0.4}^{+0.5}$ solar at $\approx 2\sigma$ significance, which agrees better with the SPS value for a Salpeter IMF. In Section 9.2 we discuss these stellar dynamical measurements and their sensitivity to the assumed DM halo concentration.

Most of the systematic errors we considered do not change M_⋆/L_H by more than the 1σ statistical error (Table 10). Of key importance is that if we use a different model to represent the stellar light profile (Section 8.6), we obtain M_⋆/L_H = 1.23 solar, very consistent with the MGE result. Much smaller values of M_⋆/L_H are given by a few models, "fixed overmassive BH" (Section 7.3), "strong AC" (Section 7.5), and to a lesser extent, "deproj" and "distance", which we believe is notable evidence disfavoring these models. The "constant Z_α/Z_Fe" test give the largest positive shift in M_⋆/L_H.

The very good agreement we find with the value of M_⋆/L_H we measure and SPS models with a Kroupa IMF is also very consistent with what we find for other fossil-like massive elliptical galaxies from X-ray HE analysis. This evidence for a Kroupa IMF is in stark contrast to the support that is typically found for a Salpeter IMF from stellar dynamics and lensing studies of massive elliptical galaxies (see discussion in Section 8.3 of B17).

Of the SMBH, stars, and DM, by far the most important component needed to describe the X-ray emission is the DM (Table 9). If the DM halo is omitted from the fiducial model, the frequentist fit gives χ² = 395.8 for 20 dof; i.e., the F-test gives a tiny probability (2 × 10⁻¹⁴) that the "no DM halo" model is preferred over the fiducial model with DM. While the DM component is clearly required, the Chandra data do not distinguish between the different DM models we investigated. The Einasto and CORELOG models give minimum χ² values that are <1 from that obtained for the fiducial model and with the same number of dof (Table 9). The weak AC model applied to the NFW and Einasto models changes the fit quality very little. The strong AC NFW model has a χ² value that is slightly larger (i.e., by 2.9), but is still a good fit.

As noted above in Section 7.4, the strong AC model gives a much smaller M_⋆/L_H, which we believe disfavors that model. The CORELOG model (not shown in Table 10) is peculiar in that the "Best Fit" and "Max Like" values for M_⋆/L_H differ significantly; i.e., M_⋆/L_H ∼ 0.7 solar for "Best Fit" and M_⋆/L_H ∼ 1.2 solar for "Max Like" with 1σ error ±0.2 solar, where the larger "Max Like" values are favored to match the SPS models (Section 7.4).

These results are very consistent with those obtained for NGC 6482 (B17).

In Paper I we obtained tentative evidence for an "overconcentrated" NFW DM halo compared to the general halo population. Our Bayesian fiducial model fitted to the Cycle 16 Chandra data yielded a Best Fit virial mass, M₂₀₀ =(9.6 ±3.7) × 10¹² M_⊙ and concentration, c₂₀₀ = 17.5 ± 6.7 considerably larger than the median ΛCDM value of ≈7 (e.g., Dutton & Macciò 2014). Intriguingly, we found that the Max Like best-fitting value suggested an even more extreme outlier from the ΛCDM relation: c₂₀₀ = 25.9 and M₂₀₀ = 5.1 ×10¹² M_⊙.

Adding the much deeper Cycle 19 Chandra observations to our analysis confirms the higher concentration solution from Paper I. In Table 10 we list the concentration and mass for overdensities Δ = 200, 500, and 2500. Now the Best Fit and Max Like values agree extremely well. The values we measure, c₂₀₀ = 30.4 ± 4.3 and M₂₀₀ = (5.3 ± 0.6) × 10¹² M_⊙, indicate an extreme outlier to the median ΛCDM relation (see Section 9.1).

The AC models give statistically significant lower values of c₂₀₀, with Strong AC giving a larger reduction (≈9) than Weak AC (≈7). As noted above in Section 7.4, the Strong AC model gives an uncomfortably low M_⋆/L_H, whereas Weak AC gives M_⋆/L_H fully compatible with the fiducial model. The Fixed Overmassive BH model gives the largest increase in c₂₀₀(+9.8), which disfavors the model as an even more extreme outlier from the median ΛCDM relation (Section 9.1). Note, however, that the frequentist fit for this model increases c₂₀₀ by only ≈4; i.e., the evidence is mostly directed against the Bayesian version of the Fixed Overmassive BH model (also see above in Section 7.3).

The Einasto model and Weak AC Einasto models behave analogously to the NFW versions, but c₂₀₀ is shifted lower by ∼3. We reiterate that these various DM models cannot be distinguished in terms of their fit quality (Section 7.5).

All of the results described in this section, in particular the evidence that Mrk 1216 is an extreme outlier in the ΛCDM c₂₀₀–M₂₀₀ relation, are very consistent with the fossil galaxy/group NGC 6482 (B17). In Section 9.1 we discuss further implications of the high value measured for c₂₀₀.

In Table 11 we list the mass-weighted slope ($\langle \gamma \rangle$) and DM fraction (f_DM) for the Bayesian fiducial model evaluated at radii for several multiples of R_e. For r = R_e we obtain $\langle \gamma \rangle$ = 2.04 ± 0.05 and f_DM = 0.34 ± 0.05, which differ significantly from the average values obtained by Y17 for their sample of 16 CEGs ($\langle \gamma \rangle =2.3$, f_DM = 0.11). We discuss the principal reason for the discrepancy between our results and those of Y17 in Section 9.2. We note that the AC models give even larger DM fractions for Mrk 1216; e.g., the Weak AC+Einasto model result shown in Table 11.

Table 11.  Mass-weighted Total Density Slope and DM Fraction

Radius	Radius	Fiducial		Weak AC + Einasto
(kpc)	(Re)	$\langle \gamma \rangle$	fDM	$\langle \gamma \rangle$	fDM
1.1	0.5	2.12 ± 0.04	0.18 ± 0.04	2.05 ± 0.04	0.27 ± 0.05
2.3	1.0	2.04 ± 0.05	0.34 ± 0.05	2.02 ± 0.04	0.43 ± 0.06
4.6	2.0	1.93 ± 0.03	0.54 ± 0.05	1.99 ± 0.02	0.60 ± 0.05
9.2	4.0	1.91 ± 0.02	0.71 ± 0.03	1.97 ± 0.02	0.74 ± 0.03
11.5	5.0	1.94 ± 0.03	0.76 ± 0.03	1.98 ± 0.02	0.78 ± 0.03
23.0	10.0	2.08 ± 0.05	0.85 ± 0.01	2.02 ± 0.04	0.87 ± 0.02

Note. The mass-weighted slope ($\langle \gamma \rangle$) is evaluated for the Bayesian HE models using Equation (2) of Dutton & Treu (2014). The DM fraction is defined at each radius r as ${f}_{\mathrm{DM}}={M}_{\mathrm{DM}}(\lt r)/{M}_{\mathrm{total}}(\lt r)$. The Weak AC Einasto model is defined in the notes to Table 10.

Download table as:  ASCIITypeset image

The mass-weighted slope we measure for Mrk 1216 within r = R_e is less than the average slopes of local massive ETGs (2.15 ± 0.03, intrinsic scatter 0.10) determined by stellar dynamics (Cappellari et al. 2015). The local ETGs also have a smaller DM fraction (0.19, accounting for the higher mass range of the CEGs –see Y17). These differences between the local "normal" ETGs and Mrk 1216 presumably reflect the unusual and very high c₂₀₀ we measure for Mrk 1216; i.e., higher DM concentration means more DM near the center (higher f_DM), which translates into a smaller slope because the NFW (or Einasto) profile is weighted higher with a flatter slope than the stellar profile.

Studies of ETGs that combine strong lensing with stellar dynamics have found a large range of DM fractions within R_e, including several with values consistent with (or larger than) what we have found for Mrk 1216 (e.g., Barnabè et al. 2011; Sonnenfeld et al. 2015). The large DM fractions in the lensing/SD studies probably do not reflect unusually high c₂₀₀ like Mrk 1216; i.e., those galaxies were not selected to be early forming objects like Mrk 1216 (or NGC 6482), and thus they should not often possess c₂₀₀ values that deviate extremely from the ΛCDM c₂₀₀–M₂₀₀ relation.

Recently, Wang et al. (2018) have used the IllustrisTNG simulation to show that massive ETGs in ΛCDM have $\langle \gamma \rangle$ = 2.003 with a scatter of 0.175 over the radial range 0.4–4 R_e, with the slope changing little for redshifts below 2. The slope we measure for Mrk 1216 is consistent with these results (Table 11). Using the same simulation, Lovell et al. (2018; see their Figure 12) obtain f_DM within 5 R_e consistent with our measurements, but within 1 R_e, they find large values (f_DM ∼0.6). Using a different simulation, Remus et al. (2017) find a large range of DM fractions within 1 R_e consistent with what we find for Mrk 1216. Interestingly, both Remus et al. (2017) and Lovell et al. (2018) obtain smaller DM fractions within 1 R_e at z = 2; e.g., f_DM ≈ 10% (Remus et al. 2017). If Mrk 1216 is truly a nearby analog of a z ≈ 2 galaxy, its higher DM fraction is in conflict with these simulations.

In Humphrey & Buote (2010) we showed that the total mass slope approximated by a power law between 0.2 and 10 R_edecreases with halo mass and R_efor halos ranging from massive galaxies (∼10¹² M_⊙) up to massive clusters (∼10¹⁵ M_⊙) with the mean relation

$$\begin{eqnarray*}&&\gamma =2.31-0.54\mathrm{log}({R}_{e}/\mathrm{kpc}).\end{eqnarray*}$$

This slope–R_e relation predicts γ = 2.11 for Mrk 1216. The values we obtain within 10R_e for the fiducial and Weak AC+Einasto models listed in Table 11 are consistent with the relation within the observed scatter (0.14 dex, Auger et al. 2010).

We have also interpreted our HE analysis in terms of modified Newtonian dynamics (MOND; Milgrom 1983) to compare to our traditional Newtonian analysis. Our method follows Sanders (1999) and Angus et al. (2008) and is most recently described in detail in Section 6.7 of B17. In Figure 7 we display the DM fraction profile within a radius of 25 kpc for our fiducial Bayesian HE model computed from the Newtonian and MOND perspectives. As is clear in the figure, MOND requires almost as much DM as in the Newtonian approach, very similar to what we found for NGC 6482 (B17).

Zoom In Zoom Out Reset image size

The MOND acceleration scale a₀ ≈ 1.2 × 10⁻⁸ cm s⁻² is reached at a radius ≈11 kpc (considering only the mass in baryons—stars+SMBH+gas, ∼42 kpc otherwise). Even considering a much smaller central DM contribution (i.e., setting c₂₀₀ = 10 in the Newtonian analysis, too small to be consistent with the data—see Section 9.2) with correspondingly much larger stellar mass contribution (M_⋆/L_H ∼ 1.6 solar) only changes this radius by ∼1 kpc; i.e., the need for DM in MOND occurs well within the Newtonian regime of Mrk 1216. Hence, Mrk 1216 and NGC 6482 provide strong evidence that MOND requires DM on the massive galaxy scale, extending to lower masses the results obtained from HE studies on the group (e.g., Angus et al. 2008) and cluster (e.g., Pointecouteau & Silk 2005) scales.

Given the close similarity between MOND and the radial acceleration relation (RAR; Lelli et al. 2017), it is also to be expected that these galaxies deviate significantly from the RAR. Indeed, as shown in Figure 7, we find that the gravitational acceleration we derive from our HE analysis (g_obs) for Mrk 1216 (and NGC 6482) significantly exceeds the acceleration from only the baryons g_bar. For example, at the MOND acceleration scale, g_bar = 1.2 × 10⁻¹⁰ m s⁻², the RAR predicts g_obs = 1.9 × 10⁻¹⁰ m s⁻², whereas we measure for Mrk 1216g_obs = (4.99 ± 0.15) × 10⁻¹⁰ m s⁻². Lelli et al. (2017) quote a scatter in the RAR of ≲0.13 dex, indicating that at this one data point the discrepancy is >3σ. In fact, as seen in Figure 7, this level of discrepancy applies over a wide range in g_bar, and thus the significance of the discrepancy with the entire RAR is in fact much larger.

Finally, we mention that in Figure 7 we also show g_obs for the c₂₀₀ = 10 model mentioned above, which has a lower central DM and a higher baryon contribution for Mrk 1216. Although the shape of the g_obs profile is different from the fiducial model, the level of discrepancy with the RAR is broadly similar. The fractional errors (not shown) for g_obs of this model are similar to the fiducial model.

We display the radial profiles of the gas (f_gas) and baryon fractions (f_b) in Figure 8 for the fiducial Bayesian HE model and list the results for these quantities within radii r₂₅₀₀, r₅₀₀, and r₂₀₀ in Table 12 along with the systematic error budget. Within r₂₅₀₀, which essentially represents the extent of the Chandra data, we measure f_b,2500 = 0.071 ± 0.006; i.e., ≈45% of the cosmic mean value f_b,U = 0.155 determined by Planck (Planck Collaboration et al. 2014), where the hot gas contributes ≈40% of the measured baryons. The gas and baryon fractions continue to increase with radius, until at r₂₀₀, we have f_b,200 = 0.120 ± 0.016, comprising nearly 80% of the comic mean. Here the hot gas contributes ≈80% of the measured baryons.

Zoom In Zoom Out Reset image size

Table 12.  Gas and Baryon Fraction

	fgas,2500	fb,2500	fgas,500	fb,500	fgas,200	fb,200
Best Fit	0.028 ± 0.003	0.071 ± 0.006	0.058 ± 0.008	0.089 ± 0.010	0.094 ± 0.015	0.120 ± 0.016
(Max Like)	(0.027)	(0.071)	(0.058)	(0.089)	(0.094)	(0.120)
1 Brk Entropy	0.001	0.001	0.004	0.004	0.007	0.008
BH Flat Prior	0.000	0.000	0.001	0.001	0.002	0.002
BH Flat Logspace Prior	−0.000	0.001	−0.000	0.000	−0.000	0.000
Fixed Overmassive BH	0.002	−0.010	0.007	−0.001	0.013	0.007
Fixed MBH–σ BH	0.001	−0.002	0.002	0.000	0.004	0.002
Sérsic Stars 2MASS	−0.001	−0.003	−0.004	−0.006	−0.008	−0.010
Einasto	−0.002	−0.009	−0.005	−0.010	−0.009	−0.012
Strong AC	−0.003	−0.022	−0.011	−0.025	−0.020	−0.033
Weak AC	−0.002	−0.008	−0.006	−0.011	−0.011	−0.015
Weak AC Einasto	−0.004	−0.016	−0.012	−0.022	−0.022	−0.031
Joint Fit of Cycle 19 Obs.	−0.000	0.001	0.001	0.002	0.003	0.004
Constant Zα/ZFe	−0.001	0.002	−0.004	−0.002	−0.007	−0.005
Annulus 10 ΔZFe	${}_{-0.002}^{+0.002}$	${}_{-0.004}^{+0.003}$	${}_{-0.006}^{+0.005}$	${}_{-0.007}^{+0.006}$	${}_{-0.010}^{+0.009}$	${}_{-0.011}^{+0.010}$
Deproj	0.010	0.007	0.023	0.022	0.040	0.038
Distance	0.003	0.006	0.005	0.007	0.007	0.008

Note. Best-fit values and 1σ error estimates for the gas and baryon fractions of the fiducial Bayesian hydrostatic equilibrium model quoted for several overdensities. See the notes to Table 10 regarding the other systematic error tests.

Download table as:  ASCIITypeset image

Thus far, we have only considered the stellar baryons associated with the MGE decomposition of the HST H-band light from Y17. Yıldırım et al. (2015) note that only two galaxies are known within a 1 Mpc radius of Mrk 1216. Using NED, we identify these galaxies as 2MASX J08284832-0704316 (PGC152584) and 2MASX J08291551-0647454 (PGC152635), which on the sky are located ≈220 kpc (≈0.9r₅₀₀) and ≈300 kpc (≈0.9r₂₀₀), respectively, from the center of Mrk 1216. There is sparse information on these galaxies. However, based on their 2MASS K-band magnitudes, they are each about 2 mag fainter than Mrk 1216, which is suggestive of a fossil group. If we were to assume a similar contribution of non-central baryons as for the fossil group NGC 6482, this would add ≈0.04 to the baryon fraction at r₂₀₀ to give f_b,200 ≈ 0.16, consistent with the cosmic mean value.

Most of the systematic errors we have considered (Table 12) shift the gas and baryon fractions by less than the statistical error. The spectral deprojection ("Deproj" test) shifts are not significant within the larger 1σ erros in the deprojection model; e.g., 1σ error is ±0.041 on f_b,200. The AC models have the largest effect, resulting in lower values. The Weak AC Einasto model, which is our favored model (see Section 9.1), yields f_b,200 lower by 0.036, which effectively cancels the expected increase from non-central baryons.

We conclude that the Bayesian HE models predict that f_b,200 is close to the cosmic mean value for Mrk 1216. In Section 9.4 we discuss this measurement in relation to previous X-ray studies of other fossil elliptical galaxies and consider the implications for the "Missing Baryons" problem.

In Figure 9 we plot t_c/t_ff, the ratio of cooling time to freefall time, as well as t_c itself, as a function of radius for the fiducial Bayesian HE model. Mrk 1216 displays a large central region of approximately constant t_c/t_ff; i.e., for r ≤ 10 kpc, t_c/t_ff ≈ 14–19. Outside of this region, t_c/t_ff increases with radius. The t_c/t_ff profile of Mrk 1216 resembles a scaled-down version of the massive cluster Hydra-A (Hogan et al. 2017), which has pronounced X-ray cavities associated with AGN radio jets. Similarly, the observed t_c/t_ff profile more closely resembles massive elliptical galaxies and groups with evidence for multiphase gas (see Figure 2 of Voit et al. 2015). The t_c/t_ff profile suggests "precipitation-regulated AGN feedback" (e.g., Voit et al. 2017) for Mrk 1216, and we discuss further the implications of Mrk 1216 for feedback models in Section 9.5.

Zoom In Zoom Out Reset image size

We also remark that the approximately constant region of t_c/t_ff ≈ 17 ends near the NFW scale radius, r_s ≈ 12 kpc (and for t_c ≈ 10⁸ yr), for the fiducial model. We obtain extremely similar t_c/t_ff profiles for the Einasto DM model (r_s ≈ 14 kpc) and other models (e.g., AC); i.e., in Mrk 1216, the DM scale radius approximately equals the feedback radius.

As noted in Section 3.1, the centroid of the X-ray brightness changes by very little within the central ∼20''. For our analysis we adopted the center position obtained by computing the centroid within a circle of radius 3'' placed initially on the nominal stellar galaxy position. The annuli used for spectral extraction (Table 6) were defined about this centroid position (8:28:47.141, −6:56:24.367).

To gauge the sensitivity of our results to this choice, we also examined using the center position (8:28:47.131, −6:56:24.047) located at the emission peak ≈035 to the NW. Using this center has no measurable effect on the results.

In Section 4.2.2 we discussed the results obtained for the spectral fitting when a deprojection analysis is used. One issue that requires clarification is the error bars reported in Table 15. Because of the correlations between annuli introduced by the deprojection procedure, we estimate errors on the gas parameters (temperature, normalization, and abundances) via Monte Carlo simulations in xspec, as we have done in previous studies of deprojected spectra (e.g., Buote 2000a; Buote et al. 2003a); i.e., we performed 100 simulations of each set of spectra, computed the gas properties (e.g., temperature) for each simulation, and then computed the standard deviations on the parameters, which we quote as the 1σ errors in Table 15.

For our analysis of the projected spectra, we have also computed parameter errors by simply finding the values of a given parameter that change the C-statistic by 1 from its minimum value; i.e., using the error command in xspec. (We quote these results by default; e.g., Table 6). As expected, for the projected spectral analysis, we find that both approaches—ΔC = 1 and Monte Carlo—give very consistent results.

As mentioned in Section 4.2.2, we performed deprojected spectral fits for two cases: (1) no gas emission was accounted for outside of the bounding annuli; and (2) gas emission predicted by the Best Fit fiducial HE model from the projected spectral analysis (Table 8) was assigned to the large background apertures exterior to the bounding annuli. The results for case (2) are presented in the "Deproj" entries in Tables 10 and 12.

For most parameters, the "Deproj" systematic test shifts the values by an amount comparable to the 1σ statistical error. As noted in Section 7.9, the listed positive shifts and the gas and baryon fractions are not statistically significant when the larger statistical uncertainties of the deprojection models are considered. Case (1) produces larger shifts that are modestly significant; e.g., f_b,200 is increased by 0.063 ± 0.038. However, we prefer not to emphasize case (1) given that it does not account for emission projected into the bounding annuli. Finally, we also mention that the minimum χ² achieved for the frequentist HE analysis for both deprojected cases is ≈8 times larger than the projected fit; i.e., the fits are formally acceptable for the deprojected cases, but the default projected analysis gives a better fit.

We conclude that the constraints on the mass profile obtained when using the spectral deprojection approach are overall very consistent with the default projection analysis.

When we do not allow the abundance ratios Z_Mg/Z_Fe and Z_Si/Z_Fe to vary with radius in the spectral fits (Section 4.2.1), the effects on the HE models are listed in the "Constant Z_α/Z_Fe" entries in Tables 10 and 12. This test shifts most of the parameters by an amount comparable to the 1σ errors, and in the case of M_⋆/L_H, by ≈1.5σ. Interestingly, the results for this test differ noticeably for the Best Fit and Max Like values. The shifts of marginal significance listed in the tables reflect only the Best Fit values. If we compare Max Like values instead, then all the shifts are <1σ.

When we allow Z_Fe to take values 0.27 Z_⊙ and 0.45 Z_⊙ for Annulus 10 of the Cycle 19 observation (see Section 4.2.1), representing the intrinsic scatter of the average group/cluster profile of Mernier et al. (2017), we list the parameter shifts for the HE models in the "Annulus 10 ΔZ_Fe" entries in Tables 10 and 12. All the quoted parameter shifts are smaller than the statistical errors.

Throughout the paper we have quoted global parameters (e.g., concentration and mass) for a series of "virial" radii representing three common overdensities; i.e., r₂₅₀₀, r₅₀₀, and r₂₀₀ (see Section 7), where r₂₅₀₀ essentially matches the radial extent of the Cycle 19 observation, while r₅₀₀ ≈ 1.9r₂₅₀₀ and r₂₀₀ ≈ 2.8r₂₅₀₀ fall outside the observed data range. Our analysis of the projected data technically constrains the projected emission for radii outside the data range in projection (i.e., r > r₂₅₀₀). In reality, the emission from radii much larger than the data extent projected onto the observed sky annuli is dominated by the emission from within the three-dimensional spherical shells corresponding to the radii of the observed sky annuli. Hence, the global parameter values quoted at radiii r₅₀₀, and particularly, r₂₀₀ are to an increasing amount extrapolations of our model outside the data range.

We do not expect significant systematic errors in the extrapolated parameter values (e.g., c₂₀₀, and M₂₀₀) provided (1) the true DM profile is accurately described by the NFW/Einasto models and (2) we have accurately measured the DM scale radius r_s. In Gastaldello et al. (2007) we have previously emphasized that accurate and precise constraints on r_s, and thus the NFW DM profile from HE X-ray studies are only possible when several radial data bins exist both above and below r_s. For Mrk 1216, we measure r_s = 12.1 ± 2.2 kpc, so that approximately six annuli lie below r_s and four annuli above r_s for the Cycle 19 data, with R_out ≈ 9r_s where R_out =110.7 kpc is the outer radius of Annulus 10. (The Cycle 16 data contribute four annuli below r_s, two annuli above r_s, and one annulus encloses r_s.) Because r_s is situated well within the data range with several radial bins above and below it, we believe it is well constrained, as is reflected in the Bayesian constraints; e.g., for the fiducial HE model, we obtain r_s = ${12.1}_{-4.4}^{+7.1}$ kpc (99% confidence).

It is instructive to consider the fossil cluster RXJ 1159+5531 for which the accuracy of some radial extrapolation (i.e., from r₅₀₀ to r₂₀₀) has been directly tested. Although it is ≈15 times more massive, RXJ 1159+5531 is similar to Mrk 1216 in that it is a highly relaxed system with an above-average halo concentration. Using a single Chandra observation with data extending out to ≈r₅₀₀ with approximately 5 (4) data bins below (above) r_s, Gastaldello et al. (2007) obtained c₂₀₀ = 8.3 ± 2.1 and M₂₀₀ = (7.9 ± 5) × 10¹³ M_⊙ with r_s ≈104 kpc. Later, adding several Suzaku observations covering the entire sky region out to r₂₀₀, in Buote et al. (2016), we obtained c₂₀₀ = 8.4 ± 1.0 and M₂₀₀ = (7.9 ±0.6) × 10¹³ M_⊙ with r_s ≈ 104 kpc, extremely consistent with the previous extrapolated measurement. Note in particular that r_s did not change between the studies.

Based on a suite of cosmological hydrodynamical cluster simulations, Rasia et al. (2013) reported that significant positive biases (up to 25%) in the inferred value of c₂₀₀ are possible if the X-ray HE analysis only fits data within ≈r₂₅₀₀. However, the largest biases occur for their simulations without AGN feedback. When AGN feedback is included, the bias is small (i.e., ≲10%, see Figure 5 of Rasia et al. 2013). The lowest mass clusters studied by Rasia et al. (2013) have M₂₀₀ ≈ 2 × 10¹⁴; i.e., ≈40 times more massive than Mrk 1216. Assuming that these results apply to lower-mass halos like Mrk 1216, then the ≲10% bias for the simulations with AGN feedback is comparable to or lower than the statistical error we measure for Mrk 1216.

In a future paper we will test the extrapolation for Mrk 1216 using a new deep XMM-Newton observation that will extend the radial range of the model fits closer to r₅₀₀.

It has been shown that spherical averaging of an ellipsoidal mass profile typically introduces only small orientation-averaged biases for global quantities such as c₂₀₀, M₂₀₀, and f_gas,200 in X-ray HE studies of galaxy and cluster masses (e.g., Buote & Humphrey 2012b, 2012c, and references therein). Because we also found such small biases to be negligible compared to the statistical errors for the massive elliptical galaxy NGC 6482 (see B17), which has DM properties very similar to Mrk 1216, we do not present a specific estimate here.

We have also considered in this paper the mass profile near R_e. Previous studies of the spherical approximation in X-ray studies of galaxy clusters do not specifically address such a small radius. However, in Buote & Humphrey (2012b; see also Churazov et al. 2008), we showed that spherically averaging an ellipsoidal scale-free logarithmic potential in an HE analysis introduces zero bias in the inferred mass for any gas temperature profile. We expect this result to apply to a good approximation near R_e, where it is well established that the total mass profile slope $\langle \gamma \rangle \approx 2$ for massive ETGs, as we have found for Mrk 1216 (see Section 7.7).

We can test this expectation for the massive elliptical galaxy NGC 720, for which we previously reported a value M_tot/L_B = 6.0 solar at R_e obtained from spherically averaging ellipsoidal HE model fits to the Chandra data (Buote et al. 2002). We compare this measurement to the results of the fully spherical HE analysis of NGC 720 we performed using the same Chandra data in Humphrey et al. (2006) after accounting for the slightly different distance used and the fact the latter study quotes M_tot/L_K; i.e., M_tot/L_B = 6.0 solar becomes M_tot/L_K = 1.1 solar, in excellent agreement with the M_tot/L_K profile displayed in Figure 5 of Humphrey et al. (2006). We therefore expect that systematic errors associated with spherical averaging should also be small near 1 R_e for Mrk 1216.

To assess the sensitivity of our results to the shape of the stellar light profile, we also considered a de Vaucouleurs model (i.e., n = 4 Sérsic model) with the half-light radius inferred from the 2MASS Extended Source Catalog (Jarrett et al. 2000). We adopted a median, circularized R_e = 2.1 kpc and normalized the model to the total H-band luminosity (Table 2). We list the results for this test in the "Sérsic Stars 2MASS" entries in Tables 10 and 12. In all cases, using the Sérsic model shifts the paramers by ≤1σ.

In Table 13 we compare the halo concentration (c₂₀₀) and mass (M₂₀₀) we have measured for Mrk 1216 using the Bayesian HE models to some theoretical models of the ΛCDMc₂₀₀–M₂₀₀ relation from the literature. Our fiducial model with an NFW DM halo exceeds the median ΛCDM relation by ∼6σ in terms of the intrinsic scatter of the theoretical relations. The size of this discrepancy implies an outlier so extreme as to be inconsistent with the ΛCDM simulations. Indeed, from visual inspection of published catalogs of c₂₀₀ and M₂₀₀ values, we do not see any ΛCDM halos consistent with the NFW values we have measured for Mrk 1216; e.g., Figure 16 of Dutton & Macciò (2014) and Figure 2 of Ragagnin et al. (2018).

Table 13.  Comparison to Theoretical ΛCDMc₂₀₀–M₂₀₀ Relations

	Mrk 1216		Dutton+14		Ludlow+16		Child+18
DM Profile	M200	c200	${\bar{c}}_{200}$	Δc200	${\bar{c}}_{200}$	Δc200	${\bar{c}}_{200}$	Δc200
NFW	(5.3 ± 0.6) × 1012 M⊙	30.4 ± 4.3	7.0	5.8σ	⋯	⋯	6.1	(11.9σ, 6.3σ)
Weak AC NFW	(6.1 ± 1.0) × 1012 M⊙	23.9 ± 4.1	6.9	4.9σ	⋯	⋯	6.1	(9.0σ, 5.5σ)
Einasto	(5.8 ± 1.0) × 1012 M⊙	27.2 ± 4.3	7.9	4.1σ	7.2	4.5σ	7.1	(8.4σ, 4.5σ)
Weak AC Einasto	(7.2 ± 1.5) × 1012 M⊙	20.8 ± 4.0	7.7	3.4σ	7.1	3.7σ	7.0	(6.1σ, 3.7σ)

Note. Results listed for Mrk 1216 pertain to the Best Fit Bayesian HE models. For each theoretical c₂₀₀–M₂₀₀ relation, we list (1) ${\bar{c}}_{200}$, the median value of c₂₀₀ predicted for the Best Fit M₂₀₀ for Mrk 1216, and (2) ${\rm{\Delta }}{c}_{200}={c}_{200}-{\bar{c}}_{200}$ expressed in terms of the intrinsic scatter of the theoretical relation. The theoretical relations are as follows. (1) Dutton+14 (Dutton & Macciò, 2014). ΛCDM Planck cosmological parameters. (2) Ludlow+16 (Ludlow et al. 2016). ΛCDM Planck cosmological parameters with the same intrinsic scatter as Dutton+14. (3) Child+18 (Child et al. 2018). ΛCDM WMAP-7 cosmological parameters. We use the results for the stacked halos quoted in their Table 1. We list two values for Δc₂₀₀. The first assumes a Gaussian scatter σ_c = c/3 preferred by Child+18. The second uses the lognormal scatter of Dutton+14.

Download table as:  ASCIITypeset image

When substituting the Einasto DM profile for NFW the value obtained for c₂₀₀ remains high, although the tension with ΛCDM is reduced modestly to 4.1–4.5σ (Table 13). The reduced tension for the Einasto model is attributed both to the smaller value of c₂₀₀ measured for Mrk 1216 and the larger theoretical intrinsic scatter compared to the NFW profile (0.13 dex, Dutton & Macciò 2014). Nevertheless, we do not see any ΛCDM halos consistent with the Einasto values we have measured for Mrk 1216 in the simulation results published by Benson et al. (2019, see their Figures 2 and 4). AC lowers c₂₀₀ more and further lessens the tension with ΛCDM to a more modest but still substantial 3.4–3.7σ.

In Figure 10 we plot the c₂₀₀ and M₂₀₀ values for Mrk 1216 and three fossil and nearly fossil systems with evidence for above-average concentrations: NGC 720 (Humphrey et al. 2011), NGC 1521 (Humphrey et al. 2012), and NGC 6482 (B17). We refer to NGC 720 and NGC 1521 as "nearly" fossil systems because within their projected virial radii, they obey the fossil classification. NGC 6482 in particular displays c₂₀₀ and M₂₀₀ behavior extremely similar to Mrk 1216–both in terms of its outlier status with respect to ΛCDM and the way the c₂₀₀ and M₂₀₀ parameters change between NFW, Einasto, and AC models.

Zoom In Zoom Out Reset image size

The very high halo concentration we have measured for Mrk 1216 has profound implications for its formation and evolutionary history. The halo formation time is a major factor determining the scatter in halo concentrations for a given mass (e.g., Neto et al. 2007; Ludlow et al. 2016), with higher concentrations reflecting halos that have formed earlier. Consequently, the high value of c₂₀₀ indicates that compared to the typical halo of its total mass, most of the mass of Mrk 1216 was in place much earlier, which is consistent with largely passive evolution for most of its existence. Thus the high concentration provides vital corroborating evidence of Mrk 1216 as a massive "relic galaxy" that is a largely untouched descendant of the red nugget population at high redshift (e.g., van den Bosch et al. 2015; Yıldırım et al. 2015, 2017; Ferré-Mateu et al. 2017).

While the connection between early formation and high c₂₀₀ is critical, other factors must also contribute. Ludlow et al. (2016) estimate that formation time accounts for a large fraction (∼50%) of the scatter in c₂₀₀, leaving about half the scatter to be likely explained by the initial conditions or the environment of a given halo. The recent study by Ragagnin et al. (2018) proposes that a measure of the "fossilness" of a halo can account for some of this scatter, although the proposed effect is too modest to explain the large c₂₀₀ of Mrk 1216. These issues are likely central to understanding how the DM and total mass profiles of Mrk 1216 and NGC 6482 are extremely similar, but the stellar component of Mrk 1216 is still more compact, obeying the size-mass relation for z ∼ 2 galaxies like other CEGs.

Finally, we note that we have emphasized comparisons with the theoretical c₂₀₀–M₂₀₀ relations assuming a lognormal intrinsic scatter. While a Gaussian interpretation of the theoretical scatter is consistent with cosmological simulations (e.g., Bhattacharya et al. 2013), the narrower tails of the Gaussian imply that Mrk 1216 would be an outlier so extreme (see Child+18 results in Table 13) for all models to be incompatible with ΛCDM. Consequently, the high c₂₀₀ values for Mrk 1216 and NGC 6482 (B17) clearly favor a lognormal distribution for the intrinsic scatter of the ΛCDMc₂₀₀–M₂₀₀ relation.

We have found significant differences in the mass properties of Mrk 1216 inferred from our HE analysis of the hot plasma within a radius $\sim {R}_{e}$ compared to the stellar dynamical (SD) constraints obtained by Y17 and Walsh et al. (2017). There is, however, a major difference in the assumed DM model between these SD studies and ours. Because the SD studies are unable to contrain the DM halo, they assume an NFW DM halo by default, with fixed c₂₀₀ = 10. Here we consider how our HE analysis changes when we make a similar assumption.

We achieve a similar DM halo in our analysis by fixing the NFW scale radius to r_s = 60 kpc, in which case we obtain c₂₀₀ = 9.8 ± 0.2 and M₂₀₀ = (2.3 ± 0.1) × 10¹³ M_⊙ from our Bayesian analysis. The fit quality, however, is very poor; i.e., the frequentist analysis gives χ² = 39.5 for 19 dof compared to χ² = 11.9 when c₂₀₀ =  is allowed to take a value c₂₀₀ = 30. In Figure 11 we plot the χ² residuals for the c₂₀₀ ≈ 10, 30 models, showing just the surface-brightness data for the Cycle 19 observation. Compared to the c₂₀₀ = 30 model, the c₂₀₀ = 10 model has much larger residuals over most of the radial range; i.e., the Chandra data clearly exclude the c₂₀₀ = 10 model in favor of the c₂₀₀ = 30 model.

Zoom In Zoom Out Reset image size

While the c₂₀₀ = 10 model does not provide a good fit to the data, the derived mass parameters agree very well with those of Y17.⁹ We obtain M_⋆/L_H = 1.60 ± 0.05 solar, which agrees with Y17's value, M_⋆/L_H = ${1.85}_{-0.40}^{+0.52}$, within their larger 1σ error bar. While Y17 do not quote values of the total mass slope and DM fraction within R_e specifically for Mrk 1216, their sample-averaged values ($\langle \gamma \rangle =2.3$, f_DM = 0.11) agree very well with what we obtain for Mrk 1216 ($\langle \gamma \rangle =2.27\,\pm 0.01$, f_DM = 0.13 ± 0.01) when we assume c₂₀₀ = 10. Furthermore, when we also fix the SMBH mass to the best-fitting "overmassive" value (M_BH = 4.9 × 10⁹ M_⊙) obtained by Walsh et al. (2017), then we obtain an even worse fit (χ² = 50.3 for 20 dof), but with M_⋆/L_H = 1.36 ± 0.05 solar, in excellent agreement with the value M_⋆/L_H = 1.3 ± 0.4 solar obtained by Walsh et al. (2017).

We conclude that the assumption of c₂₀₀ = 10 for the NFW DM halo fully accounts for the higher M_⋆/L_H measured with SD by Y17 and Walsh et al. (2017) as well as the higher $\langle \gamma \rangle$ and lower f_DM within R_e obtained by Y17. Differences in the form of the DM halo, however, have negligible impact on our constraints on M_BH, which are manifested only in the central 1'' aperture (i.e., Annulus 1 of the Cycle 19 data). Below in Section 9.3 we evaluate the SMBH constraint further and consider its implications.

In Section 7.3 we show that the HE analysis does not favor the best-fitting overmassive SMBH [M_BH = (4.9 ± 1.7) ×10⁹ M_⊙] found by the SD study of Walsh et al. (2017), but does not clearly exclude the value either when different Bayesian priors and the frequentist fits are considered. In addition, the SD error estimate on M_BH permits lower values (e.g., 1.4 × 10⁹ M_⊙ (2σ) and 0.65 × 10⁹ M_⊙ (2.5σ)) that are very consistent with our HE models. In sum, the estimated statistical and systematic errors in the X-ray HE and SD constraints allow for broad consistency and do not clearly point to whether the SMBH exceeds the M_BH–σ relation.

Based on models of the Cycle 16 data of Mrk 1216 assuming M_BH = 4.9 × 10⁹ M_⊙, we expected to measure a value of k_BT about 20% higher in Annulus 1 than we found with the Cycle 19 observation (see beginning of Section 7.3). If we attribute this underestimate to deviations from HE, it would imply ∼40% nonthermal pressure support in Annulus 1 from random turbulent motions, magnetic fields, rotational support, etc. Future observations with similar or better spatial resolution to Chandra but with much higher spectral resolution and sensitivity (i.e., Lynx)¹⁰ will be necessary to establish the hot gas kinematics within the central ≈1''. In the meantime, more precise constraints on M_BH from SD can better establish the significance of any tension with the X-ray constraints.

Our Bayesian HE models evaluated at r₂₀₀ indicate f_b,200 ≈ 80% of the cosmic mean value (Section 7.9). The evidence for near baryonic closure in Mrk 1216 is consistent with results we have obtained previously for other fossil-like massive elliptical galaxies, NGC 720 (Humphrey et al. 2011), NGC 1521 (Humphrey et al. 2012), and NGC 6482 (B17). These results suggest that in massive elliptical galaxy/small group halos, many of the "missing" baryons (Fukugita et al. 1998) can be located in the outer halo as part of the hot component.

As is clear from Figure 8, most of the baryons inferred from our HE models within a radius of r₂₀₀ are in the form of hot gas located at radii larger than the extent of the Chandra data (≳r₂₅₀₀). We have previously noted (Section 8.2 of B17) the sensitivity of the measured f_gas,200to the assumed value of Z_Fe for r > r₂₅₀₀. In fact, using a β-model analysis of the Chandra surface-brightness profile of NGC 720, Anderson & Bregman (2014) claim that the gas mass (M_gas,200) is a factor 3–4 lower than we measured in Humphrey et al. (2011).

We have reanalyzed the Chandra surface-brightness profiles using all currently available observations and performed a similar β model analysis of NGC 720. We find M_gas,200 values very consistent with those of Humphrey et al. (2011). We attribute the lower value of M_gas,200 obtained by Anderson & Bregman (2014) to primarily (1) their quoting a higher value of M_gas,200 than was actually reported in Figure 6 of Humphrey et al. (2011); i.e., 4 × 10¹¹ M_⊙ instead of the correct 3 × 10¹¹ M_⊙ at r = 300 kpc; (2) their using a different solar abundance standard (Anders & Grevesse 1989), whereas Humphrey et al. (2011) use Asplund et al. (2006); and (3) their assuming k_BT = 0.5 keV, Z_Fe = 0.6 solar for the hot plasma appropriate for the center of NGC 720 instead of quantities appropriate near r₂₅₀₀ (k_BT ≈ 0.3 keV, Z_Fe ≈ 0.3 solar).

To verify the nearly cosmic baryon fractions in these galaxies, it is necessary to measure both k_BT and Z_Fe out to radii beyond r₂₅₀₀. In future papers we will report the results of new deep XMM-Newton observations of Mrk 1216 and NGC 6482 that will allow mapping the gas out closer to ∼r₅₀₀.

The Chandra observations of Mrk 1216 suggest the hot plasma properties within the central ∼10 kpc reflect a balance between radiative cooling and episodes of gentle AGN feedback. The centrally peaked radial temperature profile itself (Figure 5) does not implicate AGN shock-heating of the gas because the continuous injection of stellar material will produce the same effect (e.g., classical wind models, David et al. 1991; Ciotti et al. 1991) and the observed central entropy is not peaked (Figure 6). Modern models of massive elliptical galaxies incorporating gentle mechanical AGN feedback episodes can produce similar centrally peaked temperature profiles (e.g., Gaspari et al. 2012a). The lack of large, significant disturbances (e.g., cavities) near the center (Section 3) or azimuthal spectral variations (Section 4.2.3) provide further evidence that episodic AGN feedback in Mrk 1216 has been gentle.

Interior to the sizable radius ∼10 kpc within which t_c/t_ff ranges from 14 to 19 (Figure 9, see Section 7.10), the "precipitation-regulated feedback" scenario (e.g., Voit et al. 2017) predicts there should be multiphase gas, but there is currently little evidence for gas cooler than the ambient hot plasma. However, Voit et al. (2017) emphasizes that condensation will be suppressed by buoyancy damping if the entropy profile is steeply rising and without a large isentropic core. We find that the entropy profile has a radial logarithmic slope α = 0.78 ± 0.05 for 0.66 < r < 16 kpc (Table 8) and no evidence of a core. This slope is larger than the critical value α_K = 2/3 for suppressing condensation derived by Voit et al. (2017) using a toy model for isothermal gas.

In fact, this result may be common for massive X-ray luminous elliptical galaxies. While the recent study by Babyk et al. (2018) finds a universal entropy profile for elliptical galaxies and galaxy clusters with an entropy slope consistent with two-thirds within 0.1r₂₅₀₀, our previous measurements within the central ∼10 kpc of massive elliptical galaxies have obtained steeper slopes α ≳ 1, consistent with suppressed condensation: NGC 4649 (Humphrey et al. 2008), NGC 1332, NGC 4261, NGC 4472 (Humphrey et al. 2009a), and NGC 6482 (B17).

Although the α_K = 2/3 criterion for suppressing condensation from the toy model of Voit et al. (2017) is consistent with the single-phase gas we observe in Mrk 1216, the three-dimensional hydrodynamical simulations of Wang et al. (2018) predict that t_c/t_ff should rise rapidly with radius for a single-phase gas, which contradicts the Chandra data (Figure 9) and the lack of evidence for multiphase gas in the amounts seen in other galaxies with t_c/t_ff < 20 (Werner et al. 2014; Voit et al. 2015).

In Section 7.10 we remarked that the region where t_c/t_ff ranges from 14 to 19 (i.e., r ≲ 10 kpc), associated with a steeply rising entropy profile with α > 2/3, also is close to the NFW DM scale radius r_s. For the galaxies listed above with good constraints on the DM halo, r_s is not far from 10 kpc; i.e., r_s ≈ 14 kpc for NGC 720 and r_s = 11.2 ± 3 kpc for NGC 6482.

We conclude by noting that the unusually quiescent evolutionary path indicated by the highly concentrated stellar and DM components of Mrk 1216 imply that the presumed episodic AGN feedback regulating the cooling in the central ∼10 kpc has not been upset by triggering from galaxy merging, but has been solely governed by accretion from the radiatively cooling hot plasma. We speculate that this has also resulted in unusually quiescent evolution in the hot plasma obeying the HE condition accurately (see below in Section 9.6).

Because the Chandra ACIS spectral data do not permit useful direct measurements of the kinematics of the hot plasma, we instead briefly summarize several indirect lines of evidence testifying to the accuracy of the HE approximation in Mrk 1216.

Azimuthal Scatter in Quadrants (Section 4.2.3): When Annuli 2,3, and 4 of the Cycle 19 observations are divided up into quadrants, we find that the small scatter of the pseudo-pressure and pseudo-entropy between quadrants in each annulus suggests HE deviations of ≲5%.

Quality of HE Model Fits (Table 9): We find that reasonable HE models provide excellent fits to the radial k_BT and Σ_x profiles of the Cycle 19 and Cycle 16 Chandra data.

Qualified Consistency with Stellar Dynamics Measurements (Section 9.2): If we fix c₂₀₀ = 10 for the NFW DM halo, as done in the stellar dynamical studies by Y17 and Walsh et al. (2017), we find excellent agreement with the measured M_⋆/L_H values as well as the total mass slope ($\langle \gamma \rangle$) and DM fraction (f_DM) obtained within 1R_e. The constraints on M_BH are uncertain and therefore ambiguous with respect to the status of HE in Annulus 1 of the Cycle 19 data (Section 9.3).

Stellar Mass-to-Light Ratio and SPS Models (Section 7.4): The value we measure for M_⋆/L_H agrees very well with SPS models with a Kroupa IMF.

High Halo Concentration (Section 9.1): The exceptionally high value of c₂₀₀ we measure for Mrk 1216 corroborates the evidence from its compact stellar size and old stellar population that it is a massive relic galaxy; i.e., Mrk 1216 has evolved passively in isolation since a redshift of ∼2 and, consequently, is highly evolved and relaxed.

We therefore expect that Mrk 1216 is a highly relaxed system, at least as relaxed as the massive Perseus cluster with its pronounced cavities and other surface-brightness irregularities for which Hitomi found with gas kinematics measurements that the pressure from turbulence is only 4% of the thermal gas pressure (Hitomi Collaboration et al. 2016). To directly verify the accuracy of the HE approximation in Mrk 1216 awaits the next generation of X-ray satellites with microcalorimeter detectors.^{11 ,12}