Galaxy metallicity scaling relations in the EAGLE simulations

TL;DR

This study uses the EAGLE simulations to explore how galaxy metallicities relate to properties like stellar mass, gas fraction, and star formation, revealing key correlations and their dependence on feedback processes.

Contribution

It provides a detailed analysis of metallicity scaling relations in simulated galaxies and links these to feedback mechanisms, improving understanding of galaxy evolution.

Findings

The $M_*-Z_{SF,gas}$ relation slope depends on resolution.

Simulated metallicity evolution increases by ~0.5 dex since z=3.

Gas fraction is the primary correlate with metallicity.

Abstract

We quantify the correlations between gas-phase and stellar metallicities and global properties of galaxies, such as stellar mass, halo mass, age and gas fraction, in the Evolution and Assembly of GaLaxies and their Environments (EAGLE) suite of cosmological hydrodynamical simulations. The slope of the correlation between stellar mass and metallicity of star-forming (SF) gas ($M_*-Z_{\rm SF,gas}$ relation) depends somewhat on resolution, with the higher-resolution run reproducing a steeper slope. This simulation predicts a non-zero metallicity evolution, increasing by $\approx 0.5$ dex at $\sim 10^9 {\rm M}_{\odot}$ since $z = 3$. The simulated relation between stellar mass, metallicity and star formation rate at $z \lesssim 5$ agrees remarkably well with the observed fundamental metallicity relation. At $M_* \lesssim 10^{10.3} {\rm M}_{\odot}$ and fixed stellar mass, higher metallicities are associated with lower specific star formation rates, lower gas fractions and older stellar populations. On the other hand, at higher $M_*$, there is a hint of an inversion of the dependence of metallicity on these parameters. The fundamental parameter that best correlates with the metal content, in the simulations, is the gas fraction. The simulated gas fraction-metallicity relation exhibits small scatter and does not evolve significantly since $z = 3$. In order to better understand the origin of these correlations, we analyse a set of lower resolution simulations in which feedback parameters are varied. We find that the slope of the simulated $M_*-Z_{\rm SF,gas}$ relation is mostly determined by stellar feedback at low stellar masses ($M_* \lesssim 10^{10} {\rm M}_{\odot}$), and at high masses ($M_* \gtrsim 10^{10} {\rm M}_{\odot}$) by the feedback from active galactic nuclei.

Tables

Table 1. Table 1: Parameters that are varied in the simulations. Columns list: simulations identifiers, the side length of the volume ( L 𝐿 L ) and the particle number per species (i.e. gas, DM) per dimension ( N 𝑁 N ), the power-law slope of the polytropic equation of state of SF gas ( γ eos subscript 𝛾 eos \gamma_{\rm eos} ), the power-law index of the star formation law ( n 𝑛 n ), the asymptotic maximum ( f th , max subscript 𝑓 th max f_{\rm th,max} ) and minimum ( f th , min subscript 𝑓 th min f_{\rm th,min} ) values of f th subscript 𝑓 th f_{\rm th} (Equation 3 ), the parameters that control the characteristic density and the power-law slope of the density dependence of the energy feedback from star formation ( n H , 0 subscript 𝑛 H 0 n_{\rm H,0} and n n subscript 𝑛 n n_{\rm n} , respectively), the subgrid accretion disc viscosity parameter ( C visc subscript 𝐶 visc C_{\rm visc} ), and the temperature increment of stochastic AGN heating ( Δ ​ T AGN Δ subscript 𝑇 AGN \Delta T_{\rm AGN} ). The upper section comprises models that have been calibrated to reproduce the z = 0.1 𝑧 0.1 z=0.1 GSMF and the lower section comprises models featuring single-parameter variations of Ref. Numbers in bold indicate variations with respect to the reference model (Ref). This table has been adapted from Table 1 1 1 in Crain et al . ( 2015 ) for the simulations used in this work.

Identifier	Side length L	$N$	${\gamma }_{\mathrm{eos}}$	$n$	$f_{\mathrm{th},\max }$	$f_{\mathrm{th},\min }$	$n_{H,0}$	$n_n$	$C_{\mathrm{visc}}/2\pi$	$\Delta T_{\mathrm{AGN}}$
	[cMpc]			[${\mathrm{cm}}^{-3}$]			[${\mathrm{cm}}^{-3}$]			${\log }_{10}$ [K]
Calibrated models
Recal-L025N0752	25	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.25$	$\mathrm{𝟏}/\ln \mathrm{𝟏𝟎}$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$9.0$
Ref-L025N0752	25	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
Ref-L025N0376	25	376	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
Ref-L050N0752	50	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
Ref-L100N1504	100	1504	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
FBconst-L050N0752	50	752	$4/3$	$1.4$	$1.0$	$1.0$	$-$	$-$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$8.5$
Ref. variations
eos1-L025N0376	25	376	$\mathrm{𝟏}$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
eos5/3-L025N0376	25	376	$\mathrm{𝟓}/\mathrm{𝟑}$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
KSLow-L025N0376	25	376	$4/3$	$1.0$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
KSHi-L025N0376	25	376	$4/3$	$1.7$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
WeakFB-L025N0376	25	376	$4/3$	$1.4$	$1.5$	$0.15$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
StrongFB-L025N0376	25	376	$4/3$	$1.4$	$6.0$	$0.6$	$0.67$	$2/\ln 10$	${10}^0$	$8.5$
NOAGN-L050N0752	50	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	$-$	$-$
AGNdT8-L050N0752	50	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$8.0$
AGNdT9-L050N0752	50	752	$4/3$	$1.4$	$3.0$	$0.3$	$0.67$	$2/\ln 10$	${10}^0$	$9.0$

Equations

$$n_{\rm H}^{*}=10^{-1}{\rm cm}^{-3}\big{(}\frac{Z}{0.002}\big{)}^{-0.64},$$

$${\dot{m}}_{\ast }=m_{\mathrm{g}}A(1{\mathrm{M}}_{\text{sun}}{\mathrm{pc}}^{-2}{)}^{-n}(\frac{\gamma }{G}{f}_{\mathrm{g}}P{)}^{(n-1)/2},$$

$$f_{\rm th}=f_{\rm th,min}+\frac{f_{\rm th,max}-f_{\rm th,min}}{1+\big{(}\frac{Z}{0.1Z_{\sun}}\big{)}^{n_{Z}}\big{(}\frac{n_{\rm H,birth}}{n_{\rm H,0}}\big{)}^{-n_{n}}},$$

$$Z=y\ln (f_{\mathrm{gas}}^{-1}),$$

$$y_{\mathrm{eff}}=\frac{y_0}{1+(M_0/M_{\mathrm{bar}}{)}^{0.57}},$$

$$Z=y(1+f{)}^{-1}\ln [1+(1+f)(f_{\mathrm{SF},\mathrm{gas}}^{-1}-1)].$$

Full text

††pagerange: Galaxy metallicity scaling relations in the EAGLE simulations–LABEL:lastpage††pubyear: 2017

1 Introduction

The metallicity properties of galaxies encode crucial information about the different physical processes (e.g. star formation, infall and outflows of gas, etc.) that drive their evolution (e.g. Tinsley, 1980; Davé, Finlator & Oppenheimer, 2011; Finlator, 2017). In this context, the determination of scaling relations between the metallicity of galaxies and other key properties such as their stellar masses ($M_{*}$) or gas fractions is a matter of great interest and the subject of an on-going debate in the community as it can help to constrain models of structure formation.

In the local Universe, there is a well-known correlation between stellar mass and metallicity, such that more massive galaxies are more metal-enriched (e.g. Lequeux et al., 1979; Tremonti et al., 2004). This mass-metallicity relation (MZR) has also been investigated at higher $z$ but with an apparent offset towards lower metallicities with respect to the local MZR (e.g. Savaglio et al., 2005; Erb et al., 2006; Maiolino et al., 2008). However, a comparison between observational works at different redshifts is not straightforward because of selection biases, aperture effects and the use of different metallicity indicators (e.g. Steidel et al., 2014). Thus, the detailed dependencies of the MZR and their level of evolution are still debated (e.g. Kewley & Ellison, 2008; Telford et al., 2016).

Recently, several authors have reported secondary dependences of metallicity at a given stellar mass, arguing that the MZR is just the projection onto two-dimensions (the mass-metallicity plane) of more fundamental relations contained in higher-dimensional parameter space (e.g. the space defined by mass, metallicity and star formation rate). In particular, Ellison et al. (2008) reported that galaxies with higher specific star formation rates (SFRs) or larger half-light radii show systematically lower gas-phase metallicities than systems with similar $M_{*}$ but lower specific SFRs or smaller sizes. More recently, Mannucci et al. (2010) have reported the existence of a ‘fundamental metallicity relation’ (FMR) between $M_{*}$, metallicity and star formation rate (SFR) that exhibits little scatter and does not show significant evolution, at least below $z\approx 2.5$ (see also Lara-López et al., 2010). According to the FMR, systems with higher SFRs tend to have lower metallicities at a given stellar mass, consistently with the behaviour reported by Ellison et al. (2008). Thus, part of the observed evolution of the MZR might be due to the fact that surveys at high $z$ tend to select systems with elevated SFRs.

More recently, Bothwell et al. (2013) suggested that the FMR might be a consequence of a more fundamental correlation between $M_{*}$, metallicity and gas fraction: systems with higher gas fractions tend to exhibit lower metallicities at a fixed mass. These studies have been extended at different $z$ by many authors in recent years (e.g. Cresci et al., 2012; Hunt et al., 2012; Henry et al., 2013; Kacprzak et al., 2016; Lara-López, López-Sánchez & Hopkins, 2013; Stott et al., 2013; Cullen et al., 2014; Maier et al., 2014; Nakajima & Ouchi, 2014; Zahid et al., 2014a; Bothwell et al., 2016a) but, as in the case of the MZR, the large uncertainties and diverse observational techniques involved in observational works prevent convergence towards a clear determination of these fundamental metallicity relations of galaxies, as discussed by Telford et al. (2016).

The tightness of the observed FMR seems to depend strongly on the adopted abundance determination method and might be affected by systematic errors such as those associated with stellar mass estimates, aperture effects, etc. Telford et al. (2016), for example, applied new abundance diagnostics to data from the Sloan Digital Sky Survey (SDSS), obtaining an anti-correlation between metallicity and SFR, at a given mass, which is weaker than that found by Mannucci et al. (2010), by 30-55%. Andrews & Martini (2013), on the other hand, reported a stronger anti-correlation between metallicity and SFR, at a fixed mass, than Mannucci et al. (2010) for nearby galaxies. According to Yates, Kauffmann & Guo (2012), the trend of increasing metallicity with decreasing SFR inverts at high stellar masses while Salim et al. (2014) suggested that this apparent turnover might be an artefact due to the signal-to-noise cuts imposed on the observational sample. Finally, other authors claimed that there is no significant dependence of the observed MZR on the SFR (e.g. Hughes et al., 2013; Sánchez et al., 2013, 2017).

From the theoretical point of view, different works have tried to address the origin and evolution of metallicity scaling relations (e.g. Tissera, De Rossi & Scannapieco, 2005; Davé, Finlator & Oppenheimer, 2012; Dayal, Ferrara & Dunlop, 2013; Yates, Kauffmann & Guo, 2012; Lilly et al., 2013; Romeo Velonà et al., 2013; Yates et al., 2013; Vogelsberger et al., 2014; Ma et al., 2016; Genel, 2016; Weinberg, Andrews & Freudenburg, 2017). These models and simulations reproduce qualitatively the observed trends but they show discrepancies regarding the exact value of the slope and level of evolution of the predicted relations. Supernova (SN)-driven outflows have often been invoked as a key ingredient for establishing a MZR (Larson, 1974; Tremonti et al., 2004; Dalcanton, 2007; Kobayashi, Springel & White, 2007). Given the shallower potential wells of low-mass galaxies, metal-enriched material can be more efficiently ejected from these systems, keeping their metallicities low. Note, however, that the mass-loading of a wind, $\dot{M}_{w}/\dot{M}_{\star}$, is not necessarily the same as the metal mass loading of a wind, $\dot{M}_{Z}/\dot{M}_{\star}$ (e.g. Mac Low & Ferrara, 1999; Creasey, Theuns & Bower, 2015), which makes the problem more complex. Besides, less efficient star formation in low-mass galaxies could also cause the lower chemical enrichment of smaller galaxies and explain the origin of the MZR (e.g. Brooks et al., 2007; Mouhcine et al., 2008; Calura et al., 2009). The evolution of the MZR could also be affected by the different SF histories associated with galaxies with different morphologies (see e.g. Calura et al., 2009). In addition, the infall of metal-poor gas onto the outer parts of galaxies or inflows triggered by mergers with other systems could play an important role (e.g. Köppen & Edmunds, 1999; Dalcanton, Yoachim & Bernstein, 2004; Finlator & Davé, 2008; Davé, Finlator & Oppenheimer, 2011). de Rossi, Tissera & Scannapieco (2007), for example, have shown that a correlation between mass and metallicity can arise naturally in a hierarchical scenario solely as a consequence of the regulation of SF by merger events. However, the lack of SN-feedback model prevented de Rossi, Tissera & Scannapieco (2007) from reproducing the observed slope of the MZR due to overcooling. According to Köppen, Weidner & Kroupa (2007), the MZR might simply reflect variations in the stellar initial mass function (IMF).

De Rossi et al. (2015a, b) and De Rossi et al. (2016) have studied the origin of the MZR and FMR by using the Galaxies-Intergalactic Medium Interaction Calculation (GIMIC, Crain et al. 2009) suite of cosmological hydrodynamical simulations. These authors found that the star-forming (SF) gas and stellar components of simulated galaxies follow local mass-metallicity relations very similar to those observed but with less scatter. The simulated relations seem to be driven mainly by infall of metal-poor gas as well as by the efficient action of SN feedback. However, the GIMIC simulations do not predict the observed level of evolution of the MZR, because of the old average stellar ages ($\sim 10$ Gyr) of simulated galaxies. In addition, the GIMIC simulations do not reproduce the observed flattening of the MZR at the high-mass end. De Rossi et al. (2015b) claimed that the latter issue is probably related to the lack of AGN feedback in GIMIC.

More recently, Schaye et al. (2015) showed that the Evolution and Assembly of GaLaxies and their Environments (EAGLE, Schaye et al. 2015) simulations are able to reproduce a correlation between stellar mass and gas-phase metallicity at $z=0$, which agrees well with observed data (Tremonti et al., 2004; Zahid et al., 2014a) for the high-resolution version of those simulations. EAGLE high-resolution simulations have also been found to predict an evolution of the MZR consistent with the observational trend reported by Zahid et al. (2013a) (e.g. Guo et al., 2016). In addition, by analysing intermediate-resolution EAGLE runs, Lagos et al. (2016) found that metallicity can be robustly determined from neutral gas fractions, or from $M_{*}$ and SFR. In particular, the strength of the stellar feedback implemented in EAGLE has an important impact on the $z=0.1$ simulated MZR (Crain et al., 2015). Furthermore, according to the results of Segers et al. (2016b), the AGN model included in EAGLE yields a relation between $M_{*}$ and stellar-$\alpha$-element-to-iron ratio ($[{\rm\alpha/Fe}]_{*}$) consistent with observations of massive ($M_{*}>10^{10.5}{\rm M}_{\sun}$) early type galaxies. Simulations can therefore play an important role in examining the relative importance of different physical processes that drive the evolution of correlations. In addition, they are valuable tools in examining observational biases by allowing the comparison of intrinsic correlations to those inferred from simulations after applying observed selections to mock observables.

In this article, we extend previous works by analysing in detail the evolution of metallicity of galaxies as a function of mass and redshift using the EAGLE suite of cosmological simulations. We focus on the analysis of the high-resolution version of the simulations which implement the so-called recalibrated model (see below). We show that the improved EAGLE subgrid prescriptions lead to a better description of the evolution of the MZR than the GIMIC simulations, preserving many key features of the observed relation that were also reproduced by GIMIC. In particular, the AGN feedback model in EAGLE yields a better description of the global metal enrichment of massive galaxies, generating the observed turn-over of the MZR at the high-mass end, as was also shown by Segers et al. (2016a).

The plan of the paper is as follows. The simulation and the sample selection are described in Section 2. The simulated mass-metallicity relation obtained from SF gas abundances is presented in Section 3 while that derived from stellar abundances are discussed in Section 4. In Section 5, we compare observed and simulated effective yields. The scatter of the simulated MZR and its dependence on secondary parameters, such as gas fraction, SFR and mass-weighted stellar age, is discussed in Section 6. In Section 7, we explore different sets of simulations corresponding to different resolutions and models in order to reveal the processes that determine the main features of the simulated metallicity scaling relations. Finally, our conclusions are summarised in Section 8.

2 The EAGLE simulations

The EAGLE suite111See http://eagle.strw.leidenuniv.nl and http://www.eaglesim.org/ for different data products, images and movies. In addition, a database with many integrated properties of EAGLE galaxies (McAlpine et al., 2016) and particle data (The EAGLE team, 2017) is publicly available. of cosmological simulations (Schaye et al., 2015; Crain et al., 2015) constitutes a set of different hydrodynamical simulations run with different resolutions, box sizes and subgrid physics models. The simulations were performed by using a modified version of the GADGET-3 Smoothed Particle Hydrodynamics (SPH) code (last described by Springel 2005). The modifications to the SPH implementation are collectively referred to as ANARCHY (Dalla Vecchia, in prep.; see also Schaller et al. 2015).

The EAGLE simulations track the joint evolution of dark matter and baryons within periodic comoving volumes of side-length up to 100 comoving megaparsec (cMpc) from $z=127$ to $z=0$. A reference model has been implemented, for which the subgrid parameters associated with energy feedback were calibrated to obtain good agreement with the $z=0.1$ galaxy stellar mass function (GSMF), whilst also reproducing the observed sizes of present-day disk galaxies. Furlong et al. (2015) show that the simulated GSMF also broadly reproduces the data up to $z\approx 7$. In addition to the reference model, other variations of subgrid parameters have been explored as discussed in detail in Crain et al. (2015).

To distinguish more easily the different runs within the EAGLE suite of simulations, the name of a given simulation includes a suffix that indicates the box length in comoving megaparsec (e.g. L100) and the cube root of the initial number of particles per species (e.g. N1504). Simulations with the same subgrid model as the primary one (the reference model) are denoted with the prefix ”Ref-” (e.g. Ref-L100N1504). As discussed in Schaye et al. (2015), for higher resolution simulations, a ”Recal-” model has been implemented in addition to the reference one. The former model uses subgrid parameters that have been recalibrated following similar procedures to those applied to the reference run to improve the fit to the $z\sim 0$ GSMF when working with the high-resolution simulations (Recal-L025N0752).

Within the EAGLE suite, the Recal-L025N0752 simulations have been found to reproduce the observed trends for the MZR (Schaye et al. 2015, see also below), because of its higher feedback efficiency. Thus, our analysis will be focused on this run. We will also perform comparisons with other simulations carried out within a fixed simulated volume (”L025” runs) but using different resolutions and model parameters to assess their effects on the obtained trends. Variations of AGN feedback parameters have only been tested within comoving volumes of side-lengths of 50 cMpc (”L050” runs); we will employ this set of simulations to analyse the impact of AGN feedback. In Table 1, we present the set of simulations analysed in this work, indicating their main subgrid parameters, which are described below.

2.1 Cosmological parameters and subgrid implementation

A flat $\Lambda$CDM cosmology is assumed with cosmological parameters consistent with those inferred by the Planck Collaboration (2014): ${\Omega}_{\Lambda}=0.693$, ${\Omega}_{\rm m}=0.307$, ${\Omega}_{\rm b}=0.04825$, ${\sigma}_{8}=0.8288$, $h=0.6777$, $n_{s}=0.9611$ and $Y=0.248$ where ${\Omega}_{\rm m}$, ${\Omega}_{\Lambda}$ and ${\Omega}_{\rm b}$ are the average densities of matter, dark energy and baryonic matter in units of the critical density at $z=0$, ${\sigma}_{8}$ is the square root of the linear variance of the matter distribution when smoothed with a top-hat filter of radius 8 $h^{-1}$ cMpc, $H_{0}\equiv h\ 100\ {\rm km\ s^{-1}\ Mpc^{-1}}$ is the Hubble parameter, $n_{s}$ is the scalar power-law index of the power spectrum of primordial adiabatic perturbations, and $Y$ is the primordial abundance of helium. Below, we briefly describe the main characteristics of EAGLE subgrid physics; more details can be found in Schaye et al. (2015).

Radiative cooling and heating rates are computed on an element-by-element basis for gas in ionization equilibrium in the presence of a Haardt & Madau (2001) ionizing UV/X-Ray background and the Cosmic Microwave Background. The total metallicity variable, $Z$, and the 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) that are important for the radiative cooling at $T>10^{4}$ K (Wiersma, Schaye & Smith, 2009) are tracked individually.

Star formation is implemented stochastically following Schaye & Dalla Vecchia (2008), but including a metallicity-dependent density threshold $n_{\rm H}^{*}$, as:

[TABLE]

that yields the Kennicutt-Schmidt relation. It exhibits a dependence on gas metallicity $Z$ that captures the transition from the warm, atomic to the cold, molecular gas phase (Schaye, 2004). A temperature floor $T_{\rm eos}({\rho}_{\rm g})$ is applied, which is associated with the equation of state $P_{\rm eos}\propto{\rho}_{\rm g}^{{\gamma}_{\rm eos}}$ (${\gamma}_{\rm eos}=4/3$, for standard runs), normalised to $T_{\rm eos}=8\times 10^{3}$ K at $n_{\rm H}=10^{-1}{\rm cm}^{-3}$, a typical temperature for the warm interstellar medium (ISM) (e.g. Richings, Schaye & Oppenheimer, 2014). In this way, cold dense gas is prevented from artificial fragmentation due to a lack of resolution. When gas particles reach densities $n_{\rm H}>n_{\rm H}^{*}$ and $\log_{10}(T/{\rm K})<\log_{10}T_{\rm eos}/{\rm K}+0.5$ , they are eligible for star formation and are assigned an SFR, ${\dot{m}}_{*}$ (Schaye & Dalla Vecchia, 2008):

[TABLE]

where $m_{\rm g}$ is the gas particle mass, $\gamma=5/3$ is the ratio of specific heats, $G$ is the gravitational constant, $f_{\rm g}$ is the mass fraction in gas and $P$ is the total pressure. For the EAGLE simulations used in this work, $f_{\rm g}=1$. The parameters $A=1.515\times 10^{-4}\ {\rm M}_{\sun}\ {\rm yr}^{-1}\ {\rm kpc}^{-2}$ and $n=1.4$ are obtained directly from the observed Kennicutt-Schmidt relation (Kennicutt, 1998), when scaled to the Chabrier (2003) IMF. Throughout this paper, when performing comparisons with observational results that assume another IMF, data are converted to a Chabrier IMF for consistency.

The chemical enrichment model follows the prescriptions of Wiersma et al. (2009). As mentioned, the simulations track the stellar mass losses of 11 elements associated with three stellar evolutionary channels: ($i)$ stellar winds and core-collapse (type II) supernovae resulting from massive stars ($M>6{\rm M}_{\odot}$), ($ii)$ type Ia supernovae assumed to result from catastrophic mass transfer in close binary stars, and ($iii$) winds from asymptotic giant branch (AGB) stars. Stellar evolutionary tracks and yields that depend on the initial metal abundance are implemented. For type II supernovae, yields from Portinari, Chiosi & Bressan (1998) were used because they consider mass loss from massive stars. In the case of AGB stars, yields of Marigo (2001) were implemented as they constitute a self consistent set with yields from Portinari, Chiosi & Bressan (1998). For type Ia supernovae, Wiersma et al. (2009) used the last version of the standard ”W7” model (Thielemann et al., 2003) (the reader is referred to Wiersma et al. 2009 for more details about yield choices and implementation). As discussed in detail in Wiersma et al. (2009), nucleosynthetic yields are uncertain by factors of a few and the abundance evolution is sensitive to the particular choice of yield tables.

Stochastic thermal feedback from star formation is applied in the EAGLE simulations. The feedback model is described by Dalla Vecchia & Schaye (2012) and is based on a stochastic selection of neighbouring gas particles that are heated by a temperature increment of $10^{7.5}$K. Taking into account the local metallicity and gas density, a fraction $f_{\rm th}$ of energy from core-collapse supernovae is injected into the ISM 30 Myr after the birth of a stellar population (Schaye et al., 2015; Crain et al., 2015). In the model, $f_{\rm th}$ is given by:

[TABLE]

where $n_{\rm H,birth}$ is the density inherited by the star particle from its parent gas particle, $Z$ is the metallicity, $f_{\rm th,min}$ and $f_{\rm th,max}$ are the asymptotic values of $f_{\rm th}$ while $n_{Z}$, $n_{n}$ and $n_{\rm H,0}$ are free parameters. It is assumed that $n_{Z}=n_{n}$. The parameters $n_{\rm H,0}$ and $n_{n}$ were chosen to reproduce the present-day GSMF and galaxy sizes.

When the halo mass of a system increases above $10^{10}h^{-1}{\rm M}_{\sun}$, seed black holes (BHs) of mass $10^{5}h^{-1}{\rm M}_{\sun}$ are placed inside them following Springel, Di Matteo & Hernquist (2005). BHs grow by subsequent gas accretion events222 When the subgrid BH mass becomes higher that its host particle mass, the BH can stochastically accrete neighbouring gas particles (see Schaye et al. 2015 for more details). and mergers at a rate computed according to the modified Bondi-Hoyle accretion rate of Rosas-Guevara et al. (2015) and Schaye et al. (2015). To regulate the Bondi rate in high-circulation flows, a viscosity parameter $C_{\rm visc}$ is introduced. AGN feedback is implemented thermally and stochastically similarly to energy feedback from star formation. Particles surrounding the BH are chosen randomly and heated by a temperature $\Delta T_{\rm AGN}$. Increasing $\Delta T_{\rm AGN}$ leads to more energetic individual feedback events, generally resulting in smaller radiative losses in the ISM. In addition, larger values of $\Delta T_{\rm AGN}$ generate a more intermittent feedback process. We note that a single feedback mode is included in EAGLE as the current implementation naturally yields an AGN feedback that mimics ’quasar-mode’ and ’radio-mode’ feedbacks at high and low accretion rates, respectively (Rosas-Guevara et al., 2016).

2.2 Set of studied simulations

Table 1 summarises the different EAGLE simulations analysed in this work and their main parameters.

As noted by Schaye et al. (2015), the correlation between $M_{*}$ and star forming (SF) gas metallicity ($Z_{\rm SF,gas}$) depends on resolution and the choice of feedback parameters. In particular, Recal-L025N0752 reproduces the slope and normalization of certain observational data sets better because of the increase of the feedback efficiency with respect to the reference model. The increase of the energy feedback from star formation required to match the observed GSMF, simultaneously decreases the metal content of the ISM of low-mass galaxies leading to a better agreement with observations. As discussed in Section 7 (see also Section 3), a decrease in resolution (which in our case also implies the use of the reference model and so, different feedback parameters) does not alter the main trends of the fundamental metallicity scaling relations but can moderately affect their detailed features.

Our analysis will be focused on the high-resolution simulation Recal-L025N0752 because it agrees better with the slope and normalization of certain observed metallicity scaling relations. For intermediate-resolution simulations, the predicted MZR is too flat at low stellar masses (see Fig. 1 and corresponding discussion in Sec. 3). Thus, unless otherwise specified, we will show results from the Recal-L025N0752 simulation in this work. In addition, for assessing resolution effects and comparison with the reference model at a fixed volume, we will use the simulations Ref-L025N0376 and Ref-L025N0752. To test the impact on our results of the slope of the $P-\rho$ relation imposed at high $\rho$, we will analyse simulations eos1-L025N0376 and eos53-L025N0376. The effects of changing the power-law index in the star formation law will be addressed by studying simulations KSLow-L025N0376 and KSHi-L025N0376. Simulations WeakFB-L025N0376 and StrongFB-L025N0376 will allow us to compare a weak and strong stellar feedback model, respectively. Simulation FBconst-L050N0752 will allow us to evaluate the effects of injecting into the ISM a fixed quantity of energy per unit stellar mass formed, independent of local conditions. We will compare the effects of varying the AGN feedback temperature, at a fixed volume and resolution, by contrasting simulations Ref-L050N0752, AGNdT8-L050N0752 and AGNdT9-L050N0752. We will also analyse the simulation NOAGN-L050N0752, for which the BH implementation is turned off. Thus, in the later simulations, BH gas accretion and AGN feedback are disabled entirely.

2.3 Identifying galaxies in EAGLE

Dark matter overdensities were identified by applying the “Friends-of-Friends” (FoF) method, assuming a linking length of 0.2 times the average inter-particle spacing (Davis et al., 1985). In the case of baryonic particles, they are assigned to the same FoF-group as their nearest dark matter neighbour. Self-bound substructures, which can contain both dark matter and baryons, are then identified by applying the SUBFIND algorithm (Springel, Yoshida & White, 2001; Dolag et al., 2009). A FoF halo can contain several SUBFIND sub-groups, or subhalos. We define the central galaxy as the subhalo with the lowest value of the gravitational potential while any remaining subhalos are classified as satellite galaxies. Unless otherwise specified, we include both central and satellite galaxies in our sample.

In order to mimic the aperture of instruments used for observations, we computed integrated quantities inside a given radius. Unless otherwise indicated, we follow Schaye et al. (2015) and report global properties of galaxies inside a sphere with radius 30 proper kpc (pkpc). In general, for the analysis carried out here, aperture effects do not change the general trends but can generate moderate variations in the slope and normalization of metallicity scaling relations (see the Appendix A) particularly for massive galaxies.

3 The $M_{*}-Z_{\rm SF,gas}$ relation

In this section, we study the correlation between stellar mass and gas-phase metallicity of simulated galaxies and compare our results with different observational findings. The mass-metallicity relation is typically inferred by measuring oxygen abundances in star-forming gas HII regions. Thus, we calculated global oxygen abundances (${\rm O/H}|_{\rm SF,gas}$) considering gas particles with ${\dot{m}}_{*}>0$ (see Eq. 2) at radius $r\leq 30$ kpc in a given subhalo. It is worth noting that, had we considered all gas particles (i.e. SF gas and non SF -NSF- gas) for these estimates, the gas metallicities would have decreased by $\approx 0.1-0.3$ dex (the exact values depend on mass and redshift). This is because the NSF gas is less metal-enriched than the SF gas (see below).

3.1 The local $M_{*}-Z_{\rm SF,gas}$ relation

Schaye et al. (2015) have compared the local $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation measured by Tremonti et al. (2004) and Zahid et al. (2014a) with results from the EAGLE intermediate- and high-resolution simulations, Ref-L100N1504 and Recal-L025N0752, respectively. The two sets of observed data are both based on SDSS data, but metallicities were obtained using different techniques. Both observed MZRs agree at $M_{*}\sim 10^{11}\ {\rm M}_{\sun}$ but the relation reported by Tremonti et al. (2004) is steeper at low masses. The intermediate-resolution simulation agrees with the flatter relation given by Zahid et al. (2014a) to better than 0.1 dex for $M_{*}>10^{9.5}\ {\rm M}_{\sun}$ but, at lower masses, this simulation does not reproduce the steep observed slopes obtained by Tremonti et al. (2004) or Zahid et al. (2014a). On the other hand, the high-resolution simulation predicts a steeper slope for the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation, in better agreement with Tremonti et al. (2004) and Zahid et al. (2014a) results. Thus, the observed MZR seems to be better reproduced by the simulation Recal-L025N0752 than by Ref-L100N1504.

The shape and normalization of the observed MZR are still a matter of extensive debate. The use of different metallicity indicators, different methods for estimating stellar masses, selection biases, aperture effects and dust can affect the slope and zero-point of the MZR. In Fig. 1, we extend the analysis by Schaye et al. (2015) studying the observed MZR and comparing it with the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation obtained from simulations with different resolutions and different simulation volumes. In the left panel, we compare the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation obtained from the EAGLE Recal-L025N0752 simulation with the SDSS best-fit relations derived by Kewley & Ellison (2008) by using a set of diverse metallicity calibrators: KK04 (Kobulnicky & Kewley, 2004), Z94 (Zaritsky, Kennicutt & Huchra, 1994), KD02 (Kewley & Dopita, 2002), M91 (McGaugh, 1991), T04 (Tremonti et al., 2004), D02 (Denicoló, Terlevich & Terlevich, 2002), PP04 O3N2, PP04 N2 (Pettini & Pagel, 2004), P01 (Pilyugin, 2001) and P05 (Pilyugin & Thuan, 2005). We see that the differences between observational results can reach $\approx 0.7$ dex in some cases. The slope of the observed MZR is also affected by the choice of the metallicity calibrator, with the steeper relation obtained by the T04 method and the shallower relation inferred when using the P05 technique.

In the right panel of Fig. 1, we compare results from simulations with observational data reported by different authors: Zahid et al. (2012) (median metallicities in bins of $M_{*}$ and standard errors on the mean), Sánchez et al. (2013) (metallicities for individual galaxies), Salim et al. (2014) (median metallicities and standard deviations in $M_{*}$ bins) and Hunt et al. (2016) (medians of metallicities with the 75% and 25% quantile levels). For the sake of clarity, only the maximum and median errors reported by Sánchez et al. (2013) are represented at the bottom right corner. Zahid et al. (2012) determined metallicities of galaxies in the SDSS by using the strong-line calibration KK04. Sánchez et al. (2013) reported data from the CALIFA survey, with metallicities derived from the strong-line PP04 O3N2-calibrator. In the case of Salim et al. (2014), they calculated metallicities from SDSS data by applying the strong-line technique of Mannucci et al. (2010). Data from Hunt et al. (2016) correspond to a compilation of different samples of galaxies at $z=0$, with metallicities calculated using the PP04 N2 strong-line method. As the normalization of the observed MZRs obtained by these authors is affected by the different metallicity calibrators used, we renormalized these observational results in order to match the median Recal-L025N0752 relation at $M_{*}\approx 10^{10.5}\ {\rm M}_{\sun}$. We added -0.10, 0.36, -0.02 and 0.27 dex to the data reported by Zahid et al. (2012), Sánchez et al. (2013), Salim et al. (2014) and Hunt et al. (2016), respectively. In this way, we avoid normalization issues and can focus on the comparison of the shapes of the relations.

In the case of simulations, we estimated median relations from mass bins containing more than 10 galaxies. The number of galaxies per bin is $N_{\rm bin}\ga 30$, for simulation Ref-L100N1504 and $N_{\rm bin}\ga 20$, for simulation Ref-L050N0752. For simulations Ref-L025N0376, Ref-L025N0752 and Recal-L025N0752, $N_{\rm bin}\ga 10$, considering the whole analysed mass range and $N_{\rm bin}\ga 20$, if we remove the highest mass bin.

Independent of the volume, intermediate-resolution simulations (Ref-L025N0376, Ref-L050N0752, Ref-L100N1504) yield similarly flat shapes for the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation, departing from observations. On the other hand, high-resolution simulations yield steeper MZR slopes more consistent with the observational trend. We can see a general good agreement between results from high-resolution simulations and observations by Zahid et al. (2012) and Salim et al. (2014). In particular, the recalibrated model predicts a MZR that reproduces encouragingly well the relation found by Salim et al. (2014), even without including the normalization adjustment (-0.02 dex). Generally, the steeper (shallower) slope obtained at low (high) masses for Recal-L025N0752 simulation improves the agreement with the observed trend.

We note that, because of their smaller volume ($25^{3}\ {\rm cMpc}^{3}$), higher resolution simulations cannot sample high-density environments. This issue leads to a dearth of galaxies, specially towards higher masses, and might affect the comparison with observations. However, intermediate resolution simulations run in different volumes show similar metallicity relations, suggesting a weak influence of environment on the obtained trends. By analysing GIMIC simulations, De Rossi et al. (2015b) also found that metallicity scaling relations are not significantly affected by the large-scale environment. Recent observations suggest a small dependence of the MZR on the environment (e.g. Wu et al., 2017), too.

Finally, it is worth recalling that the determination of abundances in the simulations are affected by the choice of the nucleosynthetic yields that are uncertain by a factor of a few (Wiersma et al., 2009). Thus, given the issues affecting observational results and the uncertainties in simulated yields, the comparison between models and observations should be taken with care.

3.2 Evolution of the $M_{*}-Z_{\rm SF,gas}$ relation

In this section, we analyse the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation as a function of $z$ and compare it to observations. For this study, the high-resolution simulation Recal-L025N0752 is used, which is the simulation run that best reproduces the slope and normalization of certain observed stellar mass-gas-phase metallicity relations at $z\approx 0$ (Schaye et al., 2015).

In Fig. 2, we show the evolution of the simulated $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation. In each of the four large panels, the simulated relation is shown at two redshifts close to the observed data, with lines denoting the median relations and error bars the 25th and 75th percentiles. As a reference, the area enclosed by the 25th and 75th percentiles corresponding to $z=0$ is reproduced in all panels as a grey shadow. We constructed the median relations considering only mass bins containing more than 10 galaxies ($N_{\rm bin}\ga 10$). We also show observational findings reported at similar redshifts to those associated with the simulations. In the following, we summarize briefly the different observational data used here (for more details, the reader is referred to the corresponding papers): de los Reyes et al. (2015) (medians and scatter at $z\approx 0.8$), Hunt et al. (2016) (medians with 75% and 25% quantile levels at $0.4<z\leq 0.7$, $0.9<z\leq 1.8$, $1.8<z\leq 2.8$ and $2.8<z\leq 3.8$), Ly et al. (2016) (medians with the 16th and 84th percentiles at $z=0.5-1.0$), Stott et al. (2013) (median values and standard errors at $z\approx 0.84-1.47$), Yabe et al. (2014) (medians and bootstrap errors at $z\approx 1.4$), Zahid et al. (2014b) (fitted metallicities and observational uncertainties at $z\approx 1.6$), Wuyts et al. (2016) (mean values and standard errors at $z\approx 0.9$ and $z\approx 2.3$), Sanders et al. (2015) (mean values and uncertainties at $z\approx 2.3$), Cullen et al. (2014) (fitted metallicities with uncertainties at $z\ga 2$), Troncoso et al. (2014) (individual measurements and their uncertainties at $z\approx 3.4$) and, Onodera et al. (2016) (results from stacking analysis at $z\approx 3.3$). Dotted lines in the large four panels indicate the fits to observations at different $z$ given by Maiolino et al. (2008). In the case of Maiolino et al. (2008) and Troncoso et al. (2014), we show the relations reported for masses estimated from Bruzual & Charlot (2003) templates.

The use of different samples, selection criteria and, mainly, the implementation of different metallicity calibrators yield discrepancies between different observations at similar $z$. In particular, as already noted for the local relation, the normalization of the MZR is still a matter of extensive debate in the literature. This problem is even worse at higher $z$ due to different selection biases and the difficulties in measuring metallicities of very distant galaxies. In particular, an unbiased comparison of the MZR at various redshifts is not guaranteed by simply using one diagnostic/calibrator. This would be only true if the implemented metallicity calibrator/diagnostic is equally valid for the HII regions of galaxies at different redshifts, and does not have any mass/SFR/redshift dependent biases. Locally calibrated metallicity diagnostics may not be appropriate for measuring metallicities at higher redshift (Steidel et al., 2014; Strom et al., 2017). Also, the N2 and $R_{23}$ diagnostics are known to become insensitive to metallicity at 12+log(O/H) $\sim$ 9.0 (e.g. Kewley & Dopita, 2002; Liang et al., 2006). Nevertheless, in order to mitigate differences in normalizations generated by the use of different metallicity calibrators (see Section 3.1) and ease the comparison of the MZR evolution, we renormalized observed relations taking into account the local MZR used by each author to compare their observed relations at $z>0$. Those local MZRs were obtained by some authors from previous works while other authors re-estimated the local MZRs using self-consistent methods to those applied at higher $z$ (for further details, see the corresponding papers). As each author compared the corresponding high-$z$ MZR with a local MZR derived from the same metallicity calibrator, the effects of the choice of the metallicity indicator on the level of evolution of the MZR should diminish. Thus, in Fig. 2, an offset is applied to observed relations at $z>0$ so that their associated local observed MZRs yields $12+{\log}_{10}({\rm O/H})=9.0$ at $M_{*}=10^{10.5}\ {\rm M_{\sun}}$ (consistent with $z=0$ results from EAGLE Recal-L025N0752).333 Offsets applied to the data: -0.02 dex (Maiolino et al., 2008; Troncoso et al., 2014; Onodera et al., 2016), +0.10 dex (de los Reyes et al., 2015), +0.20 dex (Ly et al., 2016), +0.25 dex (Sanders et al., 2015), +0.27 dex (Hunt et al., 2016), +0.30 dex (Cullen et al., 2014), +0.32 dex (Zahid et al., 2014b; Wuyts et al., 2016) and +0.35 dex (Stott et al., 2013; Yabe et al., 2014). In this way, the comparison between findings from the simulation and the level of evolution reported in observational works is more straightforward: if the shape of the MZR does not change significantly with $z$, different normalizations in Fig. 2 can be related to different levels of evolution in the simulation with respect to observations and the association is exact at $M_{*}\approx 10^{10.5}\ {\rm M}_{\sun}$.

We see that, at all $z$, the EAGLE simulation Recal-L025N0752 reproduces the observed trend of increasing ${\rm O/H}|_{\rm SF,gas}$ with increasing $M_{*}$ well, exhibiting a slope very closed to the observed one. Also, at a given $M_{*}$, metallicity tends to decrease towards higher redshifts, in agreement with the observed trend. The predicted level of evolution at $z\la 1.5$ agrees remarkably well with the observed one. At $z\approx 2$, the simulation predicts a MZR in good agreement with Zahid et al. (2014b) but other observations suggest a higher level of evolution, with the larger discrepancies towards higher $M_{*}$ (e.g. Hunt et al., 2016). At $z\approx 3$, simulations tend to produce more metal-rich galaxies than inferred by Maiolino et al. (2008) and Onodera et al. (2016). On the other hand, observational findings by Hunt et al. (2016) at $z\approx 3$ seem to yield values closer to those from the simulations, in particular, towards lower $M_{*}$.

In the small bottom panels of Fig. 2, we compare the evolution of the simulated $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation to that derived by Hunt et al. (2016). Their “MEGA” sample comprises 1000 galaxies with a common O/H calibration (PP04 N2) and spans almost two orders of magnitude in metallicity, a factor of $10^{6}$ in SFR, and a factor of $10^{5}$ in stellar mass. As explained above, we renormalized MEGA’s relations to $12+{\log}_{\rm 10}({\rm O/H})=9.0$ at $M_{*}\approx 10^{10.5}\ {\rm M}_{\sun}$ at $z=0$. In the small left panel, we compare simulation findings with multi-variable linear regressions on the MEGA dataset for 12+log(O/H) as a function of $M_{*}$ and $z$ (see Hunt et al., 2016, for details). We see that the simulation predicts an evolution of the MZR in good agreement with the fitted relations obtained for the MEGA sample at $z=0-3$. On the other hand, the small right panel shows that the MEGA median relations at different $z$ bins exhibit some discrepancies with simulations. At the low-mass end ($M_{*}\approx 10^{9.5}\ {\rm M}_{\sun}$), simulations predict the observed level of evolution at $z\approx 1-3$ but overpredict, by $\la 0.2$ dex, the evolution at $z\approx 0-1$ (observations predict a negligible evolution in this case). At the high-mass end ($M_{*}\approx 10^{10.5}\ {\rm M}_{\sun}$), simulations and observations predict a negligible evolution at $z=0-1$. Although simulations lack massive systems at high $z$, the extrapolation of simulation relations towards higher masses suggests a lower level of evolution compared to that obtained for the MEGA sample.

Ly et al. (2016) carried out the first systematic study of the evolution of the MZR to $z\approx 1$ using only the electron temperature ($T_{e}$) method. These authors reported good agreement between their findings and the evolution of the MZR obtained from the EAGLE high-resolution simulation Recal-L025N0752. Guo et al. (2016) also found a good agreement between the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation derived from this simulation and observations by Zahid et al. (2013a) at $z\la 1$. At $1\la z\la 2$, the latter authors obtained a slightly lower level of evolution (by $\approx 0.2$ dex) in simulations with respect to observations444 Note that Guo et al. (2016) computed the MZR considering all particles inside EAGLE galaxies while, in the present work, we only include particles within an aperture of 30 pkpc. (the reader is referred to Guo et al. 2016 for a comparison between EAGLE results and predictions from some semi-analytical models). On the other hand, as mentioned, the predicted metallicity evolution below $z\approx 3$ seems to be lower than what the data by Maiolino et al. (2008) suggest (e.g. $\approx 1$ dex at $M_{*}\approx 10^{9}\ {\rm M}_{\sun}$ and $\approx 0.8$ dex at $M_{*}\approx 10^{10}\ {\rm M}_{\sun}$). More recently, Onodera et al. (2016) have also reported a stronger evolution ($\approx 0.7$ dex) for the observed MZR below $z\approx 3.0-3.7$. As discussed in Section 1, because of the observed FMR, the lower abundances of observed galaxies at $z\ga 1$ might be partly explained by selection effects present in observational surveys that tend to be biased towards systems with higher SFRs at a given $M_{*}$ and, hence, lower metallicities. In Section 6, we show that EAGLE galaxies are also consistent with the existence of a FMR since $z\approx 5$; thus, at a given $z$ and $M_{*}$, there is an anti-correlation between SFR and metallicity in the simulations.

Our results regarding the increasing metallicity of the SF gas with time are consistent with Segers et al. (2016a), who show that recycling of stellar mass loss in EAGLE becomes increasingly important for fuelling star formation towards lower redshift. Also, Segers et al. (2016a) determined a characteristic mass, $M_{*}\approx 10^{10.5}{\rm M_{\sun}}$, below which, the contribution of recycled mass increases with mass and above which, it decreases with mass. The negligible evolution that we obtained for the $M_{*}=10^{10.5}-10^{11.0}{\rm M}_{\sun}$ bin in Fig. 2 and the flattening of the slope of the $z=0$ MZR at high masses (Fig. 1) might reflect the transition determined by Segers et al. (2016a). These authors claimed that this transition reflects the transition from stellar to AGN feedback.

We recall that otherwise we use 30 kpc spherical apertures. We have verified that the main trends of the scaling relations presented so far are preserved if the global properties are measured in different apertures (see Appendix A).

Finally, it is worth noting that De Rossi et al. (2015a) found a negligible evolution of the MZR in the case of the GIMIC simulations. The improved subgrid physics implemented in EAGLE leads to a better agreement with the data. A comparison of GIMIC and EAGLE MZRs is presented in Appendix B. There, we show that the implementation of energy feedback from star formation that depends on the local density and metallicity leads to an evolution of the MZR consistent with the observed trend (see, e.g. Fu et al. 2012 and Davé et al. 2017, for similar findings derived from other models and simulations). On the other hand, a constant energy feedback, as that implemented in GIMIC, leads to a negligible MZR evolution.

3.3 Impact of AGN feedback

In the previous section, we analysed the evolution of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation using the high-resolution simulation Recal-L025N0752. However, to investigate the role of AGN feedback on the MZR, we need to focus on the trends at high masses and, thus, we will need to use intermediate resolution simulations. In the overlap region at high masses, the trends followed by the simulated $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation seem to be more robust against resolution (Fig. 1) and different resolution runs agree. At $z=0$, the intermediate- and high-resolution simulations exhibit similarly flat slopes at $M_{*}\ga 10^{10}\ {\rm M}_{\sun}$, though with a decrease by $\approx 0.1$ dex in the normalization when using the recalibrated model compared to the reference one. The analysis of different simulation outputs ($z\approx$ 0.0, 0.1, 0.18, 0.27, 0.37, 0.50, 0.62 and 0.74) shows that the flattening of the predicted MZR at the high-mass end extends to $z\approx 0.7$ (see Fig. 2), in agreement with some observations (e.g. de los Reyes et al., 2015). In this section, we will show that the flattening of the simulated MZR is mainly regulated by AGN feedback.555We note that, for MZR studies, galaxies with evidence of current AGN activity are removed from observed samples. However, BH feedback may have occurred in cyclic episodes that affected also the selected observed sub-samples. In all EAGLE massive galaxies, the activity of SMBHs also varies with time.

Within the EAGLE suite of cosmological simulations, the impact of the AGN parameters was explored using the intermediate-resolution “L050N0752” version of the simulations. In addition to the Ref-L050N0752 run, three other simulations are analysed here: NOAGN-L050N0752, AGNdT8-L050N0752 and AGNdT9-L050N0752 (see Table 1). The simulation NOAGN-L050N0752 does not include AGN feedback. In the case of simulations AGNdT8-L050N0752 and AGNdT9-L050N0752, the temperature increase of the gas caused by AGN feedback has been set to $\Delta T_{\rm AGN}=10^{8}$K and $\Delta T_{\rm AGN}=10^{9}$K, respectively.666Note that the simulation AGNdT9-L050N0752 analysed by Schaye et al. (2015) assumes $C_{\rm visc}/2\pi=10^{2}$, while here $C_{\rm visc}/2\pi=10^{0}$ (Table 1, Crain et al. 2015), as we are interested in studying single-parameter variations. In Fig. 3, we compare the $z=0$ $M_{*}-M_{\rm BH}/M_{*}$ (where $M_{\rm BH}$ depicts the black hole mass), $M_{*}-{\rm SFR}$ and $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relations predicted by the four aforementioned simulations. All considered mass bins contain more than 10 galaxies inside them. We can see that BH growth and the corresponding AGN feedback sets in at $M_{*}\sim 10^{10}\ {\rm M}_{\sun}$, with $M_{\rm BH}/M_{*}$ increasing with $M_{*}$, as expected. As discussed in Crain et al. (2015), since $M_{\rm BH}$ is determined to first order by halo mass (Booth & Schaye 2010, see also Bower et al. 2017), the offset of the $M_{*}-M_{\rm BH}/M_{*}$ relation at $M_{*}\ga 10^{10}\ {\rm M}_{\sun}$ is related to the different halo mass associated with galaxies in this mass range. At a given $M_{*}$, SFR tends to decrease as $\Delta T_{\rm AGN}$ increases (middle panel of Fig. 3) because the less frequent but more energetic feedback episodes associated with a higher $\Delta T_{\rm AGN}$ are more efficient at regulating SFR in massive galaxies. In particular, if AGN feedback is suppressed, higher values of $M_{*}$ are obtained. Differences between the MZRs obtained from the NOAGN, AGNdT8, AGNdT9 and reference simulations (bottom panel of Fig. 3) are also appreciable at $M_{*}\ga 10^{10}{\rm M}_{\sun}$, reaching a metallicity offset of $\approx 0.3$ dex at $M_{*}\sim 10^{11}{\rm M}_{\sun}$. It is clear that, when AGN feedback is turned off, the MZR does not exhibit a flattening at high masses. On the contrary, the MZR slope increases towards higher masses in the case of the NOAGN model. As $\Delta T_{\rm AGN}$ increases, the MZR slope decreases and, in particular, above $\Delta T_{\rm AGN}\approx 5\times 10^{8}$K, AGN feedback can generate an inversion of the relation between stellar mass and metallicity turning it from a correlation into an anti-correlation. These trends reflect the behaviour obtained for the $M_{*}-{\rm SFR}$ relation: higher $\Delta T_{\rm AGN}$ leads to less on-going star formation. Additionally, we see that a higher $\Delta T_{\rm AGN}$ leads to larger scatter in $M_{\rm BH}/M_{*}$, SFR and ${\rm O/H}|_{\rm SF,gas}$ at a given $M_{*}$ at $10^{10}\la M_{*}/{\rm M}_{\sun}\la 10^{11}$. Taking into account the different AGN feedback histories that can take place in real galaxies, we can conclude that part of the dispersion of the observed MZR at the high-mass end may be associated with AGN feedback, in which case its study could provide important constraints on galaxy formation models (see also Section 7).

In Fig. 4, we plot ${\rm OH}|_{\rm SF,gas}$ as a function of the halo mass ($M_{200}$) for central simulated galaxies.777The halo mass, $M_{200}$ is defined as the total mass contained within the virial radius $R_{200}$, which is the radius within which the mean internal density is 200 times the critical density, $3H^{2}/8\pi G$, centred on the dark matter particle of the corresponding FoF halo with the minimum gravitational potential. More than 15 galaxies are included in each mass bin. We obtained a similar trend to that derived previously for $M_{*}$. At a given halo mass, higher $\Delta T_{\rm AGN}$ yields lower metallicities towards higher masses. When increasing $\Delta T_{\rm AGN}$ from $10^{8}$K to $10^{9}$K, the metallicity associated with central galaxies in massive halos decreases by around $0.2$ dex. If AGN feedback is suppressed, the median metallicity of massive galaxies increases significantly, reaching an offset $\ga 0.3$ dex at $M_{200}\approx 10^{13}\ {\rm M}_{\sun}$ relative to the case of $\Delta T_{\rm AGN}=10^{8}{\rm K}$. In section 7.1.2, we try to disentangle the diverse effects caused by AGN feedback in galaxies that lead to the decrease of the global metallicity of these systems. A further investigation is left for a future work.

4 The $M_{*}-Z_{\rm*}$ relation

The predicted relation between $M_{*}$ and stellar metallicity ($Z_{*}$) is shown in Fig. 5 (curve with error bars) for Recal-L025N0752. The number of galaxies per mass bin is $N_{\rm bin}\ga 10$ at all considered $z$. Observational data at $z\sim 0$ from Gallazzi et al. (2005) are represented as black circles with error bars, which indicate the median relation together with the 16th and 84th percentiles. The fitted relation to SDSS data given in Panter et al. (2008) is shown with a dotted green line. We also show findings by Gallazzi et al. (2014) at $z=0.7$ as green asterisks with error bars, depicting the median relation together with the 16th and 84th percentiles. Results at $z\approx 3$ reported by Sommariva et al. (2012) are represented with red diamonds and error bars, which denote the measured metallicities and corresponding uncertainties. For some galaxies, Sommariva et al. (2012) were not able to determine metallicity uncertainties so that no error bars are shown in those cases (the reader is referred to that paper for further details).

We see that, at a given $M_{*}$, $Z_{*}$ tends to decrease with $z$, as previously reported by Guo et al. (2016) for this simulation888For this analysis, Guo et al. (2016) considered all particles inside EAGLE galaxies while, in this work, we take into account particles within an aperture of 30 pkpc. (we referred the reader to that work for a comparison between EAGLE results and semi-analytical models). We can see that, at $z=0$, the median simulated and observed relations tend to agree at high masses. At the low-mass end, simulation predicts higher metallicities, on average, than observations. It is clear that the scatter in the data at $z\approx 0$ is larger than predicted by simulations (see also Schaye et al. 2015). This was also noted by De Rossi et al. (2015b) when using the GIMIC simulations and might be related to uncertainties in observed stellar metallicity determinations. Gallazzi et al. (2005) adopted a Bayesian statistical approach to obtain $Z_{*}$ for individual galaxies, estimating the 68 per cent confidence interval within which this parameter is constrained. Observational errors on $Z_{*}$ have an average value of 0.12 dex but can extend to $\approx 0.30$ dex (see Gallazzi et al., 2005, for further details).

At $z\approx 0.7$, observed $Z_{*}$ is only reported at $M_{*}\ga 10^{10.5}\ {\rm M}_{\sun}$ and exhibits lower median values than simulated ones by $\approx 0.2$ dex. Gallazzi et al. (2014) estimated $Z_{*}$ for individual galaxies similarly to Gallazzi et al. (2005), obtaining an average uncertainty on stellar metallicity of 0.3 dex, slightly worse than for the local observed sample. Thus, simulations and observations results agree within the error. In the mass range where simulations and observations overlap at $z\approx 3$, the simulated $M_{*}-Z_{*}$ relation is consistent with the stellar metallicities reported by Sommariva et al. (2012).

It is worth noting that the simulated $Z_{*}$ is calculated weighting by mass while, in the case of Gallazzi et al. (2005, 2014) and Sommariva et al. (2012), the observed $Z_{*}$ is a luminosity-weighted quantity. On the other hand, the approach used by Panter et al. (2008) yields a mass fraction weighted metallicity. From Fig. 5, we see that the slope of the simulation $M_{*}-Z_{*}$ relation shows a better agreement with the latter authors. To analyse this issue in more detail, we re-calculated the simulated $Z_{*}$ weighting by light in different photometric bands (u, g, r, i, z, Y, J, H, K). 999For these calculations, we used photometric data of stellar particles available in the public EAGLE database (McAlpine et al., 2016). We obtained that light-weighted metallicities tend to be higher than mass-weighted metallicities by $\la 0.2$ dex, with the larger differences towards lower masses. In addition, the scatter of the simulated $M_{*}-Z_{*}$ relation slightly increases, by $\la 0.1$ dex, when using light-weighted metallicities instead of mass-weighted metallicities. Thus, the use of light- instead of mass-weighted metallicities does not generate significant changes in the trends predicted for the simulated $M_{*}-Z_{*}$ relation.

The $M_{*}-{\rm[Fe/H]}_{*}$ relation has been widely studied in the literature for early-type galaxies (e.g. Conroy, Graves & van Dokkum, 2014). Observed ${\rm[Fe/H]}_{*}$ increases with $M_{*}$ for more massive galaxies, as expected. We have verified that the simulation reproduces this trend, however, a more detailed analysis is beyond the scope of this paper. For a study of the $M_{*}-{\rm[\alpha/Fe]}_{*}$ and ${\rm age}-{\rm[\alpha/Fe]}_{*}$ relations in EAGLE simulations, the reader is referred to Segers et al. (2016b).

A comparison between the metallicity of the different baryonic components of EAGLE galaxies is presented in Fig. 6. We show how the ratio between gas ($Z_{\rm gas}$) and stellar ($Z_{*}$) metallicity varies as a function of $M_{*}$ considering both gas components: the star forming gas and the non star forming gas (NSF gas) phases. The SF gas component of a given galaxy was calculated by summing the mass of all gas particles with ${\dot{\rho}}_{*}>0$ (see Section 2.1), while the NSF gas component has ${\dot{\rho}}_{*}=0$. Each mass bin associated to the median relations shown in Fig. 6 contains more than 10 galaxies at all considered $z$. The quantity $Z_{\rm gas}/Z_{*}$ can be considered to be a measure of the current compared to past average chemical enrichment of the ISM of galaxies. The SF and NSF gas-components present clearly distinct levels of enrichment with respect to stars at all $z$ considered.

The SF gas component is more metal-enriched than the stellar component at all $z$ considered and the opposite is true for the NSF gas. For the SF gas, the ratio $Z_{\rm SF,gas}/Z_{*}$ tends to decrease from $\approx 1.5-2.0$ at $M_{*}\sim 10^{9}{\rm M}_{\sun}$ to $\approx 1.0$ at $M_{*}\ga 10^{10.5}{\rm M}_{\sun}$, with no evidence for a dependence on $z$. Fig. 5 also shows that, at a given $M_{*}$, $Z_{*}$ tends to be lower than the metallicity of the SF gas, especially in low-mass systems (compare with Figs. 1 and 2). In addition, the $M_{*}-Z_{*}$ relation is steeper than the $M_{*}-Z|_{\rm SF,gas}$ relation, with a slope increasing towards high $z$. $Z|_{\rm SF,gas}$ can be directly associated with the recent metal enrichment of the galaxy while $Z_{*}$ reflects the average chemical evolution of the system during its periods of star formation.101010Note, however, that, in the case of observations, $Z_{*}$ uses to be luminosity-weighted and tends to be biased to the most luminous stars. Hence, star formation in smaller galaxies tends to be driven by more metal-enriched gas than that associated with their past star formation, represented by $Z_{*}$. As smaller galaxies have larger gas fractions (see Section 6), the chemical evolution at the low-mass end develops faster as can be seen from Figs. 2 and 5. On the other hand, more massive galaxies form stars from gas with a similar level of enrichment as their average stellar component. This is consistent with the flattening of the $M_{*}-Z_{\rm SF,gas}$ and $M_{*}-Z_{\rm*}$ relations at high stellar masses (see Zahid et al. 2014a for similar findings derived from analytical methods). The metallicity of high-mass galaxies saturates at early times at $12+\log_{10}({\rm O/H})\approx 8.9-9.0$, exhibiting less significant differences between SF gas and stars.

In the case of the NSF gas component, $Z_{\rm NSF,gas}/Z_{*}$ decreases towards higher stellar masses and increases with decreasing $z$. Thus, at lower $z$, the metallicities of the different baryonic phases of more massive galaxies are more homogeneous, as a result of accumulated mixing of metals over time. These results also suggest that star formation in massive galaxies might have been regulated by supply of metal-poor gas in the NSF phase that was transferred to the more enriched SF one. As most massive galaxies exhibit low-gas fractions and nearly homogeneous metallicity distributions between their different baryonic components, their subsequent chemical evolution may be driven by mergers (see Yates & Kauffmann 2014 for a similar analysis in the context of semi-analytical models).

5 Effective yields

In a closed-box model (e.g. Schmidt, 1963; Tinsley, 1980), where there are neither inflows nor outflows, the metallicity is a simple function of the stellar yield ($y$) and gas fraction ($f_{\rm gas}=M_{\rm gas}/(M_{\rm gas}+M_{*})$):

[TABLE]

where $y$ represents the ratio of the mass in metals synthesised and released to the net mass that remains locked-up in stars. If we were to unrealistically assume a constant stellar yield (e.g. Garnett, 2002; Tremonti et al., 2004), then the level of metal enrichment only depends on the gas fraction. By inverting Eq. (4), we can obtain the effective yields, $y_{\rm eff}=Z_{\rm SF,gas}/\ln(f_{\rm SF,gas}^{-1}$). Note that we calculate $y_{\rm eff}$ from the SF gas phase as this is the component usually observed. Clearly, if a galaxy evolves as a closed box, then $y_{\rm eff}=y$ while, in both cases of metal-enriched outflows and infall of primordial gas, $y_{\rm eff}\leq y$ (Edmunds, 1990). Thus, $y_{\rm eff}$ can be used to diagnose the role of inflows and/or winds on the evolution of galaxies.

In Fig. 7, we plot $y_{\rm eff}$ as a function of the mass of baryons ($M_{\rm bar}=M_{*}+M_{\rm SF,gas}$) for simulations and observations at $z\approx 0$, 2 and 3. Black, blue and red lines with error bars denote simulated median relations at $z=0$, 2 and 3, respectively. Error bars depict the 25th and 75th percentiles. Observational results by Tremonti et al. (2004) at $z\approx 0$ (black diamonds), Erb et al. (2006) at $z\approx 2$ (blue asterisks) and Troncoso et al. (2014) at $z\approx 3$ (red circles) are shown for comparison. In order to estimate $y_{\rm eff}$, Tremonti et al. (2004), Erb et al. (2006) and Troncoso et al. (2014) derived gas masses from the observed star formation rates assuming a Schmidt law (Kennicutt, 1998). The reader is referred to those papers for more details. Note that the determination of observed yields is also affected by the metallicity calibration (Eq. 4) and, as we have seen, there are significant discrepancies between observed metallicities obtained from different calibrators (Fig. 1). Erb et al. (2006) used the PP04 N2 method to determine metallicities, which yields lower metallicities that those derived by the T04 method (see Fig. 1). Troncoso et al. (2014) adopted the same method as in Maiolino et al. (2008), which gives higher metallicities than the PP04 N2 method, with the differences ranging between $\approx 0.3$ and 0.1 dex depending on mass. Thus, as explained in Section 3.2, an offset was applied to observed yields so that the corresponding $12+{\log}_{10}({\rm O/H})=9.0$ at $M_{*}=10^{10.5}\ {\rm M_{\sun}}$ at $z\approx 0$. To renormalize metallicities in the case of Erb et al. (2006) and Troncoso et al. (2014), we used the local MZR used by each author to contrast their higher-$z$ results. Each author used a local MZR derived from the same metallicity calibrator employed for obtaining the corresponding high-$z$ MZR.111111Offsets applied to observed effective yields: -0.07 dex (Tremonti et al., 2004), +0.35 dex (Erb et al., 2006) and -0.02 dex (Troncoso et al., 2014).

We see that there is negligible evolution of the simulated $y_{\rm eff}-M_{\rm bar}$ relation at $M_{*}\la 10^{10.2}\ {\rm M}_{\sun}$ below $z=3$, while more massive systems tend to present higher $y_{\rm eff}$ at higher $z$. Considering the differences between the metallicity calibrations applied by Tremonti et al. (2004) and Erb et al. (2006), observational results also suggest a decrease of $y_{\rm eff}$ for massive galaxies towards $z\approx 0$. Regarding the scatter of the $y_{\rm eff}-M_{\rm bar}$ relation, in the case of the observations, the scatter is large partly because of the significant uncertainties in the determination of gas masses. On the other hand, we have verified that in EAGLE the larger scatter of the simulated $y_{\rm eff}$ at higher masses is related to AGN feedback. By comparing the AGNdT8-L050N0752 and AGNdT9-L050N0752 models, we found that an increase of $\Delta T_{\rm AGN}$ leads to a larger scatter at $M_{\rm bar}\ga 10^{10}{\rm M}_{\sun}$. The simulation NOAGN-L050N0752 predicts a negligible scatter for the $y_{\rm eff}-M_{\rm bar}$ relation.121212 We note that, in the case of Recal-L025N0752 at $z=0$, the number of galaxies in the highest mass bin, where the scatter is larger, is $N_{\rm bin}\approx 15$. For ”L050N0752” simulations, at the same bin ($M_{\rm bar}\approx 10^{10.8}\ {\rm M_{\sun}}$) and redshift, $N_{\rm bin}\approx 100$ (NOAGN), $N_{\rm bin}\approx 100$ (AGNdT8) and $N_{\rm bin}\approx 50$ (AGNdT9).

Tremonti et al. (2004) fitted the observed data using the function:

[TABLE]

where the asymptotic value $y_{0}$ represents the true yield (if no metals are lost). These authors obtained $y_{0}=0.0104$ and $M_{0}=3.67\times 10^{9}\ {\rm M}_{\sun}$ (dashed curve in Fig. 7) for the observed sample. The simulated $z=0$ $y_{\rm eff}-M_{\rm bar}$ relation follows roughly the observed trends at $M_{*}\la 10^{10}\ {\rm M}_{\sun}$ but it is inverted at higher masses. Thus, we fitted the simulated relation at $z=0$ using the same function as Tremonti et al. (2004) but including only systems with $M_{\rm bar}\la 10^{10}\ {\rm M}_{\sun}$. We obtained $y_{0}\approx 0.013$ and $M_{0}\approx 3.4\times 10^{8}{\rm M}_{\sun}$. As the fitted curve for simulations follows tightly the median relation, for the sake of clarity, we do not include it in Fig. 7.

In the simulations, metal enrichment is implemented following the assumed stellar yields tables (Portinari, Chiosi & Bressan, 1998; Marigo, 2001; Thielemann et al., 2003) but, subsequently, metal abundances are modified by different processes such as inflows and outflows of gas. The trends predicted for $y_{\rm eff}$ provide clues about the role of gas infall and/or winds on the evolution of galaxies (e.g. Dalcanton, 2007). The decrease of the simulated $y_{\rm eff}$ towards lower masses suggests that these systems have been affected by efficient outflows of metal-enriched SN ejecta. Because small galaxies tend to have high mean gas fractions (see Sections 6), the infall of metal-poor gas could not drive significant changes in their effective yields (Dalcanton, 2007).

At high masses, $y_{\rm eff}$ tends to decrease towards $z=0$. In particular, some systems depart from the mean relation towards very low $y_{\rm eff}$ due to their very low gas fractions. A possibility is that AGN feedback drives metal-enriched outflows in more massive galaxies. In fact, we have verified that an increase of $\Delta T_{\rm AGN}$ leads to lower average $y_{\rm eff}$ at $M_{\rm bar}\ga 10^{10}{\rm M}_{\sun}$. However, infall of unenriched material probably generates the largest reductions in $y_{\rm eff}$ in gas-poor massive galaxies (Dalcanton, 2007) (see 7.1.2, for a more detailed analysis of the effects of AGN feedback on massive galaxies). Merger events can also play a significant role in determining the effective yields in massive galaxies.

6 Fundamental metallicity relation and scatter of the MZR

The scatter in the observed MZR is partially due to observational uncertainties. However, as discussed in Section 1, there is evidence for secondary dependences of the metallicity on properties such as gas fractions or SFR. In this Section, we explore secondary metallicity dependences of EAGLE galaxies that can explain part of the scatter of the MZR at fixed mass.

In Fig. 8 (top panels), we analyse the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation obtained for Recal-L025N0752 at $z=0$ by binning the galaxy sample according to star-forming gas fraction (left panels), specific SFR (sSFR$={\dot{M}_{*}}/M_{*}$, middle panels) and stellar mass-weighted mean age (right panels). The gas fractions shown in left panels were calculated using the star-forming gas component ($f_{\rm SF,gas}=M_{\rm SF,gas}/(M_{*}+M_{\rm SF,gas})$) in the galaxies but similar trends are obtained if considering the total gas ($f_{\rm gas}=M_{\rm gas}/(M_{*}+M_{\rm gas})$) instead. In each of the upper panels, we show results for three bins corresponding to the aforementioned properties, as indicated in the figure. In the bottom panels, 2D colour histograms are shown, where filled squares are coloured according to the median $f_{\rm SF,gas}$ (left panel), sSFR (middle panel) and stellar mass-weighted mean age (right panel). As a reference, the median $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation obtained from the whole galaxy sample is shown in the three bottom panels as a blue dashed line.

We can see that, at a given stellar mass, there is a clear dependence of metallicity on these three parameters. At fixed $M_{*}$, systems with higher $f_{\rm SF,gas}$ or higher sSFRs have lower values of ${\rm O/H}|_{\rm SF,gas}$, consistently with the trends found in certain observational studies (e.g. Kewley & Ellison, 2008; Salim et al., 2014; Bothwell et al., 2013). We also obtained lower values of ${\rm O/H}|_{\rm SF,gas}$ for younger galaxies, and all these trends are stronger at lower masses. Our findings are consistent with a scenario in which higher fractions of metal-poor gas drive higher sSFRs in smaller galaxies and, consequently, these subsamples of galaxies tend to be younger than more metal-enriched galaxies of similar stellar masses.

In the case of higher $M_{*}$ galaxies, as we have previously mentioned, the metallicity saturates towards an asymptotic value at early times. However, it is worth noting that there is a hint of an inversion of the dependence of metallicity on the $f_{\rm SF,gas}$ and sSFR at the massive end of the MZR. Massive galaxies with lower $f_{\rm SF,gas}$ and sSFR tend to have lower metallicities, contrary to the trend shown by smaller galaxies. This issue might be a consequence of the influence of AGN feedback which can produce a decrease in the SF activity in these galaxies and, consequently, the quenching of their chemical evolution.

In the top panels of Fig. 8, some bins used to construct the median relations are poorly sampled ($5\leq N_{\rm bin}<10$, circles). Also, in the bottom panels, the 2D histograms include bins with less that 10 galaxies in some cases. In order to check if the obtained trends are affected by sampling effects due to the limited volume of the simulation ($25^{3}\ {\rm cMpc^{3}}$), we repeated the previous analysis for the simulations Ref-L050N0752 and Ref-L100N1504, for which the number of galaxies per bin in the 2D histograms increases by up to $\approx 100$ and $\approx 700$, respectively. Although the shape and normalization of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation vary when using the simulations Ref-L050N0752 and Ref-L100N1504 (Fig. 1), the secondary dependences of metallicity (at a fixed $M_{*}$) on $f_{\rm SF,gas}$, sSFR and age are also present in these larger-volume simulations. Thus, these dependences seem not to be an artefact caused by sampling issues.

To corroborate if AGN feedback might generate a decrease in the SF activity and quench the metallicity evolution of galaxies, we compared NOAGN-L050N0752, AGNdT8-L050N0752 and AGNdT9-L050N0752 models. We obtained that for higher $\Delta T_{\rm AGN}$, the inversion of the dependence of metallicity with the $f_{\rm SF,gas}$ and sSFR is more significant (see Section 7). When AGN feedback is suppressed, we did not find a clear trend for such inversion. Moreover, we found that massive galaxies ($M_{*}\ga 10^{10.5}{\rm M}_{\sun}$) tend to exhibit lower $f_{\rm SF,gas}$ and metallicities as the ratio between the black hole mass and $M_{*}$ ($M_{\rm BH}/M_{*}$) increases.

The existence of a relation between ${\rm O/H}|_{\rm SF,gas}$, $M_{*}$ and SFR (FMR) has been investigated by many authors in recent years. However, as discussed in Section 1, the features of this relation and its strength remain poorly constrained because of the large uncertainties and systematic errors that affect observational studies. Fig. 9 clearly indicates that EAGLE galaxies follow a well-defined relation between $M_{*}$, SFR and ${\rm O/H}|_{\rm SF,gas}$ at least below $z\approx 5$. Only mass bins with $N_{\rm bin}>10$ were considered for this analysis: $N_{\rm bin}=13-296,\ 25-142$ and $23-82$ at $\log_{10}({\rm SFR}/{\rm M}_{\sun}yr^{-1})\approx-0.25,\ 0.25$ and 0.75, respectively. The predicted 3D relation agrees remarkably well with those given in Mannucci et al. (2010) for observed galaxies in the local Universe, especially at $M_{*}<10^{10.5}{\rm M}_{\sun}$ (see also Fig. 11, top left panel, for an analysis of residuals). At a given SFR, metallicity shows a strong correlation with $M_{*}$ with a flattening of the slope towards higher masses. At a given $M_{*}$, the metallicity of low-mass systems decreases with SFR while, for massive galaxies, the metallicity increases with SFR. Similar trends were reported by Yates, Kauffmann & Guo (2012) and De Rossi et al. (2015b) in the context of semi-analytical models and hydrodynamical simulations, respectively. Lagos et al. (2016) explored the relation between SF gas metallicity, SFR and $M_{*}$ in the EAGLE Ref-L100N1504 simulation. They found good agreement with the data at high masses ($M_{*}\ga 10^{10}\ {\rm M}_{\sun}$) while low-mass galaxies depart from the measured relation by up to $\approx 0.4$ dex towards higher metallicities. For the higher-resolution Recal-L025N0752 simulations used here, the agreement with the data is better because of the use of the recalibrated model (see Section 2.2), which yields more efficient stellar feedback (see Section 7).

From Fig. 9 we can see that Recal-L025N752 reproduces the observed flattening of the slope of the FMR towards high $M_{*}$. Nevertheless, in the case of massive galaxies with low SFRs, the impact of AGN feedback seems to have been very strong, generating a departure towards lower metallicity from observations by Mannucci et al. (2009). We note, however, that other observational studies (e.g. Yates, Kauffmann & Guo, 2012; Lara-López et al., 2013; Andrews & Martini, 2013; Zahid et al., 2013b) found lower metallicities in low-SFR massive systems, which is more consistent with EAGLE trends.

Mannucci et al. (2010) claimed that the FMR does not evolve significantly below $z\approx 2.5$. On the other hand, Troncoso et al. (2014) reported that observed galaxies at $z\approx 3.4$ lie below the FMR by a factor of $\sim 10$, suggesting that the observed FMR is established at $z<3$. More recent works suggested that the trend of the SFR dependence is preserved at $z<2.3$ but its strength evolves (Salim et al., 2015; Brown, Martini & Andrews, 2016; Grasshorn Gebhardt et al., 2016). EAGLE predicts a FMR that is already in place at $z\approx 5$ and does not evolve significantly with time.

The FMR might be caused by the existence of a more fundamental relation between $M_{*}$, metallicity and gas fraction (e.g. Bothwell et al., 2013; Hughes et al., 2013). De Rossi et al. (2015b) obtained that dependence in the case of GIMIC galaxies while Lagos et al. (2016) found similar results by performing a principal component analysis over different galaxy properties in the EAGLE Ref-L100N1504 simulation. At a given $M_{*}$, galaxies with higher gas fractions tend to exhibit higher SFRs and, as we have discussed, these galaxies have lower gas-phase metallicities. Fig. 10 shows the ${\rm O/H}|_{\rm SF,gas}-f_{\rm SF,gas}$ relation in Recal-L025N0752 simulations at $z=0,1,2$ and 3. The $M_{*}$ range associated to each $z$ can be seen in previous figures (e.g. Fig. 2). Again, we have calculated gas fraction using the star-forming gas component in the galaxies. We can see that systems with larger gas fractions exhibit lower metallicities. As $f_{\rm SF,gas}$ increases from 0 to 0.8, the metallicity decreases by $\approx 1$ dex. The scatter around the mean relation is relatively low ($\la 0.3$ dex), with the only exception of very gas-poor systems at $z=0$. We have verified that the scatter of the ${\rm O/H}|_{\rm SF,gas}-f_{\rm SF,gas}$ relation is generated by the different stellar masses of the systems that are included in each $f_{\rm SF,gas}$ bin. We found that, at $z\ga 0.6$ and for a fixed $f_{\rm SF,gas}$, metallicity correlates with $M_{*}$, as expected. On the other hand, at $z\la 0.6$, metallicity correlates with $M_{*}$ in the case of systems with $f_{\rm SF,gas}\ga 0.1$, while the metallicities of gas-poor systems ($f_{\rm SF,gas}\la 0.1$) tend to decrease with $M_{*}$. This behaviour is consistent with the trends obtained in Fig. 8, left panels.

Fig. 10 shows that $f_{\rm SF,gas}$ decreases and ${\rm O/H}|_{\rm SF,gas}$ increases systematically for decreasing $z$, just as in a closed-box model and most open-box models. Thus, for the sake of comparison, we estimated the correlation between ${\rm O/H}|_{\rm SF,gas}$ and $f_{\rm SF,gas}$ obtained from a closed-box model and a simple model with outflows. Following Erb et al. (2006), we implemented a simple model derived from the closed-box model to include gas outflow at a rate $\dot{M}$ (${\rm M}_{\sun}{\rm yr^{-1}}$), which is a fraction $f$ of the SFR. The metallicity of the ejected gas is assumed to be the same as the metallicity of the gas that remains in the galaxy. Under these considerations, the metallicity is

[TABLE]

Note that the closed-box model (Eq. 4) is obtained from this general expression taking $f=0$, while observations suggest that outflows eject gas from galaxies at a rate that is comparable to the SFR and may be higher ($f\ga 1$, Veilleux, Cecil & Bland-Hawthorn 2005). The curves without error bars in Fig. 10 show the results corresponding to a closed-box model with $y=Z_{\sun}$ (dotted black line) and a simple model with $y=1.5\times Z_{\sun}$ and an outflow rate $\dot{M}=1.5\times{\rm SFR}$ (dashed red line). Interestingly, the ${\rm O/H}|_{\rm SF,gas}$-$f_{\rm SF,gas}$ relation obtained from the simulations follows roughly the trend associated with the closed-box model. This can be explained considering the flat shape and negligible evolution obtained for the $y_{\rm eff}-M_{\rm bar}$ relation (Fig. 7). However, the simple outflow model reproduces better the simulated trends associated with gas-rich systems at all $z$ and also improves the description of the behaviour of gas-poor systems towards $z=0$. As discussed before, simulated systems drive metal-enriched outflows and accrete metal-poor gas as well as undergoing merger events, whose relative efficiencies might vary with time. These processes cannot be captured by means of simple analytical models but it is interesting that such models are able to roughly described the simulation results.

It is worth noting that we have analysed the metallicity evolution for all (i.e. central and satellite) galaxies. In EAGLE, satellites tend to have higher metallicities than central galaxies of similar stellar masses (Bahé et al., 2017) and this is also observed (Pasquali et al., 2010; Pasquali, Gallazzi & van den Bosch, 2012). At a given stellar mass and at $z=0$, for example, the median metallicity of satellites is higher than for central galaxies by $\approx 0.2$ dex at $M_{*}\sim 10^{9}{\rm M}_{\sun}$ and by $\approx 0.1$ dex at $M_{*}\sim 10^{10.5}{\rm M}_{\sun}$. Also, the MZR of satellites presents a scatter that is a factor $\ga 2$ larger than for central galaxies. Therefore, the inclusion of satellites in the main sample increases the scatter in the MZR. As discussed by De Rossi et al. (2015b), ram-pressure stripping could lead to higher global metallicities of satellites by removing the reservoir of metal-poor gas usually located in the outer regions of galaxies. The origin of the satellite metallicity enhancement in Ref-L100N1504 was studied by Bahé et al. (2017). These authors found an excess gas and stellar metallicities in satellites in agreement with observations, except for stellar metallicities at $M_{*}\la 10^{10.2}{\rm M}_{\sun}$ where the predicted excess is smaller than observed. Stripping of low-metallicity gas and suppression of metal-poor inflows seem to play an important role on driving the enhancement of gas metallicity.

Finally, as the predicted relation between star-forming gas fraction and metallicity shows no significant evolution and only involves two key global galaxy properties, we conclude that it is more fundamental than the FMR, at least in these simulations. It is worth noting that similar conclusions were obtained by Zahid et al. (2014a) by using empirical-constrained analytical models of chemical evolution. According to these authors, the redshift independence of the relation between metallicity and stellar-to-gas ratio suggests that it is a universal relation. The MZR would originate from this universal relation due to the evolution of the gas content of galaxies. Our results seem to support these previous findings.

7 Discussion

7.1 Subgrid parameters

We have seen that the EAGLE Recal-L025N0752 simulation reproduces the observed trends for the MZR and FMR remarkably well. In this Section, we explore the effects of varying some model parameters to assess their influence on the observed features of the fundamental metallicity relations. The different simulations used for this study and the corresponding parameter variations can be found in Table 1. In addition to checking resolution effects on the FMR (”N0376” vs ”N0752” runs) and comparing the reference and recalibrated models, we study the impact of AGN feedback (NOAGN vs AGNdT9 models) and the stellar feedback efficiency (WeakFB vs StrongFB models). We also analyse the effects of changing the slope of the equation of state imposed on the ISM (eos1 vs eos5/3 models) and the star formation law (KSLow vs KSHigh models).

7.1.1 Residuals to the observed FMR

In Fig. 11, we compare the residuals of the simulated $M_{*}-{\rm O/H|}_{\rm SF,gas}$ relation, in different SFR bins, with respect to the observed FMR given by Mannucci et al. (2010, see shaded surface in Fig. 9). In the three upper panels, we analyse the effects of changing the resolution (low resolution ”L025N0356” vs high-resolution ”L025N0752” runs) and model (reference vs recalibrated subgrid parameters, see Table 1) for a fixed simulation volume131313 We have verified that similar trends are obtained if the simulation Ref-L100N1504 is used instead of Ref-L025N0356. $(25\ {\rm cMpc})^{3}$. In Sec. 3, we have already compare the MZRs predicted by these simulations. As noted before, the best overall agreement with the observed FMR is obtained for Recal-L025N0752. The residuals of this simulation are below $\approx 0.1$ dex for smaller galaxies ($M_{*}\la 10^{10.2}{\rm M}_{\sun}$) and below $\approx 0.2$ dex for more massive ones.

Comparing the results obtained from the Recal- and Ref-L025N0752 runs, it is clear that the larger differences are obtained at low masses and for systems with SFR $\la 1\ {\rm M}_{\sun}\ {\rm yr}^{-1}$, where the residuals derived from the reference model are $\approx 0.1$ dex higher than those resulting from the recalibrated model. Thus, the averaged low-mass $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation predicted by the Ref-L025N0752 simulation is flatter than that associated with Recal-L025N0752, as already noted by Schaye et al. (2015). The stronger outflows generate a decrease in the metallicity of the ISM in the recalibrated model, producing a better match to the observed FMR.

To analyse the effects of resolution without changing subgrid parameters (the so-called strong convergence test, in the terminology introduced by Schaye et al. 2015) and volume, we can compare the simulations Ref-L025N0752 and Ref-L025N0376, shown in the middle and right upper panels of Fig. 11, respectively. The results seem to be reasonably robust against numerical resolution as differences in the metallicity residuals are below $\approx 0.1$ dex considering the whole mass range.

In Fig. 11, second row, we explore how the parameter that regulates the efficiency of AGN feedback events affect the FMR. We make use of simulations NOAGN-L050N0752 (AGN feedback suppressed) and AGNdT9-L050N0752 ($\Delta T_{\rm AGN}=10^{9}$K), which yield higher AGN feedback heating temperature than the reference model ($\Delta T_{\rm AGN}=10^{8.5}$K). We can see that variations in AGN feedback have negligible effects at low stellar masses ($M_{*}\la 10^{10.2}{\rm M}_{\sun}$). At the massive end ($M_{*}\ga 10^{10}{\rm M}_{\sun}$) of the FMR, AGN feedback produces a significant decrease of the metallicity at a given $M_{*}$ with respect to the model without AGN, leading to a better agreement with the data. A higher heating temperature drives lower oxygen abundances, especially for systems with lower SFR. This last issue might be related to the fact that lower SFRs can also be related to a more efficient heating (see Section 3.3). The increase of $\Delta T_{\rm AGN}$ reduces the overall metallicities of galaxies leading to a better agreement with the data at ${\rm SFR}\ga 1{\rm M}_{\sun}\ {\rm yr}^{-1}$. In the case of systems with low SFRs (${\rm SFR}\la 1{\rm M}_{\sun}\ {\rm yr}^{-1}$), the decrease in metallicities obtained for higher $\Delta T_{\rm AGN}$ generates a more significant departure from the observed FMR. Thus, increasing the AGN feedback heating temperature tends to strengthen the anti-correlation between metallicity and SFR at high stellar masses.

The panels in the last two rows of Fig. 11 show the effects of halving and doubling the stellar feedback (first column), of changing the power-law dependence index of the star formation law ($n$, second column) and the power law slope of the polytropic equation of state imposed on the ISM ($\gamma_{\rm eos}$, third column). See Table 1 for details. The efficiency of stellar feedback has a significant influence on the slope of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation (see also Crain et al. 2015). Weak stellar feedback yields very high metallicities towards lower stellar masses generating an important departure from the observed FMR ($\approx 0.4$ dex). Interestingly, that departure does not depend on SFR at fixed stellar mass. However, weaker stellar feedback does not significantly affect the metallicities of higher-mass galaxies (compare WeakFB-L025N0376 simulations with the Ref-L025N0376 set, in the upper right panel). As we have seen, the abundances in massive galaxies are mostly regulated by AGN feedback. On the other hand, strong stellar feedback alters the abundances of systems spanning our whole stellar mass range, which are reduced by up to $\approx 0.3$ dex with respect to those predicted by the reference model. This decrease leads to a better agreement with the data in the case of systems with lower SFR and low stellar masses. However, the MZR obtained for the intermediate and high-SFR bins moves below the observed relation. We can also deduce that an increase of the stellar feedback efficiency tends to strengthen the anti-correlation between SFR and metallicity at low stellar masses. In the case of the population of massive galaxies in the lowest SFR-bin, variations in the stellar feedback parameters do not generate important effects on the metallicities. Because of their deeper potential-wells and lower SFRs, they probably do not experience significant outflows of metal-enriched material driven by stellar feedback.

Simulations KSLow-L025N0376 and KSHigh-L025N0376 were run using power-law indices in the star formation law (Eq. 2, section 2.1) of $n=1.0$ and $n=1.7$, respectively. In the reference model, $n=1.4$. A higher value of $n$ corresponds to a higher star formation efficiency, which yields lower SF gas fractions if star formation is self-regulated (Schaye et al., 2010; Haas et al., 2013). Because of the fundamental relation between SF gas fraction and metallicity, we thus expect metallicity to increase with $n$ at fixed mass. We can compare the effects of changing $n$ in the two lower panels in the middle column of Fig. 11. Increasing $n$ generates similar effects, at low stellar masses, as those obtained by reducing the stellar feedback efficiency. For higher $n$, the metallicities of smaller systems increase with respect to the data and the slope of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation tends to flatten. In the case of a lower $n$, the main effects can be seen at the lower SFR-bins, which describe $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relations displaced towards lower abundances that those predicted by the reference model. Decreasing $n$ below the reference value also weakens the anti-correlation between metallicity and SFR at low-masses but reinforces the correlation between both variables at the massive-end.

In the two lower right panels of Fig. 11, we compare simulations eos1-L025N0376 and eos5/3-L025N0376, for which the power law slope of the polytropic equation of state imposed on the ISM ($\gamma_{\rm eos}$) has been set at $\gamma_{\rm eos}=1.0$ and $\gamma_{\rm eos}=5/3$, respectively. For the reference model, $\gamma_{\rm eos}=4/3$. It is clear that the effects of varying $\gamma_{\rm eos}$ do not generate significant variations in the metallicity relations. This is consistent with the negligible effect on the gas fraction found by Haas et al. (2013).

Finally, from Fig. 11, it is clear that the reversal of the SFR-${\rm O/H}|_{\rm SF,gas}$ correlation at $M_{*}\sim 10^{10.3}\ {\rm M}_{\sun}$ is present in all EAGLE models considered here, with the only exception of the NOAGN run. Varying the stellar feedback efficiency, the equation of state imposed on the ISM or the SF law do not seem to alter this behaviour. As discussed in section 6, our findings suggest that AGN feedback is responsible for these trends. Only when AGN feedback is turned off, the metallicity anticorrelates with the SFR at the high-mass end, at least in these simulations.

7.1.2 AGN feedback effects

As discussed above, AGN feedback plays an important role in the regulation of the chemical enrichment of massive galaxies in EAGLE simulations. The flattening of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation at $M_{*}\ga 10^{10}\ M_{\sun}$, the decrease of effective yields as well as the reversal of the SFR-${\rm O/H}|_{\rm SF,gas}$ relation at high $M_{*}$ are partly driven by the influence of AGN. In particular, in section 3.3, we have shown that the lower metallicities of more massive galaxies ($M_{*}\ga 10^{10}\ {\rm M}_{\sun}$, corresponding to $M_{200}\ga 10^{11.5}\ {\rm M}_{\sun}$) can be associated to their lower SFRs, caused by the heating of the gas by AGNs. In this section, we will try to get more insight into how AGN affects the ISM and the metal content of simulated galaxies.

In order to explore if AGN feedback may have led to a net metal depletion in simulated galaxies, we studied how efficient AGNs are at ejecting gas from these systems. In the top left panel of Fig. 12, we analyse the ratio of the total baryonic component (i.e. $M_{\rm tot,bar}=M_{*}+M_{\rm SF,gas}+M_{\rm NSF,gas}$, where the different masses were calculated inside our standard aperture of 30 pkpc) of central galaxies relative to their host halo mass ($M_{\rm 200}$) as a function of halo mass. By analysing this quantity for different AGN models (see Table 1), we can obtain a rough estimate of how efficient AGNs were at ejecting mass from galaxies. From Fig. 12, it is clear that when AGN feedback is turned off (NOAGN model), a higher percentage of baryons are present in central galaxies. In the presence of AGN feedback, the amount of baryons in central galaxies is reduced by $\approx 50$ percent for most massive systems. This percentage increases with $\Delta T_{\rm AGN}$, as expected. As model BHs inject energy into neighbouring gas particles, these particles are located in the inner regions of galaxies and, so, they are expected to be metal enriched. Thus, our results suggest that AGN feedback is efficient at driving the ejection of metal-enriched material from simulated galaxies. Therefore, part of the decrease of the metallicity obtained at the high-mass end of the simulated $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation when AGN feedback is turned on, might be caused by the ejection of metal-enriched material. A more detailed investigation regarding this issue will be the subject of a future work.

In the following, we analyse how AGN feedback affects the state of baryons that remain inside simulated galaxies. In the top right panel of Fig. 12, we calculate the fraction of gas relative to the total amount of baryons inside these systems (i.e. $M_{\rm gas}/M_{\rm tot,bar}$, with $M_{\rm gas}=M_{\rm SF,gas}+M_{\rm NSF,gas}$ and different masses calculated inside our standard aperture of 30 pkpc) as a function of halo mass. We can see that, at a given $M_{200}$, different AGN models yield similar gas fractions. However, at a given $M_{200}$, the SFR of galaxies decreases as $\Delta T_{\rm AGN}$ increases (middle left panel), leading to lower $M_{*}$ at a fixed $M_{200}$ for higher $\Delta T_{\rm AGN}$ (middle right panel). Thus, since $M_{\rm gas}/M_{\rm tot,bar}$ is not significantly affected by AGN feedback, the decrease of $M_{*}$ caused by AGNs seems to have been compensated by the ejection of $M_{\rm gas}$ that they generate, as discussed above. In this context, the chemical enrichment of galaxies is quenched by AGN feedback by decreasing their SFR and by generating the ejection of metal-enriched gas from these systems. In the bottom left panel of Fig. 12, we analyse the percentage of NSF gas relative to the total gas component (inside an aperture of 30 pkpc) of galaxies for the different AGN models. Despite the large dispersion, there is evidence for an increase of the NSF-gas fraction at higher $\Delta T_{\rm AGN}$, as expected. This is consistent with the lower SFR obtained for higher $\Delta T_{\rm AGN}$. Thus, AGNs also play a role in heating the gas remaining in the galaxies, shutting-down the SFR activity and, thus, preventing further metallicity evolution from star formation.

Finally, the flattening of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation at the massive end might be partly caused by net metal dilution after AGNs suppress SF in these systems. In order to explore this scenario, we followed Yates & Kauffmann (2014) and used the parameter $Z_{\rm SF,gas}/Z_{*}$ as a diagnostic for dilution of the ISM after the last burst of star formation. $Z_{\rm SF,gas}$ and $Z_{*}$ were calculated inside our standard aperture of 30 pkpc. When dilution takes place there is a decrease in the gas-phase metallicity without a corresponding decrease in the stellar metallicity, hence, lower values of $Z_{\rm SF,gas}/Z_{*}$ should be obtained. In the bottom right panel of Fig. 12, we can see that the median $Z_{\rm SF,gas}/Z_{*}-M_{200}$ relation is not significantly affected by changes in the AGN feedback model. However, at high masses, the scatter in the $Z_{\rm SF,gas}/Z_{*}-M_{200}$ relation spreads down to negative values of $Z_{\rm SF,gas}/Z_{*}$ as $\Delta T_{\rm AGN}$ increases, with a clearly asymmetric dispersion about the median for at least the AGNdT9p0 simulation. These findings agree with those obtained by Yates & Kauffmann (2014) (e.g. their Fig. 6), who showed that dilution does not lead to a significant decrease of $Z_{\rm SF,gas}/Z_{*}$ for massive AGN host galaxies, but rather to an increase of the scatter in the $Z_{\rm SF,gas}/Z_{*}-M_{200}$ down to low (even negative) values at old ages or high masses. These results could suggest that dilution in AGN hosts also plays a (possibly minor) role in the lower metallicities seen in massive systems, alongside shutting­down star formation and metal­rich ejection.

8 Conclusions

We have studied ”fundamental” metallicity relations in the EAGLE suite of cosmological hydrodynamical simulations. These simulations were performed using a modified version of the GADGET-3 SPH code and includes prescriptions for different baryonic processes such as radiative cooling and heating, star formation, feedback from star formation and AGNs and abundance evolution of 11 elements (see section 2). We have analysed the evolution of metal abundances of galaxies as a function of mass and redshift and compared our results with different observational works. Given the discrepancies between observed metallicities obtained from different calibrators ($\la 0.7$ dex), we focused on the shape and level of evolution of the relations when comparing to simulations. We have focused mainly on the high-resolution Recal-L025N0752 simulation, which shows better agreement with the slope and normalization of certain observed metallicity scaling relations. We have also explored other simulations to evaluate resolution effects as well as variations in the subgrid physical parameters.

Our main findings can be summarized as follows:

• •

The high-resolution Recal-L025N0752 simulation predicts a correlation between star forming gas metallicity and stellar mass in agreement with the observed trend (Fig. 1). Metallicity tends to increase with $M_{*}$, exhibiting a shallower slope towards high stellar masses ($M_{*}\ga 10^{10}{\rm M}_{\sun}$). At a given $M_{*}$, metallicity tends to decrease with $z$, consistent with the observed trend, showing an overall variation of $\approx 0.5$ dex below $z\approx 3$ (Fig. 2). The dependence of energy feedback on the local density and metallicity of the ISM seems to drive the MZR evolution, at least in these simulations.

• •

The flattening of the slope of the simulated $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation at the high-mass end is mainly the result of AGN feedback (see also Segers et al. 2016a). To explore the role of AGN feedback, we used intermediate-resolution simulations as we needed to focus on the trends at high masses. We note that, in the overlap region at high masses, different resolution runs agree. Higher values of the temperature increase ($\Delta T_{\rm AGN}$) decrease the slope of the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ and $M_{200}-{\rm O/H}|_{\rm SF,gas}$ relations at $M_{*}\ga 10^{10}{\rm M}_{\sun}$ and $M_{200}\ga 10^{11.5}{\rm M}_{\sun}$, respectively (Figs. 3 and 4). Because AGN feedback prevents further star formation, the rate of metallicity evolution decreases as the heating temperature increases. A higher $\Delta T_{\rm AGN}$ also generates a larger scatter in the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relation at high masses.

• •

The simulated $M_{*}-Z_{*}$ relation is broadly consistent with observations below $z\approx 3$ (Fig. 5). In particular, at a given $M_{*}$, $Z_{*}$ tends to be lower than the ISM metallicity, especially for low-mass systems. This is expected as current star formation uses more metal-enriched material than was used to form the older stellar components of galaxies. For a study of the recycled stellar ejecta as fuel for star formation in EAGLE, the reader is referred to Segers et al. (2016a).

• •

At a given stellar mass, metallicities associated with the non-star forming (NSF) gas components are lower than stellar ones (Fig. 6). The ratio $Z_{\rm gas}/Z_{*}$ for the NSF gas-phase increases with decreasing $z$ and decreases with stellar mass, showing a stronger level of evolution at the high-mass end. In particular, massive galaxies show the most homogeneous metallicity distribution between the different baryonic phases.

• •

Simulated effective yields tend to increase with baryonic mass ($M_{\rm bar}=M_{*}+M_{\rm SF,gas}$) at low masses ($M_{\rm bar}\la 10^{10}\ {\rm M}_{\sun}$), consistent with the observed trend (Fig. 7). At higher $M_{\rm bar}$ and $z=0$, the predicted effective yields decrease with $M_{\rm bar}$, while observations suggest a flattening of the relation. Stronger AGN feedback leads to a larger scatter and lower average values of effective yields at higher masses. In addition, the simulated relation between effective yields and $M_{\rm bar}$ does not evolve significantly with $z$ at the low-mass end but, at $M_{\rm bar}\ga 10^{10}\ {\rm M}_{\sun}$, effective yields increase with $z$.

• •

At a given $M_{*}$, there are secondary dependencies of metallicity on the star-forming gas fraction, specific SFR and mass-weighted stellar age (Fig. 8). Galaxies with higher star-forming gas fractions, higher specific SFRs or lower ages are less metal enriched, with larger variations at lower stellar masses. At the high-mass end, the dependence of metallicity on star-forming gas fraction and specific SFR turns from an anti-correlation into a correlation due to the effects of AGN feedback. Simulations reproduce the observed fundamental metallicity relation between $M_{*}$, metallicity and SFR reported by Mannucci et al. (2010) quite well (Fig. 9). However, the simulated fundamental metallicity relation seems to be in place since $z\approx 5$ while, according to some observations, the observed one is established only at $z\leq 3$. This discrepancy might be partly due to selection biases in observational studies.

• •

In the case of EAGLE, the most fundamental metallicity scaling relation seems to be the correlation between oxygen abundance in star-forming gas and star-forming gas fraction (Fig. 10) as it presents low scatter and does not evolve significantly with $z$ within the scatter. Besides, it involves only two global properties of galaxies. This is consistent with observational findings by Bothwell et al. (2016b), who reported a dependence of metallicity on the molecular gas fraction of galaxies. Note also that Lagos et al. (2016) found a strong dependence of metallicity on neutral gas fractions from intermediate-resolution EAGLE simulations. By means of empirical-constrained analytical models, Zahid et al. (2014a) also reported a universal relation between metallicity and stellar-to-gas ratio.

• •

By comparing models with different subgrid physics parameters, we found that the slope of the mass-metallicity relation is mainly modulated by stellar feedback at low-masses ($M_{*}\la 10^{10}{\rm M}_{\sun}$), while AGN feedback regulates the slope at the high-mass end. Increasing the efficiency of either type of feedback strengthens the dependence of metallicity on SFR at a given stellar mass, increasing the scatter in the mass-metallicity relation. Increasing the power-law index $n$ in the star formation law seems to generate similar effects, at low stellar masses, to those obtained by reducing the stellar feedback efficiency. On the other hand, variations in the power law slope of the polytropic equation of state (${\gamma}_{\rm eos}$) imposed on the ISM do not affect the features of the mass-metallicity relation.

• •

Our findings suggest that AGN feedback generates a decrease in the global metallicity of SF gas by quenching the star-formation activity of simulated galaxies and also by leading to ejection of metal-enriched material. Net metal dilution in AGN hosts also seems to play a (possibly minor) role in the lower metallicities obtained for massive systems.

We have seen that the EAGLE simulation Recal-L025N0752 predicts metallicity scaling relations of galaxies as a function of stellar mass and redshift that are generally consistent with observations. It also reproduces secondary dependencies of metal abundances on SFR, specific SFR, star-forming gas fraction and stellar age at a fix stellar mass. However, the reality of the detailed features of the observed metallicity relations is still a matter of debate. The use of different metallicity calibrators, different methods for estimating stellar masses and gas fractions, selection biases, redshift-dependent instrumental apertures, etc., prevent observational works from converging to a consistent picture regarding the evolution of the mass-metallicity relation and the correlation between mass, metallicity and SFR. On the other hand, in the case of simulations, the uncertainties regarding the nucleosynthetic yields, resolution issues and the different subgrid prescriptions implemented by different authors make the comparison between different models a difficult task. In this context, continued comparison between observations and theoretical models of galaxy formation provide a fruitful tool to shed light on such a challenging topic. The sensitivity of the metallicity relations to the efficiencies of stellar and AGN feedback, as well as to the star formation law, make such efforts worthwhile.

Acknowledgements

We thank the referee for constructive remarks which improved the paper. We thank John Stott, Jabran Zahid and Sara Ellison for useful suggestions and comments. M.E.D.R is grateful to María Sanz and Guadalupe Lucia for their help and support during this project. We acknowledge support from PICT-2015-3125 of ANPCyT, PIP 112-201501-00447 of CONICET, UNLP G151 of UNLP (Argentina) and STFC consolidated and rolling grants ST/L00075X/1 (Durham, UK). This work was supported by the Netherlands Organisation for Scientific Research (NWO), through VICI grant 639.043.409, and the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013)/ERC Grant agreement 278594-GasAroundGalaxies. We acknowledge support from the European Commission’s Framework Programme 7, through the Marie Curie International Research Staff Exchange Scheme LACEGAL (PIRSES-GA-2010-269264). We acknowledge the Virgo Consortium for making their simulation data available. The EAGLE simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyr̀es-le-Châtel. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

Appendix A: Aperture effects

Throughout this work we have determined global properties of galaxies within 3-D aperture radius $R_{\rm ap}=30$ proper kpc. We have checked that changing $R_{\rm ap}$ does not affect the main trends found for metallicity scaling relations. However, changes in $R_{\rm ap}$ can generate moderate variations of the slope and normalization of those relations.

In Figs. 13 and 14 we study the impact of changing $R_{\rm ap}$ on the $M_{*}-{\rm O/H}|_{\rm SF,gas}$ and $M_{*}-Z_{*}$ relations, respectively. Results for Recal-L025N0752 are analysed at different $z$. It is clear that varying $R_{\rm ap}$ between 20 and 30 proper kpc yields negligible modifications in the MZR for all considered $z$. This suggests that most of the mass in these galaxies is located at radii $<20$ kpc. On the other hand, when decreasing $R_{\rm ap}$ towards 10 kpc, there is a slight overall increase of gas and stellar metallicity at all studied $z$ and, especially, at high stellar masses.

Thus, according to our results, simulated metallicity scaling relations tend to be convergent at $R_{\rm ap}\ga 20$ kpc while moderate changes in slope and normalization are obtained when $R_{\rm ap}$ is reduced to $\approx 10$ kpc, especially at high masses.

We have also analysed the evolution of the MZR by measuring metallicities inside the apertures reported by Maiolino et al. (2008) for observed galaxies, which are redshift-dependent. Such changes tend to slightly increase the slope of the simulated MZR, especially for massive galaxies at $z\approx 2$. Nevertheless, the flattening of the MZR at $z\sim 0$ is still present when applying the aforementioned apertures. In addition, when using those redshift-dependent apertures, the median simulated MZR obtained at a given $z$ is displaced by $\approx 0.1$ dex towards higher metallicities. However, the level of evolution of the MZR between $z=0$ and $z=3$ remains at $\approx 0.5$ dex.

Appendix B: Comparison with GIMIC simulations

The subgrid physics implemented in EAGLE is based on that developed for OWLS (Schaye et al., 2010), and applied also in GIMIC (Crain et al., 2009) and cosmo-OWLS (Le Brun et al., 2014). As commented before, here we extended the work of De Rossi et al. (2015b), who studied galaxy metallicity scaling relations by analysing GIMIC simulations. The most important differences between the GIMIC and EAGLE subgrid models concern the implementation of energy feedback from star formation (which is now thermal rather than kinetic and constant and depends on the local density and metallicity), the star formation law (which now depends on metallicity) and the inclusion of an AGN feedback model.

As the subgrid physics implemented in EAGLE is based on that applied in GIMIC simulations, it is worth a comparison between the MZR derived from both models. For more details about metallicity scaling relations in GIMIC simulations, the reader is referred to De Rossi et al. (2015b, a, 2016). In the following, we will compare results from high-resolution GIMIC simulations, high-resolution EAGLE simulation Recal-L025N0752 and intermediate-resolution EAGLE simulation FBconst-L050N0752 (see Table 1). We note that high-resolution GIMIC simulations have a resolution similar to that corresponding to EAGLE intermediate-resolution simulations. We have checked that similar conclusions to the ones presented below are obtained if reference intermediate-resolution EAGLE simulations are used instead of the Recal-L025N0752 run.

In the left panel of Fig. 15, we show GIMIC (thick lines) and EAGLE Recal-L025N0752 (thin lines) $M_{*}-{\rm O/H}|_{\rm SF,gas}$ relations at different $z$, as indicated in the figure. We will not concern about the normalization of the relations as it is not well constrained observationally (see Fig. 1). We can see that, at low masses, the slope predicted by GIMIC is slightly steeper than that obtained from EAGLE Recal-L025N0752. As can be inferred from Fig. 1, intermediate-resolution EAGLE simulations predict even flatter slopes than GIMIC for less massive galaxies. The EAGLE-GIMIC differences at the low-mass end of the MZR are consistent with a weaker net effect of SN feedback in EAGLE.

Towards higher masses and lower $z$, the slope of the EAGLE Recal-L025N0752 MZR flattens, in agreement with observations. There is evidence of a saturation of the EAGLE galaxies metallicities at an asymptotic value of $12+\log_{10}({\rm O/H})\approx 8.9-9.0$. However, there is no signature of this flattening in the case of GIMIC. De Rossi et al. (2015b) claimed that this last issue is related to the lack of an AGN feedback model in GIMIC. As discussed previously in section 3.3, our findings validate that conclusion as the AGN feedback model implemented in EAGLE seems to be responsible for the shallow slope of the MZR at the high-mass end (see also Crain et al., 2015; Segers et al., 2016a). AGN feedback prevents further star formation in massive galaxies and, hence, decelerates the metal enrichment process of massive systems. AGN feedback also can drive the ejection of metal-enriched gas. In addition, De Rossi et al. (2015b) obtained a very steep slope for the FMR at high stellar masses because of the lack of AGN feedback in the GIMIC simulations. We have seen that the EAGLE simulations, which include AGN, seem to have resolved that issue too (Section 6).

As reported by De Rossi et al. (2015b), Fig. 15 shows a negligible evolution of the GIMIC MZR between $z\approx 0$ and $z\approx 3$. The simulated evolution of the EAGLE Recal-L025N0752 MZR between $z\approx 0$ and $z\approx 3$ reaches $\approx 0.5$ dex, in better agreement with observed data (see Fig. 2). However, the weaker MZR evolution derived from EAGLE at $z\ga 2$ compared to some observations (e.g. Maiolino et al., 2008; Onodera et al., 2016) suggests that stellar feedback might still be inefficient to model metal-enriched outflows in low-mass galaxies at $z>2$. However, as there is not consensus about the observed level of evolution of the MZR yet, we are not able to validate or refute findings from simulations.

We found that the evolution of the MZR in EAGLE is driven by the dependence of the energy injection from core-collapse supernovae on the local density and metallicity of the ISM. We verified this by analysing the simulation FBconst-L050N0752, which implements a model that injects into the ISM a fixed quantity of energy per unit stellar mass formed, similarly to the model implemented in GIMIC simulations. Results can be seen in the right panel of Fig. 15. Clearly, when the injected energy is independent of local conditions, we obtain a negligible evolution for the EAGLE MZR. At low masses, the MZR slope predicted by the FBconst-L050N0752 simulation is also consistent with GIMIC simulations. As mentioned, at high masses, the discrepancies between EAGLE and GIMIC are due to the implementation of an AGN feedback in the former simulations.