Distribution and Evolution of Metals in the Magneticum simulations

Metals are ideal tracers of the baryonic cycle within halos. Their composition is a fossil record connecting the evolution of the various stellar components of galaxies to the interaction with the environment by in- and outflows. The Magneticum simulations allow to study halos across a large range of masses and environments, from massive galaxy clusters containing hundreds of galaxies down to isolated field galaxies. They include a detailed treatment of the chemo-energetic feedback from the stellar component and its evolution as well as feedback from the evolution of supermassive black holes. Following the detailed evolution of various metal species and their relative composition due to continuing enrichment of the IGM and ICM by SNIa, SNII and AGB winds of the evolving stellar population reveals the complex interplay of local star formation processes, mixing, global baryonic flows, secular galactic evolution and environmental processes. We present results from the Magneticum simulations on the chemical properties of simulated galaxies and galaxy clusters, carefully comparing them to observations. We show that the simulations already reach a very high level of realism within their complex descriptions of the chemo-energetic feedback, successfully reproducing a large number of observed properties and scaling relations. Our simulated galaxies clearly indicate that there are no strong secondary parameters (like star formation rates at fixed redshift) driving the scatter in these scaling relations. The remaining differences clearly point to detailed physical processes which have to be included into future simulations.


Introduction
In cosmological hydrodynamical simulations, star formation is typically treated based on a sub-grid model. Star particles are formed when the gas locally exceeds a certain density threshold, when it is Jeans unstable, and the local gas flow converges at a gas particle. Each star particle represents a simple stellar population which is characterized by an initial mass function (IMF).
In addition, each star particle emits a feedback that is mimicking the winds which are launched by stars in reality by averaging over the whole stellar population included in each simulated stellar particle. In most models, the considered winds compose of three main stellar sources: Supernovae Type II (SNII), Supernovae Type Ia (SNIa), and asymptotic giant branch stars (AGBs). These three sources are considered to be the most important drivers of observed stellar outflows and thus the main contributers to metal enrichment of the interstellar medium: Stars with masses larger than 8 M are believed to end their lives through a so called core collapse (SNII), releasing typically an energy of 10 51 erg per supernova together with their chemical imprint into the surrounding gas. These supernovae usually occur shortly after the formation of the simulated star particle, as the most massive stars only live for a short time and thus this part of the stellar population of the star particle is emitting first. SNIa on the other hand are believed to arise from thermonuclear explosions of white dwarfs within binary stellar systems. These supernovae are delayed relatively to the SNII and occur later in a star particle's live, as the white dwarf in the according binary system has to first be formed and second needs to accrete matter from its companion until it reaches the mass threshold for the onset of thermonuclear burning. Therefore, a detailed modelling of the evolution of the stellar population has to be included in the simulations to follow the chemo-energetic imprint of the SNIa on the surrounding gas. AGB stars contribute emitting strong winds, exhibiting strong mass losses during their life and are importantly contributing to the nucleosynthesis of heavy elements.
To include the effects of these sources into the sub-grid models of simulated stellar particles raises the need to integrate a set of complicated equations describing the evolution of a simple stellar population. Such set of integral equations then allows to compute at each time the rate at which the current AGB stars pollute their environment by stellar winds and the rate at which the SNIa and SNII are exploding. Such they allow to properly treat the chemo-energetic imprint on the surrounding inter-galactic medium (IGM) and the inter-cluster medium (ICM). In the following we will give a short, schematic description of such calculations based on the description presented in Dolag et al., (2017) [1] (for a more detailed review, see Matteucci 2003 [2] and Borgani et al. 2008 [3], and references therein). Then we will introduce the cosmological hydrodynamical simulation set used in this work, and present a comparison of the metal enrichment in galaxies and clusters of galaxies caused by this model feedback with observations. Finally, we will summarize and discuss the results of this comparison in the light of model improvements needed for the future.

Chemical Enrichment
To describe the continuous enrichment of the IGM and ICM through winds from SNIa, SNII and AGB stars in the evolving stellar population of each star particle, several ingredients are needed. In the following we will shortly present the ingredients used in many state-of-the-art cosmological simulations of galaxy formation.

IMF
One of the most important quantities in models of chemical evolution is the IMF. It directly determines the relative ratio between SNII and SNIa, and therefore the relative abundance of α-elements and Fe-peak elements. The shape of the IMF also determines the ratio between low-mass long-living and massive short-living stars. This ratio directly affects the amount of energy released by the SNe as well as the present luminosity of galaxies which is dominated by low mass stars and the (metal) mass-locking in the stellar phase.
The IMF φ(m) describes the number of stars of a given mass per unit logarithmic mass interval. Historically, a commonly used form is the Salpeter IMF (Salpeter 1955 [4]) which follows a single power-law with an index of x = 1.35. However, different expressions of the IMF have been proposed more recently in order to model a flattening in the low-mass regime of the stellar mass function that is currently favoured by a number of observations. Among the newer, often used models is the Chabrier IMF (Chabrier 2003 [5]), which has a continuously changing slope and is more top heavy than the Salpeter IMF: However, the question of whether there is a global IMF or if the IMF is changing with galaxy mass, morphology or cosmological time, and which IMF has to be chosen is still an unsolved problem and a matter of heavy debate.

Lifetime functions
To model the evolution of a simple stellar population, a detailed knowledge of the lifetimes of stars with different masses is required. Different choices for the mass-dependence of the lifetime function have been proposed in the literature (e.g. Padovani and Matteucci (1993) [6], Maeder and Meynet (1989) [7], Chiappini et al., (1997) [8]), where the latest one reads: (2)

Stellar Yields
The ejected masses of the different metal species i produced by a star of mass m are called stellar yields p Z i (m,Z). These yields depend on the metallicity Z with which the star originally formed, and on the type of outflow ejected from the star. Therefore, detailed predictions for the main three sources of enrichment are needed, namely for the mass loss of AGB stars, of SNII, and of SNIa. Up to date, such predictions still suffer from significant uncertainties, mainly due to the still poorly understood mass loss through stellar winds in stellar evolution models, which depends on multiple additional physical processes.
For the mass loss through AGB stars, the most recent predictions can be found in Karakas (2007) [9]. Predictions for the mass loss from massive stars driving SNII are presented by Nomoto, Kobayashi & Tominaga (2013) [10]. The most complete table for SNIa up to date is presented by Thielmann (2003) [11].

Modelling the Enrichment Process
As summarized by Dolag (2017) [1], the assumption of a generic star formation history represented by an arbitrary function of time ψ(t) allows to compute the rates for the different contributions in form of a set of integral equations as shown in the following (for more details, see also Matteucci 2003 [2], Borgani et al. 2008 [3], and references therein). This formalism can be individually applied to the large number of particles representing the continuous star-formation process within cosmological simulations. Every star particle here represents a stellar population born in a single burst. The combination of all the stellar component within the simulated galaxies results then in a model which describes the legacy of the detailed star-formation history of any simulated galaxy.

Type Ia supernovae
SNIa occur in binary systems with masses in a mass range of 0.8-8 M . Let m B be the total mass of the binary system and m 2 the mass of the secondary companion. We can now use f (µ) as the distribution function of binary systems with µ = m 2 /m B and define A as the fraction of stars in binary systems that are progenitors of SNIa. Therefore, A has to be given or obtained by a model. Constructing such detailed models for SNIa progenitors is particular difficult, see for example Greggio & Renzini (1983) [12] or   [13]. Based on such kind of models, typical value for A are inferred to be in the range of 0.05 and 0.1, based on comparisons of chemical enrichment models with observed iron metallicities within galaxy clusters (e.g. Matteucci & Gibson 1995 [14]) or within the solar neighbourhood (e.g. Matteucci & Greggio 1986 [15]). Within the current simulations we are using a value of A = 0.1. With these ingredients and the mass dependent lifetime functions τ(m), we can model the rate of SNIa as where M Bm and M BM are the smallest and largest values, respectively, that are allowed for the progenitor binary mass m B . Then, the integral over m B runs in the range between M B,inf and M B,sup , which represent the minimum and the maximum value of the total mass of the binary system that is allowed to explode at the time t. These values in general are functions of M Bm , M BM , and m 2 (t), which in turn depend on the star formation history Ψ(t). In simulations, where the individual stellar particles are in commonly modelled as an impulsive star formation event, ψ(t) can therefore be approximated with a Dirac δ-function. The sum of all stellar particles and their individual formation time then represent the complex star-formation history of the galaxies within the simulation.

Supernova Type II, Low-, and Intermediate-mass stars
Computing the rates of SNII, low-mass stars (LMS), and intermediate-mass stars (IMS) is conceptually simpler than calculating the rates of SNIa, since they are purely driven by the lifetime function τ(m) convolved with the star formation history ψ(t) and multiplied by the IMF φ(m = τ −1 (t)). Since ψ(t) is a delta-function for the simple stellar populations used in simulations, the SNII, LMS and IMS rates read where m(t) is the mass of the star that dies at time t. For AGB rates, the above expression must be multiplied by a factor of (1 − A) if the mass m(t) falls in the range of masses which is relevant for the secondary stars of SNIa binary systems.

The equations of chemical enrichment
In order to compute the total metal release from the simple stellar population, we have to fold the above rates with the yields p SNIa|SNII|AGB Z i (m, Z) from SNIa, SNII and AGB stars for a given element i for stars born with initial metallicity Z i , and compute the evolution of the density ρ i (t) for each element i at each time t. As shown in Borgani et al., (2008) [3], this readṡ where M L and M U are the minimum and maximum mass of a star in the simple stellar population, respectively. Commonly adopted choices for these limiting masses are M L 0.1 M and M U 100 M .
In the above equation, the first line describes the locking of metals in new born stars through the currently ongoing star formation ψ(t), which for the assumed sub-grid model case vanishes as ψ(t) is

The Magneticum Simulations
The Magneticum simulation set covers a huge dynamical range, from very large cosmological volumes as shown in Fig. 1, which can be used for statistical studies of clusters and voids (e.g. Bocquet et al. 2016 [17], Pollina et al. 2017 [18]), to very high resolution simulations of smaller cosmological volumes which allow a morphological classification and a detailed analysis of galaxies and their properties (e.g. Teklu et al. 2015 [19] and Remus et al. 2017 [20]). Table 1 lists the detailed properties like size and stellar mass-resolution of the different simuations. These simulations treat the metal-dependent radiative cooling, heating from a uniform time-dependent ultraviolet background, star formation and the chemo-energetic evolution of the stellar population as traced by SNIa, SNII  and AGB stars with the associated feedback processes and stellar evolution details as described before. They also include the formation and evolution of super-massive black holes and the associated quasar and radio-mode feedback processes.

Metallicities from Magneticum in Comparison to Observations
As previously shown from re-simulations of massive galaxy clusters, the observed radial profiles of Iron within the ICM can be well reproduced in simulations: Biffi et al., (2017) [22] demonstrated that especially at high redshifts the implemented AGN feedback is the main driver in enhancing the metal enrichment in the ICM at large cluster-centric distances to the observed level.
The Magneticum simulations now allow such investigations across a much larger range of halo masses. For this study we use a large simulation volume (Box2 hr) with enough resolution to resolve mean properties of galaxies (and AGNs), resulting in the same resolution as used in Biffi et al., (2017) [22] and Hirschmann et al., (2014) [21]. For galaxies, the star formation efficiency (SFE) is a measure of the SFR per unit of gas mass. To evaluate whether the SF activity in the simulations match the observational data is a key step to proceed towards reproducing the observed mass-metallicity-relation (MZR). In Fig. 2, the evolution of the observed SFE up to a redshift of z ≈ 2 is shown in comparison to the evolution of the SFE of star forming galaxies in Magneticum. As can clearly be seen, the simulations are in excellent agreement with the observed trends up to z ≈ 2.

Galaxy Clusters: ICM Metallicities
X-ray observations of the intra cluster medium allow for a detailed study of the distribution of different metal species within the ICM via their line emissions within the X-ray band. Here, the composition of the individual metal species allows to interpret their abundances as an imprint of the relative contributions to the enrichment from SNIa and SNII, hence keeping a record on when and where these metals are injected into the ICM (e.g. De Plaa et al., 2007 [24]). The lower left panel of Fig. 3 shows a comparison of such observational data for various individual galaxy clusters with  simulations, clearly demonstrating the ability of the Magneticum simulations to reproduce the correct absolute iron abundance within the ICM self-consistently.
Here the scatter of the absolute iron abundance within the simulated clusters shows a similar spread (factor of two) as the observed values with a very mild trend of increasing iron abundance for low mass (e.g. lower ICM temperature) clusters. In addition the simulations also reproduce broadly the chemical composition footprint of various elements species, as shown in the various panels. This strongly indicates that the simulations predict the correct ratio between the contributions to the metal enrichment from SNIa and SNII and its interplay with the AGN feedback. The simulations also predict typically less spread in the composition of the metals for the indicidual clusters than they show in their absolute iron metalicity.

Galaxies: Gas Metallicities
To estimate the gas phase metallicity of galaxies, it is important consider that observationally the measurements are obtained only from star forming regions. Therefore, it can be misleading to only calculate the mean metallicity of all gas particles inside a simulated galaxy. Thus, after selecting star forming galaxies (see Teklu et al. 2017 [25]) we can either calculate their mean metallicity by averaging over all particles which are currently star-forming, or alternatively calculate the mean metallicity of the new-born stars. We tested both methods and found that the latter gives slightly better results due to the fact that it is difficult to catch the metallicity of the gas phase in the moment of star formation within simulations given the large timespan between the simulation outputs, while the young stellar population freezes the record of the metallicity of the gas from which it was formed. This leads to a MZR which is in good agreement with observations, as shown in the left panel of Fig. 4, where we use the predicted Oxygen abundances from our calculations to be consistent with observations (Sanders et al. 2015 [26] at z = 2.3 and Bresolin et al. 2016 [27] at z = 0). Interestingly, even the overall evolution of the MZR is well captured by the simulations and matches the observations, as demonstrated in the right panel of Fig. 4. However, when calculating gas-phase metallicity gradients within the galaxies, the simulations predict a steeper profile than the observations. Given that the prediction of the mean metallicities is in good agreement with observations this indicates that the simulations either still lack the resolution to properly describe the mixing of the enriched gas within galaxies, or that on these scales the diffusion of metals might have to be modelled more explicitly. Nevertheless, we clearly showed that including detailed modelling of the stellar population combined with current AGN feedback models significantly improve the predicted ICM and IGM metallicities and are needed to successfully reproduce various aspects of the observations.

Galaxies: Stellar Metallicities
In general, we assume that the mean metallicity of a stellar particle in the simulation represents the mean stellar metallicity of the stellar population represented by the stellar particle, neglecting any self-enrichment within stars. To be consistent with the observations, we based all calculations on the Iron abundance as predicted by the simulations. For this part of the study we used a smaller cosmological volume with a higher resolution, as this allows a classification of the galaxies due to their morphological type (see Teklu et al. 2015 [19]) for more details on the classification and this particular simulation). This higher resolution also allows for a more detailed resolution of the metallicity gradients with radius for individual galaxies.
Although the obtained mean metallicities of our galaxies are close to observational results, the stellar MZR obtained from the simulations is somewhat shallower that the one obtained from CALIFA observations by Gonzalez Delgado et al. (2014) [28], as shown in the left panel of Fig. 5. This again indicates that the treatment of the mixing between accreted, more pristine material from outside the galaxy and the enriched gas within the galaxies is not fully captured yet. At this point it is unclear whether that will be resolved by further enhancing the resolution of simulations or whether explicit diffusion of metals has to be taken into account.
For a proper comparison of the radial metallicity gradients between simulations and observational data, it is necessary to calculate the effective radius R eff for each galaxy. This is done by selecting all stellar particles within ten percent of the virial radius and then inferring the according half-mass radius, which we associate to the effective radius R eff . This method has already been used to demonstrate that the Magneticum simulations successfully reproduce the mass-radius relation of galaxies for different morphological types (e.g. Remus et al. 2015 [29]) and at different redshifts (e.g. Remus et al. 2017 [20]).
Interestingly, the radial stellar metallicity gradients obtained from the simulations very well match the CALIFA observations, as shown in the right panel of Fig. 5. Note also that the increasing spread towards larger radius is an intrinsic, point-by-point spread within the individual galaxies and agrees well with the behaviour measured for individual galaxies by Pastorello et al. (2014) [30]. However, while the simulations successfully reproduce the observed Iron abundances and radial gradients, they still produce a flat ratio of Oxygen (or Magnesium or Silicium) over Iron, in contrast to the observations, and here the sub-grid model clearly needs to be advanced.

Discussion and Conclusion
Metals are measured in all phases of the baryonic universe, starting from the ICM in galaxy clusters down to gas and stars in individual galaxies. As the radial distribution of different types of metals encode a record of the star formation process over the whole time of a galaxies' evolution, comparing the chemical footprint as seen in simulations with observations reflects an excellent test on the reliability of the numerical sub-grid models included in modern state-of-the-art simulations of galaxy formation in a cosmological context in reproducing the complex processed enrolled in galaxy formation.
We demonstrated that the Magneticum simulations are able to reproduce a large variety of observational findings over a large range of halo masses, indicating that the star-formation and feedback processes included start to be highly realistic. At galaxy cluster scale the simulations show an excellent agreement with both, the absolute value as well as the cluster by cluster variation of the the iron abundance when comparing to x-ray observations and also reproduce broadly the chemical composition of the ICM. On galaxy scales the tight mass-metallicity-relations found for our simulated galaxies indicate that there are no strong secondary parameters (like star formation rates at fixed redshift) driving the scatter in these relations. The remaining differences between the observed properties and the simulation results indicate, however, that the incorporation of physical processes like diffusion and mixing have to be improved within the next generation of simulations to successfully reproduce the detailed distribution of metals on individually resolved galaxy scales.