Is the stellar mass-stellar metallicity relation universal in the Milky Way satellites and beyond?

TL;DR

This study investigates whether the stellar mass-stellar metallicity relation (MZR) is universal across different galaxy types and halo masses, revealing consistent patterns with some variations influenced by halo mass and galaxy classification.

Contribution

It demonstrates that the MZR exhibits a double power law with variations depending on galaxy mass and halo environment, highlighting the universality of underlying physical processes.

Findings

MZR exponents range from 0.2 to 0.4 with scatter increasing in lower mass halos.

The MZR remains stable across different halo masses and galaxy classifications when averaged.

A double power law characterizes the MZR, with a steeper slope at higher stellar masses.

Abstract

Observations reveal a universal stellar mass-stellar metallicity relation (MZR) existing in Local Group dwarfs of different types, $Z_*\propto M_*^{\alpha}$ with $\alpha=0.30\pm0.02$. In this work, we investigate the "universality" of the MZRs for both satellites and central galaxies in a large number of different host dark matter halos covering a large mass range of $10^9$-$10^{15}h^{-1}M_\odot$, by using a semianalytical galaxy formation and evolution model. We obtain the following results. (1) The exponents ($\alpha$) for the MZRs of the satellites in halos with the same mass as the Milky Way halo but different individual assembly histories are mostly $\sim$0.2-0.4, i.e., having a scatter of $\sim 0.2$; and the scatter of $\alpha$ increases with decreasing halo masses. (2) The MZR relations are changed little by the variation of halo masses and the classification between central galaxies and satellites, if many halos with the same mass are stacked together. (3) A double power law exists in the MZR relations for both central galaxies and stacked satellites, with $\alpha\sim$0.2-0.4 at $10^3M_\odot< M_*<10^{8} M_\odot$ and a relatively higher $\alpha\sim0.5$ at $10^8M_\odot<M_*<10^{11}M_\odot$. (4) The high-mass satellites ($M_*>10^8M_\odot$) existing mostly in high-mass halos can lead to an apparent increase of $\alpha$ (from $\sim0.2$ to $\sim0.4$) with increasing host halo masses shown in the single power law fitting results of stacked satellites. The universality of the MZR suggests the common physical processes in stellar formation and chemical evolution of galaxies can be unified over a large range of galaxy masses and halo masses.

Tables

Table 1. Table 1 Simulation Runnings

$M{}_{\mathrm{halo}}$	$M{}_{\mathrm{res}}$	$N$	$M{}_{\mathrm{halo}}$	$M{}_{\mathrm{res}}$	$N$
($h^{-1}{M}_{\odot }$)	($h^{-1}{M}_{\odot }$)		($h^{-1}{M}_{\odot }$)	($h^{-1}{M}_{\odot }$)
${10}^9$	${10}^4$	10000
$5\times {10}^9$	${10}^4$	10000
${10}^{10}$	${10}^4$	10000
$5\times {10}^{10}$	${10}^4$	10000
${10}^{11}$	${10}^4$	10000
$5\times {10}^{11}$	${10}^5$	$1000$	$5\times {10}^{11}$	${10}^4$	$10$
${10}^{12}$	${10}^6$	$1000$	${10}^{12}$	${10}^5$	$10$
$2\times {10}^{12}$	${10}^6$	$1000$	$2\times {10}^{12}$	${10}^5$	$10$
$5\times {10}^{12}$	$5\times {10}^6$	$10$	$5\times {10}^{12}$	$5\times {10}^5$	$10$
${10}^{13}$	${10}^7$	$10$	${10}^{13}$	${10}^6$	$10$
${10}^{14}$	${10}^8$	$10$	${10}^{14}$	${10}^7$	$10$
${10}^{15}$	${10}^9$	$10$	${10}^{15}$	${10}^8$	$10$

Table 2. Table 2 The fitting results of the MZR

objects	$M{}_{\mathrm{halo}}$	$N=10000$ or 1000222$N=10000$ for $M{}_{\mathrm{halo}}{}^{}\le {10}^{11}{h}^{-1}{M}_{\odot }$and $N=1000$ for $M{}_{\mathrm{halo}}{}^{}=5\times {10}^{11}$– $2\times {10}^{12}{h}^{-1}{M}_{\odot }$, as listed in Table 1.				$N=10$				$N=10$
	($h^{-1}{M}_{\odot }$)					$M{}_{\mathrm{halo}}{}^{}/M{}_{\mathrm{res}}{}^{}\sim {10}^7$				$M{}_{\mathrm{halo}}{}^{}/M{}_{\mathrm{res}}{}^{}\sim {10}^6$
		$\alpha$	$b$	${\chi }^2$	$N_{\mathrm{bin}}$	$\alpha$	$b$	${\chi }^2$	$N_{\mathrm{bin}}$	$\alpha$	$b$	${\chi }^2$	$N_{\mathrm{bin}}$
	$5\times {10}^9$	$0.27\pm 0.19$	$-2.3\pm 0.2$	-	2
	${10}^{10}$	$0.24\pm 0.09$	$-2.3\pm 0.2$	$0.02$	$3$
	$5\times {10}^{10}$	$0.25\pm 0.04$	$-2.4\pm 0.1$	$0.8$	$5$
	${10}^{11}$	$0.28\pm 0.03$	$-2.5\pm 0.1$	$2.4$	$6$
satellites	$5\times {10}^{11}$	$0.31\pm 0.03$	$-2.6\pm 0.1$	$1.4$	$7$	$0.26\pm 0.13$	$-2.5\pm 0.3$	$0.09$	$3$	$0.31\pm 0.09$	$-2.5\pm 0.2$	$0.20$	$4$
(single power law)	${10}^{12}$	$0.35\pm 0.03$	$-2.7\pm 0.1$	3.2	8	$0.29\pm 0.03$	$-2.6\pm 0.1$	5.6	6	$0.27\pm 0.03$	$-2.6\pm 0.1$	2.3	6
	$2\times {10}^{12}$	$0.35\pm 0.03$	$-2.7\pm 0.1$	1.4	8	$0.30\pm 0.04$	$-2.6\pm 0.2$	0.23	6	$0.31\pm 0.02$	$-2.6\pm 0.1$	0.12	6
	${10}^{13}$					$0.39\pm 0.03$	$-2.9\pm 0.1$	4.2	8	$0.38\pm 0.03$	$-2.9\pm 0.2$	2.4	8
	${10}^{14}$					$0.39\pm 0.04$	$-3.0\pm 0.2$	1.4	8	$0.41\pm 0.03$	$-3.0\pm 0.2$	2.2	8
	${10}^{15}$					$0.41\pm 0.03$	$-3.1\pm 0.2$	2.0	8	$0.44\pm 0.03$	$-3.2\pm 0.1$	1.3	8
	$M_{\ast }/M_{\odot }$	$\alpha$	$b$	${\chi }^2$	$N_{\mathrm{bin}}$	$\alpha$	$b$	${\chi }^2$	$N_{\mathrm{bin}}$
central galaxies	${10}^3$–${10}^{11}$	$0.30\pm 0.01$	$-2.52\pm 0.05$	18.8	16
(single power law)
	${10}^3$–${10}^8$ (${10}^{13}$)333The values in the brackets are the host halo masses in unit of $h^{-1}{M}_{\odot }$.					$0.31\pm 0.05$	$-2.7\pm 0.2$
	${10}^8$–${10}^{11}$ (${10}^{13}$)					$0.54\pm 0.08$	$-3.8\pm 0.5$	0.4444The total ${\chi }^2$ in the mass range ${10}^3$–${10}^{11}{M}_{\odot }$.	8555The total bin number in the mass range ${10}^3$–${10}^{11}{M}_{\odot }$.
satellites	${10}^3$–${10}^8$ (${10}^{14}$)					$0.34\pm 0.06$	$-2.8\pm 0.2$
(double power law)	${10}^8$–${10}^{11}$ (${10}^{14}$)					$0.50\pm 0.09$	$-3.6\pm 0.6$	0.04	8
	${10}^3$–${10}^8$ (${10}^{15}$)					$0.35\pm 0.06$	$-2.9\pm 0.2$
	${10}^8$–${10}^{11}$ (${10}^{15}$)					$0.53\pm 0.09$	$-3.8\pm 0.5$	0.1	8
central galaxies	${10}^3$–${10}^8$	$0.24\pm 0.02$	$-2.38\pm 0.06$
(double power law)	${10}^8$–${10}^{11}$	$0.47\pm 0.02$	$-3.5\pm 0.1$
	$>{10}^{11}$	0	$0.26\pm 0.01$	8.3666The total ${\chi }^2$ over the whole mass range ($>{10}^3{M}_{\odot }$).	20777The total bin number over the whole mass range ($>{10}^3{M}_{\odot }$).
observation		$\alpha =0.30\pm 0.02$, $b=-2.59\pm 0.04$

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}}\ge 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}$$

$$[Fe/H]=\alpha \log (M_{\ast }/10^3{M}_{\odot })+b,$$

Full text

††affiliationtext: †yuqj@pku.edu.cn

1. Introduction

The dwarf satellites in the Milky Way (MW) are potentially powerful probes of the baryonic processes in galaxy formation occurring in the early universe. Observations reveal that there exists a universal stellar mass–stellar metallicity relation in the MW/M31 dwarf satellites and some other dwarf irregular galaxies in the Local Group, $Z_{*}\propto M_{*}^{\alpha}$ with $\alpha=0.30\pm 0.02$ (Kirby et al. 2013, K13 hereafter). The correlation between the stellar metallicity and the stellar mass was reproduced in the satellites of MW-like galaxies by using the semianalytical galaxy formation and evolution model (e.g., Font et al. 2011; Li et al. 2010; Lu et al. 2014; Hou et al. 2014, HYL14 hereafter). In this work, we investigate how universal the stellar mass–stellar metallicity relation (MZR) obtained from the semianalytical model is, that is, whether the MZRs obtained for the following different cases could fall onto the same relation: a) same halo masses but with different halo assembly histories; b) different halo masses; c) satellite galaxies or central galaxies. Here by “universal” we mean not only the “universal” existence of a correlation between the stellar metallicity and the stellar mass, but also the “universal” or roughly the same quantitative ranges for the slope ($\alpha$) and the normalization of the correlation. The exploration of whether this relation exists and how universal it is will help us to understand the common physical processes (e.g., star formation, chemical evolution) involved in the formation and evolution of different galaxies, to reveal the origin in shaping the relation, and to further provide constraints on galaxy formation and evolution models.

In this work, we employ dark matter halo merger trees and a semianalytical galaxy formation model (Cole et al. 2000; see also White & Frenk 1991; Kauffmann et al. 1993; Somerville & Primack 1999; Somerville et al. 2008) to generate galaxies and their satellites, as done in HYL14. The DM halo merger trees are generated from the modified extended Press-Schechter function (Press & Schechter, 1974; Bond et al., 1991; Lacey & Cole, 1993; Somerville et al., 2008) by the Monte Carlo method. The model provides an efficient way to explore the effects of different physical processes on observational properties. This exploration can be done by varying the model parameters of the physical processes and allows a statistical study by generating a large number of dark matter merger trees for each set of parameters. HYL14 uses this method to obtain the MZR for the dwarf satellites of MW-like central galaxies in MW-size dark matter halos, which is consistent with the observation of MW satellites; and the slope in the MZR is shown to be affected by the SN feedback strength and the reionization epoch. Note that in the semianalytical model not every central galaxy hosted by an MW-size dark matter halo is MW-like. Also note that the correlation between stellar metallicities and stellar masses exists not only in Local Group dwarfs, but also in local galaxies with stellar mass $10^{9}\lesssim M_{*}/M_{\odot}\lesssim 10^{12}$ as revealed by the SDSS spectra of over 40,000 galaxies (Gallazzi et al. 2005, 2006; K13; see also Panter et al. 2008; González Delgado et al. 2014). In this work, we generalize our study to non-MW-like galaxies and to a large range of host halo masses ($\sim 10^{9}{\,M_{\odot}}$–$10^{15}{\,M_{\odot}}$). We explore whether the MZR is universal in the central galaxies with a large range of stellar masses ($\sim 10^{3}{\,M_{\odot}}$–$10^{12}{\,M_{\odot}}$) and in their satellites.

The relation between the metallicities of galaxies and their masses resulting from semianalytical galaxy formation models has been explored extensively in the past (e.g. Yates et al. 2013; Somerville et al. 2015; Lu et al. 2017; Font et al. 2011; Guo et al. 2011; Li et al. 2010; Lu et al. 2014; HYL14). Our study is distinguished from previous works by a few aspects of the scope, the purpose, and the method, as follows.

• •

As mentioned above, the stellar mass and the halo mass of the MZR studied in this paper cover large ranges (with stellar masses $\sim 10^{3}$–$10^{12}{\,M_{\odot}}$ and halo masses $\sim 10^{9}$–$10^{15}{\,M_{\odot}}$), and the galaxies include both central galaxies and their satellites. The mass ranges covered in previous studies are smaller, and many do not put central galaxies and their satellites together in the study. For example, Yates et al. (2013) obtain the stellar mass–cold gas metallicity relation for central galaxies with relatively high stellar mass range $\sim 10^{9}$–$10^{11}{\,M_{\odot}}$. Somerville et al. (2015) present the stellar mass–cold gas metallicity relation for central galaxies with stellar mass range $\sim 10^{7}$–$10^{11}{\,M_{\odot}}$. Lu et al. (2017) present the MZR of MW dwarfs with stellar mass range $\sim 10^{3}$–$10^{9}{\,M_{\odot}}$.

• •

The previous studies do not investigate how the MZRs are different in a wide range of galaxy masses and halo masses, and they focus mainly on discussing different aspects or the effects/roles of different recipes/physical processes in semianalytical models. For example, Yates et al. (2013) discuss the effects on the element abundances obtained by updating the chemical evolution model and including delayed enrichment from stellar winds, SNe II and SNe Ia, etc., in local star forming galaxies, elliptical galaxies, and MW-like disk galaxies. Somerville et al. (2015) discuss the effects on the evolution of some fundamental galaxy properties (e.g., stellar mass functions, the relation between stellar mass and star formation rate, the relation between stellar mass and cold gas phase metallicity) obtained by partitioning cold gas into different phases and modeling the conversion of molecular gas into stars. Lu et al. (2017) discuss the importance of different types of feedback mechanisms (including both “preventive feedback” and “ejective feedback”) to explain the observational properties of MW satellites. Font et al. (2011) illustrate that chemical properties can be used to break the degeneracy in the effects of SN feedback and reionization on the luminosity function of MW dwarf satellites. As mentioned above, the physical properties obtained in those studies do not cover a wide range of galaxy masses and halo masses as done in this paper.

• •

We adopt the Monte Carlo method based on the modified extended Press-Schechter function to generate a large number (e.g., 100 or more trees) of halo merger trees, which provides an efficient way to perform a statistical study of the MZR. By contrast, the number of the merger trees obtained with sufficiently high resolutions from cosmological N-body simulations are quite limited, e.g., six MW-size halo assembly examples provided by the high-resolution Aquarius dark matter simulations(Springel et al. 2008; Starkenburg et al. 2013; see also the Via Lactea simulation, Diemand et al. 2007; Rocha et al. 2012).

Different scenarios have been proposed to explain the correlation between the metallicities of galaxies and their masses (e.g., K13; Brooks et al. 2007; Köppen et al. 2007; Dalcanton 2007; Lilly et al. 2013; Ma et al. 2016; Finlator & Davé 2008; Lu et al. 2017; Dekel & Silk 1986, and a review in Finlator 2016 and references therein). The study of semi-analytical models provides a way to reveal the common reasons that lead to a universal relation. In this paper, we mainly show the MZR obtained from semi-analytical models and illustrate the universality. In our next paper, with the tool of the semianalytical model and its results, we present our detailed chemical evolution model and the explanation for the universal MZR revealed in this paper.

This paper is organized as follows. In Section 2, we briefly describe the semi-analytical galaxy formation and evolution model used in this work. Previously we have shown that the MZR in the satellites of MW-like galaxies produced from the model 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, central galaxies, and dwarf satellites in Section 3. Remarkably, we find that the MZR is quite independent of changing these variations as long as the fiducial model parameters are fixed. A summary is given in Section 4.

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. Method

In this section, we briefly describe the semianalytical galaxy formation model that is used to explore the MZR of MW dwarf satellites and beyond in this work. The backbone of the model is the merger trees of DM halos, which may represent the hierarchical growth history of the host halo. Detailed semianalytical recipes for galaxy formation and evolution (for reference, see Cole et al. 2000; White & Frenk 1991; Kauffmann et al. 1993; Somerville & Primack 1999) are incorporated into the merger trees to obtain the observational properties of the central galaxy and its satellites.

We plant the halo merger trees using the Monte Carlo method developed by Parkinson et al. (2008; see also Kauffmann et al. 1993; Somerville & Kolatt 1999; Cole et al. 2000), which is based on a modified version of the extended Press-Schechter formula. The obtained halo mass functions are in good agreement with the $N$-body simulation results. The merger trees are built from redshift $z=0$ to 20, with 79 equal intervals in the logarithm of 1 + $z$. The minimum progenitor halo mass set in the halo merger trees is defined as the mass resolution $M{{}_{\rm res}}$, and the progenitors with mass lower than $M{{}_{\rm res}}$ cumulatively contribute to the accreted halo mass in the growth history of a halo (e.g., see eq. 3.5 in Cole et al. 2000). For each halo mass, we generate a number of merger trees (e.g., 10000, 1000, or 10) and apply the semi-analytical galaxy formation and evolution recipes to them. The related recipes, e.g., on SN feedback, reionization, and gas cooling, are summarized below in this section. More details can be found in HYL14 and references therein. In the assembly history of a halo, a satellite of the present-day central galaxy was the host galaxy of a small halo at an early time before it fell into a big halo.

As described in HYL14, the semianalytical galaxy formation model in this study is based on GALFORM (Cole et al., 2000; Benson et al., 2003; Bower et al., 2006), but with several modifications. As demonstrated by Cole et al. (2000), Kauffmann et al. (1993), Somerville & Primack (1999), Croton et al. (2006), and Bower et al. (2006), the semianalytical models can successfully reproduce a number of observations on the statistical distributions of galaxy properties, including the galaxy luminosity function, the stellar mass function, etc., and different observational constraints can provide probes of different underlying physical processes. Somerville et al. (2015) show that both the MZR and the stellar mass function at $z=0$ in the stellar mass range of $10^{7}$–$10^{11}{\,M_{\odot}}$ can be reproduced. Lu et al. (2017) discuss the status in reproducing the MZR simultaneously with other observational constraints in the stellar mass range of $10^{3}-10^{9}{\,M_{\odot}}$, and suggest that both ejective feedback (strong outflow) and preventive feedback (photoionization heating in low-mass halos) are necessary to reproduce the MZR and the stellar mass function. As shown by Figures 1 and 3 in HYL14, our model has reproduced both the luminosity function and the MZR in MW-like halos, with including the effects of both reionization and SN feedback. In this study, we focus on generalizing the same semianalytical model to see how the MZRs vary in a large range of different halo masses and galaxy masses.

Specifically, to study the properties of the MW satellites, MW-size halos with mass $\sim$1–2$\times 10^{12}h^{-1}{\,M_{\odot}}$ are set in the model and MW-like galaxies are selected from the generated galaxies, as done in HYL14. The selection criteria for the MW-like host galaxies is that the total stellar mass of a present-day host galaxy is in the range of $4-6\times 10^{10}{\,M_{\odot}}$ and its bulge mass to disk mass ratio is between 0.1 and 0.4 (e.g., Licquia & Newman 2015; McMillan 2011).

In this work, to address the universality of the MZR, we generalize the study to galaxies and DM halos with different masses, so that the properties of galaxies with different masses and their satellites in different environments can be explored.

Below are the related semi-analytical galaxy formation and evolution recipes.

• •

Gas cooling: In a newly formed dark matter halo, the initial total hot gas mass is a sum of the total hot gas mass from its progenitor halos (with mass higher than $M_{\rm res}$) and the accreted gas mass. In a halo merger tree, the accreted gas mass is the baryon mass in the accreted progenitor halos with mass lower than $M_{\rm res}$. The baryon fraction in the accreted halo mass is obtained by the cosmic average baryon fraction multiplied by a fraction reduced by the reionization of the universe. The fraction reduced by the reionization of the universe will be described below (or see eq. 1 in HYL14). The effect of reionization is also included in the estimate of the total hot gas mass that can be kept from its progenitor halos with mass higher than $M_{\rm res}$. The initial temperature of the hot gas in the newly formed halo is set to be the virial temperature of the halo.

The hot gas in the dark matter halo cools down to the halo center through atomic cooling with a rate depending on the temperature, mass density, and metallicity of the hot gas (Sutherland & Dopita, 1993). The cooling recipe of the hot halo gas in this work is adopted from HYL14, which is based on the GALFORM model and modified by including the mixing of the gas reheated by SN feedback with hot halo gas (mainly occurring in low-mass galaxies or progenitors of dwarf satellites).

As done in HYL14, the molecular hydrogen cooling process for pristine gas that occurred in the mini-halos in the early universe is modeled by using Equations (21)-(29) in Benson (2010; see also Galli & Palla 1998), and there is no molecular hydrogen cooling after the completeness of the reionization of the universe due to the suppression of the abundance of hydrogen molecules caused by the strong UV background (e.g., Wolcott-Green et al. 2017). HYL14 shows that the MZR of the satellites in MW-like galaxies is not sensitive to whether the molecular hydrogen cooling recipe is included or not.

• •

Star formation and stellar feedback: following Cole et al. (2000), stars are formed in the disk at a rate directly 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).

Star formation will re-heat the cold gas in the disk and possibly expel it out of the galaxy. The gas mass reheated by the feedback during time interval $dt$ 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 with $V{{}_{\rm hot}}=200{\rm\,km\,s^{-1}}$ and $\alpha{{}_{\rm hot}}=3.2$. Part or all of the reheated gas can escape the dark matter halo if the strength of the total feedback energy released and coupling to the intergalactic medium (IGM) is sufficiently large with

[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. For low-mass galaxies or the progenitors of satellites (e.g., $M_{*}\lesssim 10^{8}{\,M_{\odot}}$), normally the above inequality (5) cannot be satisfied, the reheated gas that is expelled out of the disk stays in the halo, with mass during time interval $dt$ given by $\beta{{}_{\rm E}}\psi dt$, where

[TABLE]

(see more details in section 2 in HYL14).

The feedback due to SN Ia explosions is included in this work. We assume that the energy released by an SN II explosion and an SN Ia explosion is the same. We use the same feedback recipe for each (Equations 3–6), but include the non-negligible time delay between the formation of SN Ia progenitors and the SN Ia explosions (See Equations 13-16 in HYL14).

As done in HYL14, in the assembly history of a halo, a satellite of the present-day central galaxy may be the host galaxy of a small isolated halo before it fell into a big halo at an early time. We apply the above 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 Equation 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. Note that 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. In Starkenburg et al. (2013), the tidal stripping and disruption of satellites was shown not to have a significant effect on the satellite luminosity functions, as it affects very few satellites.

• •

Metallicity enrichment: the metals ejected by SNe are assumed to be homogeneously and instantaneously mixed with the interstellar medium in the galaxy, and after the mixture, some metals can be ejected out of the galaxy along with the mixed interstellar medium that is ejected out by SN explosions. 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).

• •

The reionization of the universe: the reionization in the early universe reduces the baryon fraction of a DM halo. The extent of the reduction is modeled through a mass scale called the filtering mass, which is a function of redshift, as well as a function of the redshift that the first ionized bubble formed $z_{0}$, and the completion redshift of the reionization $z_{r}$. A halo with mass lower than the filtering mass loses more than 50% of the baryonic matter expected from the cosmic average. In this paper, we use Equations (B1) and (B2) in Kravtsov et al. (2004; see also Gnedin 2000; Okamoto et al. 2008) to calculate the filtering mass and model the effects of the reionization. In this work, we set $z_{0}=15$ and $z_{r}=10$ as done in the fiducial model in HYL14, which results in an MZR compatible with current observations of the MW satellites. HYL14 shows that an early reionization epoch results in a relatively high metallicity at the low-mass end and a relatively flat slope in the MZR. The relatively strong/early reionization epoch revealed in HYL14 suggests that the local universe is reionized earlier than the cosmic average; local sources may have a significant contribution to the reionization in the local region.

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

3. Results

In this section, we show our simulation results of the MZR obtained by using the method described in Section 2. Table 1 lists the simulation runnings done in this work, with the host halo mass $M{{}_{\rm halo}}$ ranging from $10^{9}h^{-1}{\,M_{\odot}}$ to $10^{15}h^{-1}{\,M_{\odot}}$ at $z=0$ and the number of the merger trees generated by the Monte Carlo method for each set of parameters. For low halo masses, a mass resolution $M{{}_{\rm res}}=10^{4}h^{-1}{\,M_{\odot}}$ is set. For high halo masses with $M{{}_{\rm halo}}\geq 5\times 10^{11}h^{-1}{\,M_{\odot}}$, although the values of $M{{}_{\rm res}}$ are set to be relatively high, we set two sets of $M{{}_{\rm res}}$ in the runnings, which shows that the simulation results are close to being convergent below.

3.1. The MZR of the satellites in MW-size halos

The MW halo mass is in the range $0.8-4.5\times 10^{12}{\,M_{\odot}}$ (e.g., Phelps et al. 2013; Kafle et al. 2014; Boylan-Kolchin et al. 2013; Callingham et al. 2018; Posti & Helmi 2018; Watkins et al. 2018, and references therein). As shown in Table 1, we use the Monte Carlo method to generate 1000 halo merger trees for both $M{{}_{\rm halo}}=10^{12}h^{-1}{\,M_{\odot}}$ and $M{{}_{\rm halo}}=2\times 10^{12}h^{-1}{\,M_{\odot}}$ at $z=0$ (with $M{{}_{\rm res}}=10^{6}h^{-1}{\,M_{\odot}}$) and incorporate the baryonic processes described in Section 2 into the halo merger trees to obtain the MZRs of the central galaxies and their satellites in the halos. An example of the MZR of the satellites in the MW-size halos obtained from the simulations is illustrated by the dots in Figure 2(b)-(c). The similar results are also shown in figure 3 in HYL14.

Figure 1 shows the fitting results for the MZR of the satellites in the 1000 host halos with mass $M{{}_{\rm halo}}=10^{12}h^{-1}{\,M_{\odot}}$ and $2\times 10^{12}h^{-1}{\,M_{\odot}}$. Note that the satellites are not uniformly distributed in the logarithm of the stellar mass, and there are more low-mass systems than high-mass ones. To avoid that low-mass satellites contributing too much weight in the fitting (which is more significant in high-mass halos), the fitting is performed in the following way. For each halo, the logarithm of the stellar mass of the satellites of the host galaxy at $z=0$ is divided into some bins starting from $\log(M_{*}/{\,M_{\odot}})=3$ with an interval of 1 dex (see an example of stacked halos in Figure 2). Noting that the stellar masses of the satellites are still not uniformly distributed within a bin and we do not pre-assume that the stellar metallicities follow a symmetric distribution or an exact Gaussian distribution at a given stellar mass, we take the medians111In practice, we have tested that our results and conclusions will be changed little if we take the means instead of the medians in each bin. of the variables ($\log M_{*}$,[Fe/H]) of all the satellites in each bin, and using the least-squares method to fit the medians as follows,

[TABLE]

where $\alpha$ is the best-fit slope and $b$ is the best-fit intercept at $\log(M_{*}/{\,M_{\odot}})=3$. In the fitting, the error bars of the medians of [Fe/H] are given by the half of the range between the 16th and 84th percentiles of the distribution of the [Fe/H] of the satellites. Each point in Figure 1 represents the best-fit result to the satellites in one host halo at $z=0$. Those systems in which the central galaxy is MW-like are marked in red, which are selected by the criteria described in Section 2. This figure demonstrates that the MZR of the satellites in MW-size halos is roughly universal, with the slopes and the intercepts of the satellites of the halos are distributed mostly within the same range. As seen from the figure, the fitted slopes $\alpha$ are scattered mostly in the range 0.2-0.4, and the fitted intercepts $b$ are scattered mostly in the range [-2.9,-2.5]. In the bottom panels with the relatively low halo mass, the slopes have relatively large scatters with large $\alpha$ up to 0.6, which is mainly because the number of the satellites are relatively small in a relatively low-mass halo and the corresponding statistical uncertainty is large. As seen from Figure 1, the MZRs of the satellites are not sensitive to the selection criteria for their central galaxies, as the red symbols representing MW-like central galaxies do not appear exceptional among the total sample.

3.2. The MZR of the satellites in a wide range of halo masses

In this subsection, we generalize the exploration of the MZR of the satellites to other halo masses (ranging from $10^{10}h^{-1}{\,M_{\odot}}$ to $10^{15}h^{-1}{\,M_{\odot}}$). Figure 2 shows the simulated MZRs of the satellites in the host halo masses ranging from $5\times 10^{11}h^{-1}{\,M_{\odot}}$ to $10^{15}h^{-1}{\,M_{\odot}}$, and each panel illustrates the stacked results of 10 merger trees of the same halo mass. The 10-tree results for $M{{}_{\rm halo}}<5\times 10^{11}h^{-1}M_{\odot}$ are not shown in this figure, as the number of their satellites is statistically small (the statistical results of their 1000 trees are listed in Table 2 below). As seen from this figure, the correlation between stellar metallicities and stellar masses of the satellites exists for the different halo masses.

As described in Section 3.1, we use the linear least-squares method to fit the MZR by Equation (7) for each tree shown in Figure 2. Note that the stellar metallicities start to become flat at the high-mass end ($M_{*}\gtrsim 10^{11}{\,M_{\odot}}$). The linear fitting is limited to the range $10^{3}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$. The best-fit parameters obtained for the satellites in 10 individual halos are shown in Figure 3(a)-(b). For high halo masses, the best-fit parameters of the individual halos do not scatter significantly. In low-mass halos, the small number statistics of the satellites in one halo can affect the best-fit parameters significantly. To view the results in a statistical way, we stack the satellites of the 10 halos together and get the best-fit results, and we show them in panels (c) and (d) of the figure. Panels (c) and (d) also include the best-fit results to the stacked 10000 or 1000 halos with low halo masses. Table 2 summarizes the MZR fitting results of the stacked halos.

According to the linear (single power law) fitting results shown in Figure 3(c)-(d) or listed in Table 2, we find the following results.

• •

Our fitting results are not affected much by the mass resolution set in the simulations, because the results obtained with the two different halo resolution masses $M{{}_{\rm res}}$ are close to being convergent for the same halo masses with the same stacked halo number $N=10$;

• •

The MZR slopes appear to be relatively lower in low-mass halos, compared to the slopes in high-mass halos, mainly because the low-mass halos have few high-mass satellites and those high-mass satellites play an important role in the fitting. For the halo masses that are not very high (e.g., $M{{}_{\rm halo}}=5\times 10^{11}$ or $10^{12}h^{-1}{\,M_{\odot}}$), an increase of the number of the stacked halos (e.g., from $N=10$ to $N=1000$) can increase the fitted MZR slopes, as the number of high-mass satellites increases. The increase of the slopes in high-mass halos should be taken to be caused by a physical effect in high-mass galaxies, as a double power law fitting to the MZR yields a higher slope in high-mass galaxies for both central galaxies and satellites in Section 3.3 and at the end of this subsection.

Figure 3(c)-(d) and Table 2 show that for the MZRs of the stacked satellites in different host halo masses, the average MZR slopes obtained by the linear (single power law) fit are generally in the same range $\sim$0.2–0.4 and slowly increase with increasing halo masses, and the average MZR intercept $b$ at $\log(M_{*}/{\,M_{\odot}})=3$ are generally in the range $\sim-2.3$ to $-3.1$.

The $\chi^{2}$ values shown in Table 2 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 linear (single power law) 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 $\sim 50\%-95\%$. However, we have the following notes on the use and the limitation of the results.

• •

As mentioned in Section 3.1, we do not pre-assume that the stellar metallicities follow a Gaussian distribution at a given stellar mass. For the scatters of data points not following an exact Gaussian distribution, the fit parameters and their errors are not the “minimum variance unbiased estimator” (MVUE) in the least-squares statistics, but the least-squares method is still helpful by being generalized to compare data with models (e.g., see Chapters 4, 6, and 11 in Bevington & Robinson 2003), which serves as a “first-order” approximation to understand the data quantitatively, as well as a mathematical way to compare the different models as done in this work.

• •

Some $\chi^{2}$ values are somewhat too small, compared with the degrees of freedom in the fitting (related with the bin number as described in the note of Table 2). Part of the reason for this is that a relatively large error has been assigned to the fitting data (i.e., using the large scatter of the [Fe/H] around the median in each bin).

• •

Understanding the fit parameters and their errors better will involve understanding the underlying physical processes/reasons leading to the data, which should be taken into account in the comparison of observational data with simulation results.

As the slopes obtained by a linear fitting to the MZR slightly increases with increasing halo masses within the range $\sim 0.2$–0.4, we perform a double power law fit for the satellites in the high-mass halos ($M{{}_{\rm halo}}=10^{13},10^{14},10^{15}h^{-1}{\,M_{\odot}}$), with a break at $M_{*}=10^{8}{\,M_{\odot}}$ in the stellar mass range $10^{3}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$. The fitting lines are continuous at the break point. The fitting is done by the least-squares method as described in Section 3.1. The fitting results are shown in Table 2 and Figure 2(d)-(f). We do not perform a double power law fit for low-mass halos, as their satellite masses are not high enough. As seen from the double power law fitting results shown in Table 2, we find that the double power law fitting yields a much smaller $\chi^{2}$ value, which indicates a better fit than the single power law fitting (though the $\chi^{2}$ values are too small compared with the degrees of freedom, which is caused by the large errors assigned to the fitting data, as mentioned above). In the low-mass range with $10^{3}{\,M_{\odot}}<M_{*}<10^{8}{\,M_{\odot}}$, the fit slopes ($\sim 0.31$–$0.35$) become smaller than their single power law fitting results and closer to the linear fit slopes in the low-mass halos; and the difference of the intercepts from those in low-mass halos also becomes smaller, which is $\lesssim 0.6$ dex at $\log(M_{*}/{\,M_{\odot}})=3$ or $\lesssim 0.2$ dex at $\log(M_{*}/{\,M_{\odot}})=6$. In the high-mass range with $10^{8}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$, the slopes obtained by the double power law fitting is roughly universal, $\sim 0.5$, which are higher than their single power law fitting results. The results of the double power law fitting to high-mass halos supports that there is no significant difference in the MZRs of the satellites in host halo masses $\sim 10^{11}$–$10^{15}h^{-1}{\,M_{\odot}}$, though the MZR slopes in lower halo masses are still slightly smaller. The results of the double power law fitting support a universal MZR in satellites.

The slight decrease of the slope with decreasing halo masses existing at the low-halo mass range of $M_{\rm halo}\sim 10^{10}$–$10^{12}h^{-1}{\,M_{\odot}}$ is associated with the monotonically increasing intercept $b$ in the same halo mass range shown in Figure 3(f) (or Figs. 4c and 4d below), which is related to the metallicity enrichment caused by SNe Ia, as well as different halo assembly histories with different halo masses and their different star formation histories at high redshifts. As discussed in section 3.1.1 in HYL14, the scatter of the simulated MZRs is generally large at the low-$M_{*}$ end, which is caused mainly through the scatter in the star formation durations of the galaxies at a given stellar mass and the difference in the chemical enrichment of SNe Ia and II. If the duration is short, the metal enrichment is mainly contributed by SN II explosions, which have a chemical pattern with a relatively low iron fraction; if the duration is long enough, SNe Ia may have a non-negligible contribution to the metal enrichment and generate more iron than SNe II. Thus, a short star formation duration would lead to a low [Fe/H], while a longer one leads to a higher [Fe/H]. High-mass satellites generally all experienced extended star formation duration; low-mass satellites are generally formed at high redshifts and some of them had relatively shorter star formation duration. The low-mass satellites with relatively shorter star formation duration have relatively lower metallicities, and this effect appears to be slightly more significant in high-mass halos. We have also done the test to obtain the MZR results by removing the metallicity enrichment due to SN Ia in the semianalytical galaxy formation and evolution model and find that the slight decrease of the slopes and the increase of the intercepts with decreasing halo masses in low-mass halos with $M_{\rm halo}\sim 10^{10}$–$10^{12}h^{-1}{\,M_{\odot}}$ disappear (where the intercepts of the satellite MZRs in low-mass halos become $\sim-2.7$ to $-2.9$).

3.3. The MZR of central host galaxies

We obtain the MZR for the simulated central galaxies with different host halo masses. Our model results and the comparison with observations are shown in Figure 4. Figure 4(a) shows the MZR of the simulated central galaxies in different host halo masses, where the stellar metallicity increases with increasing stellar masses and becomes flat at $\log(M_{*}/{\,M_{\odot}})\gtrsim 11$. The dotted lines are the fitting results to the data obtained from the simulation runnings in the left column of Table 1. The fitting is performed in these two ways: one is fit by one linear slope within the stellar mass range $10^{3}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$ (brown dotted line), and the other is fit by two slopes with a break at $M_{*}=10^{8}{\,M_{\odot}}$ in the stellar mass range $10^{3}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$ and a constant at $M_{*}>10^{11}{\,M_{\odot}}$ (black dotted line), where the fitting lines are continuous at the break points. The fitting is done by the least-squares method as described in Section 3.1. A relatively small bin size with 0.5 dex in $\log M_{*}$ (compared to 1-dex interval used for the satellites above) is used due to relatively more high-mass systems here. The best-fit slopes and intercepts are listed in Table 2. As seen from Table 2, the single power law fit gives the MZR slope $\alpha=0.30\pm 0.01$ with $b=-2.52\pm 0.05$, which is in the same range as the MZR slopes of their satellites. The double power law fit gives a steeper slope ($\alpha=0.47\pm 0.02$) in the mass range $10^{8}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$ than that ($\alpha=0.24\pm 0.02$) in the relatively low-mass range $10^{3}{\,M_{\odot}}<M_{*}<10^{8}{\,M_{\odot}}$. This tendency is consistent with the above MZR fitting results for satellites described in Section 3.2. The slope of the central galaxies in the low-$M_{*}$ range ($\alpha=0.24\pm 0.02$) is at the low end of the slope range of the satellites ($\alpha\sim 0.2$–0.4), indicating that many low-mass satellites are relatively more metal poor. In the high-mass range with $10^{8}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$, the best-fit slope is roughly the same as the corresponding values of the satellites obtained in the double power law fit. The metallicity becomes roughly constant at $M_{*}\gtrsim 10^{11}{\,M_{\odot}}$, with [Fe/H]$\sim 0.26$, which is mainly determined by the Fe yield of SNe II and Ia. As seen from Table 2, the $\chi^{2}$ value of the double power law fit for the central galaxies are smaller than that of the single power law fit, which suggests that the double power law fit is acceptable with a higher probability. Figure 4(b) shows the observational results of some Local Group dwarfs and nearby galaxies (K13; Gallazzi et al. 2005). Figure 4(c) shows that the MZRs of the simulated central galaxies (dotted line) and the simulated satellites (solid lines; same as those in Figure 2) are universally distributed within the same space. Figure 4(d) shows that our simulation results of the MZR generally agree well with the observational data.

When more observational results with good measurement quality come out and support the nonlinear behavior in the MZR or the double power law fitting does not perform well, it would be helpful to perform a nonparametric or semi-parametric regression, which requires large sample sizes, on both observational data and simulation data to uncover the underlying data structure and the contributions from different physical processes or high-order physical effects (e.g., see some modern statistical methods in Feigelson & Babu 2012; Takezawa 2005; Ruppert et al. 2003).

3.4. The MZR in Local Group dwarfs

Our simulation results obtained above support that Local Group dwarfs follow a universal MZR at least by the following points. (1) As shown in Figure 1, the MZR in the satellites of MW-size halos is universal. M31 is another big galaxy in the Local Group, and the host halo mass of M31 is also possibly in the range $\sim 1$–$2\times 10^{12}{\,M_{\odot}}$ (e.g., Kafle et al. 2018 and references therein). The results shown in Figures 1 and 2 support that the MZR of the satellites in both M31 and MW follows a universal MZR. (2) Figure 4 shows that the central galaxies of smaller halos with $M{{}_{\rm halo}}<10^{11}h^{-1}{\,M_{\odot}}$ follow the same universal MZR as that for satellites in the MW-size halos. Other dwarf galaxies in the Local Group can be explained as the central galaxies of smaller halos, so they also follow the universal MZR.

4. Summary and Discussion

We have investigated the stellar mass–stellar metallicity relations of galaxies and their satellites by using the semi-analytical models. Our study suggests a universal correlation of the stellar metallicity with stellar mass $M_{*}\sim 10^{3}{\,M_{\odot}}$–$10^{11}{\,M_{\odot}}$, [Fe/H]=$\alpha\log(M_{*}/10^{3}{\,M_{\odot}})+b$, and a roughly constant stellar metallicity for $M_{*}\gtrsim 10^{11}{\,M_{\odot}}$. The relations reproduced from our work are consistent with the observations on the MW/M31 dwarf satellites ($M_{*}\sim 10^{3}$–$10^{8.5}{\,M_{\odot}}$) and local galaxies ($M_{*}\sim 10^{9}{\,M_{\odot}}$–$10^{12}{\,M_{\odot}}$), as shown in Figure 4.

Our study shows that the slope $\alpha$ of the MZR for the satellites in a halo with mass the same as MW/M31 halo mass ([Fe/H]–$\log M_{*}$) are mostly in the same range of 0.2–0.4 (i.e., with a scatter of $\sim 0.2$ for different halo assembly histories), and the intercept at $\log(M_{*}/{\,M_{\odot}})=3$ is $b\sim-2.5$ to $-2.9$. If the satellites of many halos are stacked together, the MZR for the stacked satellites in MW/M31-size halos gives $\alpha\simeq 0.35\pm 0.03$ and the intercept $b\simeq-2.7\pm 0.1$. The slopes of the stacked satellites in the host halo masses ranging from $10^{10}h^{-1}{\,M_{\odot}}$ to $10^{15}h^{-1}{\,M_{\odot}}$ are correspondingly in the same range of 0.2–0.4, with a slight increase with increasing halo mass; and the intercept $b$ at $\log(M_{*}/{\,M_{\odot}})=3$ is in the range $-2.3$ to $-3.1$.

Our study shows that the slope of the correlation for the central galaxies with $M_{*}$ from $\sim 10^{3}$ to $10^{11}{\,M_{\odot}}$ is $\alpha\simeq 0.30\pm 0.01$ with $b=-2.52\pm 0.05$, which is in the same range of the satellite MZRs. The MZR becomes roughly constant at $M_{*}\gtrsim 10^{11}{\,M_{\odot}}$, with [Fe/H]$\sim 0.26$, which is mainly determined by the Fe yield of SNe II and Ia.

Our study also shows that a double power law provides a better fit to the stellar mass–stellar metallicity relation than the single power law for both satellites and central galaxies. The double power law fit 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 difference in the best-fit intercepts of the correlations is $\lesssim 0.6$ dex at $\log(M_{*}/{\,M_{\odot}})=3$ (with $b\sim-2.9$ to $-2.3$) or $\lesssim 0.2$ dex at $\log(M_{*}/{\,M_{\odot}})=6$. Specifically, the double power law fit to the central galaxies gives $\alpha\simeq 0.24\pm 0.02$ (close to the low end of the range 0.2–0.4) and $b\simeq-2.38\pm 0.06$ at $10^{3}{\,M_{\odot}}<M_{*}<10^{8}{\,M_{\odot}}$, and $\alpha\simeq 0.47\pm 0.02$ and $b\simeq-3.5\pm 0.1$ at $10^{8}{\,M_{\odot}}<M_{*}<10^{11}{\,M_{\odot}}$. The high-mass satellites ($M_{*}\gtrsim 10^{8}{\,M_{\odot}}$) existing mostly in high-mass halos and their relatively high $\alpha$ can lead to an apparent increase of $\alpha$ with increasing host halo masses obtained in the single power law fitting to the satellites; and after taking into account that effect, the dependence of the satellite MZRs on their host halo masses becomes little in the halo mass range $M{{}_{\rm halo}}\sim 10^{11}$–$10^{15}h^{-1}{\,M_{\odot}}$, though the slopes in lower halo masses are still slightly smaller.

For the MZRs of satellites in low-mass halos with $M_{\rm halo}\sim 10^{10}$–$10^{12}h^{-1}{\,M_{\odot}}$, the slight decrease of the slope with decreasing halo masses is associated with their monotonically increasing intercept $b$ with decreasing halo masses, which is related to the metallicity enrichment caused by SNe Ia, as well as different halo assembly histories with different halo masses and their different star formation histories at high redshifts.

Although a detailed examination of the simulated MZRs may suggest some slight differences among the central galaxies and their satellites with host halo masses spanning a wide range ($M{{}_{\rm halo}}\sim 10^{9}$–$10^{15}h^{-1}{\,M_{\odot}}$), we call the MZRs “universal” because their slopes and normalizations are within a relatively narrow range. The universal relation of the MZRs in the satellites within the large range of halo masses and their central galaxies awaits future observational tests. A precise comparison with observational results would require the consideration of observational selection effects and similar fitting methods. A double power law or non-linear fitting to the MZRs would become necessary if the number of the galaxies are sufficiently large and their masses cover a large range.

Active galactic nuclei (AGN) feedback is not included in the model. AGN feedback can have significant effects mainly in high-mass galaxies (e.g., $\gtrsim 10^{11}{\,M_{\odot}}$; Croton et al. 2006; Bower et al. 2006; Somerville & Davé 2015), which decreases the galaxy mass function at the high-mass end. As shown in this work, the stellar metallicity increases with increasing stellar mass mainly in galaxies with $M_{*}\lesssim 10^{11}{\,M_{\odot}}$, and the stellar metallicity becomes roughly constant at higher masses. AGN feedback would not have significant effects in the monotonically increasing part of the stellar metallicity with the stellar mass.

For simplicity, the reionization epoch is set to the same for the different halo masses in this work so that the effects due to other physical processes can be isolated. A late reionization epoch may result in relatively low metallicities at the low-mass end of the MZR and a steeper slope, as illustrated for MW satellites in HYL14. Future observations on the slope of the MZR for dwarfs would provide a constraint on whether or how the reionization epoch is different in different environments.

A modification to the cooling recipe in semianalytical models can result in a different luminosity function of galaxies. However, the universality in the MZR revealed in this paper and its compatibility with observations suggest that the modification itself should not be the key to the origin of the universality.

The universality of the MZR relations in satellites and host galaxies with different host halo masses provides insights or constraints on understanding the galaxy formation and evolution and their chemical evolution. One next important step of this work is to understand the origin of the universality of the MZR relation, investigate its evolution with redshift, 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 explore whether there exists a possible dependence on a third or more parameters (e.g., star formation rate, as shown in Ellison et al. 2008; Mannucci et al. 2010; Barrera-Ballesteros et al. 2017), where the multivariate Principal Components Analysis or a nonlinear extension that treats more variables would be useful. K13 discussed some possible scenarios for the MZR relations in the MW/M31 satellites. Finlator (2016) give a review on the scenarios involving inflows and outflows to explain the MZRs in massive galaxies. The semianalytical model used in this work provides an efficient way to isolate the effects of different physical processes and to see the dominant reason leading to the MZRs. The underlying reasons for the universality of the MZR obtained from the semianalytical models will be reported in a different work (in preparation).

We thank the referee, Sandra Faber, Youjun Lu, and Eric Peng for helpful comments on the manuscript. 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).