The Origin of the Stellar Mass-Stellar Metallicity Relation in the Milky Way Satellites and beyond

TL;DR

This study develops a chemical evolution model to explain the universal stellar mass-stellar metallicity relation (MZR) across various galaxy types and masses, emphasizing the roles of gas inflows, outflows, and feedback efficiency.

Contribution

It introduces a model linking the MZR to power-law scaling of cold gas and feedback efficiency, supported by semi-analytical models, highlighting the relation's universality.

Findings

The MZR is shaped by decreasing cold gas fraction with increasing stellar mass.

Stellar feedback efficiency decreases as stellar mass increases.

The model's MZR range aligns with observed universality across galaxy types.

Abstract

Observations and semianalytical galaxy formation and evolution models (SAMs) have suggested the existence of a stellar mass-stellar metallicity relation (MZR), which is shown to be universal for different types of galaxies over a large range of stellar masses ($M_*\sim 10^3$-$10^{11}M_\odot$) and dark matter (DM) halo masses ($M_{\rm halo}\sim 10^9$-$10^{15}h^{-1}M_\odot$). In this work, we construct a chemical evolution model to investigate the origin of the MZR, including both the effects of gas inflows and outflows in galaxies. We solve the MZR from the chemical evolution model, by assuming that the cold gas mass ($M_{\rm cold}$) and the stellar feedback efficiency ($\beta$) follow some power-law scaling relationships with $M_*$ during the growth of a galaxy, i.e., $M_{\rm cold}\propto M_*^{\alpha_{\rm gs}}$ and $\beta\propto M_*^{\alpha_{\beta{\rm s}}}$. We use the SAM to obtain these power-law scaling relations, which appear to be roughly universal over a large range of stellar masses for both satellites and central galaxies within a large range of halo masses. The range of the MZRs produced by our models is in a narrow space, which provides support to the universality of the MZRs. The formation of the MZR is a result caused jointly by that the cold gas fraction decreases with increasing $M_*$ and by that the stellar feedback efficiency decreases with increasing $M_*$ in the galaxy growth, and the exponent in the MZR is around $-\alpha_{\beta{\rm s}}$ or $1-\alpha_{\rm gs}$. The MZR represents an ''average'' evolutional track for the stellar metallicity of a galaxy. The comparison of our model with some previous models for the origin of MZRs is also discussed.

Figures

Tables

Table 1. Table 1 Best-fit results for the median M cold subscript 𝑀 cold M_{\rm cold} – M ∗ subscript 𝑀 M_{*} scaling relations

Objects	$M_{\text{halo}}/\left(h^{-1}{M}_{\odot }\right)$	$\log \left(M_{\ast }/M_{\odot }\right)$	${\alpha }_{\text{gs}}$	$b_{\text{gs}}$	${\chi }^2$	$N_{\text{bin}}$	$N_{\text{tree}}$	$N_{\text{galaxy}}$	$M_{\mathrm{res}}/\left(h^{-1}{M}_{\odot }\right)$
	$1\times {10}^{10}$	$3-4$	$0.24\pm 1.75$	$5.8\pm 1.6$	-	2	100	9	${10}^4$
	$5\times {10}^{10}$	$3-7$	$0.79\pm 0.05$	$5.5\pm 0.2$	0.9	8	100	62	${10}^4$
satellite	$1\times {10}^{11}$	$3-7.5$	$0.76\pm 0.06$	$5.5\pm 0.2$	2.6	9	100	165	${10}^4$
progenitors	$5\times {10}^{11}$	$3-9$	$0.68\pm 0.03$	$5.7\pm 0.1$	4.3	14	100	2163	${10}^4$
	$1\times {10}^{12}$	$3-10$	$0.61\pm 0.03$	$5.9\pm 0.1$	15.7	14	100	5970	${10}^5$
	$2\times {10}^{12}$	$3-10$	$0.66\pm 0.03$	$5.7\pm 0.1$	3.3	14	100	15265	${10}^5$
	$1\times {10}^{13}$	$3-10$	$0.65\pm 0.03$	$5.7\pm 0.1$	2.9	14	100	102164	${10}^6$
	$1\times {10}^{14}$	$3-10$	$0.65\pm 0.03$	$5.8\pm 0.1$	2.0	14	10	85563	${10}^7$
	$1\times {10}^{15}$	$3-10$	$0.61\pm 0.03$	$6.0\pm 0.1$	0.9	14	10	543679	${10}^8$
	$1\times {10}^{10}$	$3.5-6.5$	$0.69\pm 0.09$	$5.65\pm 0.22$	0.2	6	100	100	${10}^4$
	$5\times {10}^{10}$	$3.5-8.5$	$0.66\pm 0.03$	$5.82\pm 0.11$	2.4	10	100	100	${10}^4$
central	$1\times {10}^{11}$	$3.5-9$	$0.66\pm 0.03$	$5.84\pm 0.10$	1.9	11	100	100	${10}^4$
galaxies	$5\times {10}^{11}$	$3.5-10$	$0.64\pm 0.02$	$5.85\pm 0.09$	1.3	13	100	100	${10}^4$
	$1\times {10}^{12}$	$3.5-10$	$0.62\pm 0.03$	$5.89\pm 0.09$	0.9	13	100	100	${10}^5$
	$2\times {10}^{12}$	$3.5-10$	$0.61\pm 0.03$	$5.91\pm 0.10$	1.6	13	100	100	${10}^5$
	$1\times {10}^{13}$	$4-10$	$0.59\pm 0.02$	$5.97\pm 0.07$	1.1	12	100	100	${10}^6$
	$1\times {10}^{14}$	$5-10$	$0.57\pm 0.04$	$5.99\pm 0.12$	2.5	10	10	10	${10}^7$
	$1\times {10}^{15}$	$4.5-10$	$0.59\pm 0.01$	$5.94\pm 0.03$	9.2	11	10	10	${10}^8$

Table 2. Table 2 Best-fit results for the median β 𝛽 \beta - M ∗ subscript 𝑀 M_{*} scaling relations

Objects	$M_{\text{halo}}/h^{-1}{M}_{\odot }$	$\log \left(M_{\ast }/M_{\odot }\right)$	${\alpha }_{\beta s}$	$b_{\beta s}$	${\chi }^2$	$N_{\text{bin}}$	$N_{\text{tree}}$	$N_{\text{galaxy}}$	$M_{\mathrm{res}}/\left(h^{-1}{M}_{\odot }\right)$
	$1\times {10}^{10}$	$3-4$	$-0.18\pm 0.70$	$2.37\pm 0.58$	-	2	100	9	${10}^4$
	$5\times {10}^{10}$	$3-7$	$-0.23\pm 0.03$	$2.44\pm 0.12$	0.04	8	100	62	${10}^4$
satellite	$1\times {10}^{11}$	$3-7.5$	$-0.23\pm 0.03$	$2.47\pm 0.10$	0.10	9	100	165	${10}^4$
progenitors	$5\times {10}^{11}$	$3-9$	$-0.26\pm 0.03$	$2.55\pm 0.09$	0.03	14	100	2163	${10}^4$
	$1\times {10}^{12}$	$3-10$	$-0.30\pm 0.02$	$2.65\pm 0.08$	0.28	14	100	5970	${10}^5$
	$2\times {10}^{12}$	$3-10$	$-0.29\pm 0.03$	$2.62\pm 0.09$	0.07	14	100	15265	${10}^5$
	$1\times {10}^{13}$	$3-10$	$-0.30\pm 0.03$	$2.65\pm 0.09$	0.09	14	100	102164	${10}^6$
	$1\times {10}^{14}$	$3-10$	$-0.31\pm 0.03$	$2.68\pm 0.09$	0.12	14	10	85563	${10}^7$
	$1\times {10}^{15}$	$3-10$	$-0.25\pm 0.02$	$2.46\pm 0.06$	0.46	14	10	543679	${10}^8$
	$1\times {10}^{10}$	$3.5-6.5$	$-0.13\pm 0.08$	$2.30\pm 0.18$	0.04	6	100	100	${10}^4$
	$5\times {10}^{10}$	$3.5-8.5$	$-0.23\pm 0.03$	$2.49\pm 0.09$	0.48	10	100	100	${10}^4$
central	$1\times {10}^{11}$	$3.5-9$	$-0.28\pm 0.02$	$2.59\pm 0.07$	0.47	11	100	100	${10}^4$
galaxies	$5\times {10}^{11}$	$3.5-10$	$-0.34\pm 0.02$	$2.73\pm 0.05$	0.24	13	100	100	${10}^4$
	$1\times {10}^{12}$	$3.5-10$	$-0.36\pm 0.02$	$2.74\pm 0.05$	0.09	13	100	100	${10}^5$
	$2\times {10}^{12}$	$3.5-10$	$-0.40\pm 0.02$	$2.82\pm 0.06$	0.04	13	100	100	${10}^5$
	$1\times {10}^{13}$	$4-10$	$-0.45\pm 0.02$	$2.89\pm 0.06$	0.12	12	100	100	${10}^6$
	$1\times {10}^{14}$	$5-10$	$-0.58\pm 0.03$	$3.26\pm 0.10$	0.10	10	10	10	${10}^7$
	$1\times {10}^{15}$	$4.5-10$	$-0.58\pm 0.02$	$3.21\pm 0.04$	0.53	11	10	10	${10}^8$

Equations

$$\psi =M_{\mathrm{cold}}/{\tau }_{\ast }$$

$${\tau }_{\ast }={ϵ}_{\ast }^{-1}\tau {}_{\mathrm{disk}}(V{}_{\mathrm{disk}}/200\mathrm{km}{\mathrm{s}}^{-1}{)}^{{\alpha }_{\ast }},$$

$$dM{}_{\mathrm{reheat}}=\beta \psi dt,$$

$$\beta =(V{}_{\mathrm{disk}}/V{}_{\mathrm{hot}}{)}^{-\alpha {}_{\mathrm{hot}}},$$

$$d{E}_{\mathrm{SN}}-\frac{1}{2}{V}_{\mathrm{vir}}^{2}dM{}_{\mathrm{reheat}}>0,$$

$$\beta {}_{\mathrm{E}}\equiv d{E}_{\mathrm{SN}}/(\frac{1}{2}{V}_{\mathrm{vir}}^2\psi dt)=ϵ{}_{\mathrm{halo}}(V_{\mathrm{vir}}/V_{\mathrm{SN}}{)}^{-2},$$

$${\dot{M}}_{\ast }=(1-R)\psi ,$$

$${\dot{M}}_{\text{cold}}={\dot{M}}_{\text{cool}}-(1-R+\beta )\psi ,$$

$${\dot{M}}_{\ast }^Z=(1-R)Z_{\mathrm{cold}}\psi ,$$

$${\dot{M}}_{\text{cold}}^Z={\dot{M}}_{\text{cool}}{Z}_{\text{hot}}+[p-(1-R+\beta )Z_{\text{cold}}]\psi ,$$

$$Z_{\ast }(M_{\ast })$$

$$=$$

$${\dot{M}}_{\text{cool}}={\dot{M}}_{\text{cold}}+\frac{1-R+\beta }{1-R}{\dot{M}}_{\ast },$$

$$d(M_{\text{cold}}{Z}_{\text{cold}})-Z_{\text{hot}}\cdot d{M}_{\text{cold}}=\frac{p}{1-R}d{M}_{\ast }$$

$$-\frac{1-R+\beta }{1-R}(Z_{\text{cold}}-Z_{\text{hot}})d{M}_{\ast }.$$

$$d\left(M_{\text{cold}}{Z}_{\text{cold}}\right)=\frac{p}{1-R}d{M}_{\ast }-\frac{1-R+\beta }{1-R}{Z}_{\text{cold}}d{M}_{\ast }.$$

$$Z_{\text{cold}}=\frac{M_{\text{cold},0}{Z}_{\text{cold},0}}{M_{\text{cold}}}\exp \left(-{\int }_{M_{\ast ,0}}^{M_{\ast }}\frac{1-R+\beta }{1-R}\frac{d{M}_{\ast }^{\prime }}{M_{\text{cold}}}\right)$$

$$+\frac{p}{1-R}\frac{1}{M_{\text{cold}}}\exp \left(-{\int }_c^{M_{\ast }}\frac{1-R+\beta }{1-R}\frac{d{M}_{\ast }^{{}^{\prime }}}{M_{\text{cold}}}\right)$$

$$\cdot {\int }_{M_{\ast ,0}}^{M_{\ast }}\exp \left({\int }_c^{M_{\ast }^{{}^{\prime }}}\frac{1-R+\beta }{1-R}\frac{d{M}_{\ast }^{{}^{\prime \prime }}}{M_{\text{cold}}}\right)d{M}_{\ast }^{{}^{\prime }},$$

$$d{M}_{\ast }=-\frac{1-R}{1-R+\beta }d{M}_{\text{cold}}.$$

$$Z_{\text{cold}}=Z_{\text{cold,0}}+\frac{p}{1-R+\beta }\ln \frac{M_{\text{cold,0}}}{M_{\text{cold}}},$$

$$Z_{\ast }$$

$$M_{\ast }\cdot Z_{\ast }=M_{\ast ,0}\cdot Z_{\ast ,0}+\left(M_{\ast }-M_{\ast ,0}\right)Z_{\text{cold,0}}$$

$$+\left(M_{\ast }-M_{\ast ,0}\right)\frac{p}{1-R+\beta },$$

$$Z_{\ast }-Z_{\ast ,0}=(\frac{p}{1-R+\beta }+Z_{\mathrm{cold},0}-Z_{\ast ,0})$$

$$\times (1-M_{\ast ,0}/M_{\ast }).$$

$$Z_{\text{cold}}=\frac{p}{1-R+\beta }\ln f_{\text{gas}}^{-1},$$

$$Z_{\ast }=\frac{p}{1-R+\beta }\left(1+\frac{f_{\mathrm{gas}}}{1-f_{\mathrm{gas}}}\ln f_{\mathrm{gas}}\right).$$

$$Z_{\ast }=\frac{p}{1-R+\beta },$$

$$Z_{\text{cold}}=\frac{p}{1-R}\ln f_{\text{gas}}^{-1}$$

$$Z_{\ast }=\frac{p}{1-R}\left(1+\frac{f_{\mathrm{gas}}}{1-f_{\mathrm{gas}}}\ln f_{\mathrm{gas}}\right).$$

$$M_{\text{cold}}=(\frac{M_{\ast }+\Pi }{M_{\infty }+\Pi })(M_{\infty }-M_{\ast }),$$

$$\Pi =\frac{M_0{M}_{\infty }}{M_{\infty }-M_0},$$

$$Z_{\text{cold}}=\frac{p}{1-R}\frac{1}{M_{\text{cold}}}\exp \left(-{\int }_c^{M_{\ast }}\frac{d{M}_{\ast }^{\prime }}{M_{\text{cold}}}\right)$$

$$\cdot {\int }_0^{M_{\ast }}\exp \left({\int }_c^{M_{\ast }^{\prime }}\frac{d{M}_{\ast }^{\prime \prime }}{M_{\text{cold}}}\right)d{M}_{\ast }^{\prime }.$$

$$Z_{\text{cold}}=\frac{p}{1-R}{\left(\frac{M_{\infty }+\Pi }{M_{\ast }+\Pi }\right)}^2$$

Full text

††slugcomment: Accepted to The Astrophysical Journal

1. Introduction

The stellar mass and the stellar metallicity are among the most basic physical properties of a galaxy. Both observations and semi-analytical galaxy formation and evolution model (SAM) simulations reveal that the stellar metallicity correlates with the stellar mass in both dwarf satellites and normal galaxies. For example, Kirby et al. (2013) find that there exists a universal stellar mass–stellar metallicity relation in the Milky Way/M31 dwarf satellites and some other dwarf irregular galaxies in the Local Group, $Z_{*}\propto M_{*}^{0.30\pm 0.02}$. The SDSS spectra of over 40,000 galaxies show that local galaxies with stellar masses in the range of $10^{9}\lesssim M_{*}/{\,M_{\odot}}\lesssim 10^{12}$ also indicate that stellar masses have a correlation with stellar metallicities or gas metallicities (Gallazzi et al. 2005, 2006; Tremonti et al. 2004; Panter et al. 2008; Mannucci et al. 2010; González Delgado et al. 2014; Lee et al. 2006). Hereafter, for simplicity, we refer to the correlation between the stellar mass and stellar metallicity as “MZR” (see also review in Maiolino & Mannucci 2019).

The MZR in the dwarf satellites of the Milky Way ($10^{3}{\,M_{\odot}}<M_{*}<10^{8}{\,M_{\odot}}$ or $10^{9}{\,M_{\odot}}$) was reproduced in the SAMs (e.g., Font et al. 2011; Li et al. 2010; Lu et al. 2014, 2017; Hou et al. 2014, hereafter HYL14). By using SAMs, Yates et al. (2013) obtain the stellar mass–cold gas metallicity relation for central galaxies with relatively high stellar mass range $M_{*}\sim 10^{9}$–$10^{11}{\,M_{\odot}}$, and Somerville et al. (2015) present the stellar mass–cold gas metallicity relation for central galaxies with stellar mass range $M_{*}\sim 10^{7}$–$10^{11}{\,M_{\odot}}$. Xia & Yu (2019) use the SAM to produce the MZRs in galaxies with a larger range of stellar masses $M_{*}\sim 10^{3}$–$10^{11}{\,M_{\odot}}$, including both dwarf satellites and normal galaxies within a larger range of host halo masses ($M_{\rm halo}\sim 10^{9}-10^{15}h^{-1}{\,M_{\odot}}$). In this paper, we investigate the origin of the universal MZRs, by constructing an analytical chemical evolution model of galaxies and applying the reference results from the SAM. The study of the origin will help us to understand the physical processes involving stellar formation and evolution, gas inflow and outflow, and chemical evolution in galaxies.

One key to understanding the origin of the MZR is to model the galaxy metal enrichment processes. The modeling of galaxy metal enrichment has evolved from a simple “closed-box” model to relatively complex ones including the effects of gas inflow and outflow in the past fifty years. In the “closed-box” model (Schmidt, 1963; Talbot & Arnett, 1971; Searle & Sargent, 1972), an isolated galaxy formation process starts from pure metal-free gas, and no gas escapes from or is accreted onto the galaxy. Some later modified models such as the “leaky-box” model, the “accretion” model (Lynden-Bell, 1975), and the “pre-enriched” model (Pagel, 1997), include gas outflow from a galaxy, gas inflow, and initial non-zero metallicity of inflowing gas, respectively. The effects of both gas inflow and outflow are included in the models of Finlator & Davé (2008), Lilly et al. (2013), and Peng & Maiolino (2014).

Metal-enriched gas outflow is proposed to play an important role in shaping the MZR. For example, the positive correlation in the MZR can be produced in the following cases. (1) The outflows can be stronger for small galaxies due to their shallower gravity potential wells, and thus small galaxies lose more metals than larger ones (Dekel & Silk, 1986). (2) The outflow wind strength is ubiquitous for all the galaxies, but the wind metallicity is larger than the metallicity in the inter-stellar medium (ISM) of smaller galaxies (Kobayashi et al., 2007). Though gas inflows from host halos may dilute the metallicity in the ISM, Dalcanton (2007) points out that the inflow dilution in the metallicity cannot contribute much to the low effective yield of small galaxies, neither can the metal-unenriched outflow. Finlator & Davé (2008) propose an equilibrium metallicity to explain the MZR, where the metallicity is obtained by assuming that the gas inflow rate is balanced by the gas consumption rate caused by star formation plus the gas outflow rate. Lilly et al. (2013) relax the assumption of gas processing equilibrium, but assume an enrichment equilibrium to obtain the equilibrium metallicity of the gas, where metal dilution balances metal enrichment in gas metallicity (see also a review in Finlator 2016).

In addition to the above general proposals to shape the MZR by using gas outflow and inflow, there are also some models specifically involving star formation rates and histories to explain the MZR. For example, the positive correlation in the MZR is explained in Brooks et al. (2007) by relatively low star formation efficiencies in low-mass galaxies and in Vale Asari et al. (2009) by that the majority of the stars form earlier in large galaxies and the metal enrichment in large galaxies proceeds earlier. Köppen et al. (2007) suggests that the positive correlation in the MZR can be produced, if large galaxies have a top-heavy initial stellar mass function and there are more massive stars in large galaxies to produce large effective yields.

In this paper, we present a chemical evolution model to explain the origin of the MZR. Previously, we have shown in HYL14 that the MZR in the satellites of MW-like galaxies produced from the SAM matches the observation, given a specific set of parameters (i.e., the fiducial model parameters in HYL14), such as on the physical processes of feedback and the reionization of the universe. Then we fix these parameters and explore the variations as a result of changing halo masses ($M_{\rm halo}\sim 10^{9}$–$10^{15}h^{-1}{\,M_{\odot}}$), central galaxies ($M_{*}\sim 10^{3}$–$10^{11}{\,M_{\odot}}$), and dwarf satellites, and we find that the MZR is quite independent of changing these variations as long as the fiducial model parameters are fixed (Xia & Yu, 2019). In this work, by using the results obtained from the SAM to serve as the quantitative reference to some physical galaxy properties involved in the chemical evolution model, we present an explanation to the universal MZR. This paper is organized as follows. In Section 2, we give a short overview of the SAM used in this work (for more details, see HYL14 and Xia & Yu 2019). In Section 3, we present our analytic model in the galaxy metallicity evolution. In this chemical evolution model, the evolution of the cold gas mass and the stellar feedback efficiency with the growth of the stellar mass are involved and assumed to follow some power-law scaling relationships. In Section 4, we show these power-law scaling relations obtained from the SAM, and apply them to the chemical evolution model constructed in Section 3. Then we use the chemical evolution model results to give an explanation to the universal MZRs. A summary is given in Section 5.

In this paper we set the Hubble constant as $H_{0}=100\,h{\rm\,km\,s^{-1}}{\rm\,Mpc}$, and the cosmological model used is $(\Omega_{\rm m},\Omega_{\Lambda},h,\sigma_{8})=(0.25,0.75,0.70,0.90)$.

2. The semi-analytical galaxy formation and evolution model (SAM)

In this section, we briefly describe the SAM used in this work, which will provide some quantitative reference to some physical properties involved in the chemical evolution model and help to explore the origin of the MZR of MW dwarf satellites and beyond.

The SAM used in this work is based on the GALFORM (Cole et al. 2000; see also White & Frenk 1991; Kauffmann et al. 1993; Somerville & Primack 1999; Somerville et al. 2008), with several modifications by HYL14 (e.g., on gas cooling, stellar feedback, and metal enrichment). The model includes planting DM halo merger trees and then incorporating semi-analytical recipes of the baryonic physical processes into the halo merger trees, to trace the hierarchical galaxy formation and evolution until redshift $z=0$ (Cole et al. 2000; see also White & Frenk 1991; Kauffmann et al. 1993; Somerville & Primack 1999; Somerville et al. 2008. The DM halo merger trees used in this work are generated from a modified version of the extended Press-Schechter formula (Press & Schechter, 1974; Bond et al., 1991; Lacey & Cole, 1993; Somerville et al., 2008), developed by Parkinson et al. (2008) (see also White & Frenk 1991; Kauffmann et al. 1993; Somerville & Kolatt 1999; Cole et al. 2000). For each DM halo mass at $z=0$, we use the Monte Carlo method to generate a number of merger trees (e.g., from 10 to 1000) in order to do a statistical study. The progenitor halos of each tree planted in this work can be dated back to high redshift $z=20$. As shown in Xia & Yu (2019), to address the universality of the MZR, we studied galaxies and DM halos with different masses by using the SAM and the Monte Carlo method, so that the physical properties of galaxies with different masses and their satellites in different environments can be explored.

The SAM recipes are modified in HYL14 to make the model suitable for studying low-mass dwarf galaxies/satellites. Some relevant recipes are briefly described below. More details can be found in HYL14 (or Section 2 in Xia & Yu 2019).

• •

Gas cooling: The initial temperature of the hot gas in the newly formed halo is assumed to be the same as the virial temperature of the DM halo. The cooling rate of the hot gas depends on its temperature, mass density, and metallicity, and we adopt the cooling rate for atomic cooling from Sutherland & Dopita (1993). The detailed cooling recipe is based on the GALFORM model by calculating how much hot halo gas can cool down and “free-fall” to the halo center within a certain time, and modified by HYL14 by including the mixing of the gas reheated by stellar feedback with hot halo gas (mainly occurred in low-mass galaxies).

The molecular hydrogen cooling of pristine gas is also included in mini-halos before the completeness of the reionization of the universe (see Equations 21-29 in Benson 2010; see also Galli & Palla 1998). As pointed out in HYL14 that the MZR of the satellites in MW-like galaxies is not sensitive to the molecular hydrogen cooling process.

• •

Star formation and stellar feedback: The star formation rate is proportional to the mass of cold gas in the disk, given by

[TABLE]

where $\psi$ is the instantaneous star formation rate, $\tau_{*}$ is the star formation time scale, $\tau{{}_{\rm disk}}$ is the dynamic time scale of the galaxy disk, $V{{}_{\rm disk}}$ is the disk rotation velocity, and $\epsilon_{*}=0.005$ and $\alpha_{*}=-1.5$ are two parameters (see eq. 4.14 in Cole et al. 2000).

SNe and stellar winds will heat up the cold gas in the disk and possibly expel it out of the galaxy as outflows. During time interval $dt$, the reheated gas mass is given by:

[TABLE]

where $\beta$ is the feedback efficiency, and $V{{}_{\rm hot}}$ and $\alpha{{}_{\rm hot}}$ are the two parameters defining the strength of the feedback. The cold gas mass that can be expelled from the disk is constrained by the energy released by SNe and coupling to the IGM (see also Guo et al. 2011), with the following inequality:

[TABLE]

where $dE_{\rm SN}=\epsilon{{}_{\rm halo}}\times\frac{1}{2}V_{\rm SN}^{2}\psi dt$ is the total energy released by SNe and coupling to the IGM during time $dt$, $V_{\rm vir}$ is the virial velocity of the halo, $\frac{1}{2}V_{\rm SN}^{2}$ is the total energy released per unit mass by SNe with $V_{\rm SN}=630{\rm\,km\,s^{-1}}$, and $\epsilon{{}_{\rm halo}}=0.05$ is the fraction of the energy that couples to the cold gas in the disk.

If inequality (5) is satisfied, then all the reheated gas can be expelled out of the disk (part or all of the reheated gas can even escape the DM halo). In this case, $\beta$ in the chemical evolution model described in Section 3 (Equations 8 and 10) can still be obtained through Equation (4). If inequality (5) is not satisfied (normally for low-mass galaxies or the progenitors of satellites, e.g., $M_{*}\lesssim 10^{8}{\,M_{\odot}}$), the mass of the reheated gas that can be expelled out of the disk during time interval $dt$ is given by $\beta{{}_{\rm E}}\psi dt$, where

[TABLE]

and the expelled gas stays in the halo (see more details in section 2 in HYL14). In this case, $\beta$ in the chemical evolution model described in Section 3 is $\beta=\beta{{}_{\rm E}}$.

Note that for satellites, $\beta$ in the chemical evolution model described in Section 3 is also obtained through Equation (4) for the following reason. In the assembling history of a halo, a satellite of a central galaxy may be the host galaxy of a small isolated halo before it fell into the big halo at an early time. After it fell into the big halo, the original halo of the satellite can be largely tidally disrupted along its motion in the big host halo, and thus we assume that all the reheated gas from satellites (with mass expected by Equations 3 and 4) is expelled into the big host halo. The tidal stripping and disruption of the stellar and cold gas components of the satellites are not considered in our model, as they are located in a smaller central region compared with their original halo size.

The parameters in Equation (4) are set to be $V{{}_{\rm hot}}=200{\rm\,km\,s^{-1}}$ and $\alpha{{}_{\rm hot}}=3.2$, i.e. the fiducial model in HYL14, which can reproduce some observational properties of the MW dwarfs (e.g., the satellite luminosity function, the MZR) better than the other models (see also Figure 9 in HYL14 for the changes to the MZRs caused by other sets of the two parameters). Although these parameters were set to be compatible with the observational properties of the dwarf satellites in HYL14, the MZR for galaxies covering a large stellar mass range ($10^{3}{\,M_{\odot}}\lesssim M_{*}\lesssim 10^{13}{\,M_{\odot}}$) can also be reproduced with the same sets of the parameters in Xia & Yu (2019).

Through cosmological zoom-in numerical simulations for single central galaxies, Ma et al. (2016) show that a large fraction of metals ejected by stellar feedback can be retained inside a halo (as illustrated by Figure 11 therein) even for low-mass galaxies down to $M_{*}\sim 10^{6}{\,M_{\odot}}$ and the retained metals can further rain down back to galaxy disks), which provides an alternative way to meet with both the observational stellar mass function and the MZR of dwarf galaxies. This scenario has a similar effect as the application of Equation (5) or (6) in the sense of retaining metals inside halos. Satellites are not included in the work of Ma et al. (2016).

• •

Metallicity enrichment: The metals ejected by SNe are assumed to first instantaneously and homogeneously mixed with the ISM in the galaxy, and after the mixture, some metals can be ejected out of the galaxy along with the mixed ISM that is ejected out by SN explosions. The ejected gas is assumed to have the same metallicity as the cold gas. In this work, the Fe yield of SNe II is adopted from tables 2–3 in Nomoto et al. (2006), and the Fe yield of SNe Ia is from Iwamoto et al. (1999). As done in HYL14, we assume that SN II explosions are instantaneous after the formation of SN II progenitors, and SN Ia explosions have a nonnegligible time delay after the formation of SN Ia progenitors. The SN Ia event rate adopted in HYL14 is based on the observational results by Maoz et al. (2010), in which a minimum time delay of 0.1 Gyr since the birth of a stellar population is reported.

• •

The reionization of the universe: the reionization in the early universe reduces the baryon fraction of a DM halo to be below the cosmic average (e.g., Gnedin 2000). The reduced fraction can be modeled through a mass scale called the filtering mass. In this work (also in HYL14), we adopt the recipes by Kravtsov et al. (2004) (see also Gnedin 2000; Okamoto et al. 2008) to calculate the filtering mass and model the effects of the reionization.

The filtering mass is a function of redshift, as well as a function of the starting redshift $z_{0}$ (when the first HII bubble formed) and the completing redshift $z_{r}$ of the reionization history. In this work, we set $z_{0}=15$ and $z_{r}=10$, as done in the fiducial model in HYL14. This relatively earlier reionization of the universe produces a stronger effect in reducing the baryonic fraction in low-mass DM halos, resulting in a relatively shallower slope or a smaller exponent in the MZR at the low-mass end (which agrees with the observation better, according to HYL14). Throughout the paper, the word “slope” and the word “exponent” of a power-law relationship are equivalently used. Lu et al. (2017) also discuss the importance of preventive feedback in meeting with both the observational stellar mass function and the MZR of dwarf galaxies, which requires a strong prevention of baryons from collapsing into low-mass halos in the first place and has a similar effect as this early reionization scenario.

The parameters used in the recipes are chosen from the fiducial model in HYL14.

In Xia & Yu (2019), we use the above SAM to investigate the universality of the MZRs for both satellites and central galaxies in a large number of different host DM halos covering a large mass range of $10^{9}$–$10^{15}h^{-1}{\,M_{\odot}}$. We find that the satellites in those simulated host halos follow a similar relation, with the exponent $\alpha$ being in the same range $\sim$0.2–0.4. The simulated central host galaxies with $10^{3}{\,M_{\odot}}\lesssim M_{*}\lesssim 10^{11}{\,M_{\odot}}$ follow a similar relation. That study further shows that a double power law provides a better fit to the MZR than the above single power law for both satellites and central galaxies, which gives $\alpha\sim$0.2–0.4 at $10^{3}{\,M_{\odot}}\lesssim M_{*}\lesssim 10^{8}{\,M_{\odot}}$ and a relatively higher $\alpha\sim 0.5$ at $10^{8}{\,M_{\odot}}\lesssim M_{*}\lesssim 10^{11}{\,M_{\odot}}$. The stellar metallicity in massive galaxies with $M_{*}\gtrsim 10^{11}{\,M_{\odot}}$ becomes roughly constant, close to the metal yield. The difference in the best-fit normalizations of the MZRs is within a small factor (e.g., $\sim 2$ at $M_{*}=10^{6}{\,M_{\odot}}$ and $\sim 4$ at $M_{*}=10^{3}{\,M_{\odot}}$). The relatively narrow space of the fit slopes and normalizations revealed a universality in the MZRs, which suggests the common physical processes in the stellar formation and chemical evolution of the galaxies can be unified with a large range of galaxy masses and halo masses.

3. Analytic chemical evolution Models

To investigate the metallicity evolution in a galaxy, we start with the following conservation laws of the total stellar mass, cold gas mass, and metal mass within the galaxy (see eq. 4.6-4.11 in Cole et al. 2000):

[TABLE]

where $M_{*}$, $M_{\rm cold}$, $M_{*}^{Z}$, and $M_{\text{cold}}^{Z}$ are the total stellar mass, the total cold gas mass, the total metal mass in stars, and the total metal mass in cold gas, respectively, the overdot $\dot{}$ above a variable represents its derivative with respect to time $d/dt$, i.e., the changing rates of the corresponding variable, $\dot{M}_{\text{cool}}$ represents the gas mass cooling rate from the hot gas in the DM halo to the cold gas in the galaxy disk, $\psi$ is the instantaneous star formation rate, $R$ is the fraction of mass recycled by stars to the ISM (through stellar winds and SNe), $Z_{\text{cold}}\equiv M_{\text{cold}}^{Z}/M_{\text{cold}}$ is the metallicity of the cold gas, $Z_{\text{hot}}$ is the metallicity of the hot gas in the DM halo, $\dot{M}_{\text{cool}}Z_{\text{hot}}$ represents the cooling rate of the gas in metals, $p$ denotes the yield (the fraction of mass converted into stars that is returned to the ISM in the form of metals), and $\beta$ is the efficiency of stellar feedback.

The stellar mass $M_{*}$ is a monotonically increasing function with time $t$, and thus a physical variable (e.g., metallicity) as a function of time can also be expressed as a function of $M_{*}$. One advantage of using $M_{*}$ as the independent variable instead of $t$ is that some complexity in the difference of star formation histories for different galaxies (which can be expressed through the function $M_{*}(t)$) can be partly removed in understanding the origin of the MZR. Another advantage of using $M_{*}$ as the independent variable is that mass is a physical variable of a system relating to its gravitational potential and some variables (e.g., $\beta$) is constrained by the gravitational potential of the system.

The mean metallicity in stars $Z_{*}\equiv M_{*}^{Z}/M_{*}$ as a function of stellar mass can be obtained by a combination of Equations (7) and (9) and an integration over the star formation history of the galaxy as follows,

[TABLE]

where $M_{*,0}$ and $Z_{*,0}$ are the initial stellar mass and the initial metallicity of the system, respectively.

A combination of Equations (7) and (8) gives

[TABLE]

and then applying Equation (14) into Equation (10) yields

[TABLE]

Below we show how $Z_{\rm cold}$ can be solved from the above equation, given some assumption about $Z_{\text{hot}}$ or its explicit form (e.g., related to $Z_{\rm cold}$ or $M_{*}$).

3.1. Some previous analytical models

Equation (15) is a basic equation to describe the chemical evolution in a galaxy, which can be reduced to many previous well-known metallicity enrichment models, e.g. the “closed-box” model (Schmidt 1963; Talbot & Arnett 1971; Searle & Sargent 1972), the “leaky-box” model, the “pre-enriched” model (Pagel 1997), and the “accretion” model (Lynden-Bell 1975). In most of those models, a constant (or zero) hot halo gas metallicity $z_{\rm hot}$ is assumed.

For the simple case of $Z_{\text{hot}}=0$, where the hot gas in the DM halo is metal-free and its cooling and accretion onto the galactic disk is a metal dilution process in the cold gas, Equation (15) is reduced to

[TABLE]

As mentioned above, here a physical variable as a function of time can be expressed as a function of $M_{*}$, and thus the above differential equation has an analytic solution as follows,

[TABLE]

where $M_{\text{cold},0}$ is the initial cold gas mass, $Z_{\text{cold},0}$ is the initial cold gas metallicity, $c$ is an arbitrary constant and can be canceled out in the calculation, $M_{*}^{{}^{\prime}}$ and $M_{*}^{{}^{\prime\prime}}$ are two dummy variables of stellar masses, and $Z_{\text{cold}}=Z_{\text{cold},0}$ when $M_{*}=M_{*,0}$ and $M_{\text{cold}}=M_{\text{cold},0}$. Below we show that the solutions in some previous models can be obtained by reducing Equations (16) and (17) with an assumption that $p$ and $R$ are constant.

• •

Leaky-box model. In this case, we set $\dot{M}_{\text{cool}}=0$ and the feedback efficiency $\beta$ to be constant (which becomes the closed-box model below if $\beta=0$), and thus Equation (14) becomes

[TABLE]

By applying Equation (18) into Equation (17), we obtain the solution of the cold gas metallicity as follows,

[TABLE]

and an application of Equation (19) into Equation (13) yields the stellar metallicity as follows,

[TABLE]

When $M_{\text{cold}}\ll M_{\text{cold,0}}$, we have $M_{\text{cold,0}}\simeq(M_{*}-M_{*,0})(1-R+\beta)/(1-R)$, and Equation (20) can be transformed to be:

[TABLE]

where the term on the left-hand side is the total metal mass in stars, the first term on the right-hand side is the initial metal mass in stars, the second term is the metal mass contributed from the initial cold gas to the newly formed stars (with mass $M_{*}-M_{*,0}$), and the third term is the effective metal mass due to SN nucleosynthesis in the newly formed stars (which is first released into the cold gas reservoir and then returned to the newly formed stars). Equation (21) can be transformed to be:

[TABLE]

Then we have $Z_{*}\rightarrow Z_{\text{cold,0}}+\frac{p}{1-R+\beta}$ if $M_{*}\gg M_{*,0}$, and we have $Z_{*}\rightarrow Z_{*,0}$ if $M_{*}\gg\left(\frac{1-R}{1-R+\beta}\right)M_{\rm cold,0}$.

In a traditional leaky-box model with $Z_{\text{cold,0}}=Z_{*,0}=0$ and $M_{*,0}=0$, where the galaxy evolves from a pure gas cloud with a total baryonic mass $M_{0}=M_{\text{cold,0}}$, Equations (19) and (20) are reduced to

[TABLE]

[TABLE]

where $f_{\text{gas}}\equiv M_{\text{cold}}/M_{\text{0}}$ is the fraction of the gas mass relative to the initial total baryonic mass. When $f_{\rm gas}\ll 1$, the stellar metallicity shown by Equation (24) can be reduced to the following form:

[TABLE]

which is usually called the effective yield (Cole et al., 2000). If $\left(\frac{1-R+\beta}{1-R}\right)\frac{M_{*}}{M_{\rm cold}}\ll 1$, we have $Z_{\rm cold}=\frac{p}{1-R}\frac{M_{*}}{M_{\rm cold}}$ and $Z_{*}=\frac{1}{2}\frac{p}{1-R}\frac{M_{*}}{M_{\rm cold}}\propto\frac{M_{*}}{M_{\rm cold}}$ (cf., Dekel & Silk 1986).

• •

Closed-box model. If the feedback efficiency $\beta=0$ is set in the above traditional leaky-box model, Equations (23) and (24) are reduced to

[TABLE]

and

[TABLE]

• •

Accretion model. In Lynden-Bell (1975), gas leaking is not considered (i.e., $\beta=0$), but a metal-free accretion is assumed, together with a specified relation between the cold gas mass and the stellar mass (so that the mass accretion rate is implicitly given). The specified relation between the cold gas mass and the stellar mass in Lynden-Bell (1975) is shown below (see eqs. 4.4 and 4.5 therein):

[TABLE]

where $M_{0}$ and $M_{\infty}$ are the initial and final baryonic masses of the galaxy, respectively. In the above relation, when a galaxy just formed with $M_{*,0}=0$, all the baryonic mass is in the cold gas with $M_{\text{cold,0}}=M_{0}$ and $Z_{\text{cold,0}}=0$; and when star formation approaches the end, we have $M_{\text{cold}}=0$ and $M_{*}=M_{\infty}$. In this case, Equation (17) can be reduced to

[TABLE]

By applying Equation (28) into Equations (30) and (13), the cold gas metallicity (see also eq. 4.8 in Lynden-Bell 1975) and the stellar metallicity can be obtained as follows,

[TABLE]

When $M_{*}\rightarrow M_{\infty}$, we have $Z_{*}\rightarrow\frac{p}{1-R}$.

• •

Pre-enriched model. In Pagel (1997), the hot halo gas is assumed to have a constant non-zero metallicity $Z_{\text{hot}}=Z_{\rm hot,0}>0$ (pre-enriched) before its cooling onto the galaxy disk. Equation (15) becomes

[TABLE]

which has a similar form as Equation (16) if $Z_{\text{cold}}-Z_{\rm hot,0}$ is taken as one variable. The solution of $Z_{\rm cold}$ in Equation (33) can be obtained by replacing $Z_{\rm cold}$ and $Z_{\rm cold,0}$ with $Z_{\rm cold}-Z_{\rm hot,0}$ and $Z_{\rm cold,0}-Z_{\rm hot,0}$, respectively, in Equation (17).

• •

Equilibrium models including gas inflows and outflows. In the metallicity evolution model by Finlator & Davé (2008), it is assumed that $dM_{\rm cold}=0$, $dZ_{\rm cold}=0$, and $Z_{\text{hot}}=\xi Z_{\text{cold}}$ (where $\xi$ is a parameter representing the ratio of $Z_{\text{hot}}$ to $Z_{\text{cold}}$), thus Equation (15) is reduced to $Z_{\rm cold}=\frac{p}{(1-R+\beta)(1-\xi)}$.

In the metallicity evolution model by Lilly et al. (2013), it is assumed that $dM_{\rm cold}\neq 0$, $dZ_{\rm cold}=0$, and $Z_{\rm hot}=Z_{\rm hot,0}$ (a constant) to obtain an equilibrium metallicity (see eq. 25 therein). Given those assumptions, Equation (15) can be transformed to be

[TABLE]

Equation (34) is consistent with equation (25) or (26) in Lilly et al. (2013). In the model of Lilly et al. (2013), the mass ratio of cold gas to stars is further assumed not to change, i.e., $d(M_{\rm cold}/M_{*})=0$. The assumption of $d(M_{\rm cold}/M_{*})=0$ is equivalent to assuming $dM_{\rm cold}=0$ or $dM_{*}=0$ if $M_{\rm cold}$ is a function of $M_{*}$ (as shown in Equation 35 in our model below).

3.2. Our model based on some power-law scaling relations with stellar

masses

In this subsection, we introduce a model by assuming that the cold gas mass $M_{\text{cold}}$ and the stellar feedback efficiency $\beta$ follow power-law scaling relations with the stellar mass $M_{*}$ during the growth of a galaxy, where the effects of gas inflows and outflows are included and there are no specific assumptions of $dM_{\rm cold}=0$ or $dZ_{\rm cold}=0$ or $d(M_{\rm cold}/M_{*})=0$. In reality, the cold gas mass in a galaxy and the feedback efficiency are not necessarily monotonic functions of $M_{*}$ during the galaxy evolution. As to be seen in Section 4.1 below, both the cold gas mass $M_{\text{cold}}$ and the feedback efficiency $\beta$ are taken as “median values” in our model, and we will show that this power-law scaling assumption is plausible and present the quantitative relations obtained from the SAMs. In this model, we also relax the assumption of a constant $Z_{\text{hot}}$, as metals transfer among hot gas, cold gas, and stars through accretion and ejection.

We assume that $M_{\text{cold}}$ and $\beta$ have the following power-law scaling relations with $M_{*}$,

[TABLE]

We assume $Z_{\text{hot}}=\xi Z_{\text{cold}}$ ($0\leq\xi\leq 1$), and it is plausible to have $Z_{\text{cold}}\geq Z_{\text{hot}}$, since the metals produced by star formation is assumed to be returned to the ISM, not ejected directly from the stellar disc to the hot gas (i.e., $e=0$ in eqs. 4.6-4.11 in Cole et al. 2000). Thus, Equation (15) can be transformed to be the following form:

[TABLE]

If $\xi$ is a constant, $Z_{\rm cold}$ in Equation (37) can be analytically solved with $Z_{\rm cold,0}=0$ as follows,

[TABLE]

where

[TABLE]

Below we discuss the following three simplified cases with different $\xi$ values, in each of which a different term in the right-hand side of Equation (37) is simplified to be zero.

• •

$\xi=0$: Equations (37) and (38) are reduced to Equations (16) and (17) (i.e., metal-free inflows with $Z_{\text{hot}}=0$). By adopting the power-law scaling relations of Equations (35) and (36) into Equation (17), if $\beta\gg 1-R$, we have

[TABLE]

where \alpha=1+\alpha_{\text{\betas}}-\alpha_{\text{gs}}, $A=\frac{1}{1-R}\cdot\frac{1}{\alpha}\cdot\frac{\beta M_{*}}{M_{\rm cold}}$, $\gamma(a,z)\equiv\int_{0}^{z}t^{a-1}e^{-t}dt$ ($\Re a>0$) is the lower incomplete gamma function, and $\Gamma(a,z)\equiv\int_{z}^{\infty}t^{a-1}e^{-t}dt$ is the upper incomplete gamma function. By applying Equations (35) and (36) to the $\alpha=0$ case in Equation (40), we have

[TABLE]

The stellar metallicity $Z_{*}$ can be obtained by a numerical integration over $Z_{\text{cold}}$ when $\alpha\neq 0$ (i.e., Eq. 13), and $Z_{*}=\frac{1}{2-\alpha_{\rm gs}}Z_{\rm cold}$ when $\alpha=0$ by analytical integration. The condition of $\beta\gg 1-R$ can be generally satisfied at $M_{*}\lesssim 10^{8}{\,M_{\odot}}$, as to be seen from Figure 8 below.

If $\beta\ll 1-R$, the solution of $Z_{\rm cold}$ can be obtained by simply setting $\beta=1-R$ and \alpha_{\text{\betas}}=0 in Equation (40) and in the expression of $A$ above, as $\beta$ is involved in Equation (17) or (37) only through the term $\frac{1-R+\beta}{1-R}$.

Although $Z_{\text{hot}}=0$ is assumed in this case as done in many previous models, both gas inflow and outflow with $\beta\neq 0$ are allowed. Since in this case the term of the feedback efficiency $\beta$ is used in Equations (37)-(38) and shown explicitly in the solution of $Z_{\rm cold}$ (Eq. 40), we call this case the $\beta$-model.

• •

$\xi=1$: in this case, the gas inflow metallicity (hot gas metallicity) is the same as the gas outflow metallicity (cold gas metallicity). To reach this scenario, the dynamic timescales of the inflow and outflow should be both much shorter than the gas accretion timescale onto the DM halo. With $f\left(M_{*}\right)=0$, Equation (38) is reduced to

[TABLE]

By adopting the power-law scaling relations, we have

[TABLE]

Note that in this case, the relation between the gas mass and the stellar mass is explicitly used in Equations (37)-(38) (see the term with $dM_{\rm cold}/dM_{*}$) and shown in the solution of $Z_{\rm cold}$ (Eq. 42), we call this case the g-model.

• •

Specifically, if $\xi$ is set so that the second $dM_{*}$ term on the right-hand side of Equation (37) vanishes, i.e.,

[TABLE]

which means that the inflow metal mass balances the outflow metal mass, Equation (37) is reduced to the following simple form:

[TABLE]

which has the solution as follows,

[TABLE]

This model is an intermediate case between the $\beta$-model and the g-model, as we have $0<\xi=(1+\frac{1-R}{1-R+\beta}\frac{dM_{\rm cold}}{dM_{*}})^{-1}<1$ according to Equation (45). Note that in this case only the term with $p$ is used in the right-hand side of Equation (37), we call this case the $p$-model.

Though the above equations on $Z_{\rm cold}$ and $Z_{*}$ obtained in this subsection have some similarity (e.g., proportional to $M_{*}/M_{\rm cold}$) as some previous model results as reviewed in Section 3.1, they are obtained by relaxing some specific assumptions used in those previous models, and the ratio of $M_{*}/M_{\rm cold}$ is a result of the combining effects of star formation, gas inflows, and outflows.

Note that the expectation from the g-model and the $p$-models should not be taken as the reference models at $M_{*}\gtrsim M_{\rm cold}$, as the hot gas metallicity of the inflow assumed in the models would be unrealistically too high so that $Z_{*}$ expected by extrapolating Equations (44) and (48) to $M_{*}\gtrsim M_{\rm cold}$ can be even higher than $\sim p/(1-R)$.

As seen from the results obtained in this subsection, the exponents in the MZRs are closely related to the exponents in the power-law scaling relations of the $M_{\rm cold}$–$M_{*}$ relation and the $\beta$–$M_{*}$ relation. The metallicities are solely determined by the $M_{\rm cold}$–$M_{*}$ relation in the g-model and the $p$-model. According to Equations (44) and (48), the slopes of the MZRs obtained with the g- and $p$-models should be $1-\alpha_{\rm gs}$. According to the $\alpha=0$ case in Equation (41), the slopes of the MZRs obtained with the $\beta$-model should be $-\alpha_{\beta{\rm s}}=1-\alpha_{\rm gs}$.

In principle, Equations (7)–(10) are one part of the SAM, to describe the chemical evolution of a galaxy in a newly formed halo. Our analytical model presented above in Section 3.2 has the following differences from that in the SAM. In our model, we ignore the physical process of galaxy mergers in the evolution of the stellar metallicity and the stellar mass. The effects of the detailed physical processes, such as halo assembly and gas cooling, are incorporated into the evolution of gas mass with stellar mass in a statistical way, as shown in Section 4.1.1 below. The universality and the simplicity in the expression of the MZR itself suggests that it is plausible here to simplify the growth of a galaxy as the evolution of several relevant physical variables (e.g., stellar mass, gas mass, metal mass) and to use a statistical and analytical way to reveal the origin of the MZR. The combination of the chemical evolution model constructed in this work with the results obtained from the SAM provides an efficient way to isolate the effects of different physical processes and see the dominant reasons leading to the MZRs.

Note that the chemical enrichment is included in the SAM as described in Section 2, but the metallicities produced from the SAM are not used in the chemical evolution model described in this section, except that the SAM metallicity results can provide a support to the assumed range of $\xi$ ($0<\xi\lesssim 1$) in the above analysis, as well as a reference to the MZR to be shown in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond below.

4. Results

In this section, we show the power-law scaling relations between the cold gas mass and the stellar mass and between the feedback strength and the stellar mass obtained from the SAMs. Then by applying these scaling relations to the chemical evolution models presented in Section 3.2, we reproduce the MZR and explain its universality. In principle, if these scaling relations are available from observations, they can also be similarly used to the chemical evolution models to predict the MZR. The compatibility of the SAM results with observations (as shown in HYL14 and Xia & Yu 2019) suggests the reasonability to extract some other physical properties (i.e., the power-law scaling relationships here) from the SAM.

4.1. Power law scaling relations obtained from the SAM

In this subsection, the relations between the cold gas mass and the stellar mass and the relation between the feedback strength and the stellar mass during the evolution of the galaxies and satellites are extracted from the SAM presented in HYL14 and Xia & Yu (2019). In Xia & Yu (2019), we present some of our SAM simulation runnings of the SAMs for central galaxies and their satellites in DM halos ranging from $10^{9}$ to $10^{15}h^{-1}{\,M_{\odot}}$ (see table 1 therein). In this work, we add some runnings for high-mass halos for better statistics. The stellar masses of the galaxies and the satellites range from $10^{3}M_{\odot}$ to $10^{11}M_{\odot}$ at redshift $z=0$. The SAM simulations produce the evolution of the physical properties of the galaxies and their satellites.

4.1.1 The cold gas mass–stellar mass scaling relation

The gas fraction in a galaxy evolves with inflows, outflows, and star formation processes. In a closed-box model, the cold gas fraction $M_{\text{cold}}/(M_{\text{cold}}+M_{*})$ decreases monotonically as the stellar mass increases. In a realistic situation with gas inflow and outflow in the hierarchical galaxy formation and evolution, a host halo can gain new gas through accretion or merging with other halos, the newly gained gas can further cool down onto the galaxy disk, and the cold gas can also be ejected from the galaxy disk through feedback. Thus, the evolution of the cold gas fraction may not be monotonic, but fluctuate as the stellar mass increases.

Figure 1 shows the evolution of the cold gas mass as a function of stellar mass in the progenitor galaxies of the present-day satellites before their infall into the host halos, which are obtained from our SAM simulations. As mentioned above, in the assembling history of a halo, the progenitor galaxy of a present-day satellite is the host central galaxy of a small halo at an early time before it fell into a big halo. In Figure 1, each panel presents the results obtained from the progenitor galaxies of the present-day satellites in one host halo and the different panels have different host halo masses at $z=0$. Each curve represents the growth history of the progenitor of one satellite until infall, which we call a track. The number of the tracks in each panel is denoted by $N_{\rm track}$. The inset in each panel illustrates only one track.

Each track shown in Figure 1 is obtained by tracing the evolution of the progenitor galaxy back to high redshifts. Along the evolution, the progenitor may experience a galaxy merger, and we select the relatively massive one of the merging galaxies as the progenitor galaxy.

As seen from Figure 1, though fluctuating, the curves show an overall increasing tendency with increasing stellar masses at $M_{*}\lesssim 10^{10}{\,M_{\odot}}$. Below we obtain the “median” relation between the cold gas mass and the stellar mass by a linear fitting to the data. The physical properties of the progenitor galaxies are recorded at the same time interval ($\sim 10^{6}$ yr) in our simulations. Note that the data are not uniformly distributed in the logarithm of the stellar mass. To avoid some weight bias in the fitting, we divide the logarithm of the stellar mass uniformly into some bins starting from $\log(M_{*}/{\,M_{\odot}})=3$ with an interval of 0.5 dex. For each evolution track, we take the medians of their variables ($\log(M_{*}),\log(M_{\rm cold})$) in each bin; and then we show the medians of those median values of all the evolution tracks in each bin by a red square and show the half of the range between the 16th and 84th percentiles of those median values of all the evolution tracks as the error-bar of the red square in Figure 1. In that way, the different evolution tracks contribute the same weight to the red square in each bin. We use the least-squares method to fit the red squares as follows,

[TABLE]

where $\alpha_{\text{gs}}$ is the slope, and $b_{\text{gs}}$ is the intercept at $M_{*}=10^{3}M_{\odot}$ and related with $k_{\text{gs}}$ in Equation (35) by $k_{\text{gs}}=10^{b_{\text{gs}}-3\alpha_{\text{gs}}}$. The best-fit results are shown by the red solid line and the texts in each panel of Figure 1.

In order to have a statistically robust result, we stack multiple trees of the same host halo mass together, and perform a similar fitting to the $M_{\rm cold}$–$M_{*}$ relation in the progenitors of their present-day satellites. Each panel in Figure 2 shows the best-fit results to the progenitors of the present-day satellites in a stack of ten host halos with the same halo mass. The fitting method is similar as done for Figure 2. In each panel, the evolution curves of the cold gas mass as a function of stellar mass in all the progenitor galaxies in the ten host halos are used to obtain the medians of the variables ($\log(M_{*}),\log(M_{\rm cold})$) (open squares) and their error-bars. The evolution tendency of the open squares shows that $M_{\rm cold}$ increases with increasing $M_{*}$ when $M_{*}$ is low but declines at the high-$M_{*}$ end. The increase of $M_{\rm cold}$ with increasing $M_{*}$ is associated with the halo mass buildup of the satellite progenitors. The decline is associated with star formation, the exhaustion of cold gas, and the slow increase of the halo mass in the halo assembly history. We fit the median $M_{\rm cold}$–$M_{*}$ relation mainly in the increasing part, and the stellar mass ranges used in the fitting are listed in Table 1. The best-fit results are also summarized in Table 1. As seen from Figure 2 and Table 1, where the host halo masses cover a large range, the slopes $\alpha_{\text{gs}}$ and the intercepts $b_{\text{gs}}$ of the fit scaling relation are in a narrow space for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$, with $\alpha_{\text{gs}}\sim$0.61–0.79 (or $1-\alpha_{\text{gs}}\sim$0.2–0.4) and $b_{\text{gs}}\sim$5.5–6.0.

The $\chi^{2}$ values shown in Table 1 are used as a measure of the goodness of the fit, together with the number of the bins $N_{\rm bin}$ in the fit. In general the fits are acceptable, and the probability that a random set of $N_{\rm bin}$ data points drawn from the parent distribution would yield a value of $\chi^{2}$ as large as or larger than the tabulated values are mostly in the range $>$20%–95%. Some $\chi^{2}$ values are somewhat too small, compared with the degrees of freedom in the fitting ($N_{\rm bin}-2$). Part of that reason is that a relatively large error has been assigned to the fitting data (i.e., using the large scatter around the median in each bin).

To compare the satellites with the central galaxies, we also perform a similar fitting to the median $M_{\rm cold}$–$M_{*}$ relation in the central galaxies of the stacked trees. Figure 3 shows the evolution tracks of the central galaxies in the stacked host halos. The best-fit results are shown by the straight solid lines and summarized in Table 1. As seen from Figure 3 and Table 1, though the host halo masses cover a wide range, the slopes $\alpha_{\text{gs}}$ and the intercepts $b_{\text{gs}}$ of the best-fit scaling relations are in a universal range for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$, with $\alpha_{\text{gs}}\sim$0.57–0.69 (or $1-\alpha_{\text{gs}}\sim$0.3–0.4) and $b_{\text{gs}}\sim$5.7–6.0.

For a comparison of the $M_{\rm cold}$–$M_{*}$ relations between central galaxies and satellite progenitors in different halo masses, we combine the points shown in Figures 2 and 3 together into Figure 4(a). This comparison demonstrates that the cold gas mass–stellar mass relations are universal for both central galaxies and satellite galaxies, and for a wide range of host halo masses.

In Figure 4(b), we show the observational results for galaxies at redshift $z=0$ as a reference to the $M_{\rm cold}$–$M_{*}$ relations obtained during the galaxy growth in this work. The results are represented by the ratios of the cold gas mass to the stellar mass. where the blue solid line represents the results obtained from our SAMs, and the other points and lines represent the observational results for galaxies at redshift $z=0$. Note that the blue solid line obtained from our SAM simulations is a little higher than the observational results at relatively low masses ($M_{*}\lesssim 10^{8}{\,M_{\odot}}$), which is plausible as the low stellar mass range of the blue solid line normally represents the relatively early evolution stage at high redshifts.

4.1.2 The stellar feedback efficiency–stellar mass scaling relation

In the SAM, the stellar feedback efficiency $\beta$ is obtained through Equation (4) or (6). In this subsection, we obtain the best-fit results to the scaling relations between the feedback strength $\beta$ and the stellar mass, similarly as done to obtain the scaling relations between the cold gas mass and the stellar mass in Section 4.1.1.

Figures 5–8 show the $\beta$–$M_{*}$ relations for the satellite progenitors and central galaxies in different halo masses, corresponding to the same simulation examples shown for the $M_{\rm cold}$–$M_{*}$ relations in Figures 1–4(a), respectively. As seen from the figures, $\beta$ has a decreasing tendency with increasing stellar masses. We obtain the power-law scaling relation by fitting the data as follows,

[TABLE]

where the stellar mass ranges used in the fitting are the same as done for the $M_{\rm cold}$–$M_{*}$ relation, \alpha_{\text{\betas}} is the slope, and b_{\text{\betas}} is the intercept at $M_{*}=10^{3}M_{\odot}$ and related with k_{\text{\betas}} in Equation (36) by k_{\text{\betas}}=10^{b_{\text{\betas}}-3\alpha_{\text{\betas}}}. The fitting method for the $\beta$-$M_{*}$ relation is similar as done for the $M_{\rm cold}$–$M_{*}$ relations above. The best-fit parameters of \alpha_{\text{\betas}} and b_{\text{\betas}} in Figures 6–7 are listed in Table 2. Note that $\beta\gtrsim 1$ at the low-$M_{*}$ range and $\beta\lesssim 1$ at the high-$M_{*}$ range. An extrapolation of the fitting formula to higher stellar masses has little influence on the chemical evolution model results, as $\beta$ is involved in the chemical evolution model through the term $1-R+\beta$ (Eqs. 8 and 10), which is approximately $1-R$ for $\beta\ll 1$.

Figure 8 shows that the $\beta$–$M_{*}$ relation is roughly universal for satellite progenitors. Table 2 shows that the slopes \alpha_{\text{\betas}} and the intercepts b_{\text{\betas}} of the best-fit scaling relations are in a universal range, with \alpha_{\text{\betas}} being in the range from -0.31 to -0.23 and b_{\text{\betas}}\sim2.4–2.7 for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$.

The median $\beta$–$M_{*}$ relation for central galaxies in halos with $M_{\rm halo}\sim 5\times 10^{10}$–$2\times 10^{12}h^{-1}{\,M_{\odot}}$ is roughly in the same narrow space, with \alpha_{\text{\betas}} being in the range from -0.40 to -0.23 and b_{\text{\betas}}\sim2.5–2.8. In halos with higher masses ($M_{\rm halo}\sim 10^{13}$–$10^{15}h^{-1}{\,M_{\odot}}$), Table 2 shows \alpha_{\text{\betas}}\sim-0.45 to -0.58 and $b_{\beta{\rm s}}\sim$2.9–3.2.

For a comparison of the $\beta$–$M_{*}$ relations between central galaxies and satellite progenitors in different halo masses, we combine the points shown in Figures 6 and 7 together into Figure 8. This comparison demonstrates that the $\beta$–$M_{*}$ relations are roughly universal for central galaxies in $M_{\rm halo}\lesssim 10^{13}h^{-1}{\,M_{\odot}}$ and satellite galaxies in a wide range of host halo masses with $10^{10}h^{-1}{\,M_{\odot}}\lesssim M_{\rm halo}\lesssim 10^{15}h^{-1}{\,M_{\odot}}$. For high present-day halo masses ($M_{\rm halo}\gtrsim 10^{13}h^{-1}{\,M_{\odot}}$), the stellar feedback efficiencies in central galaxies are relatively lower than those in the satellite progenitors at the stellar mass range $M_{*}\gtrsim 10^{6}{\,M_{\odot}}$ (which is because the central galaxies in that stellar mass range are at higher redshifts and have more compact sizes, and their disk rotation velocities are higher.) Note that the central galaxies in $M_{\rm halo}\gtrsim 10^{13}h^{-1}{\,M_{\odot}}$ are shown to have stellar mass $M_{*}\gtrsim 10^{11}{\,M_{\odot}}$ and their stellar metallicities are close to the yield in Xia & Yu (2019), while in this paper we focus on the stellar mass range of $10^{3}{\,M_{\odot}}\lesssim M_{*}\lesssim 10^{11}{\,M_{\odot}}$ of the MZR. The stellar feedback efficiency–stellar mass relation obtained for the central galaxy in $M_{\rm halo}\gtrsim 10^{13}h^{-1}{\,M_{\odot}}$ needs not to be used to explain their present-day stellar metallicity. Thus the deviation of the stellar feedback efficiency of the central galaxies in $M_{\rm halo}\gtrsim 10^{13}h^{-1}{\,M_{\odot}}$ shown in Figure 8 will not affect our results and conclusions in this work below.

The $\beta$–$M_{*}$ power-law scaling relation revealed in this subsection implies the existence of the $V_{\rm vir}$–$M_{*}$ and the $M{{}_{\rm halo}}$–$M_{*}$ power-law scaling relations, as $\beta$ is calculated from Equation (6) with $\beta\propto V_{\rm vir}^{-2}$ for a wide range of $M_{*}$ or $V_{\rm vir}$; for example, $\beta$ is taken to be $\beta{{}_{\rm E}}$ at $V_{\rm vir}\lesssim 200{\rm\,km\,s^{-1}}$ in Figure 2 in HYL14. As to be mentioned in Section 5, testing the origin of the $\beta$–$M_{*}$ relation (and the $M_{\rm cold}$–$M_{*}$ relation) is one of the next steps of the work.

4.2. The MZRs obtained from our chemical evolution model

By applying the above best-fit $M_{\rm cold}$–$M_{*}$ and $\beta$–$M_{*}$ power-law scaling relations (listed in Tables 1–2) to the analytic chemical evolution model constructed in Section 3.2, we obtain the expected MZRs, which is done in details by substituting the scaling relations of Equations (49) and (50) into Equations (38) and (13) for the $\beta$-, g-, and $p$-models in Section 3.2. For illustration, we show one example in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond, by adopting the best-fit parameters ($\alpha_{\rm gs}$, $b_{\rm gs}$, \alpha_{\text{\betas}}, b_{\text{\betas}}) obtained for the satellite progenitors with $M_{\rm halo}=2\times 10^{12}h^{-1}{\,M_{\odot}}$. As the $M_{\rm cold}$–$M_{*}$ and $\beta$–$M_{*}$ power-law scaling relations are shown to be roughly universal, the results for other cases do not differ much. We show the results of the $\beta$-, g-, and $p$-models with different line types in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond. Note that the MZR results of the $\beta$-model are also obtained numerically from Equations (38) and (13), not from the approximations of $\beta\gg 1-R$ and $\beta\ll 1-R$ analyzed in Section 3.2, so that the presented results cover $\beta\sim 1-R$ continuously. To illustrate the effectiveness of the chemical evolution model, the MZR results obtained from both the SAM in Xia & Yu (2019) and observations are also shown in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond(d) and (e) almost all of which fall in the gap between the $g-$ and $\beta$-model results in the range of $10^{3}M_{\odot}\leq M_{*}\leq 10^{10}M_{\odot}$ as shown in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond(f). Note that the g-model and the $\beta$-model represent the upper and lower bounds on the assumption of the infalling hot gas metallicity, and the differences among the g-, $p$-, and $\beta$-model results are not larger than 1 dex, so the MZRs are located within a narrow space and appear universal as shown in Xia & Yu (2019).

Note that for the MZRs of satellites, the $M_{\rm cold}$–$M_{*}$ relation and the $\beta$–$M_{*}$ relation used to obtain the dotted lines shown in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond are obtained from the satellite progenitors (before their infall into a bigger halo), while the solid lines shown in Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond are shown for satellites obtained after the satellite progenitors fell into a bigger halo. Here we argue that the evolution of the MZRs in satellites after their infall is negligible. As mentioned in Section 2 and in HYL14, we apply the energy condition only to a galaxy before it becomes a satellite. We do not apply it to satellites, but assume that the reheated gas from satellites (with mass expected by Equations 3 and 4) is expelled into the big host halo, as the original halos of the satellites are largely tidally disrupted along their motion in the big host halo, and the tidal field induced by the big host halo also helps to keep those expelled materials out of the satellites. According to the stellar feedback model described in Section 2, after the infall, the effective stellar feedback efficiency changes from $\beta{{}_{\rm E}}$ shown by Equation (6) to $\beta$ shown by Equation (4), which increases significantly (mostly in the low-$M_{*}$ range; see also figure 1 in HYL14). The chemical evolution model of a satellite after its infall can be approximated by the leaky-box model described in Section 3. Thus the change of the stellar mass after the infall is not significant due to the significant increase of $\beta$ (see Eq. 18), and neither is the change of the stellar metallicity (see Eq. 21).

As mentioned above, the $\beta$-, g-, and $p$-models are based on some simple assumptions on the ratio of the hot gas metallicity $Z_{\rm hot}$ to the cold gas metallicity $Z_{\rm cold}$. In Figure The Origin of the Stellar Mass–Stellar Metallicity Relation in the Milky Way Satellites and beyond, we also demonstrate how the MZRs expected from the chemical evolution model with other different ratios $\xi$ are distributed between the $g$- and the $\beta$-model results. As seen from the distribution, the expected MZR when $\xi\lesssim 0.1$ is quite close to the expectation from the $\beta$-model. The realistic ratios of $\xi$ obtained for the satellite progenitors in our SAMs are mostly distributed in the range of 0.1–0.9 at their infall into bigger host halos.

As mentioned above in Section 3.2, the exponents in the MZRs are expected to be closely related to the exponents in the power-law scaling relations of the $M_{\rm cold}$–$M_{*}$ relation and the $\beta$-$M_{*}$ relation, $1-\alpha_{\rm gs}$ or $-\alpha_{\beta{\rm s}}$, which is supported by the SAM results and the exponents shown in Tables 1 and 2, as summarized below.

• •

The slopes of the MZRs obtained from the SAM in Xia & Yu (2019) are in the range of $\sim 0.2$–0.4, as mentioned at the end of Section 2.

• •

For the $\beta$-model, the slope of the MZR expected by Equation (41) with $\alpha=0$ equals to $-\alpha_{\beta{\rm s}}$, which is shown to be in the range of 0.23 to 0.31 for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$ in Table 2, consistent with the SAM results. Note that the exponents in Tables 1 and 2 are generally satisfied with $\alpha=1-\alpha_{\rm gs}+\alpha_{\beta{\rm s}}\sim-0.14$ to 0.18 and $\frac{k_{\rm gs}}{k_{\beta{\rm s}}}\cdot 10^{-3(1-\alpha_{\rm gs}+\alpha_{\beta{\rm s}})}=10^{b_{{\rm gs}}-b_{\beta{\rm s}}-3}\sim$0.5–3 for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$, and correspondingly $\frac{1-R+\beta}{1-R}\cdot\frac{M_{*}}{M_{\rm cold}}\sim$0.5–5.

• •

For the $p$-model and the g-model, the slope of the MZR expected by Equations (44) and (48) equals to $1-\alpha_{\rm gs}$, which is shown to be in the range of 0.21–0.39 for $M_{\rm halo}\sim 5\times 10^{10}$–$10^{15}h^{-1}{\,M_{\odot}}$ in Table 1, consistent with the SAM result.

5. Summary and Discussion

In this work, we investigated the origin of the stellar mass–stellar metallicity relations. We constructed a chemical evolution model, based on the continuity equations in baryonic mass and metal mass conservations in the galaxy formation and evolution model. We applied some scaling relations between the cold gas mass and the stellar mass ($M_{\rm cold}$–$M_{*}$) and between the feedback efficiency and the stellar mass ($\beta$–$M_{*}$) to the model and obtained some analytical solutions in the MZRs. The solutions agree well with the simulation results obtained from the SAMs and the MZR observations. The range between the upper and the lower bounds in MZRs predicted in our models (through the upper and lower bounds in the ratio of the hot gas metallicity to the cold gas metallicity) is narrow, which provides an explanation to the universality in the MZRs revealed in the study by Xia & Yu (2019). The exponents in the MZRs are closely connected with the exponents in the $M_{\rm cold}$–$M_{*}$ and $\beta$–$M_{*}$ power-law scaling relations.

The $M_{\rm cold}$–$M_{*}$ and $\beta$–$M_{*}$ relations in this work represent the “average” evolutional tracks with the increase of the stellar masses in a galaxy, which appears to be universal for diverse star formation histories of different galaxies. Our work shows that the formation of the MZR is a result caused jointly by that the cold gas mass fraction decreases with increasing $M_{*}$ and that the stellar feedback efficiency decreases with increasing $M_{*}$ in the galaxy growth. The MZR represents an “average” evolutional track for the stellar metallicity of a galaxy.

In the chemical evolution model constructed in this work, both the effects of gas inflows and outflows are considered, and some specific assumptions used in the previous models are removed (e.g., those in the “closed-box” model, the “leaky-box” model, cold gas equilibrium, metal mass equilibrium as discussed in Section 3). The power-law scaling relations applied to the model are obtained from the SAM in HYL14 and Xia & Yu (2019), which give how the cold gas mass and the feedback efficiency change with the increase of the stellar mass during the growth of a galaxy. The effects of halo assemblies and the physical processes of gas cooling are implicitly included in the scaling relations extracted from the SAM, and they result in star formation and the growth of the stellar mass and the gas mass in a galaxy.

Our SAM simulation results show that the median relations between the cold gas mass and the stellar mass ($M_{\rm cold}$–$M_{*}$) and between the stellar feedback efficiency and the stellar mass ($\beta$–$M_{*}$) are distributed within a narrow space during the growth of $M_{*}$. During the growth of a galaxy, the median relations between the cold gas mass and the stellar mass appear to increase when the stellar mass is small and decline when the stellar mass is sufficiently large. The increasing part of the relation appears to be universal in a power-law scaling relationship, with $\log(M_{\rm cold}/M_{\odot})=\alpha_{\rm gs}\log(M_{*}/10^{3}{\,M_{\odot}})+b_{\rm gs}$ with $\alpha_{\rm gs}\sim$0.6–0.8 and $b_{\rm gs}\sim$5.5–6.0 over a large range of halo masses, during the growth of both central galaxies and satellite progenitors. The $\beta$–$M_{*}$ relation is also roughly the same for the progenitor satellites within $10^{10}h^{-1}{\,M_{\odot}}\lesssim M_{\rm halo}\lesssim 10^{15}h^{-1}{\,M_{\odot}}$ and central galaxies in $M_{\rm halo}\lesssim 2\times 10^{12}h^{-1}{\,M_{\odot}}$, following a power-law scaling relationship with $\beta=\alpha_{\beta{\rm s}}\log(M_{*}/10^{3}{\,M_{\odot}})+b_{\beta{\rm s}}$, $\alpha_{\beta{\rm s}}\sim-0.4$ to $-0.2$ and $b_{\beta{\rm s}}\sim$2.5–2.8. The exponents of the MZRs are close to the exponents of $1-\alpha_{\rm gs}$ or $-\alpha_{\beta{\rm s}}$.

The following reasons have also been proposed in the literature to explain the increase of the stellar metallicity with increasing stellar mass revealed in nearby galaxies by using (1) outflows, where $\beta$ decreases with increasing $M_{*}$; or (2) lower specific star formation rates in high-mass systems; or (3) changing of the yield $p$ due to different IMFs in different-mass systems. In this work, the evolution of the ratio $M_{*}/M_{\rm cold}$ and the formation of the MZR have the contribution from both outflows and inflows, where inflows are associated with halo growth. Point (2) can also be a consequence of a relatively high $M_{*}/M_{\rm cold}$ ratio in high-$M_{*}$ systems. The change of $p$ is not necessary for this work to match the observational MZR.

Some assumptions in some previous equilibrium models involving both inflows and outflows have been relaxed in this work. In some special cases, e.g., $Z_{\rm hot}=0$, there exists some similarity (but not exactly the same) between the form of the solution obtained in this work and the form obtained in those previous equilibrium models. If the expectations by some equilibrium models are compatible with observations, it may imply some kind of reasonability in the assumption that a local equilibrium at a given stellar mass is close to being reached in practice. The general solution obtained in our model can be reduced to a local equilibrium model at a given stellar mass, which requires that the average stellar mass is sufficiently massive compared to the average cold gas mass, in consideration of the effect of stellar feedback. However, in general, the model in this work involves a dynamical and secular evolution of the average stellar mass and the average gas mass during the growth of a galaxy, as well as the associated evolution of inflow gas metallicity, where the mass growth of a galaxy can span a vast mass range.

Xia & Yu (2019) show that a double power law exists in the MZR relations for both central galaxies and stacked satellites. In this paper, we do not investigate the details on non-linear effects of the MZRs revealed in Xia & Yu (2019), which can be related to the non-linear relations in the $M_{\rm cold}$–$M_{*}$ relations and in the $\beta$–$M_{*}$ relations.

Some next important steps of this work are to investigate the evolution of the MZR with redshift (e.g., Ma et al. 2016, where the evolution of the MZR with redshifts is attributed to the redshift evolution of the gas mass fraction within a halo, illustrated through cosmological zoom-in simulations and a “closed-box” model), connect it to cold gas phase metallicities revealed in observations (e.g., Tremonti et al. 2004; Mannucci et al. 2010; Lee et al. 2006) and the hot gas metallicities in DM halos, test the $M_{\rm cold}$–$M_{*}$ relation and the $\beta$–$M_{*}$ relation and their origin, investigate the non-linear effects in the correlations, and explore whether there exists a possible dependence on three or more parameters (e.g., star formation rate, as shown in Ellison et al. 2008; Mannucci et al. 2010; Barrera-Ballesteros et al. 2017) etc.

We thank Sandra Faber, Youjun Lu, Filippo Mannucci, and Ying-jie Peng for helpful discussions. This work was supported in part by the National Natural Science Foundation of China under Nos. 11673001, 11273004, 10973001, 11721303, the National Key R & D Program of China (Grant No. 2016YFA0400703), and the Strategic Priority Program of the Chinese Academy of Sciences (grant No. 23040100).