Disentangling the physical parameters of gaseous nebulae and galaxies

TL;DR

This study analyzes how nebular parameters like metallicity, ionization, and density relate to galaxy properties such as stellar mass and star formation rate, using SDSS data and emission line diagnostics.

Contribution

It introduces new scaling relations linking nebular conditions to galaxy properties, derived from composite spectra and multi regression analysis.

Findings

Ionization parameter U correlates strongly with sSFR.

Metallicity Z shows a weaker inverse relation with U.

Derived relations align with theoretical models and observations.

Abstract

We present an analysis to disentangle the connection between physical quantities that characterize the conditions of ionized HII regions -- metallicity ($Z$), ionization parameter ($U$), and electron density ($n_\mathrm{e}$) -- and the global stellar mass ($M_\ast$) and specific star formation rate ($\mathrm{sSFR}=\mathrm{SFR}/M_\ast$) of the host galaxies. We construct composite spectra of galaxies at $0.027 \le z \le 0.25$ from Sloan Digital Sky Survey, separating the sample into bins of $M_\ast$ and sSFR, and estimate the nebular conditions from the emission line flux ratios. Specially, metallicity is estimated from the direct method based on the faint auroral lines [OIII]$\lambda$4363 and [OII]$\lambda\lambda$7320,7330. The metallicity estimates cover a wide range from $12+\log\mathrm{O/H}\sim7.6\textrm{--}8.9$. It is found that these three nebular parameters all are tightly correlated with the location in the $M_\ast$--sSFR plane. With simple physically-motivated ans\"atze, we derive scaling relations between these physical quantities by performing multi regression analysis. In particular, we find that $U$ is primarily controlled by sSFR, as $U \propto \mathrm{sSFR}^{0.43}$, but also depends significantly on both $Z$ and $n_\mathrm{e}$. The derived partial dependence of $U \propto Z^{-0.36}$ is weaker than the apparent correlation ($U\propto Z^{-1.52}$). The remaining negative dependence of $U$ on $n_\mathrm{e}$ is found to be $U \propto n_\mathrm{e}^{-0.29}$. The scaling relations we derived are in agreement with predictions from theoretical models and observations of each aspect of the link between these quantities. Our results provide a useful set of equations to predict the nebular conditions and emission-line fluxes of galaxies in semi-analytic models.

Tables

Table 1. Table 1: Derived physical quantities of the stacks

$M_{\ast }$ bin		sSFR bin		$N_{\mathrm{gal}}$	$\log M_{\ast }\left(M_{\odot }\right)$	$\log \mathrm{sSFR}\left({\mathrm{yr}}^{-1}\right)$	$12+\log \left(O/H\right)$	$\log U$	$\log n_e\left({\mathrm{cm}}^{-3}\right)$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
$7.4$	$7.6$	$-8.50$	$-8.25$	$2$	$7.556\pm 0.160$	$-8.348\pm 0.169$	$8.121\pm 0.057$	$-2.795\pm 0.014$	$1.913\pm 0.164$
$7.4$	$7.6$	$-8.25$	$-8.00$	$3$	$7.533\pm 0.127$	$-8.115\pm 0.132$	$7.735\pm 0.063$	$-2.631\pm 0.015$	$1.073\pm 0.810$
$7.4$	$7.6$	$-8.00$	$-7.75$	$2$	$7.512\pm 0.053$	$-7.983\pm 0.077$	$7.748\pm 0.017$	$-2.264\pm 0.010$	$0.870\pm 0.706$
$7.4$	$7.6$	$-7.75$	$-7.50$	$1$	$7.587\pm 0.082$	$-7.641\pm 0.136$	$7.744\pm 0.091$	$-2.716\pm 0.025$	$2.427\pm 0.371$
$7.6$	$7.8$	$-8.50$	$-8.25$	$1$	$7.798\pm 0.075$	$-8.476\pm 0.120$	$8.329\pm 0.073$	$-2.642\pm 0.025$	$1.306\pm 0.693$
$7.6$	$7.8$	$-8.25$	$-8.00$	$3$	$7.693\pm 0.067$	$-8.101\pm 0.082$	$7.966\pm 0.031$	$-2.245\pm 0.017$	$2.229\pm 0.677$
$7.6$	$7.8$	$-8.00$	$-7.75$	$3$	$7.650\pm 0.051$	$-7.839\pm 0.064$	$7.861\pm 0.022$	$-2.426\pm 0.010$	$1.694\pm 0.662$
$7.6$	$7.8$	$-7.75$	$-7.50$	$1$	$7.601\pm 0.092$	$-7.542\pm 0.151$	$8.043\pm 0.029$	$-2.250\pm 0.014$	$2.108\pm 0.377$
$7.8$	$8.0$	$-9.00$	$-8.75$	$1$	$7.994\pm 0.133$	$-8.754\pm 0.268$	$7.996\pm 0.063$	$-2.894\pm 0.021$	$1.166\pm 0.761$
$7.8$	$8.0$	$-8.75$	$-8.50$	$5$	$7.941\pm 0.032$	$-8.574\pm 0.060$	$8.027\pm 0.033$	$-2.732\pm 0.011$	$1.943\pm 0.175$
$7.8$	$8.0$	$-8.50$	$-8.25$	$4$	$7.973\pm 0.033$	$-8.366\pm 0.071$	$7.942\pm 0.083$	$-2.629\pm 0.029$	$2.180\pm 0.523$
$7.8$	$8.0$	$-8.25$	$-8.00$	$2$	$7.897\pm 0.082$	$-8.200\pm 0.094$	$8.103\pm 0.027$	$-2.593\pm 0.009$	$1.799\pm 0.453$
$7.8$	$8.0$	$-8.00$	$-7.75$	$1$	$7.864\pm 0.140$	$-7.958\pm 0.155$	$7.852\pm 0.022$	$-2.261\pm 0.013$	$2.214\pm 0.397$
$8.0$	$8.2$	$-9.25$	$-9.00$	$3$	$8.181\pm 0.046$	$-9.109\pm 0.159$	$8.151\pm 0.118$	$-2.925\pm 0.026$	$0.718\pm 0.714$
$8.0$	$8.2$	$-9.00$	$-8.75$	$12$	$8.171\pm 0.026$	$-8.822\pm 0.070$	$8.007\pm 0.058$	$-2.819\pm 0.019$	$1.040\pm 0.704$
$8.0$	$8.2$	$-8.75$	$-8.50$	$19$	$8.125\pm 0.030$	$-8.586\pm 0.053$	$8.177\pm 0.031$	$-2.733\pm 0.008$	$1.660\pm 0.276$
$8.0$	$8.2$	$-8.50$	$-8.25$	$16$	$8.130\pm 0.024$	$-8.331\pm 0.039$	$7.982\pm 0.021$	$-2.594\pm 0.008$	$1.825\pm 0.289$
$8.0$	$8.2$	$-8.25$	$-8.00$	$9$	$8.135\pm 0.029$	$-8.134\pm 0.041$	$8.001\pm 0.017$	$-2.525\pm 0.008$	$1.487\pm 0.683$
$8.0$	$8.2$	$-8.00$	$-7.75$	$2$	$8.164\pm 0.089$	$-7.856\pm 0.099$	$8.056\pm 0.021$	$-2.360\pm 0.012$	$1.967\pm 0.382$
$8.0$	$8.2$	$-7.75$	$-7.50$	$2$	$8.049\pm 0.129$	$-7.556\pm 0.133$	$7.617\pm 0.035$	$-2.463\pm 0.011$	$2.003\pm 0.589$
$8.2$	$8.4$	$-9.25$	$-9.00$	$30$	$8.352\pm 0.021$	$-9.073\pm 0.063$	$8.238\pm 0.064$	$-2.931\pm 0.012$	$1.378\pm 0.419$
$8.2$	$8.4$	$-9.00$	$-8.75$	$40$	$8.309\pm 0.015$	$-8.885\pm 0.042$	$8.259\pm 0.031$	$-2.791\pm 0.007$	$1.250\pm 0.472$
Notes.(1–2) lower and upper $\log M_{\ast }$ limits of the stack. (3–4) lower and upper $\log \mathrm{sSFR}$ limits of the stack. (5) number of galaxies in the stack. (6) average stellar mass of the stack. (7) average sSFR of the stack. (8–10) estimated metallicity, ionization parameter, and electron density of the stack. The table is published in its entirety in the machine-readable format.

Table 2. Table 2: Coefficients of the linear surface fitting between the parameters a

	Intercept	$\log M_{\ast }\left(M_{\odot }\right)-10$	$\log n_e\left({\mathrm{cm}}^{-3}\right)$	$4+\log \left(O/H\right)$	$\log \mathrm{sSFR}\left({\mathrm{yr}}^{-1}\right)+9$	${\chi }^2/\mathrm{dof}$b	${\sigma }_{\mathrm{int}}$c
$\log \mathrm{sSFR}+9$	$-{3.661}_{-0.125}^{+0.124}$	$-{0.627}_{-0.029}^{+0.028}$	${1.753}_{-0.064}^{+0.063}$	-	-	0.99	-
$4+\log \left(O/H\right)$	${1.238}_{-0.046}^{+0.049}$	${0.228}_{-0.011}^{+0.011}$	$-{0.401}_{-0.024}^{+0.024}$	-	-	1.61	0.0633
$\log U$	$-{2.316}_{-0.039}^{+0.043}$	-	$-{0.292}_{-0.022}^{+0.020}$	$-{0.360}_{-0.037}^{+0.037}$	${0.428}_{-0.014}^{+0.014}$	1.89	0.0585
$\log n_e$	${2.066}_{-0.012}^{+0.012}$	${0.310}_{-0.019}^{+0.021}$	-	-	${0.492}_{-0.021}^{+0.020}$	0.85	-
$4+\log \left(O/H\right)$	${0.427}_{-0.007}^{+0.007}$	${0.136}_{-0.006}^{+0.007}$	-	-	$-{0.208}_{-0.008}^{+0.008}$	5.71	0.0998
$\log U$	$-{3.073}_{-0.002}^{+0.002}$	$-{0.137}_{-0.002}^{+0.002}$	-	-	${0.372}_{-0.002}^{+0.002}$	26.07	0.0874
a We note that $M_{\ast }$, sSFR, and $n_e$ are given in units of $M_{\odot }$, ${\mathrm{yr}}^{-1}$, and ${\mathrm{cm}}^{-3}$, respectively.
b Reduced chi-square statistic.
c Intrinsic scatter with respect to the relevant objective variable.

Equations

$$\frac{\mathrm{O}}{\mathrm{H}}=\frac{{\mathrm{O}}^{+}}{{\mathrm{H}}^{+}}+\frac{{\mathrm{O}}^{++}}{{\mathrm{H}}^{+}}.$$

$$U=q_{\mathrm{ion}}/c=n_{\gamma }/n_{\mathrm{H}}$$

$$12+\log (O/H)$$

$$\log U$$

$$\log n_{\mathrm{H}}/(\mathrm{c}{\mathrm{m}}^{-3})$$

$$\log \mathrm{sSFR}+9$$

$$4+\log O/H$$

$$\log U$$

$$+r_2(4+\log O/H)+r_3(\log \mathrm{sSFR}+9).$$

$$\mathcal{L}=\prod _{i=1}^N\frac{1}{\sqrt{2\pi {\sigma }_{\mathrm{err},\mathrm{i}}^2}}\exp \left[-\frac{(y_i^{\mathrm{obs}}-y_i^{\mathrm{mod}}{)}^2}{2{\sigma }_{\mathrm{err},i}^2}\right],$$

$${\sigma }_{\mathrm{err},i}^2=\delta y_i^2+\sum _{k=1}^m{\left(p_k\delta x_{k,i}\right)}^2.$$

$${\chi }^2/\mathrm{dof}=\sum _{\mathrm{i}=1}^N\frac{(y_i^{\mathrm{obs}}-y_i^{\mathrm{model}}{)}^2}{{\sigma }_{\mathrm{err},i}^2}/(N-m-1).$$

$$\sum _{\mathrm{i}=1}^N\frac{(y_i^{\mathrm{obs}}-y_i^{\mathrm{model}}{)}^2}{{\sigma }_{\mathrm{err},i}^2+{\sigma }_{\mathrm{int}}^2}/(N-m-1)=1.$$

$$\log \mathrm{sSFR}$$

$$12+\log O/H$$

$$\log U$$

$$-0.292\log n_{\mathrm{e}}+0.428\left(\log \mathrm{sSFR}+9\right).$$

$$\mathrm{sSFR}=\mathrm{SFE}\cdot M_{\mathrm{gas}}/M_{\ast }=\mathrm{SFE}\cdot {\mu }_{\mathrm{gas}}.$$

$$\mathrm{SFE}\propto {\mathrm{sSFR}}^{0.29}{M}_{\ast }^{0.18}.$$

$${\mu }_{\mathrm{gas}}\propto M_{\ast }^{-0.63}{n}_{\mathrm{e}}^{1.3}\text{or}{M}_{\mathrm{gas}}\propto M_{\ast }^{0.37}{n}_{\mathrm{e}}^{1.3}$$

$$⟺n_{\mathrm{e}}\propto {\mu }_{\mathrm{gas}}^{0.80}{M}_{\ast }^{0.50}\propto M_{\mathrm{gas}}^{0.80}{M}_{\ast }^{-0.30}.$$

$$\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.40}{M}_{\ast }^{-0.15}.$$

$$Z\propto M_{\ast }^{0.23}{n}_{\mathrm{e}}^{-0.40}\propto {\mu }_{\mathrm{gas}}^{-0.32}{M}_{\ast }^{0.027},$$

$$Z\propto M_{\mathrm{gas}}^{-0.32}{M}_{\ast }^{0.35}.$$

$$12+\log (O/H)=8.427+0.136\left(\log M_{\ast }-10\right)-0.208\left(\log \mathrm{sSFR}+9\right).$$

$${\mu }_{\alpha }=\log M_{\ast }-\alpha \log \mathrm{SFR}.$$

$$12+\log (O/H)=(5.19\pm 0.06)+(0.344\pm 0.006){\mu }_{0.60}.$$

$$f_{\mathrm{abs}}(\mathrm{H}\alpha )=\frac{1}{2}\left[\mathrm{erf}\left(-0.373\left(x+2.331\right)\right)+1\right],$$

$$f_{\mathrm{abs}}(\mathrm{H}\beta )=\frac{1}{2}\left[\mathrm{erf}\left(-0.720\left(x+0.107\right)\right)+1\right].$$

Full text

††pubyear: 2015††pagerange: Disentangling the physical parameters of gaseous nebulae and galaxies–B

1 Introduction

Ionized nebulae in star-forming galaxies show a wide range of physical conditions. The properties of H ii regions are characterized by a set of physical quantities, including density (or pressure), gas-phase metallicity (hereafter metallicity or $Z$), and the ionization parameter $U$. These parameters are connected through various fundamental physical processes taking place in the gas, such as photoionization, collisional excitation, and radiative cooling. Interestingly, the properties of the interstellar medium (ISM) appear also strongly correlated with the global properties of the host galaxies, such as stellar mass ($M_{\ast}$) and star formation rate (SFR). This implies that the history and the on-going activities of the galaxies affect the conditions and the processes occurring in the resident nebulae. Disentangling the connections between these parameters is thus essentially important to understand the astrophysical processes govern the evolution and behavior of baryons in the galaxies.

The correlation between $M_{\ast}$ and $Z$ (the “mass–metallicity” (MZ) relation) has been well established at low redshifts ($z\sim 0.1$) based on a large sample from the Sloan Digital Sky Survey (SDSS) (e.g., Tremonti et al., 2004), and also at higher redshifts ($z\gtrsim 1$) as well (e.g., Erb et al., 2006; Zahid et al., 2011; Yabe et al., 2012; Zahid et al., 2013; Zahid et al., 2014; Sanders et al., 2015; Kashino et al., 2017, 2018). Ellison et al. (2008) pointed out the presence of an anticorrelation between $Z$ and SFR at fixed $M_{\ast}$, which has then further explored by many authors (e.g., Lara-López et al., 2010; Yates et al., 2012; Lara-López et al., 2013; Mannucci et al., 2010; Andrews & Martini, 2013). In particular, Mannucci et al. (2010) proposed the so-called fundamental metallicity relation (FMR) in the $M_{\ast}$-SFR-$Z$ space. In these studies, the scatter of the MZ relation is explained by the result of variations in the inflowing rate of circumgalactic medium into the system; inflowing gas drives formation of new stars while diluting the gas-phase metallicity. In this context, it has also been argued that the relation between $M_{\ast}$, $Z$ and the gas content is more fundamental, rather than SFR (e.g., Bothwell et al., 2013; Bothwell et al., 2016a, b). However, the shape of the FMR and its redshift evolution are still under debate (e.g., Yates et al., 2012; Andrews & Martini, 2013; Kashino et al., 2016; Telford et al., 2016).

A difficulty in measuring metallicity is that in most cases we cannot use the “direct” method and thus must rely on the so-called “strong-line” method. The direct method utilizes the flux ratio between auroral and strong lines, such as [O iii]$\lambda$4363/[O iii]$\lambda$5007, to measure the electron temperature $T_{\mathrm{e}}$ of the ionized gas. Once the temperature is determined, the absolute oxygen abundance can be estimated from the $T_{\mathrm{e}}$-dependent ratios of strong lines, such as [O ii]$\lambda$3727 and [O iii]$\lambda$5007, relative to the H$\beta$ line. On the other hand, the strong-line method is a technique that uses a relation between strong line ratios (e.g., $R_{\mathrm{23}}=(\textrm{[O\,{\sc ii}]}\lambda 3727+\textrm{[O\,{\sc iii}]}\lambda 4959,5007)/\mathrm{H\beta}$) and metallicity, calibrated theoretically or observationally (e.g., Kewley & Dopita, 2002; Pettini & Pagel, 2004; Nagao et al., 2006; Maiolino et al., 2008; Kewley & Ellison, 2008; Marino et al., 2013). Since the detection of the auroral lines is usually challenging due to their faintness, especially at high metallicity, the strong-line method has been routinely used. However, the strong-line ratios are functions not only of metallicity, but also of other parameters such as the ionization parameter. Such systematic uncertainties could result in a misleading comparison between different galaxy populations, e.g., low and high redshifts galaxies (e.g., Kewley et al., 2015; Telford et al., 2016).

The electron density $n_{\mathrm{e}}$ of the ionized gas is routinely estimated from the intensity ratio of the [S ii]$\lambda$6717,6731 or [O ii]$\lambda$3726,3729 doublet lines. The typical electron densities of H ii regions in low redshift galaxies have been measured to be $n_{\mathrm{e}}\sim 50\textrm{--}120~{}\mathrm{cm^{-3}}$ (Brinchmann et al., 2008). The densities in high redshift ($z\gtrsim 1\textrm{--}2$) galaxies, however, have found to often be elevated relative to low redshifts, ranging from $10^{2}$ to a few $\times 10^{3}~{}\mathrm{cm^{-3}}$ (e.g., Shirazi et al., 2014; Masters et al., 2014; Shimakawa et al., 2015; Sanders et al., 2016; Kashino et al., 2017; Kaasinen et al., 2017). In particular, Shimakawa et al. (2015) found that $n_{\mathrm{e}}$ is correlated with sSFR and/or SFR surface density. Such correlation could naturally explain the observed evolution in the typical $n_{\mathrm{e}}$ with redshift.

Higher ionization parameters are often measured in high redshift and/or high sSFR galaxies (e.g., Hainline et al., 2009; Nakajima & Ouchi, 2014; Shapley et al., 2015; Strom et al., 2017). These studies routinely use the [O iii]$\lambda$5007/[O ii]$\lambda$3727 ratio as a proxy of $U$. Its effectiveness has been established theoretically (e.g., Dopita et al., 2000; Kewley & Dopita, 2002), while it is known that the line ratio is also sensitive to metallicity at fixed $U$. Pioneering studies (Evans & Dopita, 1985; Dopita & Evans, 1986) claimed the presence of a strong inverse correlation between $U$ and $Z$ from comparisons between their theoretical photoionization model and observations. Dopita et al. (2006b) have theoretically shown that stellar ionizing photon flux illuminating the gas reduces at higher metallicity because of increasing absorption due to a higher opacity in the stellar winds, and due to more efficient conversion from radiation energy into kinetic energy of the winds. Recent observations of high redshift galaxies (hence higher sSFR and lower $Z$) have also found anticorrelations between $U$ and $Z$ via the strong-line method (Nakajima & Ouchi, 2014; Onodera et al., 2016). In contrast, however, Kaasinen et al. (2018) recently concluded that the ionization parameter $U$ is directly linked to sSFR, while little or probably not with metallicity, from comparison between samples of local and $z\sim 1.5$ star-forming galaxies. It remains under debate whether $U$ and $Z$ are intrinsically correlated or not, and if it is the case, what physical mechanisms work behind the relationship.

The physical quantities that characterize the conditions of galaxies and the ionized nebulae are connected with each other in complicated ways. As mentioned above, each aspect of this big puzzle has been studied by many authors. Now we need to disentangle their connection from a viewpoint beyond focusing on a correlation between two or three quantities simultaneously. In this paper, we will focus on the integrated global properties of galaxies (i.e., $M_{\ast}$ and SFR) and the basic parameters that define the conditions of the ionized gas locally; metallicity $Z$ (O/H), ionization parameter $U$, and electron density $n_{\mathrm{e}}$. To achieve the goal, accurate and precise measurements of oxygen abundance is essentially important, and thus we use the direct method that requires a detection of faint auroral lines of oxygen ions such as [O iii]$\lambda$4363, and [O ii]$\lambda$7320, 7330. To overcome the faintness of these lines, we employ techniques presented in Andrews & Martini (2013) and Curti et al. (2017), who stacked thousands of galaxy spectra from SDSS. This enables us to use the direct method even for low SFR and metal-rich galaxies.

This paper is organized as follows. We describe our sample of galaxies and the methods of spectral stacking and emission-line measurement in Section 2. Derivation of the physical properties of ionized nebulae (i.e., $Z$, $U$ and $n_{\mathrm{e}}$) is described in Section 3. In Section 4, we describe the set of simple, physically-motivated ansätze for the relations between these quantites, and analyses of the correlations of the stacked measurements. We present results in Section 5, and discussions in Section 6. Summary and the conclusions are given in Section 7. Throughout the paper, the stellar masses and SFRs are computed by using a Kroupa (2001) universal initial mass function (IMF). The metallicity $Z$ denotes the gas-phase oxygen abundance, which is defined as the abundance ratio of oxygen atoms to $10^{12}$ hydrogen atoms, as $12+\log(\mathrm{O/H})$. We adopt the solar metallicity to be $12+\log(\mathrm{O/H})_{\odot}=8.69$, or $Z_{\odot}=0.014$ for metallic mass fraction (Asplund et al., 2009).

2 Data and method

2.1 Sample selection

The parent sample of galaxies used for this study was constructed following the method described in Andrews & Martini (2013). We here present a brief overview of the sample selection. We utilize the public MPA-JHU catalog, which contains $\sim 930,000$ galaxies from the SDSS Data Release 7 (Abazajian et al., 2009). The MPA-JHU catalog contains for each galaxy the estimates of stellar mass and total SFR from the aperture-corrected H$\alpha$ luminosity, with a Kroupa IMF (Brinchmann et al., 2004; Salim et al., 2007).

The galaxies are limited to have a spectroscopic redshift between $0.027\leq z\leq 0.25$ with an uncertainty of $\sigma_{z}<0.001$. For this redshift range, both the most blue [O ii]$\lambda 3727$ line and the most red [O ii]$\lambda\lambda$7320,7330 lines fall within the observed wavelength range of 3800–9200 Å. We further imposed on the sample the detection of H$\beta$, H$\alpha$, and [N ii]$\lambda$6584 all at $S/N\geq 5$. Possible AGNs were removed by using the Baldwin-Phillips-Terlevich (BPT; Baldwin et al. 1981; see also Veilleux & Osterbrock 1987) diagram with an empirical classification curve derived by Kauffmann et al. (2003). Objects with possible poor constraints of photometry were removed, based on the warning flags (DEBLEND_NOPEAK or DEBLENDED_AT_EDGE), or visual inspection. The final sample consists of 194,005 galaxies and the median redshift is $z=0.081$. Figure 1 show the distribution of our selected galaxies in the $M_{\ast}$–sSFR plane by contours. The ridge line of the contours is well matched to the so-called main sequence of star-forming galaxies derived by Elbaz et al. (2007).

2.2 Spectral stacking

To determine the metallicity via the direct method, it is required to measure the fluxes of the faint auroral emission lines of oxygen ions, [O iii]$\lambda$4363 and [O ii]$\lambda\lambda$7320, 7330. However, the intensity of these lines rapidly declines below the noise level of the single spectrum as the metallicity increases. To probe a wide range of the parameter space, we improved the S/N ratios by creating composite spectra of galaxies.

Our motivation is to disentangle the connections between physical parameters characterizing local conditions in ionized nebulae and global properties of host galaxies. The stellar mass and SFR are both extensive variables, and hence are naturally correlated with each other. The interpretation would be more clarified by using $M_{\ast}$ as a single indicator of the size of the system, while taking sSFR, an intensive variable, rather than SFR. We therefore separated the sample galaxies into the grid in the $M_{\ast}$–sSFR plane with binsizes of 0.20 dex in $M_{\ast}$ and 0.25 dex in sSFR. Figure 1 shows the numbers of galaxies in each box that represents a single stack.

We utilized the public spectral data from the SDSS DR7, which have been processed with the SDSS Spec2D pipeline (Stoughton et al., 2002). Before creating composite spectra, the individual spectra were 1) corrected for the Milky Way reddening with the extinction values from Schlegel et al. (1998) and adopting a Cardelli et al. (1989) extinction law, 2) shifted to the rest frame based on the spectroscopic redshift, 3) resampled with a universal wavelength grid with a spacing of $\Delta\log_{10}\lambda=5\times 10^{-5}$, and 4) normalized by the H$\beta$ flux. The individual spectra have then been co-added by taking the mean flux between the 25th and the 75th percentiles at each wavelength grid, similarly to Curti et al. (2017).

2.3 Stellar continuum subtraction

To measure the fluxes of the nebular emission lines, it is important to fit and subtract the stellar component from the composite spectra. In particular, the Balmer lines and [O iii]$\lambda$4363 in the vicinity of H$\gamma$, are highly impacted by the stellar atmospheric absorption. For the stellar continuum subtraction, we generated a synthetic spectrum for each composite spectrum. We performed the fit using the IDL version of the pPXF package (Cappellari, 2017) with the built-in library of single stellar population (SSP) spectral templates. The library includes 150 model spectra based on the MILES library of stellar templates (Sánchez-Blázquez et al., 2006), computed for 25 different ages and 6 different metallicities using the code presented in Vazdekis et al. (2010). The templates cover a wavelength range of 3540–7410 Å with a resolution of 2.51Å (full width at half maximum; FWHM).

We performed the stellar template fitting first using the entire spectral range between 3650 Å and 7360 Å to determine the best-fit synthesized spectrum, generated as an arbitrary superposition of the 150 simple stellar population spectra in the library. We then fit the single best-fit template in each subrange of a few hundred angstroms around the emission lines of interest, allowing slight modifications on the systemic offset, velocity dispersion, and additive and multiplicable factors to improve the quality of the fits. In each spectral subrange, narrow wavelength regions where the emission lines present were masked out. In Figure 2, we show an example of spectral stacking and stellar continuum subtraction. The panels correspond to the spectral ranges for the key emission lines, [O ii]$\lambda\lambda$3726,3729, H$\gamma$, [O iii]$\lambda$4363, H$\beta$, [O iii]$\lambda\lambda$4959,5007, H$\alpha$, [N ii]$\lambda\lambda$6548,6584, [S ii]$\lambda\lambda$6717,6731, and [O ii]$\lambda\lambda$7320,7330, respectively.

Figure 2 shows that the absorption features are prominent at wavelengths of the Balmer series lines. These Balmer absorption lines arise from the absorption of photons by hydrogen excited at the $n=2$ level in the photospheres of stars. This effectively reduces the observed fluxes of the Balmer lines if one simply excludes the continuum based on a linear or low-order polynomial fitting around the lines. Meanwhile, it is not always possible to perform the stellar continuum subtraction based on the stellar population synthesis fitting, as done in this paper. Therefore, it would be useful to derive a simple empirical relation between the effect of Balmer absorption and the galaxy properties. In Appendix A, we present the amount of the Balmer absorption as a function of $M_{\ast}$ and sSFR, and provide empirical relations for absorption correction.

2.4 Line flux measurement

We performed emission-line fitting using the MPFIT package for IDL (Markwardt, 2009) to measure the fluxes of the emission lines in the stellar-subtracted composite spectra. Each emission line was fit with a single gaussian profile. Note that all the doublet lines ([O ii]$\lambda\lambda$3726,3729, [O iii]$\lambda\lambda$4959,5007, [N ii]$\lambda\lambda 6548,6584$, [S ii]$\lambda\lambda 6717,6731$, and [O ii]$\lambda\lambda 7320,7330$) were fit with independent amplitudes. We disregarded lines that were poorly fit with low S/N ($<5$), which is computed for each line as a ratio of measured flux and the uncertainty based on the errors on the fitting parameters returned by the MPFIT procedure. The measured fluxes of all emission lines were then corrected for dust extinction by adopting a Cardelli et al. (1989) extinction law. The absolute level of extinction was estimated based on the observed H$\alpha$/H$\beta$ ratio by assuming the intrinsic values of H$\alpha$/H$\beta$ given in Osterbrock & Ferland (2006, Table 4-4) for Case B recombination with $n_{\mathrm{e}}=10^{2}~{}\mathrm{cm^{-3}}$. The intrinsic H$\alpha$/H$\beta$ ratio is very insensitive to the electron density if $n_{\mathrm{e}}<10^{4}~{}\mathrm{cm^{-3}}$, while being slightly dependent of the electron temperature. For the range of $7000\lesssim T_{\mathrm{e}}(\mathrm{O^{+}})\lesssim 20000~{}\mathrm{K}$ of our data, the H$\alpha$/H$\beta$ ratio could vary within $\pm 5\%$, increasing as $T_{\mathrm{e}}$ decreases. Therefore we iteratively estimated the extinction as follows. We first estimate the extinction by assuming the intrinsic H$\alpha$/H$\beta$ ratio of 2.86 for $T_{\mathrm{e}}=10^{4}~{}\mathrm{K}$ and derive the electron temperature from the extinction-corrected fluxes. The estimated $T_{\mathrm{e}}$ is then used to modify the intrinsic H$\alpha$/H$\beta$ ratio, and to re-compute the level of extinction. We then repeat the extinction correction of all the line fluxes and derivation of the physical quantities. The change in the resulting $\log\mathrm{O/H}$ after the iteration is correlated with $\log\mathrm{O/H}$, and is negligible at $\log\mathrm{O/H}\lesssim 8.2$ while gradually being larger at higher $Z$. The final estimates after the iteration are reduced by $\approx 0.08~{}\mathrm{dex}$ at $12+\log\mathrm{O/H}=8.8$ relative to the initial estimates with the intrinsic H$\alpha$/H$\beta$ of 2.86 assumed.

As pointed out by Andrews & Martini (2013), an emission feature was found at $\approx$4360 Å, which is blended with the [O iii]$\lambda$4363 auroral line. Curti et al. (2017) suggested that this feature is reasonably associated with emission lines from Fe ii ions. Following their treatment, we simultaneously fit the three line components of the possible [Fe ii] lines that are centered on 4357.7, 4358.8, and 4360.6 Å, in addition to H$\gamma$ and [O iii]$\lambda$4363. As shown by Curti et al. (2017), we found that the contribution from the possible [Fe ii] lines increases as metallicity increases. Following Curti et al. (2017), we disregarded possible detection of [O iii]$\lambda$4363 when the flux measured for the $\lambda$4360 feature is $>0.5$ times the flux measured for [O iii]$\lambda$4363.

3 Physical parameters

In this work, we investigate the connection among the five physical quantities: two galaxy global properties, $M_{\ast}$ and sSFR, and three nebular parameters characterizing the conditions of the H ii regions, metallicity $Z$ ($12+\log(\mathrm{O/H})$), ionization parameter $U$, and electron density $n_{\mathrm{e}}$. These nebular properties were derived based on the stacked measurements of the emission-line fluxes as described in the following subsections. For this purpose, we used an package of IDL routines called impro111The code has been developed by John Moustakas et al. and is available here: https://github.com/moustakas/impro. This software offers functions similar to the commonly used IRAF/TEMDEN package. To evaluate the nebular parameters and their uncertainties we generated 200 random realizations for each set of the line fluxes measured from each stack while incorporating the flux errors with Gaussian distributions. We then took the median and the standard deviation for each parameter estimate from the resultant realizations.

3.1 Electron temperature and density

The electron temperatures $T_{\mathrm{e}}(\mathrm{O}^{+})$ and $T_{\mathrm{e}}(\mathrm{O}^{++})$ were determined from the line ratios [O ii]$\lambda\lambda$3726,3729/[O ii]$\lambda\lambda$7320,7330 and [O iii]$\lambda\lambda$4959,5007/[O iii]$\lambda$4363, respectively. If we detected both auroral lines, we independently estimated $T_{\mathrm{e}}(\mathrm{O}^{+})$ and $T_{\mathrm{e}}(\mathrm{O}^{++})$ using the impro routines. However, we failed to detect [O iii]$\lambda$4363 for stacks with an estimated metallicity of $12+\log\mathrm{O/H}\gtrsim 8.3$ due to the flux being faint and the contamination of the $\lambda 4360$ feature being significant. For these stacks, we estimated the temperature $T_{\mathrm{e}}(\mathrm{O}^{++})$ by utilizing the so-called flux-flux (ff) relation derived by Pilyugin et al. (2006), following the method in Curti et al. (2017). This empirically-calibrated relation predicts the [O iii]$\lambda$4363 flux relative to H$\beta$ from the strong-line ratios, [O iii]$\lambda$4959,5007/H$\beta$ and [O ii]$\lambda$3727/H$\beta$ (see Equation (14) of Pilyugin et al. 2006). Similarly, for stacks where we failed to detect the [O ii]$\lambda\lambda$7320,7330 lines, we utilized the empirical relation using the strong-line ratios given by Curti et al. (2017, see Equation 3). In this work, we use the stacks where at least we successfully detect either the [O iii] or [O ii] auroral lines, and disregarded stacks having neither detections of these. We emphasize that, at the metallicity range where the [O iii]$\lambda$4363 line is undetected, the total oxygen abundance is dominated by the singly-ionized oxygen $\mathrm{O}^{+}$. This is also the case for the opposite situation, i.e., the oxygen abundance of the stacks without the [O ii]$\lambda$7320,7330 detection is dominated by $\mathrm{O}^{++}$. Therefore, the effects of the possible uncertainties from the use of these empirical relations are negligible in the total oxygen abundance estimates.

The flux ratio of the auroral and strong lines depends not only on the electron temperature, but also weakly on the electron density. The electron density can be estimated from the [S ii]$\lambda$6717/[S ii]$\lambda$6731 ratio, with a small dependence on the temperature $T_{\mathrm{e}}(\mathrm{S}^{+})$. We therefore iteratively estimated both $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ assuming $T_{\mathrm{e}}(\mathrm{S}^{+})=T_{\mathrm{e}}(\mathrm{O}^{+})$. The [S ii] doublet ratio saturates at $n_{\mathrm{e}}<10~{}\mathrm{cm^{-3}}$. As a result, the errors is getting larger for lower $n_{\mathrm{e}}$, or cannot be constrained for some cases. In our analysis, we set the lower limit of the $n_{\mathrm{e}}$ range to be $n_{\mathrm{e}}=5~{}\mathrm{cm^{-3}}$. We substituted $\log n_{\mathrm{e}}/(\mathrm{cm^{-3}})=0.7\pm 0.7$ for the value of three stacks for which the observed [S ii] ratio is above the saturation limit. The derived $n_{\mathrm{e}}$ of the stacks ranges from $\sim 17$ to $\sim 300~{}\mathrm{cm^{-3}}$ (the central 90 percentiles), with the median $n_{\mathrm{e}}=76~{}\mathrm{cm^{-3}}$. Note that this range and the median value do not represent the electron densities of the individual SDSS galaxies since the number of galaxies in the stacks are hugely different from one to another bin. Lastly, we note that the spectral resolution of the SDSS spectra are not sufficient to obtain reasonable estimates of $n_{\mathrm{e}}$ from the [O ii]$\lambda\lambda$3726,3729 doublet ratio.

3.2 Direct abundance determination

The ionic oxygen abundance O+/H+ and O*++*/H+ were determined from [O ii]$\lambda\lambda$3726,3727/H$\beta$, and [O iii]$\lambda$5007/H$\beta$ as well as from the temperature estimates, $T_{\mathrm{e}}(\mathrm{O}^{+})$ and $T_{\mathrm{e}}(\mathrm{O}^{++})$, respectively. The procedure in the impro package includes a slight dependence on $n_{\mathrm{e}}$ as well. As described above, if we have constraint only either $T_{\mathrm{e}}(\mathrm{O}^{+})$ or $T_{\mathrm{e}}(\mathrm{O}^{++})$, we estimated the other using the empirical relations based on the strong lines. We remind that, for such cases, the contribution of the ions whose auroral line(s) is not detected is subdominant. We then assumed that the total oxygen abundance is the sum of the ionic abundances of these two ions,

[TABLE]

We neglected the minimal contribution of neutral atoms since our attention is focused on the conditions of the H ii regions, where the neutral fraction is typically small ($\sim 10^{-4}$). The contribution of $\mathrm{O^{+++}}$ and higher ionization ions was also neglected, whose ionization potential is 54.9 eV or higher, thus there is little chance to be ionized by stellar radiation.

3.3 Ionization parameter

The ionization parameter is a key quantity in this work. It is locally defined as the ratio of ionizing photon flux $S_{\mathrm{ion}}$ and the hydrogen density $n_{\mathrm{H}}$: $q_{\mathrm{ion}}=S_{\mathrm{ion}}/n_{\mathrm{H}}$. The ionization parameter can be made dimensionless by dividing by the speed of light as follows:

[TABLE]

where $n_{\gamma}$ is the volume density of the ionizing photons. Throughout the paper, we use the dimensionless $U$, instead of $q_{\mathrm{ion}}$. We note that the electron density $n_{\mathrm{e}}$ is almost equivalent to $n_{\mathrm{H}}$ in a usual H ii region.

We estimated $U$ from a sensitive indicator, the flux ratio of [O iii]$\lambda$5007/[O ii]$\lambda$3726,3729 (Dopita et al., 2000; Kewley & Dopita, 2002). Although the flux ratio is also sensitive to metallicity, we can resolve this degeneracy using the metallicity estimate from the direct method (Section 3.2). The [O iii]/[O ii] ratio also slightly varies with $n_{\mathrm{e}}$ at fixed $U$ and $Z$, as well. We consider this effect by using the $n_{\mathrm{e}}$ estimate from the [S ii] doublet ratio.

We determined the ionization parameters for each stack by interpolating the flux ratios computed for grids of $Z$, $U$, and $n_{\mathrm{H}}$, assuming $n_{\mathrm{H}}=n_{\mathrm{e}}$. We used the CLOUDY (C17.00) photoionization code to compute the emission-line fluxes for given nebular conditions. For the gas phase, we adopted the solar chemical composition, though the chemical composition does little affect the estimates of $U$ since it is estimated from the flux ratio of the same species at different ionization levels (i.e., O+ and O*++*). Photoionization calculations were executed in the pseudo plane-parallel geometry (i.e., inner radius $\gg$ thickness of the ionized shell) for the following values of $Z$, $U$, and $n_{\mathrm{H}}$.

[TABLE]

Ionizing spectra used in our CLOUDY calculations were produced by the Starburst99 population synthesis code (Leitherer et al., 1999), adopting an constant SFR history with a Kroupa (2001) IMF ($0.1\textrm{--}100~{}M_{\odot}$) and the Geneva high mass-loss evolutionary tracks with five values of stellar metallicity: $Z_{\ast}=0.001$, 0.004, 0.008, 0.020, and 0.040, corresponding to $12+\log(\mathrm{O/H})=7.54\textrm{--}9.15$. We adopted a set of spectra at $10~{}\mathrm{Myr}$, and, for a given gas-phase $Z$, the input ionizing spectrum was produced by interpolating the spectral set at $Z_{\ast}=Z$ to match the stellar and gas-phase metallicities. The ionization parameter $U$ is then derived from a $U$–[O iii]/[O ii] relation obtained by interpolating the three-dimensional grid data at the estimated $n_{\mathrm{H}}=n_{\mathrm{e}}$ and $\log\mathrm{O/H}$. The estimates of $U$ range from $\log U=-3.5$ to $-2.4$.

3.4 Derived quantities

We successfully estimated the metallicity for 138 stacks, of which we detected both [O iii]$\lambda$4363 and [O ii]$\lambda\lambda$7320,7330 for 75 stacks, only [O ii]$\lambda\lambda$7320,7330 for 61 stacks, and only [O iii]$\lambda$4363 for 2 stacks. Figure 1 shows the derived metallicity by colors in each stacking bin in the $M_{\ast}$–sSFR plane. Here a caveat is that the stacks with successful metallicity estimation are biased towards lower $M_{\ast}$ and higher SFR, resulting in being biassed to the population above the main sequence (dotted line; Elbaz et al. 2007). This is because the intensity of auroral lines rapidly declines with increasing $Z$ as well as with decreasing SFR. As a result, the detection of the auroral lines has often failed at high $M_{\ast}(\gtrsim 10^{10.6}~{}M_{\odot})$, or below the main sequence even if thousands of spectra are stacked. In contrast, the detection succeeded at lower $M_{\ast}$ and high SFRs even with a few (or only a single) galaxies in the stacks. However, we would emphasize that the successful stacks cover the whole range of the metallicity of interest ($7.6\lesssim 12+\log(\mathrm{O/H})\lesssim 8.9$), probed by most studies of the MZ relation at low and higher redshifts (references in Section 1), including both regimes of extremely metal-poor galaxy populations (e.g., Sánchez Almeida et al., 2016) and the super-solar.

In Figure 3, the ionization parameter $U$ and electron density $n_{\mathrm{e}}$ are presented in the same manner as metallicity. All these panels in Figures 1 and 3 clearly show that the nebular conditions depend strongly on $M_{\ast}$ and sSFR. It is shown that $Z$ and $U$ are more sensitive to sSFR than to $M_{\ast}$. Meanwhile, the electron density $n_{\mathrm{e}}$ increases both with $M_{\ast}$ at fixed sSFR as well as with sSFR at fixed $M_{\ast}$. Table 1 summarize the derived physical properties (full version including emission line flux measurements available online). The relations between the observed strong-line fluxes and derived physical quantities are presented in Appendix B.

4 Correlation analysis

In this section, we attempt to disentangle the connection between the parameters of the ionized H ii regions and the host galaxies. Figure 4 shows correlations for all possible combinations among the five parameters of interest, $\log M_{\ast}/M_{\odot}$, $\log\mathrm{sSFR}/(\mathrm{yr^{-1}})$, $12+\log(\mathrm{O/H})$, $\log U$, and $\log n_{\mathrm{e}}/(\mathrm{cm^{-3}})$. The stellar mass and sSFR of each stack are defined as the median values of individual galaxies within the bins. The data points are color-coded by sSFR, and the symbol size denotes the number of galaxies in each stacking bin, as indicated outside the frames.

4.1 Ansätze

We here propose a set of ansätze to figure out the connections between the parameters. First we assume that the partial correlations between any two of the parameters are expressed by power-law functions, i.e., we assume linearity for their logarithms. Although this is undoubtedly an oversimplification, it is reasonably supported by the observed relationships shown in Figure 4, and indeed a more elaborate parametrization is not justified by our limited data. Still, the power-law slopes that we obtain from the data will be helpful to qualitatively understand the physics controlling the conditions of the nebulae and galaxy evolution.

To shape a course to analyze the observed correlations, we employ further assumptions for the dependencies of the parameters. In Figure 5, we schematically show our assumption about the dependencies of the parameters. The black (white) arrows indicate the direction of possible causation which is expected to lead to a positive (negative) correlation between the two parameters at their starting and ending points. We here consider for convenience an additional parameter, the gas mass ratio $\mu_{\mathrm{gas}}=M_{\mathrm{gas}}/M_{\ast}$, which is assumed to be uniquely determined once the key parameters of interest are given.

One may naturally argue that the SFR of a galaxy increase with the amount of the fuel, i.e., the gas content in the system. Star formation consumes the gas reservoir while accumulating the stellar mass of the galaxy. Thus, the gas mass ratio is naturally expected to be anticorrelated with $M_{\ast}$. These causations are illustrated, with labels “Fueling” and “Consumption” respectively, in Figure 5.

Mathematically, the SFR integrated over a galaxy is related to its total molecular gas mass as $\mathrm{SFR}=M_{\mathrm{gas}}/\tau_{\mathrm{dep}}$ where $\tau_{\mathrm{dep}}$ is the gas depletion timescale, and the star formation efficiency is defined as its inverse $\mathrm{SFE}=1/\tau_{\mathrm{dep}}=\mathrm{SFR}/M_{\mathrm{gas}}$. On the other hand, a simple hypothesis relates the local SFR volume density $\dot{\rho}_{\ast}$ with the gas volume density $\rho_{\mathrm{gas}}$ and free-fall timescale $\tau_{\mathrm{ff}}$ of an overdense region as $\dot{\rho}_{\ast}\propto\rho_{\mathrm{gas}}/\tau_{\mathrm{ff}}\propto\rho_{\mathrm{gas}}^{1.5}$ since $\tau_{\mathrm{ff}}\propto\rho_{\mathrm{gas}}^{1/2}$ (e.g., Schmidt, 1959; Kennicutt, 1998). Although bridging local and global scales is not straightforward, one may naively argue that the global scale depletion time $\tau_{\mathrm{dep}}$ (or SFE) scales as the local free-fall time (or its inverse). We could therefore suppose that the sSFR is also controlled by density $n_{\mathrm{e}}$ through varying “SFE”, as labeled in Figure 5. We note that the electron density $n_{\mathrm{e}}$ of H ii regions is here assumed to be proportional to the density $\rho_{\mathrm{gas}}$ of star-forming molecular gas clouds, though this is not trivial.

Since volume densities are difficult to observationally measure, the relation between the SFR and gas mass densities is often expressed using the corresponding surface densities as $\dot{\Sigma}_{\ast}\approx\Sigma_{\mathrm{gas}}^{N}$, which is known as the Schmidt-Kennicutt (SK) relation. Observations have shown that the SFR and gas mass surface densities averaged over a galaxy scale obey the SK relation with $N\approx 1.4$ (Kennicutt, 1998). This indicates that the local gas density of collapsing gas clouds scales with total amount of gas content in the systems. We would therefore expect an intrinsic positive correlation between $\mu_{\mathrm{gas}}$ and $n_{\mathrm{e}}$, which is expected to be positive. Although the physical mechanism is not fully clear, this possible causation is labeled as “Compression” for convenience. In addition, the gas density may also be affected by stellar mass $M_{\ast}$ via gravitational effects (labeled as “Gravity”).

Turning to metallicity, the well-known MZ relation could be naturally expected as the consequence of the history of mass assembly and accompanying chemical enrichment. The current metallicity, however, also encloses the history of gas inflows and outflows. In particular, the inverse correlation seen between $Z$ and SFR at fixed $M_{\ast}$ has been commonly interpreted as the result of the so-called “dilution” effect: enhanced inflow of metal-poor gas in the intergalactic medium would dilute the metallicity in the ISM, while boosting star formation. This naturally predicts the anticorrelation between $Z$ and the gas content $\mu_{\mathrm{gas}}$. In Figure 5, these dependencies of $Z$ on stellar mass and the gas content are indicated, respectively, with labels “Enrichment” and “Dilution”. We note that, due to the intrinsic anticorrelation between $M_{\ast}$ and $\mu_{\mathrm{gas}}$ through the gas consuming mass assembly, $Z$ is naturally expected to be anticorrelated with $\mu_{\mathrm{gas}}$ even if we do not consider the effects of dilution due to an enhanced inflow rate. The dilution scenario implicitly assumes that galaxies of the same $M_{\ast}$ experienced the equivalent history of star formation and created the same amount of metals, and that the dispersions in the current gas content and thus metallicity are caused by recent fluctuations in the inflow rate in a short timescale. Keeping this in mind, we consider the effect of $\mu_{\mathrm{gas}}$ on $Z$, especially at fixed $M_{\ast}$, and thus named the causation just "Dilution" for simplicity.

Given $n_{\mathrm{e}}$, $Z$, and sSFR, we can then approach the ionization parameter $U$. Returning to the simple definition of $U\equiv n_{\gamma}/n_{\mathrm{H}}$, we can separate the ionization parameter into two controlling factors, the ionizing photon flux illuminating the gas and the local density of the ionized gas. It would thus be natural to suppose that $U$ is anticorrelated with $n_{\mathrm{e}}$ when fixing the ionizing radiation (labeled as “By definition” in Figure 5). Concerning the numerator, it may be natural to assume that the ionizing photon flux is primarily controlled by the sSFR (labeled as “Radiation” in Figure 5). However, the increase in the galaxy-wide sSFR would not necessarily lead an increase in the ionizing photon flux in a single H ii region, while certainly resulting in an increase in the total amount of the ionizing photons in a whole galaxy. Lastly, we expect that $U$ is also anticorrelated with $Z$. This is because an increase in metals contained in the stellar atmosphere of the central massive stars would result in a higher optical depth in the stellar winds, more efficiently absorbing ionizing photons (e.g., Dopita et al., 2006b). In addition, at higher metallicity, ionizing photons are more efficiently scattered in the photospheres and then their radiation energy are converted to kinetic energy of the stellar winds. We may also expect that an increase in the dust amount in the H ii regions, which is expected to correlate with metallicity, would reduce the ionizing photon flux contributing to ionization (Inoue et al., 2001; Inoue, 2001). The possible causation between $Z$ and $U$ is labeled as “Absorption/Conversion to kinetic energy” in Figure 5).

With these simple ansätze, we attempt to derive the dependencies of sSFR, $Z$, and $U$ on the other parameters. The ansätze denote $\mathrm{sSFR}=\mathrm{sSFR}(\mu_{\mathrm{gas}},n_{\mathrm{e}})$ and $Z=Z(M_{\ast},\mu_{\mathrm{gas}})$. However, we need to eliminate $\mu_{\mathrm{gas}}$ because we have no direct information about the gas content ($\mu_{\mathrm{gas}}$) of the galaxies. We have assumed that $n_{\mathrm{e}}$ is uniquely determined by $M_{\ast}$ and $\mu_{\mathrm{gas}}$ as shown in Figure 5, and thus $\mu_{\mathrm{gas}}$ is uniquely determined once $M_{\ast}$ and $n_{\mathrm{e}}$ are given. We thus parametrize sSFR and $Z$ as a function of $M_{\ast}$ and $n_{\mathrm{e}}$ eliminating $\mu_{\mathrm{gas}}$, as well as $U$ as a function of sSFR $Z$, and $n_{\mathrm{e}}$. In this section, we explicitly argued which positive or negative correlation is expected between two parameters. In our analysis, however, we did not impose any limits on the power-law indices and their signs in the following analysis.

4.2 Multi regression analysis

Following the ansätze we made, we derived the relations between the quantities ($M_{\ast}$, $\mathrm{sSFR}$, $\mathrm{O/H}$, $U$, and $n_{\mathrm{e}}$). The physically-motivated relations to be constrained are expressed as follows:

[TABLE]

We note that $M_{\ast}$, sSFR, and $n_{\mathrm{e}}$ are given in units of $M_{\odot}$, $\mathrm{yr^{-1}}$, and $\mathrm{cm^{-3}}$, respectively. In addition to these, we also derived the empirical functional forms for the three nebular parameters, $Z$, $U$, and $n_{\mathrm{e}}$, as a linear function of $\log M_{\ast}$ and $\log\mathrm{sSFR}$. Our fitting procedure fits a linear surface expressed by these equations to the stacked measurements, accounting for all the measurement errors on the relevant variables.

Let us consider a general form of these models as $y=p_{0}+\sum_{k=1}^{m}p_{k}x_{k}$. We refer to $y$ and $x_{k}~{}(k=1,2,...,m)$ as objective and explanatory variables, respectively. We estimate the model parameters (i.e., the normalization and coefficients of Equations (3–5) based on the maximum likelihood estimation. The likelihood $\mathcal{L}$ is computed for a given model as follows:

[TABLE]

where $N=138$ is the number of stacked measurements. The error $\sigma_{\mathrm{err,i}}$ is obtained for each stack by summing in quadrature the uncertainties expected for the objective ($\delta y_{i}$) and explanatory ($\delta x_{k,i}$) variables as follows:

[TABLE]

Note that the value $\sigma_{\mathrm{err},i}$ depends on the coefficients of the model.

We sampled the posterior probability distribution of the model parameters by employing a Markov-Chain Monte-Carlo (MCMC) technique, and using the emcee package for Python (Foreman-Mackey et al., 2013). In our fitting procedure, we adopted a uniform prior probability function for each parameter.

5 Results

In Table 2, we summarize the coefficients constrained by multi regression analysis. The upper three rows are the results for the physically-motivated models as defined in Section 4.1 (Equations 3–5). The lower three rows show the results obtained from regression on $M_{\ast}$ and sSFR. The reduced $\chi^{2}$ statistics are computed with these best-fit parameters as follows:

[TABLE]

Note that the $(m+1)$ is the number of parameters to be constrained. We also computed the residual intrinsic scatter assuming that the individual stacked measurements are intrinsically scattered with respect to the best fit following a Gaussian distribution of $\sigma_{\mathrm{int}}$. In practice, we estimated $\sigma_{\mathrm{int}}$ that satisfies the following relation:

[TABLE]

The reduced $\chi^{2}$ statistics for the fits with the physically-motivated ansätze are reasonably close to unity, with the residual scatter below $0.1$ dex. This indicates that these models well fit the data with little residual information. For $Z$ and $U$, the $\chi^{2}/\mathrm{dof}$ values are larger for the results from the regression on only $M_{\ast}$ and sSFR. Generally, the $\chi^{2}/\mathrm{dof}$ statistics can be used to evaluate which model is better fitting the same data. However, in the case here, the $\chi^{2}/\mathrm{dof}$ values do not immediately indicate that, for example for $Z$, the fit with $M_{\ast}$ and $n_{\mathrm{e}}$ is better than the fit with $M_{\ast}$ and sSFR because the different models use the different variables. The intrinsic scatters are anyway smaller than $0.1~{}\mathrm{dex}$ for the fits with $M_{\ast}$ and sSFR, similarly to the fits with the ansätze. Hereafter we focus our attention to the results for the physically-motivated models.

We found the best-fit relations for the sSFR, metallicity $Z$, and the electron density $n_{\mathrm{e}}$, as follows (the errors are given in Table 2).

[TABLE]

All these coefficients are well constrained, fully rejecting zero coefficients, i.e., the null hypothesis of no dependence on the parameter. In particular, we highlight that the scaling relation of $U(Z,n_{\mathrm{e}},\mathrm{sSFR})$ is fairly constrained even with three variables. The signs of these coefficients are all consistent with our initial assumption described in Section 4.1 (see Figure 5). Note that the magnitude of the coefficients does not immediately indicate the relative strength of the dependence on the corresponding parameter since the entire range covered by the parameters are different from each other.

In Figure 6, we show the sSFR as a function of $n_{\mathrm{e}}$, color-coded by $M_{\ast}$, in comparison with the derived liner relation of Equation (10) shown in the background. Similarly, Figure 7 shows $Z$ as a function of $M_{\ast}$, color-coded by $n_{\mathrm{e}}$ (upper panel), and as a function of $n_{\mathrm{e}}$, color-coded by $M_{\ast}$ (lower panel). For $U$, we computed the residuals after subtracting the sSFR term ($0.428\left(\log\mathrm{sSFR}+9\right)$). Figure 8 shows these residuals as a function of $Z$, color-coded by $n_{\mathrm{e}}$ (upper panel), and as a function of $n_{\mathrm{e}}$, color-coded by $Z$ (lower panel). In all these panels of Figures 6–8, the overall trends both of the correlation between the $x$- and $y$-axis variables and the color trend of the third axis variables are reasonably reproduced by the best-fit scaling relations. Meanwhile, there are several outliers possibly due to their large errors on the estimates of $n_{\mathrm{e}}$ and/or $Z$. In the next section, we attempt to interpret the global trends and discuss the physical mechanisms behind the scaling relations, not caring about minor deviations from these power-law relations.

6 Discussions

We successfully constrained the dependencies of the physical quantities, $\mathrm{sSFR}(M_{\ast},n_{\mathrm{e}})$, $Z(M_{\ast},n_{\mathrm{e}})$, and $U(\mathrm{sSFR},Z,n_{\mathrm{e}})$, based on the simplified ansätze. In this section, we attempt to interpret the results with our physically-motivated assumptions, and compare our findings to other theoretical predictions and observational results in the literature.

6.1 Specific SFR and gas content

Our result from regressing sSFR on $M_{\ast}$ and $n_{\mathrm{e}}$ indicates that sSFR is a strong function of $n_{\mathrm{e}}$ with a power-law slope of $\approx 1.75$ at fixed $M_{\ast}$. As noted above, this relatively large dependence includes the underlying correlations between $n_{\mathrm{e}}$ and $\mu_{\mathrm{gas}}$. The sSFR of galaxies can be expressed in terms of SFE ($=\mathrm{SFR}/M_{\mathrm{gas}}$) and the gas content as follows:

[TABLE]

As we mentioned in Section 4.1, it would be natural to assume that the SFE is proportional to the inverse of the free-fall timescale $\tau_{\mathrm{ff}}$ of gas clouds, thus we assume $\mathrm{SFE}\propto n_{\mathrm{e}}^{1/2}$. This simple assumption and the observed scaling of sSFR (Equation 10) yield

[TABLE]

Observationally, Sargent et al. (2014) found that SFE increases as a weak function of sSFR: $\mathrm{SFE}\propto\mathrm{sSFR}^{0.19}$, which may be in broad agreement with our result. We can also obtain scaling relations involving the gas content as follows.

[TABLE]

This indicates that the electron density $n_{\mathrm{e}}$ in H ii regions sensitively increases with the gas content at fixed $M_{\ast}$. We also obtain a relation between SFR and $M_{\mathrm{gas}}$ as

[TABLE]

The scaling relation $\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.40}$ (with minor dependence on $M_{\ast}$) agrees very well with observed integrated SK relation based on direct measurements of the gas content of galaxies (e.g., Daddi et al., 2010; Genzel et al., 2010; Tacconi et al., 2013; Santini et al., 2014), as well as the conventional SK relation that correlates the surface densities of SFR and gas mass, for which observations have measured indices $N\approx 1.4$ (e.g., Kennicutt, 1998). It is worth emphasizing that the scaling relation between sSFR and $\mu_{\mathrm{gas}}$ with a slope slightly smaller than 1.5, as supported by other observations, is derived by assuming $\mathrm{SFE}\propto\tau_{\mathrm{tt}}^{-1}\propto n_{\mathrm{e}}^{1/2}$. In doing so, $n_{\mathrm{e}}$ is regarded as being proportional to the density of the molecular clouds where stars form, though this is not trivial because $n_{\mathrm{e}}$ is the measure of the electron density of H ii regions. On the other hand, the agreement between our result ($\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.40}$) and other observations of the SK relations indicate that this assumption would be reasonable.

6.2 Metallicity

We found positive and negative dependences of metallicity $Z$, respectively, on $M_{\ast}$ and $n_{\mathrm{e}}$. In the previous subsection, we obtained the scaling relation between $n_{\mathrm{e}}$ and gas content (Equation 16) based on a simple assumption of $\mathrm{SFE}\propto t_{\mathrm{ff}}^{-1}\propto n_{\mathrm{e}}^{1/2}$ and the observed scaling of sSFR (Equation 10). Using Equation 16, we can eliminate $n_{\mathrm{e}}$ from the $Z(M_{\ast},n_{\mathrm{e}})$ relation (Equation 11), and rewrite $Z$ in terms of $M_{\ast}$ and gas content as follows

[TABLE]

or, more simply,

[TABLE]

The former equation indicates that the correlation between $Z$ and gas-to-stellar mass ratio $\mu_{\mathrm{gas}}$ would be almost independent of $M_{\ast}$ (and other parameters). This is indeed what we expect from simple, classical models of chemical evolution of galaxies (e.g., Schmidt, 1963). The latter form, however, indicates that the metallicity increases with $M_{\ast}$, as naturally expected as the result of the chemical evolution along the galaxy mass assembly history. Meanwhile, the anticorrelation with gas content is as expected from the “dilution” scenario as described in Section 4.1.

Anticorrelations between $Z$ and the gas content of galaxies have been directly measured by observations. Bothwell et al. (2013) analyzed a galaxy sample with the measurement of the H i gas content, and found that the metallicity (via the strong-line method) decreases with the H i gas mass $M_{\mathrm{HI}}$ at fixed $M_{\ast}$, approximately as $Z\propto M_{\mathrm{HI}}^{-0.17}$ (see Figure 2 of the reference). Furthermore, Davé et al. (2013) attempted to model the global properties of the H i gas content in galaxies using cosmological simulations. Their preferable outflow model, which is able to reasonably reproduce the observed H i mass function, predicts that the metallicity is inversely correlated with the H i gas mass with a slope of $-0.26$ at a given $M_{\ast}$. Our result is in reasonable agreement with theses dependencies found by both observations and numerical simulations.

The inverse correlation between $Z$ and sSFR (at fixed $M_{\ast}$) has been also characterized by regressing $Z$ on the $M_{\ast}$–sSFR plane (Table 2 lower part) as follows:

[TABLE]

We note that this relation is nearly equivalent to what is obtained by eliminating $n_{\mathrm{e}}$ from $Z(M_{\ast},n_{\mathrm{e}})$ and $n_{\mathrm{e}}(M_{\ast},\mathrm{sSFR})$ (see Table 2), with a slight difference in the coefficient of $\log M_{\ast}$. The relation between these three parameters ($M_{\ast}$, $Z$, sSFR), i.e., the FMR, has been commonly characterized by introducing a variable $\mu_{\alpha}$ defined as follows:

[TABLE]

where the parameter $\alpha$ is determined to minimize the scatter in terms of metallicity around the average relation. Our result from a simple plane fitting provides $\alpha=0.208/(0.136+0.208)=0.60\pm 0.01$. This is in good agreement with the value ($\alpha=0.66$) found by Andrews & Martini (2013). Meanwhile, such values of $\alpha$ are larger than those ($\alpha\approx 0.3$) as reported by studies relying on the strong-line metallicity indicators (e.g., Mannucci et al., 2010). Note that Equation (20) can be rewritten using $\mu_{0.61}$ as follows:

[TABLE]

The slope of $0.34$ is lower than the slope of 0.43 found by Andrews & Martini (2013), and the difference is not distinguished if we adopt $\alpha=0.66$.

6.3 Ionization parameter

Our results indicate that the ionization parameter $U$ varies as a strong function of the sSFR of galaxies, but also depends significantly on both $Z$ and $n_{\mathrm{e}}$. The former dependence implies that the ionization photon flux in a single H ii region increases on average with increasing total (i.e., galaxy-wide) sSFR of the host galaxies. To directly examine this hypothesis, we need to resolve the individual H ii regions in galaxies across wide ranges of $M_{\ast}$ and sSFR.

There are other ideas to account for the correlations between sSFR and $U$. The observed [O iii]/[O ii] ratio, a proxy of $U$, could effectively increase for a “density-bounded” H ii region, which is fully ionized and whose size is determined by the cloud size, instead of the Strömgren radius. Since O*++* dominates at the inner region of the ionized sphere, close to the radiation source, whereas O*+* becomes more dominant at the outer region, density-bounded H ii regions are expected to have a significantly smaller [O ii] zone, resulting in higher [O iii]/[O ii]. Specially, Binette et al. (1996) considered the effects of density-bounded (“matter-bounded” in the reference) ionized gas clouds to the changes of the emission-line ratios, and showed that the line ratios (especially, [O iii]/[O ii]) vary along the increasing fraction of the density-bounded clouds in a very similar way to the variation with increasing $U$ (see Figure 6 in the reference)222This study specially focused on the conditions in the narrow-line region of active galactic nuclei, but can also applied to star-forming H ii regions.. A larger contribution of density-bounded H ii regions has proposed to account for observed higher [O iii]/[O ii] ratios (Brinchmann et al., 2008; Nakajima & Ouchi, 2014). In addition, it is recently thought that extremely high [O iii]/[O ii] ratios are often incident to non zero escape fraction, i.e., the fraction of ionizing photons that escape into intergalactic space (e.g., Nakajima & Ouchi, 2014; Nakajima et al., 2016; Izotov et al., 2016, 2018). Similarly, if an H ii region overlaps an adjacent H ii region, the effective ionization parameter would increase locally because the overlap region is illuminated by both ionizing sources. It could be naturally expected that such situations tend to happen more frequently in systems with higher sSFRs because the number density of H ii regions may increase with increasing sSFR. We note that our estimate of the ionization parameter assumes a simple ionization-bounded plane-parallel geometry in order to limit the number of free parameters, and that our estimation could thus be inaccurate for some cases. We defer further investigations of the relations between $U$ and the emission-line indicators, such as [O iii]/[O ii], including more realistic geometric conditions to future work.

We now turn to the dependence on metallicity. There have been observational studies that indicate the presence of an anticorrelation between $U$ and $Z$ based on a strong-line method (Nakajima & Ouchi, 2014; Onodera et al., 2016). Shapley et al. (2015) and Kojima et al. (2017) mentioned the anticorrelation by using the stacked measurements based on the direct method from Andrews & Martini (2013). In Figure 4, our data show that $U$ is apparently strongly anticorrelated with $Z$: a simple power-law fit yields $U\propto Z^{-1.52\pm 0.04}$. Compared with this apparent relation, however, our result from the four-dimensional fitting indicates a shallower slope for the partial correlation ($U\propto Z^{-0.36}$; Equation 12) at fixed sSFR and $n_{\mathrm{e}}$. The apparent stronger $Z$-dependence can be attributed to the anticorrelation between sSFR and $Z$ (Section 6.2). Our result thus indicates necessity of physical mechanisms to generate the intrinsic moderate inverse correlation between $U$ and $Z$. Some explanations have been suggested to account for the $U$-$Z$ anticorrelation as mentioned in Section 4.1 (e.g., Dopita & Evans, 1986; Dopita et al., 2006a; Kewley et al., 2013).

Dopita et al. (2006b) have theoretically found that $U$ correlates with $Z$ as $U\propto Z^{-0.8}$ including the following two effects. First, at higher metallicity, ionizing photons emitted by a star are more absorbed by the stellar winds due to an increased optical depth. Second, at higher metallicity, ionizing photons are more efficiently scattered in the photospheres and their radiation energy is more converted to kinetic energy in the stellar wind base region (Leitherer et al., 1992). It may be worth noting that the line-blanketing effects could also be considered as an origin of the $U$–$Z$ anticorrelation: metal elements stored in the stellar photosphere block the radiated ionizing photons, reducing the emergent flux at wavelengths of the metallic lines. However, the back-warming effect occurs immediately following the blocking, where metals re-emit photons at lower frequencies, raising the adjacent continuum. Kudritzki (2002) demonstrated that these effects have little influence on the number of hydrogen ionizing photons ($\geq 13.6~{}\mathrm{eV}$) because the two effects almost entirely cancel each other.

Another explanation is that the $U$–$Z$ correlation is in part caused by dust absorption in the H ii region. Inoue et al. (2001) established a theoretical framework for the dust absorption of ionizing photon in an H ii region, and found observationally an anticorrelation between the fraction $f_{\gamma}$ of the ionizing photons contributing to ionization and metallicity (or dust-to-gas mass ratio) for individual H ii regions in the Milkey Way. Similarly, Inoue (2001) found signs of such anticorrelation also for nearby galaxies, with $f_{\gamma}$ ranging from $\sim 40\textrm{--}70\%$.

The $U$–$Z$ anticorrelation could also be related to the possible correlation between the IMF and metallicity (Evans & Dopita, 1985): a top-heavy IMF would yield more ionizing photons per SFR. Recently, Martín-Navarro et al. (2015) found that the local IMF is tightly related to the local metallicity, becoming top-heavier at lower metallicity. Such correlation may be explained if the fragmentation scales of the molecular clouds are determined by the Jeans mass, which increases with increasing temperature (i.e., towards lower metallicity). A change in the IMF, however, may also change the estimations of $M_{\ast}$ and SFR, as well as the ionization parameter because the shape of the ionizing spectrum changes. We defer further investigation including the possible dependence of the IMF to future work.

Turning to the $n_{\mathrm{e}}$-dependence, our four-dimensional fitting found a negative partial correlation between $U$ and $n_{\mathrm{e}}$ as $U\propto n_{\mathrm{e}}^{-0.29}$, at fixed sSFR and $Z$, as opposed to their apparent positive correlation (with the Spearman’s rank correlation coefficient $\rho=0.476$ with $p$-value $<10^{-8}$) as shown in the relevant panel in Figure 4. The apparent correlation is caused by the positive correlation between sSFR and $n_{\mathrm{e}}$ and the tight correlation between $U$ and sSFR. Indeed, the bottom-leftmost panel in Figure 4 shows the presence of an inverse correlation between $U$ and $n_{\mathrm{e}}$ at fixed sSFR (i.e., at fixed color). The observed dependence is weaker than what is expected from the definition of $U$, i.e., $U\propto n_{\mathrm{e}}^{-1}$. This may indicates that the differential electron densities between relatively inner regions, where the ionization parameter is defined, and outer regions, which may dominate the [S ii] fluxes, in the H ii regions.

The weak negative dependence of $U$ on $n_{\mathrm{e}}$ is in good agreement with what was derived by Dopita et al. (2006b), who found that the theoretical ionization parameter is weakly related to the density $n_{0}$ of the ambient ISM as $U\propto n_{0}^{-1/5}$ by modeling the evolution of an expanding H ii region, pressurized by the stellar winds from the central star cluster. In their model, the internal density $n_{\mathrm{in}}$ of an isothermal spherical H ii region is related to the ambient density, as $n_{\mathrm{in}}\propto n_{0}^{3/5}$, as well as the inner radius $R_{\mathrm{in}}$ is also linked to it as $R_{\mathrm{in}}\propto n_{0}^{-1/5}$ (Oey & Clarke, 1997). For a fixed ionizing luminosity of the central source, the ionization parameter depends on these quantities as $U\propto 1/(n_{\mathrm{in}}R_{\mathrm{in}}^{2})$, yielding the relation between $U$ and the ambient density as above.

Lastly, we note that we recently obtained a conclusion that the ionization parameter $U$ is directly linked to sSFR, while little or probably not with metallicity (Kaasinen et al., 2018), from comparison between the high-redshift sample of the main-sequence star-forming galaxies at $z\sim 1.6$ (a subset from Kashino et al. 2017, 2018) and the “sSFR-matched” sample of the SDSS galaxies. These two samples that have the equivalent average sSFRs present similar levels of the ionization parameter, whereas the sSFR-matched low-redshift sample shows relatively lower metallicities. However, our results in this paper are indeed not against what is found by Kaasinen et al. (2018) within the uncertainties. The difference in the metallicity of these two samples is found to be $\sim 0.1~{}\mathrm{dex}$. From our result, this metallicity difference induces only a very small difference ($\lesssim 0.05~{}\mathrm{dex}$) in $U$ at fixed sSFR, which is equivalent to the level of statistical uncertainties presented in Kaasinen et al. (2018). We also note that our analysis has been applied only for the low redshift galaxies, whereas the attention in Kaasinen et al. (2018) was more focused on the evolution over a wide range of cosmic time. The investigation of redshift evolution of the relationships that we found between the nebular parameters is beyond the scope of this paper.

6.4 Application of the scaling relations

Lastly we mention possible application of the scaling relations we derived in this paper. Our results can be implemented into a subgrid recipe of galaxy evolution simulations to compute the conditions of ionized nebular from the properties of host galaxies, and predict the fluxes of emission lines of interest. Predicting strong-line fluxes is useful to prepare mock catalogs for future spectroscopic surveys as well as to test the models by comparing observations. Our results would be also implemented into SED fitting codes to include the contribution from emission lines for given galaxy properties.

Here we present an application of our scaling relations for a semi-analytic model galaxy catalog. In semi-analytic models, the evolution of galaxies embedded in collapsed dark halos, is traced following the halo merger trees while solving a simplified model that controls the star formation involving relevant physical phenomena such as supernova and/or AGN feedback. As a result, some fundamental quantities of the galaxies, such as $M_{\ast}$, $M_{\mathrm{gas}}$, SFR, metal abundance and the stellar population, are obtained as functions of cosmic time. Once $M_{\ast}$ and $M_{\mathrm{gas}}$ are given, our scaling relations can be utilized to compute the conditions of the H ii regions, $n_{\mathrm{e}}$ and $U$, from Equations (16) and (12), respectively. With the gas-phase $Z$ traced in the simulation, it is then possible to compute the emission line fluxes by interpolating a prepared table of the line fluxes for grids of ($n_{\mathrm{e}},Z,U$).

For a demonstration, we use the New Numerical Galaxy Catalog ($\nu^{2}$GC; Ishiyama et al. 2015; Makiya et al. 2016)333http://hpc.imit.chiba-u.jp/ nngc/index.html. From the public light-cone catalog, we extracted a sample of 4855 galaxies limiting to those at $z<0.3$ and having $\mathrm{sSFR}\geq 10^{-11}~{}\mathrm{yr^{-1}}$. We define the stellar mass as the sum of disk (Mstard) and bulge (Mstarb) masses, and take the gas mass ($M_{\mathrm{gas}}=$Mcool) and metal mass ($M_{Z}=$MZc) in the cold phase. The metallicity is defined as $12+\log(\mathrm{O/H})=8.69+\log(M_{Z}/M_{\mathrm{gas}}/Z_{\odot})$ where $Z_{\odot}=0.014$. In Figure 9, we show the correlations among $M_{\ast}$, sSFR, $Z$, given by the semi-analytic model, and $n_{\mathrm{e}}$ and $U$ computed with our scaling relations. In the $\nu^{2}$GC sample, galaxies with $\mathrm{sSFR}\geq 10^{-8.5}~{}\mathrm{yr^{-1}}$ are missing, while a space of low $M_{\ast}$ ($<10^{9}~{}M_{\odot}$) and low sSFR ($<10^{-9.5}~{}\mathrm{yr^{-1}}$) is filled up, where we have no SDSS stacks (see Figure 4). Moreover, a population of very metal-rich galaxies ($12+\log\mathrm{O/H}>9$) exits. Regardless of these discrepancies, the trends of the SDSS sample seen in Figure 4 are globally reproduced here, manifesting the effectiveness of the application to semi-analytic model galaxy samples. Detailed comparison of the $\nu^{2}$GC catalog with observations and/or predictions of any concrete observables are beyond the scope of this paper.

7 Conclusions

We have investigated the connection between the physical parameters characterizing the conditions of star-forming ionized gas and the global properties of the host galaxies, especially focusing on five key parameters, i.e., stellar mass, sSFR, oxygen abundance, ionization parameter, and electron density of the ionized gas. In order to obtain accurate estimates of the oxygen abundance, we have utilized the direct method based on the detection of auroral lines ([O iii]$\lambda$4363, and [O ii]$\lambda\lambda$7320,7330) in the set of stacked spectra of the SDSS galaxies at $z\sim 0.1$, which are binned into the $M_{\ast}$–sSFR grids. Our key conclusions can be summarized as follows.

The gas-phase oxygen abundance was successfully measured for a wide range of $M_{\ast}$ and sSFR by using stacked spectra. The resultant values range between $7.6\lesssim 12+\log(\mathrm{O/H})\lesssim 8.9$, tightly correlating with sSFR and $M_{\ast}$. 2. 2.

The electron density $n_{\mathrm{e}}$ was measured from the [S ii] doublet ratio, accounting for the dependence on the electron temperature, and found to increase with increasing sSFR at fixed $M_{\ast}$ as well as with increasing $M_{\ast}$ at fixed sSFR. 3. 3.

The ionization parameter $U$ was measured from the [O iii]/[O ii] ratio, accounting for the dependencies on $Z$ and $n_{\mathrm{e}}$. The resultant values of $U$ were found to be tightly correlated with the sSFR, and also to depend on $Z$ and $n_{\mathrm{e}}$. 4. 4.

The partial anticorrelation between $U$ and $Z$ ($U\propto Z^{-0.36}$ at fixed sSFR and $n_{\mathrm{e}}$) was found to be weaker than the apparent correlation ($U\propto Z^{-1.52}$). The observed moderate (i.e., not as strong as its apparent correlation) anticorrelation between $U$ and $Z$ is qualitatively consistent with theoretical predictions of the intrinsic $U$–$Z$ relation. 5. 5.

The ionization parameter is also found to negatively correlate with $n_{\mathrm{e}}$ at fixed sSFR and $Z$ as $U\propto n_{\mathrm{e}}^{-0.29}$. This is qualitatively consistent with the definition of the ionization parameter, though the dependence on $n_{\mathrm{e}}$ is weaker than the inverse proportion. This may indicates that the electron density probed by the [S ii] doublet scales with the inner surface density more sensitively than proportionality.

Our results presented in this paper have firmly provided the global picture on the behavior of the nebulae at the galaxy-wide scales. The scaling relationships that we found between the fundamental properties of galaxies and the nebulae provide a set of boundary conditions which are to be reproduced by high-resolution simulations. On the other hand, our results provide a convenient tool to compute the nebular conditions and to predict fluxes of strong emission lines with little computational cost in semi-analytic models and/or SED fitting codes. Such application is thus useful to construct mock catalogs for forthcoming large spectroscopic surveys with the next-generation spectrographs, such as the the Subaru/Prime Focus Spectrograph (PFS), Dark Energy Spectroscopic Instrument (DESI), Multi-Object Optical Near-infrared Spectrograph (MOONS), and Euclid and WFIRST in space. These next-generation surveys will then be enable us to extend our analysis beyond the low redshift Universe in the near future, and to understand the physical processes behind the connections of galaxy properties in more detail.

Acknowledgements

We are grateful Brett Andrews for kindly providing us with their catalog and data. We also thank K. Yabe and T. Kojima for useful discussions and comments. This work has been supported by KAKENHI (17H01114) and by NAOJ ALMA Scientific Research Grant number 2016-01 A.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

Appendix A Balmer absorption

Here we present empirically-calibrated formulae for correcting the observed fluxes of the Balmer lines (H$\alpha$ and H$\beta$) for the reduction due to the stellar atmospheric absorption as a function of $M_{\ast}$ and sSFR. To evaluate the effects of the Balmer absorption, we compare the flux measurements on the composite spectra in which the stellar continuum is subtracted using pPXF (Section 2.3) with those based on the continuum subtraction using a simple linear fitting, as usually done for continuum-faint or low-S/N spectra. From these measurements, we defined the fraction of the absorbed flux $f_{\mathrm{abs}}$ for each stack, $f_{\mathrm{abs}}=(F_{\mathrm{int}}-F_{\mathrm{obs}})/F_{\mathrm{int}}$, where $F_{\mathrm{obs}}$ and $F_{\mathrm{int}}$ are the flux measured, respectively, with the linear-fit based and pPXF-based continuum subtraction.

In Figure 10, we show the absorbed fraction of the H$\beta$ line in the $M_{\ast}$–sSFR plane. Stacks are limited to those with $S/N(\mathrm{H\beta})\geq 40$. Each box represents a single stack, color-coded by $f_{\mathrm{abs}}(\mathrm{H\beta})$. It is clear that a strong trend exists between $f_{\mathrm{abs}}$ and a quantity, which can be defined as a combination of $M_{\ast}$ and SFR. The absorption fraction of the H$\alpha$ flux exhibits a similar trend to H$\beta$, but the absolute values are much smaller. Considering a linear combination as $x_{\alpha}=\log\mathrm{SFR}-\alpha\log M_{\ast}$, we found that ,when $\alpha=1.32$, the correlation becomes the tightest for both H$\alpha$ and H$\beta$. The relation $x_{\mathrm{\alpha=1.32}}=0$ is shown by a red solid line. In Figure 11, we show the absorption fraction of the H$\alpha$ and H$\beta$ fluxes as a function of $x_{\alpha=1.32}$. To express these trends, we fit an empirical functional form to data, and obtained the forms as follows,

[TABLE]

Appendix B Nebular diagnostic diagrams

It would be useful to see the relation between the strong-line flux ratios and the galactic and nebular properties. In Figure 12, we show four routinely-used diagrams – $R_{23}$ vs. $O_{32}$, $O_{2}$ vs. $O_{3}$, $N_{2}$ vs. $O_{3}$, and $S_{2}$ vs. $O_{3}$ diagrams (the so-called BPT diagrams; Veilleux & Osterbrock 1987). Here the line ratio indices are defined as follows: $R_{23}=(\textrm{[O\,{\sc ii}]}\lambda 3727+\textrm{[O\,{\sc iii}]}\lambda\lambda 4959,5007)/\mathrm{H\beta}$, $O_{32}=\textrm{[O\,{\sc iii}]}\lambda 5007/\textrm{[O\,{\sc ii}]}\lambda 3727$, $O_{2}=\textrm{[O\,{\sc ii}]}\lambda 3727/\mathrm{H\beta}$, $O_{3}=\textrm{[O\,{\sc iii}]}\lambda 5007/\mathrm{H\beta}$, $N_{2}=\textrm{[N\,{\sc ii}]}\lambda 6584/\mathrm{H\alpha}$, and $S_{2}=\textrm{[S\,{\sc ii}]}\lambda\lambda 6717,6731/\mathrm{H\alpha}$. The contours show the distribution of the individual SDSS galaxies in the parent catalog. To draw the contours, we additionally impose on the sample the detection of [O iii]$\lambda$5007 at $>3\sigma$ (89 percent of the sample), although the original selection does not care about the detection of [O iii]$\lambda$5007 (Section 2.1). The line ratios of the individual sources are corrected for extinction using the Cardelli et al. (1989) reddening curve and assuming the intrinsic H$\alpha$/H$\beta=2.86$. The SDSS stacks with the successful metallicity estimation are shown by circles, color-coded by $M_{\ast}$, SFR, sSFR, $12+\log(\mathrm{O/H})$, $U$, and $n_{\mathrm{e}}$ in each panel.

We can see clear trends of the location in the diagrams against all the average physical properties considered here. The trends with $M_{\ast}$ are consistent with those of the individual galaxies (Masters et al., 2016, see Figure 2). At a fixed SFR, the stacks cover a wide range of the emission-line ratios, almost from the one to the other tail-end in all four diagrams (second row from top in Figure 12), forming a narrow sequence parallel to the ridge line of the contours. The nebular conditions are also tightly correlated with the location in the diagram. The metallicity varies along the sequence, from the lowest $Z$ at the high end-tail of the $y$-axis values to the highest $Z$ at the opposite side (forth row in Figure 12). The trend with $U$ is very similar to sSFR (fifth row). This is corresponding to the tight correlation between $U$ and sSFR which we have seen in Figure 4. The electron density $n_{\mathrm{e}}$ is correlated with the location in a similar way to SFR: stacks of fixed $n_{\mathrm{e}}$ covers a wide range of these strong-line ratios and form a narrow sequence parallel to the ridge line of the contours. The trend in the $N_{2}$ vs. $O_{3}$ diagram was reported by Brinchmann et al. (2008).