© ESO, 2015

1. Introduction

The deep gravitational potential of clusters of galaxies retains large amounts of hot (~10⁷–10⁸ K) gas, mainly visible in X-rays, which accounts for no less than 80% of the total baryonic mass. This so-called intra-cluster medium (ICM) contains not only H and He ions, but also heavier metals. Iron (Fe) was discovered in the ICM with the first generation of X-ray satellites (Mitchell et al. 1976); then neon (Ne), magnesium (Mg), silicon (Si), sulfur (S), argon (Ar), and calcium (Ca) were measured with ASCA (e.g. Mushotzky et al. 1996). Precise abundance measurements of these elements have been made possible thanks to the good spectral resolution and the large effective area of the XMM-Newton (Jansen et al. 2001) instruments (e.g. Tamura et al. 2001). Nickel (Ni) abundance measurements and the detection of rare elements like chromium (Cr) have been reported as well (e.g. Werner et al. 2006; Tamura et al. 2009). Finally, thanks to its low and stable instrumental background, Suzaku is capable of providing accurate abundance measurements in the cluster outskirts (e.g. Werner et al. 2013).

These metals clearly do not have a primordial origin; they are thought to be mostly produced by supernovae (SNe) within cluster galaxy members and have enriched the ICM mainly around z ~ 2−3, i.e. during a peak of the star formation rate (Hopkins & Beacom 2006). However, the respective contributions of the different transport processes required to explain this enrichment are still under debate. Among them, galactic winds (De Young 1978; Baumgartner & Breitschwerdt 2009) are thought to play the most important role in the ICM enrichment itself. Ram-pressure stripping (Gunn & Gott 1972; Schindler et al. 2005), galaxy-galaxy interactions (Gnedin 1998; Kapferer et al. 2005), active galactic nucleus (AGN) outflows (Simionescu et al. 2008, 2009b), and perhaps gas sloshing (Simionescu et al. 2010) can also contribute to the redistribution of elements. Studying the metal distribution in the ICM is a crucial step in order to understand and quantify the role of these mechanisms in the chemical enrichment of clusters.

Another open question is the relative contribution of SNe types producing each chemical element. While O, Ne, and Mg are thought to be produced mainly by core-collapse SNe (SNcc, including types Ib, Ic, and II, e.g. Nomoto et al. 2006), heavier elements like Ar, Ca, Fe, and Ni are probably produced mainly by type Ia SNe (SNIa, e.g. Iwamoto et al. 1999). The elements Si and S are produced by both types (see de Plaa 2013, for a review). The abundances of high-mass elements highly depend on SNIa explosion mechanisms, while the abundances of the low-mass elements (e.g. nitrogen) are sensitive to the stellar initial mass function (IMF). Therefore, measuring accurate abundances in the ICM can help to constrain or even rule out some models and scenarios. Moreover, significant discrepancies exist between recent measurements and expectations from current favoured theoretical yields (e.g. de Plaa et al. 2007), and thus require further investigation.

The temperature distribution in the ICM is often complicated and its underlying physics is not yet fully understood. For instance, many relaxed cluster cores are radiatively cooling on short cosmic timescales, which was presumed to lead to so-called cooling flows (see Fabian 1994, for a review). However, the lack of cool gas (including the associated star formation) in the core revealed in particular by XMM-Newton (Peterson et al. 2001; Tamura et al. 2001; Kaastra et al. 2001) leads to the so-called cooling-flow problem and argues for substantial heating mechanisms, yet to be found and understood. For example, heating by AGN could explain the lack of cool gas (see e.g. Cattaneo & Teyssier 2007). Studying the spatial structure of the ICM temperature in galaxy clusters may help to solve it.

Abell 4059 is a good example of a nearby (z = 0.0460, Reiprich & Böhringer 2002) cool-core cluster. Its central cD galaxy hosts the radio source PKS 2354-35 which exhibits two radio lobes along the galaxy major axis (Taylor et al. 1994). In addition to ASCA and ROSAT observations (Ohashi 1995; Huang & Sarazin 1998), previous Chandra studies (Heinz et al. 2002; Choi et al. 2004; Reynolds et al. 2008) show a ridge of additional X-ray emission located ~20 kpc south-west of the core, as well as two X-ray ghost cavities that only partly coincide with the radio lobes. Moreover, the south-west ridge has been found to be colder, denser, and with a higher metallicity than the rest of the ICM, suggesting a past merging history of the core prior to the triggering of the AGN activity.

In this paper we analyse in detail two deep XMM-Newton observations (~200 ks in total) of A 4059, obtained through the CHEERS¹ project (de Plaa et al., in prep.). The XMM-Newton European Photon Imaging Camera (EPIC) instruments allow us to derive the abundances of O, Ne, Mg, Si, S, Ar, Ca, Fe, and Ni not only in the core, but also up to ~650 kpc in the outer parts of the ICM. The XMM-Newton Reflection Grating Spectrometer (RGS) instruments are also used to measure N, O, Ne, Mg, Si, and Fe. This paper is structured as follows. The data reduction is described in Sect. 2. We discuss our selected spectral models and our background estimation in Sect. 3. We then present our temperature and abundance measurements in the cluster core, as well as their systematic uncertainties (Sect. 4), measured radial profiles (Sect. 5), and temperature and Fe abundance maps (Sect. 6). We discuss and interpret our results in Sect. 7 and conclude in Sect. 8. Throughout this paper we assume H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3, and Ω_Λ = 0.7. At the redshift of 0.0460, 1 arcmin corresponds to ~54 kpc. The whole EPIC field of view (FoV) covers R ≃ 0.81 Mpc ≃0.51r₂₀₀ (Reiprich & Böhringer 2002, where r₂₀₀ is the radius within which the density of cluster reaches 200 times the critical density of the Universe). All the abundances are given relative to the proto-solar values from Lodders et al. (2009). The error bars indicate 1σ uncertainties at a 68% confidence level. Unless mentioned otherwise, all our spectral analyses are done within 0.3–10 keV by using the Cash statistic (Cash 1979).

2. Observations and data reduction

Two deep observations (DO) of A 4059 were taken on 11 and 13 May 2013 with a gross exposure time of 96 ks and 95 ks respectively (here after DO 1 and DO 2). In addition to these deep observations, two shorter observations (SO; see also Zhang et al. 2011) are available from the XMM-Newton archive. The observations are summarised in Table 1. Both DO and SO datasets are used for the RGS analysis while for the EPIC analysis we only use the DO datasets. In fact, the SO observations account for ~20% of the total exposure time, and consequently the signal-to-noise ratio S/N would increase only by $\sqrt{\mathrm{1.20}}\mathrm{\simeq }\mathrm{1.10}$, while the risk of including extra systematic errors and unstable fits due to the EPIC background components (Sect. 3 and Appendix B) is high. The RGS extraction region is small, has a high S/N, and its background modelling is simpler than using EPIC; therefore, combining the DO and SO remains safe.

The datasets are reduced using the XMM-Newton Science Analysis System (SAS) v13 and partly with the SPEX spectral fitting package (Kaastra et al. 1996) v2.04.

2.1. EPIC

In both DO datasets the MOS and pn instruments were operating in Full Frame mode and Extended Full Frame mode respectively. We reduce MOS 1, MOS 2 and pn data using the SAS tasks emproc and epproc. Next, we filter our data to exclude soft-proton (SP) flares by building appropriate Good Time Intervals (GTI) files (Appendix A.1) and we excise visible point sources to keep the ICM emission only (Appendix A.2). We keep the single, double, triple, and quadruple events in MOS (pattern ≤ 12). Owing to problems regarding charge transfer inefficiency for the double events in the pn detector², we keep only single events in pn (pattern=0). We remove out-of-time events from both images and spectra. After the screening process, the EPIC total net exposure time is ~150 ks (i.e. ~80% of the initial observing time). In addition to EPIC MOS 1 CCD3 and CCD6 which are no longer operational, CCD4 shows obvious signs of deterioration, so we discard its events from both datasets as well.

Figures 1 and 2 show an exposure map corrected combined EPIC image of our full filtered dataset (both detectors cover the full EPIC FoV). The peak of the X-ray emission is seen at ~23h 57′0.8′′ RA, −34°45′34′′ Dec.

We extract the EPIC spectra of the cluster core from a circular region centred on the X-ray peak emission and with a radius of 3 arcmin (Fig. 2). Using the same centre we extract the spectra of eight concentric annuli, together covering the FoV within R ≤ 12 arcmin (Fig. 1). The core region corresponds to the four innermost annuli. The RMFs and ARFs are processed using the SAS tasks rmfgen and arfgen, respectively. In order to look at possible substructures in temperature and metallicity, we also create EPIC maps. We divide our EPIC observations in spatial cells using the Weighted Voronoi Tesselations (WVT) adaptive binning algorithm (Diehl & Statler 2006). We restrict the size of our full maps to R ≤ 6 arcmin. The cell sizes are defined in such a way that in every cell S/N = 100. The relative errors of the measured temperature and Fe abundance are then expected to be not higher than ~5% and ~20%, respectively (see Appendix C for more details). Because SAS does not allow RMFs and ARFs to be processed for complex geometry regions, we extract them on 10 × 10 square regions covering together our whole map and we attribute the raw spectra of each cell to the response files of its closest square region. The spectra and response files are converted into SPEX format using the auxiliary program trafo.

2.2. RGS

Reflection Grating Spectrometer data of all four observations are used (see Table 1 and also Pinto et al. 2015, for details). The RGS detector is centred on the cluster core and its dispersion direction extends from the north-east to the south-west. We process RGS data with the SAS task rgsproc. We correct for contamination from SP flares by using the data from CCD9, where hardly any emission from the source is expected. We build the GTI files similarly to the EPIC analysis (Appendix A.1) and we process the data again with rgsproc by filtering the events with these GTI files. The total RGS net exposure time is 208.4 ks. We extract response matrices and RGS spectra for the observations. The final net exposure times are given in Table 1.

We subtract a model background spectrum created by the standard RGS pipeline from the total spectrum. This is a template background file, based on the count rate in CCD9 of RGS.

We combine the RGS 1 and RGS 2 spectra, responses and background files of the four observations through the SAS task rgscombine obtaining one stacked spectrum for spectral order 1 and one for order 2. The two combined spectra are converted to SPEX format through trafo. Based on the MOS 1 image, we correct the RGS spectra for instrumental broadening as described in Appendix A.3. We include 95% of the cross-dispersion direction in the spectrum.

3. Spectral models

The spectral analysis is done using SPEX. Since there is an important offset in the pointing of the two observations, stacking the spectra and the response files of each of them may lead to bias in the fittings. Moreover, the remaining SP component is found to change from one observation to another (see Appendix B). Therefore, the better option is to fit simultaneously the single spectra of every EPIC instrument and observation. This has been done using trafo.

3.1. CIE

We assume that the ICM is in collisional ionisation equilibrium (CIE) and we use the CIE model in our fits (see the SPEX manual³). Our emission models are corrected from the cosmological redshift and are absorbed by the interstellar medium of the Galaxy (for this pointing N_H ≃ 1.26 × 10²⁰ cm^-2 as obtained with the method of Willingale et al. 2013). The free parameters in the fits are the emission measure $\mathit{Y}\mathrm{=}{}^{\mathrm{\int }}\mathit{n}_{\mathrm{e}}{\mathit{n}}_{\mathrm{H}}\mathrm{d}\mathit{V}$, the single-temperature kT, and O, Ne, Mg, Si, S, Ar, Ca, Fe, and Ni abundances. The other abundances with an atomic number Z ≥ 6 are fixed to the Fe value.

3.2. GDEM

Although CIE single-temperature models (i.e. isothermal) fit the X-ray spectra from the ICM reasonably well, previous papers (see e.g. Peterson et al. 2003; Kaastra et al. 2004; Werner et al. 2006; de Plaa et al. 2006; Simionescu et al. 2009b) have shown that employing a distribution of temperatures in the models provides significantly better fits, especially in the cluster cores. The strong temperature gradient in the case of cooling flows and the 2D projection of the supposed spherical geometry of the ICM suggest that using multi-temperature models would be preferable. Apart from the CIE model mentioned above, we also fit a Gaussian differential emission measure (GDEM) model to our spectra. This model assumes that the emission measure Y follows a Gaussian temperature distribution centred on kT_mean and as defined by $\mathit{Y}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{Y}}_{\mathrm{0}}}{{\mathit{\sigma }}_{\mathit{T}}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\frac{{\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{mean}}\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{T}}^{\mathrm{2}}}\right)\mathit{,}$(1)where x = log (kT) and x_mean = log (kT_mean) (see de Plaa et al. 2006). Compared to the CIE model, the additional free parameter from the GDEM model is the width of the Gaussian emission measure profile σ_T. By definition σ_T = 0 is the isothermal case.

3.3. Cluster emission and background modelling

We fit the spectra of the cluster emission with a CIE and a GDEM model successively, except for the EPIC radial profiles and maps, where only a GDEM model is considered.

Since the EPIC cameras are highly sensitive to the particle background, a precise estimate of the local background is crucial in order to estimate ICM parameters beyond the core (i.e. where this background is comparable to the cluster emission). The emission of A 4059 entirely fills the EPIC FoV, making a direct measure of the local background impossible. Some efforts have been made in the past to deal with this problem (see e.g. Zhang et al. 2009, 2011; Snowden & Kuntz 2013), but a customised procedure based on full modelling is more convenient in our case. In fact, an incorrect subtraction of instrumental fluorescence lines might lead to incorrect abundance estimates.

For each extraction region, several background components are modelled in the EPIC spectra in addition to the cluster emission. This modelling procedure and its application to our extracted regions are fully described in Appendix B. We note that we do not explicitly model the cosmic X-ray background in RGS (although we did in EPIC) because any diffuse emission feature would be smeared out into a broad continuum-like component.

4. Cluster core

4.1. EPIC

Our deep exposure time allows us to get precise abundance measurements in the core, even using EPIC (Fig. 3 left). Moreover, the background is very limited since the cluster emission clearly dominates. Table 2 shows our results, both for the combined fits (MOS+pn) and independent fits (either MOS or pn only).

Using a multi-temperature model clearly improves the combined MOS+pn fit. Nevertheless, even by using a GDEM model, the reduced C-stat value is still high because the excellent statistics of our data reveal anti-correlated residuals between MOS and pn, especially below ~1 keV (Fig. 3, right).

When we fit the EPIC instruments independently, the reduced C-stat number decreases from 1.87 to 1.40 and 1.78 in the MOS and pn fits, respectively. Visually, the models reproduce the spectra better as well. We also note that the temperature and abundances measurements in the core are different between the instruments (Table 2). While temperature discrepancies between MOS and pn have been already reported and investigated (Schellenberger et al. 2015), here we focus on the MOS-pn abundance discrepancies. Figure 4 (left) illustrates these values and shows the absolute abundance measurements obtained from our GDEM models. Except for Ne, Ar, and Ca (all consistent within 2σ), we observe systematically higher values in MOS than in pn. Assuming (for convenience) that the systematic errors are roughly in a Gaussian distribution, we can estimate them for different abundance measurements Z_MOS and Z_pn, having their respective statistical errors σ_MOS and σ_pn, ${\mathit{\sigma }}_{\mathrm{sys}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathrm{tot}}^{\mathrm{2}}\mathrm{-}\frac{{\mathit{\sigma }}_{\mathrm{MOS}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{pn}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{,}$(2)where ${\mathit{\sigma }}_{\mathrm{tot}}\mathrm{=}\sqrt{\mathrm{\left(}\mathrm{\right(}{\mathit{Z}}_{\mathrm{MOS}}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{Z}}_{\mathrm{pn}}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{2}}$ and μ = (Z_MOS + Z_pn)/2. We obtain absolute O, Si, S, and Fe systematic errors of ±25%, ±30%, ±34%, and ±14% respectively. The MOS-pn discrepancies in Mg and Ni are too big to be estimated as reasonable systematic errors (Fig. 4). No systematic errors are necessary for the absolute abundances of Ne, Ar, and Ca.

If we normalise the abundances relative to Fe in each instrument (Fig. 4, right panel), O/Fe is consistent within 2σ and Si/Fe and S/Fe within 3σ. Inversely, the discrepancies on Ar/Fe measurements slightly increase, but their statistical uncertainties are quite large because the main line (~3.1 keV) is weak. We note that the discrepancies in Mg and Ni measurements remain huge and almost unchanged. Based on the same method as above, we find that systematic errors of O/Fe, Si/Fe, and S/Fe are reduced to ±8%, ±15%, and ±20% while the systematic errors of Ar/Fe increase to ±27%.

4.1.1. Equivalent widths

One way of determining the origin of the discrepancies in the fitted abundance from different instruments is to derive the abundances using a more robust approach. Instead of fitting the abundances using the GDEM model directly, we model each main emission line/complex by a Gaussian and a local continuum (hereafter the Gauss method). The GDEM model is still used to fit the local continuum; however, only the Fe abundance is kept to its best-fit value and the other abundances are set to zero⁴. We then check the consistency of this method by comparing it with the abundances reported above (hereafter the GDEM method) in terms of equivalent width (EW), which we define for each line as $\mathit{EW}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{line}}}{{\mathit{F}}_{\mathrm{c}}\mathrm{\left(}\mathit{E}\mathrm{\right)}}\mathit{,}$(3)where F_line and F_c(E) are the fluxes of the line and the continuum at the line energy E, respectively. Since the EW of a line is proportional to the abundance of its ion, in principle both methods should yield the same abundance. We compare them on the strongest lines of Mg, Si, S, Ca, Fe, and Ni in MOS and pn spectra (Table 3) and we convert the average MOS+pn EWs into abundance measurements (Fig. 4). While we find consistency between the Gauss and GDEM methods for Ca and Fe-K lines both in MOS and pn, the other elements need to be further discussed.

The EW of Mg obtained in pn using the Gauss method is, significantly, ~9 times higher than when using the GDEM method. In the latter case, the pn continuum of the model is largely overestimated around ~1.5 keV, making the Mg abundance underestimated. The elements Si and S also show significantly larger EWs in pn using the Gauss method. In terms of abundance measurements, they both agree with the MOS measurements (Fig. 4). We also note that beyond ~1.5 keV the MOS residuals ratio are known to be significantly higher than the pn ones (Read et al. 2014), and peak near the Si line. This might also partly explain the discrepancies found for S, Si, and maybe Mg.

When using the GDEM method for pn, the Ni-K line is poorly fitted. The large difference in EW obtained when fitting it using the Gauss method emphasises this effect. In fact, when fitting the pn spectra using a CIE or GDEM model, a low Ni abundance is computed by the model to compensate the issues in the calibration of the effective area around 1.0–1.5 keV (i.e. where most Ni-L lines are present). For this reason and because of large error bars for the Ni-K line, the fit in pn ignores it.

If we fit the spectra only between 2–10 keV, after freezing kT, σ_T, O, Mg, and Si obtained in our previous fits, we obtain Ni abundances of 1.61 ± 0.35 and 1.37 ± 0.26 for MOS and pn, respectively, making them consistent between each other. This clearly favours the Ni abundance measured with MOS in our previous fits. Interestingly, we also measure Fe abundances of 0.752 ± 0.019 and 0.676 ± 0.017 for MOS and pn, respectively; their discrepancies are then reduced, but still remain. Finally, we note that the pn data are shifted by ~ − 20 eV compared to the model around the Fe-K line; this shift does not affect the abundance measurements though.

Our results on the abundance analysis in the core are summarised in Table 4 and Fig. 4 and are briefly discussed in Sect. 7.1. Because their uncertainties are too large, we choose not to consider Mg and Ni abundances in the rest of the paper. Moreover, although the MOS-pn discrepancies are sometimes large and make some absolute abundance measurements quite uncertain, in the following sections we are more interested in their spatial variations. By comparing combined MOS + pn measurements only, the systematic errors we have shown here should not play an important role in this purpose.

4.2. RGS

Our RGS analysis of the core region focuses on the 7–28 Å (0.44–1.77 keV) first and second order spectra of the RGS detector; RGS stacked spectra are binned by a factor of 5. We test single-, two-temperature CIE models, and a GDEM model for comparison.

The models are redshifted and, to model the absorption, multiplied by a hot model (i.e. an absorption model where the gas is assumed to be in CIE) with a total N_H = 1.26 × 10²⁰ cm^-2 (Willingale et al. 2013), kT = 0.5 eV, and proto-solar abundances.

In order to take into account the emission-line broadening due to the spatial extent of the source, we have convolved the emission components by the lpro multiplicative model in SPEX (Tamura et al. 2004; Pinto et al. 2015).

The RGS order 1 and 2 stacked spectra have been fitted simultaneously (Fig. 5) and the results of the spectral fits are shown in Table 5 and Fig. 4. The 2-T CIE and GDEM fits are comparable in terms of Cash statistics and the models are visually similar. Although there might be some residual emission at temperature below 1 keV that can be reproduced by the 2-T CIE model (Frank et al. 2013), using a GDEM model is more realistic regarding temperature distribution found in the core of most clusters. The abundances are in agreement between the different models because they depend on the relative strength of the lines.

5. EPIC radial profiles

We fit the EPIC spectra from each of the eight annular regions mentioned in Sect. 2 using a GDEM model. We derive projected radial profiles of the temperature, temperature broadening, and abundances (Table 6). In our measurements, all the cluster parameters (Y, kT, σ_T, and abundances) are coupled between the three instruments and the two datasets. Since we ignore the channels below 0.4 keV (MOS) and 0.6 keV (pn) in the outermost annulus to avoid background contamination (Appendix B), we restrict our O radial profile within 9′. For the same reason, the O abundance measurement between 6′–9′ might be biased up to ~25% (i.e. our presumed MOS-pn systematic uncertainty for the O measurement).

In order to quantify the trends that appear in our profiles, we fit them with simple empirical distributions. For temperature and abundance profiles, $\begin{array}{ccc}& & \mathit{kT}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{\infty }}\mathrm{+}\mathit{\alpha }\exp \mathrm{(}\mathrm{-}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}\mathrm{)}\\ & & \mathit{Z}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{\infty }}\mathrm{+}\mathit{\alpha }\exp \mathrm{(}\mathrm{-}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}\mathrm{)}\end{array}$and for σ_T radial profile, ${\mathit{\sigma }}_{\mathit{T}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{\infty }}\mathrm{+}\mathit{\alpha }{\mathit{r}}^{\mathit{\gamma }}\mathit{.}$(6)Table 7 shows the results of our fitted trends. Figure 6 shows the radial profiles and their respective best-fit distributions.

The temperature profile reveals a significant drop from ~2.5′ to the innermost annuli, confirming the presence of a cool-core. Beyond ~2.5′, the temperature stabilises around kT ~ 4.2 keV. More surprisingly, after a plateau around 0.22 from the core to ~2.5′, σ_T increases up to 0.33 ± 0.04 in the outermost annulus. This increase is significant in our best-fit distribution. In this outer region, we show that kT and σ_T are slightly correlated (Fig. 7); however, the radial profiles of kT and σ_T show different trends. Moreover, constraining σ_T = 0 in the outermost annulus clearly deteriorates the goodness of the fit (Fig. 7), meaning that the σ_T increase is probably genuine.

Our analysis reveals a slightly decreasing O radial profile. Even if fully excluding a flat trend is hard based on our data, the exponential model (Eq. (5)) gives a better fit than a constant model Z(r) = D_∞ (Table 7). A decrease from 0.54 ± 0.06 to 0.29 ± 0.06 is observed between 0.5′–6′ as well. Finally, O is still strongly detected in the outermost annuli. We note, however, that additional uncertainties should be taken into account (see above). In fact, the O measurement near the edge of the FoV may also be slightly affected by the modelling of the Local Hot Bubble (Appendix B) through its flux and its assumed O abundance.

As mentioned earlier, Ne is hard to constrain, but is detected. Its abundance drops to zero outside the core while it is found to be more than half its proto-solar value within 0.5 arcmin. Profiles of Si and S abundances also decrease, typically from ~0.8 to very low values in the outermost annuli. In every annulus the Si and S measurements are quite similar; this is also confirmed by the best-fit trends which exhibit consistent parameters between the two profiles. The Ar radial profile is harder to interpret because of its large uncertainties, but the trend suggests the same decreasing profile as observed for Si and S.

The Ca radial profile shows particularly high abundances in general, significantly peaked toward the core where it reaches 1.8 ± 0.3 times the proto-solar value and 2.0 ± 0.3 times the local Fe abundance. Finally, we show that Fe abundance is also significantly peaked within the core and decreases toward the outskirts, where our fitted model suggests a flattening to 0.14 ± 0.03.

We note that our radial analysis focuses on the projected profiles only. Although deprojection can give a rough idea about the 3D behaviour of the radial profiles, they are based on the assumption of a spherical symmetry, which is far from being the case in the innermost parts of A 4059 (Sect. 6). Moreover, the deprojected abundance radial profiles are not thought to deviate significantly from the projected ones (see e.g. Werner et al. 2006). Based on the analysis of Kaastra et al. (2004), we estimate that the contamination of photons into incorrect annuli as a result of the EPIC point-spread function (PSF) changes our Fe abundance measurements by ~2% and ~4% in the first and second innermost annuli, respectively, which is not significant regarding our 1σ error bars. The choice of a GDEM model should take into account both the multi-temperature features due to projection effects and the possible PSF contamination in the kT radial profile.

6. Temperature, σ_T, and Fe abundance maps

Using a GDEM model, we derive temperature and abundance maps from the EPIC data of our two deep observations. The long net exposure time (~140 ks) for A 4059 allows the distribution of kT, σ_T, and Fe abundance to be mapped within 6′. As in the radial analysis, all the EPIC instruments and the two datasets are fitted simultaneously.

In order to emphasise the impact of the statistical errors on the maps and to possibly reveal substructures, we create so-called residuals maps following the method of Lovisari et al. (2011). In each cell, we subtract from each measured parameter the respective value estimated from our modelled radial profile (Fig. 6) at the distance r of the geometric centre of the cell. The significance index is defined as being this difference divided by the error on the measured parameter. The kT, σ_T, and Fe abundance maps and their respective error and residuals maps are shown in Fig. 8.

The kT map reveals the cool core of the cluster in detail. It appears to be asymmetric and to have a roughly conic shape extending from the north to the east and pointing toward the south-west. Along this axis, the temperature gradient is steeper to the south-west than to the north-east of the core. Most of the relative errors obtained with the CIE model (not shown here) are within 2–5%, which is in agreement with our expectations (Appendix C); however, they slightly increase with radius. This trend is stronger when using the GDEM model, and the errors are somewhat larger. A very local part (~5 cells) of the core is up to 8σ cooler than our modelled temperature profile. This coldest part is offset ~25′′ SW from the X-ray peak emission. This contrasts with the western part of the core, which shows a significantly hotter bow than the average ~55′′ away from the X-ray peak emission. We also note that some outer cells are found significantly (>2σ) colder or hotter than the radial trend.

The σ_T map confirms the positive σ_T measurements in most of the cells outside the core, typically within 0.1–0.4. Globally, σ_T is consistent with that measured from the σ_T radial profile. We note that outside the core the errors are inhomogeneous and are sometimes hard to estimate precisely.

The Fe map also shows that the core is asymmetric. As it is in the kT map, the abundance gradient from the core toward the south-west is steeper than toward the north-east. The highest Fe emitting region is found to be ~25′′ SW offset from the X-ray peak emission and coincides with the coldest region. In this offset SW region, Fe is measured to be more than 7σ over-abundant.

We note that the smallest cells (~12′′) have a size comparable to the EPIC PSF (~6′′ FWHM); a contamination from leaking photons between adjacent cells might thus slightly affect our mapping analysis. However, the PSF has a smoothing effect on the spatial distributions, and gradients may be only stronger than they actually show in the map. This does not affect our conclusion of important asymmetries of temperature and Fe abundance in the core of A 4059.

7. Discussion

We determined the temperature distribution and the elemental abundances of O, Ne, Si, S, Ar, Ca, and Fe in the core region (≤3′) of A 4059 and in eight concentric annuli centred on the core. In addition, we built 2D maps of the mean temperature (kT), the temperature broadening (σ_T), and the Fe abundance. Because of the large cross-calibration uncertainties, Mg and Ni abundances are not reliable in these datasets using EPIC, and we prefer to measure the Mg abundance using RGS instead.

7.1. Abundance uncertainties and SNe yields

As shown in Table 2, the Ne abundance measured using EPIC depends strongly on the choice of the modelled temperature distribution. The main Ne lines are hidden in the Fe-L complex, around ~1 keV. This complex contains many strong Fe lines and is extremely sensitive to temperature. A slight change in the temperature distribution will thus significantly affect the Ne abundance measurement, making it not very reliable using EPIC (see also Werner et al. 2006). For the same reason, Fe abundances of single- and multi-temperature models might change slightly but already cause a significant difference between both models.

Most of the discrepancies in the abundance determination between the EPIC instruments come from an incorrect estimation of the lines and/or the continuum in pn (Sect. 4.1). Cross-calibration issues between MOS and pn have been already reported (see e.g. de Plaa et al. 2007; Schellenberger et al. 2015), but their deterioration has probably increased over time despite current calibration efforts (Read et al. 2014). Our analysis using the Gauss method (Table 3 and Fig. 4) suggests that in general MOS is more reliable than pn in our case, even though MOS might slightly overestimate some elements as well (e.g. Mg, S, or even Fe). In all cases, this latest method is the most robust one with which to estimate the abundances in the core using EPIC.

Another interesting result is our detection of very high Ca/Fe abundances in the core. This trend has been already reported by de Plaa et al. (2006) in Sérsic 159-03 (see also de Plaa et al. 2007). Within 0.5′ the combined EPIC measurements give a Ca/Fe ratio of 2.0 ± 0.3. This is even higher than measured within 3′ (Ca/Fe = 1.45 ± 0.14). Following the approach of de Plaa et al. (2007) and assuming a Salpeter IMF, we select different SNIa models (soft deflagration versus delayed-detonation, Iwamoto et al. 1999) as well as different initial metallicities affecting the yields from SNcc population (Nomoto et al. 2006). We fit the constructed SNe models to our measured abundances in the core (O, Ne, Mg, and Si from RGS; Ar and Ca from EPIC; Fe from the Gauss method). We find that a WDD2 model, taken with Z = 0.02 and a Salpeter IMF, reproduce our measurements best, with (χ²/ d.o.f.)_WDD2 = 4.28/6 (Fig. 9). Although the fit is reasonable in terms of reduced χ², it is unable to explain the high Ca/Fe value that we found. Based again on de Plaa et al. (2007), we also considered a delayed-detonation model that fitted the Tycho SNIa remnant best (Badenes et al. 2006). The fit is improved ((χ²/ d.o.f.)_Tycho = 1.77/5), but the model barely reaches the lower error bar of our measured Ca/Fe. Assuming that the problem is not fully solved even by using the latest model, we can raise two further hypotheses that might explain it:

• 1.

  Calcium abundance measurements might suffer fromadditional systematic uncertainties. Our analysis(Sects. 4.1 and 5) shows, however,that MOS and pn Ca/Fe measurements are consistent within theentire core (3′). Moreover, the continuum and EW of Ca lines (~3.9 keV) are correctly estimated by our CIE models. Because of current efforts to limit them, uncertainties in the atomic database can contribute only partly. Finally the effective area at the position of this line is smooth and no instrumental-line feature is known around ~3.9 keV.

• 2.

  Some SNe subclasses, so far ignored, might contribute to the metal enrichment in the ICM. For example, the so-called calcium-rich gap transients as a possible subclass of SNIa, are expected to produce a large amount of Ca even outside galaxies, making the transportation of Ca in the ICM much easier (Mulchaey et al. 2014).

7.2. Abundance radial profiles

All the abundance radial profiles decrease with radius. Interestingly, O shows a slight decrease (confirmed by our empirical fitted distribution), even though a flat profile cannot be fully excluded. This decreasing trend has been observed in other clusters, such as Hydra A (Simionescu et al. 2009a), A2029, and Centaurus (Lovisari et al. 2011). However, the observations of A 496 (Lovisari et al. 2011) and A 1060 (Sato et al. 2007) suggest a flatter profile. The O distribution is less clear in Sérsic 159-03 (de Plaa et al. 2006; Lovisari et al. 2011).

Moreover, only O and Fe profiles show abundances significantly higher than zero in the outermost annuli. The Fe profile is clearly peaked to the core, and agrees with typical slopes found in many other clusters (e.g. Simionescu et al. 2009a; Lovisari et al. 2011). Moreover, its apparent plateau in the outer regions may suggest a constant Fe abundance in the ICM even outside r₅₀₀, as recently observed by Suzaku in Perseus (Werner et al. 2013) and other clusters (e.g. Leccardi & Molendi 2008; Matsushita 2011). As seen in Fig. 6, the Fe abundance found in the outskirts of Perseus (0.303 ± 0.012, in proto-solar abundance units) is higher than what we find for A 4059, even when accounting for the systematic uncertainties estimated from the core in Sect. 4.1. This constant Fe abundance found in other cluster outskirts and this work suggest that the bulk of the enrichment at least by SNIa started in the early stages of the cluster formation.

In the previous cluster analyses where O appeared to be flat, the increase of O/Fe with radius is usually justified by arguing a very early population of SNIa and SNcc, starting after an intense star formation around z ~ 2–3 (Hopkins & Beacom 2006) and undergoing a very efficient mixing all over the potential well, followed by a delayed population of SNIa responsible for the Fe peaked profile, and produced preferably in the central galaxy members in which a strong ram-pressure stripping is assumed (see also discussion for Sérsic 159-03 from de Plaa et al. 2006). It has also been suggested that ram-pressure stripping could shape the Fe peak profile between z = 1 and z = 0 (Schindler et al. 2005). However, De Grandi et al. (2014) suggest that the bulk of the Fe peak was already in place before z = 1 in most clusters, meaning that at least SNIa type products started to get a centrally peaked distribution early on in the cluster formation. In fact, Fe seems to follow the near-infrared light profile of the central cD galaxies much better at z = 1 than at z = 0, suggesting that most of the current mixing mechanisms tend to spread out the metals in the ICM.

The decreasing O radial profile measured in this work suggests that the same kind of scenario is likely for SNcc type products. Although its best-fit slope of the profile appears to be flatter than the slope of the Fe radial profile (Table 7), the O/Fe radial values are still compatible with a constant distribution (except possibly for the 6′–9′ annulus, where systematics might affect the O measurements). Consequently, it is not necessary to invoke a delayed population of SNIa and/or SNcc occurring after z = 1, although it might contribute to a minor part of the metals found in the core. At z ~ 2–3 the central cD galaxy and its surrounding galaxy members were already actively star-forming and could have produced the bulk of all metals observed in the core, probably injected into the ICM through galactic winds. More recently, ram-pressure stripping could have also played a minor role in the enrichment of the core, for example to explain the asymmetry found on the maps (see below).

Assuming a flat and positive distribution of Fe and O beyond the FoV, the mixing of the metals is likely very efficient in the outskirts, where the entropy is high. In the core however, the entropy was already very stratified early on without any major mergers to disturb it, and the mixing mechanisms could be less efficient there.

While O and Fe are detected far from the core and this favours an early initial enrichment from SNIa and SNcc types, puzzlingly we do not detect significant abundances of Ne and Si in the outermost annuli. This result is less striking in the S and Ar radial measurements, even though our fitted trends give small upper limits for D_∞. Nevertheless, abundance measurements in the outer parts of the FoV can also suffer from additional systematic uncertainties related to the background contribution. These uncertainties may explain our lack of clear detection of Ne, Si, S, and Ar in the outermost annuli. Finally, we note the similarity between the Si and S profiles, already reported in the cD galaxy M 87 by Million et al. (2011).

In addition to these radial trends, our maps show local regions of anomalously rich Fe abundance in the core. This is particularly striking in the south-west ridge, where the Fe abundance is >7σ higher than the average trend from its corresponding radial profile. Since no galaxy can be associated with this particular region, it is hard to explain its enrichment with galactic winds. As previously reported and discussed by Reynolds et al. (2008), it is possible that an important part of the metals in the core comes from one early starburst galaxy that passed very close to the cD central galaxy before the onset of the central AGN. In this case ram-pressure stripping could probably have played a dominant role in the enrichment within ~0.5 arcmin after the initial enrichment seen through the radial profiles. This possible scenario is also discussed in the next section.

7.3. Temperature structures and asymmetries

Although the ICM appears homogeneous and symmetric at large scale, the inner part appears to be more asymmetric (Fig. 2). As already observed in the past by Chandra (Heinz et al. 2002; Reynolds et al. 2008), the south-west ridge is clearly visible as an additional peaked X-ray emission near the core, and a diffuse tail from the core toward the north-east can also be detected.

Evidence of asymmetries are also found in our spectral analyses. Although our radial kT profile looks similar to other cool-core clusters, our kT and Fe abundance maps show clear inhomogeneities in the ICM structure of A 4059. Compared to the 2D maps previously measured using Chandra (Reynolds et al. 2008), the S/N of the cells in our EPIC maps are ~3.3 and ~2.5 times greater for kT and the Fe abundance, respectively, allowing us to confirm these substructures with a higher precision and over a larger FoV.

First, like the Fe abundance, the temperature gradient is steeper within the south-west ridge than north-east of the core. The central core (including the south-west ridge) is also significantly colder (~2.3 keV) and the south-west ridge has a higher Fe abundance (~1.5) than the rest of the core within 0.5′. These results confirm the previous study by Reynolds et al. (2008) who also found strong asymmetry in the core of A 4059 using Chandra. Their pressure map shows neither asymmetry nor discontinuity in the core, even around the south-west ridge. From both Chandra and XMM-Newton studies, it is clear that this ridge plays a role in the metal enrichment of the core (see also Sect. 7.2) and must be closely linked to the history of the cluster (Reynolds et al. 2008). The hotter bow region found W of the core is likely related to it. Based on the Chandra images (Reynolds et al. 2008), sloshing seems an unlikely explaination for the origin of the ridge. Indeed, it appears to be a second brightness peak separated from the core, and its particular morphology is very different from the typical spiral regular pattern of sloshing fronts (see e.g. Paterno-Mahler et al. 2013; Ichinohe et al. 2014). Another scenario is that this local cool, dense, and Fe-rich asymmetry was already present before the triggering of central AGN radio-activity; it was formed by a gas-rich late-type galaxy that plunged very close to the central cD galaxy. An intense starburst caused by its interactions with the dense local ICM occurred and it lost an important part of its metals as a result of the strong gravitational interaction coupled with intense ram-pressure stripping.

Reynolds et al. (2008) estimated that such a galaxy should be within 300v₃ kpc of the cluster core. They suggested the bright spiral galaxy ESO 349-G009 as being a good candidate, although they were not sure whether this object belongs to A 4059. Looking at the caustic taken from Zhang et al. (2011, see individual galaxy redshifts in the references therein, e.g. Andernach et al. 2005, we can confirm that this is indeed the case (Fig. 10). The galaxy is located in the front part of the cluster and moves with a high radial velocity compared to the cD galaxy (Δv ≃ 1800 km s^-1). Assuming that this scenario is correct and that the movement of this galaxy near the central cD galaxy was essentially along the line of sight, the absence of an obvious metal tail from ram-pressure stripping on the plane of the FoV is naturally explained. Moreover, the X-ray isophotes joining ESO 349-G009 and the cluster ICM (Fig. 11) show an interaction between them and might suggest that the galaxy is escaping from the core. The UV light detected in its arms using the XMM-Newton OM instrument (e.g. UVM2 filter) reveals that the galaxy still has a high star formation rate. The gas mass of the ridge (5 × 10⁹M_⊙) is a small percentage of the total stellar mass of ESO 349-G009 (Reynolds et al. 2008).

Finally, both the radial profile and map reveal a constant or increasing trend of σ_T with radius. This is likely explained by projection effects such as the increased effective length along our line of sight. For cooling core clusters, this effective length increases as a function of radius, and a longer effective length will mix more temperatures along the line of sight. A still broad range of temperatures in the local ICM beyond the core cannot be fully excluded, but seems more unlikely. Indeed, although the few outer local colder or hotter cells found in the kT residuals map (Fig. 8, top-right) might argue in favour of this second explanation, the temperature (and thus σ_T) measurements in the outer map cells are very sensitive to the background modelling, and are thus affected by these additional systematic uncertainties.

8. Conclusions

In this paper we have studied a very deep XMM-Newton observation (~140 ks of net exposure time) of the nearby cool-core cluster A 4059. Several temperature and abundance parameters have been derived from the spectra both in the core and in eight concentric annuli; moreover, we were able to derive kT, σ_T, and Fe abundance maps. We conclude the following:

• The temperature structure shows the cool-core and in additionincreasing deviations from apparent isothermality in and out ofthe core.

• The abundances of O, Ne, Si, S, Ar, Ca, and Fe all are peaked toward the core, and we report the presence of Fe and O beyond ~0.3 Mpc from the core. This suggests that the enrichment from SNIa and SNcc started early on in the cluster formation, probably through galactic winds in the young galaxy members.

• The EPIC image as well as the temperature and Fe abundance maps reveal strong asymmetries in the cluster core. We confirm a colder and Fe-richer ridge south-west of the core, previously found by Chandra, perhaps due to an intense ram-pressure stripping event. Therefore, in addition to an early enrichment through galactic winds, ram-pressure stripping might have greatly contributed to a more recent enrichment of the inner core.

• The Ca/Fe abundance ratio in the core is particularly high (1.45 ± 0.14 using a combined EPIC fit), even accounting for systematic uncertainties. If we assume the Ca/Fe abundance of the entire core to be genuine, it is unlikely explained by current standard SNe yields models. Recently proposed calcium-rich gap transient SNIa might be an interesting alternative with which to explain the high Ca abundance generally found in the ICM.

• Because of cross-calibration issues, the EPIC MOS and pn detectors measure significantly different values of temperature and most abundances. Although this leads to systematic uncertainties on their absolute values, the discrepancies are generally smaller when considering abundances relative to Fe. Moreover, it should not affect relative differences between spectra from different regions if the same instrument(s) are used. Fitting a Gaussian line and a local continuum instead of CIE models is a robust method to measure more reliable abundances.

Online material

Appendix A: Detailled data reduction

Appendix A.1: GTI filtering

In order to reduce the soft-proton (SP) background, we build good time intervals (GTI) using the light curves in the 10–12 keV band for MOS and 12–14 keV band for pn. We fit the count-rate histograms from the light curves of each instrument, binned in 100 s intervals, with a Poissonian function and we reject all time bins for which the number of counts lies outside the interval μ ± 2σ (i.e. $\mathit{\mu }\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{2}\sqrt{\mathit{\mu }}$), where μ is the fitted average of the distribution. We repeat the same screening procedure and threshold (so-called 2σ-clipping) for 10 s binned histograms in the 0.3–10 keV band because De Luca & Molendi (2004) reported episodes of particularly soft background flares. In order to get a qualitative estimation of the residual SP flare contamination, we use the Fin_over_Fout algorithm which compares the count rates in and out of the FoV of each detector (De Luca & Molendi 2004). We found that in both observations MOS 1 displays a F_in/F_out ratio higher than 1.3, meaning that the observations have been significantly contaminated by SP events. This value is still reasonable though, and a look at the filtered light curve lead us to keep the MOS 1 datasets. Furthermore, a careful modelling of convenient SP spectral components are used in our spectral fittings as well (see Appendix B).

Appendix A.2: Resolved point sources excision

The point sources in our FoV contribute to the total flux and may bias the astrophysical results that we aim to derive from the cluster emission. Therefore, they should be discarded. We detect all the resolved point sources (RPS) with the SAS task edetect_chain and we proceed with a second check by eye in order to discard erroneous detections and possibly include a few missing candidates. It is common practice in extended source analysis to excise bright point sources from the EPIC data. We note, however, that many sources have fluxes below the detection limit S_cut and an unresolved component might remain (Appendix B.2).

A remaining problem is how to choose the excision radius in the best way. A very small excision radius may leave residual flux from the excised point sources while a very large radius may cut out a significant fraction of the cluster emission leading to decreased S/N. We define A_eff as the extraction region area for the cluster emission when the point sources are excised with a radius r_s, ${\mathit{A}}_{\mathrm{eff}}\mathrm{=}\mathit{A}\left(\mathrm{1}\mathrm{-}\mathit{\pi }{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{S}}_{\mathrm{cut}}}^{\mathrm{\infty }}\left(\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\right)\mathrm{d}\mathit{S}\right)\mathit{,}$(A.1)where N is the number of sources, S is the flux, and A is the full detection area.

Since we are dealing with a Poissonian process, S/N can be estimated as $\mathit{S}\mathit{/}\mathit{N}\mathrm{=}\frac{\mathit{C}}{\sqrt{\mathit{C}\mathrm{+}\mathit{B}}}$, where C and B are the number of counts of the cluster emission and the total background, respectively. The value of C depends on the extraction area and can thus be written C = C^∗A_eff, where C^∗ is the local surface brightness of the cluster (counts/′′), and B can be divided into the instrumental or hard particle (HP) background I, an unresolved point sources (UPS) component, and the remaining excised point source flux outside the excision region. The total background can be thus written as $\mathit{B}\mathrm{=}\mathit{I}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{S}}_{\mathrm{cut}}}\mathit{S}\left(\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\right)\mathrm{d}\mathit{S}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{EEF}\mathrm{(}{\mathit{r}}_{\mathrm{s}}\mathrm{\left)}\mathrm{\right)}{\mathrm{\int }}_{{\mathit{S}}_{\mathrm{cut}}}^{\mathrm{\infty }}\mathit{S}\left(\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\right)\mathrm{d}\mathit{S,}$(A.2)where EEF(r_s) is the encircled energy fraction of the PSF as a function of radius. We can finally write the total S/N as $\begin{array}{ccc}\mathit{S}\mathit{/}\mathit{N}\mathrm{=}\frac{{\mathit{C}}^{\mathrm{\ast }}\sqrt{\mathit{A}\left(\mathrm{1}\mathrm{-}\mathit{\pi }{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{S}}_{\mathrm{cut}}}^{\mathrm{\infty }}\mathrm{(}{\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{)}}\mathrm{d}\mathit{S}\right)}}{\sqrt{{\mathit{C}}^{\mathrm{\ast }}\mathrm{+}\mathit{I}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{S}}_{\mathrm{cut}}}\mathit{S}\mathrm{(}{\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{)}}\mathrm{d}\mathit{S}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{EEF}\mathrm{(}{\mathit{r}}_{\mathrm{s}}\mathrm{\left)}\mathrm{\right)}{\mathrm{\int }}_{{\mathit{S}}_{\mathrm{cut}}}^{\mathrm{\infty }}\mathit{S}\mathrm{(}{\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{)}}\mathrm{d}\mathit{S}}}\mathrm{·}& & \end{array}$(A.3)The optimum S/N can be then computed as a function of r_s and S_cut (Eq. (A.3)). In Appendix B.2 we discuss the origin of dN/ dS. We find and adopt an optimised radius for RPS excision in our dataset of ~10′′.

Appendix A.3: RGS spectral broadening correction from MOS 1 image

Because the RGS spectrometers are slitless and the source is spatially extended in the dispersion direction, the RGS spectra are broadened. The effect of the broadening of a spectrum by the spatial extent of the source is given by $\mathrm{\Delta }\mathit{\lambda }\mathrm{=}\frac{\mathrm{0.138}}{\mathit{m}}\mathrm{\Delta }\mathit{\theta }\AA \mathit{,}$(A.4)where m is the spectral order and θ is the offset angle in arcmin (see the XMM-Newton Users Handbook).

The MOS 1 DET Y direction is parallel to the RGS dispersion direction. Therefore, we extract the brightness profile of the source in the dispersion direction from the MOS 1 image and use this to account for the broadening following the method described by Tamura et al. (2004). This method is implemented through the Rgsvprof task in SPEX. As an input of this task, we choose a width of 10′ around the core and along the dispersion axis, in which the cumulative brightness profile is estimated. In order to correct for continuum and background, we use a MOS 1 image extracted within 0.5–1.8 keV (i.e. the RGS energy band). This procedure is applied to both observations and we average the two spatial profiles obtaining a single profile that will be used for the stacked RGS spectrum.

Appendix B: EPIC background modelling

We split the total EPIC background into two categories, divided further into several components:

• 1.

  Astrophysical X-ray background (AXB), from the emission ofastrophysical sources and thus folded by the response files. TheAXB includes the Local Hot Bubble (LHB), the galactic thermalemission (GTE), and the UPS.

• 2.

  Non-X-ray background (NXB), consisting of soft or hard particles hitting the CCD chips and considered as photon events. For this reason, they are not folded by the response files. The NXB contains the SP and the HP backgrounds.

In total, five components are thus carefully modelled.

Appendix B.1: Hard particle background

High energy particles are able to reach the EPIC detectors from every direction, even when the filter wheel is closed. Besides continuum emission, they also produce instrumental fluorescence lines which should be carefully modelled. Moreover, for low S/N areas, we observe a soft tail in the spectra due to the intrinsic noise of the detector chips. A good estimate of the HP background can be obtained by using Filter Wheel Closed (FWC) data which are publicly available on the XMM-Newton SOC webpage⁵. We select FWC data that were taken on 1 October 2011 and 28 April 2011 with an exposure time of 53.7 ks and 35.5 ks for MOS and pn, respectively. We removed the MOS 1 events from CCD3, CCD4, and CCD6 to be consistent with our current dataset.

Instead of subtracting directly the FWC events from our observed spectra, modelling the HP background directly allows a much more precise estimate of the instrumental lines fluxes, which are known to vary across the detector (Snowden & Kuntz 2013).

We fit the individual FWC MOS and pn continuum spectra with a broken power law F(E) = YE^{− Γ}e^η(E) where η(E) is given by $\mathit{\eta }\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\frac{\mathit{r\xi }\mathrm{+}\sqrt{{\mathit{r}}^{\mathrm{2}}{\mathit{\xi }}^{\mathrm{2}}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathrm{)}}}{\mathrm{1}\mathrm{-}{\mathit{r}}^{\mathrm{2}}}$(B.1)with ξ = ln(E/E₀) and $\mathit{r}\mathrm{=}\frac{\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Delta \Gamma }{\mathrm{)}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}}{\mathrm{\left|}\mathrm{\Delta \Gamma }\mathrm{\right|}}$ (see SPEX manual). In this model, the independent parameters are A, Γ (spectra index), ΔΓ (spectral index break), E₀ (break energy), and b (break strength). Unlike the instrumental lines, this continuum does not vary strongly across the detector. Tables B.1 and B.2 show the best-fit parameters that we found for the entire FoV extraction area and the modelled instrumental lines, respectively. In addition to the broken power-law, each instrumental line is modelled with a narrow (FWHM ≤ 0.3) Gaussian function. Although a delta function is more realistic, in this case allowing a slight broadening optimises the correction for the energy redistribution on the instrumental lines.

Appendix B.2: Unresolved point sources

An important component of the EPIC background is the contribution of UPS to the total X-ray background. Its flux can be estimated using the so called log N–log S curve derived from blank field data. This curve describes how many sources are expected at a certain flux level. The source function has the form of a derivative (dN/ dS) and can be integrated to estimate the number of sources in a certain flux range, $\mathit{N}\mathrm{(}\mathit{<}\mathit{S}\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathit{S}}^{\mathrm{\infty }}\left(\frac{\mathrm{d}{\mathit{N}}^{\mathrm{\prime }}}{\mathrm{d}{\mathit{S}}^{\mathrm{\prime }}}\mathrm{d}{\mathit{S}}^{\mathrm{\prime }}\right)\mathit{,}$(B.2)where N is the number of sources and S is the low-flux limit.

The most common bright UPS are AGNs, but galaxies and hot stars contribute as well. Based on the Chandra deep field, Lehmer et al. (2012) find that AGNs are the most dominant in terms of number counts, but in the 0.5–2 keV band the galaxy counts become higher than the AGN counts below a few times 10^-28 W m^-2 deg^-2. The assumed spectral model of this component is a power-law with a photon index of Γ = 1.41 (see e.g. Moretti et al. 2003; De Luca & Molendi 2004). In reality, the power-law index may vary slightly between 1.4–1.5, given the uncertainties in the different surveys and estimations (Moretti et al. 2009). Based on the Chandra Deep Field South (CDF-S) data, Lehmer et al. (2012) define the (dN/ dS) relations for each source category as follows: $\begin{array}{ccc}& & {\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{AGN}}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}\mathit{,}\\ & & {\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{gal}}\mathrm{=}{\mathit{K}}^{\mathrm{gal}}\mathrm{(}\mathit{S}\mathit{/}{\mathit{S}}_{\mathrm{ref}}{\mathrm{)}}^{\mathrm{-}{\mathit{\beta }}^{\mathrm{gal}}}\mathit{,}\\ & & {\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}}^{\mathrm{star}}\mathrm{=}{\mathit{K}}^{\mathrm{star}}\mathrm{(}\mathit{S}\mathit{/}{\mathit{S}}_{\mathrm{ref}}{\mathrm{)}}^{\mathrm{-}{\mathit{\beta }}^{\mathrm{star}}}\mathit{.}\end{array}$Each relation describes a power law with a normalisation constant K and a slope β. Since the (dN/ dS) relation of AGNs shows a break, there is an additional β₂ parameter and a break flux f_b. The reference flux is defined as S_ref ≡ 10^-14 erg cm^-2 s^-1. The best-fit parameters for the studied energy bands are listed in Table 1 of Lehmer et al. (2012).

The relations above can be used to estimate the flux from sources that are not detected in our EPIC observations. The UPS component also holds for the deepest Chandra observations. Hickox & Markevitch (2006) found a detection limit of 1.4 × 10^-16 in a 1 Ms CDF-S observation and estimated the unresolved flux to be (3.4 ± 1.7) × 10^-12 erg cm^-2 s^-1 deg^-2 in the 2–8 keV band. Since Chandra has a much lower confusion limit and a narrower PSF, we do not expect EPIC to reach this detection limit even in a deep cluster observation. It is therefore not necessary to know the log N–log S curve below this flux limit to obtain a reasonable estimate for the unresolved flux.

In the flux range from 1.4 × 10^-16 up to the EPIC flux limit, we can calculate the flux using the log N–log S relation. The total unresolved flux Ω_UPS for the 2–8 keV band is then calculated using ${\mathrm{\Omega }}_{\mathrm{UPS}}\mathrm{=}\mathrm{3.4}\mathrm{\times }{\mathrm{10}}^{-12}\mathrm{+}{\mathrm{\int }}_{\mathrm{1.4}\mathrm{\times }{\mathrm{10}}^{-16}}^{{\mathit{S}}_{\mathrm{cut}}}{\mathit{S}}^{\mathrm{\prime }}\left(\frac{\mathrm{d}\mathit{N}}{\mathrm{d}{\mathit{S}}^{\mathrm{\prime }}}\right)\mathrm{d}{\mathit{S}}^{\mathrm{\prime }}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}\mathrm{de}{\mathrm{g}}^{-2}\mathit{.}$(B.6)Using the Eqs. (B.3)–(B.5)for $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}$ in the integral above, the unresolved flux calculation is straightforward. Given the detection limit of our observations S_cut = 3.83 × 10^-15 W m^-2, we find a total UPS flux of 8.07 × 10^-15 W m^-2 deg^-2. This value can be used to constrain the normalisation of the power-law component describing the AXB background in cluster spectral fits. We note that this method does not take the cosmic variance into account (see e.g. Miyaji et al. 2003), which means that the normalisation may still be slightly biased.

Appendix B.3: Local Hot Bubble and Galactic thermal emission

The LHB component is thought to originate from a shock region between the solar wind and our local interstellar medium (Kuntz & Snowden 2008), while the GTE is the X-ray thermal emission from the Milky Way halo. At soft energies (below ~1 keV), the flux of these two foreground components is significant. They are both known to vary spatially across the sky, but we assume that they do not change significantly within the EPIC FoV. Both components are modelled with a CIE component where we assume the abundances to be proto-solar. Both temperatures are left free, but are expected to be within 0.1–0.7 keV. The GTE component is absorbed by a gas with hydrogen column density (N_H = 1.26 × 10²⁰ cm^-2), while the LHB component is not.

Appendix B.4: Residual soft-proton component

Even after filtering soft flare events from our raw datasets, a quiescent level of SP remains that might affect the spectra, especially at low S/N and above ~1 keV. It is extremely hard to precisely estimate the normalisation and the shape of its spectrum since SP quiescent events strongly vary with detector position and time (Snowden & Kuntz 2013). They may also depend on the attitude of the satellite. For these reasons, blank sky XMM-Newton observations are not good enough for our deep exposures. The safest way to deal with this issue is to model the spectrum by a single power law (Snowden & Kuntz 2013). Using a broken power law might be slightly more realistic, but the number of free parameters is then too high to make the fits stable. The spectral index Γ of the power law is unfortunately unpredictable and may be different for MOS and pn instruments and between different observations. Since Snowden & Kuntz (2013) reported spectral indices between ~0.1–1.4, we allow the Γ parameter in our fits to vary within this range.

Appendix B.5: Application to our datasets

We apply the procedure described above for each component on our two observations of A 4059. We extract an annular region with inner and outer radii of 6′ and 12′, respectively, and centred on the cluster core (Fig. 1, the outer two annuli), assuming that all the background components described above contribute to the detected events covered by this area. In order to get a better estimation of the foreground thermal emission (GTE and LHB), we fit a ROSAT PSPC spectrum from Zhang et al. (2011) simultaneously with our EPIC spectra. This additional observation covers an annulus centred to the core and with inner and outer radii of 28′ (~r₂₀₀) and 40′ (~r₂₀₀ + 12′), respectively, avoiding instrumental features and visible sources. We note that in this fit we also take the UPS contribution into account. Depending on the extraction area, all the normalisations (except for the UPS component, Appendix B.2) are left free, but are properly coupled between each observation and instrument.

Table B.3 shows the different background values that we found for the extracted annulus. Figure B.1 shows the result for the MOS 2 spectrum at the first observation, its best fit model, and the contribution of every modelled component. As expected, the NXB contribution is more important at high energies. Above ~5 keV, the cluster emission is much smaller than the HP background. Consequently and as already reported, the temperature and abundances measured by EPIC are harder to estimate in the outer parts of the FoV.

Finally, we apply and adapt our best background model to the core region (Sect. 4) and the eight concentric annuli (Sect. 5). The normalisation of every background component has been scaled and corrected for vignetting if necessary. From the background parameters, only the normalisations of the HP component (initially evaluated from the 10–12 keV band, where negligible cluster emission is expected), as well as those of the instrumental fluorescent lines, are kept free for all the spectra. In the outermost annulus (9′–12′) we ignore the channels below 0.4 keV (MOS) and 0.6 keV (pn) to avoid low energy instrumental noise. For the same reason we ignore the channels below 0.4 keV (MOS) and 0.5 keV (pn) in the second outermost annulus (6′–9′). The background is also applied to and adapted for the analysis of the spectra of each map cell (Sect. 6).

Appendix C: S/N requirement for the maps

Despite their good statistics, we want to optimise the use of our data and find the best compromise between the required spatial resolution of our maps (Sect. 6) and S/N. The former is necessary when searching for inhomogeneities and kT/metal clumps (i.e. the smaller the better), the latter to ensure that the associate error bars are small enough to make our measurement significant. Clearly, these variables depend on the properties of the cluster and on the exposure time of our observations.

We perform a set of simulations to determine what the best combination of S/N and spatial resolution is for the case of A 4059. For every annulus (i.e. the ones determined in Sect. 5) we simulate a spectrum with input parameters (i.e. kT, O, Ne, Mg, Si, S, Ca, Fe, Ni, and the normalisation) corresponding to the ones determined in the radial profiles analysis.

The AXB and the HP background are added to the total spectrum by using the properties derived in Appendix B. We allow their respective normalisations to vary within ±3% in order to take into account spatial variations on the FoV. Starting from the value we derived for the radial profile, we rescale the normalisation of the simulated spectrum to the particular spatial resolutions we are interested in (here we test 15′′, 20′′, 25′′, 30′′, 40′′, 50′′, and 60′′). We then fit the spectrum as done for the real data and for all the annuli and spatial resolutions we calculate the relative errors on the temperature and Fe abundance as a function of S/N. The median values of 300 realisations are shown in Fig. B.2 with their 1σ errors.

A S/N of 100 is required to measure the abundance with a relative error lower than ~20%. With this choice the temperature will be also determined with a very good accuracy, i.e. relative errors always lower than ~5%.

---

Acknowledgments

This work is a part of the CHEERS (Chemical Evolution RGS cluster Sample) collaboration. We would like to thank its members and the anonymous referee for their feedback and discussions. F.M. thanks Huub Röttgering and Darko Donevski for useful discussions. L.L. acknowledges support by the DFG through grant LO 2009/1-1. Y.Y.Z. acknowledges support by the German BMWi through the Verbundforschung under grant 50 OR 1304. 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 the USA (NASA). The SRON Netherlands Institute for Space Research is supported financially by NWO, the Netherlands Organisation for Space Research.

References

All Tables

All Figures