© ESO, 2016

1. Introduction

For more than 30 years, cosmological models based on the Λ cold dark matter(CDM) paradigm have explained the origin and evolution of structures in the Universe remarkably well (Blumenthal et al. 1984; Fu et al. 2008). However, the description of the baryonic processes acting at smaller scale (galaxies) remains problematic.

Gas flows and star formation are supposed to be the main actors of galaxy evolution (e.g. Dekel et al. 2009a). These mechanisms leave an imprint in the amount of heavy elements that are produced by stars and progressively enrich the interstellar medium (ISM). A better understanding of both the gas phase and stellar metallicities is essential to better constrain the galaxy formation scenario.

In the early 2000s, the analysis of hydrodynamic simulations allowed the formulation of a bimodal accretion scenario (e.g. Birnboim & Dekel 2003; Kereš et al. 2005; Ocvirk et al. 2008; Dekel et al. 2009a,b; Khochfar & Silk 2009; Brooks et al. 2009; Faucher-Giguère et al. 2011). According to this bimodal scenario, the galaxies can grow through a cold and a hot mode. The smaller galaxies grow through a very efficient cold mode. At high redshift (z> 2), this mode is based on filamentary streams that feed the galaxy directly. In the hot accretion mode, the cosmological gas is shock heated in the halo, and then this hot gas has to cool before feeding the galaxy. The hot mode is predicted to complete and progressively replace the cold mode in more massive structures (Faucher-Giguère et al. 2011). However, even when a hot isotropic atmosphere is present around a galaxy, some residual cold accretion can exist, in particular for massive halos at high redshifts (e.g. Thom et al. 2011; Ribaudo et al. 2011). The recent works of Nelson et al. 2013, based on a new generation of hydrodynamic codes, call into question this bimodal accretion scenario. According to this work, all the accreted gas could follow the hot mode. In this context, the existence of cold streams is questioned.

In parallel to hydrodynamic simulations, semi-analytic models (SAMs) are good tools for studying the impact of these different prescriptions. Since the pioneering work of White & Frenk (1991), the assumptions of i) a virial shock and ii) an isotropic hot atmosphere have been used in SAM (e.g. Cole et al. 2000; Hatton et al. 2003; Somerville et al. 2008; Guo et al. 2011; Henriques et al. 2013). A rapid and a slow cooling accretion mode can be identified in semi-analytical prescriptions, even if all the accreted gas is assumed to be shock heated (Benson & Bower 2011). These two cooling regimes are linked to the dynamical timescale of dark matter structures. In the case of low-mass structures the cooling rate is very efficient and the gas is rapidly accreted onto the galaxy. With this kind of model, the gas accretion cannot be followed in detail. In particular it is impossible to explicitly follow metals in the two accretion modes. Observations (e.g. Thom et al. 2011; Ribaudo et al. 2011; Crighton et al. 2013) and simulations (e.g. Fumagalli et al. 2011) indicate that the two accretion modes are associated with different metallicity contents. The cold mode corresponds to metal-poor gas, log (Z/Z_⊙) ≤ − 2.0. The hot gas is found with a higher metallicity, log (Z/Z_⊙) ≥ − 1.0. The metals, observed in the hot gas surrounding galaxies are initially formed by stars in the galaxy. The enriched gas is then ejected by a set of feedback mechanisms in the form of supernovae (SN) and active galaxy nuclei (AGN).

These feedback processes are also essential to regulate the star formation activity in galaxies but they are still misunderstood (e.g. Bouché et al. 2012). Indeed, standard scenarios of regulation of star formation lead to strong excesses of low-mass structures, mainly at high redshift (z> 3). This standard scenarios based only on galactic fountains are gradually abandoned in favour of models using long re-accretion cycles (e.g. Guo et al. 2011; Henriques et al. 2013) or an intermediate gas reservoir of non-star-forming gas that strongly regulates the star formation activity Cousin et al. (2015a).

In this work, we compare two accretion scenarios and two different star formation recipes. The first accretion scenario is described in detail in Cousin et al. (2015b). It assumes a bimodal accretion as initially proposed by Khochfar & Silk (2009) or Benson & Bower (2011). This model is based on a cold and a hot reservoir. Both are fed by the metal-free cosmological accretion but only the hot reservoir receives the galactic metal-rich ejecta in addition. As the metallicity of the wind phase depends on the galaxy metal enrichment process, the metallicity of the hot reservoir evolves with time. In the second accretion scenario, we assume a unique gas reservoir that receives both the metal-free cosmological accretion and metal-rich gas ejected by the galactic feedback processes. In addition to these two accretion scenarios, we compare two different star formation recipes. First, we assume a standard star formation process where the gas is progressively converted into stars by assuming a given efficiency. Second, we use the recipe proposed by Cousin et al. (2015a): the freshly accreted gas is assumed to be non-star forming. It is progressively converted into star-forming gas and then into stars.

Our objectives are to characterise the impacts of these scenarios onto the galaxy metallicity signature in both stars and gas measured at z ≃ 0. Our work is based on a new chemodynamical package that we present here. Similar to its predecessors, Nagashima et al. (2005), Arrigoni et al. (2010), Yates et al. (2013), or more recently De Lucia et al. (2014), our chemodynamical model is based on stellar evolution models. We assume a Chabrier (2003) IMF¹. The new library is built with a large set of initial stellar metallicties, including metal-free stellar populations. We also introduce a time-deleted SN feedback prescription, based on SNII and SNIa event rates.

Our analysis is performed on the whole sample of our simulated galaxies and on a sub-sample of star-forming (SF) galaxies. Observational measurements show that star formation rates (SFRs) and stellar masses are strongly correlated for SF galaxies (e.g. Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007). This tight correlation is called main sequence. Our sub-sample of SF galaxies is therefore defined in the context of this main sequence.

The paper is organised as follows. In Sect. 2 we describe the two accretion scenarios and give an overview of the other components of the model. In Sect. 3 we describe the new chemodynamical model used to follow the metal-enrichment process. In Sect. 4 we present the four different models and we describe the set of rules used to define our star-forming galaxy sample. In Sect. 5 we present the metallicity signatures in both stars and gas phase measured at z ≃ 0 for our different scenarios. This first analysis allows us to define a best model among the four different prescriptions used. With this best model, we explore the link between the average metallicity of the gas and star formation activity (e.g. Ellison et al. 2008; Mannucci et al. 2010; Lara-López et al. 2010). In Sect. 6 we focus on a sub-sample of galaxies: the Milky Way (MW) like galaxies. We analyse the impact of a set of selection criteria onto several properties such as the stellar and gas metallicity, the gas fraction, and the time elapsed since the last major merger. Then we detail their star formation and their metal enrichment histories. Finally, we replace the properties of the progenitors of MW-like galaxies in the context of the main sequence evolution. The conclusion of this work is presented in Sect. 7.

2. The eGalICS model

We use the new eGalICS model described in detail in Cousin et al. (2015b). This SMA is applied on dark matter merger trees extracted from a pure N-body simulation. This simulation is based on a WMAP-3yr cosmology (Ω_m = 0.24, Ω_Λ = 0.76, f_b = 0.16, h = 0.73) and describes a volume of [100 /h]³ Mpc. In this volume, 1024³ particles evolve. Each particle has an elementary mass of m_p = 8.536 × 10⁷ M_⊙. We use the HaloMaker code described in Tweed et al. (2009) to identify haloes and their respective sub-structures (satellites) and build the corresponding merger trees. In the merger tree evolution, we only consider halos containing, at least 20 dark matter particles. This limit leads to a minimal dark matter halo mass of 1.707 × 10⁹ M_⊙.

2.1. Dark matter and baryon accretion

Dark matter haloes are fed following a smooth accretion mode. The dark matter accretion rate Ṁ_dm associated with a halo is built between two time-steps with the particles that are newly detected in the halo ℋ, and that have never been identified in an other halo before. The baryonic mass is then added progressively following this smooth, dark matter accretion, (1)In this relation, is the effective baryonic fraction that depends on the photoionisation model².

From this global baryonic accretion (Ṁ_b), we apply two different models. The first model is based on a single hot isotropic accretion. The second model assumes a bimodal accretion. The baryonic cosmological accretion is therefore divided into two phases. At a given redshift and for a given halo ℋ, the fraction of shock-heated gas (f_sh(M_vir,z)) is computed following the Lu et al. (2011a) prescription (their Eqs. (24) and (25)) (2)We use two distinct reservoirs (cold and hot). The cold reservoir is fed by metal-free cosmological gas when the hot reservoir receives both metal-free cosmological gas and metal-enriched gas ejected by the galaxy. We assume that the gas ejected by the galaxy has no impact on the cold mode. While the cold reservoir feeds the galaxy directly, the hot gas is assumed to reach the virial temperature and has to support the radiative cooling process before feeding the galaxy.

The cooling rate is computed using the classical model proposed by White & Frenk (1991). The condensed mass enclosed in the cooling radius r_cool is estimated assuming

• a hydrostatic equilibrium gas profile ρ_g(r);

• a mean constant temperature ; and

• a temperature and metal dependent cooling function (Sutherland & Dopita 1993).

The cooling rate of the hot atmosphere is then computed as (3)where the dynamical time is adapted to the cooling radius³ of the hot gas structure. The efficiency of the cooling process is governed by the free parameter ε_cool (see Table 1).

2.2. Star formation in disks

The accreted gas feeds the galaxy evolving in the centre of the dark matter halo. We assume that this cold gas initially forms a thin exponential disk. Gas acquires angular momentum during the mass transfer (Peebles 1969). After its formation, the disk is supported by this angular momentum (e.g. Blumenthal et al. 1986; Mo et al. 1998).

We apply two different models for star formation. A standard model described in Cousin et al. (2015b) and the recipes defined in Cousin et al. (2015a). In the standard prescription, the star formation is based on the empirical Schmidt-Kennicutt law (Kennicutt 1998). In this context we use a critical surface density threshold (Σ_crit = 10 M_⊙ pc^-2, Schaye 2004) and assume that only gas lying at surface density larger than this limit is available to form new stars⁴. In the second model, the accreted gas is assumed to be initially non-star forming. It is then progressively converted into star-forming gas and concentrates in the centre of the disk. We call r_s the radius that encloses the star-forming gas. The SFR is given by (4)where ε_⋆ is a free efficiency parameter (see Table 1), M_g^⋆ is the mass of star-forming gas. The star-forming timescale t_dyn is then given by(5)where V_c is the circular velocity of the disk computed at the radius r_s by taking into the dark matter, the bulge, and the disk potential wells into account.

2.3. Supernovae feedback

We compute the ejected mass rate due to SN using kinetic energy conservation (Dekel & Silk 1986) (6)We use E_sn = 10⁴⁴ joules. The wind velocity is fixed and set to V_wind = 250 km s^-1, and ε_ej is a free efficiency parameter. Unlike the previous version of the model, the ejecta rate is not directly linked to the SFR but to more realistic SNII and SNIa event rates η_sn, thus integrating all the stellar populations formed previously. The explicit computation of η_sn is given in Sect. A.3. The realistic SNII and SNIa event rates η_sn, used in this new version, has modified the average mass-loading factor produced by SN feedback. We increased the efficiency of the process to keep the same total SN feedback impact (ε_ej = 0.6 >ε_snf_kin,sn = 0.3).

2.4. Mergers

We define the merger type (major or minor) of two galaxies 1 and 2 using the following mass ratio, (7)where M_{1 / 2} is the total mass (galaxy and dark matter halo) inside the galaxy half-mass radius. We consider the merger event as minor when η_merger< 1 / 3.

In the case of major merger (η_merger> 1 / 3), the galaxy structure and dynamics are completely modified. The stellar population of the remnant galaxy is a spheroid modelled by a Hernquist (1990) mass distribution. The gas components of the two progenitors are added and a new disk is formed.

During a minor merger, we add the stellar population of the bulges and disks separately to form the remnant bulge/disk. All cold gas reservoirs are added to the remnant disk.

To take the strong modifications impacting galaxies during a merger into account, Cousin et al. (2015a) defined a boost factor, MAX[1,ε_boost(Δt)], which increases the efficiencies of the star formation process. The boost factor is defined as (8)This boost factor simulates a gas compression (decreasing with time) as a function of the gas content of the two progenitors, where η_gas is the gas fraction in the post-merger structure, t is the time elapsed since the last merger event, τ_merger = 0.05 Gyr is a characteristic merger timescale, and η_merger is given in Eq. (7).

3. The chemodynamical model

The impact of stars on the ISM-enrichment process is a function of their stellar mass (m_⋆) and of their metallicity (Z_⋆). We take stellar masses into account between m_⋆,min = 0.1 M_⊙ and m_⋆,max = 100 M_⊙. This scale has been chosen to match the stellar mass range of the BC03 spectrum library.

This stellar mass range can be divided into three bins (Romano et al. 2005): i) the very low-mass stars (VLMS) [0.1 − 1 [ M_⊙; ii) the low- and intermediate-mass stars (LIMS) [1 − 9]M_⊙; and iii) the massive stars (MS) ] 9 − 100]M_⊙. For each mass bin, the metal production and the ISM re-injection processes are different.

The new chemodynamical library is built by following six of the main ISM elements ¹H, ⁴He, ¹²C, ¹⁴N, ¹⁶O, and ⁵⁶Fe. The production and re-injection processes are associated with stars with initial mass between m_⋆,min and m_⋆,max, and for seven metallicity⁵ bins listed in Table 2. These bins are able to cover gas composition from metal-free to super-solar metal fraction. We built this metallicity scale by merging different sets of stellar evolution models: Portinari et al. (1998), Campbell & Lattanzio (2008), Karakas (2010) and Heger & Woosley (2010), and Gil-Pons et al. (2013) (see Table 2). Some more recent library of massive stars yields are currently available (e.g. Geneva group). These works have shown that C, N, O, for example are affected by mass loss and rotation. The current work is dedicated to the chemical evolution of galaxies, but we have also included in parallel the prediction of stellar population spectra via the BC03 spectra library. To stay in agreement with the BC03 library, we based our yield library on the Padova model and we have therefore used the same metallicity scale.

3.1. Very low-mass stars (VLMS)

The VLMS contains stars with mass between m_⋆ = 0.1 M_⊙ and m_⋆ = 1.0 M_⊙. We assume that VLMS are eternal. They simply remove their mass from the gas evolution cycle.

3.2. Low- and intermediate-mass stars (LIMS)

The LIMS are stars with a initial mass between m_⋆ = 1 M_⊙ and m_⋆ = 9 M_⊙. They end their life as white dwarfs (WD) after passing through a phase of double-shell burning, known as the thermally pulsing asymptotic giant branch (TP-AGB) phase. During this phase He, C, and N are convected to the outermost layers and re-injected into the ISM. We assume that all the mass ejected by a LIMS is injected in the ISM at the end of the stellar life⁶ of a star of initial mass m_⋆. Ejecta processes coming from LIMS are taken into account following Karakas & Lattanzio (2007) and Karakas (2010). All their models were evolved from the zero-age main sequence to the tip of the TP-AGB. The details of this procedure is described in Karakas et al. (2009) and references therein. In these works, ejected masses are given for a range of initial stellar masses, m_⋆ = [1 − 6]M_⊙, and a set of metallicities, Z_⋆ = [0.0001, 0.004, 0.008, 0.02]. These metallicity bins correspond to the , , , and bins of our chemodynamical library. To complete the metallicity set, we interpolated(extrapolated) the total ejected mass produced by LIMS on a metallicity log scale for metallicity bin (). To complete the initial stellar mass range given by Karakas (2010), [1 − 6]M_⊙, we take ejected mass coming from stars with intermediate masses m_⋆ = 7 M_⊙ and m_⋆ = 9 M_⊙ into account following Portinari et al. (1998).

3.3. Massive stars (MS)

The MS correspond to m_⋆> 9 M_⊙. They represent only a small fraction of the stellar mass, i.e. 21% for a Chabrier (2003) IMF (e.g. Romano et al. 2005, 2010). Nevertheless, MS are very important contributors to the ISM enrichment process. Their nucleosynthesis produces heavy elements that are directly re-injected into the ISM through stellar winds during their entire life. MS end their life in a type II supernovae episode (SNII) during which a large amount of energy and the majority of the newly synthesised metals are injected into the ISM. Ejecta from MS are taken into account following Portinari et al. (1998). This choice is based on the large available initial stellar mass range they provide, i.e. m_⋆ = [12, 120]M_⊙. Portinari et al. (1998) also take the mass loss into account through winds prior to the SNII episode. Stellar ejecta have been computed for the following metalicities: Z_⋆ = [0.0004, 0.004, 0.008, 0.02, 0.05]. To complete our metallicity scale, ejecta associated with MS with Z = 0.0001 were extrapolated on a metallicity log scale. As for LIMS we consider that the elements released during the SNII phase are instantaneously injected into the ISM at a time t corresponding to the lifetime of the star t = τ_m. Before the SNII episode, we also assume a continuous wind ejecta process acting during a time lapse equal to the lifetime of the star Δt = τ_m.

3.4. Total ejected masses, the overall mass range

Figure 1 shows the evolution of the IMF-weighted total ejected mass for six main ISM elements (¹H, ⁴He, ¹²C, ¹⁴N, ¹⁶O, and ⁵⁶Fe). The metal ejected mass and total ejected mass are also plotted.

In the LIMS range, [1, 9]M_⊙, a good complementarity is found between ejected masses coming from Karakas (2010)[1 − 6]M_⊙ and Portinari et al. (1998)[7,9]M_⊙.

In the different panels, dedicated to different metallicities, we note a gradual increase of the metal ejected mass with the metallicity. For example, solar metallicity (Z_⋆ = 0.02) LIMS eject, on average, ten times more oxygen than LIMS with metallicity Z_⋆ = 0.0004.

3.5. Event rates

To convert the mass of gas ejected by each star (according to its initial mass) in a continuous ejection rate over time, we have to build the event rates associated with each specific ejection process.

3.5.1. TP-AGB and Type II supernovae events

For LIMS we assume that the mass is ejected during the TP-AGB phase. For MS the majority of the mass is ejected during the type II supernovae episode. The remaining part is continuously ejected by the stellar winds during the whole life of the MS. The TP-AGB and the Type II SN event rates are therefore given by the death rate of their associated stars. This death rate is given by the derivative of the inverse function of the stellar lifetime τ_m, (9)where ϕ(m) is the IMF. For stars with masses larger than 1 M_⊙ the stellar lifetime is a decreasing function of the stellar mass. With this formulation, the injection process associated with a single stellar population (SSP) begins with the contribution of the most massive stars (m_⋆ = 100 M_⊙). It is then progressively fed by less and less massive stars (m_⋆ → 1 M_⊙).

3.5.2. Type Ia supernovae contribution

In parallel to the TP-AGB phases of LIMS and SNII explosion of MS, the ISM is also enriched by type Ia supernovae explosions coming from binary systems. The observed properties of SN Ia suggest that these systems may be produced by the thermonuclear explosion of accreting WD, even if the exact binary evolution has not yet been identified e.g. Greggio & Renzini 1983; Matteucci & Greggio 1986; Nomoto et al. 1995; Wheeler et al. 1995; Matteucci & Recchi 2001).

We compute the SNIa contribution according to the model proposed by Greggio & Renzini (1983). The mass ejected for each of the six elements that we followed are given by the model of Iwamoto et al. (1999) (W7).

The SNIa event rate following a burst of star formation and leading to a stellar mass of 1 M_⊙ (SSP) is given by (10)where the integration limits M_B,inf and M_B,sup are given by (11)and (12)In these previous expressions, M_B = M₁ + M₂ is the total mass of the binary system. We assume that this mass is enclosed between M_BM = 16 M_⊙ and M_Bm = 3 M_⊙, and M₂ is the mass of the companion. At a given time t after the formation, M₂ is the mass of the LIMS that died at t and is therefore given by the inverse function of the stellar lifetime m(τ_m = t). As the stellar lifetime is a function of the stellar metallicity, the SNIa event rate is impacted by the stellar metallicity. In Eq. (10), μ = M₂(t) /M_B is the companion mass ratio and f(μ) = 24μ² is the distribution function of binary systems. Finally, is the fraction of binary systems with the right characteristics to become SNIa.

3.6. Ejecta rates for a SSP

We introduce here the expression of the time dependent ejecta rate, ė(t), generated by a SSP. This formulation is the elementary step used in the computation of the evolution of a more complex stellar population.

The time dependent ejecta rate is given by the sum of three different terms. Each of these terms corresponds to a specific enrichment process, i.e. In this expression, the first term corresponds to the mass of element elt ejected by stars of initial mass m_⋆ and initial metallicity Z_⋆ at the end of their life. Stars of smaller and lower mass m_⋆ die progressively, each ejecting gas into the ISM. m_ej(m_⋆,Z_⋆,elt) represents the contribution of LIMS during their TP-AGB phase or the contribution of the MS during the SNII explosion. The second term represents the contribution of SNIa. It is based on m_ej(snIa,elt), which is the mass of the element elt ejected by a type Ia supernovae. The third term of Eq. (15) represents the MS wind contribution. We assume that the mass ejected through these winds are continuously injected in the ISM during the stellar lifetime, Δt = τ_m. m_ej,w(m_⋆,Z_⋆,elt) is the mass of the element elt ejected through stellar winds by a MS of initial mass m_⋆ and initial metallicity Z_⋆. η_⋆ is the instantaneous stellar death rate and η_snIa the instantaneous SN rate.

Figure 2 shows the evolution of the ejecta rates associated with a SSP of 1 M_⊙ for different metallicity bins. These rates are given for the main ISM elements. The initial ejecta rate corresponds to the contribution of winds from massive stars. This plateau ends at the death of the most massive star (m_⋆ = 100 M_⊙).

3.7. Metal-free primordial interstellar medium

During the first steps of a galaxy formation process, the gas accreted onto the disk is metal-free. To take this feature into account, we have to build a metal-free () bin in our chemodynamical model. We use stellar evolution models from Campbell & Lattanzio (2008) and Gil-Pons et al. (2013) for LIMS and from Heger & Woosley (2010) for MS. The instantaneous ejecta rate is given by Eq. (15) but without the MS wind term; the Heger & Woosley (2010) model does not take an explicit wind contribution into account. For metal-free MS, all the ejected mass is injected into the ISM during the SNII phase.

The upper and lower panels of Fig. 3 show the IMF-weighted total ejected mass and the time evolution of the instantaneous ejecta rates computed for a metal-free SSP of 1 M_⊙. Without the MS wind contribution, there is no plateau at the beginning of the time evolution of the ejecta rates. Gas and metals are injected only after the death of the most massive stars. For SNIa, the ejected mass is still deduced from the W7 model of Iwamoto et al. (1999).

3.8. From single stellar populations to realistic star formation histories

During its evolution a galaxy forms stars in a continuous process. All these populations have to be followed in all the galaxies. The formation and evolution of the stellar population are taken into account using a set of history tables described in Appendix A.

4. Models and sample selection criteria

4.1. The models

We generated four different models, noted m₁ to m₄. We based m₁ and m₃ on the isotropic hot accretion paradigm, and m₂ and m₄ on bimodal accretion in which small structures are essentially fed by cold metal-free gas. We based m₁ and m₂ on the standard star formation process (Cousin et al. 2015b) while, in m₃ and m₄, the new stars are formed following the recipe of Cousin et al. (2015a, see Sect. 2.2). The different configurations are listed in Table 3.

The values of the star formation efficiency, ε_⋆, and the ejecta efficiency, ε_sn, are the same in all models (Table 1).

4.2. Sample selection: star-forming galaxy sample

Star-forming galaxies (SF) are the subject of many observational studies. These objects are key elements of the stellar mass assembly. The relation found between their stellar mass and their SFR (main sequence) argues for a secular evolution. This relation is the main topic of a large number of studies (Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007).

The tight correlation between the stellar mass and the SFR, leads to the measure of a specific SFR (sSFR = SFR/M_⋆) that slowly evolves with the stellar mass. We define our sub-sample of SF galaxies using a sSFR range that evolves with the stellar mass. We use the average relation measured by Peng et al. (2010) in the SDSS-zCOSMOS survey (their Fig. 1), i.e. (16)From this average relation, at a given stellar mass, the sSFR range is defined by applying a scatter of 0.3 dex corresponding to the standard dispersion measured around the main sequence by, for example Salim et al. (2007), Elbaz et al. (2007), or Peng et al. (2010).

We add a set of rules to the sSFR criterion to define our sub-sample of SF galaxies. To be included in the sample, a galaxy must:

• be hosted by a dark matter halo with M_vir> 10¹⁰ M_⊙; and

• not have suffered a major merger during the last 0.05 Gyr.

The first criterion reduces the effects of resolution in mass. We recall that the minimum halo mass is .

During a (major-)merger, a SF-galaxies leaves the main sequence for a short time. During this period the SFR of a merger galaxy is larger than that of a main sequence galaxy of similar stellar mass. The last criterion allows us to remove those galaxies that have recently suffered a major merger.

Figure 4 shows the M_⋆-SFR relation extracted at z = 0.01 and z ≃ 0.09 from our model m₃ for the whole sample and for our SF galaxy sample. These results are compared to observational data of Peng et al. (2010) based on the SDSS-zCOSMOS survey (0.02 <z_obs< 0.09). As illustrate in Fig. 4 the predictions of our model extracted at the two limits of the Peng et al. (2010) sample are very close and almost indistinguishable.

Comparing to the observed relation, we note that our SFRs are lower in low-mass objects (M_⋆< 10¹⁰ M_⊙), whereas the star formation activity is higher in massive objects (M_⋆> 10¹¹ M_⊙). The gap is strongly reduced for our sub-sample of SF galaxies. Concerning the high-mass range (M_⋆> 10¹¹ M_⊙), it is possible to invoke the activity of a super-massive black hole (SMBH) to reduce gas accretion onto the galaxy and reduce star formation activity (e.g. Croton et al. 2006; Bower et al. 2006; Somerville et al. 2008). As explained in Cousin et al. (2015b), the impact onto the cooling process of the SMBHs in our model is weaker than in other models and does not sufficiently reduce star formation activity in massive galaxies. In the low-mass regime (M_⋆< 10¹⁰ M_⊙), the galaxies have star formation activity that is too low. The difference between the mean trend of our whole sample and the observed trend is close to 0.4 dex. These low SFRs are due to a very efficient regulation of the star formation process (in this case, the non-sfg reservoir). This regulation must be strong to reproduce the faint end of the stellar mass function. We point out here the difficulty of simultaneously matching the stellar mass contents of the galaxies and the rhythm of this stellar mass assembly. This discrepancy is also highlighted by Henriques et al. (2015). Indeed, even if the last version of the Munich model leads to an excellent agreement on the stellar mass function, the cosmic SFR density is lower than the observed data points. Their main sequences also seem to show a SFR that is too low for the lowest galaxies.

5. Galaxy metal enrichment

In this section, we compare the behaviour of metallicity in our model to local (z ≃ 0.09) data on the stellar mass to stellar metallicity relation Sect. 5.1, on the stellar mass to gas-phase metallicity, Sect. 5.2, and on the relation with the SFR Sect. (5.3).

5.1. Stellar metallicity

We present the M_⋆ − Z_⋆ relation for our four different models in Fig. 5. The luminosity-weighted stellar metallicity is Z_⋆ and is given in terms of mass fraction relative to the solar metallicity (Z_⊙ = 0.02). Our model predictions are compared to the observational relation measured by Peng et al. (2015). This measure is based on more than 3900 nearby star-forming galaxy spectra extracted from the SDSS-DR7 survey.

Our predictions are also compared to the Gallazzi et al. (2005) measurements. Based on the previous SDSS-DR2 survey, these investigators measured both stellar mass and light-weighted stellar metallicity for more than 53 000 nearby galaxies (0.005 <z< 0.22). This work is still a reference and is largely used for comparisons with models (e.g. Yates et al. 2013; Lu et al. 2015).

We also compare our M_⋆ − Z_⋆ relations with a set of predictions from different SAMs (Croton et al. 2006; Somerville et al. 2012; Lu et al. 2011b). The mean trend predicted by these SAMs was extracted from the recent paper of Lu et al. (2014). In these three SAMs the metal enrichment process is computed by assuming the instantaneous recycling approximation and fixed values for metal yields.

5.1.1. General behaviours

The general behaviour our M_⋆ − Z_⋆ relations are in general agreement with those observed by Gallazzi et al. (2005). For stellar masses above 10⁹ M_⊙, the M_⋆ − Z_⋆ relations exhibit very similar trends for the four different models. In this range, the average M_⋆ − Z_⋆ relation, computed from both the whole sample and the SF galaxy sample, are indistinguishable. These strong similarities between models are mainly due to the computation of a luminosity-weighted stellar metallicity. This quantity is mostly sensitive to young stellar populations: the youngest stellar populations, formed at z ≃ 0., have similar metallicity properties, even if they evolved previously in different conditions.

In our four models, the majority of our modelled galaxies in the range 10⁹<M_⋆< 10^10.5 M_⊙ are located on a plateau around log (Z_⋆) ≃ − 0.5. This average value is consistent with those measured by Gallazzi et al. (2005) around 10⁹ M_⊙. After this plateau, our models show a strong and prompt increase of the metallicity in the high-mass domain as measured by Gallazzi et al. (2005). After this strong increase in metallicity, the M_⋆ − Z_⋆ relation again becomes constant with a level close to solar metallicity. This average level is slightly smaller (≃0.1 dex) than that measured by Gallazzi et al. (2005) in the SDSS data.

Unlike our models, the other SAMs do not reproduce the general trend observed in the SDSS-DR2 survey by Gallazzi et al. (2005). However, their average trend is similar to that measured by Peng et al. (2015) in the SDSS-DR7 survey.

In our models, the increase starts above 10^10.5 M_⊙. This transition mass is higher than that obtained by Gallazzi et al. (2005) with a rapid increase in metallicity around 10^9.5 M_⊙. In the mass range 10^9.5<M_⋆< 10^10.5 M_⊙, our simulated galaxies stay on the plateau and exhibit a metallicity (−0.3 dex) that is too low. In this specific range of mass, the rate of the metal enrichment is slower for our modelled galaxies than for the SDSS sample. The behaviour is similar in all of our models, meaning that the causes of the discrepancy are not linked to the accretion scenario or to the star formation process. The causes of this delay might be related to the set of hypotheses used to model the link between the stars and ISM. At any time, all the stellar populations previously formed re-inject enriched gas in the ISM. This ejected gas is instantaneously mixed with the pre-existing gas in the ISM. The gas ejected is therefore diluted and all the ISM has the same metallicity. In real galaxy discs, the stars are formed and evolve in regions of limited size. The gas ejected by the stellar population is not diluted with the total gas of the galaxy and remains distributed around the star-forming region. Some stars are formed in regions with a higher metallicity than average. In our models, the metallicity signature of the gas ejected by the SN or AGN feedback process is similar to that of the ISM. We can presume that it is more difficult to eject metals, condensed in dust, than hydrogen. Even if this second point is essentially linked to the average gas metallicity, it could also impact the stellar metallicity: a higher metallicity of the ISM leads obviously to a higher stellar metallicity.

5.1.2. Impact of model prescriptions

Even if the general trend of the M_⋆ − Z_⋆ relation is very similar in the four models, some differences can be observed.

The impact of the non-star-forming gas recipe is visible in the low-mass domain (M_⋆< 10⁹ M_⊙). In this range, the models m₃ and m₄ form a lower number of galaxies with log (Z/Z_⊙) > − 0.5. A decrease of the M_⋆ − Z_⋆ relation is observed for these models, while the average metallicity of low-mass objects is roughly constant in the two standard models (m₁ and m₂). With the non-star-forming gas reservoir, the star formation activity is strongly reduced and therefore the metal enrichment cycle is also strongly impacted.

The two different accretion scenarios lead to different M_⋆ − Z_⋆ relations below 10¹⁰ M_⊙. The M_⋆ − Z_⋆ relation predicted by models m₁ and m₃ (isotropic accretion) is above the prediction found for the models m₂ and m₄ (bimodal accretion). As expected the direct feeding with metal-free gas introduced in m₂ and m₄ reduces the average stellar metallicity (≃0.1 dex).

5.2. Gas metallicity

Oxygen is the most abundant element formed into stars. It can therefore be used as a proxy for all heavy elements and gas metallicity. We define the gas metallicity Z_g as the number of oxygen atoms to hydrogen atoms with a logarithmic scale: Z_g = 12 + log (O / H). In this formalism we adopt Z_⊙ = 8.94 (Karakas 2010).

Figure 6 shows the evolution of gas metallicity as a function of stellar mass. Our predictions are compared to various observational results:

• Tremonti et al. (2004) measured both stellar mass and gas-phase metallicity for more than 53 000 nearby (z ≃ 0.1) star-forming galaxies from the SDSS-DR2 survey (York et al. 2000). Owing to our large statistical set, this work has been the main observational set to be compared with models in more than ten years;

• Lee et al. (2006) added infrared (IR) data coming from the Spitzer Space Telescope to extend the M_⋆ − Z_g relations measured by Tremonti et al. (2004) to the low stellar mass range (down to ≃10⁶ M_⊙);

• Andrews & Martini (2013) measured the gas-phase metallicity of more than 200 000 nearby (z ≃ 0.078) star-forming galaxies from the SDSS-DR7 survey.

Tremonti et al. (2004) and Andrews & Martini (2013) used different methods to measure metallicity: i.e., the strong line and the direct method, respectively. The strong line method uses only the flux ratio of the strong lines. The metallicity of the gas is deduced using empirical or theoretical calibrations. The direct method uses the flux ratio of auroral to strong lines to obtain a direct measure of the electron temperature of the gas. This temperature enables a better determination of the gas-phase metallicity.

As discussed by Kewley & Ellison (2008) or Moustakas et al. (2010), metallicities measured using theoretical calibrations based on strong lines are systematically larger than those obtained with the direct method. A gap between the metallicity measurements performed by Tremonti et al. (2004) and Andrews & Martini (2013) is visible in Fig. 6 but found lower than the maximum systematic difference inferred by the different methods of measurements (up to ≃0.7 dex; Kewley & Ellison 2008). Given all these uncertainties on metallicity measurements, we focus on the general trends instead of absolute values to constrain our models.

The SF galaxy sample and the whole sample exhibit very similar trends for all the models with, on average, a slightly smaller metallicity for the SF galaxy sample than for the whole sample.

While the M_⋆ − Z_⋆ relation has a very similar trend for all of the models, the M_⋆ − Z_g relation shows very different behaviour for the four different models. We observe a good agreement between the non-sfg scenario (m₃, m₄) and the Andrews & Martini (2013) measurements.

5.2.1. The low-mass range: M_⋆< 10⁹M_⊙

For the stellar masses below M_⋆< 10⁹ M_⊙, the two different sets of star-forming models show distinct behaviours. The oxygen abundance is clearly over-predicted in the two classical models m₁ and m₂. The non-star-forming hypothesis used in the models m₃ and m₄ leads to a good agreement with the Lee et al. (2006) and Andrews & Martini (2013) data points.

The oxygen excess found in m₁ and m₂ is directly linked to the low-mass galaxy excess observed in standard galaxy formation models and caused by a star formation activity that is too efficient. If the star formation process is not strongly regulated, the stellar mass assembly is too fast in low-mass structures and leads to a systematic overproduction of stars. An excess of stars leads to an excess of metals. These results confirm the need for a strong regulation at the low-mass scale.

For a given star formation process ([m₁,m₂] or [m₃,m₄]), the two accretion scenarios lead to similar trends. In the non-sfg scenario (m₃,m₄), the gas-phase metallicity increases faster with an isotropic accretion than in the bimodal accretion mode. In the low-mass domain (M_⋆< 10⁸ M_⊙), the slope of the distribution is slightly larger for model m₃ than for model m₄. In addition, at 10^8.5 M_⊙ the average gas-phase metallicity predicted by m₃ is larger than that predicted by m₄ by 0.2 dex. The bimodal accretion mode allows feeding of the galaxy with metal-free cosmological gas directly. The metal-enriched gas produced by the stars is diluted with the metal-free cosmological gas. In the isotropic accretion scenario, the newly accreted metal-free gas is mixed into the hot atmosphere with the metal-enriched gas coming from ejecta. The average gas-phase metallicity is therefore higher in the isotropic accretion scenario. With the classical star formation models (m₁,m₂), the gas-phase metallicity is also slightly higher in the isotropic case.

5.2.2. The intermediate-mass range: 10⁹<M_⋆< 10^10.5M_⋆

In this mass range, the predicted gas-phase metallicities are very similar for the four models. The bimodal accretion mode leads again to a lower average metallicity for the gas phase as compared to the isotropic accretion mode. In all cases, the average metallicity of the models is consistent with that measured by Andrews & Martini (2013).

5.2.3. The high-mass regime: M_⋆> 10^10.5M_⋆

In this range of mass, the four models give also similar results. In all cases, oxygen abundances of our SF galaxies are consistent with Andrews & Martini (2013) metallicity measurements. However our relations still increase when a saturation of the metallicity is observed by Andrews & Martini (2013) or Tremonti et al. (2004). This effect is likely to be due to our AGN feedback mechanism, which is not efficient enough (see Sect. 4.4 in Cousin et al. 2015b, for more details).

5.2.4. Our best model

In the low-mass range (M_⋆< 10⁹ M_⊙), standard star formation recipes lead to an overproduction of metals. These two models cannot be selected as our best model. To identify this best model, we perform a reduced χ² test between both the Andrews & Martini (2013) and our SF galaxy sample mean trends. The reduced χ² values for each model are given in Fig. 6. As mentioned previously, the two standard models obtain high χ² values, they are therefore rejected. The model m₃, based on the non-sfg reservoir and the isotropic accretion process, is selected as the best model. Despite its better agreement in the low-mass regime, the m₄ model is rejected by the reduced χ² test⁷. The gap with Andrews & Martini (2013) measurements in the intermediate-mass range is larger than for model m₃.

In the following sections, we use the model m₃ as our best model.

5.3. SFR, M_⋆–Z_g relation

In a seminal work, Mannucci et al. (2010; see also Ellison et al. 2008; Lara-López et al. 2010) introduced a fundamental relation between stellar mass, gas metallicity, and SFR in galaxies at low and high redshift as a generalisation of the mass metallicity relation. For a given stellar mass, metallicity globally decreases with increasing SFR. This anti-correlation between star formation activity and metallicity can be because of a modification of the gas content (e.g. Bothwell et al. 2013; Lara-López et al. 2013). The main actor could be low metallicity gas inflows which dilute metals in the gas phase, fuel star formation, and increase SFR (e.g. Zahid et al. 2014).

Recently, Yates et al. (2012) and Andrews & Martini (2013) revisited these relations with the SDSS-DR7 sample, their results are reported in Fig. 7. With a sample of more than 177 000 emission line galaxies Yates et al. (2012) confirmed the dependence of metallicity on SFR at a fixed stellar mass, but with an opposite trend for low- and high-mass regimes with a transition around M_⋆ ≃ 10¹⁰ M_⊙. The relations of Andrews & Martini (2013) also shows a transition in their three first SFR-bins, which is similar to that reported by Yates et al. (2012). However, without measurements for stellar masses higher than 10^10.5 M_⊙ in Andrews & Martini (2013), it is difficult to give a definitive conclusion. Andrews & Martini (2013) use the direct method to compute gas-phase metallicity while in Yates et al. (2012), metallicities are calculated using a grid of photo-ionisation models (Charlot & Longhetti 2001). The difference in absolute levels of metallicity (Fig. 7) could be assigned to the different methods used.

From the whole galaxy sample generated with our best model m₃, we extract galaxies with − 1.5 < log (SFR) < 2.0M_⊙/ yr and 8 < log (M_⋆) < 11.5 M_⊙. These ranges of SFR and stellar masses are similar to those used in Yates et al. (2012) and Andrews & Martini (2013). Our sample is then divided in (SFR, M_⋆) bins corresponding the samples defined by Yates et al. (2012) and Andrews & Martini (2013) for the SFRs and with a width of 0.5 dex in M_⋆, respectively. For each of these M_⋆-SFR bins, the median metallicity of the gas is computed. In Fig. 7, our M_⋆ − Z_g relations are compared to the measurements of Yates et al. (2012) and Andrews & Martini (2013).

The M_⋆ − Z_g relation extracted from our best model m₃ shows a clear variation with SFR at fixed stellar mass. As in Sect. 5.2 the absolute level measured by Andrews & Martini (2013) is in good agreement with our best model m₃. The dependence with the SFR is opposite for low- and high-mass regimes, as observed by Yates et al. (2012). The transition occurs around 3 × 10¹⁰ M_⊙ for our model, which is slightly higher than for the Yates et al. (2012) sample (10¹⁰ M_⊙). This shift of our model towards high stellar masses also appears when we compare our model with Andrews & Martini (2013) results. Our model shows a good agreement with their average relation in the first SFR bin corresponding to log (SFR) ∈ [ − 1.0: − 0.5]. However we also observe a systematic shift for the bins of higher SFR (log (SFR) ∈ [ − 0.5:1.0]) of the order of 0.5 dex. For log (SFR) > 1.0 the shift between our best model and the observational relation is smaller.

Galaxies extracted from our best model appear too massive as compared to observations. This is consistent with the discrepancy found between the observed main sequence of galaxies and that extracted from our best model (see Fig. 4). As in the M_⋆ − Z_g-SFR relation, the discrepancy for the main sequence is smaller at higher stellar mass (M_⋆> 10¹⁰).

In summary, our SFR-M_⋆ − Z_g relation indicates that in a given SFR bin our galaxies are too massive as compared to observations but, it shows a clear variation with SFR at fixed stellar mass. As observed by Yates et al. (2012), the low-SFR galaxies that are the most metal rich at low masses become the least metal rich at high masses.

6. Focus on Milky Way-like galaxies

In this section we study a sample of Milky Way-like galaxies generated with our best model m₃. The Milky Way is a SF galaxy. As a result of the number of constraints available, i.e. stellar mass distribution, ages, metallicities, abundance gradient, a large number of works are dedicated to the evolution of the MW in term of stellar mass assembly and metallicity enrichment process. In this context, MW-like galaxies are good references to test our new model. In this section, we base our study on two groups of galaxies: the MW cousins and the MW sisters. The first group is selected using selection criteria based on the dark matter halo and total stellar masses. Then, by applying additional criteria we select the MW sisters as a sub-sample of the MW cousins.

6.1. Milky Way cousins

The measurements compiled by Yin et al. (2009) or recently performed by Licquia & Newman (2015) indicate that the Milky Way has a total stellar mass between 3 × 10¹⁰ M_⊙ and 7 × 10¹⁰ M_⊙. In this stellar mass range, the agreement between the predicted M_⋆-SFR relation and the observed relation is satisfying (Fig. 4). Abundance matching methods (e.g. Behroozi et al. 2010) allow us to define the following range: 7 × 10¹¹ M_⊙<M_h< 3 × 10¹² M_⊙ for the dark matter halo hosting this range of stellar masses.

We randomly select 1440 halos in our simulation at z = 0 satisfying this dark matter halo mass criterion. For each one, the main branch of the merger tree is followed and saved with a high temporal resolution⁸. The whole merger tree is also saved⁹ to follow dark matter and baryonic properties of all progenitors.

From our initial selection of 1440 halos, we extract the MW-cousins sample using the stellar mass criterion: 3 × 10¹⁰<M_⋆< 7 × 10¹⁰ M_⊙. This selection reduces our initial sample from 1440 to 336 objects. The majority of the rejected galaxies have a stellar mass between 1 × 10¹⁰ M_⊙ and 3 × 10¹⁰ M_⊙. In the range of dark matter halo mass selected, our modelled galaxies are therefore on average less massive than predicted by the abundance matching method.

6.2. Milky Way sisters

In addition to its total stellar mass and dark matter halo mass, the Milky Way can be characterised by some other properties such as its SFR or its morphology. Based as previously on observational measurements we apply the following rules (Yin et al. 2009; Licquia & Newman 2015):

• star formation activity: 1 <SFR< 5 M_⊙/ yr; and

• morphology, defined as the ratio of the stellar mass hosted in disk to the total stellar mass of the galaxy: .

With these additional criteria, our MW-sisters sample contains 56 objects. The majority of the rejected galaxies have a smaller disk to total mass ratio ≃0.5 and a SFR lower than 1 M_⊙/ yr.

6.3. Global properties of Milky Way-like galaxies

To analyse our selection, we focus on six quantities: the time elapsed since the last major merger, stellar metallicity, stellar age, gas fraction, and the oxygen to hydrogen and iron to hydrogen relative abundances in the gas phase. The distributions of these quantities are shown in Fig. 8 for both MW cousins and MW sisters.

The panel a of the Fig. 8 shows the distribution of the time elapsed since the last major merger. While for the MW-cousins sample the width of the distribution is large, from 0.5 Gyr to 12.5 Gyr, the distribution associated with the MW sisters is tightened around 10.5 Gyr with a small scatter of 2 Gyrs. The morphology criterion that we used to select MW sisters induces a very high portion of the total stellar mass to be located in the disk of the galaxy and a small central bulge. Bulges are formed during major mergers; disk-dominated galaxies with a small bulge, such as MW-sister galaxies, encountered a major merger event a long time ago. This kind of morphology is rare. The majority of our MW-cousins galaxies are rejected from the MW-sisters sample because they have a spheroidal morphology or a smaller disk than observed in our Milky Way. This old major merger leads to an old stellar population in the central bulge. This is confirmed by the distribution of the luminosity-weighted age and the metallicity of the stellar populations shown in panels b and c of Fig. 8, respectively. Stars in the central bulge of MW-sisters galaxies have a luminosity-weighted ages between 10 and 12 Gyr. Their luminosity-weighted metallicity is also lower than that found for bulge stars of the MW-cousins sample, which is consistent with an earlier formation for the MW − sisters bulges. Stellar metallicities and ages of the stars located in the disk are similar in the two samples. The distributions associated with the MW cousins have a larger scatter than those predicted for MW-sisters sample, while their mean values are very close. For both MW-cousins and MW-sisters samples, the stellar age and metallicity in the disk are representative of the recent story, and the past evolution has only a very small impact on them. We compare our distributions of stellar metallicities with those of Bensby et al. (2014) and Bensby et al. (2013) for the disk and bulge stellar populations, respectively. In comparison to these reference values our stellar metallicities are too low by 0.2 dex in the bulge and by 0.5 dex in the disk. The discrepency in terms of stellar metallicity in the bulge component is small given the uncertainty about yields (up to a factor 2). The low metallicity of our stellar disk, and therefore the gap between the measure of Bensby et al. (2013) and our prediction, was already visible in Fig. 5. Indeed, in Fig. 5 our model predictions are offset by 0.5 dex as compared to the measurements of (Gallazzi et al. 2005).

In the three lower panels of Fig. 8, the distributions of the gas fraction are plotted together with those of the oxygen and iron abundance relative to the hydrogen in the gas phase of the disk. The distributions associated with the MW-cousins and MW-sisters sample exhibit similar average values and scatters. As for the stellar component of the discs, the past history does not have a strong impact on the gas properties.

We also compare our predictions to observed quantities. The measure of the gas fraction performed by Yin et al. (2009) in the Milky-Way is reported in panel d). Our modelled galaxies have a smaller gas fraction than observed. The instantaneous SFR of our MW-sisters galaxies is one of the criteria used to define the sample. The gas is therefore consumed with a rate similar to our Galaxy. A predicted gas fraction that is lower than observed seems therefore to indicate a lack of fresh gas. Our MW-sisters galaxies are hosted by dark matter halo of ≃10¹² M_⊙ in which accretion is dominated by the cooling process of the hot gas phase. Our cooling efficiency has to be increased to reproduced the observed gas fraction.

In panel e and f, we compare our predicted [O / H] _disk and [Fe / H] _disk abundances, respectively, with a compilation of measures performed in recently formed open clusters (<3 Gyr) (Pancino et al. 2010; Balser et al. 2011; Maderak et al. 2013). The stars in these open clusters are young. Their chemical composition are therefore close to that of the gas in which they formed.

Our predicted [O / H] _disk abundance shows an average value that is slightly higher than observed. In contrast, the average value of our [Fe / H] _disk abundance distribution is slightly lower than the observed value. The scatter of our predicted distribution is also smaller than observed. However, in the two cases, considering the dispersion and errors in measurements of the formation time of these open clusters, the predicted and observed abundances are in reasonable agreement.

In summary, MW sisters represents only 16% of the MW-cousins sample. Even if they host the same total stellar mass, the disk dominated morphology criterion strongly reduces the sample.

6.4. Comparison with detailed chemical evolution models

In this section we present both the star formation and the metal enrichment histories of our 56 MW-sisters galaxies. Based on their merger trees, we extract the median trends of both SFR and metal enrichment histories. The SFR history of each galaxy is computed by including the contribution of all its progenitors. Oxygen and iron abundances are followed only in the main disk’s progenitor of the galaxies.

Our median trends are compared to results coming from different detailed chemical evolution models aimed at reproducing MW properties in Fig. 9.

• Mollá & Díaz (2005) assume a universal rotation curve hypothesis and a multiphase chemical model and predict the radial mass distribution of 44 theoretical galaxies. They also provide the evolutions with time of oxygen and iron abundances. We compare our result with their model referenced as V = 200 and N = 4. This model assumes intermediate efficiencies for both: i) the conversion of the diffuse gas into molecular gas cloud and ii) the SFR;

• Boissier et al. (2010) give the SFR and chemical evolution of a set of disk galaxies. The grid of models is function of both the spin and maximal rotational velocity of the dark matter host halo. The properties of these modelled disks are consistent with observed properties for redshift included in 0 <z< 1;

• Kubryk et al. (2015) base their models on a multiphase chemical evolution. They also include a description of radial migration based on the analysis of an N-body plus hydrodynamic simulation. This radial flow mixes gas of metal-poor regions into metal-rich regions. Their model reproduces current values of main global observables of the MW disk and bulge.

This kind of models follows a classical approach for the mass assembly of the MW. These models do not follow the hierarchical formation of structures (as we do) but are constrained by additional observations (e.g. radial profile).

6.4.1. Star formation rate history

The median trend and individual SFR histories of our 56 MW-sisters are shown in the upper panel of Fig. 9. We observe strong short-time variations around the median trend. These variations are due to merger events or strong instantaneous baryonic accretion events that are not taken into account in Mollá & Díaz (2005), Boissier et al. (2010), or Kubryk et al. (2015).

Beyond these variations, our median SFR history reaches a maximum after roughly 1 Gyr of evolution (3. <z< 4.). Then the SFR slowly decreases until z = 0.

This median SFR history is in agreement with the results of chemical evolution models. This history is especially consistent with that of Kubryk et al. (2015), which takes radial transport inside the disk into account. Their predictions are compared favourably with a large suite of constraints in the Milky Way, such as age-metallicity relations.

In the first 5 Gyr of evolution, Boissier et al. (2010) predict a lower SFR than both Kubryk et al. (2015) and our model. In contrast, in the last 4 Gyr of evolution, Boissier et al. (2010) predict a slightly larger SFR than our model even if their values are still in the dispersion generated by all our individual SFH histories.

6.4.2. Metal enrichment history

In our model, the strong increase of the SFR, predicted just after the formation of the MW-sisters galaxies, leads to a large production of metals. After only ≃3 Gyr of evolution, both our median oxygen and iron abundances reach an asymptotic level. This saturation is not predicted by the three theoretical models, which instead predict a monotonic increase of both the oxygen and iron abundances. These different trends can be linked to the different feedback prescriptions and ejection processes assumed. In Boissier et al. (2010) and Kubryk et al. (2015) the new metals stay in the galaxy. These new metals are diluted with new fresh accreted gas, but their amount increases continuously. In our models the feedback processes produce winds that eject metals in the surrounding hot gas phase. Star formation and its associated feedback processes lead to an equilibrium between metals production and ejection. These two mechanisms explain the saturation in metal content and the asymptotic trend predicted by our models.

In the first Gyr of evolution, the SFR predicted by Boissier et al. (2010) is lower than the SFR predicted by Kubryk et al. (2015) or by our models. Consequently the oxygen and iron enrichment process is also slower during this period. However, the absolute value reached by Boissier et al. (2010) at z = 0 stays in good agreement with our results.

About Kubryk et al. (2015) results, despite a SFR history similar to our, the oxygen and iron abundances predicted reach larger values than in our MW-sisters galaxies. As discussed previously, this result is consistent with a scenario with an accretion of fresh gas without gas ejection.

Some other models have implemented a radial migration of gas and have tested the impact of these phenomenon on gas metallicity (e.g. Minchev et al. 2014). Their predictions for the evolution with time of [Fe/H] are fully compatible with our results (Minchev et al. 2014 their Fig. 12).

The chemical evolution predicted by Mollá & Díaz (2005) shows the best agreements with our predictions for both oxygen and iron abundances. The behaviours (shapes and slopes) and the absolute levels at z = 0 are especially consistent.

In the next section, we extend the comparison of our predicted star formation histories of MW-like galaxies to observational measurements.

6.5. Star formation rate history of our MW-like galaxy samples

In this section, we present the median star formation histories of both MW-cousins and MW-sisters samples. Our results are compared in Fig. 10 to the following observational measurements:

• van Dokkum et al. (2013) measure the stellar mass and the SFR of Milky Way-like galaxy progenitor out to z ≃ 3. Their analysis is based on 3D-HST and CANDELS Treasury surveys;

• Papovich et al. (2015) use ZFOURGE, a deep medium-band, near-infrared imaging survey coupled with HST, Spitzer, and Herschel images, to measure the evolution in stellar mass, SFR, and also morphologies of Milky Way progenitors.

In the two studies, the properties in stellar mass are deduced from abundance-matching techniques. As a prior, Papovich et al. (2015) and van Dokkum et al. (2013) adopt a final stellar mass of 5 × 10¹⁰ M_⊙ for the Milky Way. We recall that our stellar mass range criterion is 3 × 10¹⁰<M_⋆< 7 × 10¹⁰.

In Fig. 10, the two median trends extracted from our two samples show very similar behaviours. We note a slight deviation in the four last Gyrs of evolution between the two samples. Our MW-sisters galaxies are the most active; the SFR criterion generates a bias towards high SFR in the last four Gyrs of evolution. This result is consistent with the findings of Boissier et al. (2010). To match the stellar mass distribution of SF galaxy, these authors suggest that roughly 50% of the Milky Way-like galaxies¹⁰ have quenched their star formation since redshift 1 and therefore have a smaller SFR than the Milky Way. If this picture is correct, the Milky Way is more active than the majority of galaxies with similar stellar masses.

The star formation history measured by van Dokkum et al. (2013) is similar to what we obtain. Even if their maximum of activity is reached at lower redshift (z ≃ 2.), the absolute value of this peak is fully consistent with ours. Between z = 4 and z = 0, the slope of their star formation history is steeper than ours. Their measures are always included in our dispersion but reach the lowest values of our distribution at z = 0.

The results of Papovich et al. (2015) show a maximum for the star formation history around 1 <z< 2, which is lower than the maximum that we find. This maximum is higher (0.3 dex) than both the van Dokkum et al. (2013) measurements and our predicted value. Between z = 2. and z = 0.5, they measure a strong decrease of the SFR, which is steeper than the SFR we predict. The high SFR measured by Papovich et al. (2015) at z ≃ 2 is also in disagreement with van Dokkum et al. (2013) results, our model, and all the chemical evolution models presented previously. Indeed, this strong SFR peak would lead to a very strong metal enrichment process, which cannot be consistent with local (z = 0) measurements.

6.6. Redshift evolution of the metallicity, SFR, and stellar mass of MW-sisters progenitors

6.6.1. Metallicity of MW-sisters progenitors

The central panel of the Fig. 9 shows the evolution with time of the oxygen abundance predicted in the disk of the main progenitor of MW sisters. In this section, we extend this analysis at all the progenitors of our MW-sisters galaxies.

Figure 11 shows the relation M_⋆ −SFR relation extracted from our best model m₃ at four different redhifts. In this figure, all individual points show a MW-sister galaxy or a progenitor.The oxygen abundance in the disk’s gas phase is colour-coded. For each redshift, we can see that for a very large percentage of the sample the oxygen abundance is high (12 + log (O / H) ≥ 8.5). The metal enrichment process is therefore very efficient. We note that at z ≃ 3 the gas in the disks is already metal rich. By extending to all the progenitors, we confirm here the strong enrichment process observed in the main progenitor’s disk in Fig. 9. At z ≃ 2 and z = 3, we also confirm that on average at a given stellar mass, the higher the SFR, the lower the oxygen abundance. However, we can see some extreme points with low metallicities (12 + log (O / H) ≤ 7.7). These progenitors exhibit a low or very low star formation activity. These galaxies were extracted just after a merger event. During these events, the star formation activity is strongly increased (boost factor) and the gas is strongly consumed. After these violent events, in some cases, the star formation is strongly reduced because of a lack of fresh gas. These disks, without gas, accrete new gas, coming from the surrounding environment of the galaxy. This freshly accreted gas is generally poorer in metals in comparison to that previously consumed. After a merger, the metallicity can therefore strongly decrease.

6.6.2. Evolution onto the main sequence

In this section we address the evolution on the main sequence of the progenitors of our MW-sister galaxies. The evolution of the M_⋆ −SFR relation extracted from our best model m₃ is shown in Fig. 11. In this figure, we also report the relation measured by Peng et al. (2010) at z ≃ 0.09 and the evolution of the main sequence with redshift predicted by the empirical model of Béthermin et al. (2012) (for z> = 1), which reproduces the observed galaxy main sequence well between 1. <z< 3. (Karim et al. 2011; Rodighiero et al. 2011; Whitaker et al. 2012).

Measurements performed by Daddi et al. (2007), Karim et al. (2011), and Salmi et al. (2012) indicate that, at fixed stellar mass, the average SFR increases with redshift. This trend is also predicted by our model. However, as in Fig. 4, we observe that, at all redshifts and in the low-mass regime, the observed SFR is lower than that obtained by Béthermin et al. (2012). This is because of the strong regulation process applied on the star formation activity in low-mass structures. The non-sfg reservoir allows us to reproduce the stellar mass function well (Cousin et al. 2015a) and gives satisfying results for the gas metallicity (Fig. 5.2), but at fixed stellar mass and in the low-mass regime the star formation activity is lower than observed.

6.6.3. Milky Way-sister galaxies and the main sequence

At z = 0 a significant fraction of our MW sisters is above the mean trend of our predicted main sequence. As explained previously, SFRs of our MW-sisters galaxies are higher than the majority of our simulated galaxy with a similar stellar mass (Fig. 9). We also note that our MW sisters are close to the observed main sequence (Peng et al. 2010), but are still mostly below the observed main sequence.

At higher redshifts z ≥ 1, the majority of the MW-sisters progenitors are close to the M_⋆ −SFR relation extracted from our whole sample. The progenitors of our MW-sisters sample are therefore standard galaxies of our best model m₃. We note however that they are located below the observed main sequence at the corresponding redshift.

These two previous points are consistent with the result of Heinis et al. (2014). Indeed they found that galaxies would have a stellar mass that is too large if they stay on the observed main sequence from high redshift to z = 0. Therefore the MW-sisters progenitors cannot stay their entire life on the observed main-sequence relation, but reach this relation only at redshift close to 0. Currently the Milky Way is on the observed local main sequence. This is based on the currently observed main sequence. If in reality the low-mass galaxies (M_⋆< 10⁹) have lower SFR than currently measured (e.g. Whitaker et al. 2014), the observed main sequence should be closer to our predicted main sequence at high redshift. In these conditions the MW-sisters progenitors could have evolved mainly onto the main sequence.

7. Conclusions

In this work, we have focused our attention on the metallicity signature of the galaxies in the context of SAMs. To achieve this goal, we added a new chemodynamical extension to our SAM, which follows the metal enrichment process into both the stellar population and gas phase. We analysed the results of four different models based on two star formation scenarios and two accretion scenarios.

We show that the four models lead to similar M_⋆ − Z_⋆ relations at z ≃ 0. The different scenarios of gas accretion and star formation have therefore no significant impact on this relationship. The shape of our average M_⋆ − Z_⋆ relations is in good agreement with observations, but we observe a systematic shift to the high masses. This shift should be linked to the set of simple rules used in the gas-stars coupling.

In the M_⋆ − Z_g relationship, the two different models of star formation regulation lead to very different behaviours. In the low-mass regime, the oxygen abundance predicted by our classical model of star formation is too high compared with the measurements carried out in dwarf galaxies. Conversely, our model based on a strong regulation of the star formation shows better agreements. For a given star formation prescription, we note a slightly lower metallicity in the model based on bimodal accretion (i.e, metal-free cold-streams).

We then analysed the dependence of the M_⋆ − Z_g relation with the star formation activity. The absolute levels of our M_⋆ − Z_g −SFR relations, extracted from our best model m₃, are in good agreement with the recent measurements of Andrews & Martini (2013). In addition as observed by Yates et al. (2012) our low-SFR modelled galaxies, which are the most metal-rich at low masses, become the least metal-rich at high masses.

Based on dark matter halo mass and stellar mass, we selected a group of MW cousins. From this first sample, we extracted the MW sisters by applying both SFR and morphological criteria. This last criterion strongly reduce the sample. The MW-sisters sample represents only 16% of the initial MW-cousins sample. The SFR criterion also has an impact, we note that MW sisters are more active than their cousins. The median SFR history of our MW-sisters sample is consistent with observational measurements and also shows very good agreement with detailed chemical evolution models. The high values of oxygen and iron abundances predicted by our best model at high redshifts z ≤ 3.0 shows that the metal enrichment process is very efficient.

We finally addressed the evolution onto the main sequence of the progenitors of our MW-sister galaxies. At fixed stellar mass, our model reproduces the observed increase with the redshift of the SFR. However our predicted main sequence is lower (by 0.5 dex) than the observed main sequence in the low-mass regime. Our MW-sisters progenitors appear at high redshift, as the standard galaxies of our best model. However we have shown that our MW-sisters galaxies start their evolution below the observed main sequence and progressively reach this observed relation at z ≃ 0.

---

Acknowledgments

Authors thank the Centre National d’Études Spatiales for its financial support. We thank Laura Portinari for very useful discussions. We thank the referee for his helpful report and his comments.

References

Appendix A: From SSPs to realistic star formation histories: The history tables

Appendix A.1: The SSP ejecta rate table: Ej - TAB

Based on Eq. (15) we build a 3D (age, metallicity, and element) table: Ej - TAB. We use 130 different age bins and 7 initial stellar metallicity bins. The chosen stellar ages values correspond to a sub-sample (Δt_min = 1 Myr) of the SSP-age list available in the BC03 stellar spectrum library. We assume that ejecta rates are constant between two successive stellar ages. The Ej - TAB is computed by assuming a given IMF. We use Chabrier (2003) but this table is available for Salpeter (1955), Scalo (1998), and Kroupa (2001) IMFs. In order to be consistent with the mass conservation, in each age bin, the instantaneous ejecta rate is normalised to m_⋆ + m_remnent = 1 M_⊙.

The Ej - TAB is the first elementary block used to compute the stellar evolution and the metal enrichment process associated with the complex star formation history that occurs in our simulated galaxies.

Appendix A.2: The stellar mass history table: SMH - TAB

For each galaxy, the stellar mass assembly is followed using a 2D (age and metallicity) table: SMH - TAB. At each time of the evolution process, the cell gives the mass of stars (m = m_⋆ + m_remnent) with an age: a_i<age<a_{i + 1} and with an initial stellar metallicity . A clock clk is associated with each cell of the SMH - TAB. This clock runs as long as the cell contains stellar mass. The SMH - TAB is the second elementary block used to compute the stellar evolution and the metal enrichment process associated with the complex star formation history.

Appendix A.2.1: Increasing stellar mass

During a star formation episode, a total mass of gas M_g with a metallicity Z_g (metal mass fractions) is transformed into new stars. The mass of each main-ISM element contained in this mass of gas is also known. The mass M_g has to be added to the first age bin of SMH - TAB and distributed into the different available initial stellar metallicity bins. We apply a set of rules:

• i)

  No stars can be formed with an initial stellar metallicity higher than the input average gas metallicity.

• ii)

  To be consistent with the composition of the total stellar ejecta associated with a given initial stellar metallicity, no interpolation between stellar metallicity bins is performed.

• iii)

  For each initial stellar metallicity bin, the initial chemical composition is fixed and given by Karakas (2010).

To be consistent with the amount of each element contained in the star-forming gas, the distribution of the gas between the different metallicity bins is carried out by maximising the mass affected in the bin that is closest to (but smaller than) the average gas metallicity Z_g. By following the previous rules, the chemical composition of the stellar mass actually formed is different from the input star-forming gas. To respect the mass conservation for each element followed, a residual mass is returned to the ISM¹¹.

Appendix A.2.2: Evolution of the stellar population

The evolution of a given stellar population is then performed following a adaptive time-steps scheme applying time-steps Δt smaller or equal to 1 Myr. For each bin , the evolution is carried out as follows:

• i)

  For each element tracked by the library, the massm_elt, ejected by the stellar population hosted in the bin is re-incorporated into the ISM. This mass is given by (A.1)

• ii)

  In a second time-step, the mass contained in a given bin of can be (or not) transferred to the next (older) age bin: a_{i + 1}. Depending on the clock linked to this bin, two cases are possible:

  • i) If clk, no mass transfer is performed and we just add the current time-step to the clock, (A.2)

  • ii) If , all the mass in the bin is transferred into the next age bin . This transfer of mass towards the next age bin leads to a change of its clock. The new value is set to (A.3)

The SMH - TAB is updated starting with the oldest age bin containing mass and, afterwards, by browsing all younger bins. The evolution procedure ends by the younger age bin: 0 Myr in which a new star formation episode can be taken into account.

Appendix A.3: Instantaneous SN rate: SN - TAB

We introduce a realistic SNII + SNIa even rate. This rate is used to compute SN feedback.

Based on Eqs. (9) and (10), we build a 2D (age and metallicity) table: SN - TAB: (A.4)In this definition the first term is the death rate of a massive star applied for a massive star with a lifetime τ_m = a_i. The second term, is the instantaneous SNIa rate computed at a time t = a_i. We assume that the global SN rate is constant between two successive ages. As for the instantaneous ejecta rate, at a given age, the SN rate is normalised to m_⋆ + m_remnent = 1 M_⊙.

Figure A.1 shows the instantaneous SNII and SNIa event rates for different stellar metalicities. This SN rates are associated with a SSP of 1 M_⊙ distributed through a Chabrier IMF.

By using SN - TAB and SMH - TAB the instantaneous SN rate, η_sn, used in Eq. (6) is simply given by (A.5)

Appendix B: Data repository

Appendix B.1: Model outputs

Outputs from all models presented here are available at the CDS. Data are distributed under fits-formatted files and are therefore compatible with the TOPCAT software¹².

We also provide the median star formation history, median gas accretion history, and metal enrichment histories (12 + log[O/H] and [Fe/H]) associated with our MW-sisters sample (Fig. 10). These data are saved in ascii-formated files.

Appendix B.2: Chemodynamical library

The study presented here is mainly based on a new chemodynamical library. We provide a set of two files associated with this library:

• mass_loss_rates_IMF.fits files, which contain the ejecta rate [M_⊙/Gyr] associated with a SSP of 1 M_⊙ for the six different main-ISM elements followed: ¹H, ⁴He, ¹²C, ¹⁴N, ¹⁶O, and ⁵⁶Fe. The total ejected mass and metal ejected mass are also provided (Eq. (15) and Fig. 2);

• SN_rates_IMF.fits files, which contain the total SN rate (SNII + SNIa [nb/Gyr]) associated with a SSP of 1 M_⊙ and with the individual contribution of SNII and SNIa (Fig. A.1).

These files are available for four different IMFs: Salpeter+55, Chabrier+03, Kroupa+93 and Scalo+98. Both ejecta rates and SN rates are computed for the complete list of stellar ages provided in the BC03 spectra library. They are saved in fits-formated files and structured with different extensions corresponding to the different initial stellar metallicity bins.

All Tables

All Figures