EDP Sciences
Free Access
Issue
A&A
Volume 607, November 2017
Article Number A98
Number of page(s) 16
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201629926
Published online 21 November 2017

© ESO, 2017

1. Introduction

Line cooling of chemical elements from C to Fe plays an important role in the formation of galaxies, stars, and planets. Most of the elements in the Universe today are thought to have formed in star bursts at z ≈ 2−3 (Hopkins & Beacom 2006; Madau & Dickinson 2014). The hot intracluster medium (ICM) in groups and clusters of galaxies is an excellent probe of this chemical evolution in the dense regions of the Universe. Metals are accumulated over very long times (>5 Gyr) in the cluster centers, and the total mass of metals in the hot plasma in the core is about a factor of 2–6 higher than the metal mass locked up in the galaxies (Renzini & Andreon 2014). The abundances measured in the plasma thus provide a “fossil” record of the integral yield of all the different stars (releasing metals in supernova explosions and winds) that have left their specific abundance patterns in the gas before and during cluster evolution (see, e.g., de Plaa 2013; Werner et al. 2008, for a review).

X-ray spectroscopy provides a precise measure of metal abundances in the ICM. The observations with the European Photon Imaging Camera (EPIC, Strüder et al. 2001; Turner et al. 2001) provide highly significant measurements of the abundances of Si, S, Ar, Ca, Fe, and Ni. The high-resolution X-ray spectra obtained with the Reflection Grating Spectrometer (RGS, den Herder et al. 2001) on board XMM-Newton, which resolves the Fe-L complex into individual lines, allow for precise abundance measurements of O, Ne, Mg, and Fe. In cooler clusters (3 keV), RGS also detects lines from N (Sanders & Fabian 2011). Because non-equilibrium ionization and optical depth effects in the ICM are very weak, these abundances are more reliable than abundances measured in stars or in the cold low-ionization interstellar medium. However, a thorough study of systematic uncertainties in abundance measurements with RGS has not been performed to date. The error bars in Fig. 1 show the expected statistical error bars on cluster abundances based on previous studies (de Plaa et al. 2007; de Plaa 2013).

Most of the elements detected with XMM-Newton are produced by supernovae. Core-collapse supernovae (SNcc) produce large amounts of O, Ne, and Mg (e.g., Woosley & Weaver 1995; Nomoto et al. 2006), while type Ia supernovae (SNIa) produce large quantities of Fe, Ni, and relatively little O, Ne, and Mg (e.g., Iwamoto et al. 1999; Bravo & Martínez-Pinedo 2012). The Si-group elements (Si, S, Ar, and Ca) are produced by both supernova types (see Fig. 1). N is produced mainly by asymptotic giant branch (AGB) stars (Karakas 2010; Werner et al. 2006a; Grange et al. 2011) and by winds of massive stars, especially in rotating stars and at low metallicity (e.g., Romano et al. 2010). In this paper, we focus mainly on the O/Fe abundance ratio. The O/Fe, Ne/Fe, and Mg/Fe ratios are good indicators for the relative contribution of SNIa with respect to SNcc. The knowledge of these ratios is important for determining the amount of Si-group elements produced by SNIa.

thumbnail Fig. 1

Expected abundances measured in a typical long XMM-Newton observation of 120 ks (bottom panel). The estimates for the SNIa, SNcc, and AGB contribution are based on a sample of 22 clusters (de Plaa et al. 2007) and two elliptical galaxies (Grange et al. 2011). The top panel shows the typical range in SNIa and IMF models with respect to the statistical error bars in the observation. Figure adapted from de Plaa (2013).

Open with DEXTER

SNIa have likely produced a substantial fraction of the Fe, Ni, and Si-group elements observed in cluster cores. The statistical precision of the abundance ratios derived from X-ray observations of nearby clusters and groups of galaxies in a typical XMM-Newton orbit of 120 ks is typically better than 10–20%, while the spread in yields (see Fig. 1, top panel) obtained from simulations assuming different SNIa explosion mechanisms can be up to a factor of a few for elements such as Ca and Ni (e.g., Iwamoto et al. 1999; Badenes et al. 2003). Accurate cluster abundances therefore allow us to constrain supernova models.

Much of the uncertainty in SNIa yields is due to the variety in possible type Ia supernova progenitors and the subsequent explosion mechanism. In recent years, the search for type Ia supernovae in galaxies has become much more efficient. Large samples of SNIa observed in mainly optical, infrared, and UV wavelengths revealed variations in SNIa properties that appear to correlate with the properties of the host stellar populations (see, e.g., Howell 2011; Wang & Han 2012, for a review). In the single-degenerate (SD) supernova scenario, a carbon-oxygen white dwarf accretes matter from a non-degenerate companion star before it reaches the critical temperature for explosive carbon ignition. It has become clear that the properties of the companion star are important for the properties of the type Ia explosion that follows after the accretion phase and is one of the origins of the variety of SNIa that is observed. In addition, in the double-degenerate (DD) scenario, two white dwarfs merge and disintegrate in a supernova explosion, creating yet another variety of type Ia supernovae.

In addition to our lack of knowledge about the explosion mechanism, it is also unclear how the progenitor systems form. Attempts have been made to explain the observed type Ia rate theoretically through simulations of the evolution of binary populations (e.g., Claeys et al. 2014). This study showed that if the SNIa rate is due to the standard SD channel, the SNIa rate can be explained only under the assumption that the accretion onto the white dwarf is not limited (e.g., that the Eddington limit does not hold). The result of this and similar studies makes clear that type Ia supernovae are still a poorly understood phenomenon. Abundance ratios determined from clusters are therefore a key test for binary population synthesis and SNIa explosion models.

A simple test has been performed by, for example, de Plaa et al. (2007) using a sample of 22 clusters observed with XMM-Newton. The authors analyzed the abundances of Si, S, Ar, Ca, Fe, and Ni within a radius of 0.2 R500 from the cluster center. A good fit was obtained with a one-dimensional delayed-detonation model from Badenes et al. (2003), while models from Iwamoto et al. (1999) were unsuccessful because they underestimated the Ca abundance. The model from Badenes et al. (2003) that fitted the cluster abundances in de Plaa et al. (2007) also fits the abundances of the Tycho supernova remnant (Badenes et al. 2006), which is thought to have been a fairly typical SNIa with an average luminosity. Recently, Mulchaey et al. (2014) suggested that a subclass of supernovae, called Ca-rich gap transients, may provide enough calcium to explain the high calcium abundance found in the ICM of clusters.

However, the work of de Plaa et al. (2007) only used abundances of elements heavier than Si determined from EPIC, and their O, Ne, and Mg measurements were determined from only two clusters for which they analyzed deep RGS spectra (de Plaa et al. 2006; Werner et al. 2006b). In order to distinguish the SNcc and SNIa contribution and place stronger constraints on the SNIa explosion mechanism, accurate knowledge of the abundances of these elements from a large number of clusters is necessary. Accurate measurements of O, Ne, and Mg, which are primarily products of SNcc, require the unique capabilities of the RGS.

The O, Ne, and Mg yields from SNcc and N from AGB stars depend strongly on the progenitor mass. The difference in the total population yields between a top-heavy or Salpeter initial mass function (IMF) corresponds to a spread of 80% in the N abundance, and a 3040% spread in O, Ne, and Mg (see Fig. 1). Unfortunately, N is not only produced in AGB stars, but also ejected in winds of massive stars before the SNcc explosion, especially in rotating stars and at low metallicity, which makes the models for the origin of N very uncertain (see, e.g., Romano et al. 2010). Since the accuracy of the measured abundances is higher, a large sample of clusters with a broad range of masses, cool-core properties, and optical characteristics of the central dominant (cD) galaxy may provide constraints on these models. Ultimately, if we are also able to measure carbon and sodium, the abundance sets may provide a test of the IMF universality as well.

This paper introduces the CHEmical Enrichment RGS Sample (CHEERS) project, which mainly aims to obtain reliable chemical abundances in the intracluster medium of galaxy clusters through deep XMM-Newton observations of 44 clusters. We introduce the sample and describe the selection of the clusters, which is optimized to exploit the RGS to the best of its abilities. The observations required to complete this sample were performed in AO-12 as part of an XMM-Newton Very Large Program. This paper reports the abundance results for oxygen and iron obtained from the RGS spectra. Because of the high spectral resolution in the soft X-ray band, the RGS is better capable of resolving the oxygen lines than EPIC. We aim to provide a thorough discussion about the reliability of the measurements and the robustness of the cluster to cluster variations. This sample also provides very high quality XMM-Newton EPIC data. In two companion papers, we describe the EPIC abundance measurements of the other common elements using this sample (Mernier et al. 2016a) and the interpretation of the combined RGS and EPIC abundances (Mernier et al. 2016b). The radial abundance profiles are studied in Mernier et al. (2017). In another companion paper, we report detections of nitrogen in a subset of the RGS observations (Mao et al. 2017). Our measured abundances are relative to the proto-solar abundances by Lodders et al. (2009) unless stated otherwise. Error bars are given at the 1σ (68%) confidence level.

2. Sample selection

In order to study the chemical enrichment history of individual clusters of galaxies and the differences in enrichment between clusters, we need a moderately large sample of clusters with deep exposure times per cluster. Because the sample of de Plaa et al. (2007) lacked sufficient RGS coverage, we need to expand this sample to be able to divide it into subsamples of different cooling properties and study the spatial distribution of the elements. Our aim is to have a “complete” sample of high-quality RGS cluster spectra that can be obtained within a reasonable exposure time of 200 ks each. With “complete”, we mean that we aim to have observations of all suitable RGS cluster targets within a redshift of z = 0.1. Obviously, many clusters already have deep RGS spectra, but the XMM-Newton archive did not contain deep observations of all the suitable targets. We obtained deep XMM-Newton observations of 11 clusters in AO-12 as part of a very large program to complete the sample.

Table 1

XMM-Newton observations that define the complete CHEERS sample.

2.1. Selection of the proposed targets

A substantial sample of suitable clusters is needed to study chemical enrichment in different cluster environments. O, Ne, Mg, Ca and Ni in particular are key elements for constraining the SNIa/SNcc contributions and the SNIa explosion mechanism. The nitrogen abundance, which is sensitive to the IMF, can then be measured in the subsamples that contain cool clusters. We need the RGS to measure the O and N abundance accurately. However, the varying spatial extent and brightness of clusters means that not all clusters are suitable RGS targets because the spatial surface brightness distribution of a cluster determines the spectral line width in the RGS (see Sect. 3.2). The clusters need to be bright and centrally peaked to resolve at least the brightest spectral lines. To select the brightest and the best-suitable clusters for the RGS, we selected this sample mostly from the HIFLUGCS sample (Reiprich & Böhringer 2002).

The SNIa/SNcc contribution ratio is mostly sensitive to the ratio between O and Fe, since they are the best-determined elements. Therefore, we required a statistical significance of about 10σ on O in a single observation. With this criterion, we expect to obtain significances for Ne and Mg of 6σ and 4σ with the RGS, respectively, and ~7σ for N in cool systems (kT ≲ 1 keV). This would in principle be enough to constrain SNcc and AGB models, for which N/Fe, O/Fe, Ne/Fe, and Mg/Fe vary between 30–80% (see Fig. 1). A 10σ signal-to-noise ratio for O is a reasonable requirement for the selection of proposed clusters.

A second criterion is to measure accurate abundances of less abundant elements. A key element is calcium (de Plaa et al. 2007). An accurate Ca abundance guarantees even more precise values for the other elements. The Ca/Fe ratio for different type Ia models varies between 0.33–0.97 times solar (Werner et al. 2008), so an accuracy of 10% solar on the Ca abundance would in principle be sufficient to distinguish between different SNIa models or at least rule out certain models.

Our final criterion for selecting the proposed clusters in AO-12 is an expected uncertainty of 0.036 (10σ) solar for oxygen with the RGS and 0.10 times solar for Ca with EPIC. Clusters were selected to be proposed if this criterion could be reached within an exposure time of ~200 ks. In 85% of all cases, oxygen gives the most stringent selection. The exposure times were increased by 40% to account for possible loss of data due to soft-proton flares. The observations that were performed based on this selection are marked in boldface in Table 1.

2.2. Selection of archival observations

In addition to the proposed targets that were observed by XMM-Newton in AO-12, we also considered archival cluster observations with high-quality RGS data. Sanders et al. (2011), for example, presented a list of high-quality RGS cluster observations. We reprocessed these data, and in contrast to our proposed clusters, selected the clusters for which the oxygen abundance is detected at the 5σ level to ensure that we obtained a reasonably large sample with sufficient spectral quality. Since we obviously do not have control over the exposure time of archival data, we selected the clusters based on the bare minimum statistical quality data that we required. The final list of cluster observations is shown in Table 1.

2.3. Other applications of the selected sample

In this paper, we focus on the O/Fe abundance as measured with the RGS. Recently, the CHEERS sample also yielded science results other than abundance measurements. Some examples are the discovery of cool (~0.2 keV) gas in the CHEERS RGS spectra of elliptical galaxies (Pinto et al. 2014) and constraints on turbulent velocities measured with RGS using line broadening (e.g., Pinto et al. 2015). The RGS data have also shown for NGC 4636 that spatially resolved resonant scattering analysis is capable of revealing velocity structure in the ICM (Ahoranta et al. 2016). This technique will soon be applied to more members of the sample in follow-up papers.

3. Data analysis

We used both the archival and new XMM-Newton exposures listed in Table 1. The observations were processed with the XMM-Newton Science Analysis Software (SAS) version 14.0.0. For each observation, we extracted the event files from the ODF data files using calibration (CCF) files available on 2016/01/31. We used high-resolution spectra from the RGS and data from the MOS1 instrument to extract the spatial line profile used in the spectral fit of the RGS data.

3.1. RGS spectral extraction

We processed the RGS data with the SAS task rgsproc following the standard procedures. In order to decrease the contamination from the soft-proton flares, we extracted RGS light curves from CCD number 9, where hardly any source emission is expected. We binned the light curves in 100 s intervals and fit a Poissonian distribution to the count-rate histogram. We rejected all the time bins for which the number of counts lies outside the interval μ ± 2σ, where μ is the fitted average of the distribution. We used the resulting good time intervals (GTI) files to obtain the filtered event files.

thumbnail Fig. 2

RGS extraction regions and MOS1 stacked image of M 87.

Open with DEXTER

We extracted the RGS source spectra in a region centered on the peak of the source emission with a width of 0.8. We used the model background spectrum created by the standard RGS pipeline, which is a template background file, based on the count rate in CCD 9 of the RGS. In Fig. 2 we show the 0.8 RGS extraction region overlaid on the MOS1 image of M 87.

thumbnail Fig. 3

Example stacked RGS spectrum of NGC 5846. 1T, 2T, and gdem model fits are shown. The colored line labels indicate the most probable origin of the element, e.g., SNIa, SNcc, or AGB stars. The residuals for the three models are shown at the right side.

Open with DEXTER

We also combined the RGS1 and 2 source spectra, the response matrices, and the background files extracted within the 3.4 region (see Fig. 2) through the XMM-SAS task rgscombine. These stacked spectra were only used for plotting purposes. The spectral fits were performed simultaneously on the individual spectra. The stacked RGS spectrum of NGC 5846 is shown in Fig. 3 as an example. We converted the spectra into the SPEX format because we used the SPEX1 spectral fitting package version 3.02.00 for the spectral fitting (Kaastra et al. 1996).

Since we analyzed RGS spectra in units of counts, the errors on the data points are Poisson distributed. Therefore, we minimized the C-statistic (Cash 1979) when we fit models to the RGS spectra.

3.2. RGS spectral broadening

Since RGS is a spectrometer without a slit, the spatial extent of the source causes the measured spectral lines to be broadened (see Davis 2001, for a discussion about grating responses). Photons originating from a region near the cluster center but offset in the direction along the dispersion axis (ΔΘ in arcmin) will be slightly shifted in wavelength (Δλ) with respect to line emission from the cluster center. The wavelength shift is calculated using the following relation: (1)where m is the spectral order (see the XMM-Newton Users Handbook). We corrected for this effect by carefully constructing spatial profiles in corresponding spectral bands from MOS1 data. The MOS1 detector coordinate DETY direction is parallel to the dispersion direction in RGS1 and RGS2, which allows a direct extraction of the surface brightness profile from a MOS1 image. We extracted MOS1 images in detector coordinates for each exposure in the 0.5–1.8 keV (7–25 Å) energy band. For each image, we extracted the surface brightness profile in the dispersion direction through the Rgsvprof task, which is part of the SPEX spectral fitting package. From a MOS1 detector image, this task derives the cumulative spatial profile along the dispersion direction for a certain width in the cross-dispersion and the dispersion direction. For the width in cross-dispersion, we chose widths of 3.4 and 0.8. The width in the dispersion direction is set to 10, since the bulk of the cluster emission is contained within this radius.

The spatial profile from Rgsvprof was convolved with the model spectrum during spectral fitting using the lpro model component in SPEX. The main point of this procedure is to include the line broadening that is due to the spatial extent of the source in the spectral model. It allows us to fit the broadening of the spectral lines and propagate the uncertainty in the spatial broadening into the uncertainties of the other fit parameters, such as the O and Fe abundances. Pinto et al. (2015) showed some examples of spatial profiles in their Fig. 2.

3.3. Spectral modeling

In order to model the multi-temperature structure in clusters (see, e.g., Frank et al. 2013), we used and compared several models available in the SPEX package. In addition to the simple one-temperature (1CIE) and two-temperature (2CIE) models, we also used differential emission-measure (DEM) models. In these models, emission measures are assumed for a range of temperatures on a grid that follow a model or empirical parametrization of the emission measure distribution. The empirical parametrizations are either a truncated power-law distribution (wdem, Sect. 3.3.2) or a Gaussian distribution (gdem, Sect. 3.3.3). For spectral simulations, we also used the classical cooling-flow model. In the models, we fixed the redshift to the most accurate value from optical observations, and we used the Galactic column densities estimated using EPIC (Mernier et al. 2016a), unless stated otherwise. We did not use literature values for the NH, because we found that NH is a source of systematic uncertainty (see Sect. 5.3).

We note that by fitting X-ray spectra, it is difficult to distinguish between the different DEM model parametrizations, like wdem and gdem. Kaastra et al. (2004b) showed that different temperature distributions that share the same emission-weighted average temperature and the same total emission measure produce very similar X-ray spectra that are usually statistically indistinguishable from each other. These DEM models do yield somewhat different abundances when fitted to spectra, therefore multi-temperature structure is a source of systematic uncertainty for abundances that needs to be addressed (see Sect. 5.2).

In all the DEM models, it is implicitly assumed that the abundances in the plasma are the same for all temperature components in the region where the spectrum was extracted. With the current spectral resolution, it is in most cases very hard or even impossible to resolve individual thermal components and to uniquely determine abundances for each temperature. Therefore, we need this assumption to obtain stable fit solutions. This means that the abundances that we measure are essentially emission-weighted average abundances in the fitted region.

3.3.1. 1CIE and 2CIE modeling

All spectra were initially fit using two temperature components (in collisional ionization equilibrium, CIE). The temperatures and emission measures of the two components were left to vary. When one of the components was poorly constrained, a single-temperature or gdem model was chosen. The abundances of both components were coupled to each other in order to be consistent with the DEM models, which assume that all temperature components have the same abundance.

3.3.2. wdem model

One of the differential emission measure models we used is the so-called wdem model, where the emission measure, , of a number of thermal components is distributed as a truncated power law. This is shown in Eq. (2) adapted from Kaastra et al. (2004a): (2)This distribution is cut off at a fraction of Tmax that is βTmax. The value of β was set to 0.1 in this study, which roughly corresponds to the lowest temperatures that are typically detectable with the RGS. The model above is an empirical parametrization of the DEM distribution found in the cores of cool-core clusters (Kaastra et al. 2004a; Sanders et al. 2010). In this form, the limit α → 0 yields the isothermal model at Tmax.

3.3.3. gdem model

Another DEM model that we used is a Gaussian differential emission measure distribution, gdem, in log  T (de Plaa et al. 2006): (3)In this equation, x = log  T and x0 = log  T0, where T0 is the average temperature of the distribution. The width of the Gaussian is σT. Compared to the wdem model, this distribution contains more emission measure at higher temperatures. This model usually yields very similar C-statistic values as the wdem model when fitted to cluster spectra. It resembles the overall shape of the isobaric cooling-flow model, but without the strong emission measure around 0.4 keV, which was not detected with XMM-Newton (see, e.g., Peterson et al. 2001; Tamura et al. 2001).

3.3.4. Cooling-flow model

The isobaric cooling-flow model (see, e.g., Fabian 1994) is the only physical DEM model we used. For this DEM model, the emission measure distribution (dY/ dT) is described by (4)where is the cooling rate in M/yr, k is the Boltzmann constant, μ is the mean molecular weight, and mH is the mass of a hydrogen atom. The Λ(T) stands for the cooling function, which is pre-calculated for an abundance of 0.5 times the solar abundance (Kaastra et al. 2004a). In the model, the dY/ dT is calculated for a grid truncated at a low temperature T1 and a high temperature Tn. The temperature grid typically contains 16 bins.

3.3.5. Updated atomic data and radiation processes

We used thermal plasma models that were developed from the original MEKAL code (Mewe et al. 1985, 1986), with a major update of the Fe-L complex lines by Liedahl et al. (1995). This original code has been included in XSPEC, but has not been updated since then. The MEKAL development continued as an integral part of the SPEX code (Kaastra et al. 1996), where it has been available as the default CIE model. Over the years, it received some updates. Although it is not available separately from SPEX, the model is built up from an atomic database and a set of routines that calculate the emission processes and the resulting model spectrum. The database and the related routines are called SPEXACT2. Between 1996 and 2016, the CIE model in SPEX was updated regularly and was the default SPEX CIE model. We refer to this model as SPEXACT version 2.06.

With the release of SPEX version 3.0 early in 2016, a newly developed spectral emission code became publicly available in the SPEX package. This code contains newly calculated atomic data and more accurate approximations of the emission processes in hot plasmas. For example, the radiative recombination (RR) component of the line emissivity was approximated by a power law in SPEXACT v1 and v2, while the true relation is slightly curved, which causes the oxygen abundance to be biased in certain temperature ranges (see Mernier et al. 2016a). In SPEXACT v3.02, the RR rates are updated and now produce a much more accurate oxygen abundance values. In this paper, we mainly used SPEXACT version 3.02 to fit the spectra. The iron lines, however, were still calculated using SPEXACT 2.06 because the Fe xvii lines are very uncertain in the models (de Plaa et al. 2012). For this paper, we used the calculation by Doron & Behar (2002), which appears to describe the observed Fe xvii line ratios reasonably well.

4. Results

We show the final choice of models that were fit to the spectra in Table 2. For M 87 and Perseus, it was necessary to include a power-law component to account for emission from a central active galactic nucleus (AGN). Most objects clearly needed (at least) two temperatures because we observed lines from Fe xvii and Fe xx. In these cases, a two-temperature fit provides the lowest C-statistics value. For some, mainly cool, objects like M 89, we do not have enough statistics to probe the multi-temperature structure, therefore we chose a single-temperature (1CIE) model. In Abell 85, the two-temperature model provides a slight improvement to a single-temperature model, while the gdem model does not. In the fits, we used the best-fit NH from the EPIC analysis (Mernier et al. 2016a), except for Perseus, which benefits from a free NH value in the fit (see Sect. 5.3).

Table 2

Best-fit (multi-)temperature model for each cluster.

Table 3

Oxygen and iron abundances measured with the RGS.

The final fit results for oxygen and iron obtained from the RGS are listed in Table 3. Since absolute abundances can be more sensitive to systematic effects than relative abundances, we also calculated the O/Fe ratio for comparison. The weighted mean abundances for O and Fe are 0.551 ± 0.010 and 0.556 ± 0.007, respectively. The O and Fe values show considerable scatter. A calculation of the variance yields a value of 0.22 for O and 0.52 for Fe. The scatter in the ratio O/Fe is 0.34.

We can assume that the statistical errors on the measured O and Fe abundances are approximately normally distributed, which means that we can use χ2 statistics when we fit a model to these abundances. When we assume that the parent population of clusters has a constant O/Fe ratio, a fit with a constant value to the abundances yields a χ2 of 102/43 d.o.f., which is formally not acceptable, but much smaller than the χ2 of the individual O and Fe abundances.

thumbnail Fig. 4

O/Fe ratio plotted against the emission weighted temperature determined from RGS. The dark grey line shows the error weighted mean of the sample.

Open with DEXTER

thumbnail Fig. 5

Histogram of the measured O/Fe ratios in the CHEERS sample. The blue line shows a Gaussian fit to the distribution.

Open with DEXTER

The weighted average of the measured O/Fe ratio is 0.853 ± 0.018, which is not the same as the ratio between the mean absolute oxygen and iron abundance. These are just two different estimators, and it is not expected that they yield the same result. Since we used linear abundance ratios, it cannot be expected either that the average O/Fe is the same as Fe/O-1. When we plot the O/Fe ratio as a function of the dominant plasma temperature (see Fig. 4), the points with the smallest error bars appear to cluster around a O/Fe value of 0.8. The points above 1.0 have larger error bars, which is partly due to the error propagation. This means that objects with a lower O/Fe are assigned a slightly higher weight than clusters with a high O/Fe. In a histogram of the O/Fe values (see Fig. 5), a hint of a tail toward higher abundance ratios is visible. When we fit a single Gaussian to the histogram, we find the center of the Gaussian at 0.95 ± 0.04 and a width of 0.19 ± 0.03. The χ2/d.o.f. of the distribution is 9/6, which is formally acceptable. The number of objects in our sample is too low to allow us to statistically distinguish between more complicated shapes of the O/Fe distribution, for example, a bimodal distribution.

5. Fitting biases

The accuracy of the measured abundances need to be studied carefully because several systematic effects can influence the value measured in the spectral fit. First, we study the effect of spatial line broadening on the abundance measurement that is due to the slitless nature of the RGS (see Sect. 5.1). Since the line profile shows the spatial distribution of the emissivity of the ion, each line has in principle a different broadening, which is not corrected for in the spectral fit. Second, the choice of the multi-temperature model fixes the shape of the emission measure distribution, which may not reflect the true distribution in the gas and bias the abundance measurement (see Sect. 5.2). In Sect. 5.3 we study the influence of the assumed NH value on the abundance measurement. Finally, thermal plasma codes also contain uncertainties that affect the abundance measurement (see Sect. 5.4). In Sects. 5.5 and 5.6 we attempt to estimate the effect of the systematic biases on the abundance result.

For each bias effect, we simulate spectra using the RGS response matrices. Poisson noise is not added to the simulated spectra, which means that the value of each data point is the exact mean number of counts expected for that bin in a 100 ks observation. The error on each data point is set to be the square root of the expected number of counts. We checked that the results of the fit are the same as long as the simulated spectra have a minimum quality level (10 000 counts). We do not need many Monte Carlo simulations in this case because we know the expected value for each simulated spectral bin exactly: it is the model value for that bin. We also know the expected variance for each bin, which is the square root of the model value for a given exposure time. If Poisson noise were added, we would need to average out the statistical fluctuations through many Monte Carlo runs to obtain an approximation of the input value, which we already know exactly. Since we wish to compare models with each other, statistical fluctuations are not relevant. The average difference between the best-fit model values and the consequences for the fitted parameters are important, but do not depend on the statistical noise. Only the statistical weight of a data point in the fit minimization matters for the comparison, which is included by the error bars on the simulated data points. We chose this method over a Monte Carlo method with Poisson noise because it saves much unnecessary computation time.

When we compared the SPEX CIE model to APEC version 3.0.1, we used XSPEC to simulate spectra using the APEC model. Since we need to use the same solar abundance table in both SPEX and XSPEC for the comparison and it does not really matter which one, we chose the Lodders (2003) solar photospheric abundances that are available in both packages. In the simulations, solar abundances are assumed for both O and Fe (always with respect to the same solar abundance table as used in the corresponding fits). In the cooling-flow model simulation, the low-temperature cutoff is set at 0.5 keV, and in the wdem models, the cutoff is set at 0.2 keV and α to 2.0.

5.1. Bias that is due to spatial line-broadening

Owing to the nature of the RGS, spectral lines of extended sources are broadened with a width that depends on the spatial emissivity distribution of the respective ion on the sky (see Sect. 3.2 for an explanation of spatial line broadening). This means that essentially two factors determine the width of a line in the RGS: the radial abundance distribution, which sets the amount of elements in the gas as a function of radius, and the radial temperature distribution, which sets the abundance of each ion following the ionization balance at that temperature. These two parameters govern the line emissivity as a function of radius. Temperature and abundance gradients in the cores of clusters therefore can cause lines of different ions to show different line broadening profiles in RGS spectra. Because abundance profiles of different elements are usually not very different from each other (e.g., Mernier et al. 2017), the main cause for differences in line broadening is a steep temperature profile that causes the ionic fractions of ions to vary strongly with radius. This was confirmed by Pinto et al. (2016), who found smaller line widths in cool-temperature components and broader line widths in hotter temperature components in the cores of elliptical galaxies.

In some clusters, the variation in broadening between lines of different ions can be a significant effect. In extreme cases, the O viii line width can be about a factor of 3 larger than the average widths of the iron lines (de Plaa et al. 2006). Since current spectral fitting programs cannot fully model the widths of each line, the spectra are fit with an average width. The question is whether this introduces a systematic bias in the abundance determination.

thumbnail Fig. 6

O, Fe, and O/Fe abundance results for fits to simulated spectra of a range of temperatures. In the simulated spectra, we set the width of the oxygen lines to be twice the width of the Fe lines. In the fit, this difference in width is not fit and is assumed to be the same for all lines. The squares and stars show the absolute O and Fe abundance, respectively. The circles show the O/Fe ratio. All measurements are compared to the input value of 1.

Open with DEXTER

To estimate the effect of line broadening on the measured O/Fe abundance, we simulated RGS spectra for a range of temperatures from 0.6 keV to 6.0 keV. The simulated spectra consisted of the addition of two model spectra with the same temperature and normalization, but with a different broadening profile. For the first model spectrum, we set the O, Ne and Mg abundance to twice solar and the Fe abundance to zero. We broadened this spectrum with a typical spatial profile for a cluster. In the second model spectrum, the O, Ne and Mg abundance were set to zero and Fe was set to twice solar. This model spectrum was broadened with a profile with the same overall shape, but scaled to half the width of the first spectrum. When the two spectra were added, the total spectrum mimicked a spectrum with twice the normalization of the individual components and with the average abundance of both components (i.e., once solar). We chose an average abundance of once solar in the simulations to facilitate recognizing the relative difference between the input and output abundance.

We divided the detectable elements in the RGS into two groups (we assigned O, Ne, and Mg with a broadening twice larger than the broadening of the Fe lines) because the broadening is mainly determined by the strongest lines of Fe and O. The Ne and Mg lines are not that strong and have a lower weight in the determination of the broadening. Given the similar origin of O, Ne, and Mg (core-collapse supernovae) and their comparable atomic weight, we assumed that the line widths of O, Ne, and Mg are very similar to each other and coupled their widths in this simulation to the width of oxygen. The factor of two difference in line width is based on the typical ratio between the line widths of the observed O and Fe lines in different clusters, which in practice varies from 1 (not significantly different) up to a factor of 3 in an extreme case (de Plaa et al. 2006).

We fit the simulated spectrum for each temperature bin with a single-temperature model, but now with a single line-broadening component. The resulting O, Fe, and O/Fe abundance measurements are shown in Fig. 6. We find a bias of about −10% in the O/Fe abundance ratio, which weakly depends on the plasma temperature. Because the fit to the simulated spectra was now constrained to a single broadening profile, it tried to find a weighted-average width between the O and Fe line widths. For O, the model line width was smaller than the simulated width, which means that some of the line flux in the wings of the line was effectively attributed to the continuum. A lower line/continuum ratio leads to an underestimation of the O abundance. For Fe, this works in the opposite way. Because the line model is broader than the simulated line, continuum photons are attributed to the wing of the line, which leads to a higher model line flux and hence to an overestimation of the Fe abundance. Around 1 keV, the Fe line flux is increased through the higher contribution of Fe xvii lines. Therefore, Fe has a relatively large weight in the determination of the line width, and the measured Fe abundance is therefore closer to the input value (at the cost of a larger bias in O). For higher temperatures, the balance between the Fe and O line fluxes changes in favor of O, which means that the bias in O decreases while the bias in Fe increases.

For verification purposes, we also simulated and fit spectra without line broadening. The differences between the input and output values were then smaller than 1%. Therefore, we can fully attribute the biases observed in Fig. 6 to line broadening effects. In reality, the effect of line broadening may be slightly different. In the simulation, the lines are broadened by a convolution, while in RGS observations, the core of the line is as strong as it should be and the wings of the lines are enhanced by the line emission from the regions around the core. Our simulations show the typical magnitude of the systematic differences in abundances that are to be expected, but in a real case, the systematic difference may be in the other direction.

5.2. Bias that is due to multi-temperature structure

It is clear that because of the relatively large field of view of the RGS and the thermally complicated nature of cooling cores, the measured spectra might likely not consist of a single-temperature spectrum. Their true multi-temperature structure is, however, not uniquely constrained. The only physical multi-temperature model that we have is the cooling-flow model (see Sect. 3.3.4), but it was found to describe the first XMM-Newton spectra of clusters not very well (e.g., Peterson et al. 2001; Tamura et al. 2001). Therefore, we used other (empirical) models described in Sect. 3.3 in an attempt to approximate the emission-measure distribution with a power law (wdem) or a Gaussian function (gdem). Since the true distribution is not known and the emissivity of lines depends on the temperature, imperfections in the multi-temperature approximation may bias the measured abundances. We therefore tried to determine biases in abundances for different combinations of multi-temperature models to estimate the typical magnitude of the bias that is due to multi-temperature effects.

thumbnail Fig. 7

Results from two-temperature fits to simulated RGS cooling-flow spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER

In Fig. 7 we show how O and Fe abundance measurements are biased when we simulate an input spectrum using a cooling-flow model and fit it with two single-temperature models. The input cooling-flow model had a low kT limit of 0.5 keV because otherwise we would have created strong O vii and Fe xvii lines that are not observed at this high emissivity (Peterson et al. 2001). At low temperatures, the two-temperature model reproduced the O and Fe abundances well. For higher temperatures, however, the O/Fe was overestimated because the Fe abundance was underestimated.

thumbnail Fig. 8

Results from gdem fits to simulated RGS cooling-flow spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER

When we performed the same experiment, but fit the simulated spectra from the cooling-flow model with a gdem model, we detected a bias in the other direction. Again, the input cooling-flow model had a low kT limit of 0.5 keV. Figure 8 shows the results. The gdem model fits the O/Fe abundance well below 1 keV, but above this temperature, the measured abundances deviate from each other. The resulting O/Fe abundance is underestimated by about 10% in this case.

thumbnail Fig. 9

Results from gdem fits to simulated RGS wdem spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER

In a similar way, we also compared the wdem and gdem models. The input DEM parameters for the wdem simulated spectra were α = 2 and β = 0.2, which are typical observed values for clusters. Figure 9 shows results from the simulated wdem spectra fit with gdem models. In this case, the variation in the results is much larger with temperature. The most significant variations are seen at the low-temperature end, where O and Fe are biased in opposite directions. Around 1 keV, the bias is about 20% in the O/Fe ratio. However, above 2 keV, the bias in the O/Fe drops to a few percent.

For the line-broadening bias (see Sect. 5.1), it is relatively easy to qualitatively explain the observed biases. For the multi-temperature biases that we have estimated in this section, however, it is far more difficult. The reason is that normalization, temperature, and abundance can partly compensate for each other, especially for temperature components that are not dominant. When Fe xvii lines are detected, for example, which indicates the presence of gas with a temperature around 0.50.7 keV, the fit can either try to increase the normalization of the low temperatures in the DEM model, lower the central temperature of the DEM, or increase the Fe abundance. Changing these parameters also affects other bands in the spectrum through the continuum. Therefore, the fit result is the result of a complicated interplay between the assumed temperature distribution and the abundances.

The experiments above show that the bias in the O/Fe ratio and the individual abundances are diverse. We only performed a small subset of tests, which provides a general idea of the accuracy, but not a precise measure. It is therefore difficult to know the bias that is due to multi-temperature structure exactly. Based on this experiment, we can only estimate that the accuracy of the O/Fe abundance ratio is about 10−20%.

Since we do not observe a strong relation between the O/Fe ratio and temperature in the CHEERS sample, it does not appear that we systematically used an inappropriate DEM distribution to fit the spectra. However, we do observe a significant scatter in the measured abundances values. We speculate that the multi-temperature structure varies between objects and that we are unable to model this structure precisely, which leads to random biases in the abundances. Random biases of about 10−20%, as we typically find in our simulations, could be one of the main sources of the scatter that we observe in the O and Fe abundances.

5.3. Bias that is due to uncertainties in the broad-band continuum

Since line strengths are measured with respect to the continuum, uncertainties in the broad-band continuum can significantly affect the measured abundances. The continuum consists of multiple continua from the CIE model, background components, possible AGN power-law emission, and is absorbed by the Galactic interstellar medium (ISM). For most of the objects in the CHEERS sample, AGN power-law emission and background components do not play a major role because we focus on the bright inner cluster cores where the thermal emission in most cases dominates the AGN emission and the background. In this section, we concentrate on the uncertainties that are due to the Galactic NH and discuss the biases that are due to the CIE model in Sect. 5.4.

The absorption column that is due to the ISM is usually quantified using the column density of atomic hydrogen (NH), which is determined using radio surveys (e.g., Kalberla et al. 2005). However, in high-density regions, molecular hydrogen and dust can also significantly contribute to the absorption (Willingale et al. 2013). The X-ray determined NH in clusters of galaxies is not always consistent with the estimates derived from radio and other wavelength bands. This is also partly due to calibration uncertainties in the X-ray instruments and uncertainties in the solar abundance table (Schellenberger et al. 2015). In the EPIC analysis of the CHEERS sample (Mernier et al. 2016a), the X-ray measured NH sometimes deviates from the NH estimated by the tool provided by Willingale et al. (2013), which affects the abundances.

thumbnail Fig. 10

Relative change in O/Fe as a function of difference in NH. Fits with NH values of Willingale et al. (2013) are compared to best-fit NH values using EPIC (Mernier et al. 2016a).

Open with DEXTER

To check the effect of varying NH on the RGS measured abundances, we performed the spectral fitting twice. The first fit assumed the NH using the NH tool by Willingale et al. (2013), and in the second fit, the best-fit NH from EPIC (Mernier et al. 2016a) was assumed. For the Perseus cluster, the NH was left free because the EPIC-determined NH does not fit the RGS data well. In Fig. 10 we show the relative difference between the O/Fe ratio measured using the Willingale et al. (2013)NH values and the ratio measured using the NH determined with EPIC. The plot shows a negative trend with increasing absolute NH. The bias in the O/Fe abundance ratio could increase to ~40% at maximum in rare cases where the absorption column is high (1021 cm-2), like for Perseus. The reason that we see a decreasing trend with increasing differences between the Willingale et al. (2013) and EPIC determined NH values is that a higher NH forces the fit to increase the continuum of the spectral model to fit the data, which affects the line/continuum ratio such that it becomes lower. Since the strongest O line is located at a lower energy than the Fe-L complex, it is more strongly affected by NH variations, and thus we observe a decreasing trend.

Table 4

Results of two gdem fits to the EPIC pn spectrum of the Perseus cluster using different effective areas.

The combination of uncertainty in NH and calibration uncertainty in the soft X-ray band can set a significant bias on the abundance determination. The effective area calibration of the RGS and EPIC is estimated to be accurate within ~5%3 and the systematic uncertainty in the NH values of Willingale et al. (2013) are estimated to be between ~8–16%. As a test of the effect of the calibration uncertainty and the NH value on the O/Fe abundance, we fit the EPIC pn spectrum of Perseus with a gdem model. The Perseus cluster shows the largest deviation in Fig. 10 and has one of the highest temperatures of the sample. Since the effect of the effective area calibration increases with cluster temperature (Schellenberger et al. 2015), the Perseus cluster is a conservative choice for this test. We fit the EPIC spectrum of Perseus twice: first, using the original effective area file of EPIC pn, and second, using a modified effective area that assumes that the Chandra ACIS calibration is correct (modified using the MODARF tool, Schellenberger et al. 2015). In the fits, the NH was left free. The results are shown in Table 4. The change in effective area between EPIC pn and ACIS appears to mainly affect the measured temperature structure. The ACIS temperature is 10% higher than the pn temperature. The NH is hardly affected by the change in effective area, and its value is close to the H i value from radio observations (1.38 × 1021 cm-2, Kalberla et al. 2005). The effect on the O/Fe ratio is modest with a difference of ~8%. From this test, we conclude that uncertainties in the calibration can be compensated for by changing the model parameters, for example, the temperature in this case, which in turn can affect the abundance determination. It is, however, difficult to draw more general conclusions from this test because the compensation effects in the spectral fit can be different for each cluster or instrument.

The test shows that fixing NH to the value of Willingale et al. (2013; 2.12 × 1021 cm-2) would bias the fit substantially, since both the pn and ACIS calibration favor the value of Kalberla et al. (2005). Especially for abundance determinations, it is most important that the continuum is estimated properly. Therefore, it is advisable to fit the NH, instead of fixing it to a literature value, which could limit the freedom of the fit to accommodate for a mismatch between the observed and modeled continuum. In practice, we limited the NH during the fit to the range between the H i value of Kalberla et al. (2005) and the Willingale et al. (2013) value, to avoid the fit to optimize to unphysical values. The range is slightly extended on both the low and high end with 5 × 1019 cm-2 and 1 × 1020 cm-2, respectively, to account for the uncertainties in the NH values (see Mernier et al. 2016a, for details).

5.4. Bias that is due to atomic database and spectral modeling accuracy

The accuracy of spectral models is generally hard to assess because it requires expert knowledge of atomic physics and radiation processes. In X-ray astrophysics, there are two groups actively developing codes to model soft X-ray emission from thermal plasmas in collisional ionization equilibrium (CIE): the group developing APEC/ATOMDB4, and the group developing SPEX. The bases for these codes were originally developed in the 1970s and 1980s, but in recent years, they have been upgraded in preparation for new instruments dedicated to high-resolution spectroscopy. Although both groups rely on atomic data from similar sources, the radiation processes can be complicated, and different assumptions or approximations can result in differences in measured abundances.

thumbnail Fig. 11

Results from one-temperature CIE fits with SPEXACT v2 to simulated APEC spectra for a range of temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER

In this section, we compare the spectrum as calculated by the APEC v3.0.1 code with the default SPEXACT v2 code in SPEX. We also compare APEC to the new SPEXACT v3 code (see Sect. 3.3.5). In this comparison, we simulated spectra for a range of temperatures using APEC (and Lodders 2003, solar photospheric abundances). These spectra were subsequently fit to SPEX models using the same abundance set, and the normalization, temperature, and the abundances of the relevant elements were left free to vary. In Fig. 11 we show the results for the comparison of the APEC spectra with SPEXACT v2 spectra. There appears to be a strong bias in the O/Fe abundance, especially above 1 keV, where the O/Fe ratio is about 50% of the original value. The main origin of this difference is the crude approximation of the radiative recombination process in the original MEKAL model (Mao & Kaastra 2016). Our RGS results for oxygen are corrected for this effect.

In the SPEXACT v3 code, the issue with the radiative recombination is fixed and we can compare the APEC/ATOMDB v3.0.1 and SPEXACT v3.02 to see what the remaining uncertainties are in O/Fe. Figure 12 shows the fit results for the same simulated APEC spectra, but now fit with the latest SPEX model. The strong trend with temperature that was visible in the previous plot for the SPEXACT v2 code has become less pronounced. However, a 10–20% discrepancy between APEC v3.0.1 and SPEX v3.02.00 remains. The origin of these differences are subject of study by the APEC and SPEX teams and include biases that are due to, for example, differences in included atomic data and differences in modeled radiation processes. We expect that only new high-resolution X-ray spectrometers like SXS on board Hitomi (Takahashi et al. 2010) and the X-ray Astronomy Recovery Mission (XARM) or the X-IFU instrument on board the future X-ray observatory ATHENA (Nandra et al. 2013) will be able to provide a sufficient benchmark to improve these spectral models.

Thanks to this study, we were able to find and correct a source of bias in the old thermal models of MEKAL (Mao & Kaastra 2016). This reduced the difference in the O/Fe ratio from 50% to at most 1020%. Differences between the SPEXACT v3 CIE model and APEC/ATOMDB v3.0.1 remain, which can likely be attributed to differences in the amount of spectral lines in the database and differences in the implementation of radiative processes. This case shows that deep observations help to identify and solve systematic biases in the results.

thumbnail Fig. 12

Results from one-temperature CIE fits with SPEXACT v3 to simulated APEC spectra for a range of temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER

5.5. Effect of biases on a simulated sample

The systematic errors on O/Fe that we estimated in the previous sections cannot simply be added because the bias in O/Fe has a highly nonlinear dependence on the observation and cluster parameters. For each cluster, the total bias may turn out to be different. One possible way to estimate the total effect of the biases is to simulate a number of clusters, add as many biases as possible, and measure the effect on O/Fe. Therefore we simulated a set of 50 RGS spectra of 100 ks using SPEXACT v3. In the simulated spectra, the input abundances for O and Fe are solar, but we varied the line widths of O and Fe randomly. For O, we drew a random line width from a distribution with an average of once the spatial broadening profile and with a random Gaussian distributed variance of 0.25. For Fe, the random width is on average 0.75 with the same variance of 0.25. Using these numbers, we obtain random O/Fe line width ratios with an average of ~1.3 and a variance of ~0.4. Most of these random values will fall in the typically observed range of the O/Fe width ratios between 1–2. For each simulated spectrum, we generated a random variety of an asymmetric gdem multi-temperature model. We drew a random number for the sigma parameter for the low-temperature tail with an average of 0.2 and for the high-temperature tail of 0.1, both with a variance of 0.05. Each object was assigned a random redshift between 0.01 and 0.1 and an NH between 1019 and 1022 cm-2. These numbers roughly represent the typical numbers obtained from RGS fits to cluster spectra. The simulated spectra were subsequently fit with a wdem model, a single line width, and a free NH. With this setup, we simulated the effect of bias in the multi-temperature structure, line width, and NH fitting on a sample of clusters. The simulated biases were random with respect to each other and would average out for large numbers of simulated clusters.

thumbnail Fig. 13

O/Fe versus the maximum temperature of a wdem model (kTmax) for a simulated set of 50 RGS cluster spectra. Three O/Fe ratios are not shown, since the O/Fe ratio was not significantly measured.

Open with DEXTER

In Fig. 13 we show the O/Fe ratio versus the kTmax temperature from a wdem model fit to the simulated spectra. The O/Fe ratios do not appear to show a significant trend with temperature, but there is significant intrinsic scatter in the results. The weighted average of O/Fe in the simulated sample is 0.936 ± 0.005, which is significantly lower than the input value of 1. The variance of the O/Fe abundances is 0.14. Although this test does not include all possible biases, it shows that the biases combined are likely to appear as additional systematic scatter on the O/Fe measurements and a bias in the average value, similar to what we observe in the CHEERS sample. It is therefore likely that a considerable fraction of the scatter in O/Fe observed in the CHEERS sample is due to systematic effects, such as inaccurate multi-temperature structure parametrizations or insufficient modeling of line broadening effects.

5.6. An alternative estimation of the total systematic bias

The tests we performed in Sect. 5 basically show the sensitivity of our result, the O/Fe abundance ratio, to different initial assumptions of the spatial broadening, temperature structure, absorption column, and spectral code. This only provides us with a rough indication of what the bias in the O/Fe ratio may be if our initial assumptions are incorrect. The problem is that we do not know which assumptions are right. We try to compare initial assumptions that are generally considered reasonable in the field in order to obtain a typical bias on the O/Fe ratio. As explained in Sect. 5.5, combining these biases into a total bias estimate cannot be calculated exactly. The biases have a direction and magnitude that also depends on parameters such as temperature and absorption column, which makes it hard to attach a general number on it.

However, we also approached the problem from a different direction because we are able to use constraints from our CHEERS sample. Assuming the variation in abundances that we measure are only due to systematic uncertainties and that all objects have the same intrinsic abundances. Then, following the variance in Table 3, the typical systematic error for O and O/Fe would be about 40% and for Fe nearly 100%. Because of the limited bandwidth of the RGS and thus the lack of much line-free continuum, the absolute abundance of Fe is strongly degenerate with respect to the normalization of the CIE component. It is therefore better to consider the O/Fe ratio for now, which is not affected by this degeneracy. When we do this, the maximum systematic uncertainty would typically be 40% on O/Fe. It is unlikely, however, that all local clusters would have exactly the same chemical enrichment history and O/Fe ratio in their cores. The 40% can thus be considered to be an upper limit.

When we consider the typical biases in the O/Fe ratio as shown in our simulations, the differences range between ~10–40%, where the 40% differences only occur in a few extreme cases with large discrepancies in the assumed NH value. When the biases are combined linearly, the combined bias typically ranges between 20−30% if we exclude some exceptional cases. We note that this estimate holds when the latest spectral models (APEC v3.0.1 or SPEXACT v3) are used. The older spectral codes show much larger discrepancies. We therefore recommend using the most recent versions of the spectral codes and being cautious with using literature values for the contribution of molecules to the Galactic absorption.

6. Discussion

Through the analysis of X-ray spectra of a sample of 44 clusters, groups, and elliptical galaxies, we measured the O, Fe, and O/Fe abundance ratios in the cores of these objects. We also estimated the accuracy of our results by testing the sensitivity of the abundance measurements to our initial assumptions in the data analysis. In this section, we discuss the astrophysical implications of the O/Fe distribution as a function of temperature and the observed scatter in the abundance measurements.

Based on the O/Fe ratio alone, we cannot estimate the fraction of SNIa versus SNcc that enriched the ICM because the theoretical yields of supernova models vary too much. We need to combine the measured RGS abundances with abundances measured using EPIC for a more meaningful comparison with supernova models. This comparison is reported in Mernier et al. (2016b).

6.1. O/Fe ratios and the star formation history of the central BCGs

The central metal abundance peaks seen in giant elliptical galaxies, groups, and clusters of galaxies are thought to be dominated by the metal enrichment from the brightest cluster galaxy (BCG) in the cluster center (De Grandi et al. 2004). Most BCGs appear red and dead, indicating that they stopped forming stars approximately 10 billion years ago (Serra & Oosterloo 2010). Today, their enrichment is thought to be dominated by SNIa with a long delay time, and to a lesser extent, by metal-rich stellar winds. While metal uplift by AGN will only redistribute both SNcc and SNIa products locked up in the stars and the gas in the BCG, SNIa explosions are thought to currently dominate the Fe and Ni enrichment. The measured O/Fe ratio is therefore expected to decrease with time since the last significant episode of star formation.

The fact that the O/Fe ratio that we measure does not show a significant trend as a function of temperature could indicate that the star formation history of BCGs does not change as a function of cluster mass, which is consistent with an earlier conclusion by de Grandi & Molendi (2009). The constant O/Fe may mean that most of the metals were already formed before the cluster ICM was assembled around z ≲ 2 (Werner et al. 2013; Simionescu et al. 2015) and the contribution from SNIa with long delay times from the BCG is comparatively small. The constant SNIa/SNcc ratio as a function of cluster radius found by Mernier et al. (2017) suggests that late enrichment by the BCG either produces similar amounts of SNIa and SNcc, or more likely, that this late enrichment does, on average, not contribute much to the ICM abundance. The flat radial abundance profiles in the very core of clusters could be formed by redistribution of “old” enrichment products from the starburst period of the BCG (z ~ 2−4) into the ICM by processes such as AGN feedback and sloshing.

6.2. Intrinsic scatter in the O/Fe abundances

The measured O/Fe ratios are not consistent with a simple constant value, but there is an additional scatter of 0.34 in the O/Fe abundance values. In Sect. 5.5 we showed that systematic effects, such as effects due to uncertainties in the multi-temperature structure or line broadening, contribute to the scatter in the O/Fe abundance. Although these biases are largely due to inaccuracies in the spectral modeling, a part of the scatter may be due to real variations in cluster cores. It is likely that the true multi-temperature structure in the cores of clusters and groups depends on the level of AGN activity. The chemical enrichment history may also show slight variations between objects. In this section, we explore possible astrophysical origins for the remaining fraction of the observed scatter.

In Fig. 5 we find a weak hint of a tail of objects showing O/Fe ratios that are higher than the weighted average, above about 0.9. Although this tail is not formally significant, the distribution may be broadened and skewed due to intrinsic scatter. This may be partly the result of a late star-forming period in some of the BCGs. This result could be consistent with optical observations of BCGs, which indicate that even though most systems stopped forming stars at high rates a long time ago (current star formation rates are on the order of ~0.1 M yr-1), a few systems show star formation rates of a few solar masses per year, very exceptionally up to a few 100 solar masses per year (O’Dea et al. 2008; McDonald et al. 2012). In order to test this hypothesis, we compared our O/Fe measurements to estimates of the age of the stellar population in elliptical galaxies by Serra & Oosterloo (2010). Only eight clusters and groups appear in both samples, which limits the comparison. Interestingly, the single stellar population (SSP) equivalent age (tSSP) from Serra & Oosterloo (2010) and our O/Fe measurements show no correlation. A larger overlapping sample would be necessary to determine whether a significant correlation exists. Thus, the test of our hypothesis is unfortunately inconclusive.

The scatter in the absolute abundances that we measure appears to be larger for iron than for oxygen. Although we should be careful with interpreting absolute abundances with the RGS because of the systematic effects, a comparison between the scatter in oxygen relative to the scatter in iron can be made. The larger scatter in the iron abundance suggests that the object-to-object variation in SNIa enrichment of the core of the ICM exceeds the SNcc variation, either because of variations in the size or age of the stellar populations in the core or variations in the enrichment mechanisms that transport the iron from the BCG into the ICM. Ehlert et al. (2011) showed that in extreme cases, the metallicity peak in the cores of clusters can be mixed into the surrounding gas through violent AGN feedback. Metal uplift by AGN-blown bubbles has been observed in many clusters (Simionescu et al. 2008, 2009; Kirkpatrick et al. 2009, 2011; Kirkpatrick & McNamara 2015). Hydrodynamic simulations by Planelles et al. (2014) show that AGN feedback is likely to play a major role in the radial metal profile. Late-time SNIa enrichment in the BCG may thus be more peaked in one cluster and more extended in the other through AGN feedback. This may explain a part of the scatter in O/Fe that we observe in the CHEERS sample. However, when we plot the O/Fe abundance ratio and the Fe versus the radio luminosity from Bîrzan et al. (2012), we do not see any relation. Since the radio mode activity only shows the recent AGN feedback, an episode of major feedback in the past could have redistributed the metals in the core a long time ago and therefore it does not show a correlation between O/Fe and radio luminosity today. Hence, we cannot draw conclusions from the lack of correlation.

7. Conclusions

We have measured the O/Fe abundances in 44 clusters and groups of galaxies with the RGS and summarize our results below.

  • The O/Fe ratio does not show a significant trend as a function oftemperature in the 0.66 keV range, which suggests that the enrichment of the ICM does not depend on cluster mass and that most of the enrichment took place before the ICM was formed.

  • We estimate that the systematic biases in the O/Fe ratio that are due to the spatial broadening of spectral lines in RGS, the multi-temperature structure, and the spectral models are about 2030%. A thorough analysis of the systematic uncertainties proved to be essential to reach this accuracy.

  • We find a significant scatter in the O/Fe ratio from cluster to cluster. A considerable fraction of the scatter is most likely due to systematic uncertainties originating from the assumed multi-temperature structure, line broadening effects, NH uncertainties, and uncertainties in spectral models. A fraction of the scatter is due to true object-to-object variations in the abundance and multi-temperature structure, for instance. However, using the current data, we are not able to separate these components and estimate the true intrinsic scatter in cluster and group abundances accurately and quantitatively.

  • Thanks to the high statistics of the observations, we have identified and corrected a bias in the old MEKAL code (Mewe et al. 1985) that caused a bias of about 50% in the O/Fe ratio. This shows that deep observations of samples can help reduce the systematic biases in our analyses.

  • Before new missions such as XARM and ATHENA fly, it is essential to invest in improving atomic line databases and spectral models to reduce their biases to a level that makes them comparable to the expected level of statistical precision for these instruments.

For the systematic biases in the determination of abundances with the RGS, we summarize our findings here.

  • The spatial line broadening bias is about 10%. This can be reducedby using multiple lpro models to fit the widths.

  • The multi-temperature structure bias is about 10−20%. Minimizing it requires better observations and modeling of cluster cores.

  • The broadband continuum shape bias is typically ~20%. It can be reduced by fitting the continuum more carefully, for instance, free NH and/or local continuum fitting.

  • The spectral modeling bias is about 10−20%. This can only be reduced by improving atomic data and spectral modeling.


2

SPEX atomic code and tables.

3

See XMM-SOC-CAL-TN-0018 and XMM-SOC-CAL-TN-0030 at https://www.cosmos.esa.int/web/xmm-newton/ calibration-documentation

Acknowledgments

We dedicate this paper to our dear colleague and coauthor Yu-Ying Zhang, who passed away on December 11, 2016. We are deeply grateful for her important scientific contributions and her generous support of the CHEERS project. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. Part of the data is obtained through a successful XMM-Newton VLP proposal by the CHEERS collaboration (proposal ID 072380). The Netherlands Institute for Space Research (SRON) is supported financially by NWO, The Netherlands Organisation for Scientific Research. N.W. has been supported by the Lendület LP2016-11 grant awarded by the Hungarian Academy of Sciences. C.P. and A.C.F. acknowledge funding through ERC Advanced Grant number 340442. P.K. thanks Steve Allen and Ondrej Urban for support and hospitality at Stanford University. Y.Y.Z. acknowledged support by the German BMWi through the Verbundforschung under grant 50 OR 1506. L.L. acknowledges support from the Transregional Collaborative Research Centre TRR33 The Dark Universe and from the HRC contract number NAS8-03060. H.A. acknowledges support from NWO via a Veni grant. G.S. acknowledges support from NASA through XMM-Newton grant NNX15AG25G. T.H.R. and G.S. acknowledge support from the DFG through Heisenberg grant RE 1462/5, grant RE 1462/6, and from TRR 33. The authors thank Liyi Gu for helpful suggestions regarding the manuscript.

References

All Tables

Table 1

XMM-Newton observations that define the complete CHEERS sample.

Table 2

Best-fit (multi-)temperature model for each cluster.

Table 3

Oxygen and iron abundances measured with the RGS.

Table 4

Results of two gdem fits to the EPIC pn spectrum of the Perseus cluster using different effective areas.

All Figures

thumbnail Fig. 1

Expected abundances measured in a typical long XMM-Newton observation of 120 ks (bottom panel). The estimates for the SNIa, SNcc, and AGB contribution are based on a sample of 22 clusters (de Plaa et al. 2007) and two elliptical galaxies (Grange et al. 2011). The top panel shows the typical range in SNIa and IMF models with respect to the statistical error bars in the observation. Figure adapted from de Plaa (2013).

Open with DEXTER
In the text
thumbnail Fig. 2

RGS extraction regions and MOS1 stacked image of M 87.

Open with DEXTER
In the text
thumbnail Fig. 3

Example stacked RGS spectrum of NGC 5846. 1T, 2T, and gdem model fits are shown. The colored line labels indicate the most probable origin of the element, e.g., SNIa, SNcc, or AGB stars. The residuals for the three models are shown at the right side.

Open with DEXTER
In the text
thumbnail Fig. 4

O/Fe ratio plotted against the emission weighted temperature determined from RGS. The dark grey line shows the error weighted mean of the sample.

Open with DEXTER
In the text
thumbnail Fig. 5

Histogram of the measured O/Fe ratios in the CHEERS sample. The blue line shows a Gaussian fit to the distribution.

Open with DEXTER
In the text
thumbnail Fig. 6

O, Fe, and O/Fe abundance results for fits to simulated spectra of a range of temperatures. In the simulated spectra, we set the width of the oxygen lines to be twice the width of the Fe lines. In the fit, this difference in width is not fit and is assumed to be the same for all lines. The squares and stars show the absolute O and Fe abundance, respectively. The circles show the O/Fe ratio. All measurements are compared to the input value of 1.

Open with DEXTER
In the text
thumbnail Fig. 7

Results from two-temperature fits to simulated RGS cooling-flow spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER
In the text
thumbnail Fig. 8

Results from gdem fits to simulated RGS cooling-flow spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER
In the text
thumbnail Fig. 9

Results from gdem fits to simulated RGS wdem spectra for a range of (maximum) temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER
In the text
thumbnail Fig. 10

Relative change in O/Fe as a function of difference in NH. Fits with NH values of Willingale et al. (2013) are compared to best-fit NH values using EPIC (Mernier et al. 2016a).

Open with DEXTER
In the text
thumbnail Fig. 11

Results from one-temperature CIE fits with SPEXACT v2 to simulated APEC spectra for a range of temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER
In the text
thumbnail Fig. 12

Results from one-temperature CIE fits with SPEXACT v3 to simulated APEC spectra for a range of temperatures. The measured O, Fe, and O/Fe abundances are shown and compared to their input value of once solar.

Open with DEXTER
In the text
thumbnail Fig. 13

O/Fe versus the maximum temperature of a wdem model (kTmax) for a simulated set of 50 RGS cluster spectra. Three O/Fe ratios are not shown, since the O/Fe ratio was not significantly measured.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.