Coevolution of metallicity and star formation in galaxies to z ≃ 3.7 – II. A theoretical model

Abstract

Recent work suggests that galaxy evolution, and the build-up of stellar mass (M_*) over cosmic time, is characterized by changes with redshift of star formation rate (SFR) and oxygen abundance (O/H). In a companion paper, we have compiled a large data set to study Metallicity Evolution and Galaxy Assembly (MEGA), consisting of ∼1000 galaxies to z ≃ 3.7 with a common O/H calibration. Here we interpret the MEGA scaling relations of M_*, SFR, and O/H with an updated version of the model presented by Dayal et al. This model successfully reproduces the observed O/H ratio of ∼80 000 galaxies selected from the Sloan Digital Sky Survey to within 0.05–0.06 dex. By extending the model to the higher redshift MEGA sample, we find that although the specific mass loading of outflows does not change measurably during the evolution, the accretion rate and gas content of galaxies increase significantly with redshift. These two effects can explain, either separately or possibly in tandem, the observed lower metal abundance of high-z galaxies.

galaxies: abundances, galaxies: evolution, galaxies: high-redshift, galaxies: star formation

1 INTRODUCTION

Galaxy evolution takes place through the build-up over time of stellar mass (M_*) through various episodes of star formation. As stars evolve, metals are produced, so that metallicity and its relation with M_* and star formation rate (SFR) are important gauges of star formation history (SFH). Thus, metal content, typically measured through the gas-phase oxygen abundance (O/H), the most abundant heavy element produced by massive stars, can be used to assess a galaxy's evolutionary state. The integrated star formation activity over the lifetime of the galaxy is measured by stellar mass, while the SFR measures how much gas is being converted into stars over the time-scale of the SFR indicator. The gas-phase metallicity (Z) not only quantifies the production of metals from high-mass stars, but also indicates the degree to which galaxies interact with their environment through gas accretion and outflows in the form of galactic winds.

Several well-established scaling relations reflect the intimate link between M_*, SFR, and O/H. A mass–metallicity relation (MZR) connecting stellar mass and Z is clearly present in the Local Universe (e.g. Tremonti et al. 2004), and it apparently extends to the highest redshifts examined so far (e.g. Erb et al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Cresci et al. 2012; Xia et al. 2012; Yabe et al. 2012; Henry et al. 2013; Cullen et al. 2014; Steidel et al. 2014; Troncoso et al. 2014; Wuyts et al. 2014; Zahid et al. 2014; de los Reyes et al. 2015; Ly et al. 2015; Onodera et al. 2016). Stellar mass and SFR are also related in a SF ‘main sequence’ (SFMS) both locally (Brinchmann et al. 2004; Salim et al. 2007) and at high redshift (e.g. Noeske et al. 2007; Elbaz et al. 2011; Karim et al. 2011; Wuyts et al. 2011; Speagle et al. 2014). It is generally agreed that the high-redshift relations of both the MZR and the SFMS are similar in form to the local ones, but they show different normalizations: at a given M_*, SFR (and sSFR) increases with increasing redshift while metallicity decreases. Moreover, both relations show an inflection at high M_*, with a decrease in O/H and SFR (Tremonti et al. 2004; Noeske et al. 2007; Wyder et al. 2007; Whitaker et al. 2014; Gavazzi et al. 2015; Lee et al. 2015).

In order to observationally constrain the evolution of metallicity with redshift, we have compiled a new data set of ∼1000 star-forming galaxies from z ≃ 0 to z ∼ 3.7 with nebular oxygen abundance measurements; we have dubbed this compilation the ‘MEGA’ data set, corresponding to Metallicity Evolution and Galaxy Assembly. This data set presents a marked improvement over the data set used by Hunt et al. (2012) because of the inclusion of several more high-z samples and, more importantly, because of a common metallicity calibration. The coevolution of SFR and O/H with redshift in the MEGA data set is presented in a companion paper (Hunt et al. 2016, hereafter Paper I), and summarized in Section 2. Here, we describe a physical model for understanding the coevolution of M_*, SFR, and O/H. This is an extension of the model presented by Dayal, Ferrara & Dunlop (2013), and also a re-assessment of the model parameters with different O/H calibrations, different stellar yields, and a different initial mass function (IMF). Section 3 describes this analytical approach to physically interpret the scaling relations locally, and Section 4 extends the analytical formulation to quantify evolution of the scaling relations with redshift. We discuss the results of our modelling in Section 5, and summarize our conclusions in Section 6. As in Paper I, we use a Chabrier (2003) IMF throughout.

2 METALLICITY AND SFR COEVOLUTION

Paper I presents common O/H calibrations for 990 galaxies from z ∼ 0 to z ∼ 3.7 culled from 19 different samples. The MEGA data set spans a range of 10⁵ in M_*, 10⁶ in SFR, and almost two orders of magnitude in O/H. The conversion to common metallicity calibrations was performed according to the transformations given by Kewley & Ellison (2008). We chose three O/H calibrations, because of their similarity to electron-temperature or ‘direct’ methods (e.g. Andrews & Martini 2013): Denicoló, Terlevich & Terlevich (2002, hereafter D02) and the nitrogen and oxygen+nitrogen-based abundances by Pettini & Pagel (2004, hereafter PP04N2 and PP04O3N2, respectively).

We hypothesized that the inflections in the MZR and SFMS at high M_* were compensating one another, and also that the decrease in O/H at high redshift was somehow quantitatively associated with the increase in SFR. This hypothesis was essentially aimed at assessing the possibility that the mutual correlations among M_*, SFR, and O/H could be expressed as a planar relation, i.e. with only two of the three pseudo-observables being required to describe the variation. From a Principal Component Analysis (PCA), we showed that this was indeed possible, and found a two-dimensional relation, namely the Fundamental Plane of Metallicity (FPZ):

\begin{equation} 12+\log ({\rm O/H}) = -0.14\,{\rm log(SFR)} + 0.37\,{\rm log(M_{\ast })} + 4.82. \ \ \end{equation}

(1)

The FPZ is an accurate reflection of the relations among M_*, SFR, and O/H at least to within 0.16 dex in O/H (equation 1 represents the best-behaved O/H calibration, PP04N2, Pettini & Pagel 2004). As described in Paper I, such a level of residuals is comparable to trends found in other galaxy samples but with a smaller range in M_*, SFR, and O/H.

We also analysed whether or not the FPZ varied over time, and found a significant correlation of the residuals with redshift, but with a residual standard error of 0.16 dex, the same as the 0.16 dex uncertainty of the FPZ itself. Thus, we concluded that the FPZ is approximately redshift invariant, since any redshift variation of the FPZ is within the noise of the current data.

As a comparison to the MEGA data set, here and in Paper I we also analyse the emission-line galaxies from the Sloan Digital Sky Survey (SDSS) defined by Mannucci et al. (2010, hereafter the SDSS10 sample). Although its parameter-space coverage is relatively limited, because of its large size this sample adds superb statistical power to our analysis, and also enables a comparison with Dayal et al. (2013) who analysed the same galaxies. As for the MEGA data set, we have converted (Kewley & Dopita 2002, hereafter KD02) the metallicities reported by Mannucci et al. (2010) to the D02, PP04N2, and PP04O3N2 calibrations, as described in Paper I.

3 COEVOLUTION: PHYSICAL INSIGHTS

Much theoretical effort has been devoted to understand the evolution of the MZR and its dependence, if any, on SFR and the surrounding environment (e.g. Peng et al. 2010; Davé, Finlator & Oppenheimer 2011; Peeples & Shankar 2011; Hopkins, Quataert & Murray 2012; Krumholz & Dekel 2012; Yates, Kauffmann & Guo 2012; Dayal, Ferrara & Dunlop 2013; Lilly et al. 2013; Schaye et al. 2015). Here, following the analytical formalism of Dayal et al. (2013), we present a physical basis for the observed coevolution of SFR and O/H. With a nebular oxygen abundance X ≡ M_O/M_g where M_O and M_g are the galaxy oxygen and total gas mass, respectively, SFR, O/H, and M_* are related through a set of evolutionary differential equations (see Dayal et al. 2013):

\begin{equation} \frac{{\rm d}M_{\rm star}}{{\rm d}t}\,\equiv \,\psi \,=\,\epsilon _{\ast }\,M_{\rm g} \end{equation}

(2)

\begin{equation} \frac{{\rm d}M_{\rm g}}{{\rm d}t}\,=\,-(1-R) \psi + (a-w) \psi \end{equation}

(3)

\begin{equation} M_{\rm g}\,\frac{{\rm d}X}{{\rm d}t}\,=\,y(1-R)\psi - a X \psi . \end{equation}

(4)

Here |$\epsilon _{{\ast }}^{-1}$| is the star formation time-scale and we have renamed SFR as ψ; we assume that the SFR is proportional to the gas mass, and that the infalling gas is much less metal enriched than the ambient ISM, X_i ≪ X. The two constants, R and y, represent the returned fraction from stars and the yield per stellar generation, respectively, and depend on the IMF; we have adopted (R, y) = (0.46, 0.05), taken from Vincenzo et al. (2016) using the stellar yields from Nomoto, Kobayashi & Tominaga (2013) assuming a Chabrier (2003) IMF.¹ The lower (upper) mass limit for these yields is 0.1 M_⊙ (100 M_⊙).

Like Dayal et al. (2013), we further assume that gas outflow through winds, w(M)ψ, is proportional to the SFR, where M is the total (dark+baryonic) galaxy mass. There is some observational evidence of such a proportionality, at least in the hot and/or ionized gas components of galaxies (e.g. Martin 2005; Veilleux, Cecil & Bland-Hawthorn 2005). For convenience in the analytical formulation, we also take gas accretion, a(M)ψ, to be proportional to SFR. Although there is little evidence to directly support such a proportionality, it may be that gas is replenished in galaxies by cooling of pre-existing ionized gas in expanding supershells or galactic fountains (e.g. Fraternali & Binney 2008; Hopkins et al. 2008; Fraternali & Tomassetti 2012).

The solution of these equations can be written as:

\begin{equation} X = \frac{y(1-R)}{a} \left( 1 - \mu ^{-\alpha } \right) \end{equation}

(5)

where μ ≡ M_g/M_g0 = ψ[ε_*M_g0]^{− 1}, and α ≡ a (R − 1 + a − w)⁻¹. The initial gas mass can be expressed as:

\begin{equation} M_{\rm g0}\,\simeq \,M_{\rm star}\,(1 + w -a) + M_{\rm g}\,=\,(\Omega _{\rm b}/\Omega _{\rm M})\,M. \end{equation}

(6)

Altogether there are five free parameters in the (z ≃ 0) model: the inverse of the star formation time-scale, ε_*; and the coefficients and power-law indices describing the mass dependence of infall through gas accretion (a), |$a(M_{\rm star})\,=\,a_{\rm coeff}\,(M_{\rm star}/M_0)^{a_{\rm pow}}$|⁠; and of outflow by galactic winds (w), |$w(M_{\rm star})\,=\,w_{\rm coeff}\,(M_{\rm star}/M_0)^{w_{\rm pow}}$|⁠.

The solution to these equations results in an implicit SFH assumed by our models, namely an exponentially declining one. Such a SFH is observationally supported by some studies (e.g. Eyles et al. 2007; Stark et al. 2009), but contested by others (e.g. Pacifici et al. 2015; Salmon et al. 2015). In our case, the exponential time-scale depends on halo mass (here assumed to be proportional to M_*). As described by Dayal et al. (2013), this means that galaxies of different masses build up their metal contents at different rates: galaxies start from low X, low M_*, and high ψ values and then move towards higher metallicities as their M_* increases. However, depending on their mass, they move at different velocities along the track: the most massive observed galaxies are very evolved with a large specific age (i.e. in units of the star formation time-scale) and move essentially along constant metallicity curves; smaller objects are younger and are still gently building up their metal content.

We can accommodate this model (hereafter the FPZ model) to the MEGA data set and the SDSS10 sample for the three O/H calibrations considered here by fitting for these five parameters. Fig. 1 illustrates the resulting best-fitting FPZ model overplotted on the binned SDSS10 data for all three O/H calibrations (D02, PP04N2, PP04O3N2). Since this is a local calibration, the FPZ model fit to the MEGA data set included only galaxies with z ≤ 0.1; Fig. 2 shows this best fit to the MEGA data set of the FPZ model. Table 1 reports the best-fitting parameters.

Figure 1.

SDSS10 galaxies: 12+log(O/H) plotted against (log of) SFR; the data points correspond to the average 12+log(O/H) (with the three O/H calibrations) for SDSS galaxies in a given SFR bin, and the error bars to the 1σ spread of the data. The curves show three FPZ models with different O/H calibrations (D02, PP04N2, PP04O3N2) as discussed in the text. The colour coding is by (log) M_*: 7.5–8 (magenta), 8–8.5 (purple), 8.5–9 (blue), 9–9.5 (cyan), 9.5–10 (green), 10–10.5 (orange), 10.5–11 (red), 11–11.5 (brown). The model and the data are in excellent agreement with a residual scatter of 0.05−0.07 dex, according to the O/H calibration (see Table 1). A colour version of this figure is available in the online version.

Figure 2.

MEGA galaxies: 12+log(O/H) plotted against (log of) SFR; the data points correspond to the average 12+log(O/H) (with the three O/H calibrations) for galaxies in a given SFR bin, and the error bars to the 1σ spread of the data. The curves show three FPZ models with different O/H calibrations (D02, PP04N2, PP04O3N2) as discussed in the text; here we show the FPZ best-fitting MEGA model, rather than the FPZ best-fitting SDSS model given in Fig. 1. As in Fig. 1, the colour coding is by (log) M_*: 7.5–8 (magenta), 8–8.5 (purple), 8.5–9 (blue), 9–9.5 (cyan), 9.5–10 (green), 10–10.5 (orange), 10.5–11 (red), 11–11.5 (brown). Only bins with more than two galaxies are shown; in some cases the standard deviations within the bins are smaller than the data point. The residual scatter of the model is 0.24−0.26 dex, according to the O/H calibration (see Table 1). A colour version of this figure is available in the online version.

Table 1.

FPZ model best-fitting parameters.^a

Sample	Degrees of	|$\sigma _{\rm fit}^{\mathrm{b}}$|	O/H	acoeff	apow	wcoeff	wpow	ε*
	freedom		calibrationc					(Gyr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Local Universe
MEGAd	252	0.264	D02	0.47 ± 0.45	0.22 ± 0.60	15.55 ± 4.0	−0.34 ± 0.04	2.0 ± 1.3	−	−
MEGAd	252	0.243	PP04N2	0.88 ± 0.32	−0.27 ± 0.06	9.9 ± 4.7	−0.29 ± 0.06	1.1 ± 0.8	−	−
MEGAd	252	0.251	PP04O3N2	0.73 ± 0.29	−0.30 ± 0.06	8.6 ± 3.9	−0.32 ± 0.06	0.9 ± 0.6	−	−
SDSS10	78536	0.052	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	−	−
SDSS10	78536	0.063	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	−	−
SDSS10	78536	0.068	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	−	−
z > 0: |$\epsilon _{{\ast }}(z)\,=\epsilon _{{\ast }}(0)\,({M}_{\rm star}/{10^{9.5}})^{\gamma (\epsilon _{\ast })} \ (1+z)^{\delta (\epsilon _{\ast })}$|
									γ(ε*)	δ(ε*)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
MEGAe	25	0.060	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	−0.23 ± 0.07	−0.16 ± 0.09
MEGAe	25	0.038	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	−0.49 ± 0.05	−0.17 ± 0.08
MEGAe	25	0.046	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	−0.44 ± 0.05	−0.10 ± 0.08
Schechter f	28	0.079	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	−0.40 ± 0.06	−0.68 ± 0.10
Schechter f	28	0.025	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	−0.54 ± 0.03	−0.23 ± 0.05
Schechter f	28	0.025	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	−0.56 ± 0.03	−0.18 ± 0.05
z > 0: a(z) = a(0) (1 + z)γ(a), w(z) = w(0) (1 + z)δ(w)
									γ(a)	δ(w)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
MEGAe	25	0.046	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	0.46 ± 0.06	0.05 ± 0.06
MEGAe	25	0.039	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	0.68 ± 0.04	−0.07 ± 0.04
MEGAe	25	0.049	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	0.71 ± 0.05	−0.08 ± 0.06
Schechter f	28	0.035	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	0.48 ± 0.04	0.00 ± 0.04
Schechter f	28	0.042	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	0.73 ± 0.04	−0.19 ± 0.04
Schechter f	28	0.049	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	0.75 ± 0.05	−0.19 ± 0.05
z > 0: μ(z) = μ(0) [ (Mstar/109.5)γ(μ) (1 + z)δ(μ) ]−1
									γ(μ)	δ(μ)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
MEGAe	25	0.062	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	−0.16 ± 0.05	−0.10 ± 0.08
MEGAe	25	0.039	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	−0.44 ± 0.05	−0.14 ± 0.08
MEGAe	25	0.048	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	−0.36 ± 0.05	−0.05 ± 0.08
Schechter f	28	0.035	D02	0.994 ± 0.001	0.036 ± 0.002	11.25 ± 0.11	−0.449 ± 0.002	6.3 ± 0.2	−0.19 ± 0.03	−0.25 ± 0.04
Schechter f	28	0.026	PP04N2	1.136 ± 0.005	0.104 ± 0.005	20.55 ± 0.26	−0.413 ± 0.003	13.1 ± 0.8	−0.50 ± 0.03	−0.20 ± 0.05
Schechter f	28	0.032	PP04O3N2	1.020 ± 0.004	0.096 ± 0.005	18.09 ± 0.21	−0.451 ± 0.003	8.2 ± 0.4	−0.42 ± 0.03	−0.13 ± 0.05

^a The free parameters in the z ∼ 0 FPZ models are: |$a\,=\,a_{\rm coeff}\,({M_{\rm star}}/{10^{10.75}})^{a_{\rm pow}}$|⁠; |$w\,=\,w_{\rm coeff}\,({M_{\rm star}}/{10^{9.0}})^{w_{\rm pow}}$|⁠; ε_*. At z > 0, there are three classes of best-fitting models as described in the text: one that depends on the variation of the inverse of the SF time-scale (ε_*); one that depends on the variations of the inflow (a) and outflow (w) amplitudes; and one that depends on μ, the relative gas content.

^b Residual standard error of global fit.

^cD02 corresponds to D02; PP04N2 and PP04O3N2 to nitrogen-based and oxygen+nitrogen-based calibrations by Pettini & Pagel (2004).

^d These fits were obtained for the MEGA galaxies with z < 0.1.

^e We have adopted the SDSS10 z ∼ 0 best-fitting model parameters at z ∼ 0 in order to fit FPZ(z > 0); only galaxies with z > 0.1 are included in the binned data.

^f These are fits of the model to mean M_* and mean SFR using integrals of double Schechter functions taken from Ilbert et al. (2013), and assuming that sSFR varies with redshift as the star-forming galaxies observed by Karim et al. (2011); there are 28 dof in these fits (six redshift bins × five mass bins, since we do not use the lowest-mass bin, and two of the mass bins are empty at the given redshift bin). The redshift bins and the metallicity behaviour with redshift to emulate are taken from mean trends of the MEGA data set. We have adopted the SDSS10 z ∼ 0 best-fitting model parameters at z ∼ 0 in order to fit FPZ(z > 0). See text for more details.

The model is an excellent representation of the SDSS10 data (see also Dayal et al. 2013); over ∼80 000 individual data points, the best-fitting model with only five free parameters gives a residual error in 12+log(O/H) predicted versus observed of 0.05–0.07 dex, according to the O/H calibration. The model is somewhat less successful at fitting the MEGA data set, with mean O/H residuals of predicted versus observed of ∼0.24–0.26 dex, larger than the FPZ residuals of ∼0.16 dex (see Paper I). However, these galaxies are only a subset of the entire MEGA sample (because of the redshift limitation), and the best-fitting MEGA parameters are generally consistent, within the errors, to SDSS (see below for more details).

Despite the differences (O/H calibrations, stellar yields, IMF), our new best-fitting FPZ model parameters are generally similar to those of Dayal et al. (2013). The small M_* dependence of the accretion term and its small amplitude are consistent with previous values; for the galactic winds, however, we find a slightly larger coefficient and a slightly steeper M_* dependence (∼− 0.4 versus ∼− 0.33). The inverse SF time-scales |$\epsilon _{{\ast }}^{-1}$| are also smaller than in Dayal et al. (2013), corresponding to ∼100 Myr versus ∼600 Myr previously. Part of the differences in the best-fitting parameters may be due to the choice of the model IMF and our assumptions for the O/H calibration; here we use the Asplund et al. (2009) values of 12+log(O/H) = 8.69 for solar oxygen abundance, and a metal mass fraction Z/H of 0.0198.² Nevertheless, despite this and the difference in the fitting methods, the results are in reasonably good agreement.

The mass loading factor for galactic winds, w, found by the SDSS10 best-fitting model is roughly consistent with those found by the cosmological hydrodynamic simulations of Davé et al. (2011), despite their different underlying assumptions about the relation of SFR, gas accretion, and galactic winds. They assume an ‘equilibrium’ model (further described below, Section 3.1) in which the infall rate a is balanced by the sum of the SFR ψ and the outflow w. The results of their simulations show a trend with stellar mass of w|$\propto {M}_{\rm star}^{-0.33}$|⁠, or momentum-driven winds, consistent with the MEGA best fit, and only slightly flatter than the power-law index a_pow ∼ −0.4 given by the SDSS10 best fit. This power-law M_* dependence of ∼− 0.4 corresponds roughly to the low-mass linear portion of SDSS10 MZR with a slope of ∼0.38 (Paper I).

In contrast, the accretion or infall parameter a is significantly different from what is found in some equilibrium models. While the infall amplitude a is required to be larger than the outflow w in such models, a = (w + 1) ψ (see e.g. Davé et al. 2011), we find a much smaller a value, roughly a factor of 10 smaller than w, and with a much shallower M_* dependence. These differences will be discussed further in Section 3.1. In any case, our model is able to reproduce O/H of the ∼80 000 SDSS10 galaxies to within ∼ 0.05–0.07 dex, lending credibility to our approach.

The quantity |$\epsilon _{{\ast }}^{-1}$| corresponds approximately to a gas depletion time, with some caveats. Observationally, gas depletion times, τ_depl, are inferred from the measured ratio of SFR and gas mass, τ_depl = ψ/M_g, either locally as surface densities, or globally integrated over the entire galaxy (e.g. Bigiel et al. 2011; Saintonge et al. 2011; Huang & Kauffmann 2014; Genzel et al. 2015; Hunt et al. 2015). For spiral galaxies, typical depletion times τ_depl for the H₂ component are ∼2.4 Gyr (e.g. Bigiel et al. 2011) but can be as small as 50 Myr in metal-poor starbursts (Hunt et al. 2015) or as large as ∼10 Gyr in more quiescent systems (e.g. Saintonge et al. 2011; Bothwell et al. 2014). On the other hand, atomic gas H i depletion times are relatively constant, ∼3.4 Gyr (Schiminovich et al. 2010). The values we find for |$\epsilon _{{\ast }}^{-1}$| (∼ 100–200 Myr for SDSS10, ∼ 400–700 Myr for MEGA) are shorter than either of these.

The reason for this disagreement can be understood as follows: the constant of proportionality ε_* relating gas mass and SFR (see equation 2) is in reality the ratio of a (dimensionless) star formation efficiency ε_ff and a SF time-scale, t_ff (e.g. Krumholz, McKee & Tumlinson 2009; Krumholz, Dekel & McKee 2012), where the subscript ‘ff’ refers to a free-fall or dynamical time-scale. Physically, ε_ff gives the fraction of gas converted into stars over a dynamical time-scale, in a process that is usually highly inefficient with ε_ff∼ 0.01–0.05 (e.g. Krumholz & Tan 2007). With such efficiencies and assuming that ε_* = ε_ff/t_ff the SDSS10 fitted values of |$\epsilon _{{\ast }}^{-1}$| would give t_ff∼ 1–8 Myr. These times are unrealistically small for typical giant molecular clouds or spiral discs, which instead have typical dynamical times t_ff∼15–20 Myr (Krumholz et al. 2012).

However, what we call M_g is in fact only the gas that enriches the ISM with metals as it is recycled through star formation; it is not the total gas available that would be observed. In other words, M_g(model) = M_g(observed) ε_ff, or from equation (2):

\begin{equation} \epsilon _{{\ast }}\,=\,\frac{\psi }{M_{\rm g}}\,=\frac{\psi }{M_{\rm g}{\rm (observed)}\,\epsilon _{\rm ff}}\,=\,\frac{1}{\tau _{\rm depl}\,\epsilon _{\rm ff}}. \end{equation}

(7)

We can use this equation to calculate typical gas depletion times as estimated from our fits of the SDSS10 galaxies. Assuming ε_ff∼0.03, we find t_ff∼3 Gyr, similar to typical gas depletion times for both the molecular and the atomic component.³

Comparison of the MEGA results with those of SDSS10 shows that the amplitude and power-law index for the wind parameters are consistent with SDSS10, but the accretion parameters and the inverse SF time-scales, ε_*, differ. As can be seen in Fig. 2, the high-mass end of the SFR dependency of O/H falls off more steeply for the MEGA galaxies than for SDSS10. For the PP04N2 and PP04O3N2 calibrations, gas accretion for the MEGA galaxies seems to increase with decreasing M_*, in contrast to the behaviour of the SDSS10 sample; because of the highly star-forming nature of the MEGA data set, this could be telling us that the lower-mass galaxies are more gas rich than higher-mass ones. On the other hand, the SF time-scales implied by the FPZ model fit to the MEGA galaxies are longer than for SDSS10, although with relatively large uncertainties. More generally, the SDSS10 best-fitting parameters are in all cases much better determined than those for the MEGA data set, probably because of the much larger number of galaxies in the former. Therefore, in what follows, we will adopt the best-fitting FPZ parameters for the SDSS10 z ≃ 0 galaxies to extend the model to high redshift. Table 1 gives the best-fitting FPZ model parameters for the SDSS10 and MEGA data sets.

3.1 Comparison with previous work at z ≃ 0

As discussed in Paper I, the intention of the FPZ is similar to that of the ‘Fundamental Metallicity Relation’ (‘FMR’; Mannucci et al. 2010), namely to relate the metal content of galaxies to their star formation activity (see also Ellison et al. 2008). Nevertheless, the two formulations are different; the planar formulation of the FPZ, accurate to ∼0.16 dex, extends to galaxies with M_* as low as ∼10⁶ M_⊙ and to redshifts ≲ 3.7 (see Section 2 and Paper I). In contrast, the extension to lower stellar masses of the FMR is quadratic (Mannucci, Salvaterra & Campisi 2011), and, as shown in Paper I, has a larger mean dispersion than the FPZ and a significant offset for the MEGA sample which spans a much broader parameter space than the original SDSS data set at z ∼ 0 for which the FMR was developed.

Due to its success, which by including SFR reduced significantly the dispersion for 12+log(O/H) of the original data set relative to only stellar mass, many models have focused on reproducing the FMR at z ≃ 0 by Mannucci et al. (2010). Because of the underlying notion that metallicity is somehow connected to SFR, common to both the FPZ and the FMR formulations, before pursuing our model for the FPZ, here we discuss previous attempts to theoretically understand the FMR.

Despite using simulation runs calibrated to a number z ≃ 0 observables (stellar mass function, stellar mass–size relation and the stellar mass–black hole relation) the current gold-standard EAGLE simulations are unable to reproduce the local FMR; whilst reproducing observations to within 0.15 dex for high-mass (M_* ≲ 10¹⁰ M_⊙) galaxies, they severely over-estimate the metal content of lower mass haloes by as much as 0.4 dex (Lagos et al. 2016). Given the sub-grid limitations of simulation, most effort has thus been diverted to developing an analytic understanding of the z ≃ 0 FMR.

We briefly discuss the two previous works most close in spirit to ours. The Davé, Finlator & Oppenheimer (2012) model assumes every galaxy to be in ‘equilibrium’ with the infall rate exactly balancing the gas lost in outflows and star formation. While this model correctly yields a metallicity that is independent of star formation for high-mass galaxies, it fails to capture the metallicity downturn observed at low mass. Interestingly, however, this model finds outflows to be momentum driven with a mass power-law dependence of ∼− 0.33, which however we obtain only for the MEGA sample; the SDSS10 slopes are slightly steeper than this, ∼− 0.4. Peeples & Shankar (2011) find that reproducing the local MZR requires a high metal-expulsion efficiency (compared to expelled gas, ‘preferentially metal-enriched winds’) that scales steeply with increasing M_*. Although this might yield the observed FMR (or FPZ) for low-mass galaxies, it is doubtful if such an assumption could reproduce the constant metallicity observed for the most massive systems, independently of SFRs ranging over five orders of magnitude presented in the MEGA sample.

Lilly et al. (2013) use the ‘gas-regulator’ model wherein stars are embedded in a gas reservoir whose mass is increased by infall, and the SFR being proportional to the gas mass available at any time. Assuming both the star formation efficiency and outflow rate to be independent of the halo mass (the ‘ideal regulator’), these authors find the metallicity to be largely independent of the evolutionary path because of gas continually ‘flushing’ the system. This is quite contrary to our model that requires outflows to scale with the halo mass, and finds the metallicity to explicitly depend on the gas evolution history.

3.2 Testing the z ≃ 0 model at higher redshift

In order to extend the FPZ model to z > 0, we first test whether the SDSS10 best-fitting z ≃ 0 parameters apply to the data at higher redshift. Such behaviour would be expected if the lower metallicities at high redshift were compensated for merely by the more intense star formation activity, which could power larger outflow of metals through stronger galactic winds. The left panel of Fig. 3 shows the behaviour of the MEGA sample if we infer metallicity through the FPZ(z ≃ 0) model parameters. The curves are the fits to the observed O/H of the MEGA data set as a function of redshift as reported in Paper I; the symbols correspond to the values of 12+log(O/H) that would be inferred from the FPZ z = 0 model applied to the median SFR and M_* in each redshift bin. In all M_* bins, up to z ∼ 0.7 the FPZ correctly predicts the observed O/H. However, for higher redshifts, given the observed (median) SFRs, the FPZ model overpredicts the observed O/H, in particular for the highest M_* mass bins where the discrepancy is ≳ 0.4 dex. On the other hand, in the lowest mass bin [log(M_*/M_⊙) =8.5–9.0] O/H is underestimated by ∼0.1 dex. Similar tensions relative to the z = 0 model by Dayal et al. (2013) are seen also in the sample of z ∼ 2 galaxies by Grasshorn Gebhardt et al. (2016) where they find that galaxies with M_* ≳ 10⁹ M_⊙ have higher O/H relative to their local counterparts, while lower-mass systems have lower abundances.

Figure 3.

Binned measurements of 12+log(O/H) estimated from SFR and the FPZ model at z ≃ 0 as a function of redshift for the MEGA data set (left panel) and the simulated data set (right). In both panels, the symbols correspond to the values of 12+log(O/H) that would be inferred using the FPZ(z ≃ 0) model applied to the median SFR and M_* within each redshift bin (according to binned M_* values); the curves show the average behaviour of the MEGA data set as reported in Paper I. See text for a description of the simulated galaxy populations shown in the right panel. A colour version of this figure is available in the online version.

Because the high SFRs of the MEGA data set may not be fully representative of typical high-z ‘main-sequence’ galaxy populations (see Paper I), as a check, we have simulated truly main-sequence galaxy populations as a function of redshift using the Galaxy Stellar Mass Functions (GSMFs) by Ilbert et al. (2013) for z > 0 and for z ≃ 0 the GSMF from Baldry et al. (2012). For a given M_*, the SFR has been inferred by using the SFMS formulation of Speagle et al. (2014) as a function of cosmic time. For each redshift bin, integration over the appropriate weighted GSMFs was performed to obtain the mean M_* and the mean SFR at that redshift for main-sequence populations. Using these values of M_* and SFR, we then calculated O/H using the FPZ(z ≃ 0) model parameters as for the MEGA data. The result is shown in the right panel of Fig. 3, where the metallicities predicted from the FPZ(z ≃ 0) model are compared to the mean observed O/H behaviour (shown by curves) of the MEGA data set as reported in Paper I. As in the left panel, symbols correspond to the metallicities that would be inferred for the simulated galaxies of a given M_* and SFR(M_*,z) from the FPZ(z ≃ 0) model within each M_* and redshift bin. Despite the two very different approaches, both the simulated galaxies (right panel of Fig. 3) and the MEGA galaxies (left panel) show similar behaviour. The FPZ(z ≃ 0) parameters are successful to z ∼ 0.7, but fail in the same way as for the MEGA data set at higher z. The consistency of the simulations and the MEGA galaxies is encouraging, and in the next section we use both methods to extend the FPZ model to high redshift.

4 EXTENDING THE MODEL TO HIGHER REDSHIFT

As illustrated in the previous section, the higher SFRs at z > 0 are insufficient by themselves through the FPZ(z ≃ 0) model to lower the metallicities to the levels observed. We have thus investigated three avenues of adapting the FPZ model to z > 0: (a) changes in ε_*, since time-scales and/or star formation efficiencies might be expected to change with redshift (Section 4.1); (b) redshift variations of accretion and wind parameters (a and w) (Section 4.2); and (c) higher gas fractions through possible changes in the μ parameter (Section 4.3, see equation 5). Since the discrepancies in Fig. 3 seem to depend on M_*, we considered separate dependencies on M_* and redshift (or cosmic time). Thus, for each approach to establish the z-dependent FPZ model, we fix FPZ(z ≃ 0) to the SDSS10 parameters and introduce a scaling factor for either ε_*, (a, w), or μ. Thus only two parameters are to be fit: one for a possible M_* scaling for z > 0 (ε_*, μ) and one for a scaling with redshift (ε_*, μ, a, w). Although these are only a few of the many possible adaptations of the FPZ model to high z, we have limited the possibilities for simplicity.

To constrain the model, we have assumed that the true metallicities of high-z galaxy populations are approximated by the mean behaviour of the MEGA data set as reported in Paper I. The best-fitting model parameters for z > 0 are obtained by minimizing the residuals of the observed mean behaviour of 12+log(O/H) relative to the model over the various mass and redshift bins, not including z ≃ 0; because of small numbers at high redshift, we also do not consider the lowest mass bin [log(M_*/M_⊙)<8.5]. The resulting degrees of freedom in the fits are 25 for the MEGA data set (since some M_* bins are missing at z > 0), and 28 for the simulated galaxies (five M_* bins over six redshift bins, two free parameters). As mentioned above, in all high-z FPZ models, we adopted the z ≃ 0 SDSS10 parameters and adjusted ε_* (first approach), a and w, the accretion and wind coefficients (second), or μ (third). Table 1 gives the best-fitting FPZ(z) (D02, PP04N2, PP04O3N2) model parameters for the three approaches which are described in the following sections.

4.1 Redshift variation of star formation time-scales and/or efficiencies

We first investigated whether scaling ε_* with M_* and z could have the desired effect of lowering O/H at high redshift. Several methods were explored for introducing a scaling factor e_scale applied as a multiplicative factor to ε_* in the FPZ models. Although we studied numerous ways to scale ε_*, including a mass-quenching scenario as formulated by Peng et al. (2010), only two gave reasonably low mean O/H residuals, ∼0.03 − 0.04 dex for both MEGA and simulated galaxies. These correspond to ε_*(z) = e_scale ε_*(0) where e_scale is given by:

• (M_*/|$10^{9.5})^{\gamma (\epsilon _{\ast })}$||$(1+z)^{\delta (\epsilon _{\ast })}$|

• (M_*/|$10^{9.5})^{\gamma (\epsilon _{\ast })}$||$({\rm age}_z/{\rm age}_{\rm H})^{\delta (\epsilon _{\ast })}$|

where ε_*(0) gives the value of ε_* at z = 0, age_z and age_H are the ages of the universe at redshift z and z = 0, respectively, and γ(ε_*) and δ(ε_*) are the parameters to be fit. Since the two fits are equivalent in terms of quality (low residuals), we chose (a) since it is generally easier to formulate numerically.

The values of 12+log(O/H) predicted by this model as a function of redshift are shown in Fig. 4; as in Fig. 3 the left panel gives the MEGA data set and the right panel the simulated galaxies. The fit is excellent for the simulated galaxies, with a mean residual σ ∼ 0.025 dex (except for D02), but for the MEGA data set it is slightly worse, σ = 0.04–0.06 dex. There are some problems in the intermediate-z regime where the metallicities inferred from the FPZ model by varying ε_* with mass and z tend to be overestimated. The predictions for low-mass galaxies also tend to be overestimated at low redshift.

Figure 4.

Binned measurements of 12+log(O/H) estimated from SFR and the FPZ(z) model with ε_*(z) added to the FMZ(z = 0) model as a function of redshift for the MEGA data set (left panel) and the simulated data set (right). In both panels, the symbols correspond to the (PP04N2) values of 12+log(O/H) that would be inferred using the FPZ(z ≃ 0) model applied to the median SFR and M_* within each redshift bin (according to binned M_* values); the curves show the average behaviour of the MEGA data set as reported in Paper I. See text for a description of the simulated galaxy populations shown in the right panel. A colour version of this figure is available in the online version.

4.2 Redshift variation of gas accretion and galactic winds

The second avenue of investigation relied on simultaneously scaling the gas accretion coefficient a and the galactic wind coefficient w, in order to quantify whether the scaling with SFR changes with redshift. For simplicity, we used one coefficient γ(a) to define the redshift variation of a, and another δ(w) for that of w:

\begin{eqnarray*} a(z)&=&a(0)\,(1+z)^{\gamma (a)} \nonumber \\ w(z)&=&w(0)\,(1+z)^{\delta (w)} \nonumber \end{eqnarray*}

where a(0) and w(0) are the values of a_coeff and w_coeff at z = 0. The implicit assumption is that the M_* dependence of a(0) and w(0) given by the SDSS10 fit in Table 1 holds also at z > 0.

Fig. 5 gives the results of this fit, where, as before, the left panel shows the MEGA data set and the right the simulated galaxies. The fit for the MEGA data set is good, with a mean residual between data and predicted O/H of σ = 0.04–0.05 dex, similar to the previous ε_* approach described above. For the simulated galaxies, the fit is still good (σ = 0.04–0.05 dex), although slightly worse than the ε_* model (except for D02 where it is better).

Figure 5.

Binned measurements of 12+log(O/H) estimated from SFR and the FPZ(z) model with a(z), w(z) added to the FMZ(z = 0) model as a function of redshift for the MEGA data set (left panel) and the simulated data set (right). In both panels, the symbols correspond to the (PP04N2) values of 12+log(O/H) that would be inferred using the FPZ(z ≃ 0) model applied to the median SFR and M_* within each redshift bin (according to binned M_* values); the curves show the average behaviour of the MEGA data set as reported in Paper I. See text for a description of the simulated galaxy populations shown in the right panel. A colour version of this figure is available in the online version.

4.3 Redshift variation of gas fraction

The last successful formulation involved scaling the model parameter μ which basically defines the gas-mass fraction relative to its initial value in the galaxy's evolution (see equations 5 and 6). For computational reasons, the scaling factor was introduced in the denominator of μ:

\begin{equation*} \mu (z)\,=\,\mu (0)\ [\,({M}_{\rm star}/10^{9.5})^{\gamma (\mu )}\,(1+z)^{\delta (\mu )})\,]^{-1}. \nonumber \end{equation*}

The FPZ(z) model predictions of O/H with this approach are shown in Fig. 6 where the MEGA data set is shown in the left panel and the simulated galaxies in the right. Again, the fit for the MEGA data set is good, with mean residuals σ = 0.04–0.06 dex, consistent with the previous two formulations. The behaviour of the MEGA galaxies is similar to the ε_* approach described in Section 4.1 where O/H is over-predicted by the model in the intermediate-z redshift range, although the discrepancies at the low-mass end are slightly less pronounced. On the other hand, this formulation of the FPZ(z) model for the simulated populations is excellent, with a mean residual of σ ∼ 0.03 dex.

Figure 6.

Binned measurements of 12+log(O/H) estimated from SFR and the FPZ(z) model with M_g(z) added to the FMZ(z = 0) model as a function of redshift for the MEGA data set (left panel) and the simulated data set (right). In both panels, the symbols correspond to the (PP04N2) values of 12+log(O/H) that would be inferred using the FPZ(z ≃ 0) model applied to the median SFR and M_* within each redshift bin (according to binned M_* values); the curves show the average behaviour of the MEGA data set as reported in Paper I. See text for a description of the simulated galaxy populations shown in the right panel. A colour version of this figure is available in the online version.

5 DISCUSSION

This is one of the first works that aims to study the physics governing the redshift-evolution of the FPZ (or FMR). Lilly et al. (2013) predicted an un-evolving FPZ unless the parameters governing inflow/outflow/star formation evolve with redshift; our approach suggests that indeed these parameters evolve. Yates et al. (2012) study the evolution of the MZR from z ∼ 3 to z ≃ 0 by comparing the semi-analytic models of Guo et al. (2011) with an SDSS-selected sample similar to SDSS10. However, their models were unable to predict the observed decrease of metallicity for a fixed M_*, and they concluded that ‘chemical enrichment in the model galaxies proceeded too rapidly at early times’.

Here, we have shown that to correctly reproduce the coevolution of SFR and O/H, i.e. correctly fit the observations, the model describing the co-dependence of M_*, SFR, and metallicity must change with redshift. All three FPZ(z) fits to both the MEGA data set and the simulated galaxies as described above successfully predict the O/H trends with redshift to within 0.03–0.05 dex. Although the three formulations are apparently describing different physical manifestations of the coevolution of metallicity and SFR, they are generally painting a similar picture, namely an increase of gas content with redshift.

5.1 Redshift variation of the model

According to Table 1, the scaling of the FPZ(z) model with ε_* suggests that:

\begin{equation} \epsilon _{\ast }(z)\,\approx \,\epsilon _{\ast }(0)\,({M}_{\rm star}/10^{9.5})^{-0.5}\,(1+z)^{-0.2}. \end{equation}

(8)

The above γ(ε_*), δ(ε_*) coefficients are approximate, taken as a rough average of the (PP04N2) simulated and MEGA data sets. Equation (8) would imply that at a given M_*, |$\epsilon _{{\ast }}^{-1}$|⁠, roughly the product of gas depletion time τ_depl and SF efficiency ε_ffis increasing slightly with redshift (roughly 30 per cent larger at z ∼ 3). Assuming that ε_ff remains constant, this would be contrary to observational evidence that suggests that molecular depletion times are shorter at higher redshift (e.g. Genzel et al. 2010, 2015; Silverman et al. 2015). Galaxies are also more gas-rich at higher redshift, with as much as 40–60 per cent of their dynamical mass in H₂ (e.g. Daddi et al. 2010; Tacconi et al. 2010), a fraction ∼5–10 times higher than typical spiral discs at z ≃ 0. However, as discussed in Section 3, depletion times are defined observationally as the ratio of observed gas mass to SFR. Our models do not know about observed gas mass; they only consider the gas that is enriching the ISM with metals, i.e. the gas that is effectively converted into stars. Thus, of necessity, our models cannot distinguish between the observationally defined τ_depl and the product of τ_depl and ε_ff, or equivalently the ratio of t_ff and ε_ff. Given that |$\epsilon _{{\ast }}^{-1}$| = M_g(observed) ε_ff/ψ (see equation 7) it is likely that the ε_* FPZ model scaling is pointing to increasing gas content with z, rather than longer depletion times.

The FPZ(z) formulation with the accretion and wind coefficients is consistent with this. As shown in Table 1, the scaling of the model with a and w:

\begin{eqnarray} a(z)&\approx &a(0)\,(1+z)^{0.7}, \nonumber \\ w(z)&\approx &w(0)\,(1+z)^{0.0}, \end{eqnarray}

(9)

indicates that the mass loading relative to SFR in the galactic outflows powered by star formation is roughly the same at high z as at z ≃ 0. In contrast, the mass loading in accretion is significantly increased, almost three times higher at z ∼ 3.5 than locally. Thus, to explain the drop in metallicity at high redshift, more (pristine) gas is needed, in the model acquired through accretion. Our result is roughly consistent with that of Papovich et al. (2011) who find that gas fractions increase roughly as (1 + z)^0.9 from z ∼ 3 to z ∼ 8, and after this ‘accretion epoch’, at lower redshifts both accretion and SFR are reduced.

A similar indication is given by the FPZ(z) model based on the redshift variation of μ, equal to the ratio of M_g and the initial gas mass (see equation 5). The coefficients in Table 1 indicate that:

\begin{equation} \mu (z)\,\approx \,\mu (0)\,\frac{1}{({M}_{\rm star}/10^{9.5})^{-0.4}\,(1+z)^{-0.15}}.\quad \end{equation}

(10)

Both the M_* and redshift dependencies are similar to (although slightly smaller than) those for the ε_* formulation; this is not unexpected given that |$\mu \,=\,M_{\rm g}/M_{\rm g0}\,\propto \epsilon _{\ast }^{-1}$|⁠. The μ approach is perhaps more direct, but is essentially indicating, as above, that more gas is needed with increasing redshift in order to form stars at the necessary levels, recycle the ISM, and achieve the observed reduced metallicities.

The increase in gas content with redshift roughly ∝(1 + z)^0.2 given by equations (8) and (10) is similar to (the inverse of) that found observationally for molecular depletion times τ_depl of main-sequence galaxies by Genzel et al. (2015). Since, as discussed above, our models are unable to separate gas content from τ_depl and ε_ff, this is an encouraging consistency. Moreover, for massive galaxies with M_* ≈10¹¹ M_⊙, we would predict gas fractions seven times higher at z ∼ 2 than at z ≃ 0, roughly consistent with the CO observations by Geach et al. (2011), and only slightly lower than the increase of a factor of 10 found by Bothwell et al. (2013b) for luminous sub-millimetre galaxies.

Nevertheless, a limitation of our model is that it does not explicitly distinguish between accreted gas, presumably H i, and the gas which forms stars that, at these high SFRs, must be molecular. Indeed, the changes with redshift of relative gas content predicted by the FPZ(z) model are slightly lower than the results derived observationally for H₂; Tacconi et al. (2013) and Saintonge et al. (2013) find that, as z goes from 2 to 1, gas fractions decrease by a factor of ∼1.4, while we would predict a change of only ∼1.1–1.2 (for a galaxy of fixed stellar mass). The gas content of our model comprises not only the molecular component but also H i; although it is not yet possible to observe H i at high z, it is likely that H i content shows a smaller redshift variation than H₂ (e.g. Lagos et al. 2011; Popping et al. 2012; Lagos et al. 2016). Thus, total gas fractions are expected to increase less rapidly with redshift than H₂ alone. More detailed comparisons with observations will require observations of atomic gas to cosmological redshifts which should be possible with the Square Kilometre Array (SKA).

5.2 Redshift variation of the stellar mass scaling

It is well established that at a given redshift the gas fraction in massive galaxies is lower than in less massive ones (e.g. Saintonge et al. 2011; Huang et al. 2012; Popping et al. 2012; Boselli et al. 2014; Bothwell et al. 2014; Popping et al. 2015). Moreover, as discussed above, it is clear that gas fraction increases with redshift. The FPZ(z) model with the ε_* and μ approaches also predicts that the amount of increase in gas content with redshift varies with stellar mass. In particular, the M_* scaling of the FPZ(z) models (see equations 8 and 10) suggests that the gas content of high-mass galaxies relative to less massive ones increases with increasing redshift. For a galaxy with M_* = 10¹⁰ M_⊙ at z = 2, the growth in gas content is predicted to be ∼2.2 times that at z = 0; for a more massive galaxy, M_* = 10¹¹ M_⊙ at the same redshift the increase would be approximately seven times the gas content of a galaxy of the same mass at z = 0. A less massive galaxy, with M_* = 10⁹ M_⊙, would be expected to have roughly the same gas content at z = 2 as at z = 0. Such a result would be consistent with ‘downsizing’ in which massive galaxies evolve more rapidly than lower-mass ones, thus consuming their gas at earlier epochs (e.g. Cowie et al. 1996; Bundy et al. 2006; De Lucia et al. 2006; Thomas et al. 2010). Specifically, Thomas et al. (2010) found that SF activity in galaxies with M_* ∼10¹¹ − 4 × 10¹¹ M_⊙ peaked at increasingly higher redshifts, z ∼ 1.5–2, well within the redshift and mass ranges probed by the MEGA and simulated galaxies.

Such a trend of gas fraction with stellar mass and redshift is also in accord with some observational estimates of the mass variations of gas fraction with redshift, either through indirect determinations of gas content by inverting gas-SFR scaling relations (e.g. Popping et al. 2012, 2015), or with gas masses inferred from measured dust masses up to z ∼ 2.5 after correcting for metallicity (Santini et al. 2014). Santini et al. (2014) find that from z ∼ 2.5 to z ∼ 1, the gas fraction in massive galaxies decreases more sharply than in lower-mass ones; less massive galaxies show a shallower decrease in the gas content.

However, observations of H₂ alone are in disagreement with this. Morokuma-Matsui & Baba (2015) estimated redshift variations of H₂ fraction through a compilation of CO observations and find that more massive galaxies show less evolution in H₂ fraction than lower-mass ones, in direct contrast to our results. The same contradiction emerges in the study by Dessauges-Zavadsky et al. (2015) who studied five lensed galaxies and compared them with CO observations from the literature. Genel et al. (2014) Illustris simulations also find that less massive galaxies at higher redshift have a larger change in gas fraction with redshift. Observationally, our result for galaxies with masses ≲ 10¹⁰ M_⊙ is very difficult to verify given that virtually all observational studies so far of gas content at high redshift are limited to M_* ≳ 10¹⁰ M_⊙ (although see Saintonge et al. 2013; Dessauges-Zavadsky et al. 2015, who studied lensed galaxies). Nevertheless, the results of the FPZ(z) models suggest that massive galaxies at z ∼ 3 relative to z ≃ 0 must have a higher gas content in order to achieve the relatively lower metallicities. Incorporating a multi-phase approach as in Magrini et al. (2012) might help interpret the different gas components in our models.

One limitation of the FPZ model was discussed above, namely the inability to distinguish between H₂, H i, and total gas content. At least one other possible limitation to the model is the lack of mass scaling of ε_*, or equivalently the SF time-scale. Such an omission can be justified for the SDSS10 galaxies because, as discussed in Section 3, in typical spiral discs both molecular depletion times τ_depl and H i depletion time are roughly constant, 2–3 Gyr (Catinella et al. 2010; Schiminovich et al. 2010; Bigiel et al. 2011). Nevertheless, this could be a problem for the MEGA data set because of its highly star-forming nature, and large variation in M_*. It could also be troublesome for the high-z main-sequence simulated galaxies because of their relatively higher SFRs compared to local ones. It is possible that the M_* scaling we find for the ε_* and μ adaptations of the FPZ model to z > 0 is a consequence of neglecting such a treatment at z = 0.

Our treatment of trends of FPZ parameters with z and M_* as simple power laws is also highly over-simplified. As found by Lagos et al. (2011, 2014) and Morokuma-Matsui & Baba (2015), the redshift and M_* dependence of gas-mass fractions is more complicated than this, as there are inflections and slight curvatures in the behaviour of both quantities. Thus, our models are almost certainly not exact, but rather a simplified scaling aimed at a general representation of the trends necessary to explain the coevolution of SFR and metallicity.

Recently, it has been proposed that scaling laws based on gas fraction, SFR, and M_* are more fundamental than those based on metallicity (e.g. Bothwell et al. 2013a; Santini et al. 2014; Bothwell et al. 2016; Lagos et al. 2016). The physical basis of such a scaling is clear, since gas provides the fuel for star formation and the reservoir of gas-phase metals. However, the role of H i versus H₂ is still under debate, and the possibility of verifying through observations the amount of total gas in H i is currently not possible, although SKA will be able to help shed light on this problem.

6 SUMMARY

In a companion paper, we constructed a new MEGA data set for studying the coevolution of metallicity and SFR; the data are compiled from 19 different samples up to z ≃ 3.7, spanning a factor of ∼10⁵ in M_*, ∼10⁶ in SFR, and almost two orders of magnitude in O/H. As a comparison sample, we include the SDSS10 galaxies studied by Mannucci et al. (2010).

• Here we update the model of Dayal et al. (2013) for the z ≃ 0 behaviour of M_*, SFR, and O/H and apply it to the SDSS10 and MEGA data sets; the new model relies on the stellar yields from Vincenzo et al. (2016) and Nomoto et al. (2013) assuming a Chabrier (2003) IMF. With only five free parameters, the FPZ model at z ≃ 0 is able to reproduce the observed oxygen abundance of the ∼80 000 SDSS10 galaxies to within 0.05−0.06 dex, and of the ∼250 z ≃ 0 MEGA galaxies to within 0.24–0.26 dex.

• We have extended the FPZ(z ≃ 0) model to higher redshift by exploring three possibilities: (1) a redshift (and mass) scaling of SF time-scale ε_* (Section 4.1); (2) a redshift scaling of accretion and galactic wind mass loading (Section 4.2); and (3) a redshift (and mass) scaling of gas fraction (Section 4.3). These approaches give similar results and are able to reproduce the observed metallicity trends in the MEGA sample to within ∼ 0.03–0.05 dex.

• The extension of the FPZ(z) model allows us to quantify how gas accretion and outflow depend on redshift. Although the specific mass loading of outflows does not change measurably during the evolution, the accretion rate and gas content of galaxies increase significantly with redshift. These two effects explain, either separately or possibly in tandem, the observed lower metal abundance of high-z galaxies.

Acknowledgments

We thank Filippo Fraternali for interesting discussion and insightful comments. This work has benefited from the DAVID network (http://wiki.arcetri.astro.it/DAVID/WebHome) for fostering a fruitful collaborative environment. PD gladly acknowledges funding from the EU COFUND Rosalind Franklin programme, and LH from a national research grant, PRIN-INAF/2012.

1

The stellar yields from Romano et al. (2010) are similar, with (R, y) = (0.45, 0.06). Our values for (R, y) are roughly the mean values of the yields over the range of metallicities given by Vincenzo et al. (2016) since our model does not contemplate time variation of the yields.

2

This corresponds to an oxygen abundance of 4.9 × 10⁻⁴ while Dayal et al. (2013) used 7.9 × 10⁻⁴.

3

To explore the molecular gas fraction necessary to form stars would be an important addition to this discussion, but is beyond the scope of this work.