The uniformity and time-invariance of the intra-cluster metal distribution in galaxy clusters from the IllustrisTNG simulations

The distribution of metals in the intra-cluster medium encodes important information about the enrichment history and formation of galaxy clusters. Here we explore the metal content of clusters in IllustrisTNG - a new suite of galaxy formation simulations building on the Illustris project. Our cluster sample contains 20 objects in TNG100 - a ~(100 Mpc)^3 volume simulation with 2x1820^3 resolution elements, and 370 objects in TNG300 - a ~(300 Mpc)^3 volume simulation with 2x2500^3 resolution elements. The z=0 metallicity profiles agree with observations, and the enrichment history is consistent with observational data going beyond z~1, showing nearly no metallicity evolution. The abundance profiles vary only minimally within the cluster samples, especially in the outskirts with a relative scatter of ~15%. The average metallicity profile flattens towards the center, where we find a logarithmic slope of -0.1 compared to -0.5 in the outskirts. Cool core clusters have more centrally peaked metallicity profiles (~0.8 solar) compared to non-cool core systems (~0.5 solar), similar to observational trends. Si/Fe and O/Fe radial profiles follow positive gradients. The outer abundance profiles do not evolve below z~2, whereas the inner profiles flatten towards z=0. More than ~80% of the metals in the intra-cluster medium have been accreted from the proto-cluster environment, which has been enriched to ~0.1 solar already at z~2. We conclude that the intra-cluster metal distribution is uniform among our cluster sample, nearly time-invariant in the outskirts for more than 10 Gyr, and forms through a universal enrichment history.

2014). These first results were followed by much deeper and better resolved observations with XMM-Newton, Chandra, Suzaku and Hitomi, which allowed a characterisation of the spatial distribution of metals and their evolution in much greater detail (e.g., Werner et al. 2006;Sanderson et al. 2009;Matsushita 2011;Bulbul et al. 2012a;Molendi et al. 2016;McDonald et al. 2016;. Those studies have revealed a striking universality and constancy of the metallicity and elemental abundance distributions within the ICM. Except for the innermost parts of the ICM most clusters have similar metallicity profiles with shallow negative gradients. The detailed metallicity profile near the cluster center depends on the thermodynamic profile. For cool core clusters the metallicity peaks at the center whereas non-cool core clusters exhibit shallow gradients all the way towards the center without such a peak (e.g., De Grandi & Molendi 2001;De Grandi et al. 2004;Baldi et al. 2007;Leccardi & Molendi 2008;Johnson et al. 2011;Elkholy et al. 2015). Consequently, the level of core entropy anti-correlates with the central enrichment level (Leccardi et al. 2010). The central excess in cool core clusters could, for example, have been produced by the stellar population of the brightest cluster galaxy (Böhringer et al. 2004). Most observations also find that there is only a weak dependence of the average ICM metallicity on the cluster temperature; e.g. Baldi et al. (2012) find Z ∝ kT 0.06±0.16 ; i.e. the average metallicity does not depend strongly on cluster mass or X-ray luminosity although Yates et al. (2017) recently found a slight anti-correlation for the temperatureiron abundance relation.
More recent observations also quantified the cluster enrichment history to higher redshifts finding that the median average metallicity since z ∼ 1 changes by less then 40% compared to its present-day value (e.g., Baldi et al. 2012;Ettori et al. 2015;Mc-Donald et al. 2016;; i.e. there is no strong evolutionary trend visible in metallicity going back in time. Observations of the outskirts of galaxy clusters (see Reiprich et al. 2013, for a review) also revealed remarkable uniformity with azimuth (Werner et al. 2013;Simionescu et al. 2015;Urban et al. 2017) supporting the idea that those regions should have been enriched early on.
The iron abundance within the ICM can be constrained rather well through Fe-K and Fe-L emission measurements. However, quantifying the abundance of elements other than iron is more challenging and more uncertain. The most recent studies tend to find typically rather flat X/Fe ratios (Simionescu et al. 2015;Mernier et al. 2017;Ezer et al. 2017;Simionescu et al. 2017) pointing towards a similar distribution of iron and lighter elements in the ICM supporting the idea of an early ICM enrichment.
Observations of the detailed distribution of metals in the ICM provide important information to help understand the metal production and cluster formation processes. Core-collapse supernovae (SNcc) are mainly responsible for the production of light elements such as oxygen, neon, magnesium and silicon. SNIa, on the other hand, produce large amounts of heavier elements like iron and nickel (for a review, see Nomoto et al. 2013). These two distinct production channels furthermore operate on different timescales in releasing the metals. Observations of relative abundance ratios can therefore help to constrain the various enrichment processes in the ICM. However, clusters provide a rather complex astrophysical and dynamical environment, where lots of different physical processes like galactic winds, active galactic nuclei (AGN) feedback, and gas stripping lead to mixing and a redistribution of metals in the ICM (Churazov et al. 2001;Rebusco et al. 2005;Simionescu et al. 2008Simionescu et al. , 2009Kirkpatrick & McNamara 2015). Despite this complexity, a key interpretation of the uniformity of the iron distribution and the flat abundance ratio profiles within the ICM is that early enrichment (z > 2) might be required to achieve this level of homogeneity throughout the ICM. Otherwise, more recent astrophysical events within the ICM should cause more significant metal abundance variations. Indeed, some clusters show a metallicity distribution indicative of such recent events with fluctuations in the abundance patterns. For example, Kirkpatrick & McNamara (2015) found higher metallicities along the axes of X-ray cavities, which could be caused by AGN feedback.
Some simple theoretical attempts of chemical and population synthesis modelling failed to simultaneously match the constraints on ICM metallicity and on the stellar mass to light ratios (Renzini & Andreon 2014). These authors claim that the iron yield in the largest clusters (M500,crit > 10 14 M ) must be a few times the solar value. Other studies came to similar conclusions (Nagashima et al. 2005;Arrigoni et al. 2010). Various solutions to this modelling problem have been suggested; e.g., changes to the shape of the initial stellar mass function, increasing the efficiency of iron production in SNIa, or more efficient metal ejection from galaxies. Hydrodynamical simulations on the other hand provide a more selfconsistent and more complete approach for studying the growth and enrichment of galaxy clusters. Over the last decade those simulations had have mixed success in explaining some of the observed metallicity relations and abundance patterns (Tornatore et al. 2007;Cora et al. 2008;Rasia et al. 2008;Fabjan et al. 2010;McCarthy et al. 2010;Planelles et al. 2014;Rasia et al. 2015;Martizzi et al. 2016;Biffi et al. 2017). Martizzi et al. (2016), for example, find too low metallicities in the ICM by up to factors of about ∼ 2 − 5 depending on cluster-centric distance. Rasia et al. (2015) and Biffi et al. (2017), on the other hand, find reasonable agreement with the observed iron and metallicity profiles. However, they have to resort to artificial thermal conduction, and this heat transport can flatten the metal distribution as Kannan et al. (2017) recently pointed out. It is therefore important to gauge the effects of thermal conduction properly before drawing conclusions. Other simulations, for example those presented in Tornatore et al. (2007), correctly capture aspects of the ICM enrichment, but do not convincingly demonstrate that they also predict the correct stellar masses of galaxies within clusters; i.e. they do not provide a consistent picture. Besides the actual galaxy formation model, results might also depend on the numerical scheme and resolution as discussed by Martizzi et al. (2016). So far no hydrodynamical simulation could simultaneously achieve the following goals: produce a reasonable galaxy population from dwarf to cluster scales (e.g., correct stellar masses, galaxy colors, clustering); produce the right amount of ICM enrichment (e.g., correct metallicity gradients and overall enrichment in clusters, abundance ratios, cool core/non-cool core dichotomy, enrichment history); and demonstrate the latter two points over a large unbiased sample of numerically well-resolved galaxy clusters.
The goal of this paper is to fill this gap and present a first assessment of the metal content of the ICM of galaxy clusters in Illus-trisTNG -an ongoing follow-up simulation campaign of the Illustris project containing various large-scale hydrodynamical galaxy formation simulations. IllustrisTNG allows us to study the chemical composition of the ICM of a large sample of galaxy clusters, which is much larger than any previously studied cluster sample at that numerical resolution. Our study is therefore distinctly different from previous theoretical attempts to study the ICM enrichment for several reasons: (i) We employ a state-of-the-art galaxy formation model that passes a large number of observational tests from dwarf to cluster scales. (ii) The numerical resolution of our cluster sam- (Ngas + N DM = 2 × 2160 3 , m b = 5.74 × 10 4 h −1 M ), TNG100 (Ngas + N DM = 2 × 1820 3 , m b = 9.44 × 10 5 h −1 M ), and TNG300 (Ngas + N DM = 2 × 2500 3 , m b = 7.44 × 10 6 h −1 M ). The table presents the basic numerical parameters of the IllustrisTNG simulations studied in this paper (TNG100 and TNG300): simulation volume, number of gas cells (Ngas), number of dark matter particles (N DM ), number of tracer particles (Ntracer), baryon mass resolution (m b ), dark matter mass resolution (m DM ), and Plummer-equivalent gravitational softening length ( ). For each volume, we have run different numerical resolutions spaced by a factor of eight in mass resolution. The gravitational softening lengths refer to the maximum physical softening length of dark matter and star particles. The softening of gaseous cells is tied to their radius and allowed to fall below this value.
ples is higher than in most of the previous studies. (iii) The cluster sample size of our simulation is one to two orders of magnitude larger than any previous attempt studying the ICM enrichment. (iv) Our cluster sample is large enough to probe significantly different assembly and merging histories thereby minimising potential biases. (v) Our simulation code has been demonstrated to be efficient and accurate in solving the hydrodynamical equations. The paper is structured as follows. In Section 2 we briefly describe our methods and the simulation suite that we are using for our analysis. In Section 3 we present a detailed analysis of the present-day metallicity distribution in the ICM of our clusters. We inspect the redshift evolution of the metal content in Section 4, where we also analyse the origin of the metals in the ICM. Finally, we present our conclusions in Section 5.
The simulations were carried out with the moving-mesh code AREPO (Springel 2010). The IllustrisTNG simulations employ a comprehensive module for galaxy formation physics, which is an updated version of the Illustris model Torrey et al. 2014;Vogelsberger et al. 2014a,b;Genel et al. 2014). Details of the updated model are described in Weinberger et al. (2017a) and Pillepich et al. (2017b). The new model features a new radio mode AGN feedback scheme (Weinberger et al. 2017a), numerical improvements for the convergence properties (Pakmor et al. 2016), the inclusion of ideal magnetohydrodynamics (Pakmor & Springel 2013), a retuned SN wind model, and refinements/extensions to the chemical evolution scheme (Pillepich et al. 2017b).
The main differences between the IllustrisTNG and Illustris galaxy formation models affecting the metal production and distribution are: the AGN feedback implementation, the galactic wind scalings, and the details of the stellar evolution and chemical enrichment model. IllustrisTNG employs a kinetic AGN feedback model in the low-accretion state, which produces black hole-driven winds, whereas Illustris models the low-accretion state with a radio bubble scheme (Sijacki et al. 2007;Vogelsberger et al. 2013). The SN-driven galactic wind characteristics have also been modified for IllustrisTNG: wind material has now a thermal component; it is injected isotropically instead of bipolar; injection occurs with a gas-metallicity dependent wind energy; winds have a minimum wind velocity and are injected with a redshift dependent wind velocity implying that the wind velocity and the growth of the virial halo mass have the same scaling with redshift. The stellar evolution and chemical enrichment schemes also differ in a few ways. Both models employ a Chabrier (2003) IMF, but the IllustrisTNG model assumes that stars pass through an AGB phase in the mass range 1 − 8 M , while stars with masses between 8 and 100 M end as SNcc. The transition in the Illustris model was 6 M . Furthermore the yield tables have also been updated going from the Illustris to the IllustrisTNG model. The AGB yields are extended from their original Illustris 1 − 6 M mass range to 1 − 7.5 M . For masses between 1−6 M IllustrisTNG uses the Karakas (2010) yields. For masses of 7 and 7.5 M IllustrisTNG employs the Doherty et al. (2014) and Fishlock et al. (2014) yields. The largest update is to the SNII tables, which for IllustrisTNG extend from 8 to 120 M : for 13 − 40 M the yields are taken from Kobayashi et al. (2006), which are then complemented to lower and higher masses with the yields presented in Portinari et al. (1998). IllustrisTNG also employs the Nomoto et al. (1997) yields for SNIa, which lead only to minor differences with respect to the SNIa yields of the Illustris model. However, the IllustrisTNG SNIa rates also differ from the 5.0 5.5 6.0 6.5 7.0 7.5 Figure 1. Overview of two of the three IllustrisTNG simulations. The main panel shows a slice of TNG300 with a thickness of 22.13 Mpc and an extent of the full TNG300 box (302.63 Mpc) at z = 0. We show the DM density field (outer part) and iron distribution (inner part). The maps are centered on the most massive cluster of TNG300 (M 200,crit = 1.54 × 10 15 M ). The iron distribution traces the large scale matter distribution with gas outflows also enriching the intergalactic medium. The upper right inset shows the total metal distribution in a 3 Mpc region around the most massive cluster in TNG300. The lower right inset shows the stellar light of the central galaxy of the corresponding cluster. The bottom inset shows the X-ray emission of the hot gas around this most massive cluster. In the upper left panel we present the DM density fields of TNG100 and TNG300 in relative scale, and the iron map of TNG100.
Illustris model based on the changes proposed in Marinacci et al. (2014). The combined changes in the stellar evolution model and yields can cause quite significant changes in the integrated metal production. For example, the amount of produced iron over one Hubble time differs by a factor of ∼ 2.
In the following we will present our metallicity results in units of the solar abundances of Anders & Grevesse (1989), which is the most commonly used system of units for presenting observational metallicity and abundance measurements. Specifically, we assume Z = 0.0194, Fe = 1.853 × 10 −3 , O = 9.593 × 10 −3 , Si = 6.998 × 10 −4 gas mass fractions when presenting metallicities, iron, oxygen, and silicon abundances. In the following we calculate based on the simulation data the total metal mass of a cell or the mass of a given element in that cell and then divide this by the  In all clusters the metallicity is increasing towards the center and reaches values around ∼ 0.1 Z at the virial radius. In some clusters stream like features are also visible due to metal stripping from infalling cluster galaxies.

Fe
Fe/Z  total gas mass in a cell to derive a mass fraction. The resulting mass fractions are the divided by the corresponding solar mass fractions. We note that more recent measurements find lower solar metallicities; e.g. Z = 0.0133 (Lodders 2003), Z = 0.0141 (Lodders et al. 2009), and Z = 0.0134 (Asplund et al. 2009). These solar values can be used to easily convert our results into other systems of units.

Metallicity and abundance ratio maps
Before describing the metal content of individual clusters, we show in Fig. 1 the large-scale distribution of DM and iron of TNG300. The iron abundance is presented in solar units as described above. This map has a side-length of 302.63 Mpc and a thickness of 22.13 Mpc. The central part shows the iron distribution, whereas the outer parts visualise the DM mass distribution. The map is focused on the most massive cluster of TNG300 (M200,crit = 1.54 × 10 15 M ). As expected the overall iron abundance follows the large-scale structure. We note that the extent of the iron distribution differs between Illustris and IllustrisTNG due to the different AGN feedback; i.e. the new kinetic AGN feedback of Illus-trisTNG does not transport metals out as far as the Illustris bubble radio mode (Weinberger et al. 2017a). Fig. 1 also contains four insets showing the metal distribution, X-ray emission and stellar light associated with the most massive cluster at the center of the projection. The fourth, upper left inset compares the size of the TNG300 volume with TNG100, for which we also show the iron distribution. We begin our detailed analysis of the ICM metal enrichment and chemical composition by visually inspecting a few massive clusters of TNG300 and TNG100. In Fig. 2 we present metallicity maps for twelve (eight) massive clusters of TNG300 (TNG100) at z = 0. These clusters span a mass range from M500,crit = 4.17 × 10 14 M to M500,crit = 1.16 × 10 15 M . The corresponding radii cover a range from r500,crit = 1 Mpc to r500,crit = 1.63 Mpc. The circles in the different panels of Fig. 2 denote r500,crit, and the virial masses (M200,crit) of the clusters are indicated in each panel. The total horizontal and vertical extent of the projection region is 3 × r500,crit. The maps in Fig. 2 demonstrate that the metal distribution in the ICM is not fully homogeneous. However, the average metallicity roughly approaches the typical values of about ∼ 0.1 − 0.5 solar, in agreement with observational findings. One can also see how the metallicity increases towards the centers of clusters.
Despite strong observational evidence for this nearly universal metallicity of the ICM, not every simulation has been able to actually reproduce those values and the observed shallow metallicity gradients in the ICM. For example Martizzi et al. (2016) find significantly too low metallicities in their simulations. We note however that the actual amount of metals in the ICM and the overall amplitude of the corresponding profiles depend on the employed stellar yields that enter the chemical evolution model. Uncertainties in these yields directly propagate into uncertainties in the amplitude of ICM metal profiles. One should therefore not over-interpret the normalisation of the profiles given the large yield uncertainties.
Inspecting the maps of the different clusters in Fig. 2 also reveals that the average metallicity seems to be rather similar across this small cluster sample. No cluster has a particularly low or high average ICM metallicity compared to the other clusters in the sample. We will quantify this similarity in more detail below. Again, this finding is in good agreement with observational data, which indicates that, for example, the kBT − Z relation of clusters has a rather shallow slope with little scatter; i.e. there is only a weak dependence of the average cluster metal content on its mass, temperature or X-ray luminosity. For example, Baldi et al. (2012) find a slope of (0.06 ± 0.16) for the mean trend of the temperaturemetallicity scaling relation. Despite the uniformity and similarity of the different maps, one can also appreciate that detailed features in each map differ from cluster to cluster. These variations are driven by the specific formation history of the cluster as well as internal processes like AGN feedback, SN winds, stripping processes, and self-enrichment within the cluster environment. Hints of these irregularities are also seen in some recent cluster observations (Kirkpatrick & McNamara 2015;Vagshette et al. 2016).
Such small-scale features and fluctuations should also be present in the abundance patterns of other elements in the ICM. Our chemical evolution model individually tracks a variety of chemical elements (H, He, C, N, O, Ne, Mg, Si, Fe) beyond pure metallicity. Maps of some of these elements are presented in Fig. 3, where we show four different views of the most massive cluster of TNG300 with M200,crit = 1.54 × 10 15 M . The maps present from upper left to lower right: the iron (Fe) abundance, the iron over total metallicity (Fe/Z) ratio, the silicon over iron (Si/Fe) ratio, and the oxygen over iron (O/Fe) ratio as indicated. These maps can directly be compared to the first map of Fig. 2, which shows the corresponding distribution of all metals. Any visible differences in these maps must be attributed to the different production and distribu-tion mechanisms for iron, oxygen, silicon compared to the general metal distribution mapped by Z. Silicon and oxygen have different production channels compared to iron, with their main contribution coming from SNcc, whereas iron is mostly produced through SNIa. SNIa produce a large amount of iron and nickel, while SNcc are the main contributors of oxygen and silicon. Lighter elements (carbon, nitrogen) are also produced by low-and intermediate-mass stars during their asymptotic giant branch (AGB) phase. Despite these differences the maps show overall rather similar structures as the metallicity maps in Fig. 2. The abundance of the elements is also quite close to solar. However, the ratio maps demonstrate that the actual abundance ratios of elements can fluctuate within the ICM. For example, the silicon over iron and oxygen over iron ratios are smaller at the cluster center compared to its outskirts at r500,crit. Furthermore, silicon over iron and oxygen over iron ratios are low at the centers of strongly emitting regions. On average the Si/Fe values are larger than O/Fe. The fact that the Si/Fe and O/Fe ratios slightly decrease towards the center can be interpreted as a slight decrease in enrichment due to SNcc compared to SNIa for smaller radii. We will show below that this is indeed the case for our cluster samples in the IllustrisTNG simulations.

Radial metallicity profiles
Observations with XMM-Newton, Chandra and Suzaku have established a variety of insights into the metal content and metal distribution within galaxy clusters. A general result is that radial profiles of the iron abundance show shallow negative gradients, steeper for relaxed cool core clusters, with central metallicity values approaching nearly solar abundances and with a global average enrichment at a level of about ∼ 1/3 of the solar value. It is therefore interesting to compute metallicity and iron profiles of our cluster samples of TNG300 and TNG100 and contrast our simulation predictions with these observational findings.
We present those profiles along with observational data in Fig. 4. For definiteness, we present here and in the following massweighted three-dimensional averages of non-star forming gas in the clusters. We note that we explicitly tested that those profiles are very similar to projected emission-weighted profiles (see also Biffi et al. 2017), and that any differences are significantly smaller than uncertainties introduced due to our choice of elemental yields. Observationally CCD-resolution X-ray spectra from Chandra and XMM-Newton put constraints on the ICM metal abundance profiles, which are inferred from equivalent width measurement of the Fe Kα emission line at 6.7 keV. In Fig. 4 we show both the metallicity and iron profile since the latter is observationally the more fundamental quantity. The lines in Fig. 4 show the median over all mass-weighted cluster profiles for all clusters with M500,crit > 10 13.75 M . In TNG100 our cluster sample contains 20 clusters, and in TNG300 we find 370 clusters fulfilling the M500,crit cut in cluster mass. TNG300 offers significantly more statistics due to the increased simulation volume, which is about 20 times larger than the TNG100 volume. The top panel of Fig. 4 presents the results for the small cluster sample of TNG100, while the middle panel shows the significantly larger cluster sample of TNG300. We note that we apply the same cluster mass cut for both simulations implying that the mean mass of the cluster samples of TNG100 and TNG300 differ due to the different simulation volumes. The shaded colored regions indicate in each panel the 1σ spread between the different individual cluster profiles by calculating the 16% and 84% percentiles. All quoted 1σ ranges in the following will be computed using this percentile range.  For TNG300-1 we also present the resolution corrected (corr) profile. The shaded regions denote the 1σ spread among the different cluster profiles of the samples. Observational data bands are taken from Leccardi & Molendi (2008) and Ettori et al. (2015); data points are mean values taken from Leccardi & Molendi (2008). The TNG100 profiles agree with the observational data within the error bars. At all radii, especially in the outer parts, the scatter between the cluster sample is small demonstrating that the metallicity profiles are rather uniform across the full cluster sample. The grey dashed line shows a scaled version of the empirical iron profile fit in Mernier et al. (2017). The bottom panel shows the negative logarithmic slopes of the metallicity and iron profiles (solid lines), and the relative 1σ metallicity scatter as a function of radius (dashed line). The grey line shows the slope of the empirical iron profile fit in Mernier et al. (2017).
For the comparison to observations we take data points from the Leccardi & Molendi (2008) XMM-Newton sample of 50 objects, which contains a mixture of cool core and non-cool core clusters with redshifts between 0.1 and 0.4 and average temperatures kT > 3.3 keV. We compare this observational sample to our simulation data at z = 0.2, which is the mean redshift of the observational sample. We convert the r180,crit data points of Leccardi & Molendi (2008) to r500,crit radii by assuming r500,crit 0.6 × r180,crit, which on average holds for our cluster sample. Besides this observational sample we also compare to more recent observational data from Ettori et al. (2015). They provide a combined analysis of the metal content of 83 objects in the redshift range 0.09 − 1.39 and with a gas temperature between 2 keV and 12.8 keV. Their data is spatially-resolved in three radial bins with (0 − 0.15, 0.15 − 0.4, > 0.4) × r500,crit, and obtained with a similar analysis technique using XMM-Newton observational data as in Leccardi & Molendi (2008) and Baldi et al. (2012). For the Ettori et al. (2015) data we take the range of metallicity values for z < 0.2 in each radial bin. The grey shaded region in Fig. 4 shows the expected range of metallicities for clusters within that redshift range.
The metallicity profiles of TNG100-1 agree with the observational data of Leccardi & Molendi (2008) within the observational error bars and also lie well within the metallicity ranges of Ettori et al. (2015). We note that IllustrisTNG (TNG100-1) correctly predicts the amplitude and overall shape of the observed metallicity profiles. More specifically, the metallicity and iron profiles presented in Fig. 4 typically slowly rise towards the center of the clusters; i.e. they show a negative radial gradient following the trends also seen in observations. At the smallest radii they reach values around ∼ 0.6 solar in agreement with observational results. At a radial distance of about r500,crit we find average metallicities of around ∼ 0.2 solar for the TNG100 sample and ∼ 0.1 for the TNG300 sample. We note that the simulated profiles are rather steep beyond r500,crit, which is not fully consistent with observations. We will show below that the detailed metal distribution is affected by the details of the galaxy formation model even at these radii. The difference in normalisation between the TNG100-1 and TNG300-1 profiles is due to numerical resolution; i.e. the numerical resolution of TNG300 is not sufficient to achieve metallicity profile convergence at the level of TNG100-1. Based on convergence studies, presented below, we expect that the TNG300-1 profiles should be shifted towards larger metallicities by a factor of 1.6 to match the numerical resolution of TNG100-1. The corrected profiles are shown as solid lines in the middle panel of Fig. 4, and the uncorrected profiles as dashed lines. We will in the following correct the TNG300-1 results by this factor to extrapolate to TNG100-1 resolution. The bottom panel of Fig. 4 shows the logarithmic slopes of the metallicity and iron profile (solid lines). The inner region shows a flattening of the iron profile with a slope approaching ∼ −0.2 at a percent of r500,crit, whereas the profiles steepen in the outer parts reaching slopes of ∼ −0.5 around ∼ r500,crit. Interestingly, the iron profile is at all radii steeper than the total metallicity profile. We will find a similar result below when inspecting the radial profiles of silicon and oxygen. Those are also shallower than the iron profile, and contribute to the shallower behaviour of the overall metallicity profile. The gray line represents the iron slope derived from the empirical iron profile fits to stacked iron measurements presented in Mernier et al. (2017) based on measurements of 44 nearby cool core galaxy clusters, groups, and ellipticals taken from the CHEERS catalogue. We also include a scaled version of this fit, with functional form . Median metallicity profiles for TNG100-1 divided into cool core and non-cool core clusters at z = 0.1, where we require a cool core cluster to have a central entropy of less than 30 keV cm −2 . Observational data bands are taken from Ettori et al. (2015). The median cool core and noncool core metallicity profiles of the simulation are in broad agreement with the observational results and show a similar trend; i.e. cool core clusters show rising metallicity profiles towards the center. In agreement with observations the profiles agree beyond 0.15 × r 500,crit . Towards the center the cool core profiles approach metallicities nearly twice as large as those of the non-cool core systems. Dashed lines show the 1σ scatter for both the cool core and non-cool core cluster samples.
massive cluster of TNG300 measured at 1.4 × r500,crit. Observational data points are taken from eight azimuthal directions of the Perseus cluster based on 84 Suzaku observations (Werner et al. 2013). The overall fluctuation strength of the iron abundance agrees between the simulation prediction and observations. For the comparison we have identified the north direction of the Perseus observational data with 130 degree as indicated by the vertical dashed line. The shaded regions shows the envelope of all profiles for the ten most massive clusters, which all show a similar level of fluctuations at that radial distance. Both observations and our simulation consistently find a rather low fluctuation level in the outskirts of galaxy clusters. Typically, this is interpreted as a common origin of the metals in the outskirts due to some early enrichment within the proto-cluster environment. We will get back to this point below. So far we have not yet divided our cluster samples into subsets according to mass or cool core versus non-cool core for the analysis of the metallicity or iron profiles. Observationally cluster metallicity profiles, average metallicities and iron abundances do not seem to be strongly dependent on cluster mass or temperature, which we will also demonstrate below. However, observations do find that the metallicity profile depends on the thermodynamic profile of clusters. Specifically, cool core clusters seem to have increased metallicities at their centers compared to noncool core clusters. These differences in metallicity profiles for cool core versus non-cool core clusters are only found for the innermost regions of clusters (r < 0.15 × r500,crit), where central enrichment increases the metallicity and element abundances for cool core clusters. Observationally cool core and non-cool core clusters seem to have similar metallicity and iron profiles beyond that radius. To explore this dichotomy in our simulations, we show in Fig. 6 the median metallicity profiles for TNG100-1 divided into cool core and non-cool core clusters. We require a cool core cluster to have a central entropy of less than 30 keV cm −2 (Hudson et al. 2010). The detailed thermodynamical profiles of clusters require sufficient numerical resolution. We therefore do not present a similar analysis for TNG300 due to its eight times lower mass resolution. We compare the cool core and non-cool core subsample profiles of TNG100-1 to observational data taken from Ettori et al. (2015), where we select the low-redshift observations, z < 0.2, which we then compare to the z = 0.1 predictions of the simulation. We note that their selection criterion for cool core clusters is based on the pseudo-entropy ratio. The cool core and non-cool core profiles of the simulation are in broad agreement with the observational results and show a similar trend; i.e. cool core clusters show rising metallicity profiles towards the center. Cool core clusters approach 0.8 Z towards their centers whereas non-cool core clusters have central metallicities below 0.5 Z . We also find that the median metallicity profiles of the cool core and non-cool core subsamples agree beyond 0.15 × r500,crit similar to the trends seen in observations. Dashed lines show the 1σ scatter for both the cool core and non-cool core cluster samples. We note that we have explored here only one cool core versus non-cool core criterion based on the central electron density. Other criteria like central cooling time or entropy along with a more detailed census of cool core versus non-cool core systems TNG300-1 is discussed in Barnes et al. (2017a). We further note that the cluster sample in Fig. 6 is rather small, but the simulation profiles indicate that the sample shows, at least qualitatively, the observed dichotomy in the metallicity profiles.
Besides total metallicity and iron our galaxy formation model also traces other elements as mentioned above. Here we specifically focus on silicon and oxygen since they are, compared to iron, also produced by SNcc on rather short timescales, whereas iron is mostly produced by SNIa that operate on much longer time scales. Differences in the iron versus lighter element profiles therefore encodes information on the metal production process and cluster formation history in general. Observationally it appears that the typical fraction of SNIa (SNcc) contributing to the enrichment lies within ∼ 20 − 45% (55 − 80%), depending on the selected yield models (Mernier et al. 2017). Bulbul et al. (2012c) finds a 30 − 37% SNIa fraction in the core of A3112. Observations further find that the silicon over iron and oxygen over iron ratio distributions seem to be rather flat, which suggests a flat proportion and a similar metal contribution of SNcc and SNIa (Mernier et al. 2017).
Given the wealth of information encoded in the different element distributions, we present in the left two panels of Fig. 7 the median profiles for silicon and oxygen. For TNG300 we show the raw profiles (dashed) and the resolution corrected profiles (solid lines). The general shape of the silicon and oxygen profiles coarsely follow those of iron and the total metallicity; i.e. the abundance of these elements is increasing towards the center. At around ∼ 0.1 × r500,crit we find typical silicon abundances of ∼ 0.4 solar and oxygen abundances of ∼ 0.3 solar. Those values drop towards larger radii, ∼ r500,crit, down to ∼ 0.1 − 0.2 solar for oxygen and silicon. The general similarity of the oxygen and silicon distribution compared to the overall metals and iron could already be seen in Fig. 3. We compare the detailed silicon and oxygen median profiles to the recent stacking analysis of Mernier et al. (2017) who studied 44 nearby cool core galaxy clusters, groups, and ellipticals for their chemical composition using deep XMM-Newton/EPIC observations to extract the average abundance profiles of oxygen, magnesium, silicon, sulfur, argon, calcium, and iron. We find that our results lie within the scatter of measurements of the Mernier et al. (2017) observations. However, it seems that the slope of the simulated oxygen and silicon abundances differs slightly from those of the observational sample.
To quantify this in more detail, we present in the two right For TNG300-1 we present also the resolution corrected profiles (corr). Observational data is taken from Mernier et al. (2017). Right two panels: Abundance ratios of silicon and oxygen with respect to iron. The figure includes results for the small TNG100-1 (top) sample and the significantly larger TNG300-1 (bottom) cluster sample. Those ratios are resolution-independent and we therefore do not apply any resolution corrections for TNG300-1. We include observational Suzaku data points for AWM7 (Sato et al. 2008), Centaurus (Sakuma et al. 2011) and Coma (Matsushita et al. 2013) (only Si/Fe). The error bars on those measurements are large and vary with cluster-centric distance. We indicate a typical abundance error in the lower right of the panels. We also include more recent observational data from the stacking analysis of Mernier et al. (2017). The simulation ratios are slightly too high compared to the observational data. Furthermore, we also find that these ratios are declining towards the cluster center. Both ratio decrease nearly by 0.5 from r 500,crit to 10 −2 × r 500,crit . panels of Fig. 7 the actual silicon over iron and oxygen over iron ratio profiles for TNG100 (top) and TNG300 (bottom) in solar units. Since these are ratios, we do not have to apply a resolution correction for TNG300. As in the other profile plots we present here the median of the cluster samples. The ratios are rather flat with a slight positive gradient, which shows that the overall silicon and oxygen gradients are shallower than those of iron. This is consistent with Fig. 4, where we found that the iron profiles are steeper than those of the total metallicity. We compare our results to observational data from Suzaku for the massive galaxy clusters: AWM7 (Sato et al. 2008), Centaurus (Sakuma et al. 2011) andComa (Matsushita et al. 2013) (only Si/Fe). We note that the observational error bars are large, and the data is consistent with our simulation findings within those errors. We do not show the error bars for each data point, but instead show the typical abundance ratio error in the lower right of each panel. Mernier et al. (2017) also recently presented abundance ratio profiles, which we show for comparison. In general, we find that oxygen is suppressed compared to silicon. Interestingly, the silicon over iron and oxygen over iron ratios decrease slightly towards the cluster center in our simulations. This is different from the ratio profiles of Mernier et al. (2017), which are typically flatter in the inner part of the cluster. In the simulation both ratios decrease by nearly 0.5 from r500,crit to 10 −2 ×r500,crit. We note that some observational studies also find peaked oxygen distributions (e.g. Werner et al. 2006;Bulbul et al. 2012b). However, it is difficult to reliably measure oxygen abundances. Oxygen lines are located at 0.65 and 0.56 keV, and this band is heavily contaminated by soft X-ray foreground emission coming from our own Galaxy.

Average metal content
So far we have mainly studied radial abundance profiles. However, we can also average the metallicity within clusters. We have already mentioned above that this average metallicity of the ICM is observationally found to be quite independent of the overall cluster properties like mass and temperature. Specifically, the temperaturemetallicity and X-ray-metallicity scaling relations are found to be rather flat. Furthermore, the scatter in the average metallicities is also quite small. We confirm these observational findings in Fig. 8 where we show best-fit Gaussians to the distribution of mass-weighted average ICM metallicities, iron, silicon, and oxy- For metallicity and iron we also show the actual simulation data demonstrating that Gaussians indeed describe those distributions well. The lower right panel presents the relation between average cluster temperature and metallicity both for TNG100-1 and TNG300-1 (with and without resolution correction) where averages are calculated within r 500,crit . In agreement with observational data, there is no strong dependence of the iron abundance on temperature over the full cluster sample; i.e. most clusters show actually a similar amount of enrichment within r 500,crit . However, there seems to be a slight anti-correlation which has also been pointed out by Yates et al. (2017) with a similar slope as indicated (gray line). We also show the actual data taken from Yates et al. (2017) (gray square symbols).
gen within different radial cuts (first three panels) for TNG300 with the resolution correction discussed above. For total metallicity and iron we also show the actual data to demonstrate that the Gaussians provide a good fit. The averages for each cluster were computed in a mass-weighted way for r < r500,crit (upper left panel), 0.15 × r500,crit < r < r500,crit (upper right panel), and r < 0.15 × r500,crit (lower left panel). In each panel we can see a clear ordering of the abundance of elements with iron being the least abundant and silicon being the most abundant in solar units. We also find that the innermost abundances within 0.15 × r500,crit are higher than those in the outer parts of the ICM. The different distributions in Fig. 8 also demonstrate that the spread in the element abundances over the full cluster sample of TNG300 is rather small showing the uniformity of the enrichment process for the different clusters. The best fits yield the following 1σ dispersions in solar units for the different radial cuts: In the lower right panel of Fig. 8 we present the temperature-iron abundance relation for TNG100-1 and TNG300-1, which demonstrates that the average iron abundances of our simulation indeed do not depend strongly on the average temperatures of the clusters in agreement with observational findings. We calculate both the iron abundance and temperature as mass-weighted averages within r500,crit. Actually there seems to be some slight anti-correlation, which can also be seen in the recent compilation of observational data presented in Yates et al. (2017). The gray line shows the best-fit to this anti-correlation found by Yates et al. (2017) using 79 lowredshift (z ∼ = 0.03) systems resulting in a slope of −0.26 ± 0.03, whereas points show the original data points following the averaging procedure in Barnes et al. (2017b) taking mean values if more than one estimate for the observed mass-weighted metallicity or temperature was present in the Yates et al. (2017) sample. The resolution corrected values have been scaled up by 1.6 to match the TNG100-1 resolution. However, the simulation metallicities are slightly lower, ∼ 0.2 dex, than the observational values at our fiducial resolution. A similar discrepancy has also been found by Barnes et al. (2017b) TNG100-1  TNG100-2  TNG100-3   TNG300-1  TNG300-2  TNG300-3 Illustris-1 TNG100-2 x M 1 /M 2 TNG300-1 x M 1 /M 2 Figure 9. Comparison of median cluster metallicity profiles extracted from Illustris and IllustrisTNG, and TNG300/TNG100 at three different numerical resolutions. Illustris and IllustrisTNG differ in the employed yields and details of the feedback parametrisation, which impacts the ICM enrichment resulting in rather different metallicity profiles.The Illustris galaxy formation model produces central metallicity drops, which are absent in Illus-trisTNG. We speculate that this is related to the Illustris bubble radio-mode AGN feedback, which leads to more efficient uplifting of central gas and metals compared to the kinetic low-accretion rate AGN feedback of Illus-trisTNG. The profiles of the different simulation resolutions demonstrate that convergence for the ICM enrichment is rather slow and it is typically difficult to achieve fully converged results. However, a scaled up version of the TNG100-2 profile with a scaling factor corresponding to the difference in stellar mass at level 1 and level 2 resolution (M 1 /M 2 ) agrees with the TNG100-1 profile. This demonstrates that the reason for the slow convergence rate is ultimately due to the slow convergence rate in stellar mass. Therefore, also the TNG300-1 metallicity profile overlaps with the median TNG100-2 profile, and a scaled up version of TNG300-1 agrees approximately also with the profile of TNG100-1. This resolution dependence motivates a simple resolution correction by scaling up ICM abundances by M 1 /M 2 ∼ = 1.6 for TNG300-1 to the resolution level of TNG100-1.

Model variations and convergence
The results presented so far and the enrichment of the ICM depend on the actual galaxy formation model that has been employed in the simulation. The IllustrisTNG galaxy formation model is an updated version of the previous Illustris model as described above, and furthermore the TNG100 simulation retains all initial condition phases of the original Illustris volume, but with an updated cosmology. We can therefore directly compare the TNG100-1 results with those of Illustris-1 to see how strongly the ICM enrichment varies between the two simulations. Such a comparison provides important insights into the sensitivity of the ICM enrichment on the galaxy formation model. In Fig. 9 we therefore present a comparison between the median metallicity profiles derived from Illustris (ILL-1) and IllustrisTNG (TNG100-1). We note that Illustris contains only 10 clusters fulfilling the mass cut criterion M500,crit > 10 13.75 M compared to the 20 cluster within TNG100 of IllustrisTNG, most importantly due to the different σ8 normalisation. The Illustris and IllustrisTNG simulations furthermore differ in various aspects: (i) the yields have been changed between the two simulations, (ii) the low accretion state AGN feedback works differently between the two simulations with Illustris injecting thermal energy in inflating bubbles in the ICM, and IllustrisTNG employing a kinetic feedback model injecting momentum at the center of the cluster modelling a black hole driven wind, (iii) the cosmol-ogy has been updated to the more recent Planck results for Illus-trisTNG, (iv) the details of the SN winds have been modified, (v) the metal advection has been improved in IllustrisTNG. All these changes affect the median abundance profiles in the ICM, such that Illustris predicts nearly a factor ∼ 1.5 − 2 higher metallicity values compared to IllustrisTNG. Aside from the difference in normalisation we also see a difference in the shape of the metallicity profile in the inner region of the ICM. The median profile of the Illustris simulation is decreasing towards the center with a positive gradient. IllustrisTNG, on the other hand, predicts a rising metallicity profile towards the center with a negative metallicity gradient over the full radial range. We speculate that this difference is related to the modelling of AGN radio mode feedback. The radio mode feedback of Illustris drives gas out of the center of cluster ) carrying along also metals. On the contrary, the kinetic radio mode feedback of IllustrisTNG (Weinberger et al. 2017a) does not result in this gas removal and therefore also does not expel metals from the central regions of the cluster. Besides this difference in the inner part, we also find that the slope of the metallicity profile in the outer part of the cluster, around and beyond r500,crit, is shallower for the Illustris model. In fact, the IllustrisTNG model shows a kink, which is absent in Illustris. This could also be related to the difference in the AGN feedback since the Illustris model distributes metals out to significantly larger radii than the IllustrisTNG model.
Interestingly, there is also some observational evidence for decreasing metallicities towards the center of clusters, and this has been interpreted as uplifting of metals through AGN bubbles (e.g., Rafferty et al. 2013;Kirkpatrick & McNamara 2015;Panagoulia et al. 2015;Mernier et al. 2017). This seems also to be consistent with the Illustris metallicity profile, where the profile indicates that central metals have been moved to larger radii in the ICM. Mernier et al. (2017) found that ∼ 32% of the clusters in their sample exhibit these kind of metallicity drops towards the centers. For Illustris the median metallicity profiles peak at 0.03 × r500,crit reaching ∼ 0.7 Z and drop to ∼ 0.6 Z at 10 −2 × r500,crit. This is comparable to the central oxygen drops seen in Mernier et al. (2017), but larger than the observed iron drop therein. Metallicity drops could also be caused by other physical effects, for example, dust-induced metal depletion (Vogelsberger et al., in prep) or thermal conduction in the ICM (Kannan et al. 2017). We note that the observed central metallicity drops might also be a result of fitting simplified singletemperature models to a multi-phase gas. Indeed, those drops often disappear if multi-temperature (2T or 3T) plasma models are employed in the innermost regions (e.g., Buote et al. 2003).
We discussed above a few times the difficulty of numerical convergence of the ICM element abundance profiles. To make this point clearer, we also present in Fig. 9 the median metallicity profiles at all three numerical resolution levels for TNG100 and TNG300. The mass resolution differs by a factor of ∼ 64 between level 3 and level 1. Higher numerical resolution leads to more enrichment in the ICM reflected by an overall increase of the normalisation of the metallicity profiles. Convergence for the ICM enrichment is rather slow and is typically difficult to achieve (see also the discussion in Martizzi et al. 2016). Strictly speaking our results therefore represent lower limits for the amount of metals in the ICM. Interestingly, the median metallicity profile of TNG300-1 overlaps with the median profile of TNG100-2, and TNG300-2 overlaps with TNG100-3. It seems therefore that the profiles at different resolutions can be scaled to each other with a constant factor. This is indeed possible because the normalisation offsets between different resolution levels are mostly driven by differences in the amount of stellar mass that is formed at the different resolu- tion levels. This can be seen by looking at the profile of TNG100-2 scaled up with a constant factor corresponding to the ratio between the stellar masses formed at the two resolution levels (M 1 /M 2 ). To calculate this ratio, we have taken all clusters in the sample and summed up the stellar mass of each associated friends-of-friends group. This then results in total stellar masses M 1,2 for TNG100-1,2, respectively. The scaled version of TNG100-2 then overlaps with the metallicity profile of TNG100-1, which demonstrates that the reason for the low convergence rate of the metallicity profiles is ultimately due to the low convergence rate in stellar mass.
We can use the scale factor M 1 /M 2 ∼ = 1.6 between TNG100-1 and TNG100-2 to also scale up the profile of TNG300-1 to the mass resolution of TNG100-1 as we have done above. This works reasonably well as Fig. 9 demonstrates. The difference in the profile shapes between TNG100-1 and the scaled up version of TNG300-1 are most likely related to the fact that the cluster sample of TNG300 is much larger and also includes more massive clusters than the small sample of TNG100, which leads to a slightly different shape of the median profile. Nevertheless, the fact that the scaled up ver-sion of TNG300-1 agrees so well with TNG100-1 demonstrates that stellar mass is the main driver for non-convergence, and the uniformity of metal enrichment among different clusters is the reason for the approximate similarity of the profile shapes that allows a resolution scaling solely based on the stellar mass ratio. One might ask whether the normalisation differences between Illustris and Il-lustrisTNG are also mostly driven by different amounts of stellar mass. We have checked that this is not the case since the stellar mass within the cluster sample is only ∼ 10% higher in Illustris compared to IllustrisTNG, which is not sufficient to explain the offset. Furthermore, the profile shapes are obviously different such that a simple scaling cannot map one profile to the other.
The results presented in this section demonstrate that the ICM enrichment of IllustrisTNG clusters is in good agreement with observational data. However, we note that the silicon over iron and oxygen over iron abundance ratios are in slight tension with observations since they follow a shallow positive gradient which is disfavored by current observations pointing towards flatter profiles.  Figure 11. Median metallicity, iron (left), silicon and oxygen (right) profiles of TNG300 at z = 0.5, 1, 2 compared to the z = 0 profiles. The outer parts of the ICM (r ∼ r 500,crit ) show a nearly constant level of enrichment back to z ∼ 2 changing by only about ∼ 10%. The slope of the inner profiles (r < 0.1 × r 500,crit ) becomes shallower towards z = 0, and the central values decrease by nearly a factor of two from z = 2 to z = 0. Therefore, the outer cluster abundance profiles beyond 0.4 × r 500,crit are strikingly constant over time with nearly no evolution over the last ∼ 10 Gyr.

PAST ICM COMPOSITION AND ITS EVOLUTION
In the previous section we have focused on the present-day metal, iron, silicon and oxygen abundances within the ICM. One key result is the uniformity of these profiles across the large cluster samples that we studied. Additionally, observations also find only little evolution of the metal content of galaxy clusters as a function of time (e.g., Ettori et al. 2015;McDonald et al. 2016;. Therefore, both the variability from cluster to cluster and also the variability as a function of time seem to be rather small in terms of the ICM metal content. Observationally obtaining reliable radial metallicity profiles for nearby clusters is challenging. Performing this exercise at higher redshifts is even more uncertain. However, average iron abundances and metallicities can be extracted for reasonably large samples of clusters spanning significant lookback times. McDonald et al. (2016) analysed cluster observations out to z ∼ 1.5 finding nearly no evolution of the averaged cluster metallicities within r500,crit as a function of redshift, in contrast to results of Ettori et al. (2015) who found a mild redshift dependence. However, in both of these studies the evolution seems to be mostly driven by changes in the actual core region with little evolution in the outer cluster parts beyond 0.15 × r500,crit. Specifically, both studies do not find a significant redshift dependence once the core region is excluded, in which case they also do not find any significant trend with cool core versus non-cool core thermodynamic profiles. More recently Mantz et al. (2017) analysed the redshift evolution of the iron content of 254 massive clusters based on Chandra and Suzaku observations also finding essentially no evolution especially in the outer parts of the ICM. In summary it seems that observationally the outer ICM metallicity profiles do not strongly evolve with redshift, while the inner parts show some evolution, although not all observational results are fully consistent with each other. In this section we will study the enrichment evolution of our cluster sample in more detail.
To begin we first show qualitatively the evolution of the metallicity maps of the most massive cluster of TNG300 in Fig. 10. At each redshift the extent of the projection is 3 × r500,crit, and as in the maps in the previous section the circles denote r500,crit. The mass-weighted average metallicity (resolution corrected) in solar units within this radius changes from 0.127 (0.203) at z = 0, to 0.124 (0.198) at z = 0.5, and to 0.138 (0.221) at z = 1. Therefore, the overall metal budget is not strongly evolving with redshift for more than ∼ 7 Gyr. However, local changes in the actual distributions of the metals in the ICM are visible in the maps. We note that this is in agreement with the observational findings of McDonald et al. (2016); i.e. the average metal content does not evolve below z ∼ 1.5. At some redshifts, for example at z = 0.5, one can also see how mergers bring in metal rich gas. The second and third row of Fig. 10 show the corresponding time evolution of the silicon over iron and oxygen over iron abundance ratios.

Redshift evolution of metallicities
To quantify the redshift evolution of the ICM enrichment in more detail, we present in Fig. 11 the median metallicity, iron, silicon and oxygen profiles of the cluster samples of TNG300 at different redshifts (z = 0.5, z = 1, z = 2) compared to the present-day profiles by showing their ratios. Here, and in the following, we trace back the clusters in our sample and look at their progenitors at different redshifts. We find that the outer parts of the ICM do not evolve much since z = 2, and the profiles agree well with those at z = 0. Therefore the outer parts (∼ r500,crit) of the ICM are not only very universal among different clusters, reflected by the small scatter at larger radii found in the previous section, but they are also invariant over a long period of time. Hence, enrichment does not seem to occur in those outer parts and any evolution is restricted to the more central parts of the cluster. In fact, the outer cluster abundance profiles beyond 0.4 × r500,crit are strikingly constant over time with nearly no evolution over the last ∼ 10 Gyr. This finding is in agreement with other theoretical studies (e.g., Biffi et al. 2017) and is also consistent with observational results as described above. The median profiles in Fig. 11 also demonstrate that there is some evolution of the profiles occurring at the centers of the clusters, where they tend to become flatter towards later times. The central amplitudes decrease by nearly a factor of 2 at the cluster centers from z = 2 to z = 0. We note that Martizzi et al. (2016) found a similar effect in their simulations, but this significant evolution seems to be absent in observations (e.g., . It is unclear whether this points to a systematic uncertainty in the observations, or whether the central enrichment in simulations, and related to this potentially the AGN feedback, are inconsistent with the observed central metallicity evolution at cluster centers. Besides studying median abundance profiles we can also aver- . Redshift evolution of the mass-weighted average metallicities for clusters of TNG300 with resolution correction. We present the redshift evolution within two different radial cuts: inner (r < 0.1 × r 500,crit ), and total (r < r 500,crit ). For the total cut the median over the cluster sample does not evolve significantly with redshift. For the inner regions we find that the metallicity slightly decreases towards lower redshifts consistent with the findings in Fig. 11. The constancy of the total metallicity is in agreement with recent observational findings as indicated by the gray median line and shaded region taken from McDonald et al. (2016) also measured within r 500,crit . We also compare the evolution of the central average iron abundance to the best-fit model of . The bands for observations and simulation data show the 1σ scatter. This is much larger for the observational data since the cluster to cluster variation is dominated by measurement uncertainties. Based on our simulations we find that the true cluster to cluster variation is significantly smaller, ∼ 0.03 Fe , for all redshifts for the total cut (r < r 500,crit ).
age the metallicity over the full cluster within r500,crit as we have done in the previous section. We can then calculate medians and percentiles of all these averages to inspect the redshift evolution of the median metallicity of our cluster sample. This procedure mimics the observational analysis of, for example, McDonald et al. (2016) and . We present the resulting evolution in Fig. 12. At each redshift we neglect cluster progenitors below M500,crit < 2 × 10 14 M when calculating the medians to be con- In Fig. 12 we present the median evolution of the metallicity averages within two different radial regions within the cluster: inner (r < 0.1 × r500,crit), and total (r < r500,crit). Our findings for the redshift evolution of the cluster sample median is in qualitative agreement with the results of McDonald et al. (2016) who found for 153 galaxy clusters between 0 < z < 1.5 only a very weak redshift dependence for the average metallicities within r500,crit. Specifi-cally, McDonald et al. (2016) fit their median metallicity redshift dependence within r500,crit with a linear model Z(z ) = a + bz where z = z − 0.6 is chosen to roughly minimise the covariance between a and b for the subsample of clusters observed with Chandra (z ∼ 0.6). For the global metallicity within r500,crit they find a = 0.23 ± 0.01 Z and b = −0.06 ± 0.04 Z . Within the error bars this is consistent with nearly no evolution. We note that this overall redshift dependence is weaker than the one reported in Ettori et al. (2015). However, the analysis of McDonald et al. (2016) presented the first study of the ICM metallicity evolution in a sample of galaxy clusters selected via the Sunyaev-Zel'dovich effect resulting in a selection that represents a mass-limited, redshiftindependent sample that is relatively unbiased with respect to the presence or lack of a cool core. This sample is therefore significantly less biased than previously employed cluster samples to explore the redshift evolution of the average metallicity of clusters. Correspondingly, the X-ray follow-up of galaxy clusters at 0 < z < 1.2 provides therefore excellent constraints on the global metallicity evolution and the radius-dependent evolution over this redshift range. Mantz et al. (2017)  There is a slight redshift evolution in the simulation data, which might be caused by the limited number of progenitors towards higher redshifts. We also note that the redshift evolution in the central region (r < 0.1 × r500,crit) is stronger, showing an increase of about ∼ 20% from z = 0 to z ∼ 1, in agreement with the evolution of the radial profiles presented in Fig. 11. However, this trend is not seen in the observational data although the average enrichment level agrees with observations. Fig. 12 also includes bands showing the 1σ scatter over the simulated cluster samples. This scatter is much larger for the observational data since the cluster to cluster variation is dominated by measurement uncertainties. Based on our simulations we find that the true physical scatter for r < r500,crit, ∼ 0.03 Fe at all redshifts, is significantly smaller than the observational variation.

Origin of metals
The metal and element budget presented above originates from three different stellar population sources: AGB stars, SNIa, and SNcc. We can trace back the origin of metals and iron at different radii within the ICM using the metal tagging technique of AREPO. The resulting fractional mass profiles with the different metal and iron sources are presented in Fig. 13 for TNG300. At all radii the contribution from SNIa to the overall metal budget is not very large and metal production is largely dominated by SNcc. This is expected and also in agreement with observations that find based on the relative metal abundances in the outskirts of clusters that the enrichment should have been dominated by SNcc with only a small metal fraction coming from SNIa (Simionescu et al. 2015). Furthermore, the contribution of AGB stars to the iron budget is negligible and SNIa dominate the iron production as expected. We also note that the scatter between the different source fractions is small, indicated by the contours. Especially, in the outer parts we find a very low dispersion within our cluster sample for the fractional metal and iron contributions from different sources. There is, therefore, also a strong uniformity among all 370 clusters regarding the origin of metals and iron at most radii within the ICM. The three panels also demonstrate that fractional source profiles for metal and iron  Figure 13. Median radial profiles of the production sources of iron and metals split in three different production channels: mass fraction resulting from AGB, SNIa, SNcc for TNG300 at different times (z = 0, 1, 2). Shaded regions show the 1σ scatter over the cluster sample. The iron budget is at all radii dominated by production through SNIa, whereas SNcc dominate the overall metal budget in the ICM. The observationally inferred SNIa mass fraction is based on the SNIa versus SNcc rates presented in Ezer et al. (2017). Those are nearly consistent with the simulation results within the observational error bars. We find at all redshifts that the metal contribution from SNcc is slowly decreasing towards the cluster center, whereas the SNIa is very slowly increasing. As a consequence the silicon over iron and oxygen over iron ratios are also decreasing towards the center as demonstrated in Fig. 7.  Figure 14. The fractional mass contribution of SNIa and SNcc to the total iron budget within r 500,crit at different redshifts. At early times SNcc contribute most of the iron until z ∼ 2.5 since SNcc enrichment occurs on shorter timescales than SNIa enrichment. At z ∼ 2.5 the SNIa contribution begins to dominate the overall iron budget. The transition is rather quick such that below z ∼ 2 the contribution of SNIa versus SNcc is quite constant in agreement with the profiles presented in Fig. 13 below z = 2.
are rather flat in the outskirts of the clusters, and similar for all redshifts back to z ∼ 2. Towards the center the iron contribution from SNIa turns slightly up while at the same time the contribution from SNcc goes down. We see similar trends for the total metal content: in the outer parts of the cluster SNcc dominate the overall metal budget. However, towards the center the contribution of AGB stars is increasing. Biffi et al. (2017) performed a similar analysis for their zoomin cluster sample. However they inspected only a small number of clusters for the origin of metals, so a statistically meaningful comparison with our findings is not possible. Nevertheless, the ratio of metals produced from SNcc over metals produced by SNIa decreases towards the cluster center in our simulations, whereas the ratio of metals from SNcc to metals from SNIa slightly increases at smaller radii for the cluster sample in Biffi et al. (2017). This directly impacts the radial profiles of the silicon over iron and oxygen over iron ratios since silicon and oxygen are produced by SNcc, whereas iron is produced by SNIa. We therefore expect that those abundance ratios should follow the trend seen for the metal sources in Fig. 13. Indeed, as demonstrated in the previous section, we find that the oxygen over iron and silicon over iron ratios slightly decrease towards the cluster centers, which is a consequence of the enrichment pattern seen in Fig. 13. On the contrary, Biffi et al. (2017) find slightly increasing silicon over iron and oxygen over iron ratios towards the cluster centers in agreement with their different enrichment history reflected by the different contributions from SNIa and SNcc as a function of radius in their analysis. In both cases, and as expected, those abundance ratios are good tracers to differentiate between different enrichment histories. However, it also seems that the different galaxy formation models and enrichment schemes of the two simulations lead to different trends for the detailed enrichment history within the ICM.
Recent observations favor rather flat abundance ratios for silicon over iron and oxygen over iron (e.g., Mernier et al. 2017) as shown in Fig. 7. These observational abundance ratios can also be fit with a combination of SNIa and SNcc yield models to quantify the radial contribution of SNIa and SNcc products. Mernier et al. (2017) shows that this leads to quite constant contributions from  In total more than ∼ 80% of all metals within the virial radius at z = 0 have been accreted from the IGM; i.e. has an ex-situ origin. The shaded region indicate the variation of this fraction within the cluster sample. This scatter is small towards z = 0 demonstrating that most clusters in the sample have accreted a similar amount of metals. We also show the metal mass within r 200,crit at each redshift, and the average metallicity of accreted gas. The latter does not evolve much since z ∼ 2, i.e. despite the fact that significant metal mass is accreted afterward, this gas has a rather uniform metallicity which causes the outer metallicity profiles to be nearly constant as a function of time. The black line shows the growth of the median total mass of the cluster sample. Upper right panel: Median of the average metallicities of gas and stars, and star formation rate density as a function of redshift measured in spheres with a comoving radius of 5 Mpc around the cluster progenitors. The proto-cluster environment is enriched at the time when the star formation rate peaks. This enriched material is then later accreted onto the clusters. The enrichment of stars is delayed compared to the gas since those stars have to form from enriched gas first. Lower panels: Metallicity maps for the most massive cluster of TNG300 at different redshifts. The extent of the maps are 10 Mpc comoving, i.e. covering the region over which averaged quantities are measured in the upper right panel. Circles indicate r 500,crit .
SNIa and SNcc, which seems to be slightly inconsistent with the trends seen in current simulations. We also include some observationally derived SNIa fraction in the top panel of Fig. 13 based on Ezer et al. (2017). They find for Abell 3112 that the ratio of the number of SNIa explosions to SNcc has a uniform distribution at a level of 12 − 16% out to the virial radius. We convert this ratio to an overall metal mass production ratio to be consistent with the results presented in Fig. 13. Our simulations predict a lower SNIa metal contribution compared to these observations although the outer data points are marginally consistent with the simulation predictions within the error bars. In general, this nearly constant SNIa contribution found in observations suggests that the ICM in cluster outskirts has been enriched by metals at an early stage before the cluster itself was formed during the period of intense star formation activity. Overall this finding agrees with our results, although our SNIa contribution seems to be slightly smaller. We can also use the metal tagging technique to understand the origin of iron in clusters as a function of redshift as presented in Fig. 12. In Fig. 14 we show the median fraction of iron mass originating from SNIa and SNcc as a function of time for the cluster sample of TNG300. As expected the SNcc contribution dominates at early times due to the short SNcc timescales. At later times the contribution from SNIa dominate since SNIa operate on much longer timescales. The figure also reveals that we expect an equal contribution from SNIa and SNcc at around z ∼ 2.5. At earlier times SNcc dominate the iron budget, whereas SNIa take over at later times. The transition is rather quick such that below z ∼ 2 the contribution of SNIa versus SNcc is quite constant in agreement with the profiles presented in Fig. 13.
To get more insights into the origin of ICM metals and the reason for the uniformity and time-invariance of the corresponding distributions, we track the amount of metal mass accreted onto the cluster from beyond its virial radius using the tracer particles ). In the upper left panel of Fig. 15 we present the cumulative fraction of accreted metal mass as a function of redshift for the cluster sample of TNG300. We find that more than ∼ 80% of all metals within r200,crit at z = 0 have been accreted during the evolution of the cluster sample. This number is quite universal and robust across the cluster sample as can be seen by the small scatter of the accretion fraction at low redshifts. We note that we only count metals entering the virial radius in the form of gas; i.e. metal enriched gas returned by infalling stars is not accounted for. Therefore, the ∼ 80% value should be interpreted as a lower limit. We also include in the upper left panel of Fig. 15 the normalised total mass of metals within r200,crit at each redshift. At z ∼ 2 the cluster has assembled only about 10% of its final metal mass, which is consistent with the fact that at that time only a small fraction of metals have been accreted. The black line shows the growth of the median total mass of the cluster sample.
Although these results seem to be consistent with each other they seem to contradict the results we have presented in Fig. 11, where we have demonstrated that the outer metallicity profiles are time-invariant since around z ∼ 2. This inconsistency can be understood by looking at the average metallicity of the gas that is accreted onto the clusters, which we also show in in the upper left panel of Fig. 15 (dashed line). Below z ∼ 2 the average metallicity of accreted gas is nearly constant and close to Z/Z ∼ 0.1, which is consistent with the constant metallicity profiles in the outer cluster regions discussed above. The accreted gas has therefore been enriched before accretion, so that the average metallicity within the cluster stays rather constant in the outer parts in agreement with the findings presented above. We note that the metallicity of the ac-creted gas varies only between 0.05 Z Z 0.2 Z since z ∼ 2 over the cluster sample. Therefore, the metallicity of accreted gas is also quite universal for our cluster sample.
In the upper right panel of Fig. 15 we quantify the early enrichment of the proto-cluster environment in more detail. Here we present the median mass-averaged gas and stellar metallicity within a comoving radius of 5 Mpc around the cluster progenitors at various redshifts along with the comoving star formation rate density in the volume. We exclude the gas within r200,crit at each redshift. The median star formation rate peaks around z ∼ 3 and leads to significant enrichment of the nearby gas below z ∼ 2. This median metallicity is consistent with the accreted metallicity seen in the upper left panel. At the same time the average stellar metallicity is also increasing towards lower redshifts. However, this occurs only after some delay since the stars have to form from the enriched gas first. In the lower panels we present metallicity maps for the most massive halo of TNG300 at different redshifts. The extent of the maps is 10 Mpc comoving, i.e. covering the region over which averaged quantities are measured in the upper right panel. Those maps also demonstrate the high average metallicity in the protocluster environment already at early times. This enriched material is then accreted towards lower redshift and is the reason for the time-invariant metallicity profiles in the outer parts of the clusters.
We conclude that our galaxy formation model is able to reproduce broadly also the high redshift metallicity measurements. Combined with the low-redshift results presented in the previous section, we arrive at a self-consistent enrichment history which seems to favor an early enrichment. However, we also note that the detailed contribution of SNIa and SNcc to the metal budget and consequently the abundance ratios are slightly inconsistent with most recent observational data. A natural extension of the analysis presented here would be a study to investigate the exact origin of the metals within the proto-cluster environment using the tracer particles.

CONCLUSIONS
We have presented a first analysis of the metal content of galaxy clusters of the new IllustrisTNG simulations Marinacci et al. 2017;Naiman et al. 2017;Pillepich et al. 2017a;Nelson et al. 2017). IllustrisTNG is the follow-up project of Illustris (Vogelsberger et al. 2014a,b;Genel et al. 2014;Sijacki et al. 2015), and consists of three different main simulations TNG50, TNG100, TNG300 each at three different resolution levels. The side length of TNG100 (2 × 1820 3 resolution elements) is 110.72 Mpc with a mass resolution of 7.46 × 10 6 M and 1.39 × 10 6 M for the DM and baryonic components, respectively. TNG300 (2 × 2500 3 resolution elements) has a side length of 302.63 Mpc and a mass resolution of 5.88 × 10 7 M and 1.1 × 10 7 M for the DM and baryonic components, respectively. TNG50 (2 × 2160 3 resolution elements) has a side length of 51.67 Mpc and a mass resolution more than an order of magnitude higher than TNG100.
Here we have focused on the analysis of two cluster samples (M500,crit > 10 13.75 M ) of the TNG100 (20 clusters) and TNG300 (370 clusters) simulations, which contain sizeable cluster populations. We find that the low-redshift metallicity profiles in those simulations agree with X-ray inferred observed abundance profiles. Furthermore, the simulations also reproduce high-redshift evolutionary trends of the ICM metal content. Our main findings are: • The radial median metallicity profiles agree with observational X-ray data, and the clusters in the different simulations show similar profiles across a wide mass range. The median metallicity, iron, silicon and oxygen profiles follow shallow negative gradients with very little scatter in the outskirts of the ICM (σZ/Z ∼ 15% relative scatter). This scatter increases to σZ/Z ∼ 40% towards the centers of the clusters (∼ 0.01 × r500,crit). The central metallicity reaches ∼ 0.4 − 1.0 Z , and drops to ∼ 0.2 Z at larger radii (∼ r500,crit). The central logarithmic slope of the metallicity profile is close to −0.1 at one percent of r500,crit and −0.5 at r500,crit. The azimuthal iron scatter around r500,crit has roughly the same amplitude as the iron scatter found through Suzaku observations of the Perseus cluster (Werner et al. 2013). We note that we have focused here only on total metallicity, iron, silicon, and oxygen abundance profiles and distributions. However, our galaxy model also traces other elements beyond those. We present profiles for the other abundance ratios in Fig. 16. All profiles are rather flat with a slight positive gradient. The highest ratio is predicted for neon, whereas nitrogen has the lowest values at all radii.
• Cool core clusters in our simulations, classified by their central being less than 30 keV cm −2 (Hudson et al. 2010), exhibit radial metallicity distributions, which are steeper and peaked towards the cluster centers compared to those of non-cool core clusters in agreement with observational data (Ettori et al. 2015). Cool core clusters approach 0.8 Z towards their centers whereas non-cool core clusters have central metallicities below 0.5 Z . The metallicity profiles of both overlap beyond 0.15 × r500,crit in agreement with observational data (Ettori et al. 2015).
• Silicon over iron and oxygen over iron ratios show small positive gradients, and are therefore in slight tension with most recent observations. The ratios increase by about 50% from 10 −2 × r500,crit to r500,crit. The increasing ratios of these two elements are consistent with the radial distribution of the production sites of metals in our simulations. Towards the center the relative contribution of SNIa compared to SNcc is slightly rising resulting in a decrease of the X/Fe ratios towards the centers.
• The distribution of the mass-weighted averaged metallicities within r500,crit can be described by narrow Gaussians with a dispersion of 0.027 Z . We find a small anti-correlation for the temperature-iron abundance relation with a logarithmic slope, ∼ −0.26, in agreement with the compiled analysis of Yates et al. (2017). However, the overall normalisation of the temperature-iron abundance relation is about 0.2 dex too low in the simulation compared to observational data.
• The outskirts, r ∼ r500,crit, of the metallicity, iron, silicon and oxygen profiles do not evolve since redshift ∼ 2. Therefore, the outer cluster abundance profiles beyond 0.4 × r500,crit are constant over time with nearly no evolution over the last ∼ 10 Gyr. The inner profiles, r < 0.1 × r500,crit, are flattening towards lower redshift and decrease in their amplitude. The central amplitudes decrease by nearly a factor of 2 at the cluster centers from z = 2 to z = 0. This effect is not seen in observations (e.g., . It is unclear whether this points to a systematic uncertainty in the observations, or whether the central enrichment in simulations, and related to this the AGN feedback, are inconsistent with the observed behaviour at cluster centers. • The averaged metallicities within the ICM (r < r500,crit) do not evolve with redshift since z ∼ 1 in agreement with recent ob- . Based on our simulations we find that the true physical scatter of iron abundances at various redshifts is significantly smaller (∼ 0.03 Fe ) than the observational variation, which is dominated by measurement uncertainties. The iron is produced at early times by SNcc (beyond z ∼ 2.5) and at late times through SNIa below z ∼ 2.5. The transition of the two regimes is rather quick such that the contributions of SNcc and SNIa do not evolve significantly below z ∼ 2. The average metallicity within 0.1 × r500,crit is evolving from z = 0 to z = 1 with an increase of about ∼ 20% towards z ∼ 1. We compare the evolution of this central average iron abundance to the best-fit model of , and find that the average enrichment level agrees with the simulation data, but we see a different redshift evolution.
• We find no central metallicity drops (Mernier et al. 2017) in the clusters of the IllustrisTNG simulations. However, we observe those for clusters of the Illustris simulation where a different radio mode AGN model has been employed. This difference might point towards bubble-driven uplifting of metals through AGN feedback in Illustris, but not in IllustrisTNG since the latter simulation suite uses a kinetic AGN feedback model at the cluster center not supporting such a mechanism at the same strength. For Illustris the median metallicity profiles peak at 0.03 × r500,crit reaching ∼ 0.7 Z and drop to ∼ 0.6 Z at 10 −2 × r500,crit. This is comparable to the central oxygen drop seen in Mernier et al. (2017), but larger than the observer iron drop therein.
• At least ∼ 80% of the metals within the virial radius (r200,crit) are accreted. This fraction is universal across the cluster sample with very little scatter. The accreted gas is enriched to ∼ 0.1 Z at around z ∼ 2 slightly after the peak of star formation in the protocluster environment. At z ∼ 2 the cluster has assembled only about 10% of its final metal mass, which is consistent with the fact that at that time only a small fraction of metals have been accreted yet. The proto-cluster environment for the clusters is already enriched to > 0.1 Z below z ∼ 2 briefly after the peak of the local star formation. This proto-cluster gas is then accreted onto the cluster and leads to the constancy of the outer metallicity profiles from z ∼ 2 to z = 0.
• The detailed shapes of the metallicity profiles and metal distri-butions are converged. The normalisation however depends on the degree of convergence of stellar mass. Any relative ratios derived from the simulation, e.g., abundance ratios profiles, are therefore converged. Furthermore, low resolution results can be extrapolated to higher resolutions through a scaling with the stellar mass ratio of the two resolutions. Specifically, the metal content of clusters of TNG300-1 can be scaled to TNG100-1 resolution by applying a correction factor of 1.6 which follows from the stellar mass ratio between the two simulations.
We conclude that the predictions of IllustrisTNG are generally in good agreement with observational findings at low and high redshifts. Overall we find that the abundance profiles of clusters and the enrichment exhibit a strikingly uniform behaviour over our large cluster sample. This is true for the radial metallicity profiles, and also for the metal accretion history. Besides this uniformity, cluster average metallicities and metal profiles are also remarkably constant as a function of time. The TNG300 cluster sample with 370 objects demonstrate these characteristics in a convincing way through the large sample size, which is currently unparalleled.
Further progress in understanding the origin and evolution of metals in clusters has to come from two fronts: observationally and theoretically. Observationally, McDonald et al. (2016) and  presented the so far most detailed accounting of the enrichment history of the ICM. A more refined analysis will require next-generation observatories, such as XARM, Athena and X-ray Surveyer, combined with samples of clusters at z ∼ 2 which will be available with the next generation of Sunyaev-Zel'dovich experiments. Significant progress in measuring metallicity profile will also come with next generation X-ray observatories. The larger effective area of future missions will also allow studies of the outskirts of clusters in more detail. This combined with increased spectral resolution will lead to new insights into the distribution of metals in the ICM of clusters. From the theory side there still seems to be no perfect consensus on reproducing all details of observed abundance trends among different simulations. The detailed metal distribution, especially at cluster centers, also depends sensitively on details of the galaxy formation model, like the AGN feedback implementation, as we have demonstrated. Ultimately, more detailed galaxy formation models should be employed going beyond the coarse sub-resolution models used in simulations like IllustrisTNG. Those models should include more refined schemes for stellar winds and AGN feedback. Both of them affect the gas enrichment and should therefore be modelled more faithfully at resolutions beyond those of IllustrisTNG. In fact, even at the TNG100-1 resolution more detailed AGN feedback models are applicable (e.g., Weinberger et al. 2017b), which might lead to differences in the central metal distribution. Higher resolution cluster simulations should therefore definitely explore those new models and contrast the results with those of coarser models like those of Illustris and IllustrisTNG. Last, also other physical processes like thermal conduction (Kannan et al. 2017) or dust formation (McKinnon et al. 2016 affect the metallicity profiles in and around galaxies and should therefore be considered in future models.