Gas-phase Metallicity as a Diagnostic of the Drivers of Star Formation on Different Spatial Scales

The basic idea of the gas-regulator model is that the formation of stars is instantaneously determined by the mass of a cold gas reservoir, which is regulated by the interplay between inflow, outflow, and star formation (Lilly et al. 2013). The instantaneous SFR can be written as

$$\begin{eqnarray}&&\mathrm{SFR}(t)={M}_{\mathrm{gas}}(t)\cdot \mathrm{SFE}(t),\end{eqnarray} \tag{ 1 }$$

where SFE(t) is the instantaneous SFE. We note that the SFE is the inverse of the gas depletion time (τ_dep) by definition, i.e., SFE ≡ 1/τ_dep. Following the work of Lilly et al. (2013) and W19, we assume that the mass loss due to outflow is scaled by the instantaneous SFR(t) with a mass-loading factor λ, i.e., λSFR(t). The fraction of stellar mass that is returned to the interstellar medium through winds and supernova explosions is denoted as R. We utilize the instantaneous return assumption and take R = 0.4 from stellar population models (e.g., Bruzual & Charlot 2003). We define the effective gas depletion timescale (τ_def,eff) as the gas depletion timescale due to not only star formation but also the wind-loading outflow, i.e., τ_dep,eff = τ_def/(1 − R + λ). We denote the inflow rate as Φ(t) and the metallicity of infalling gas as Z₀. The metal mass of the gas reservoir is denoted as M_Z(t). The yield, i.e., the mass of metals returned to the interstellar medium per unit mass of instantaneously formed stars, is denoted as y.

The basic continuity equations for gas and metals are (see Equations (9) and (20) in Lilly et al. 2013)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dM}}_{\mathrm{gas}({\rm{t}})}}{{dt}} & = & \ {\rm{\Phi }}(t)-(1-R)\cdot \mathrm{SFR}(t)-\lambda \cdot \mathrm{SFR}(t)\\ & = & \ {\rm{\Phi }}(t)-(1-R+\lambda )\cdot \mathrm{SFE}(t)\cdot {M}_{\mathrm{gas}}(t)\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\frac{{{dM}}_{{\rm{Z}}}(t)}{{dt}} & = & y\mathrm{SFR}(t)-Z(t)\cdot (1-R+\lambda )\cdot \mathrm{SFR}(t)+{\rm{\Phi }}(t){Z}_{0}\\ & = & \ y\mathrm{SFE}(t)\cdot {M}_{\mathrm{gas}}(t)-(1-R+\lambda )\mathrm{SFE}(t)\cdot {M}_{{\rm{Z}}}(t)\\ & & +\ {\rm{\Phi }}(t)\cdot {Z}_{0},\end{array}\end{eqnarray} \tag{ 3 }$$

where Z(t) = M_Z(t)/M_gas(t) by definition.

We apply the instantaneous recycling approximation in this work, i.e., we assume that the metal enhancement by star formation is instantaneous. We are ignoring the timescale for supernova ejecta to mix with the interstellar medium and start forming a new generation of stars (see more discussion in Section 6.3), which can be a few tens of megayears (Roy & Kunth 1995).

In Equation (2) and (3), there are a total of five quantities driving the solution: the possibly time-varying Φ(t), SFE(t), and the (assumed constant) λ, Z₀, and y. The M_gas(t), SFR(t), M_Z(t), and Z(t) are then the response of the regulator system to these five quantities. In this work we assume that λ, Z₀, and y are time independent. The mass-loading factor is tightly correlated with the stellar mass of galaxies (Hayward & Hopkins 2017) and is not likely to change significantly on timescales of gigayears given the current relatively low sSFR ∼ 0.1 Gyr⁻¹ of local SFMS galaxies. The yield y is a physical parameter reflecting the nuclear synthesis activity and the relative number of massive stars, which is expected to be tightly correlated with the IMF but is not likely to change significantly on gigayear timescales. However, we note that λ, Z₀, and even y may well change from galaxy to galaxy and even from region to region within individual galaxies.

Equations (2) and (3) are the two basic continuity equations in the gas-regulator model. The following analysis in Section 2 will focus on the solution of these two equations by driving the regulator system with sinusoidal Φ(t) or sinusoidal SFE(t) and investigate the correlation between the resulting instantaneous SFR(t) and the gas-phase metallicity Z(t). This correlation will be the main diagnostic used in our later analysis of observational data.

In Sections 2.2 and 2.3, we will investigate the properties of the gas-regulator model using a sinusoidally varying inflow rate or SFE. Although this is undoubtedly artificial in isolation, we can argue that any given inflow rate or SFE can be expressed as a linear combination of sinusoidal functions with different frequencies via the Fourier transform. The gas-regulator system will respond to these individual sinusoidal components independently. For any given time-varying inflow rate, the resulting M_gas (or M_Z) should be the superposition of the resulting M_gas (or M_Z) obtained from these sinusoidal variations of inflow rate. This statement is also true for SFE, when the variation of SFE is small (see details in Section 2.3).

We now first drive the gas-regulator system with time-varying inflow and time-invariant SFE. As in W19, we input the simple case that inflow rate is a linear sinusoidal function of time with a period of T_p:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)={{\rm{\Phi }}}_{0}+{{\rm{\Phi }}}_{{\rm{t}}}\cdot \sin (2\pi t/{T}_{{\rm{p}}}).\end{eqnarray} \tag{ 4 }$$

Then, we look for solutions in which M_gas(t) and M_Z(t) are sinusoidal functions phase-shifted from the inflow rate:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{gas}}(t) & = & \ {M}_{0}+{M}_{{\rm{t}}}\cdot \sin (2\pi t/{T}_{{\rm{p}}}-\delta )\\ {M}_{{\rm{Z}}}(t) & = & \ {M}_{{\rm{Z}}0}+{M}_{\mathrm{Zt}}\cdot \sin (2\pi t/{T}_{{\rm{p}}}-\beta ).\end{array}\end{eqnarray} \tag{ 5 }$$

In W19, by substituting M_gas(t) into Equation (2) and equalizing the various time-dependent terms in the usual way, we have the solution of M_gas(t):

$$\begin{eqnarray}\begin{array}{rcl}{M}_{0} & = & \ {{\rm{\Phi }}}_{0}{\tau }_{\mathrm{dep},\mathrm{eff}}\\ \delta & = & \ \arctan (\xi )\\ \displaystyle \frac{{M}_{{\rm{t}}}}{{M}_{0}} & = & \ \displaystyle \frac{1}{{\left(1+{\xi }^{2}\right)}^{1/2}}\times \displaystyle \frac{{{\rm{\Phi }}}_{t}}{{{\rm{\Phi }}}_{0}},\end{array}\end{eqnarray} \tag{ 6 }$$

where τ_def,eff is the effective gas depletion time, defined as τ_dep · (1 − R + λ)⁻¹ or SFE⁻¹ · (1 − R + λ)⁻¹ in Section 2.1, and ξ is the ratio of the effective gas depletion timescale to T_p(2π)⁻¹, i.e.,

$$\begin{eqnarray}&&\xi \equiv 2\pi {\tau }_{\mathrm{dep},\mathrm{eff}}/{T}_{{\rm{p}}}.\end{eqnarray} \tag{ 7 }$$

As we discussed in W19, the amplitude and phase delay of the output M_gas(t) strongly depends on the parameter ξ, the relative timescale of gas depletion time to the driving period. At fixed T_p, galaxies or regions with shorter gas depletion time are more able to follow the changes of the inflow rate, leading to a smaller phase delay and larger amplitude, and vice versa (see more discussion in Section 4 of W19).

In the similar way, we substitute Equations (5) and (6) into Equation (3) and equate the various time-dependent terms to find the solution of M_Z(t):

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{Z}}0} & = & \ ({y}_{\mathrm{eff}}+{Z}_{0}){{\rm{\Phi }}}_{0}{\tau }_{\mathrm{dep},\mathrm{eff}}\\ \beta \ & = & \ \arctan \left[\displaystyle \frac{2{y}_{\mathrm{eff}}\xi +{Z}_{0}\xi (1+{\xi }^{2})}{{y}_{\mathrm{eff}}(1-{\xi }^{2})+{Z}_{0}(1+{\xi }^{2})}\right]\\ \displaystyle \frac{{M}_{\mathrm{Zt}}}{{M}_{{\rm{Z}}0}} & = & \ \displaystyle \frac{{\left(1+{\eta }^{2}\right)}^{1/2}}{1+{\xi }^{2}}\times \displaystyle \frac{{{\rm{\Phi }}}_{{\rm{t}}}}{{{\rm{\Phi }}}_{0}},\end{array}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{y}_{\mathrm{eff}}\equiv y\cdot {\left(1-R+\lambda \right)}^{-1}\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\eta =\xi {Z}_{0}\cdot {\left({y}_{\mathrm{eff}}+{Z}_{0}\right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

If β is less than zero, then β should be equal to β+π. The shorthand η is defined for convenience. We remind readers that the effective yield y_eff defined in this way is different from some previous papers (e.g., Edmunds 1990; Garnett 2002). We prefer this definition because we believe that it is more fundamental.

If we assume that the inflow gas is in a pristine state, i.e., Z₀ ∼ 0, then the solution of M_Z(t) can be simplified further to be

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{Z}}0} & = & \ {y}_{\mathrm{eff}}{{\rm{\Phi }}}_{0}{\tau }_{\mathrm{dep},\mathrm{eff}}\\ \beta \ & = & \ \arctan \left[\displaystyle \frac{2\xi }{1-{\xi }^{2}}\right]\\ \displaystyle \frac{{M}_{\mathrm{Zt}}}{{M}_{{\rm{Z}}0}} & = & \ \displaystyle \frac{1}{1+{\xi }^{2}}\times \displaystyle \frac{{{\rm{\Phi }}}_{{\rm{t}}}}{{{\rm{\Phi }}}_{0}}.\end{array}\end{eqnarray} \tag{ 11 }$$

Interestingly, in this specific case with Z₀ ∼ 0, the phase delay of M_Z(t) is twice that of M_gas(t), i.e., β = 2δ. Similar to M_gas(t), the phase delay β and relative amplitude M_Zt/M_Z0 of M_Z(t) strongly depend on the parameter ξ. At fixed T_p, galaxies or regions with shorter (effective) gas depletion time can more easily follow the change of inflow rate and gas mass, resulting in smaller β and larger M_Zt/M_Z0. Specifically, if ξ is much less than unity, then both δ and β are close to zero, and both M_t/M₀ and M_Zt/M_Z0 are close to Φ_t/Φ₀. In other words, when the (effective) gas depletion time is much less than the driving period, i.e., ξ ≪ 1, the change of mass of the gas reservoir and of the mass of metals in the gas-regulator system can nearly follow the change of inflow rate, with little phase delay and with nearly the same relative amplitude of variation. If, however, ξ is much larger than 1, then δ is close to π/2, β is close to π, and both M_t/M₀ and M_Zt/M_Z0 are close to zero. This means that when the (effective) gas depletion time is much longer than the driving period, i.e., ξ ≫ 1, the gas-regulator system is unable to follow the relatively fast changes in the inflow rate, resulting in little variation in either M_gas(t) or M_Z(t).

The dependence of M_gas(t) and M_Z(t) on ξ can be clearly seen in the top left panel of Figure 1, where we show examples of the evolution of M_gas (blue), SFR (red), M_Z (orange), and Z (purple) when driving the gas-regulator system with periodic sinusoidal inflow. For illustrative purposes, we set Z₀ = 0, Φ_t/Φ₀ = 0.1, T_p = 1 Gyr, and $\mathrm{log}\xi$ = −0.5, 0.0, and 0.5.

Zoom In Zoom Out Reset image size

Given the solutions of M_gas(t) and M_Z(t) in Equations (6) and (8), the resulting SFR(t) and Z(t) can be easily obtained. Since the SFE is assumed here (in this subsection) to be time invariant, the change of SFR will exactly follow the change of cold gas mass. Therefore, the blue and red lines in the top left panel of Figure 1 are overlapped together. However, the Z(t), i.e., the ratio of M_Z(t) to M_gas(t), has a more complicated behavior than SFR(t), because it is not a sinusoidal function. The variation of the metallicity depends on the amplitude of variations in M_Z(t) and M_gas(t), as well as the phase delay between the two.

To clarify the correlation between the instantaneous SFR(t) and Z(t), we plot the $\mathrm{log}(\mathrm{SFR}(t)/\langle \mathrm{SFR}\rangle )$ versus $\mathrm{log}(Z(t)/\langle Z\rangle )$ and $\mathrm{log}({M}_{\mathrm{gas}}(t)/\langle {M}_{\mathrm{gas}}\rangle )$ versus $\mathrm{log}(Z(t)/\langle Z\rangle )$ for a set of different ξ in the top middle and right panels of Figure 1, where 〈Z〉, 〈SFR〉, and 〈M_gas〉 are the average metallicity, SFR, and cold gas mass, respectively. Since the $\mathrm{log}(\mathrm{SFR}(t)/\langle \mathrm{SFR}\rangle )$ is a relative quantity, i.e., $\mathrm{logSFR}(t)-\mathrm{log}\langle \mathrm{SFR}\rangle$, we also denote $\mathrm{log}(\mathrm{SFR}(t)/\langle \mathrm{SFR}\rangle )$ as ${\rm{\Delta }}\mathrm{log}$ SFR. In the same way, we denote $\mathrm{log}(Z(t)/\langle Z\rangle )$ as ${\rm{\Delta }}\mathrm{log}Z$ and $\mathrm{log}({M}_{\mathrm{gas}}(t)/\langle {M}_{\mathrm{gas}}\rangle )$ as ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$.

As shown, at all the different ξ shown here, the gas-regulator model predicts that ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ are negatively correlated when the system is driven with a sinusoidal inflow rate. The slope and the tightness of the ${\rm{\Delta }}\mathrm{log}$ SFR−${\rm{\Delta }}\mathrm{log}Z$ correlation strongly depend on ξ. Generally speaking, at fixed T_p, the correlation of ${\rm{\Delta }}\mathrm{log}$ SFR–${\rm{\Delta }}\mathrm{log}Z$ becomes weaker and steeper with increasing effective gas depletion time. The slope of the ${\rm{\Delta }}\mathrm{log}$ SFR−${\rm{\Delta }}\mathrm{log}Z$ relation is always steeper than −1. This means that the gas-regulator model requires that the scatter of ${\rm{\Delta }}\mathrm{log}Z$ is always less than or equal to the scatter of ${\rm{\Delta }}\mathrm{log}$ SFR. We will come back to this point later.

In addition to the sinusoidal inflow rate, we also for completeness explored the effect of a periodic step function in the inflow rate. We solve Equations (2) and (3) numerically for this case. The bottom panels of Figure 1 show the resulting M_gas(t), SFR(t), M_Z(t), and Z(t), as well as the resulting correlation between ${\rm{\Delta }}\mathrm{log}$ SFR (and ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$) and ${\rm{\Delta }}\mathrm{log}Z$. In generating the plots, we set the period of step function as 2 Gyr and change the upper-state duration of inflow rate (τ_s). We allow the τ_s to vary from 0.1τ_dep,eff to τ_dep,eff.

As shown in the bottom left panel of Figure 1, a sudden increase of inflow rate causes an increase of the SFR (or cold gas mass) and a decrease of gas-phase metallicity, and vice versa. This therefore also leads to a negative correlation between SFR (or cold gas mass) and metallicity, i.e., between ${\rm{\Delta }}\mathrm{log}$ SFR (or ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$) and ${\rm{\Delta }}\mathrm{log}Z$, consistent with the result for the sinusoidal variation in the top panels of Figure 1.

Of course, observationally we cannot follow the temporal evolution of a single galaxy. In W19 we therefore explored the scatter of the instantaneous SFR in a population of gas regulators, σ(SFR), when they are driven with simple sinusoidal Φ(t), and showed that within this population σ(SFR)/σ(Φ) is a monotonically decreasing function of ξ:

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\mathrm{SFR})}{\sigma ({\rm{\Phi }})}=\displaystyle \frac{1}{{\left(1+{\xi }^{2}\right)}^{1/2}}\cdot {\left(1-R+\lambda \right)}^{-1}.\end{eqnarray} \tag{ 12 }$$

In Equation (12), the scatter of SFR and Φ are calculated in linear space, while in the observations, the scatter of SFR is usually measured in logarithmic space. Here we present an approximate analytical solution of $\sigma (\mathrm{log}\,\mathrm{SFR})/\sigma (\mathrm{log}\,{\rm{\Phi }})$, which can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\mathrm{log}\,\mathrm{SFR})}{\sigma (\mathrm{log}\,{\rm{\Phi }})}\approx \displaystyle \frac{1}{{\left(1+{\xi }^{2}\right)}^{1/2}}.\end{eqnarray} \tag{ 13 }$$

As can be seen in Equation (13), if the scatter is measured in logarithmic space, the factor 1 − R + λ vanishes, since this can be viewed as a constant "inefficiency" in the star formation (see Lilly et al. 2013). The details to derive Equation (13) are given in Appendix A.

The left panel of Figure 2 shows the numerical solution (black solid curve) and the approximate analytical solution (gray dashed curve) of $\sigma (\mathrm{logSFR})/\sigma (\mathrm{log}{\rm{\Phi }})$ as a function of $\mathrm{log}\xi$, which are in excellent agreement. Specifically, in obtaining the numerical solution, we first solve Equation (2) to obtain the SFR(t) at a set of different ξ. Then, we calculate the σ(logSFR) within a single period for each ξ. The $\sigma (\mathrm{log}Z)$ is calculated in a similar way. The analytic solution provides the physical insight to how ξ determines the response of the gas-regulator model, while the numerical solution provides a double-check for the analytic solution and provides the judgment for the validation of the approximate analytic solution.

Zoom In Zoom Out Reset image size

As pointed out in W19 and Wang & Lilly (2020a), when driving the gas regulator with time-varying inflow rate, the amplitude of the variations in SFR is reduced from the amplitude of variations in the inflow rate by a frequency-dependent (ν = 1/T_p) response curve, i.e., Equation (13) or Equation (9) in W19. In other words, for a given inflow rate with any given ${\mathrm{PSD}}_{\mathrm{log}{\rm{\Phi }}}$(ν), the PSD of the resulting logSFR, PSD_logSFR(ν), can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{PSD}}_{\mathrm{logSFR}}(\nu ) & \approx & \displaystyle \frac{{\sigma }^{2}(\mathrm{logSFR})}{{\sigma }^{2}(\mathrm{log}{\rm{\Phi }})}\cdot {\mathrm{PSD}}_{\mathrm{log}{\rm{\Phi }}}(\nu )\\ & \approx & \ \displaystyle \frac{1}{1+{\xi }^{2}}\cdot {\mathrm{PSD}}_{\mathrm{log}{\rm{\Phi }}}(\nu )\\ & = & \ \displaystyle \frac{1}{1+{\left(2\pi {\tau }_{\mathrm{dep},\mathrm{eff}}\nu \right)}^{2}}\cdot {\mathrm{PSD}}_{\mathrm{log}{\rm{\Phi }}}(\nu ).\end{array}\end{eqnarray} \tag{ 14 }$$

Equation (14) is established under the condition that the variation of inflow rate is small, and therefore logSFR(t) and logΦ(t) are close to sinusoidal functions. However, as examined in Appendix B, we directly input the sinusoidal Φ(t) in logarithmic space with a large amplitude (0.5 dex), and find that Equation (14) is still valid. Therefore, we conclude that Equation (14) theoretically predicts the connection between the PSDs of $\mathrm{log}$ SFR and $\mathrm{log}{\rm{\Phi }}$, when driving the gas regulator with time-varying inflow rate (see also W19; Wang & Lilly 2020b; Tacchella et al. 2020). However, from the observations, the inflow rate history of galaxies is not, of course, a directly observable quantity.

Finally, we come to the scatter of the gas-phase metallicities in the population. As shown in top middle panel of Figure 1, the ratio of the scatter in Z(t) to the scatter of SFR(t) (i.e., what would be observed in a given population of regulators at a fixed time) is predicted to be strongly dependent on ξ. Here we present the approximate analytical solution of $\sigma (\mathrm{log}Z)/\sigma (\mathrm{logSFR})$, which can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\mathrm{log}Z)}{\sigma (\mathrm{logSFR})}\approx \displaystyle \frac{\xi }{{\left(1+{\xi }^{2}\right)}^{1/2}}\cdot \displaystyle \frac{1}{1+{Z}_{0}/{y}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 15 }$$

The detailed derivation of Equation (15) is shown in Appendix A. Similar to Equation (13), Equation (15) provides the link between the variability of $\mathrm{log}$ SFR and $\mathrm{log}$ Z. However, one cannot write the correlation of PSD_logSFR and PSD_logZ in a similar way to Equation (14), because the resulting metallicity is not a sinusoidal function. Since the σ(logZ)/σ(logSFR) only depends on ξ at fixed Z₀/y_eff, we therefore present both the numerical (blue, green, and red curves) and the analytic solutions (gray dashed curves) of $\sigma (\mathrm{log}Z)/\sigma (\mathrm{logSFR})$ in the left panel of Figure 2 at three different Z₀/y_eff, which are again in excellent agreement. For different Z₀/y_eff, σ(logZ)/σ(logSFR) monotonically increases with $\mathrm{log}\xi$, which is opposite to $\sigma (\mathrm{logSFR})/\sigma (\mathrm{log}{\rm{\Phi }})$. Intriguingly, if Z₀ = 0, Equations (13) and (15) are strictly symmetrical across the axis of $\mathrm{log}\xi =0$. Unlike Equation (13), Equation (15) relates two readily observable quantities, the instantaneous SFR and the instantaneous gas-phase metallicity.

In the previous subsection, we looked at the behavior of the gas-regulator system with time-varying inflow and time-invariant SFE. In this subsection, we will explore the behavior when the regulator is driven with a constant inflow rate but experiences a time-varying SFE.

Similar to Section 2.2, we first input a time-invariant inflow, i.e., Φ(t) = Φ₀, and a sinusoidally varying SFE(t):

$$\begin{eqnarray}&&\mathrm{SFE}(t)={\mathrm{SFE}}_{0}+{\mathrm{SFE}}_{t}\cdot \sin (2\pi t/{T}_{{\rm{p}}}).\end{eqnarray} \tag{ 16 }$$

The driving period of SFE(t) is again denoted as T_p. As before, we look for the solution of M_gas(t) and M_Z(t) in terms of Equation (5). However, we note that, unlike in Section 2.2, the solutions of M_gas(t) and M_Z(t) are not exact sinusoidal functions but can be represented as approximations to sinusoids. We assume that the variation of the input SFE(t) is small, i.e., SFE_t ≪ SFE₀. Therefore, the variations of the resulting M_gas(t) and M_Z(t) are also small, i.e., M_t ≪ M₀ and M_Zt ≪ M_Z0.

Actually, based on Equations (2) and (3), varying the SFE is mathematically (but, of course, not physically) equivalent to varying inflow rate as SFE_t = −Φ_t/Φ₀ × SFE₀ and SFE_t ≪ SFE₀. Therefore, the solution of Equations (2) and (3) can be directly written as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{0} & = & \ {{\rm{\Phi }}}_{0}{\tau }_{\mathrm{dep},\mathrm{eff}}\\ \delta \ & = & \ \arctan (\xi )\\ \frac{{M}_{{\rm{t}}}}{{M}_{0}} & = & \ -\frac{1}{{\left(1+{\xi }^{2}\right)}^{1/2}}\times \frac{{\mathrm{SFE}}_{{\rm{t}}}}{{\mathrm{SFE}}_{0}}\end{array}\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{Z}}0} & = & \ ({y}_{\mathrm{eff}}+{Z}_{0}){{\rm{\Phi }}}_{0}{\tau }_{\mathrm{dep},\mathrm{eff}}\\ \beta \ & = & \ \arctan \left[\displaystyle \frac{2{y}_{\mathrm{eff}}\xi +{Z}_{0}\xi (1+{\xi }^{2})}{{y}_{\mathrm{eff}}(1-{\xi }^{2})+{Z}_{0}(1+{\xi }^{2})}\right]\\ \displaystyle \frac{{M}_{\mathrm{Zt}}}{{M}_{{\rm{Z}}0}} & = & \ -\displaystyle \frac{{\left(1+{\eta }^{2}\right)}^{1/2}}{1+{\xi }^{2}}\times \displaystyle \frac{{\mathrm{SFE}}_{{\rm{t}}}}{{\mathrm{SFE}}_{0}}.\end{array}\end{eqnarray} \tag{ 18 }$$

We emphasize that the definitions of ξ and η are the same as in Section 2.2.

Comparing Equations (17) and (18) with Equations (6) and (8), the difference in sign indicates the fact that, when driving the regulator with time-varying Φ(t), the increase of Φ(t) produces an increase of both M_gas(t) and M_Z(t) with some time delays, while when driving the regulator with time-varying SFE(t), the increase of SFE(t) leads to decreases of both M_gas(t) and M_Z(t), again with some time delay. We note that the difference in sign cannot be simply treated as an additional phase delay of π, i.e., a time lag of T_p/2, because an inverse connection between the mass of cold gas and the metallicity is physically expected.

The solutions in Equations (17) and (18) are illustrated in the top left panel of Figure 2, where we show examples of the evolution of M_gas (blue), SFR (red), M_Z (orange), and Z (purple). As previously, we also investigate the correlation between ${\rm{\Delta }}\mathrm{log}$ SFR (and ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$) and ${\rm{\Delta }}\mathrm{log}Z$ at a set of different ξ, as shown in the top middle (and top right) panel of Figure 3.

Zoom In Zoom Out Reset image size

In contrast with the result in Section 2.2, ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ now show a strong positive correlation. However, we note that ${\rm{\Delta }}\mathrm{log}\,{M}_{\mathrm{gas}}$ and ${\rm{\Delta }}\mathrm{log}\,Z$ remain negatively correlated, as would be expected, independent of whether the gas regulator is driven with time-varying inflow or time-varying SFE. In this sense, metallicity variations fundamentally reflect variations in total gas mass in the gas-regulator reservoir (see Lilly et al. 2013; Bothwell et al. 2013, 2016).

As in Section 2.2, we also look at periodic step functions for the time variation of SFE. Such changes in SFE may well be more realistic than sinusoidal changes. The bottom panels of Figure 3 demonstrate the resulting M_gas(t), SFR(t), M_Z(t), and Z(t), as well as the correlation between ${\rm{\Delta }}\mathrm{log}$ SFR (and ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$) and ${\rm{\Delta }}\mathrm{log}Z$. In generating the plots, we set the period of the step function to be 2 Gyr and change the upper-state duration (τ_s). We allow the τ_s to vary from 0.1τ_dep,eff to τ_dep,eff. Here the τ_dep,eff is calculated based on the SFE in its lower state.

A sudden increase of SFE causes a sudden increase, but a following decrease, of the SFR, a following decrease of M_gas, and a following increase of the gas-phase metallicity. As a whole, it is not as immediately clear what the sign of the ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ correlation will be, given the fact that the lower branch has many more data points (reflecting the longer interval of time) than the upper branch in the bottom middle panel of Figure 3. Therefore, it appears that the sign of the correlation between SFR and metallicity depends on the detailed forms of input time-varying SFE. There is a strong asymmetry in the distribution of SFR through the cycle, indicated by the number density of the data points in the bottom middle panel of Figure 3. Specifically, SFR(t) stays close to its median value for most of the time but shows a strong increase for a short period. The asymmetry becomes more significant as the relative duration of the increased phase of SFE is decreased. However, one thing is clear: the states with strongly enhanced SFR are always metal enhanced with respect to the mean metallicity. These phases are represented in the upper locus of points in the figure that have Z > 〈Z〉.

Consistent with the top right panels of Figure 3, for ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{gas}}$ and ${\rm{\Delta }}\mathrm{log}Z$ we conclude that there will always overall be a negative correlation that will be most clearly seen in the highest SFR points.

Similar to Section 2.2, we again present the approximate analytical solution of σ(logSFR)/σ(logSFE), and $\sigma (\mathrm{log}Z)/\sigma (\mathrm{logSFR})$ when driving the gas regulator with sinusoidal SFE. These quantities can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\mathrm{logSFR})}{\sigma (\mathrm{logSFE})}\approx \displaystyle \frac{\xi }{{\left(1+{\xi }^{2}\right)}^{1/2}}\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\mathrm{log}Z)}{\sigma (\mathrm{logSFR})}\approx \displaystyle \frac{1}{{\left(1+{\xi }^{2}\right)}^{1/2}}\cdot \displaystyle \frac{1}{1+{Z}_{0}/{y}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 20 }$$

In a similar way to that in Section 2.2, according to Equation (19), we can write the correlation between the PSDs of $\mathrm{log}$ SFR(t) and $\mathrm{log}$ SFE(t) as

$$\begin{eqnarray}&&{\mathrm{PSD}}_{\mathrm{logSFR}}(\nu )\approx \frac{1}{1+{\left(2\pi {\tau }_{\mathrm{dep},\mathrm{eff}}\nu \right)}^{-2}}\cdot {\mathrm{PSD}}_{\mathrm{logSFE}}(\nu ).\end{eqnarray} \tag{ 21 }$$

Equation (21) is established when the variation of SFE is small, i.e., SFE_t ≪ SFE₀. The right panel of Figure 2 shows the numerical solution (solid curves) and the approximate analytical solution (gray dashed curves) of σ(logSFR)/σ(logSFE) and $\sigma (\mathrm{log}Z)/\sigma (\mathrm{logSFR})$ as a function of $\mathrm{log}\xi$. As shown, the numerical solution is well matched by the analytical solution. Intriguingly, Equations (19) and (20) are strictly symmetrical to Equations (13) and (15), respectively, at the axis of $\mathrm{log}\xi =0$.

When driving the gas regulator with a time-varying SFE, the σ(logSFR)/σ(logSFE) is predicted to increase with ξ, while $\sigma (\mathrm{log}Z)/\sigma (\mathrm{logSFR})$ is predicted to decrease with ξ. Although Equations (13), (15), (19), and (20) are only approximate solutions obtained in the limit of small variations of inflow rate or SFE, we have verified numerically that these equations remain reasonable approximations even when the variations of inflow rate or SFE are quite significant. In Appendix B, we examine the correlation of ${\rm{\Delta }}\mathrm{log}$ SFR versus ${\rm{\Delta }}\mathrm{log}$ Z, when there is a large variation (0.5 dex) of inflow rate or SFE. We confirm that the model predictions of ${\rm{\Delta }}\mathrm{log}$ SFR versus ${\rm{\Delta }}\mathrm{log}Z$ in the top middle panels of Figures 1 and 3 are scalable to the case in which the variation of driving factor is quite significant.

In Sections 2.2 and 2.3, we have explored the behavior of the gas-regulator system when it is driven by time-varying inflow and time-varying SFE, respectively. In this subsection, we will explore how the metallicities are modified by changes in the wind mass-loading factor and the metallicity of inflow gas. The assumed yield enters as a simple scaling factor throughout.

Following Section 2.2, we drive the gas regulator with a sinusoidally varying inflow rate for a set of different mass-loading factors and Z₀. First, we set Z₀ = 0 and show the ${\rm{\Delta }}\mathrm{log}$ SFR versus $\mathrm{log}\,Z-\mathrm{log}\,y$ in the top panel of Figure 4 for different mass-loading factors. To eliminate the dependence on the value of y, we set the x-axis to be $\mathrm{log}Z-\mathrm{log}y$ rather than $\mathrm{log}Z$. At Z₀ = 0, the relative changes in metallicity and star formation, i.e., the ${\rm{\Delta }}\mathrm{log}$ SFR−${\rm{\Delta }}\mathrm{log}Z$ relation, stay the same with varying λ, while the mean metallicity decreases with increasing λ. We then set Z₀ = 0.5y and obtain the ${\rm{\Delta }}\mathrm{log}$ SFR versus $\mathrm{log}Z-\mathrm{log}y$ shown in the middle panel of Figure 4. As shown, the basic negative correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ still remains, but the relative change in metallicity, i.e., the scatter of metallicity in the population, decreases with increasing Z₀ and increasing λ (if Z₀ ≠ 0), compared to the top panel of Figure 4. Finally, we set Z₀ = 0.5y_eff, and the ${\rm{\Delta }}\mathrm{log}$ SFR versus $\mathrm{log}Z-\mathrm{log}y$ is shown in the bottom panel of Figure 4. The correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ is the same when Z₀/y_eff is fixed.

Zoom In Zoom Out Reset image size

These results follow Equation (15), where it is clear that the ratio of $\sigma (\mathrm{log}Z)$ to σ(logSFR) only depends on Z₀/y_eff once ξ is fixed. At a given ξ, the scatter of $\mathrm{log}Z$ is scaled by the factor ${(1+{Z}_{0}/{y}_{\mathrm{eff}})}^{-1}$ with respect to $\sigma (\mathrm{log}Z)$ at Z₀ = 0. Another interesting point is that the mean metallicity does depend on Z₀, λ, and y. As a whole, the mean gas-phase metallicity increases with increasing Z₀ and/or decreasing λ. Actually, the mean SFR and metallicity can be solved analytically based on Equations (6) and (8) (or Equations (17) and (18)), which can be written as (see also Sánchez Almeida et al. 2014)

$$\begin{eqnarray} &&\langle \mathrm{SFR} \rangle ={M}_{0}\cdot \mathrm{SFE}={{\rm{\Phi }}}_{0}\cdot {\left(1-R+\lambda \right)}^{-1}\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&\langle Z\rangle =\displaystyle \frac{{M}_{{\rm{Z}}0}}{{M}_{0}}={Z}_{0}+{y}_{\mathrm{eff}}.\end{eqnarray} \tag{ 23 }$$

For completeness, we also look at driving the gas-regulator system with sinusoidal SFE. The results are shown in Figure 5. As shown, the effects of Z₀, λ, and y on the correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ follow Equation (20), and the effects of them on mean SFR and metallicity follow Equations (22) and (23). Although here we do not try to set a different yield, we argue that the effect of varying yield would also follow Equation (15) or Equation (20).

Zoom In Zoom Out Reset image size

Based on Equations (22) and (23), we find that the mean SFR is determined by the mean inflow rate and mass-loading factor regardless of the SFE, simply because the gas reservoir adjusts to maintain long-term balance, while the mean metallicity only depends on Z₀ and y_eff regardless of how large is the inflow or SFE. We conclude that, under the gas-regulator framework, the metallicity is primarily determined by Z₀ and y_eff with a secondary effect on SFR (or cold gas mass) due to the time variation of the inflow rate or SFE.

It is important to note that the metallicity does not depend on the absolute values of inflow rate and SFE, but rather on the change in them. Therefore, when investigating the question whether the SFR is a secondary parameter to determine the metallicity, one should look at the correlation of the relative values of SFR and Z (or residuals like ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ in Figures 1 and 3), rather than the absolute values, in order to eliminate different 〈SFR〉 and 〈Z〉 for different galaxies or regions.

In Sections 2.2 and 2.3, we have investigated the properties of the gas regulator by driving the gas-regulator system with (1) time-varying Φ and time-invariant SFE and (2) time-varying SFE and time-invariant Φ. The most important result is that there are opposite correlations (negative and positive, respectively) between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ for these two modes.

However, in the real universe, the inflow rate and SFE could both be time-varying. If the variation of inflow rate is much more dominant than the variation of SFE, a negative correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ is expected according to the analysis of Section 2.2, and vice versa. This implies that the correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ (at least for regions of significantly enhanced star formation) is a powerful and observationally accessible diagnostic of the dominant mechanism for the variation of SFR, changes in inflow, or changes in SFE, even if both of them may vary with time.

In the whole of Section 2, we have used ${\rm{\Delta }}\mathrm{log}$ SFR to indicate the temporal elevation and suppression of star formation. We note that one could also use the displacement of sSFR to its smooth cosmic evolution in logarithmic space, ${\rm{\Delta }}\mathrm{log}$ sSFR. These are essentially equivalent provided that the timescales of interest (essentially a single period of oscillation) are much shorter than the inverse sSFR, so that the mass changes negligibly through the oscillation cycle. We have examined that the difference is negligible in the case of sSFR ∼ 10⁻¹⁰ yr⁻¹ for typical local SF galaxies when replacing ${\rm{\Delta }}\mathrm{log}$ SFR with ${\rm{\Delta }}\mathrm{log}$ sSFR. This is true for the present work, since the galaxies used in the observational part are low-redshift galaxies. As discussed below, there are good reasons to use ${\rm{\Delta }}\mathrm{log}$ sSFR as the observational measure of the relative SFR rather than ${\rm{\Delta }}\mathrm{log}$ SFR.

In addition, the metallicity we measure from the observation is the oxygen abundance, defined as the number ratio of oxygen to hydrogen. We argue that in Equation (3) the mass of metals can be replaced by the number of oxygen, if the y and Z₀ are defined in the terms of the number fraction of oxygen. Therefore, the model predictions in the whole of Section 2 in terms of Z are also valid in terms of O/H, including Figures 1–5, as well as the corresponding equations.

The models predict the correlation of the temporal changes in SFR and Z within a given period for an individual gas-regulator system. Observationally, it is, of course, not possible to monitor a single galaxy or region for the timescales of interest. Instead, we have the SFR and metallicity for a large population of galaxies or regions. Therefore, we must make the assumption of ergodicity in order to compare the models with the observational data. In other words, we must assume that the variations that we see of the SFR and Z across the observed galaxy population (or between different regions within a galaxy) observed at a single epoch do indeed reflect the typical temporal changes of the SFR and Z within a given system (see more details in Section 3.2 of Wang & Lilly 2020a).

The final released sample ¹ of the MAD survey includes 38 weakly inclined, spiral galaxies, spanning a large range of stellar mass from 10^8.5 to 10^11.2 M_⊙. These galaxies were observed during the MUSE Guaranteed Time Observing runs on the Very Large Telescope (VLT) between 2015 April and 2017 September. The on-source exposure time was 1 hr for each galaxy, and the seeing ranged between 04 and 09. These galaxies are very nearby, with z < 0.013, leading to an average spatial resolution of ∼100 pc or better. The MUSE field of view is 1 arcmin², and the wavelength coverage of the spectra is from 4650 to 9300 Å. The coverage for MAD galaxies is in the range of ∼0.5–3 effective radius (R_e), with a median value of 1.3 R_e. The data were reduced using the MUSE pipeline (Weilbacher et al. 2012), including bias and dark subtraction, flat-fielding, wavelength calibration, and so on.

The data downloaded from the data release are not the original two-dimensional spectra, but the derived measurements from the spectra. These include the strengths of the emission lines, such as the flux map of Hβ, [O iii] λλ4959, 5007, Hα, [N ii] λ λ6548, 6583, and [S ii] λλ6717, 6731. The emission lines are modeled by single Gaussian profiles. The fluxes of emission lines are corrected for the dust attenuation by using the Balmer decrement with case B recombination. The intrinsic flux ratio of Hα to Hβ is taken to be 2.86, assuming the electron temperature of 10⁴ K and the electron density of 10² cm⁻³. The adopted attenuation curve is the CCM (Cardelli et al. 1989; O'Donnell 1994) extinction curve with R_V = 3.1. In addition to the maps of emission lines, the released MAD data also include the maps of derived quantities, such as SFR surface density (Σ_SFR) and stellar mass surface density (Σ_*). The SFR is derived from the Hα luminosity (Kennicutt 1998; Hao et al. 2011), assuming the Kroupa (2001) IMF. The stellar mass density is derived by fitting stellar population models to the stellar continuum spectra twice. The first continuum fitting is performed spaxel by spaxel, and the second fitting was performed on the stellar Voronoi tessellation (signal-to-noise ratio of 50 on the continuum) using the stellar templates from MILES (Sánchez-Blázquez et al. 2006) with pPXF (Cappellari & Emsellem 2004). Then, the resulting Σ_* map is a spaxel-by-spaxel-based map, assuming that the Σ_* is the same for all spaxels within a single Voronoi bin.

We note that in the released maps, spaxels that are located within the Seyfert and LINER regions of the BPT diagram (e.g., Baldwin et al. 1981; Kewley et al. 2001) are masked out. The SFR and gas-phase metallicity cannot be well measured in these masked regions. Their exclusion should not affect our analysis; therefore, in this work we only focus on the SF and "composite" regions. The SF and composite regions are identified by the demarcation lines of Kauffmann et al. (2003) and Kewley et al. (2001) on the BPT diagram (see Figure 1 in Erroz-Ferrer et al. 2019). This also means that for each galaxy we only use some fraction of the spaxels within the MUSE field. The fraction of valid spaxels varies from galaxy to galaxy, with a median value of 0.33. We have checked that our result does not depend on the fraction of spaxels that are used.

The highly resolved MAD data enable us to investigate the correlation between star formation and metal enhancement down to the scale of giant molecular clouds (GMCs). However, the MAD sample only includes 38 galaxies, which limits the statistical power when examining the galaxy population as a whole. Therefore, we utilize complementary data on the integrated spectra of galaxies from the MaNGA survey. We do not use the individual spaxel data for the MaNGA galaxies because the resolution is so much worse than for MAD.

MaNGA is the largest IFS survey of nearby galaxies up to now, and it aims at mapping the two-dimensional spectra for ∼10,000 galaxies with redshifts in the range of 0.01 < z < 0.15. Using the two dual-channel BOSS spectrographs at the Sloan Telescope (Gunn et al. 2006; Smee et al. 2013), MaNGA covers the wavelength of 3600–10300 Å at R ∼ 2000. The spatial coverage of individual galaxies is usually larger than 1.5 R_e with a resolution of 1–2 kpc. The flux calibration, including the flux loss due to atmospheric absorption and instrument response, is accurate to better than 5% for more than 89% of MaNGA's wavelength range (Yan et al. 2016).

In this work, we use the well-defined sample of SF galaxies from W19. Here we only briefly describe the sample definition, and we refer the reader to W19 for further details. This galaxy sample is originally selected from the SDSS Data Release 14 (Abolfathi et al. 2018), excluding quenched galaxies, mergers, irregulars, and heavily disturbed galaxies. The quenched galaxies are identified and excluded based on the stellar mass and SFR diagram. For each individual galaxy, the stellar mass and SFR are measured within the effective radius, i.e., M_*(<R_e) and SFR(<R_e), based on the MaNGA two-dimensional spectra. The final MaNGA sample includes 976 SF galaxies and is a good representation of normal SFMS galaxies.

Similar to the measurements of SFR for the MAD galaxies, the map of SFR surface density of MaNGA galaxies is also determined by the extinction-corrected Hα luminosity (Kennicutt 1998). The correction of dust attenuation follows the same approach for MAD galaxies, as described in Section 3.1, by using the Balmer decrement and adopting the CCM extinction curve. The maps of stellar mass surface density for MaNGA galaxies are obtained by fitting the stellar continuum using the public fitting code STARLIGHT (Cid Fernandes et al. 2004), using single stellar populations with Padova isochrones from Bruzual & Charlot (2003). However, we note that in determining the Σ_SFR and Σ_* for MaNGA galaxies, the Chabrier (2003) IMF is assumed, which is different from the one adopted for MAD galaxies. We argue that the two IMFs are quite close to each other, with only a small overall shift on SFR and M_* (or Σ_SFR and Σ_*), which does not change any of our conclusions in this work.

The T_e-based method is widely understood to represent the "gold standard" in determining the gas-phase metallicity (e.g., Skillman et al. 1989; Garnett 2002; Bresolin 2007; Berg et al. 2015; Bian et al. 2018). However, the measurement of the weak [O iii] λ4363 emission line is needed, which is often not detected in the available spectra. Therefore, a set of empirical recipes have been proposed to derive the gas-phase metallicity based only on the strong emission lines (e.g., Kobulnicky & Kewley 2004; Pettini & Pagel 2004; Maiolino et al. 2008; Pérez-Montero & Contini 2009; Pilyugin et al. 2010; Marino et al. 2013; Vogt et al. 2015), such as [O ii] λ3227, Hβ, [O iii] λ5007, Hα, [N ii] λ6584, and [S ii] λλ6717, 6731. However, systematic offsets of 0.2 dex or more are found between these different empirical calibrations, even using the same line measurements (Kewley & Ellison 2008; Blanc et al. 2015). Not only are there systematic offsets, but the ranges of the derived metallicities are also different. This is unfortunately important since we have argued that the variation of gas-phase metallicity (or the residuals of metallicity) as the SFR varies is an important diagnostic. The dispersion in metallicities must be considered in the context of the different ranges of the oxygen abundance measurements using the different methods.

Unfortunately, the wavelength coverage of MUSE does not include the [O ii] λ3727 and [O iii] λ4363 lines, leading to a limited number of usable strong-line prescriptions. In this work, we adopt the empirical relations from Dopita et al. (2016, hereafter DOP16) and Pilyugin & Grebel (2016, hereafter PG16). These two empirical relations are the most recently constructed and have advantages over the previous methods, although ultimately it is not known whether they are more accurate than the other methods. In the following, we briefly introduce these two calibrators.

3.3.1.  N2S2H α

The N2S2H α is a remarkably simple diagnostic proposed by DOP16, which can be written as

$$\begin{eqnarray}&&{\mathtt{N}}{\mathtt{2}}{\mathtt{S}}{\mathtt{2}}{\mathtt{H}}\alpha =\mathrm{log}([{\rm{N}}\,{\rm\small{ii}}]/[{\rm{S}}\,{\rm\small{ii}}])+0.264\mathrm{log}([{\rm{N}}\,{\rm\small{ii}}]/{\rm{H}}\alpha ),\end{eqnarray} \tag{ 24 }$$

where [N ii] is the flux of [N ii] λ6584 and [S ii] is the total flux of [S ii] λλ6717, 6731. Then, the empirical relation of the metallicity can be written as

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.77+{\mathtt{N}}{\mathtt{2}}{\mathtt{S}}{\mathtt{2}}{\mathtt{H}}\alpha +0.45{\left({\mathtt{N}}{\mathtt{2}}{\mathtt{S}}{\mathtt{2}}{\mathtt{H}}\alpha +0.3\right)}^{5}.\end{eqnarray} \tag{ 25 }$$

This simple empirical relation is calibrated in the range of $7.7\lt 12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 9.2$, which is valid for both MAD and MaNGA galaxies. Hα, [N ii] λ6584, and [S ii] λλ6717, 6731 are located close together in wavelength, limiting the spectral range needed, and making the N2S2H α diagnostic insensitive to reddening. The [N ii]/Hα term provides a correction for the weak dependence on the ionization parameter and gas pressure. DOP16 argued that this diagnostic should be a fairly reliable metallicity estimator with only small residual dependence on other physical parameters, and that it can be used in a wide range of environments.

However, this metallicity estimator strongly depends on the relative N/O ratio. In the calibration of N2S2H α, DOP16 assumed a universal correlation between N/O and O/H, which is determined based on a mixture of both stellar and nebular sources (Nicholls et al. 2017). This means that any galaxies or regions that deviate from the adopted N/O–O/H relation would have an extra uncertainty in metallicity measurement when using N2S2H α as the metallicity estimator. The metallicity also depends on the relative S/O ration, while DOP16 have adopted the [N ii/S ii] term in the metallicity indicator N2S2H α, which likely accounts for some of the dependence of metallicity on S/O.

3.3.2.  Scal

The S-calibration (Scal) metallicity estimator has been proposed by PG16, based on three standard diagnostic line ratios:

$$\begin{eqnarray}\begin{array}{rcl}{\mathtt{N}}{\mathtt{2}} & =\ & [{\rm{N}}\,{\rm\small{ii}}]\lambda \lambda 6548,6584/{\rm{H}}\beta ,\\ {\mathtt{S}}{\mathtt{2}} & =\ & [{\rm{S}}\,{\rm\small{ii}}]\lambda \lambda 6717,6731/{\rm{H}}\beta ,\\ {\mathtt{R}}{\mathtt{3}} & =\ & [{\rm{O}}\,{\rm\small{iii}}]\lambda \lambda 4959,5007/{\rm{H}}\beta .\end{array}\end{eqnarray} \tag{ 26 }$$

The Scal diagnostic is defined separately for the upper and lower branches of $\mathrm{log}{\mathtt{N}}{\mathtt{2}}$. The Scal indicator for the upper branch ($\mathrm{log}{\mathtt{N}}{\mathtt{2}}$ ≥ −0.6) can be written as

$$\begin{eqnarray}\begin{array}{rcl}12+\mathrm{log}({\rm{O}}/{\rm{H}}) & = & 8.424+0.030\mathrm{log}({\mathtt{R}}{\mathtt{3}}/{\mathtt{S}}{\mathtt{2}})+0.751\mathrm{log}{\mathtt{N}}{\mathtt{2}}\\ \ & & +\ (-0.349+0.182\mathrm{log}({\mathtt{R}}{\mathtt{3}}/{\mathtt{S}}{\mathtt{2}})\\ \ & & +\ 0.508\mathrm{log}{\mathtt{N}}{\mathtt{2}})\times \mathrm{log}{\mathtt{S}}{\mathtt{2}},\end{array}\end{eqnarray} \tag{ 27 }$$

and the Scal indicator for the lower branch ($\mathrm{log}{\mathtt{N}}{\mathtt{2}}$<−0.6) can be written as

$$\begin{eqnarray}\begin{array}{rcl}12+\mathrm{log}({\rm{O}}/{\rm{H}}) & = & 8.072+0.789\mathrm{log}({\mathtt{R}}{\mathtt{3}}/{\mathtt{S}}{\mathtt{2}})+0.726\mathrm{log}{\mathtt{N}}{\mathtt{2}}\\ \ & & +\ (1.069\mbox{--}0.170\mathrm{log}({\mathtt{R}}{\mathtt{3}}/{\mathtt{S}}{\mathtt{2}})\\ \ & & +\ 0.022\mathrm{log}{\mathtt{N}}{\mathtt{2}})\times \mathrm{log}{\mathtt{S}}{\mathtt{2}}.\end{array}\end{eqnarray} \tag{ 28 }$$

The Scal prescription is calibrated to metallicity measurements, from the T_e-based method, of several hundred nearby H ii regions. PG16 found that the Scal indicator gives metallicities in very good agreement with the T_e-based methods, with a scatter of only ∼0.05 dex across the metallicity range of $7.0\lt 12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 8.8$. Furthermore, the Scal indicator takes advantage of three emission-line ratios, which is an improvement over previous strong-line methods.

In principle, given the wavelength coverage of MUSE, the N2 and O3N2 diagnostics, calibrated by Marino et al. (2013), are also applicable. As pointed out by Kreckel et al. (2019), using O3N2 (or N2) and Scal can result in qualitatively different results from the MUSE data. In this work, we prefer to use metallicity indicators of N2S2H α and Scal, rather than N2 and O3N2. Actually, Marino et al. (2013) calibrated the O3N2 and N2 diagnostics to the T_e-based method and found that the O3N2 and N2 diagnostics result in uncertainties in oxygen abundance of 0.18 and 0.16 dex, respectively. Given the similar ranges of the metallicity determined by Scal, O3N2, and N2 for a given data set, the smaller uncertainty of the Scal diagnostic indicates a significant improvement over the previous O3N2 and N2 diagnostics. This improvement may come from the fact that (1) Scal use more emission-line ratios and (2) these break the degeneracy between metallicity and ionization parameter (see also Kreckel et al. 2019). The latter is also true for the N2S2H α indicator.

3.3.3. The Contamination from Diffuse Ionized Gas

In Section 2, we investigated the behavior of the SFR(t), the cold gas mass M_gas(t), and the gas-phase metallicity Z(t) in the gas-regulator system in response to variations in the inflow rate Φ(t) and SFE(t), as well as how this response depends on the (assumed constant) wind mass-loading factor λ and metallicity of the inflowing gas Z₀. Specifically, we found a negative correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ (i.e., $\mathrm{log}(\mathrm{SFR}/\langle \mathrm{SFR}\rangle$) vs. $\mathrm{log}(Z/\langle Z\rangle )$) when driving the gas regulator with time-varying inflow rate, and a positive correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ when driving with time-varying SFE(t). Therefore, one can in principle identify the driving mechanism of star formation activity by looking at the sign of the correlation between SFR and gas-phase metallicity in observational data.

However, as pointed out above in Section 2.4, one should look at the correlations of the relative values of SFR and Z (i.e., the residuals ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ in Figures 1 and 3), rather than the absolute values, in order to take out the effects of different 〈SFR〉 and 〈Z〉 for different galaxies or of different regions within them, e.g., the overall mass–metallicity or mass–sSFR relations or radial gradients in metallicity or sSFR within galaxies.

In this section, we will therefore construct radial profiles of sSFR and log(O/H) for the MAD galaxies and use these to construct localized ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) data points from the observations.

We here clarify why the sSFR is better than the SFR to characterize variations in the SFR within and between galaxies. Ideally, we would wish, in order to compare observations with the model in Section 2, to follow a given galaxy, or a region within a galaxy, as it changes temporally. This is, of course, not possible. Instead, we must compare different galaxies in the population, or different regions at the same radius within a galaxy, and, invoking ergodicity, assume that these reflect the temporal variations that we are interested in. For each observational point, we therefore need the best estimate of the average (long-term) state of that location. It is an observational fact that, both from galaxy to galaxy and within galaxies, the range of sSFR is smaller than the range of SFR. For this reason, the average sSFR will be better defined than the average SFR. Related to this, normalizing by the local Σ_*, even for spaxels at the same galactic radius, also removes the effect of azimuthal variations in the density of gas to the extent that both gas and stars vary in the same way.

We emphasize again that, in the model predictions, the ${\rm{\Delta }}\mathrm{log}$ SFR is replaceable by ${\rm{\Delta }}\mathrm{log}$ sSFR, and the ${\rm{\Delta }}\mathrm{log}Z$ is replaceable by ${\rm{\Delta }}\mathrm{log}$(O/H), as mentioned in the end of Section 2.

Figure 6 shows an example of the radial profiles and two-dimensional maps of sSFR and 12 + log(O/H) for one individual MAD galaxy, NGC 1483. This typical galaxy is not specially selected in any way and is shown for illustration purposes as being representative of the general sample. The top three panels of Figure 6 show the sSFR(r) and the 12 + log(O/H) as estimated by N2S2H α and Scal as a function of radius for individual spaxels. The log(O/H) based on the N2S2H α diagnostic is denoted as log(O/H)-DOP16, and the log(O/H) based on the Scal approach is denoted as log(O/H)-PG16. The radius used here is the deprojected radius scaled by the effective radius of the galaxy. In computing this deprojected radius, we use the disk inclination based on the measured minor-to-major-axis ratio and the position angle taken from the S4G photometry (Salo et al. 2015), assuming an infinitely thin disk. In each of the top panels, the blue line is a running median of 201 spaxels.

Zoom In Zoom Out Reset image size

As shown, for NGC 1483 the distribution of the sSFR at a given radius is quite strongly asymmetric, over nearly the whole range of galactic radius. While the sSFRs of most spaxels are close to the median profile (or slightly less), a small fraction of spaxels have sSFR that is enhanced by up to an order of magnitude relative to the median profile. This is due to the fact that star formation activity is not uniform across the disk, but happens mostly in spiral arms or other star formation regions. The regions with strong or enhanced sSFR can clearly be seen in the sSFR map in the middle left panel of Figure 6. In addition, the sSFR profile shows a positive radial gradient, which is consistent with the inside-out growth expected in disk galaxies (e.g., Pérez et al. 2013; Li et al. 2015; Ibarra-Medel et al. 2016; Lilly & Carollo 2016; Goddard et al. 2017; Rowlands et al. 2018; Wang et al. 2018b).

The impression of strong asymmetry is reduced in 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$, for both metallicity indicators. In addition, the overall 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$ profiles for both indicators have a negative radial gradient, consistent with the previous studies of disk galaxies (e.g., Pilkington et al. 2012; Belfiore et al. 2017; Erroz-Ferrer et al. 2019). This feature also can be seen in the maps of 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$, shown in the middle row of Figure 6. The measurements of $\mathrm{log}({\rm{O}}/{\rm{H}})$ based on N2S2H α and Scal are qualitatively consistent. However, for a given data set, $\mathrm{log}({\rm{O}}/{\rm{H}})$-DOP16 is usually larger than $\mathrm{log}({\rm{O}}/{\rm{H}})$-PG16, and the range of $\mathrm{log}({\rm{O}}/{\rm{H}})$-DOP16 is nearly twice that of $\mathrm{log}({\rm{O}}/{\rm{H}})$-PG16. We note that the particular galaxy is typical of the sample, and most of the SF disk galaxies show similar features in sSFR and $\mathrm{log}({\rm{O}}/{\rm{H}})$ to those shown for this galaxy.

As pointed out in Section 2.4, the gas regulator predicts that the average SFR will only depend on the average inflow rate Φ₀ and mass-loading factor λ, and the average Z is determined by the effective yield, defined from the wind mass-loading as y(1 − R + λ)⁻¹, and the metallicity of the inflowing gas Z₀. However, the λ (and possibly also the Z₀) may well change radially within individual galaxies, because λ should be directly related to the gravitational potential well. If so, the fitted average $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ profile should therefore reflect the radial dependence of the wind mass loading and/or Z₀ (see further discussion in Section 4.2.2). We are not interested in this effect if it is caused by time-invariant factors (λ and Z₀). Rather, we are interested in the spaxel-by-spaxel variations that are superposed on this smooth radial profile, because these are presumably caused by shorter-term temporal variations.

Therefore, we focus on the values of individual spaxels relative to these underlying trends, i.e., on the ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) residuals when the underlying sSFR(r) and $\mathrm{log}({\rm{O}}/{\rm{H}})(r)$ profiles are subtracted from each spaxel (compare to the relation of ${\rm{\Delta }}\mathrm{log}$ SFR vs. ${\rm{\Delta }}\mathrm{log}Z$ in Figures 1 and 3). In this way, we are effectively not allowing a variation of λ and Z₀ between spaxels at a given radius for an individual galaxy, but allowing these quantities to vary radially, and removing the effect of this from our analysis.

To achieve this, we first perform a linear fit to the sSFR(r) and 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})(r)$ profiles based on all the individual spaxels. These fits are shown as red lines in the top panels of Figure 6. As shown, for both sSFR and 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$, the linear fit is quite a good representation of the median profile, although it is not perfect. We then define the ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ for each individual spaxel as the deviation of each spaxel from the fitted profile of sSFR(r) and 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})(r)$, respectively. In this way, we can eliminate the overall radial dependences of sSFR and $\mathrm{log}({\rm{O}}/{\rm{H}})$, as well as global differences between galaxies, such as the overall sSFR or effects due to the mass–metallicity relation. As shown in Section 2, in the gas-regulator framework, these changes in 〈SFR〉 or 〈Z〉 will reflect differences (radially within galaxies or from galaxy to galaxy) in the overall inflow rate, mass-loading factor, and Z₀.

The bottom three panels of Figure 6 show the maps of ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ that are obtained for NGC 1483 after removing the radial gradients in $\mathrm{log}$ sSFR and $\mathrm{log}({\rm{O}}/{\rm{H}})$. It is immediately apparent that regions with enhanced SFR, indicated by red bumps, nearly always show enhanced metallicity for both of the metallicity indicators. This is consistent with the previous analysis of Kreckel et al. (2019) and Erroz-Ferrer et al. (2019), at a similar spatial resolution of ∼100 pc. It should be noted that the color scale used for ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$-DOP16 has twice the range of that used for ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$-PG16.

Figure 7 shows the fitted linear profiles of sSFR, 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$-DOP16 and 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$-PG16 for all 38 MAD galaxies. In displaying these profiles, we separate galaxies into three mass bins: $\mathrm{log}({\text{}}{M}_{* }/{\text{}}{M}_{\odot })\lt 10.0$ (blue lines), $10.0\lt \mathrm{log}({\text{}}{M}_{* }/{\text{}}{M}_{\odot })\lt 10.8$ (green lines), and $10.8\lt \mathrm{log}({\text{}}{M}_{* }/{\text{}}{M}_{\odot })$ (red lines). As shown, almost all of the MAD galaxies have positive radial gradients in sSFR with only four exceptions. It can be seen that there is no strong dependence of the sSFR profile on the stellar mass of galaxies. While the profiles in 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$ have similar slopes, the overall values show a very strong dependence on global stellar mass, reflecting the MZR.

Zoom In Zoom Out Reset image size

The definition of ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ for the individual spaxels, based on Figure 7, enables us to investigate the correlations of ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ for small-scale regions (∼100 pc) within individual MAD galaxies. However, one can well imagine that the physical processes driving the small-scale star formation may be very different from those driving the star formation on galactic scales, and this motivates us to define analogous quantities, ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$, to reflect the global properties of galaxies with the MAD population. For this purpose, we choose the fitted values of sSFR and 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$ at 0.5 R_e (see the red lines in the top panels of Figure 6) as representative of the global properties of individual galaxies, ² because the spatial coverage for MAD galaxies generally extends to at least 0.5 R_e.

Figure 8 shows the sSFR${}_{0.5\mathrm{Re}}$, 12 + $\mathrm{log}{({\rm{O}}/{\rm{H}})}_{0.5\mathrm{Re}}$-DOP16, and 12 + $\mathrm{log}{({\rm{O}}/{\rm{H}})}_{0.5\mathrm{Re}}$-PG16 as a function of the overall stellar mass of the galaxies. The overall stellar mass is obtained by broadband spectral energy distribution fitting (Erroz-Ferrer et al. 2019). The sSFR${}_{0.5\mathrm{Re}}$ decreases slightly with stellar mass, and 12 + $\mathrm{log}{({\rm{O}}/{\rm{H}})}_{0.5\mathrm{Re}}$ increases significantly with stellar mass, for both metallicity indicators. Both of these trends are, of course, well established for SF galaxies in the literature. As pointed out above, in the framework of the gas-regulator model, the dependence of SFR and Z on stellar mass is due to the stellar mass dependence of the inflow rate, λ, and Z₀ (see Equations (22) and (23) and discussion in Lilly et al. 2013). To eliminate the mass dependence, we perform a linear fit for each of these two relations, as shown with the black lines in Figure 8. In a similar way, for each individual MAD galaxy, we can then define the ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{0.5\mathrm{Re}}$ or ${\rm{\Delta }}\mathrm{log}{({\rm{O}}/{\rm{H}})}_{0.5\mathrm{Re}}$ to be the deviation of an individual galaxy from the linearly fitted relation. This is useful to study the driving mechanisms of star formation from galaxy to galaxy within the population.

Zoom In Zoom Out Reset image size

In Appendix C, we repeat our basic analysis by defining ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ for individual spaxels using average relations with Σ_* rather than with galactic radius, and we find that the basic results remain the same.

To compare with our theoretical expectations of the gas-regulator model constructed in Section 4.1, we have defined in the previous section ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) both on the 100 pc scale of individual spaxels in MAD galaxies and for the much larger galactic scale of MAD galaxies as a whole (defined at 0.5 R_e). In this section, we will explore the correlations between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) on these scales and interpret the results in terms of the predictions of the gas-regulator model.

4.2.1. 100 pc Scales

We first look at the ${\rm{\Delta }}\mathrm{log}$ sSFR−${\rm{\Delta }}\mathrm{log}$(O/H) relation on scales of ∼100 pc. Panel (a) of Figure 9 shows the distribution of all valid MAD spaxels on the ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H)-DOP16 diagram. The gray scale shows the number density of spaxels in logarithmic space. In constructing panel (a), we give each MAD galaxy the same weight, by weighting the contribution of spaxels from each galaxy in panel (a) of Figure 9. In other words, for a given MAD galaxy, we weight the individual spaxels of that galaxy by ${N}_{\mathrm{spaxel}}^{-1}$, where N_spaxel is the total number of valid spaxels for that particular galaxy. This ensures that the sum of the weights of all the spaxels in each galaxy is the same, and therefore that each galaxy is equally represented in the figure.

Zoom In Zoom Out Reset image size

As a whole, we find a significant positive correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H)-DOP16 for all the individual spaxels in the 38 MAD galaxies. Furthermore, we also investigate the correlation of ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H) for each individual MAD galaxy. Panel (b) of Figure 9 shows bisector fits of the ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H)-DOP16 relation of individual spaxels for each of the 38 MAD galaxies. Here we adopt a bisector fitting (Isobe et al. 1990) because there is no reason for us to prefer regression of ${\rm{\Delta }}\mathrm{log}$(O/H) on ${\rm{\Delta }}\mathrm{log}$ sSFR or ${\rm{\Delta }}\mathrm{log}$ sSFR on ${\rm{\Delta }}\mathrm{log}$(O/H). As can be seen, consistent with the result of panel (a), ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H)-DOP16 show positive correlation for 37 of the 38 MAD galaxies with only one exception. The exception is NGC 4030, the most massive galaxy in MAD sample. Inspection of the color-coding suggests that this result is not dependent on the mass of the galaxy.

The same analysis is repeated for log(O/H)-PG16, in panels (a) and (b) of Figure 10. Overall, we find that ${\rm{\Delta }}\mathrm{log}$ sSFR and Δlog(O/H)-PG16 still show a positive correlation, although not as significant as in panel (a) of Figure 9. Consistent with this, panel (b) shows that 32 of the MAD galaxies show positive correlations of ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H)-PG16, but now six MAD galaxies show negative correlations. We note that when using O3N2 the result is qualitatively different from that when using Scal (or N2S2H α) for galaxies with stellar mass below ∼10^10.5 M_⊙. This is likely due to the fact that the O3N2 indicator has a larger uncertainty (∼0.18 dex) than Scal (∼0.05 dex) in determining the metallicity (see more discussion in Sections 3.3.2 and 6.3).

Zoom In Zoom Out Reset image size

For comparison purposes, we also show the model predictions of ${\rm{\Delta }}\mathrm{log}$ SFR versus ${\rm{\Delta }}\mathrm{log}Z$ on panels (a) of Figures 9 and 10. The two white roughly elliptical shapes are based on the model predictions shown in the top middle panel of Figure 3. As discussed in Section 2, the correlation of ${\rm{\Delta }}\mathrm{log}Z$ versus ${\rm{\Delta }}\mathrm{log}$ SFR depends on both ξ and Z₀/y_eff. For simplicity, we only consider here the model predictions for Z₀ = 0. We then choose the $\mathrm{log}\xi$ (i.e., the ellipse in Figure 3) that has the same ratio of the dispersion in ${\rm{\Delta }}\mathrm{log}Z$ (x-axis) and ${\rm{\Delta }}\mathrm{log}$ SFR (y-axis) as the observed data points in panels (a), and then we scale this model ellipse so as to match the 1σ and 2σ dispersions of the observed quantities (see further discussion in Section 4.3). Given the simplicity of the models in Section 2, and especially given the different ranges of metallicity returned by the two metallicity indicators, we do not believe that a quantitative fitting would be meaningful. Nominally, we get $\mathrm{log}\xi =0.6$ for log(O/H)-DOP16 and $\mathrm{log}\xi =0.96$ for log(O/H)-PG16. Within the limitations of model and metallicity estimators, this qualitative matching approach suggests that the prediction of time-varying SFE in the gas-regulator model is in good agreement with the data in predicting a positive correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) at 100 pc scale.

We emphasize that, given the highly idealized nature of the model, any precise comparison of the model with the data would not be very useful. Rather, the model was intended to establish two things: the sign of the correlation between changes in the SFR and Z (positive or negative), and the approximate relative variation of these two quantities. Throughout the paper, the comparisons between the model and the data should be treated in this spirit.

To try to quantify the significance of this positive correlation, we calculate the distribution of Pearson correlation coefficients for the 38 MAD galaxies in Figure 11. For each individual MAD galaxy, the coefficient is calculated in two ways, both by weighting the individual spaxels equally and by weighting them by the SFR of the spaxels. There are two reasons for doing the SFR weighting. First, it ensures that the inner parts of galaxies, which have fewer but brighter pixels, are not swamped by the outer parts. Second, it is clear in Figure 3 that while the sinusoidal variation in SFE produces symmetric changes in SFR and metallicity, the (possibly more realistic) step function changes produce asymmetric changes, in which the positive correlation is most clearly established by the regions of highest SFR. As shown, for both approaches of computing the correlation coefficient and for both metallicity indicators, we find that most of the coefficients are positive, with only a few less than zero. This is consistent with the results in panels (a) and (b) of Figure 9 and Figure 10. The correlation of ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H) becomes more significant when weighting the spaxels with their SFR. This is due to the fact that regions with strongly enhanced star formation always show enhanced gas-phase metallicity (see also panel (a) of Figures 9 and 10). Both the positive correlation of ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H) and the increased significance of the correlation for regions of enhanced star formation are also visible to the eye in Figure 6. This may reflect the point made in the context of Figure 3 in Section 2.3, that the positive correlations caused by step function changes to the SFE are most clearly seen when the SFR is highest.

Zoom In Zoom Out Reset image size

We have also examined the significance of the correlation based on the p-value, ³ which is the probability of obtaining the real observed result under the assumption of no correlation. A smaller p-value means that the correlation is more significant. For both metallicity indicators, we find that the p-values are very close to zero (p < 0.001) for all except one or two (for different metallicity indicators) of the 38 MAD galaxies. This test shows that the positive correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) is present also within individual galaxies, as well as in the ensemble analysis.

In the gas-regulator framework, we showed that ${\rm{\Delta }}\mathrm{log}$ SFR (or ${\rm{\Delta }}\mathrm{log}$ sSFR) and ${\rm{\Delta }}\mathrm{log}Z$ will be positively correlated when the gas-regulator system is driven by a time-varying SFE (see Figure 3). Comparing the model predictions with panel (a) of Figure 9 or Figure 10, we can conclude that at ∼100 pc scales within galaxies the variation of SFR (and Z) is due to a time-varying SFE experienced by a particular packet of gas. This conclusion implies an assumption that different spatial regions of galaxies, around an annulus of given galactic radius, constitute different temporal phases of a gas packet, as modeled in Figure 3. This assumption is not unreasonable and is supported by other strong observational evidence in favor of a time-varying SFE at ∼100 pc scale.

Leroy et al. (2013) and Kreckel et al. (2018) have found that the dispersion of the SFE based on molecular gas measurements increases significantly toward smaller scales (see also Kruijssen & Longmore 2014). Specifically, the scatter of the SFE is ∼0.4 dex at ∼100 pc (Kreckel et al. 2018). Consistent with this, Kruijssen et al. (2019) and Chevance et al. (2020) have shown that molecular gas mass and SFRs are spatially decorrelated on the scales of individual GMCs in nearby disk galaxies, contrary to the tight correlation seen on kiloparsec and galactic scales (e.g., Shi et al. 2011; Bigiel et al. 2008). This decorrelation implies rapid temporal cycling between GMCs and star formation, likely involving feedback processes. By using this decorrelation, Chevance et al. (2020) constrain the properties of GMCs and claimed that the GMC lifetimes are typically 1030 Myr, which consists of a long inert phase without massive star formation and a short phase (of less than 5 Myr) with a burst of star formation. These observational results indicate a strong spatial variation of SFE at ∼100 pc scales, simultaneously suggesting a strong temporal variation of SFE for a given packet of gas.

In Section 2.3, we constructed two models of time-varying SFE(t) in the form of both a sinusoidal function and a periodic step function. According to the above discussion, we find that the SFE(t) at ∼100 pc scale may be better characterized by a step function with a short duration of top phase and a long duration of bottom phase, rather than the sinusoidal model. Consistent with this, we find that the model prediction with the step function is in good qualitative agreement with the observational results. Specifically, the distribution of ${\rm{\Delta }}\mathrm{log}$ sSFR is strongly asymmetric, with a very small fraction of spaxels showing strongly enhanced SFR, and with these regions also showing strongly enhanced metallicity, as found in panel (a) of both Figures 9 and 10. These features can also be found in the model prediction shown in the bottom middle panel of Figure 3.

It goes without saying that the models explored in Section 2 are simple heuristic models, which cannot be expected to explain all the details of the situation. Not least, timescales of star formation of only 10–30 Myr (Chevance et al. 2020) are comparable to the timescales for chemical enrichment, an effect neglected by our use of the instantaneous recycling approximation in Section 2. Nevertheless, we can conclude qualitatively that the variation of SFR and gas-phase metallicity on 100 pc scales is primarily due to a time-varying SFE experienced by the gas.

4.2.2. Subgalactic Scales

4.2.3. Galactic Scales

We showed in Section 2 that the relative strength of variations in the SFR and metallicity of the gas regulator should depend on the response timescale of the system, set by the gas depletion timescale, relative to the timescales of any driving variations in the inflow or in the SFE. Equation (20) describes the relative amplitude of these variations, characterized by $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogSFR) across the population, as a function of ξ when driving the gas-regulator system with sinusoidal SFE(t). According to Equation (20), the $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogSFR) should decrease with increasing ξ, i.e., decrease with increasing effective gas depletion time if we ignore the possible variation of the driving period of SFE(t) (see Figure 5).

We have therefore calculated the dispersions across the spaxels of the residual quantities ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ and ΔlogsSFR that were plotted in the two panels (a) in Figures 9 and 10 for the two metallicity indicators, respectively. We calculate this dispersion for each of the MAD galaxies, respectively.

Figure 12 shows the resulting ratios of $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$ to σ(ΔlogsSFR) for each of the 38 MAD galaxies for the two metallicity indicators. It is obvious that the $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$ based on N2S2H α is overall greater than that based on Scal. This is not a noise issue, but is due to the fact that the range of log(O/H)-DOP16 is nearly twice the range of log(O/H)-PG16 for a given data set, as mentioned earlier in this paper. This systematic uncertainty hinders the quantitative interpretation of these dispersions, although trends established within a single estimator (i.e., within a single panel of Figure 12) should have some validity.

Zoom In Zoom Out Reset image size

As can be seen, the $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogsSFR) based on N2S2H α is in the range of 0.12–0.42, with a median value of 0.24. The $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogsSFR) based on Scal is about half this value, mostly in the range of 0.06–0.17 with a median value of 0.11. We do not find a significant dependence of $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogsSFR) on stellar mass, except possibly for an increase at the lowest stellar masses below 10⁹ M_⊙. The MAD galaxy with the largest $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogsSFR) is ESO499-G37. A small fraction of spaxels in ESO499-G37 have $\mathrm{log}$ N2 < −0.6, resulting in very low metallicity in these regions (see Equation (28)), and producing a large dispersion in ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$-PG16.

Although we suspect that a sinusoidally time-varying SFE is not likely to be realistic in the universe at ∼100 pc, Equation (20) (see also Figure 2) does permit a rough order-of-magnitude estimate of whether these relative dispersions are reasonable and of the approximate timescales involved.

If we examine what values of ξ in Equation (20) produce the observed $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/σ(ΔlogsSFR) for typical MAD galaxies, i.e., 0.24 and 0.12 for the two metallicity estimators, then these are $\mathrm{log}\xi =0.607$ for N2S2H α and $\mathrm{log}\xi =0.956$ for Scal at Z₀ = 0. If we take 1 Gyr as a reasonable estimate for the effective gas depletion timescale for the overall galaxy population (see W19), we get rough estimates of T_p = 1.5 Gyr for N2S2H α and T_p = 0.7 Gyr for Scal as the nominal period of a time-varying SFE.

Intriguingly, in the Milky Way, a periodic star formation history with a period of ∼0.5 Gyr has been suggested in the solar neighborhood, from analysis of the resolved stellar population (Hernandez et al. 2000; de la Fuente Marcos & de la Fuente Marcos 2004). As further pointed out by Egusa et al. (2009), this periodic star formation history may be associated with the passage of the spiral potential in the density wave theory (Lin et al. 1969). Assuming that the potential has a two-armed pattern as suggested by Drimmel (2000), the pattern speed can be calculated as Ω_P = Ω(r = R_⊙) − π/(0.5 Gyr) = 21 km s⁻¹ kpc⁻¹. This is also consistent with the result from numerical simulations of the stellar and gaseous response to the spiral potential, presented by Martos et al. (2004).

It is quite suggestive that this same sort of timescale emerges in our own analysis of the metallicities in terms of a periodically varying SFE, and it suggests that the periodic (or time-varying) SFE(t) that is relevant on 100 pc scales may be explained by the passage of orbiting packets of gas through the spiral density wave. It should be noted in the context of the gas regulator that changes in the density of the gas due to local flows associated with the passage of a spiral density wave have no effect, except (quite possibly) to change the SFE because of the changed density of gas, and should not be considered as "inflows" or "outflows," which are explicitly concerned with changes in the total (baryonic) mass of the gas-regulator "system" being considered. Both the metallicity and the sSFR are formally intrinsic quantities that make no statement about the size of the system, i.e., whether the boundary of the reservoir is of fixed physical size or not. Put another way, neither the measurements of metallicity nor the sSFR will change (per se) owing to a compression of the whole system of gas and stars (unless there are associated physical changes, e.g., in the SFE).

However, given the many steps and uncertainties in our analysis, we certainly caution against overinterpretation of our analysis.

In Section 4.2.3 we examined the relation between the measure of the overall, i.e., "galactic-scale," ${\rm{\Delta }}\mathrm{log}$ sSFR−${\rm{\Delta }}\mathrm{log}$(O/H) for MAD galaxies, taking as measures of these quantities the values from the radial profiles at a fiducial radius of 0.5 R_e, as shown for the 38 MAD galaxies in the two panels (e) in Figures 9 and 10. Since the MaNGA coverage usually extends farther than 1.5 R_e for individual galaxies, we can now measure the integrated sSFR and metallicity within 1.5 R_e. The metallicity within 1.5 R_e is computed as the Hα luminosity-weighted 12 + log(O/H) for all the spaxels within 1.5 R_e. This is probably more representative of the global quantities than the sSFR (or metallicity) that was measured in the MAD sample at one particular radius.

As before, we first construct the sSFR${}_{\lt 1.5\mathrm{Re}}$ versus mass and log(O/H)${}_{\lt 1.5\mathrm{Re}}$ versus mass relations, and we use these to normalize the measurements of individual galaxies. The first is obtained with a linear fit of the sSFR${}_{\lt 1.5\mathrm{Re}}$ versus stellar mass relation. For the metallicity, we do a polynomial fit to the log(O/H)${}_{\lt 1.5\mathrm{Re}}$ versus stellar mass relation, since it clearly flattens at the high-mass end. For each individual MaNGA galaxy, we thereby define the ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt 1.5\mathrm{Re}}$ and ${\rm{\Delta }}\mathrm{log}$(O/H)${}_{\lt 1.5\mathrm{Re}}$ to be the (logarithmic) deviation from these relations.

The correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt 1.5\mathrm{Re}}$ and ${\rm{\Delta }}\mathrm{log}$(O/H)${}_{\lt 1.5\mathrm{Re}}$ for the MaNGA sample is shown for the two metallicity indicators in the two panels of Figure 13. As shown, enlarging the sample size by a factor of about 25, the negative correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) is very clearly seen for both metallicity indicators. The Pearson correlation coefficients for the MaNGA sample are −0.23 and −0.36 for N2S2H α and Scal indicators, respectively. The p-values are zero (<10⁻¹⁰) for both metallicity indicators, substantially strengthening the results found in panels (e) of Figures 9 and 10 for the much smaller MAD sample.

Zoom In Zoom Out Reset image size

For each metallicity indicator, the linear slopes obtained in MaNGA by bisector fitting are very similar to those seen in the MAD sample (panels (e) of Figures 9 and 10). However, we note again that the slopes for the two metallicity indicators are significantly different, due to the range of log(O/H)-DOP16 being nearly twice of that of log(O/H)-PG16. As also noted above, this clear inverse correlation is a restatement of the existence of SFR as a (negative) second parameter in the mass–metallicity relation. For comparison, we show the model predictions for a sinusoidally time-varying inflow rate as black curves scaled to the MaNGA data using the same procedure as described in Section 4.2.1. These same curves were plotted in panels (e) of Figures 9 and 10 for comparison with the MAD data. As can be seen, a single model prediction is broadly consistent with both the observed MaNGA and MAD data points.

We now turn to look more closely at the radial variations of SFR and log(O/H) in galaxies and how these vary across the galaxy population, using again the MaNGA sample. This analysis may be regarded as an extension of the analysis of SFR profiles that was presented in W19.

In W19, we investigated the dependence of the normalized star formation surface density profiles, Σ_SFR(R/R_e), of star-forming MaNGA galaxies on the overall stellar mass and on ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}$, the vertical deviation of the galaxy from the "nominal" SFMS. ⁴ We showed (Figure 5 of that paper) that star-forming MaNGA galaxies that lie significantly above (or below) the overall SFMS show an elevation (or suppression) of SFR at all radii. In addition, we showed that whereas at low stellar masses this elevation (suppression) of star formation is more or less uniform with radius, for the more massive galaxies the elevation (or suppression) of star formation becomes more pronounced in the central regions of the galaxies. As a direct consequence of this, the dispersion in the (normalized) Σ_SFR across the galaxy population, which we designated $\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$, was found to be correlated with the local surface mass density, which is equivalent to an inverse correlation with the gas depletion timescale in the extended Schmidt law (Shi et al. 2011), since in that local surface density and the gas depletion timescale are directly related. This result was shown in Figure 9 of W19.

In this subsection, we will carry out an analogous study of the metallicity profiles for the same MaNGA galaxies used in W19. We will investigate the dependence of these metallicity profiles on stellar mass and the deviation from the SFMS, analogously to the analysis of Σ_SFR in that paper.

For consistency with that earlier work, we use ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}$ to separate galaxies. Galaxies are classified into four subsamples, $0.33\lt {\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}$, $0.0\lt {\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}\lt 0.33$, −0.33 < ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}\lt 0.0$, and ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}$<−0.33, and we further split the sample into five bins of overall stellar mass, as in W19. For each individual galaxy, we first compute a radial profile of 12 + log(O/H), determined as the median 12 + log(O/H) of spaxels located within each radial bin. In each of these subsamples we then construct the median 12 + log(O/H)(r/R_e) radial profiles, using the two metallicity estimators in turn, and estimating uncertainties by bootstrapping the sample galaxies.

Figure 14 shows the median 12 + log(O/H)-DOP16 profiles for these subsamples of galaxies. In each stellar mass bin, the metallicity profiles are displayed in blue, green, yellow, and red in descending order of their overall ${\rm{\Delta }}\mathrm{log}$ sSFR${}_{\lt \mathrm{Re}}$. Figure 15 is the same as Figure 14, but for log(O/H)-PG16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first-order result is clear: independent of which metallicity indicator is used, low-mass galaxies that lie significantly above (or below) the SFMS in their overall sSFR have systematically lower (or higher) log(O/H) over the whole galactic radii. This dependence of log(O/H) on ${\rm{\Delta }}\mathrm{log}$ sSFR, however, decreases (and even vanishes and possibly reverses) with increasing stellar mass. There is some evidence that it also decreases to the centers of the galaxies (see, e.g., the higher-mass bins in Figure 14).

The result of Figure 14 is broadly consistent with the result of Figure 15, except for the highest-mass bins. In the highest-mass bin of Figure 15, a positive correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) can be seen, which is clearly opposite to the result of other mass bins. This might be due to the failure of the N2S2H α indicator that the assumed N/O–O/H relation does not hold for the most massive SF galaxies.

The overall negative correlation between the overall ${\rm{\Delta }}\mathrm{log}$ sSFR and log(O/H) shown by the sets of profiles in Figures 14 and 15 is another manifestation of the inverse correlations in panels (e) of Figures 9 and 10 and in Figure 13. We can then argue that these indicate that time-varying inflows are the primary drivers of variations of sSFR and log(O/H), across the galaxy population (see also W19; Wang & Lilly 2020b). As noted above, this result is a manifestation of the general presence of SFR as a (negative) second parameter in the overall mass–metallicity relation (e.g., Mannucci et al. 2010). The fact that we see the range of log(O/H) (at given R/R_e) decreasing with stellar mass is also consistent with previous studies of the overall Z(M_*, SFR) relation (e.g., Mannucci et al. 2010; Curti et al. 2020).

In Figure 5 of W19, the dispersion in the (normalized) Σ_SFR increases slightly with increasing stellar mass and also increases toward the centers of galaxies for galaxies of a given stellar mass. It is quite striking how the dispersion of $\mathrm{log}({\rm{O}}/{\rm{H}})$ (at a given mass and radius) behaves in the opposite way to the dispersion in (normalized) Σ_SFR shown in our previous work. Whereas the former decreases with increasing stellar mass (and possibly toward the centers of galaxies), the dispersion of Σ_SFR (or sSFR) increases slightly with mass and increases toward the centers of galaxies. We will discuss this in detail in the next section.

As discussed in Section 2, the simple gas-regulator model predicts not only the sign of the correlation between ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ but also the variation of ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}Z$ as a function of ξ (see Figure 2), defined as the ratio of the driving period to the effective gas depletion timescale in Equation (7). In this subsection, we will investigate more quantitatively the variation of ${\rm{\Delta }}\mathrm{log}$ SFR and ${\rm{\Delta }}\mathrm{log}$(O/H) in MaNGA galaxies. This will inevitably run up against the systematic uncertainties arising from the significantly different outputs of our two chosen metallicity indicators. We can, however, try to minimize the effects of these by looking at relative effects across the galaxy population, expecting that systematic effects will thereby cancel out.

In W19, we constructed the Σ_SFR profiles for MaNGA galaxies at five different stellar mass bins and defined the parameter ${\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$ as the deviation of a given galaxy from the median Σ_SFR profile at a given galactic radius and in a given stellar mass bin. We then computed the dispersion of this quantity across the population within five stellar mass bins and within five radial bins of r/R_e. It was found that the scatter of ${\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$, which we denoted as $\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$, increases significantly with stellar surface mass density, Σ_*. We interpreted this trend in terms of a decreasing gas depletion time, linking the stellar surface mass density to the gas depletion timescale via the extended Schmidt law (Shi et al. 2011), i.e., ${\tau }_{\mathrm{dep}}\propto {{\rm{\Sigma }}}_{* }^{-1/2}$. The trend with the inferred gas depletion time was found to be consistent with the model prediction of the time-varying inflow-rate-driven scenario, provided that the driving period of the inflow was more or less the same for all galaxies. In this subsection, we look at the analogous variation of ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ in the population and compare the result with the model predictions.

Similar to the definition of ${\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$, we define the ${\rm{\Delta }}\mathrm{log}$(O/H) to be the deviation from the median $\mathrm{log}$(O/H) at a given specific radius and a given stellar mass bin. We then compute the dispersion of this quantity, σ(${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$), and compare this to the dispersion of the completely independent (normalized) star formation surface density, $\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$, used in W19.

Figure 16 shows the ratio of σ(${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$) to σ(${\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$) as a function of the inferred gas depletion time, for both metallicity indicators. The gas depletion time is again derived from the Σ_* on the basis of the extended Schmidt law (Shi et al. 2011). To be in line with W19, we here use the scatter of ${\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}}$, rather than the scatter of ΔlogsSFR. We have checked that using the latter gives a consistent result to that in Figure 16. In each panel of Figure 16, different colors are used for different stellar mass bins. The different data points of the same color represent the different radial bins, evenly spaced with a bin width of 0.2 R_e. At a fixed stellar mass bin, the τ_dep monotonically increases with galactic radius as Σ_* decreases. As in W19, data points at galactic radii larger than R_e are indicated in gray, since these may be more easily affected by environmental effects (see also W19; Wang & Lilly 2020b). Readers are invited to consult Figure 9 in W19 and Figures 15 and 17 in Wang & Lilly (2020b) for further insights into the variations of Σ_SFR.

Zoom In Zoom Out Reset image size

As a whole, the ratio $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$/$\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$ increases quite sharply with the inferred gas depletion time for both metallicity indicators. This is individually true for each the five stellar mass bins (where it reflects radial changes), except for the highest stellar mass bin for log(O/H)-DOP16, and it is true when comparing galaxies of different stellar masses. These trends reflect quantitatively a combination of the effects that are seen in Figures 14 and 15 of this paper and Figure 5 of W19. It should, however, be noted that the dispersions σ are calculated using the individual galaxies, whereas these figures plot the median values within the four bins of ${\rm{\Delta }}\mathrm{log}$ SFR and so are not directly comparable. We have discussed the particular case of the highest-mass bin for log(O/H)-DOP16 in Section 5.2.

It is again clear that the dispersions $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ obtained using the different metallicity estimators differ by a factor of two, reflecting their twofold difference in range within the sample. In the absence of a reliable reason to prefer one over the other, this makes any precise quantitative comparison with the predictions of the simple gas-regulator model (see Figure 2) impossible. However, three points may be made, independent of the choice of estimator.

First, both metallicity estimators show about a factor of four increase in the ratio $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$/$\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$ as the inferred gas depletion timescale increases by an order of magnitude from $\mathrm{log}{\tau }_{\mathrm{dep}}$ ∼ 8.8 to 9.8. This is as expected for variations in inflow rate (around $\mathrm{log}\xi \sim 0$) and quite opposite to the trend expected for variations in SFE, which would have predicted a decrease in this ratio as the gas depletion timescale increases.

By horizontally shifting the model prediction in the left panel of Figure 2 (according to Equation (15)) so as to match the result of Figure 16 (see the black lines), we get (1 + λ)T_p/2π ∼ 10^9.5 yr for log(O/H)-DOP16 and ∼10^9.8 yr for log(O/H)-PG16. These are equivalent to driving periods of the inflow of T_p of a few to several gigayears. This is broadly consistent with the independent argument presented in Wang & Lilly (2020b) that the variation of ${\rm{\Delta }}\mathrm{log}$ sSFR is mainly produced by a variation of inflow rate on relatively long timescales. The driving period of changes in the inflow appears to be considerably longer than the period of the temporal variations in SFE discussed in Section 4.3.

Second, it can be seen that at a given gas depletion time (stellar surface mass density) more massive galaxies tend to have a lower value of $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$/$\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$. This could possibly reflect the quite plausible expectation that more massive galaxies might well have higher Z₀/y_eff than less massive galaxies as a result of either a lower wind mass loading λ or a higher inflow metallicity Z₀, leading to a reduction in $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}}))$/$\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$, as shown in Figure 2.

Finally, it is noticeable that, even with the larger range of log(O/H)-DOP16, the ratio of $\sigma ({\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$/$\sigma ({\rm{\Delta }}\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SFR}})$ never exceeds unity, the maximum value permitted by the gas-regulator model, as shown in Figure 2.

In this work, we find that on ∼100 pc scales within individual galaxies the local ${\rm{\Delta }}\mathrm{log}$(O/H) appears to be positively correlated with the local ${\rm{\Delta }}\mathrm{log}$ sSFR, while when looking into the integrated quantities of the same galaxies across the galaxy population, the ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})$ is found to be negatively correlated with ${\rm{\Delta }}\mathrm{log}$ sSFR. These results are quite consistent with previous findings, as discussed in Section 1. Specifically, based on the ∼1000 SAMI galaxies, Sánchez et al. (2019) found that ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) (defined in a similar way to that in the present work) show a negative correlation across the galaxy population for a wide range of metallicity indicators, with the correlation coefficient between −0.32 and −0.14. At highly resolved scale (∼100 pc), many authors have found that regions with strongly enhanced star formation show enhanced gas-phase metallicity (e.g., Ho et al. 2018; Erroz-Ferrer et al. 2019; Kreckel et al. 2019). Our results are consistent with these previous results.

In the context of our simple gas-regulator framework, the opposite sign of the correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) on 100 pc (GMC-scales) and on galactic scales and large subgalactic scales indicates that different physical processes regulate the star formation and chemical enhancement on these different scales. As a whole, a positive ${\rm{\Delta }}\mathrm{log}$ sSFR−${\rm{\Delta }}\mathrm{log}$(O/H) relation arises driving the gas regulator with time-varying SFE(t), and a negative ${\rm{\Delta }}\mathrm{log}$ sSFR–${\rm{\Delta }}\mathrm{log}$(O/H) relation is the result of driving it with a time-varying inflow rate. A time-varying SFE(t) at ∼100 pc scale (see Kruijssen et al. 2019; Chevance et al. 2020) and a time-varying inflow rate at galactic scale (see Wang et al. 2019; Wang & Lilly 2020a) are also suggested by other recent works.

In this work, we have not examined the intermediate scales of, say, 1 kpc. However, it is not difficult to infer that, as the scale is increased, the effect of time-varying SFE(t) becomes weaker and that of the time-varying inflow rate becomes stronger. This is likely the reason that at ∼1 kpc scale the correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) is weaker or even has disappeared, as seen by previous works (Moran et al. 2012; Barrera-Ballesteros et al. 2017; Berhane Teklu et al. 2020).

As shown in Equation (23), the metallicity of the steady state (i.e., constant inflow rate and constant SFE) is only determined by the metallicity of the inflow gas Z₀ and the effective yield including the effects of any wind, y(1 − R + λ)⁻¹. These two parameters are expected to be strongly correlated with the global stellar mass. The Z₀ may be expected to increase with stellar mass because the circumgalactic medium of galaxies was enriched by the outflow of enriched material driven by star formation in the past. The wind loading λ is expected to decrease with stellar mass because more massive galaxies have deeper gravitational potential wells. This is probably the origin of the observed mass–metallicity relation. In addition, we emphasize that Z does not depend on the absolute value of SFE or inflow rate, but depends on the change (with time) of it, under the gas-regulator model frame.

Based on the analysis of the present work (and also some previous works), the so-called FMR is clearly not valid at subkiloparsec scales. In the gas-regulator framework, we predict that the mass of the gas reservoir is always negatively correlated with metallicity (see Figures 1 and 3), regardless of the driving mechanisms by time-varying SFE or time-varying inflow rate. Observationally, Bothwell et al. (2013) and Bothwell et al. (2016) found that the cold gas mass appears to be more fundamental in determining the MZR than SFR, for both atomic and molecular gas, consistent with this picture. In this sense, the mass of the gas reservoir is a better secondary parameter than SFR in determining the metallicity across the GMC scale to galactic scale.

The importance of the SFR is that, in studying the correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) (as in this paper), we can distinguish the underlying physical processes that are driving variations in gas content, SFR, and metallicity, even for individual galaxies. Specifically, even though the negative ${\rm{\Delta }}\mathrm{log}$ sSFR versus ${\rm{\Delta }}\mathrm{log}$(O/H) relation across the galaxy population indicates that time-varying inflow is the main driver of variations in SFR, it is completely possible that in some particular galaxies existing reservoirs of cold gas are undergoing strong gravitational instability, leading to a temporary increase in their SFE.

For instance, Jin et al. (2016) identified 10 SF galaxies with kinematically misaligned gas and stars from the MaNGA survey. These galaxies have intense ongoing star formation and high gas-phase metallicity in their central regions with respect to normal SF galaxies, which can be easily interpreted as evidence of a temporary increase in SFE due to gas−gas collision between the preexisting gas and the misaligned inflowing gas.

In Section 2.2 (or Section 2.3), we always explore the effects of time-varying inflow rate (or SFE) while assuming the other to be time invariant. However, we note that in the real universe both inflow rate and SFE could vary with time simultaneously. On galactic scales, the negative correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) in the observations indicates that time-varying inflow rates are dominant. However, this does not mean that the SFE is fully time invariant at galactic scales for all SF galaxies. Actually, as also mentioned in Section 6.2, the SFE may be temporally enhanced in some galaxies, due to some physical processes, such as a merger or interaction with close companions. Mergers and interactions may trigger gravitational instabilities in the cold gas and further enhance the SFE. At the small 100 pc scale, although the variation of SFE is dominant, we can probably ignore the possible variation of inflow rate for different regions within individual galaxies. In addition, we do not consider other feedback processes of star formation in the model, except for the outflow, which is assumed to be simply proportional to the SFR.

The dispersal and mixing of enriched supernova ejecta with the interstellar medium are assumed for simplicity in this paper to be instantaneous. The mixing timescale is expected to be strongly dependent on physical scale (Roy & Kunth 1995; de Avillez & Mac Low 2002). Specifically, on scales of 100 pc or less, the mixing timescale should be a few tens of megayears, of order 10 times shorter than the few-hundred-megayear variation timescales of the SFE that we proposed in this paper. Therefore, the simplifying assumption of instantaneous mixing is unlikely to invalidate the conclusions drawn here.

In the model, we also assume that the yield y of metals is uniform both within galaxies and across the galaxy population. The yield y is closely related to the relative number of Type II supernovae, and therefore to the IMF. In the real universe, the IMF may be different from galaxy to galaxy, or even different in different parts of the same galaxy. Indeed, by using a sensitive index of the IMF, ¹³CO/C¹⁸O, Zhang et al. (2018) found that the IMF in dusty starburst galaxies at redshift ∼2–3 may be more top-heavy with respect to the Chabrier (2003) IMF. The top-heavy IMFs would result in larger y than the bottom-heavy IMF, which increases the complexity of understanding the metal enhancement process by star formation.

In Section 4, we presented the observational results obtained using the N2S2H α and Scal metallicity indicators. The two indicators produce broadly consistent results. As discussed in Section 3.3, these two indicators, proposed by DOP16 and PG16, respectively, offer significant improvements and advantages over the previous indicators, like N2 and O3N2. However, we note that when using O3N2 the derived results differ in part from those presented here. Specifically, a negative correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) across the MaNGA galaxy population can still be seen with O3N2, while the ${\rm{\Delta }}\mathrm{log}$ sSFR−${\rm{\Delta }}\mathrm{log}$(O/H) relation of individual spaxels for MAD galaxies is found to depend quite strongly on stellar mass. For galaxies with stellar mass above ∼10^10.5 M_⊙, the ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) of spaxels show positive correlation, which is similar to the results based on N2S2H α and Scal. But for galaxies with stellar masses below ∼10^10.5 M_⊙, the correlation between ${\rm{\Delta }}\mathrm{log}$ sSFR and ${\rm{\Delta }}\mathrm{log}$(O/H) of spaxels becomes negative, different from the results of N2S2H α and Scal shown in this paper. This may be due to the fact that both the N2S2H α and Scal indicators break the degeneracy between metallicity and ionization parameter. Although we prefer to use the N2S2H α and Scal indicators, we mention this alteration of the results with using O3N2 here for those readers who may prefer that metallicity indicator.

In the present work, we have only compared the observational results on low-redshift galaxies with the model predictions. However, we note that our model prediction is also suitable for high-redshift galaxies. One may expect to push the analysis in the current work to high redshift in the near future, based on near-infrared spectroscopic galaxy surveys with the JWST.