© ESO 2021

1. Introduction

Star-forming galaxies evolve and assemble their mass by accreting gas from the intergalactic medium (IGM), converting it in stars, then expelling gas and metals, thus changing the gas metallicity in the interstellar medium (ISM). This baryon cycle is mainly regulated by internal processes¹ (e.g., Wu et al. 2017), which are responsible for metal and dust production, and the subsequent expulsion of material. Outflows driven by stellar winds and supernovae (SNe) explosions dominate in the lower-mass regime, while at higher masses, outflows tend to be powered by feedback from active galactic nuclei (AGN, e.g., Tremonti et al. 2004; Dekel & Birnboim 2006; Dalcanton 2007; Erb 2008; Finlator & Davé 2008; Zahid et al. 2014; De Rossi et al. 2017; Lara-López et al. 2019).

The efficiency of these processes is predicted to change with galaxy properties, and in particular with the virial mass (M_vir) and stellar mass (M_star) of the galaxies. These variations lead to the most important relations among the main galaxy parameters driving galaxy evolution, namely the “main sequence” correlation between M_star and the star-formation rate (SFR, e.g., Brinchmann et al. 2004, the SFMS); the mass–metallicity relation (e.g., Tremonti et al. 2004, the MZR); and the consequent inverse correlation between metallicity and specific SFR (sSFR ≡ SFR/M_star, e.g., Mannucci et al. 2010; Lara-López et al. 2010; Hunt et al. 2012, 2016). Baryonic processes are thought to shape the SFMS and the MZR, in addition to many other correlations among other stellar properties and dark matter: the virial-to-stellar mass (e.g., Moster et al. 2010) and the baryonic Tully–Fisher (e.g., Lelli et al. 2016a) relations; the correlation between color and stellar population gradients with mass (e.g., Tortora et al. 2010); and the relation of dark matter fraction and mass density slopes with stellar mass (e.g., Tortora et al. 2019).

These scaling relations are not necessarily straight power laws, but can also show curvatures, bending, breaks, and U-shapes, suggesting the existence of typical mass scales which regulate the efficiency of the physical processes involved in baryonic cycling (e.g., Dekel & Birnboim 2006). Kannappan et al. (2013) identified two transitions in galaxy morphology, gas fractions, and fueling regimes. The first transition was considered a “gas-richness” threshold; galaxies below a transition mass of log(M_star/M_⊙) ∼ 9.5 were found to be “accretion-dominated”, characterized by overwhelming gas accretion. The second M_star transition corresponds, instead, to a bimodality threshold, roughly indicating a “quenched” regime at log(M_star/M_⊙) > 10.5, populated by spheroid-dominated, gas-poor galaxies.

Both M_star thresholds, and the mass regimes they delineate, are really tracing gas content, and its intimate link with star-formation activity (Kannappan et al. 2013). To further explore this link, in the first paper of this series (Ginolfi et al. 2020, hereafter Paper I), we constructed a sample of galaxies in the Local Universe with homogeneously derived global measurements of the main parameters that quantify galaxy evolution: namely M_star, SFR, metallicity O/H, and gas content including molecular gas mass M_H2 and atomic gas mass M_HI. The sample comprises ∼400 galaxies, and was dubbed Metallicity And Gas for Mass Assembly (MAGMA).

In a second paper (Hunt et al. 2020, hereafter Paper II), we used MAGMA to explore how molecular and atomic gas properties vary as a function of M_star and SFR. The main focus of Paper II was on feedback processes, and how they govern scaling relations across the mass spectrum. Essentially three processes play a role in regulating galaxy evolution²: (P1) preventive feedback, the (lack of) availability of cold baryons from the host halo; (P2) inefficiency of the star-formation process, the conversion of the available gas into stars; and (P3) ejective feedback, the production of energy and momentum, and the expulsion of material including gas, dust, and metals. In Paper II, we investigated two of the three main feedback mechanisms behind the inefficiency of low-mass halos to form stars (e.g., Behroozi et al. 2013; Moster et al. 2013; Graziani et al. 2015, 2017): namely the availability of cold gas to form stars (P1), and the inefficiency of the star-formation process itself (P2).

Here, in the third paper of the series, we examine the remaining mechanism, (P3), the efficiency of metal enrichment and ejective feedback. Understanding the processes that shape galaxy scaling relations with metallicity has proved to be a challenging task. It has been historically important to compare observations with the “closed-box” model, assuming that the total baryon content in the galaxy is fixed and the gas and stellar content are only determined by star formation and the stellar yields (e.g., Pagel & Patchett 1975). However, it is widely recognized that galaxies do not evolve as closed boxes (e.g., Tremonti et al. 2004), and more complex models have considered inflows and outflows of material from the galactic ecosystem (e.g., Larson 1972; Pagel & Patchett 1975; Tinsley 1980; Pagel & Edmunds 1981; Edmunds & Pagel 1984; Edmunds 1990; Tremonti et al. 2004; Dalcanton 2007; Erb 2008; Finlator & Davé 2008; Recchi et al. 2008; Peeples & Shankar 2011; Dayal et al. 2013; Zahid et al. 2014; Peng & Maiolino 2014). More recently, hydrodynamical simulations have tested these predictions (e.g., Ma et al. 2016; De Rossi et al. 2017; Torrey et al. 2019).

Our approach here is mainly empirical, requiring a sample with measurements of M_star, SFR, O/H, and gas content, homogeneously determined over a wide range of stellar mass and gas fractions; MAGMA is such a sample. In this paper, Sect. 2 describes briefly the MAGMA sample and how it was assembled. Then, Sect. 3 presents the theoretical formalism for metallicity evolution and the scaling of metallicity with M_star, SFR, and gas content. The computation of outflow parameters and metal-loading factors for MAGMA galaxies is described in Sect. 4, and their impact on shaping the mass–metallicity relation in Sect. 5. We recast our results in the context of the effective yields in Sect. 6, and discuss the physical implications for wind metallicity and mass loading of the winds in Sect. 7. Our conclusions are drawn in Sect. 8.

2. MAGMA sample

Paper I presented the MAGMA sample, comprising 392 local galaxies having homogeneously derived global quantities including M_star, SFR, metallicities [12 + log(O/H)], and gas masses (both atomic, M_HI, and molecular, M_H2, and total M_g = M_HI + M_H2). Molecular gas masses were obtained by measurements of the CO luminosity, $L_{\text{CO}}^{\prime }$. We assembled this sample by combining a variety of previous surveys at z ∼ 0 with new observations of CO in low-mass galaxies. Only galaxies with robust detections of physical parameters were retained, that is to say we required clear detections of CO and H I. Galaxies were eliminated if they showed optical signatures of an AGN in the parent samples (obtained from line-ratio diagnostics, similar to Baldwin et al. 1981). We also eliminated potentially H I-deficient galaxies by applying a cut based on H I-deficiency measurements, where available (e.g., Boselli et al. 2009, 2014a). The MAGMA sample spans more than 5 orders of magnitude in M_star, SFR, and M_g, and a factor of ∼50 in metallicity [Z, 12 + log(O/H)].

Because of potential systematics that could perturb our results, in Paper I we homogenized the sample by recalculating both M_star and SFR in a uniform way. To calculate SFR, we adopted the hybrid formulations of Leroy et al. (2019), based on combinations of GALEX and WISE luminosities. New M_star measurements were computed using 3.4 μm (or IRAC 3.6 μm when not available) luminosities from fluxes in the ALLWISE Source Catalogue (Wright et al. 2010) and a luminosity-dependent M/L from Hunt et al. (2019); the continuum contribution from free-free emission was subtracted from these luminosities before the M_star computation. Both SFR and M_star assume a Chabrier (2003) IMF.

For metallicity, 12 + log(O/H), we adopted either “direct” calibrations, based on electron temperatures (e.g., Izotov et al. 2007), or the [N II]-based strong-line calibration by Pettini & Pagel (2004, PP04N2). When PP04N2-derived metallicities were not available, they were converted from the original calibration to PP04N2 according to Kewley & Ellison (2008). PP04N2 is the metallicity calibration closest to the direct one (e.g., Andrews & Martini 2013; Hunt et al. 2016), which makes it appropriate for a sample like MAGMA that spans a wide range of M_star and SFR.

The molecular gas mass was determined from $L_{\text{CO}}^{\prime }$ based mainly on ¹²CO(1 − 0); only in a few cases did we have to convert from ¹²CO(2 − 1) to ¹²CO(1 − 0) using “standard” prescriptions (e.g., Leroy et al. 2009; Schruba et al. 2012). The conversion factor α_CO from CO luminosity to molecular gas mass was calculated according to the piecewise relation presented in Paper II, with α_CO ∝ Z^1.55, for 12 + log(O/H) < 8.69 and ${\alpha }_{\mathrm{CO}}^{\odot }=3.2M_{\odot }$ (K km s⁻¹ pc²)⁻¹ (see also Saintonge et al. 2011). The impact of a different metallicity dependence is discussed in Paper II. Here we do not include a factor of 1.36 to account for helium. For more details about the parent surveys, as well as parameter determinations and calibrations, quality checks, or the derivation of scaling relations, we refer the reader to Paper I and Paper II.

3. Metallicity evolution and scaling relations: formalism

Much effort has been spent over the last fifty years to understand the theoretical framework behind the baryon exchange cycle and metallicity evolution in galaxies (e.g., Larson 1972; Pagel & Patchett 1975; Tinsley 1980; Pagel & Edmunds 1981; Edmunds & Pagel 1984; Edmunds 1990; Tremonti et al. 2004; Dalcanton 2007; Erb 2008; Finlator & Davé 2008; Recchi et al. 2008; Spitoni et al. 2010; Peeples & Shankar 2011; Dayal et al. 2013; Zahid et al. 2014; Peng & Maiolino 2014). It has been fairly well established that bulk outflows and galactic winds driven by star formation are necessary to explain the MZR and other scaling relations involving metallicity. Here we follow these pioneering works, and explore the behavior of the MZR, SFR, and gas content in the MAGMA galaxies. MAGMA provides the important constraints of gas measurements together with SFR and metallicity, thus enabling, for the first time, a quantification of the mass and metallicity loading in stellar outflows, and their role in shaping the MZR in the Local Universe.

3.1. Galactic chemical evolution

In what follows, we rely heavily on the formulations for galactic chemical evolution by Pagel (2009). We first make the assumption of instantaneous recycling, namely that the chemical enrichment provided by stellar evolution, nucleosynthesis, and expulsion of metals into the ISM take place “instantaneously” relative to the timescales of overall galaxy evolution.

There are two main quantities that this assumption makes time-invariant: the first is the return fraction R, the mass fraction of a single stellar generation that is returned to the ISM; the second is the (true) stellar yield, y, defined as the mass of an element newly produced and expelled by a generation of stars relative to the mass that remains locked up in long-lived stars and compact remnants. The metal mass produced relative to the initial mass of a single stellar generation is then q = α y where α = 1 − R, (or “lock-up” fraction), with the implicit assumption that stellar lifetimes can be neglected so that R is a net mass return fraction (see e.g., Vincenzo et al. 2016).

These parameters are governed by the IMF, and in particular the mass fraction of a stellar generation above a certain limit, since these are the main contributors to metal production. We further assume that metals are homogeneously mixed at all times, and that initally there are no stars and no metals, only the initial gas M_g(t = 0) with mass M_i.

There are four main physical quantities involved in and constrained by chemical evolution: the mass in stars, M_star; the star-formation rate SFR, here referred to as ψ; the mass in gas, M_g; and the mass in metals, M_Z, or the metal abundance in the gas, Z_g = M_Z/M_g. The mass of stars that have been born up to time t can be written as:

$$\begin{array}{c}\hfill S(t)={\displaystyle {\int }_0^t\psi (t\prime )\mathrm{d}{t}^{\prime }.}\end{array}$$(1)

Consequently, the mass M_star remaining in stars and long-lived remnants is:

$$\begin{array}{c}\hfill M_{\star }(t)=\alpha S(t).\end{array}$$(2)

The total baryonic mass in a galactic system is the sum of gas and stars:

$$\begin{array}{c}\hfill M_{\mathrm{tot}}(t)=M_{\star }(t)+M_{\mathrm{g}}(t)=M_{\mathrm{i}}-M_{\mathrm{w}}(t)+M_{\mathrm{a}}(t),\end{array}$$(3)

where M_w is the mass expelled in bulk winds and outflows, and M_a is the mass of newly accreted gas. This is essentially a mass continuity equation that defines the mass conservation between infall, outflows, and star formation. Thus, the time derivative of the gas mass can be written as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{g}}={\dot{M}}_{\mathrm{a}}-{\dot{M}}_{\mathrm{w}}-{\dot{M}}_{\star }.\end{array}$$(4)

The time variation of the mass of metals M_Z in the ISM gas depends on ψ, R, the stellar yield q, Ṁ_a, and Ṁ_w as follows:

$$\begin{array}{c}\hfill {\dot{M}}_Z=Z_{\mathrm{g}}R\psi +q\psi -Z_{\mathrm{g}}\psi +Z_{\mathrm{a}}{\dot{M}}_{\mathrm{a}}-Z_{\mathrm{w}}{\dot{M}}_{\mathrm{w}},\end{array}$$(5)

where Z_a and Z_w are the metallicities of the accreted and outflow gas, respectively. The first two terms on the right-hand side of Eq. (5) correspond to the addition of metals of ejecta from stars; the third term to the loss to the ISM from astration; and the remaining terms to the metals gained from accretion and lost in bulk stellar outflows from winds and SNe.

Equation (5) can be simplified by two common assumptions: we take Ṁ_a and Ṁ_w to be directly proportional to ψ (e.g., Erb 2008; Peeples & Shankar 2011; Dayal et al. 2013; Lilly et al. 2013):

$$\begin{array}{cc}& {\eta }_{\mathrm{a}}={\dot{M}}_{\mathrm{a}}/\psi ,\hfill \\ \hfill & {\eta }_{\mathrm{w}}={\dot{M}}_{\mathrm{w}}/\psi ,\hfill \end{array}$$(6)

where we have introduced the mass-loading factors η_a and η_w. It is fairly well established that the latter relation is true, namely that the mass in outflows is proportional to SFR, at least to its integral over time (e.g., Muratov et al. 2015; Christensen et al. 2016). In contrast, the former proportionality is more subject to speculation, although evidence is mounting that gas accretion powers star formation (e.g., Fraternali & Tomassetti 2012; Sánchez Almeida et al. 2014; Tacchella et al. 2016). In any case, this assumption provides the noteable advantage of analytical simplicity; we plan to investigate more completely its ramifications in a future paper. Following Peeples & Shankar (2011), we also introduce the metal-loading factors ζ_a and ζ_w relating Z_g, Z_a, η_a, Z_w, and η_w:

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{a}}\equiv \left(\frac{Z_{\mathrm{a}}}{Z_{\mathrm{g}}}\right)\frac{{\dot{M}}_{\mathrm{a}}}{\psi }=\left(\frac{Z_{\mathrm{a}}}{Z_{\mathrm{g}}}\right){\eta }_{\mathrm{a}},\end{array}$$(7)

and

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{w}}\equiv \left(\frac{Z_{\mathrm{w}}}{Z_{\mathrm{g}}}\right)\frac{{\dot{M}}_{\mathrm{w}}}{\psi }=\left(\frac{Z_{\mathrm{w}}}{Z_{\mathrm{g}}}\right){\eta }_{\mathrm{w}}.\end{array}$$(8)

ζ_w and ζ_a essentially describe the efficiency with which gas accretion and star-formation-driven outflows change metal content within the galaxy. Thus, we can rewrite Eqs. (4) and (5) as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{g}}=\psi ({\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}-\alpha ),\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\dot{M}}_Z=\psi [q-\alpha Z_{\mathrm{g}}+Z_{\mathrm{g}}({\zeta }_{\mathrm{a}}-{\zeta }_{\mathrm{w}})].\end{array}$$(10)

Our observations do not directly constrain M_Z, but rather Z_g, so we need to reformulate Eq. (10):

$$\begin{array}{c}\hfill {\dot{M}}_Z=\frac{\mathrm{d}}{\mathrm{d}t}\left(M_{\mathrm{g}}{Z}_{\mathrm{g}}\right)=M_{\mathrm{g}}{\dot{Z}}_{\mathrm{g}}+Z_{\mathrm{g}}{\dot{M}}_{\mathrm{g}}\end{array}$$(11)

in order to solve for Ż_g:

$$\begin{array}{c}\hfill {\dot{Z}}_{\mathrm{g}}=\frac{\psi }{M_{\mathrm{g}}}[q+Z_{\mathrm{g}}({\zeta }_{\mathrm{a}}-{\zeta }_{\mathrm{w}}-{\eta }_{\mathrm{a}}+{\eta }_{\mathrm{w}})],\end{array}$$(12)

where we have used Eq. (9) for Ṁ_g and Eq. (10) for Ṁ_Z.

The necessary ingredients are now in place to derive the solution for the dependence of Z_g on the remaining parameters. First, we know that SFR is related to gas mass by:

$$\begin{array}{c}\hfill \psi ={ϵ}_{\star }{M}_{\mathrm{g}}\end{array}$$(13)

where ϵ_⋆ is the inverse depletion time or time-scale for star formation. Equations (9) and (13) are coupled, so must be solved simultaneously. Assuming that ϵ_⋆ does not vary with time, we can thus integrate Eq. (9), to obtain the following solution for the gas mass:

$$\begin{array}{c}\hfill M_g(t)=M_{\mathrm{i}}{e}^{t{ϵ}_{\star }({\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}-\alpha )}.\end{array}$$(14)

Substituting this solution into Eq. (12) and integrating, we obtain:

$$\begin{array}{c}\hfill Z_{\mathrm{g}}=\left(\frac{q}{{\zeta }_{\mathrm{w}}-{\zeta }_{\mathrm{a}}+{\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}}\right)[1-{\left(\frac{M_{\mathrm{g}}}{M_{\mathrm{i}}}\right)}^{\frac{{\zeta }_{\mathrm{a}}-{\zeta }_{\mathrm{w}}-{\eta }_{\mathrm{a}}+{\eta }_{\mathrm{w}}}{{\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}-\alpha }}],\end{array}$$(15)

where it can be seen that the time dependence is encapsulated in the evolution of the gas (M_g/M_i, e.g., Recchi et al. 2008). Equation (15) is essentially telling us that galaxies evolve along a hypersurface constrained by M_star, M_g, gas metallicity, and accretion and wind mass and metal loadings. The galaxy’s position on this hypersurface can change with time as gas is accreted, stars are formed, and mass and metals are ejected in winds. However, at each time, the relation constraining these fundamental variables is approximately invariant.

The ratio of the system gas mass M_g and the initial gas mass M_i can be written in terms of the baryonic gas mass fraction, μ_g = M_g/(M_g + M_⋆):

$$\begin{array}{c}\hfill \frac{M_{\mathrm{g}}}{M_{\mathrm{i}}}=\frac{{\mu }_{\mathrm{g}}}{1+({\mu }_{\mathrm{g}}-1)\left(\frac{{\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}}{\alpha }\right)}\end{array}$$(16)

where we have assumed that the total baryonic mass is described by Eq. (3). This is again an expression of mass conservation, meaning that no gas has been expelled altogether from the galaxy, but rather has been recycled into stars as described in Eq. (13).

Similar solutions for Z_g were first formulated by Edmunds & Pagel (1984), and later by Recchi et al. (2008), Erb (2008), Spitoni et al. (2010), Peeples & Shankar (2011), Dayal et al. (2013), and many others. These earlier works mostly assumed that Z_w/Z_g = 1 in a so-called homogeneous wind (e.g., Pagel 2009), leading to ζ_w = η_w. This assumption simplifies Eq. (15) placing the dependence of η_w only in the exponent of M_g/M_i. However, as has been shown previously (e.g., Peeples & Shankar 2011) and as we show below for MAGMA, this is almost certainly an over-simplified assumption because the metallicity of the outflowing material Z_w is not generally the same as the ISM gas, Z_g, especially at low masses. We further discuss this simplification and its implications later on in the paper (see also Appendix C).

3.2. Limitations of the model

First, Eq. (15) is highly sensitive to the value of η_a − η_w − α in the denominator of the exponent of M_g/M_i. For convenience in what follows, we define:

$$\begin{array}{c}\hfill \mathrm{\Delta }\eta ={\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}.\end{array}$$(17)

When Δη = α, the exponent for M_g/M_i in Eq. (15) becomes infinite. This equality, Δη ≈ α, is exactly what would be expected for the “equilibrium scenario” for gas and metallicity evolution, in which gas accretion is balanced by star formation and winds (e.g., Davé et al. 2012; Mitra et al. 2015), such that Ṁ_g ≈ 0. From Eq. (9), this implies that either ψ ≈ 0 which is not physically relevant here, or Δη ≈ α. In other words, for “equilibrium”, the difference between the accretion η_a and wind η_w coefficients is balanced by the lock-up fraction α. However, this formalism gives for the equilibrium scenario a solution for Z_g that is apparently ill-conditioned³. When Δη = α, Eq. (14) (and Eq. (16)) are irrelevant, as M_g = M_i.

More generally, the trend of gas growth in Eq. (14) depends on the sign of the exponent Δη − α. This exponent is a critical component of Eq. (15) because of the sensitivity to the sign of the exponent of M_g/M_i. Where Ṁ_g > 0, with Δη − α slightly positive, a finite Z_g solution is obtained only when this exponent is negative, and conversely, if Δη − α is slightly negative, then this exponent should be positive.

Equation (15) is also sensitive to the gas fraction μ_g through Eq. (16). As also noticed by Recchi et al. (2008), this approach leads to constraints on the possible values of gas fraction μ_g. Depending on the value of Δη − α, very small μ_g leads to infinite solutions for Z_g.

In any case, despite the complexity of possible solutions, Eq. (15) is a useful tool to assess the interplay of gas and metallicity in galaxies. In the next section, we explore this by letting the data guide the possible solutions.

4. Mass and metal loading in MAGMA

MAGMA combines homogeneous measurements of M_star, SFR, Z_g, and M_g; so, for the first time, we can infer from observations the mass loading factors, η_a and η_w, as well as the metal-loading factor ζ_w. To do this, we first make another common assumption, namely that the accreted gas is pristine (Z_a = 0, ζ_a = 0, e.g., Erb 2008; Peeples & Shankar 2011; Dayal et al. 2013; Lilly et al. 2013; Creasey et al. 2015). Although this assumption may not be exactly true for low-redshift galaxy populations, there is considerable observational evidence that the gas accreted by local galaxies is significantly metal poor (e.g., Tosi 1982; Sánchez Almeida et al. 2014; Fraternali 2017).

To apply the formalism described above to the observations, we convert the MAGMA oxygen abundance 12 + log(O/H) to Z_g according to:

$$\begin{array}{c}\hfill log(Z_{\mathrm{g}})=log(\mathrm{O}/\mathrm{H})+0.956,\end{array}$$(18)

where the underlying assumption is that oxygen measures the enrichment of a primordial mix of gas comprising approximately 75% hydrogen by mass (e.g., Peeples & Shankar 2011; Chisholm et al. 2018). We adopt the yields and return fractions for the Chabrier (2003) IMF from Vincenzo et al. (2016) using the Romano et al. (2010) stellar yields: (y, R) = (0.037, 0.455). We have purposely not considered a conversion of O/H to total metal abundance Z_g, because of the inherent uncertainties in this conversion (e.g., Vincenzo et al. 2016); thus our value for y is the oxygen yield y_O from Vincenzo et al. (2016).

From Eq. (15), and fixing ζ_a to 0, because of our assumption that inflowing gas is pristine, it can be seen that there are two variables that shape the Z_g − M_star–SFR–M_g hypersurface: Δη and ζ_w. With this formulation, it is impossible to determine separately the accretion and outflow mass-loading factors, η_a and η_w. Physically, this means that the important driving factor is the difference Δη, that is to say the dominance (or not) of gas accretion over bulk outflows. Below we describe the results from a Bayesian approach for inferring mass- and metallicity-loading factors in MAGMA.

4.1. Bayesian approach

For each galaxy with four observables (M_star, SFR, Z_g, M_g), we need to determine two parameters, Δη and ζ_w. However, only two observables enter into Eq. (15): Z_g and μ_g. Thus, the problem is severely underdetermined. To overcome this limitation, we have divided MAGMA into M_star bins, and assumed that Δη and ζ_w are the same for each mass bin. The M_star bins were defined to guarantee a sufficient number of galaxies in each bin, but also to adequately emphasize the lowest-M_star MAGMA galaxies.

To find the best-fit combination of (Δη, ζ_w) for each M_star bin, we first calculate χ² value for each of N galaxies in the bin⁴:

$$\begin{array}{c}\hfill {\chi }_{i,j}^2={\left(\frac{Z_{\mathrm{g}}^{\mathrm{obs}(j)}-Z_{\mathrm{g}}^{\mathrm{mod}(i)}}{\sigma }\right)}^2,\end{array}$$(19)

where $Z_{\mathrm{g}}^{\mathrm{obs}(j)}$ for the jth galaxy is given by Eq. (18) and $Z_{\mathrm{g}}^{\mathrm{mod}(i)}$ by Eqs. (15) and (16) for the ith parameter pair (Δη, ζ_w); we have assumed a constant error σ = 0.1, similar to the deviation of the MZR fit discussed in Paper I. χ² is computed over the grid of Δη and ζ_w parameters, in order to obtain N values (indexed by j in Eq. (19)) of χ² for each parameter pair (Δη, ζ_w).

Since only one value of χ² can be associated with each unique parameter pair, a decision must be made for the assignment of this single value for the N galaxies in the mass bin. This can sometimes be accomplished through stacking (e.g., Belfiore et al. 2019), but we cannot do that here, because of the way $Z_{\mathrm{g}}^{\mathrm{mod}}$ is defined (Eqs. (15) and (16)). We experimented with several methods; the approach that best captures the sensitivity of the χ² distribution to the variation in the parameters (Δη, ζ_w) is to assign the single χ² value to the mean of the lowest quartile in the χ² distribution of the N galaxies in the mass bin. This is a rather arbitrary choice but the numbers of galaxies in each M_star bin are sufficient to sample in a meaningful way the galaxies with the best fits. Thus for each mass bin, assigning this single χ² value to each parameter pair, we obtain the distribution of χ² values for that bin.

We then apply the standard Bayesian formulation:

$$\begin{array}{c}\hfill P(\mathit{\theta }|\mathit{D})\propto P(\mathit{\theta })P(\mathit{D}|\mathit{\theta })\end{array}$$(20)

to derive the full posterior probability distribution P(θ|D) of the parameter vector θ = (Δη, ζ_w), given the data vector D = (Z_g, μ_g). This posterior is proportional to the product of the prior P(θ) on all model parameters (the probability of a given model being obtained without knowledge of the data), and the likelihood P(D|θ) that the data are compatible with a model generated by a particular set of parameters. The data are assumed to be characterized by Gaussian uncertainties, so the likelihood of a given model is proportional to exp(−χ²/2). We assume uniform priors: Δη ∈ [ − 3.4, 3.4] in linearly sampled steps, and log₁₀(ζ_w) ∈ [ − 0.1, 2.5] in logarithmic steps of 0.05 dex.

The best-fit parameter vector (Δη, ζ_w) is evaluated by constructing the probability density functions (PDFs), weighting each model with the likelihood exp(χ²/2), and normalizing to ensure a total probability of unity. Specifically, this was achieved by marginalization over other parameters:

$$\begin{array}{c}\hfill P(\mathit{\theta }|\mathit{D})={\displaystyle \int }P(\mathit{\theta },Y|\mathit{D})\mathrm{d}Y\end{array}$$(21)

where Y is either Δη or ζ_w, excluding the parameter of interest (ζ_w or Δη, respectively).

The results for the different mass bins are shown in Fig. 1. The χ² surfaces show the relation between Δη and ζ_w. For each M_star bin, the PDF median is shown as a vertical dashed line, and the 1σ confidence intervals (CIs) by a shaded violet region. We have taken the maximum-likelihood estimate (MLE) to be the median of the PDF for ζ_w and the mode for Δη. This choice is discussed more fully below.

Fig. 1. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA binned by Mstar. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF median) is shown by a vertical dashed line.

Figure 1 shows that Δη is degenerate; low values of χ² can be achieved for a broad range of Δη. However, strikingly, the spread of allowed values for Δη at low M_star (in the two upper left-most panels in Fig. 1) is much broader than for higher stellar masses. In the higher M_star bins, the allowed distribution of Δη is truncated to increasingly small values of Δη. This implies that η_a, with equal probability, can be significantly larger than η_w at low masses, as also shown by the PDFs in each panel. Another clear result that emerges from Fig. 1 is that higher ζ_w occur at lower M_star.

Figure 1 also shows that while the PDF of ζ_w shows a form similar to a Gaussian distribution, and is thus fairly well constrained, the PDF for Δη is not. It peaks at a given value, but then falls off toward lower values gradually, leading to the conclusion that the best-fit value may depend on the interval over which the uniform prior is sampled. To investigate this, we have performed the Bayesian calculations using two different intervals in Δη:Δη ∈ [ − 10, 3.4] and Δη ∈ [0, 3.4], again in linearly sampled steps. These cornerplots are shown in Appendix A, and the best-fit values for Δη and ζ_w for Δη ∈ [0, 3.4] are given in the last two columns of Table 1.

Results are shown graphically in Fig. 2, where Δη is plotted against V_circ in the upper panel and ζ_w vs. V_circ in the lower one. Because the outflow strengths probably depend on the depth of the potential well, we assess them in terms of the virial velocities, V_circ, as proposed by Peeples & Shankar (2011). We have calculated virial velocities in several ways, but ultimately adopt V_out, velocities from the baryonic Tully–Fisher relation (BTFR) predicted by APOSTLE/EAGLE simulations from Sales et al. (2017), focused on the low-mass regime. Other alternatives are also viable, such as the power-law BTFR from Lelli et al. (2016b), or velocities obtained from the M_vir − M_star relation via halo abundance matching (Moster et al. 2010). However, at low M_star, and low baryonic masses, the inferred virial velocities change significantly according to whether the relation is a straight power law as in Lelli et al. (2016b) or inflected at low masses as in Sales et al. (2017). The simple power-law relation gives V_circ values that are apparently too low relative to the rest of the sample, while the V_circ values by Sales et al. (2017) are more consistent. Moreover, because the low-mass end of our sample is somewhat extreme, there is also some doubt about the applicability of the usual halo-abundance matching techniques. In any case, our conclusions are independent of the formulation used to translate the mass scale to the kinematic one.

Fig. 2. Difference of mass-loading factors Δη (top panel) and metal-loading factors ζw (bottom panel) as a function of virial velocities obtained from the BTFR from Sales et al. (2017). In both panels, filled symbols show the Bayesian results, and in the top panel, open ones illustrate the upper bound of the PDF 1σ confidence intervals (CI). Different symbol types correspond to the interval of Δη in the Bayesian priors, as described in the text; the symbols for Δη ∈ [0, 3.4] and Δη ∈ [ − 10, 3.4] are (arbitrarily) offset along the abscissa for better visibility. The gray regions give the ±1σ uncertainties in the best-fit values. Also shown is the “equilibrium” asymptote, namely Δη = α (see Sect. 3.1). In the top panel, the two estimates of Δη from different priors are exactly coincident. In the bottom panel, the formulations by Peeples & Shankar (2011) for two O/H calibrations, redone here for our value of y, are given by (purple) short-dashed and (blue) long-dashed curves for the PP04N2 and D02 calibrations, respectively. Also shown are the metal-loading factors computed by Chisholm et al. (2018) from COS data. The text gives more details.

As shown in Fig. 2, the different intervals used for the Bayesian priors give fairly consistent results both for Δη and for ζ_w. The best-fit Bayesian value for Δη is ≈α for all the prior intervals sampled, and there is a clear trend for ζ_w to increase toward lower V_circ (or stellar mass, as we show later). Over the mass range sampled by MAGMA, ζ_w is almost a factor of ten higher at lower M_star (V_circ) relative to the most massive galaxies.

Also shown in the bottom panel of Fig. 2 are the estimates of ζ_w obtained by Chisholm et al. (2018) with Cosmic Origins Spectrograph (COS) observations of ultraviolet absorption lines. By analyzing the underlying stellar continuum, and using the observed line profiles to constrain optical depth and covering fraction, Chisholm et al. (2018) were able to quantify the mass outflow rate, Ṁ_w. Then, through photoionization modeling, they obtain column densities of each element that are then used to constrain the metallicity of the outflow Z_w. As shown in Chisholm et al. (2018), their results are well fit by the analytical predictions of Peeples & Shankar (2011), with yields in the range 0.007−0.038, which agree with our reference value of y = 0.037.

4.2. Comparison with a closed-form approximation

Figure 2 also illustrates the predictions of Peeples & Shankar (2011) using their formulation for the MZR (taken from Kewley & Ellison 2008, and converted to a Chabrier 2003 IMF). Peeples & Shankar (2011) inferred ζ_w from the behavior of M_g and Z_g with M_star, using known derivatives derived empirically; we discuss this approach more fully below. Figure 2 shows their formulation relating ζ_w and V_circ:

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{w}}={\left(\frac{V_0}{V_{\mathrm{circ}}}\right)}^b+{\zeta }_0,\end{array}$$(22)

where, for consistency, we have recomputed their fits using our higher adopted yield (for y_O) of 0.037 taken from Vincenzo et al. (2016) for a Chabrier (2003) IMF. As shown in Fig. 2, the PP04N2 results are in good agreement with MAGMA, while the D02 calibration (Denicoló et al. 2002) is slightly discrepant, and more consistent with Chisholm et al. (2018), despite the higher yield. This illustrates how the computed ζ_w values depend on the adopted yields, and underscores the necessity of a self-consistent framework. We have performed our entire metallicity analysis with lower (and higher) yields than the ones used here, and the results do not change.

Peeples & Shankar (2011) assessed Δη and ζ_w based on known derivatives calculated from observed scaling relations of M_g and Z_g with M_star. In some sense, this is complementary to our Bayesian approach, but may also skew results as we show below. Following Pagel (2009), the time dependence in Eq. (9) can be eliminated by noting that Ṁ_⋆ = α ψ. Thus, by changing variables we obtain:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\mathrm{g}}}{\mathrm{d}{M}_{\star }}=\frac{{\dot{M}}_{\mathrm{g}}}{{\dot{M}}_{\star }}=\frac{{\dot{M}}_{\mathrm{a}}-{\dot{M}}_{\mathrm{w}}}{\alpha \psi }-1=\left(\frac{{\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}}{\alpha }\right)-1.\end{array}$$(23)

The same can be done for Eq. (10), and eliminating the time dependence as above gives:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_Z}{\mathrm{d}{M}_{\star }}=y+Z_{\mathrm{g}}(\frac{{\zeta }_{\mathrm{a}}-{\zeta }_{\mathrm{w}}}{\alpha }-1)·\end{array}$$(24)

Our observations do not directly constrain M_Z, but rather Z_g, so we need to formulate dZ_g/dM_⋆ as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_Z}{\mathrm{d}{M}_{\star }}=Z_{\mathrm{g}}\frac{\mathrm{d}{M}_{\mathrm{g}}}{\mathrm{d}{M}_{\star }}+M_{\mathrm{g}}\frac{\mathrm{d}{Z}_{\mathrm{g}}}{\mathrm{d}{M}_{\star }}·\end{array}$$(25)

The left-hand side is given by Eq. (24), and dM_g/dM_⋆ by Eq. (23), so we can solve for dZ_g/dM_⋆:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{Z}_{\mathrm{g}}}{\mathrm{d}{M}_{\star }}=\frac{1}{\alpha M_{\mathrm{g}}}[q+Z_{\mathrm{g}}({\zeta }_{\mathrm{a}}-{\zeta }_{\mathrm{w}}-{\eta }_{\mathrm{a}}+{\eta }_{\mathrm{w}})].\end{array}$$(26)

As they should be, Eqs. (23) and (26) are related to Eqs. (9) and (12) by the simple derivative, Ṁ_⋆ = α ψ. Equations (23) and (26) are slightly different from the ones derived by Peeples & Shankar (2011), because of the different assumptions they made for the stellar enrichment given by Eq. (5).

For a given M_g vs. M_star relation, we know dM_g/dM_star, the left side of Eq. (23); for a given MZR, we know dZ_g/dM_star, the left side of Eq. (26). Thus, in the spirit of Peeples & Shankar (2011), Δη = η_a − η_w and ζ_w can be directly constrained observationally. Inverting Eqs. (23) and (26) gives Δη and ζ_w explicitly:

$$\begin{array}{cc}& {\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}}=\alpha (F_{\mathrm{g}}\frac{\mathrm{d}log{M}_{\mathrm{g}}}{\mathrm{d}log{M}_{\star }}+1),\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\zeta }_{\mathrm{w}}=\frac{q}{Z_{\mathrm{g}}}-\alpha [F_{\mathrm{g}}(\frac{\mathrm{d}log{Z}_{\mathrm{g}}}{\mathrm{d}log{M}_{\star }}+\frac{\mathrm{d}log{M}_{\mathrm{g}}}{\mathrm{d}log{M}_{\star }})+1]\hfill \\ \hfill & =\frac{q}{Z_{\mathrm{g}}}-({\eta }_{\mathrm{a}}-{\eta }_{\mathrm{w}})-\alpha F_{\mathrm{g}}\frac{\mathrm{d}log{Z}_{\mathrm{g}}}{\mathrm{d}log{M}_{\star }}·\hfill \end{array}$$(28)

In Eq. (28), we have defined F_g ≡ M_g/M_star. Peeples & Shankar (2011) were the first to derive such a formulation for ζ_w; as mentioned in Sect. 3, the formalism here differs slightly from theirs because of the way they dealt with metal production in stellar ejecta.

The problem with Eq. (27) is that for positive d log M_g/d log M_star, Δη is by definition > 0, because F_g is ≥0. In Paper II, the robust least-squares power–law relation⁵ between M_g and M_star is found to be:

$$\begin{array}{c}\hfill log(M_{\mathrm{g}})=(0.72\pm 0.02)log(M_{\star })+(2.51\pm 0.22),\end{array}$$(29)

thus with d log M_g/d log M_star ≈ 0.7, roughly consistent with earlier works (e.g., Leroy et al. 2008). Because of this positive derivative, to infer Δη by relying on Eq. (27) and d log M_g/d log M_star would give results that are inherently skewed toward higher values of Δη. Qualitatively, we would expect that as a galaxy evolves, the stellar mass continues to grow, while the gas mass diminishes; such a scenario would give d log M_g/d log M_star < 0, and possibly Δη < α. Consequently, from a physical point of view, it may be inaccurate or unrealistic to invoke the sample’s overall d log M_g/d log M_star to calculate Δη. What we require is an estimate of Δη for individual galaxies which may or may not be related to the statistical properties of the parent sample at z ≈ 0.

In any case, Fig. 2 illustrates that ζ_w is not particularly sensitive to Δη. The predictions of ζ_w from Peeples & Shankar (2011) using d log M_g/d log M_star are fairly consistent with those we find with the Bayesian approach, even though the former requires Δη > 0 and the latter suggests that Δη ≈ 0.

5. Shaping the mass–metallicity relation with outflows

In the previous section, we estimated Δη and ζ_w through a Bayesian approach. Here we parameterize Δη and ζ_w as a function of both M_star and μ_g to predict quantitatively how mass loading Δη and metal loading ζ_w shape the MZR.

5.1. Stellar mass, metallicity, and outflows

First, we incorporate M_star as the defining variable. Figure 3 (left panel) shows the trends of Bayesian-inferred Δη and ζ_w with observed Z_g as a function of M_star (in logarithmic space). The bottom left panel shows the predictions for Z_g that are inferred from the functional forms we have found for Δη and ζ_w (see Table 2). The middle left panel shows ζ_w together with the robust cubic best-fit polynomials with coefficients as given in Table 2. The cubic polynomial has no physical significance; rather our objective is to parameterize Δη and ζ_w relative to M_star in order to compute Z_g.

Fig. 3. Difference of mass-loading factors Δη and metal-loading factors ζw inferred from the Bayesian analysis, and observed Zg in the bottom panel, as a function of Mstar on the left, and as a function of μg on the right. In the two upper panels, the best-fit Bayesian values are shown as filled symbols (circles or squares according to the Δη prior), with error bars corresponding to the ±1σ excursion, which is also illustrated by light-gray shaded regions. The rightmost point (magenta circle) in the top and middle right panels corresponds to the “extra” high-μg bin shown in Fig. 6. In the top panels, a horizontal dashed line illustrates the “equilibrium solution” with Δη = α. In the top left panel, the (purple) dashed curve is a linear power-law fit of Δη to log(Mstar), while the dashed curve in the top right panel shows Eq. (27), with d log Mg/d log Mstar = 0.72, as found in Paper II. The fits to the curves for ζw in the middle panels are cubic: log(ζw) as a function of log(Mstar) in the left panel, and ζw as a function of Fg in the right. In the lower panels, observed MAGMA Zg are shown as small filled circles, color coded by gas fraction, while the Zg medians and 1σ deviations are shown by large (black) filled circles and error bars. The dashed curves in the lower panels correspond to the calculation of Zg based on the fits in the upper panels (as a function of μg on the left, and log(Mstar) on the right): the long-dashed curves (dark orange) correspond to the Δη ≥ 0 priors, while the short-dashed ones (purple) to the “symmetric” priors. The green curves in the left panels are described in the text, obtained by fitting individual galaxy Δη and ζw calculated using Eqs. (27) and (28).

The slight differences (∼0.15 dex for M_star ≳ 10⁹ M_⊙, long-dashed orange- vs. short-dashed purple curves) in ζ_w according to the Δη prior are evident in the middle left panel in Fig. 3 and the lower one. When the constraint on the prior Δη > 0, the Bayesian inference of ζ_w gives lower values, and is a better approximation of observed Z_g (orange long-dashed curve). On the other hand, when the Δη priors include negative values ζ_w is larger (see also Fig. 2), but Z_g is underestimated (purple short-dashed curve). Since the curves shown in the lower panel are not fits to Z_g, but rather the parameterization of Z_g as a function of Δη and ζ_w, this may be telling us that positive values over negative Δη are preferred, at least in the current evolutionary states of the MAGMA galaxies.

The green curves in the top two left panels of Fig. 3 illustrate the polynomial fits (quadratic for Δη and cubic for log ζ_w) that would be obtained had we calculated Δη and ζ_w from Eqs. (27) and (28), respectively, using d log Z_g/d log M_⋆ and d log M_g/d log M_⋆ derived from the MAGMA scaling relations. The corresponding green curve in the lower left panel is the approximation of Z_g that would be obtained from these fits. It is clear that for log M_star/M_⊙ ≳ 8.5, Z_g is equally well estimated by the priors of Δη ≥ 0 and the fits to individual parameters computed with Eqs. (27) and (28). Wind metal loading ζ_w estimated with the symmetric prior of Δη falls short (∼0.15 dex) of observed Z_g in this mass range. In this mass range, Fig. 3 also shows that Δη is unimportant; the main factor that shapes Z_g is ζ_w.

On the other hand, the fit to Δη in the top panel (green curve, Eq. (27)) is significantly rising toward lower M_star, unlike the Bayesian values for Δη that, with both sets of priors, is roughly constant ≈α. The influence of Δη at low M_star (log M_star/M_⊙ ≲ 8.5) is clearly shown by the green curve in the lower left panel of Fig. 3. The fall-off of Z_g toward low M_star is less severe than would be predicted by a flat Δη, and is more consistent with the observations.

5.2. Gas fraction, metallicity, and outflows

In Paper II (see also Eq. (29)), we found that gas content F_g grows with decreasing stellar mass (see also e.g., Haynes et al. 2011; Huang et al. 2012; Boselli et al. 2014b; Saintonge et al. 2017; Catinella et al. 2018). Figure 4 shows this comparison explicitly for baryonic gas fraction μ_g and M_star in MAGMA galaxies. There is a clear trend for higher μ_g as M_star decreases, but in the very lowest M_star regime, below log(M_star/M_⊙) ≲ 8.5, the scatter is large. In MAGMA, most galaxies in this regime have μ_g ≳ 0.5, but there are two exceptions, Sextans A and WLM, at very low μ_g despite low global stellar mass. This apparent inconsistency probably arises from underestimates of the total gas content corresponding to the localized regions where CO was detected (e.g., Shi et al. 2015, 2016; Elmegreen et al. 2013; Rubio et al. 2015), and indicates the inherent difficulties of probing the gas, stellar, and metal contents in these extreme dwarf systems. Although the median trend for less-massive galaxies to be more gas rich is evident down to M_star ∼ 3 × 10⁸ M_⊙, the large overall scatter makes it virtually impossible to infer trends of Δη and ζ_w with μ_g from their trends with M_star.

Fig. 4. Baryonic gas fraction μg plotted against (log) Mstar for MAGMA galaxies. Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the gray region. As noted in the text, the two low-mass galaxies with low μg are Sextans A and WLM, with probable substantial underestimates of the total gas content.

Thus, to assess the dependence of Δη and ζ_w on μ_g, we cannot use the analysis of Sect. 4.1. Consequently, we repeated the Bayesian estimation by binning in μ_g, rather than in M_star. As before, bin boundaries were chosen such that a similar number of galaxies is found in each of the nine bins; the priors for Δη and ζ_w are as before (see Sect. 4.1). The results for the symmetric prior on Δη ∈ [ − 3.4, 3.4] are shown in Fig. 5, where the layout of the figure is the same as in Fig. 1. Here there is a trend for larger Δη excursions at high μ_g (see lower panels), similar to what we found for low M_star. There is also a tendency for the distribution of ζ_w to be broader at high μ_g, that is also reflected in the broader range of Δη.

Fig. 5. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA binned by μg. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line.

At high μ_g, there is some evidence for a bimodal distribution of ζ_w, as shown in Fig. 6, where we have raised the boundary of the highest μ_g bin. In this μ_g bin, the rightmost one plotted in Fig. 5, the best-fit (mode) Δη is higher than any. Nevertheless, as can be seen in the lower-right panel, higher values of ζ_w are possible for virtually every Δη sampled, so that the dependence of ζ_w on Δη is apparently degenerate, as we concluded in Sect. 4.1.

Fig. 6. Corner plot of χ2 surface as a function of the model parameters (Δη, ζw) for the highest μg bin MAGMA, here with a higher boundary than in Fig. 5: μg > 0.75. As in previous figures, the violet contours correspond to the minimum χ2 value, and the top and right panels report the probability density distributions for the marginalized parameters. Confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line.

To parameterize Δη and ζ_w, here we adopt F_g as the independent variable because it is statistically better behaved than μ_g, as a result of the simpler dependence on M_star. However, it is more convenient to compare quantities to μ_g as shown in the right panel of Fig. 3 where Bayesian estimates for Δη, ζ_w, and observed Z_g are plotted as a function of μ_g. The curve in the top right panel is not a fit, but rather the prediction of Eq. (27) using d log M_g/d log M_star = 0.72 from Eq. (29). The middle right panel shows ζ_w together with the robust cubic best-fit polynomial with coefficients as given in Table 2. Metallicity vs. μ_g is given in the bottom right panel, where the long-dashed curve corresponds not to a fit, but rather to the prediction of Eq. (15) with the functional forms for ζ_w (cubic polynomial) and Δη. The difference in ζ_w for the two different priors for Δη in the Bayesian analysis are analogous to those for M_star (left panel Fig. 3). At low μ_g, Z_g is better approximated by ζ_w inferred from Δη > 0 relative to ζ_w from the symmetric prior.

The influences of ζ_w and Δη are clearly seen in the trend for Z_g with μ_g. The steep increase of ζ_w and Δη toward high μ_g is reflected in the dramatic reduction of metallicity toward high μ_g. The shape of the MAGMA Z_g curve relative to μ_g resembles that found by Erb (2008, see her Fig. 3), for z ∼ 2 galaxies, but with some important differences. The first is that Erb (2008) assume that Z_w/Z_g = 1, that is to say that ζ_w = η_w. Appendix C presents the trivial derivation of this result. Under this assumption, Z_g in Eq. (15) is essentially governed by η_a, with η_w entering only into the M_g/M_i exponent (see also Appendix C). Consequently, there is less variation at low μ_g, with the prediction by Erb (2008) assuming a constant asymptotic value for μ_g as high as ∼0.35. By considering the change in Z_w relative to Z_g, as we do here, Z_g never reaches a truly asymptotic value, although there is flattening toward lower μ_g.

The second difference is the assumption for R (thus α) in the metallicity treatment of Erb (2008). She uses the same metallicity calibration as adopted here (PP04N2), but assumes that α, the fraction of gas retained in stars, is unity, thus neglecting the gas that is expelled in SNe and stellar winds. This would imply a difference −0.34 dex in her and our absolute metallicity scale. An additional factor is the different yield adopted here (y = 0.037) relative to y = 0.019 used by Erb (2008). Considering these two numerical differences, the absolute scale in Erb (2008) would be reduced by ∼0.6 dex. Such a difference is consistent with the galaxies shown by Erb (2008), considering that her sample is at z ∼ 2, while MAGMA is at z ∼ 0. Thus, as observed, lower Z_g (∼ − 2.7 dex) at low gas fractions would be expected for her sample relative to ours (e.g., Erb et al. 2006; Maiolino et al. 2008).

5.3. Considerations

Interpretation of our results is not altogether straightforward. On the one hand, the Bayesian analysis shows that Δη as a function of M_star tends toward α which would indicate a condition of “gas equilibrium” (e.g., Davé et al. 2012; Mitra et al. 2015). On the other, Δη as a function of μ_g does not show this behavior, but rather resembles what would be predicted by Eq. (27). This distinction stems from the highly imperfect correlation of μ_g with M_star as shown in Fig. 4. Figures 1 and 5 also illustrate that a wide range of Δη is possible (Δη < 0, Δη < α, Δη ≥ α), but that this only weakly affects ζ_w.

Our analysis also shows that within a 1σ range, Δη can be < α, over all mass ranges, but possibly more importantly, Δη can be significantly > α toward low M_star and high μ_g (see Fig. 3). We take this to imply that at low M_star and high μ_g there is more possibility for increased overall gas accretion (η_a) relative to wind outflows (η_w). This could be related to the stochastic nature of star-formation histories in low-mass galaxies (e.g., Weisz et al. 2012; Madau et al. 2014; Krumholz et al. 2015; Emami et al. 2019), and to the availability of gas.

Figure 3 shows that Δη ≈ α almost independently of M_star, thus implying gas equilibrium at all masses; however, Δη increases as a function of μ_g roughly as predicted by Eq. (27). Because we would expect that younger, less evolved, galaxies contain more gas, thus higher μ_g, this would seem to suggest that Δη could depend on the evolutionary phase of a galaxy. Since we cannot know the star-formation histories (SFHs) of the individual galaxies, this suggestion is difficult to test. However, we can fix a simple SFH, such as a falling or rising exponential, and experiment with what such a “toy model” would predict for d log M_g/d log M_star as a function of time. If we consider that the timescale of star formation τ_SF is tightly linked to the gas depletion time τ_g (${\tau }_{\mathrm{g}}={{ϵ}_{\star }}^{-1}$), we can explicitly calculate the growth of M_star and M_g using Eqs. (2) and (14), respectively.

These toy-model calculations are shown in Appendix B. As expected, d log M_g/d log M_star depends clearly on Δη; both positive and negative values of d log M_g/d log M_star are possible depending on the SFH, the age of the galaxy t_gal, the timescale for star formation τ_SF, and most importantly on Δη. The value of Δη, relative to α in the exponent of Eq. (14), determines whether gas mass is increasing or diminishing with time. Since stellar mass is monotonically increasing, the value of Δη sets the value of d log M_g/d log M_star. Nevertheless, our results show that Δη is almost impossible to constrain, demonstrating behavior that varies according to whether the independent variable is M_star or gas content μ_g.

As discussed in Sect. 3.2, when Δη = α, as for the “gas-equilibrium” scenario, Eq. (15) is numerically intractable because of the zero denominator in the exponent for M_g/M_i. This same condition gives Ṁ_g = 0 and M_g = M_i (see Eq. (16)). Because of the relatively long total-gas depletion times, ∼a few Gyr, even without a true equilibrium solution, essentially M_g/M_i remains unity throughout most of the lifetime of the galaxy (see Eq. (14)).

An important aspect of our formalism is that, in essence, through Eq. (13), we assume that the SFH corresponds also to the time evolution of gas mass. Moreover, because we assume that total system mass is conserved and that the time derivative of accreted gas mass Ṁ_a is proportional to SFR through η_a, it is possible that the formalism itself requires a rough equality between Δη and α, the lock-up fraction.

Despite the degeneracy of Δη and the inability of our formalism to determine unambiguously its behavior, ζ_w is fairly well determined (see Figs. 2 and 3). After discussing effective yields in Sect. 6, in the remainder of the paper, we focus on ζ_w, and what can be learned from our results about metal- and mass-loading in the stellar-driven outflows of star-forming galaxies.

6. Effective yields

Above we have quantified a scenario for the relation between Z_g and gas content in galaxies that requires significant gas accretion and stellar-driven outflows. In this section, we recast our results in the context of the effective yield, as measured against the predictions from a different, “closed-box”, paradigm. A galaxy evolving as a closed box obeys a simple relationship between the metallicity and the gas mass fraction (Searle & Sargent 1972; Pagel & Patchett 1975):

$$\begin{array}{c}\hfill Z_{\mathrm{g}}=yln(1/{\mu }_{\mathrm{g}})\end{array}$$(30)

where μ_g is the baryonic gas mass fraction and y is defined as the “true” yield (Sect. 3.1). As gas is converted into stars, metal content grows and the gas mass diminishes. If the galaxy evolves as a closed box, then the yield y in Eq. (30) should be equal to the nucleosynthetic yield. However, galaxies do not evolve as closed systems, so we expect the overall metal yield to differ from the nucleosynthetic value. The effective yield y_eff can be defined by inverting Eq. (30) (e.g., Lequeux et al. 1979):

$$\begin{array}{c}\hfill y_{\mathrm{eff}}=\frac{Z_{\mathrm{g}}}{ln(1/{\mu }_{\mathrm{g}})}·\end{array}$$(31)

The effective yield is expected to be lower than the nucleosynthetic value y because metals in galaxies are lost from the system through outflows, and the metal content of the ISM gas will be diluted with infall of metal-poor gas (e.g., Edmunds 1990).

The effective yields y_eff for MAGMA, computed using Eq. (31), are shown in Fig. 7, plotted against baryonic gas mass fraction μ_g. The model prediction from Eq. (15) through our parameterizations of Δη and ζ_w is shown as two curves, corresponding to the two different Δη priors. For reference, our adopted value of y is shown as a horizontal short-dashed line. As also evident in previous figures, the results from the Δη ≥ 0 prior (long-dashed orange curve) are in slightly better agreement with the data. The inflection of y_eff at high μ_g is due mainly to the strong downward curvature of Z_g as seen in Fig. 3. Moreover, y_eff is always (with a single possible exception at high μ_g, NGC 4731) significantly smaller than the true yield y = 0.037 adopted here. Edmunds (1990) has shown that this behavior for y_eff would be expected if metals were lost through stellar- and SNe-driven bulk outflows, or diluted by accretion of metal-poor gas. MAGMA shows a strong confirmation of this scenario. It is clear from Fig. 7 that our parameterized model is a good approximation to the data (see also Fig. 3).

Fig. 7. Effective yield as a function of μg. Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the surrounding gray regions. The two dashed curves correspond to the predictions from the Δη and ζw parameterizations (see also Fig. 3), with the symmetric Δη priors given as a short-dashed (purple) curve, and the Δη ≥ 0 priors as a long-dashed (orange) one. The horizontal dashed line corresponds to our adopted true yield y = 0.037 (see text).

Models that do not take into account the metal-enriched bulk outflows at low M_star, namely assuming that Z_w = Z_g, predict that at high μ_g, y_eff would approach the nucleosynthetic value y (e.g., Erb 2008). This can be seen by expanding Eq. (31) in a Taylor’s series (see Dalcanton 2007; Filho et al. 2013) so that for high μ_g, we have:

$$\begin{array}{c}\hfill y_{\mathrm{eff}}\approx \frac{M_Z}{M_{\star }}·\end{array}$$(32)

Thus, in very gas-rich galaxies, the effective yield would be approximately independent of gas fraction. As pointed out by Dalcanton (2007), metal-poor accreted gas will lower the metallicity of the system, but it will also augment μ_g, thus maintaining the system in a pseudo-closed-box equilibrium. The linear dependence of y_eff on metal mass M_Z, as shown by Eq. (32), also means that metal-enriched outflows are very effective at lowering the effective yield in galaxies with high gas-mass fractions (see Dalcanton 2007). Indeed, Fig. 7 shows that even at high μ_g, y_eff is significantly lower than the true yield y.

It is important to remember that Eq. (31) measures Z_g in the ionized gas, but the gas content μ_g is derived from neutral gas (H I + H₂). There is some evidence that the metal contents of these two gas phases are different; in UV absorption studies of metal-poor dwarf galaxies (e.g., Thuan et al. 2005; Lebouteiller et al. 2013), the metallicity of the neutral gas in these systems was always found to be lower than in the ionized gas. A similar conclusion was reached by Filho et al. (2013) and Thuan et al. (2016) by comparing y_eff with y, and attributing the difference to the lower metallicity in the neutral gas.

Figure 8 plots the effective yield y_eff against stellar mass M_star (left panel) and circular velocity, V_circ (right). Even though there is much scatter at low M_star and V_circ (M_star ≲ 3 × 10⁸ M_⊙ and V_circ ≲ 50 km s⁻¹), y_eff is increasing with mass and velocity up to 3 × 10⁹ M_⊙ and ∼100 km s⁻¹, in qualitative agreement with similar trends reported in the literature (e.g., Garnett 2002; Pilyugin et al. 2004; Tremonti et al. 2004; Dalcanton 2007; Ekta & Chengalur 2010). For more massive galaxies, there is an inflection in the relations with a mild decline at masses and velocities larger than 3 × 10⁹ M_⊙ and ∼100 km s⁻¹.

Fig. 8. Effective yield as a function of Mstar (left panel) and Vcirc (right panel). Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the surrounding gray regions. Also shown in the left panel as two dashed curves are the predictions from the Δη and ζw parameterizations (see also Fig. 3), with the symmetric Δη priors given as a short-dashed (purple) curve, and the Δη ≥ 0 priors as a long-dashed (orange) one. The right panel illustrates galaxies with high star-formation efficiency (short gas depletion times) with a blue dashed curve, and those with low efficiency (long depletion times) as a red one. The horizontal short-dashed line corresponds to the yield of 0.037 adopted in this paper. See text for more details.

The right panel of Fig. 8 also shows the MAGMA sample divided into two regimes in gas depletion times, τ_g, defined as M_g/SFR (see Paper II); galaxies with long τ_g (τ_g > 4 × 10⁹ yr) are shown as a dashed red curve, and short times τ_g < 4 × 10⁹ yr as a dashed blue one. The gas depletion time τ_g is the inverse of the so-called star-formation efficiency (SFE, see Paper II), that also enters into Eq. (13), so short τ_g would correspond to higher SFE, and vice versa. There can be as much as a factor of ten in y_eff between these two categories of galaxies, with a mean difference of ∼2 − 3 in y_eff. At fixed mass, larger yields correspond to longer τ_g or lower SFE. In these galaxies with longer τ_g, the peak in y_eff at M_star ∼ 3 × 10⁹ M_⊙ seems to disappear, while there is a steady decline of y_eff toward lower M_star for shorter depletion times (see also Lara-López et al. 2019). Longer τ_g are also associated with higher gas fractions (see color coding in Fig. 8); thus Eq. (32) may be telling us that the higher effective yields are simply a function of the enhanced impact of metal-enriched outflows in gas-rich galaxies.

7. Characterizing stellar feedback through winds

By adopting the formalism for galactic chemical evolution formulated by Pagel (2009), we have quantified with MAGMA the metal-weighted mass loading factor ζ_w (=η_wZ_w/Z_g), and with somewhat less success, the difference in mass loading of gas accretion and outflows, Δη. However, we are interested in the individual mass-loading factor relative to the SFR, η_w, and the metal content of the winds, Z_w. In this section, we explore the constraints for mass loading η_w and Z_w provided by our results.

7.1. The entrainment fraction

Assuming that Z_w = Z_g, as done in much previous work (e.g., Erb 2008; Finlator & Davé 2008; Dayal et al. 2013; Hunt et al. 2016; Belfiore et al. 2019), implies that virtually all of the material in the outflow is ISM gas. The gas metallicity could be only a lower limit to the metallicity of the wind (e.g., Peeples & Shankar 2011), because the metals produced by SNe would be expected to enrich the gas above that of the galaxy (e.g., Vader 1986). Thus, to estimate Z_w, we also need to consider the fraction of gas that has been entrained in the wind, f_ent, and the metallicity of the SNe ejecta, Z_eject (see Dalcanton 2007; Peeples & Shankar 2011):

$$\begin{array}{cc}& Z_{\mathrm{w}}=f_{\mathrm{ent}}{Z}_{\mathrm{g}}+(1-f_{\mathrm{ent}})Z_{\mathrm{eject}}\hfill \\ \hfill & \frac{Z_{\mathrm{w}}}{Z_{\mathrm{g}}}=f_{\mathrm{ent}}+(1-f_{\mathrm{ent}})\frac{Z_{\mathrm{eject}}}{Z_{\mathrm{g}}}·\hfill \end{array}$$(33)

In low-mass galaxies, f_ent is found to be low, both observationally (e.g., Martin et al. 2002) and from simulations (e.g., Mac Low & Ferrara 1999). The entrainment fraction is closely related to the mass-loading factor of the wind η_w, and probably varies with galaxy mass (gravitational potential, see e.g., Hopkins et al. 2012). The maximum oxygen abundance of SNe ejecta Z_eject is known to within a factor of 2−3, Z_eject ∼ 0.04 − 0.1 (e.g., Woosley & Weaver 1995), and can constrain f_ent and η_w (e.g., Martin et al. 2002). The increase of ζ_w with decreasing virial velocity and stellar mass (see Fig. 2) is possibly due to diminishing f_ent to raise Z_w closer to Z_eject. To verify this, we need to isolate η_w, the mass-loading factor in the product ζ_w = η_wZ_w/Z_g.

7.2. The driving mechanisms behind stellar outflows

One way to do this is to compare the V_circ scaling of ζ_w with expected V_circ scalings of η_w (see Eq. (22)). Galactic outflows powered by stellar feedback are thought to be driven either by momentum, imparted to the ISM by massive stellar winds and SNe through radiation pressure (e.g., Murray et al. 2005, 2010; Hopkins et al. 2011; Zhang & Thompson 2012), or by energy, injected into the ISM by massive stars and core-collapse SNe (e.g., Dekel & Silk 1986; Heckman et al. 1990; Dalla Vecchia & Schaye 2012; Christensen et al. 2016; Naab & Ostriker 2017).

For outflows that conserve momentum, we would expect a power-law scaling ${\eta }_{\mathrm{w}}\propto V_{\mathrm{circ}}^{-1}$ while for energy-conserving winds, ${\eta }_{\mathrm{w}}\propto V_{\mathrm{circ}}^{-2}$ (e.g., Dekel & Silk 1986; Murray et al. 2005; Hopkins et al. 2012). The only way such scalings can be reconciled with the (modified power-law) form of Eq. (22) is by a complementary scaling of Z_w. Examples of such complementary scaling are shown in Fig. 9 where the ratio of wind and ISM metallicities Z_w and Z_g is plotted against virial velocity V_circ. In the left panel of Fig. 9, we have taken the normalization of Z_w/Z_g to be the same as that as the best-fit of Eq. (22) using our Bayesian estimates, simply V_circ/V₀, where V₀ = 32 km s⁻¹ for symmetric Δη priors, and V₀ = 35 km s⁻¹ for priors with Δη ≥ 0. Both fits have power-law index b ∼ 2. The assumptions of momentum-conserving or energy-conserving outflows lead to different estimates for the trend of Z_w/Z_g.

Fig. 9. Ratio of inferred wind and ISM metallicities Zw/Zg plotted against virial velocity Vcirc. Also shown (with arbitrary scaling) is the shape of Bayesian ζw as also plotted in Fig. 2; it is clearly not a power law, thus the different power-law forms of ηw introduce inflections in Zw/Zg. Left panel: as indicated in the legend, different trends of Zw/Zg are inferred by dividing the fitted form of ζw (see Eq. (22)) by various scalings of ηw (see text for more details). The shaded regions correspond to the trends allowed by the two different Bayesian estimates of ζw according to the different Δη priors. Right panel: broken-power law trend for ηw found by Muratov et al. (2015) is shown as a solid heavy (gray) curve, where the inflection at Vcirc = 60 km s−1 corresponds to a sharp decrease of Zw/Zg (see text for more details). The short-dashed (dark-blue) curve shows the relation found by Mitchell et al. (2020), in the low-mass regime where stellar feedback dominates, and the long-dashed (blue) curve shows the observed relation derived by Chisholm et al. (2018).

Because dwarf galaxies are found to have a low entrainment fraction (e.g., Mac Low & Ferrara 1999; Martin et al. 2002), we would expect Z_w/Z_g to increase toward low V_circ because low f_ent would enrich the outflows toward the higher super-Solar abundance of the SNe ejecta, Z_eject (see Eq. (33)). Energy-conserving winds result in a flattening of Z_w/Z_g toward low V_circ (blue curves), and a cubic trend of η_w with V_circ gives a strongly declining Z_w/Z_g (gray curves). The trend of ${\eta }_{\mathrm{w}}\propto V_{\mathrm{circ}}^{-1}$ is the only one that gives a rising trend of Z_w/Z_g at low V_circ. Thus, momentum-conserving winds would be favored by our results.

Quantifying this conclusion is, however, dependent on the metallicity calibration. For the PP04N2 calibration used here, the roughly quadratic low-mass V_circ dependence is shallower than the cubic (or quartic) dependence found by Peeples & Shankar (2011) for different O/H calibrations. Other calibrations such as those by Tremonti et al. (2004) or Denicoló et al. (2002) give differing trends (with power-law slopes b > 3 at low V_circ), that are inconclusive for determining the nature of the mass-loading factor η_w. The PP04N2 calibration adopted here results in the shallowest increase of ζ_w toward low V_circ (power-law slope b ≈ 2), relative to other O/H calibrations, and this is reflected in a more incisive low-mass (low V_circ) increase in Z_w/Z_g.

The right panel of Fig. 9 shows the predictions for η_w from the sub-grid FIRE (Feedback in Realistic Environments) simulations (Muratov et al. 2015), and the EAGLE simulations (Mitchell et al. 2020). At low V_circ (< 60 km s⁻¹), Muratov et al. (2015) find the power-law dependence on V_circ to be steep (∼ − 3.2), while at higher masses (V_circ) the winds are momentum-driven with η_w ∝ $V_{\mathrm{circ}}^{-1}$. A similar conclusion was reached by Davé et al. (2013), namely that at low masses, winds are energy-driven (or even more steeply with V_circ), while momentum-driven winds dominate at higher masses. Judging from Fig. 9, these trends would be incompatible with our results. In contrast, EAGLE simulations (e.g., Mitchell et al. 2020) indicate a slope midway between, namely ∼ − 1.5, leading to a very gradual increase in metallicities toward lower V_circ. Comparing the left- and right-hand sides of Fig. 9 shows that such a dependence better reflects our results at low V_circ. At higher V_circ (masses), the situation is more complex because the flattening of ζ_w is difficult to interpret in terms of a single power-law dependence of η_w on V_circ. In any case, we conclude that at low masses, winds in stellar outflows cannot be energy driven, but rather are more consistent with momentum-driven outflows ($V_{\mathrm{circ}}^{-1}$), or a slightly steeper “combined” slope as found by the EAGLE models with $V_{\mathrm{circ}}^{-1.5}$.

There is also the question, not touched upon here, of whether or not outflows are able to expel material from the galaxy halo, thus enriching the IGM. The entrainment fraction f_ent would be key to relating the strength of the outflows with metallicity in the winds (e.g., Martin et al. 2002; Chisholm et al. 2018) and deserves further study. In addition, entrainment is also affected by the gas phases involved in the outflow. Our analysis traces the metallicity in the ionized gas (considered the “cool” phase by most simulations, T < 2 × 10⁴ K), thought to carry most of the mass, whereas physically realistic winds also contain a “hot phase” (T > 5 × 10⁵ K) that harbors most of the energy (Kim et al. 2020). Consequently, with MAGMA, we are missing the diagnostic of the hot X-ray phase.

8. Summary and future perspectives

In this paper, we explore the impact of ejective feedback on galaxy evolution using data from the MAGMA sample of 392 galaxies with H I and CO detections, which collects a homogeneous sample of local field star-forming galaxies, covering the widest range of stellar and gas masses, metalliticies and star formation compiled up to now. After presenting the sample, and the correlations among metallicity, SFR, M_star and (cool) gas-phase components in Paper I and Paper II, in this third paper of the series, we investigate outflow metal loading factors and effective yields in terms of gas fraction, mass and circular velocity, using the formalism of Pagel (2009). The main results of the paper follow.

• We model the chemical enrichment of MAGMA galaxies using a formalism based on Pagel (2009), without any constraint on the M_g time evolution and including the assumption that wind metallicity Z_w is different from the ISM metallicity Z_g. Wind metallicity Z_w can be greater than the ambient Z_g if it is primarily comprised of supernova ejecta or depressed if a sufficient amount of metal-poor gas is entrained as the wind propagates out of the galaxy (Chisholm et al. 2015; Creasey et al. 2015; Muratov et al. 2017; Pillepich et al. 2018). Using this formalism we have observationally constrained two parameters important for stellar feedback: (a) the metal-loading factor of the winds ζ_w; and (b) the difference of the mass loading of the accretion and winds, Δη = η_a − η_w. This is done through a Bayesian approach applied to subsamples of galaxies in different M_star and μ_g bins.

• In terms of trends with M_star, we find that Δη is degenerate, and exceedingly difficult to constrain. However, our Bayesian estimates based on μ_g bins show that Δη increases significantly with increasing gas fraction, implying that Δη is a strong function of gas content and possibly evolutionary state (see Appendix B).

• The wind metal-loading factors ζ_w are found to increase with decreasing stellar mass and circular velocity, with a flattening toward larger M_star. This indicates that the efficiency of metal-loading in the winds depends on the potential depth and is strongly enhanced in low-mass systems, in agreement with previous work based on different methodologies for estimating ζ_w (e.g., Peeples & Shankar 2011; Chisholm et al. 2018).

• We have also constrained effective yields y_eff with MAGMA data, and find that y_eff increases with stellar mass and V_circ up to a certain threshold, followed by a flattening toward larger M_star, and possibly a mild inversion at the highest masses.

• Finally, we have inferred possible trends with V_circ for the wind mass-loading factors η_w and find that our results for low V_circ (M_star) favor momentum-driven winds, rather than energy-driven ones.

Various quantities in our analysis, ζ_w, y_eff, and Z_g, present inflections at M_star ∼ 3 × 10⁹ M_⊙, corresponding to virial velocities ∼100 km s⁻¹. This corresponds to a “gas-richness” mass scale, associated with transitions of physical processes, separating two different mass regimes (Paper II; see also Kannappan 2004). It delineates an accretion-dominated regime at low masses, where stellar outflows are relevant, and an intermediate mass, gas-equilibrium regime, where outflows are less effective, and galaxies are able to consume gas through SF as fast as it is accreted. There would also be an additional, larger “bimodality” mass threshold at M_star ∼ 3 × 10¹⁰ M_⊙ (Paper II, Kannappan et al. 2013) but this is only weakly evident in MAGMA, mainly because of the exclusion of passive or “quenched” gas-poor galaxies and AGN.

These results can also be interpreted within a broader context, since the mass thresholds defining metallicity and gas properties are associated with changes of correlations involving stellar populations and stellar and dark matter distributions (Tortora et al. 2010, 2019). In particular, inflections with mass of metal-loading and metal yields at the “gas-richness” threshold seem to define changes of the slopes of total mass density profiles including dark matter (e.g., Tortora et al. 2019); this threshold approximatively separates galaxies with mass density slopes shallower and steeper than ∼ − 1 (see their Fig. 1). More investigations on this mass scale will follow.

As outlined in Paper II, further improvement of the MAGMA coverage, with more detections of total gas content in dwarf galaxies, will strengthen the conclusions of our series of papers. In addition, it is important to reassess the impact on our conclusions of assumptions made here, namely that the accretion rate Ṁ_a is proportional to the SFR and the chemically unenriched nature of the accreted gas. In this paper, we have also pursued a first step to combine information from completely independent datasets and observables. Only by including gas, stellar, kinematical and dynamical observables, will we be able to determine the efficiency of the physical phenomena which shape and regulate galaxy evolution.

Acknowledgments

We thank the anonymous referee and the editor for profound insight that helped to greatly improve the paper. We are grateful for funding from the INAF PRIN-SKA 2017 program 1.05.01.88.04.

References

Appendix A: Additional cornerplots for Bayesian best fits

Here we show the additional corner plots of the χ² surface as a function of the model priors (Δη, ζ_w) for MAGMA, for Δη ∈ [ − 10, 3.4] (Fig. A.1) and Δη ∈ [0, 3.4] (Fig. A.2).

Fig. A.1. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line. Here, Δη ∈ [ − 10, 3.4] (see Fig. 1 for the “symmetric” Δη intervals).

Fig. A.2. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF median) is shown by a vertical dashed line. Here, Δη ∈ [0,3.4] (see Fig. 1 for the “symmetric” Δη intervals).

Appendix B: SFHs and toy-model predictions for gas and stellar growth

In order to investigate the time derivatives of M_g and M_star and their ratio, we introduce an arbitrary SFH. As a simple example, we have taken an exponential with arbitrary scaling defined by β and the star-formation time scale τ_SF:

$$\begin{array}{c}\hfill \psi (t)={\psi }_{0}exp\left(\frac{\beta t}{{\tau }_{\mathrm{SF}}}\right).\end{array}$$(B.1)

Equation (B.1) can be integrated to obtain M_star [see Eq. (2)], at the current age of the galaxy t_gal:

$$\begin{array}{c}\hfill M_{\star }(t=t_{\mathrm{gal}})=\frac{{\psi }_0\alpha {\tau }_{\mathrm{SF}}}{\beta }[exp\left(\frac{\beta t_{\mathrm{gal}}}{{\tau }_{\mathrm{SF}}}\right)-1].\end{array}$$(B.2)

If we take the time-scale of star formation τ_SF to be equal to the gas depletion time τ_g [the inverse of ϵ_⋆, see Eq. (13)], there are two unknowns in Eq. (B.2): t_gal and ψ₀ (β is fixed arbitrarily for rising and declining SFHs). However, ψ(t_gal) = ψ₀ exp(βt_gal/τ_SF) = SFR and M_⋆(t_gal) = M_star, are the observables that we measure, so that Eq. (B.2) can be solved for ψ₀:

$$\begin{array}{c}\hfill {\psi }_0=\frac{\beta M_{\mathrm{star}}}{\alpha {\tau }_{\mathrm{SF}}}+\mathrm{SFR}.\end{array}$$(B.3)

This toy model also gives the initial gas mass M_i = ψ₀τ_SF.

The remaining variable needed to define the galaxy’s evolutionary state is t_gal, which can be established with ψ₀ through Eq. (B.1):

$$\begin{array}{c}\hfill t_{\mathrm{gal}}=\frac{{\tau }_{\mathrm{SF}}}{\beta }ln\left(\frac{\mathrm{SFR}}{{\psi }_0}\right)(\beta \ne 0).\end{array}$$(B.4)

In case of a constant SFH (β = 0), with ψ(t) = ψ₀ = SFR, the solutions are quite simple: M_⋆(t) = M_star, M_gas = M_i = SFR τ_SF, and Δη = α.

Figure B.1 shows the values inferred for t_gal with the simple SFH of Eq. (B.1) with β = −1 (declining exponential) plotted against M_star (top panel) and ψ₀ plotted against M_star (middle). Although there is significant scatter, galaxies with M_star ≲ 10⁹ M_⊙ are apparently younger, with ages as low as 10⁹ yr, rather than the overall inferred age of more massive galaxies, 10¹⁰ yr. The initial SFR ψ₀ is very tightly linked with M_star, over a dynamic range of ∼10³. The bottom panel of Fig. B.1 illustrates that the total gas depletion times (i.e., including both molecular and atomic components) are relatively constant with M_star.

Fig. B.1. Inferred galaxy age tgal plotted against Mstar for MAGMA galaxies (top panel); initial SFR at t = 0, ψ0 plotted against Mstar (middle); τg (assumed equal to τSF) plotted against Mstar (bottom). In the two top panels, we assumed a simple declining exponential SFH with β = −1 [Eq. (B.1)], and that ${\tau }_{\mathrm{SF}}={\tau }_g={{ϵ}_{\star }}^{-1}$.

Now the only missing ingredient to describe M_g is Δη [see Eq. (14)]. Ultimately, the toy model presented here requires through Eq. (13) that the SFH given by Eq. (B.1) is related to the time evolution of M_g [Eq. (14)]. If ${{\tau }_{\mathrm{SF}}}^{-1}$ is the same as ϵ_⋆, as we have assumed, the formalism requires:

$$\begin{array}{c}\hfill \mathrm{\Delta }\eta =\alpha +\beta .\end{array}$$(B.5)

Thus, fixing β allows only certain values of Δη. Figure B.2 shows the time evolution of M_g and M_star over the last 50% of the galaxies’ lifetime, assuming that t_gal is given by Eq. (B.4), and ψ₀ by Eq. (B.3). The different sets of curves correspond to different M_star medians and the corresponding medians of τ_SF, t_gal, and ψ₀, while the vertical separation is achieved through the assignment of different Δη as required by β = ( − 1, 0, 0.1). The rising exponential SFH cannot have β arbitrarily large because otherwise Eq. (B.3) has negative values; β = 0.1 is close to the largest value allowed by the data. Figure B.2 clearly shows that according to the SFH, and the value of β and consequently Δη, dlog M_g/dlog M_star can be positive or negative. Thus, it is difficult to assign Δη through the application of Eq. (27).

Fig. B.2. Time evolution of Mg vs Mstar over the last 50% of the galaxies’ lifetime, according to the toy-model based on three exponential SFHs: β = ( − 1, 0, 1), and also assuming that the SFR timescale τSF is the same as the gas depletion time, τg = ${{ϵ}_{\star }}^{-1}$. Galaxy age tgal is given by Eq. (B.4), ψ0 by Eq. (B.3), and the Mg evolution by Eq. (14). The median curves in different Mstar bins are plotted, where the (median observed) gas mass Mg is taken to be Mg(tgal). For Δη = α, the equilibrium solution, Mg = Mi throughout the lifetime of the galaxy.

Appendix C: The formalism for the “homogeneous wind”

In this paper we have determined solutions for Z_g in the most general case of Z_w ≠ Z_g. We demonstrate (as already shown by Peeples & Shankar 2011) that in general this is a valid assumption since the outflowing gas Z_w is not the same as the ISM gas, Z_g, especially at low masses (e.g., Chisholm et al. 2015; Creasey et al. 2015; Muratov et al. 2017).

However, much earlier work assumed, instead, that Z_w = Z_g, adopting the so-called homogeneous wind model (e.g., Pagel 2009; Erb 2008). This assumption leads to ζ_w = η_w, and modifies Eq. (15) as follows:

$$\begin{array}{c}\hfill Z_g=\left(\frac{q}{{\eta }_a}\right)[1-{\left(\frac{M_g}{M_i}\right)}^{\frac{{\eta }_a}{\alpha -{\eta }_a+{\eta }_w}}].\end{array}$$(C.1)

This solution, already reported by Erb (2008) in her Eq. 11, places the dependence of η_w only in the exponent of M_g/M_i, while η_a also appears in the denominator in the Z_g expression.

The very interesting outcome of this exercise is that the η_w derived in this way reproduce the values of ζ_w shown in Fig. 3 and obtained using the more general formalism (using Z_w ≠ Z_g) of Eq. (15). Thus, metal loading ζ_w and mass loading η_w in the winds can be confused, depending on the assumptions made in the application of the formalism.

All Tables

All Figures

Fig. 1. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA binned by Mstar. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF median) is shown by a vertical dashed line.

In the text

Fig. 2. Difference of mass-loading factors Δη (top panel) and metal-loading factors ζw (bottom panel) as a function of virial velocities obtained from the BTFR from Sales et al. (2017). In both panels, filled symbols show the Bayesian results, and in the top panel, open ones illustrate the upper bound of the PDF 1σ confidence intervals (CI). Different symbol types correspond to the interval of Δη in the Bayesian priors, as described in the text; the symbols for Δη ∈ [0, 3.4] and Δη ∈ [ − 10, 3.4] are (arbitrarily) offset along the abscissa for better visibility. The gray regions give the ±1σ uncertainties in the best-fit values. Also shown is the “equilibrium” asymptote, namely Δη = α (see Sect. 3.1). In the top panel, the two estimates of Δη from different priors are exactly coincident. In the bottom panel, the formulations by Peeples & Shankar (2011) for two O/H calibrations, redone here for our value of y, are given by (purple) short-dashed and (blue) long-dashed curves for the PP04N2 and D02 calibrations, respectively. Also shown are the metal-loading factors computed by Chisholm et al. (2018) from COS data. The text gives more details.

In the text

Fig. 3. Difference of mass-loading factors Δη and metal-loading factors ζw inferred from the Bayesian analysis, and observed Zg in the bottom panel, as a function of Mstar on the left, and as a function of μg on the right. In the two upper panels, the best-fit Bayesian values are shown as filled symbols (circles or squares according to the Δη prior), with error bars corresponding to the ±1σ excursion, which is also illustrated by light-gray shaded regions. The rightmost point (magenta circle) in the top and middle right panels corresponds to the “extra” high-μg bin shown in Fig. 6. In the top panels, a horizontal dashed line illustrates the “equilibrium solution” with Δη = α. In the top left panel, the (purple) dashed curve is a linear power-law fit of Δη to log(Mstar), while the dashed curve in the top right panel shows Eq. (27), with d log Mg/d log Mstar = 0.72, as found in Paper II. The fits to the curves for ζw in the middle panels are cubic: log(ζw) as a function of log(Mstar) in the left panel, and ζw as a function of Fg in the right. In the lower panels, observed MAGMA Zg are shown as small filled circles, color coded by gas fraction, while the Zg medians and 1σ deviations are shown by large (black) filled circles and error bars. The dashed curves in the lower panels correspond to the calculation of Zg based on the fits in the upper panels (as a function of μg on the left, and log(Mstar) on the right): the long-dashed curves (dark orange) correspond to the Δη ≥ 0 priors, while the short-dashed ones (purple) to the “symmetric” priors. The green curves in the left panels are described in the text, obtained by fitting individual galaxy Δη and ζw calculated using Eqs. (27) and (28).

In the text

	Fig. 4. Baryonic gas fraction μg plotted against (log) Mstar for MAGMA galaxies. Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the gray region. As noted in the text, the two low-mass galaxies with low μg are Sextans A and WLM, with probable substantial underestimates of the total gas content.
In the text

Fig. 5. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA binned by μg. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line.

In the text

Fig. 6. Corner plot of χ2 surface as a function of the model parameters (Δη, ζw) for the highest μg bin MAGMA, here with a higher boundary than in Fig. 5: μg > 0.75. As in previous figures, the violet contours correspond to the minimum χ2 value, and the top and right panels report the probability density distributions for the marginalized parameters. Confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line.

In the text

Fig. 7. Effective yield as a function of μg. Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the surrounding gray regions. The two dashed curves correspond to the predictions from the Δη and ζw parameterizations (see also Fig. 3), with the symmetric Δη priors given as a short-dashed (purple) curve, and the Δη ≥ 0 priors as a long-dashed (orange) one. The horizontal dashed line corresponds to our adopted true yield y = 0.037 (see text).

In the text

Fig. 8. Effective yield as a function of Mstar (left panel) and Vcirc (right panel). Galaxies are coded by gas fraction, Fg, as shown in the vertical color bar. MAGMA medians are shown as a heavy gray line, with ±1σ excursions as the surrounding gray regions. Also shown in the left panel as two dashed curves are the predictions from the Δη and ζw parameterizations (see also Fig. 3), with the symmetric Δη priors given as a short-dashed (purple) curve, and the Δη ≥ 0 priors as a long-dashed (orange) one. The right panel illustrates galaxies with high star-formation efficiency (short gas depletion times) with a blue dashed curve, and those with low efficiency (long depletion times) as a red one. The horizontal short-dashed line corresponds to the yield of 0.037 adopted in this paper. See text for more details.

In the text

Fig. 9. Ratio of inferred wind and ISM metallicities Zw/Zg plotted against virial velocity Vcirc. Also shown (with arbitrary scaling) is the shape of Bayesian ζw as also plotted in Fig. 2; it is clearly not a power law, thus the different power-law forms of ηw introduce inflections in Zw/Zg. Left panel: as indicated in the legend, different trends of Zw/Zg are inferred by dividing the fitted form of ζw (see Eq. (22)) by various scalings of ηw (see text for more details). The shaded regions correspond to the trends allowed by the two different Bayesian estimates of ζw according to the different Δη priors. Right panel: broken-power law trend for ηw found by Muratov et al. (2015) is shown as a solid heavy (gray) curve, where the inflection at Vcirc = 60 km s−1 corresponds to a sharp decrease of Zw/Zg (see text for more details). The short-dashed (dark-blue) curve shows the relation found by Mitchell et al. (2020), in the low-mass regime where stellar feedback dominates, and the long-dashed (blue) curve shows the observed relation derived by Chisholm et al. (2018).

In the text

Fig. A.1. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF mode for Δη, median for ζw) is shown by a vertical dashed line. Here, Δη ∈ [ − 10, 3.4] (see Fig. 1 for the “symmetric” Δη intervals).

In the text

Fig. A.2. Corner plots of χ2 surface as a function of the model parameters (Δη, ζw) for MAGMA. The violet contours correspond to the minimum χ2 value. The top and right panels of each corner plot report the probability density distributions for the marginalized parameters; confidence intervals (±1σ) are shown as violet-tinted shaded rectangular regions, and the MLE (PDF median) is shown by a vertical dashed line. Here, Δη ∈ [0,3.4] (see Fig. 1 for the “symmetric” Δη intervals).

In the text

	Fig. B.1. Inferred galaxy age tgal plotted against Mstar for MAGMA galaxies (top panel); initial SFR at t = 0, ψ0 plotted against Mstar (middle); τg (assumed equal to τSF) plotted against Mstar (bottom). In the two top panels, we assumed a simple declining exponential SFH with β = −1 [Eq. (B.1)], and that ${\tau }_{\mathrm{SF}}={\tau }_g={{ϵ}_{\star }}^{-1}$.
In the text

Fig. B.2. Time evolution of Mg vs Mstar over the last 50% of the galaxies’ lifetime, according to the toy-model based on three exponential SFHs: β = ( − 1, 0, 1), and also assuming that the SFR timescale τSF is the same as the gas depletion time, τg = ${{ϵ}_{\star }}^{-1}$. Galaxy age tgal is given by Eq. (B.4), ψ0 by Eq. (B.3), and the Mg evolution by Eq. (14). The median curves in different Mstar bins are plotted, where the (median observed) gas mass Mg is taken to be Mg(tgal). For Δη = α, the equilibrium solution, Mg = Mi throughout the lifetime of the galaxy.

In the text