© ESO 2022

1 Introduction

As one of the most widely observed and studied classes of ionized regions in galaxy studies, H II regions are vital for our understanding of galaxy evolution. By studying the emission line spectra of H II regions, one can learn about their chemical compositions as well as their ionization states, which are related to the past and current evolution statuses of their host galaxies. Various methods have been developed to derive the properties of H II regions based on their emission line spectra, among which the so-called strong line method is widely used in the literature. The strong line method makes use of strong optical emission lines that are easy to observe to predict the gas-phase metallicity and ionization parameter. There are in general two ways of applying this method, including empirical calibration and theoretical modeling. The former takes advantage of the observations of nearby H II regions that have very high signal-to-noise ratios (S/N), of which the metallicities can be determined by the direct method that uses faint optical auroral lines. The derived values can then be compared with the strong line ratios to obtain an empirical relation, which could later be applied to distant H II regions or star-forming (SF) galaxies without measurements of auroral lines (e.g., Pagel et al. 1979, 1980; Edmunds & Pagel 1984; Pettini & Pagel 2004; Pilyugin & Thuan 2005; Marino et al. 2013). Although the empirical method does not rely on any assumption on the physical conditions of H II regions, there are unavoidable intrinsic scatters in most of the relations introduced by the variations in parameters other than metallicity (which we refer to as “secondary parameters" hereafter). In addition, the ionization parameter, which is defined as the ratio between the flux of the hydrogen ionizing photons to the hydrogen density, is almost impossible to be measured directly. It cannot be inferred from empirical calibrations and has to depend on certain assumptions on the internal structures of H II regions.

Theoretical modeling, on the other hand, can self-consistently predict both the metallicity and ionization parameter, provided that the secondary model parameters are properly set and match the realistic H II regions (e.g., Charlot & Longhetti 2001; Kewley & Dopita 2002; Tremonti et al. 2004; Kobulnicky & Kewley 2004; Dors et al. 2011). These parameters include the stellar spectral energy distribution (SED), the chemical abundance pattern, the dust composition, the density structure, the geometry, etc. State-of-the-art photoionization codes such as CLOUDY and MAPPINGS can be used to compute 1D models of ionized clouds based on existing atomic data and to predict emission-line ratios that are in good agreement with the observations (e.g., Dopita et al. 2000, 2013; Kewley et al. 2001, 2019, and references therein). Overall, photoionization models have been shown to be able to successfully reproduce several important emission-line ratios in observations. For example, the model predictions on [N II]λ6583/Hα, [S II]λλ6716, 6731/Hα, and [O III]λ5007/Hβ are largely consistent with observations, which explains the observed SF sequence in standard optical diagnostic diagrams (Stasińska et al. 2006; Dopita et al. 2013). However, it is important to note that when the model-predicted line ratios are compared with the observations to find the metallicity and ionization parameter, the results have a strong dependence on the input parameters. Unfortunately, many of the secondary model parameters are hard to constrain, and their importance is often overlooked in practice. This can lead to large discrepancies among the strong-line relations based on different models, as we show in this paper. A prime example of strong model-dependence can be shown by the results on the correlation between the metallicity and ionization parameter.

The metallicity and ionization parameter are often assumed to be the main contributors to the variations of line ratios among SF regions (for a counterargument, however, see Pellegrini et al. 2020). An interesting question is whether they are correlated with each other or independent. Surprisingly, there are a considerable number of contradictory findings in the literature on this issue. Some early works show the two quantities are clearly anti-correlated (e.g., Dopita & Evans 1986; Maier et al. 2006; Nagao et al. 2006), which is further supported by the theoretical calculation of a wind-driven bubble model for H II regions by Dopita et al. (2006) (hereafter D06). D06 argued that there are two effects driving this anticorrelation. First of all, at higher metallicities, stellar winds become more opaque and absorb more ionizing photons. In addition, the stellar atmosphere scatters photons more effectively, leading to stronger winds. This enlarges the inner shocked wind region and dilutes the ionizing flux received by the outer H II region. The anticorrelation has also been confirmed by a number of more recent studies using a relatively large sample of H II regions or SF galaxies (e.g., Pérez-Montero 2014; Morisset et al. 2016; Thomas et al. 2019).

Interestingly, some other studies found that the metallicity and ionization parameter either have no obvious correlation, or have a positive correlation in observed H II regions (e.g., Dors et al. 2011; Dopita et al. 2013, 2014; Poetrodjojo et al. 2018; Kreckel et al. 2019; Zinchenko et al. 2019; Mingozzi et al. 2020), in contrast to the theoretical prediction of D06. Dopita et al. (2014) (hereafter D14) propose that a positive correlation between the metallicity and ionization parameter only exists in starburst galaxies and not in normal SF galaxies. They argue that the underlying reason for the positive correlation in starburst galaxies could be the relation between the star-formation rate (SFR) and ionization parameter. A higher SFR could increase the ionization parameter either due to the higher mass of the star cluster or the change of the geometry of the ionized cloud. However, according to Telford et al. (2016) and Mingozzi et al. (2020), such a relation between the SFR and ionization parameter can also be found in normal SF galaxies. In addition, their analyses show that the correlation strongly depends on the stellar masses of the galaxies. Furthermore, Mingozzi et al. (2020) found that the ionization parameter is actually more tightlycorrelated with the specific star-formation rate (sSFR, which is traced by the equivalent width of the Hα line in their work) than the SFR. To make the matter even more convoluted, Poetrodjojo et al. (2018) compared the slopes of the metallicity versus ionization parameter relations (MI relations) and the SFR versus ionization parameter relations (SFR-I relations) for H II regions in different galaxies. They found that the slopes of these two relations are very different for different galaxies and the overall correlations between these quantities are not significant.

The tension between the two contradictory kinds of findings is puzzling and this remains unresolved. Kewley et al. (2019) comment that the evidence against the expected anticorrelation mainly came from spatially resolved studies. However, as we show in this paper, the correlation between the metallicity and ionization parameter in observed H II regions or SF galaxies is less related to the scale of the observations, and it is rather a result of model assumptions and methodology. In this work, we reexamine this long-standing issue by focusing on the choice of input parameters for photoionization models that are used to fit the metallicity and ionization parameter of H II regions. We provide a self-consistent model that resolves this problem, which could help reveal the potential physical process that leads to the correlation. The layout of the paper is as follows. In Sect. 2, we describe the observational data we use and the input parameters for the photoionization models. In Sect. 3, we compare the model predictions on the correlations under different assumptions and evaluate the models based on the consistency of their predictions on different combinations of emission line ratios. We discuss the physical interpretations for the correlation indicated by our best-fit model in Sect. 4 and examine the robustness of our analyses in Sect. 5. Finally, we summarize our findings and draw our conclusions in Sect. 6.

Throughout this work, we use the following abbreviations for some of the frequently mentioned emission line ratios. We denote log([N II]λ6583∕Hα), log([S II]λλ6716, 6731∕Hα), log([O III]λ5007∕Hβ), log([N II]λ6583∕[O II]λλ3726, 3729), and log([O III]λ5007∕[O II]λλ3726, 3729) as N2, S2, R3, N2O2, and O3O2, respectively.

2 Data and models

In this work, we use observational data from the Mapping Nearby Galaxy at Apache Point Observatory survey (MaNGA, Bundy et al. 2015). As one of the three major experiments of SDSS-IV (Blanton et al. 2017), MaNGA was designed to obtain spatially resolved spectroscopic data of ~10 000 nearby galaxies andit finished its observation in the summer of 2020. Its observing strategy and survey execution are detailed in Law et al. (2015) and Yan et al. (2016a). With a median redshift of 0.03, MaNGA’s targets form a primary sample made up of galaxies observed out to 1.5 R_e and a secondary sample made up of galaxies observed out to 2.5 R_e. These targets are designed to have a nearly flat stellar mass distribution with M_* = 10⁹ ~ 10¹¹ M_⊙ (Wake et al. 2017). MaNGA uses the 2.5 m Sloan telescope for its observation (Gunn et al. 2006). Lights from galaxies are fed through IFU fiber bundles with fields of view ranging from 12′′ to 32′′ (Drory et al. 2015) to the BOSS spectrographs (Smee et al. 2013). The collected spectra have a median spectral resolution of R ~ 2000 and cover a wavelength range from 3622 Åto 10 354 Å. These raw spectra data are then reduced and calibrated by the Data Reduction Pipeline (DRP, Law et al. 2016, 2021b; Yan et al. 2016b) and eventually fed to the Data Analysis Pipeline (DAP, Belfiore et al. 2019; Westfall et al. 2019). The DAP uses the PENALIZED PIXEL-FITTING software (PPXF, Cappellari & Emsellem 2004; Cappellari 2017) at its core and fits the stellar continua and emission line spectra simultaneously. The final data products include spatially resolved models and measurements of emission lines and stellar continua. A python-based toolkit, MARVIN, further facilitates steamlined access and visulization of the DAP products (Cherinka et al. 2019).

Our sample is drawn from the MaNGA products in the 15th public data release of SDSS (DR15), which includes a total of 4639 unique galaxies. This data set is equivalent to the seventh product launch of MaNGA (MPL-7). We selected a sample of H II regions by using the 3D optical diagnostic diagram introduced by Ji & Yan (2020), which combines the [N II]- and [S II]- Baldwin, Phillips & Terlevich (BPT) diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987) and selects spaxels in a 3D space. Basically, Ji & Yan (2020) used two model surfaces to describe the SF and active galactic nucleus (AGN) loci in 3D space, which are shown in Fig. 1. The two model surfaces are well separated in 3D, allowing thedecomposition of spatially mixed ionized regions into pure ionized regions, of which the line ratios are indicated by the locations of the model surfaces. With this method, Ji & Yan (2020) computed quantitatively a demarcationsurface in 3D that is able to isolate SF regions with an AGN contribution of less than 10% to the total Hα fluxes¹. The selection function is given by $$P_1<-1.57P_2^2+0.53P_2-0.48,$$(1)

where $$P_1=0.63N2+0.51S2+0.59R3,$$(2)

and $$P_2=-0.63N2+0.78S2.$$(3)

Figure 1 shows the density distribution of our sample in the 3D space spanned by N2, S2, and R3. It is clear that there is a continuous mixing sequence connecting the two model surfaces, and our selected sample is well traced by the fiducial SF model surface. The sample includes a total of 2782 galaxies and 1.65 million spaxels. Figure 2 shows the data distribution in a specific 2D projection, which we denote as the P₁ –P₂ diagram, where the photoionization model surfaces appear edge-on. Since both the data and the (relevant) model surface cover a very narrow region in this projection, one can use it to visualize how well the model surfaces fit the center of the data distribution. In addition, the projected AGN model grid is clearly separated from the SF model grid, indicating that AGN or composite regions are unlikely to be projected into the SF locus and contaminate the sample significantly. Although the two models seem to overlap with each other at the upper part of the diagram, the corresponding part of the AGN model actually describes very low metallicity AGNs, which are rare in the MaNGA sample, as indicated by the spaxels in the central regions of MaNGA galaxies (Ji & Yan 2020).

Since we have tried to carry out a statistical study on the distribution of both the metallicity and ionization parameter in “typical” H II regions for this work, whether the selected sample is representative enough has a non-negligible impact on the result. For our main MaNGA sample, the mass distribution is relatively flat, as can be seen in Fig. 3. This ensures that the high-mass regime is well populated and helps to constrain photoionization models at high metallicities. We defer a more detailed discussion on the selection effect to Sect. 5.

Our fiducial photoionization model for H II regions (hereafter JY20 model) is generated by the photoionization code CLOUDY (v17.00, Ferland et al. 2017), and it is described in detail by Ji & Yan (2020). In short, this model simulates an isobaric H II region with plane-parallel geometry. The initial hydrogen density is set to be 14 cm⁻², which is derived from the median [S II]λ6716/[S II]λ6731 of H II regions in MaNGA (Ji et al. 2020). The ionizing SED is produced by the code STARBURST99 (v7.01, Leitherer et al. 1999), assuming a continuous star-formation history (SFH) of 4 Myr and a Kroupa initial mass function (IMF, Kroupa 2001). When computing this SED, we used the Pauldrach et al. (2001) and Hillier & Miller (1998) stellar atmospheres, and a standard Geneva evolutionary track. To account for secondary nitrogen, we used the N/O versus 12 + log(O/H) relation given by Dopita et al. (2013), which was derived by fitting a set of H II-region measurements derived by van Zee et al. (1998). For the gas-phase chemical abundances, we chose the solar abundance set of Grevesse et al. (2010) as the reference abundance and used the default dust depletion factors in CLOUDY, which are based on measurements by Cowie & Songaila (1986) and Jenkins (1987). For each gas-phase abundance used in the model, we tried to match it with the stellar SED with the same metallicity, assuming the abundance ratio among heavy elements (including α-elements but excluding C and N) is the same as the Sun. Since STARBURST99 only computes stellar metallicities up to 2 Z_⊙, we extrapolated the stellar SED when the gas-phase oxygen abundance was larger than two times the solar value. Basically, we assumed that the logarithms of the stellar fluxes change linearly with the logarithmic oxygen abundances².

Besides this fiducial model, we also examined three other SF-ionized models in the literature, computed by Levesque et al. (2010), Dopita et al. (2013), and Byler et al. (2017), respectively (we denote these models as L10, D13, and B17 hereafter). The key input parameters of these models, including the hydrogen density, geometry, ionizing SED, solar abundance set, nitrogenprescription, and dust depletion, are listed in Table 1. The relative importance of these parameters is further discussed in Sect. 3. Throughout this work, we express the gas-phase metallicity as 12 + log(O/H), which refersto the “pre-depletion” abundance unless otherwise specified. We chose 12 + log(O/H) instead of [O/H] (≡ log((O/H)/(O/H)_⊙)) adopted by these models to avoid any confusion about the value of the solar oxygen abundance. For the stellar metallicity, we denote it as Z. For all photoionization models in Table 1, the relative stellar metallicity, log (Z∕Z_⊙), and the relative gas-phase metallicity, [O/H], were set to be the same and were varied together. The ionization parameter is defined as $$U=\frac{{\Phi }_0}{n_{\text{H}}c}=\frac{Q_0}{4\pi r^2{n}_{\text{H}}c},$$(4)

where Φ₀ is the ionizing flux from the central stars, Q₀ is the number of ionizing photons per unit time, r is the radius at which the ionization parameter is measured³, n_H is the hydrogen density, and c is the speed of light.

Fig. 1 Density distribution of MaNGA MPL-7 spaxels in the 3D space spanned by N2, S2, and R3. Our sample H II region spaxels are colored from yellow to purple, while the rest of the spaxels are colored from white to black. Two photoionization model surfaces are shown. The cyan model and the red model are our fiducial SF model and AGN model, respectively. We show three viewing angles that lead to different 2D projections. Specifically, panel a and panel b correspond to the [S II]- and [N II]- BPT diagrams. Panel c shows that the two model surfaces are separate in 3D and that they are connected by a continuous mixing sequence.

Fig. 2 Density distribution of the MaNGA data in the 2D P1–P2 diagram, where the relevant parts of the (interpolated) SF model (cyan grid) and AGN model (red grid) appear edge-on and well separated. The sample H II regions are colored from yellow to purple. Here, we only plotted the parts of the model surfaces that cover the middle 98% of the data along the P3 axis (line of sight), which is perpendicular to the P1 versus P2 plane. This projection corresponds to a line of sight at (θ, ϕ) = (36°, 219°) in the 3D space of (N2, S2, R3), where θ and ϕ are the polar angle and the azimuthal angle, respectively (Ji & Yan 2020).

Fig. 3 Redshift and stellar mass distributions of our sample galaxies. The histogram shows a relatively flat distribution in stellar mass. The sample galaxies include both primary and secondary MaNGA galaxies. The former cover spatial regions out to 1.5 Re, while the latter cover spatial regions out to 2.5 Re.

3 Dependence of the MI correlation on model parameters

In this section, we use photoionization models to fit gas-phase metallicities as well as ionization parameters for our sample. The basic idea ofthis approach is to compare the model predictions on multiple emission line ratios with the observed values. Depending on the number of line ratios involved, the calculation could be 2D (minimal dimensions required to fit both the metallicity and ionization parameter), or it could have three or more dimensions.

Starting with the 2D diagnostics: the standard 2D BPT diagrams ([N II]-, [S II]-, and [O I]-based diagrams) can provide useful constraints on the ionizing sources for the ionized regions, but they are not suitable for fitting metallicities and ionization parameters. The photoionization model grids tend to wrap around in these diagrams, which results in significant degeneracies between the fitted metallicities and ionization parameters. Another 2D diagnostic diagram that is frequently used to fit metallicities and ionization parameters simultaneously is the [N II]/[O II] versus [O III]/[O II] diagram (Dopita et al. 2000). The merit of this diagnostic diagram is that photoionization model grids do not overlap with themselves in this diagram⁴. However, a few caveats should be noted. First, this diagram has an extra dependence on the dust extinction correction due to the large wavelength separation of the emission lines used. Second, it relies on a good understanding of how the N/O ratio changes with the metallicity. The fact that there is a significant degeneracy between models with different N/O versus O/H relations makes this diagram alone not very useful in constraining the model parameters or providing unique fitting results.

To overcome the above difficulties associated with the 2D diagrams, one could take more line ratios into consideration and constrain the secondary model parameters before fitting the metallicity and ionization parameter. This approach effectively brings us to the regime of high-dimensional analysis. Ji & Yan (2020) show that by simply combining [N II]/Hα, [S II]/Hα, and [O III]/Hβ and forminga 3D diagnostic diagram (N2–S2–R3 diagnostic), one can place strong constraints on secondary parameters, including the shape of the ionizing SED and the N/O abundance pattern. This is achieved by simply requiring that the model surface goes through the densest part of the data surface in 3D. The N2–S2–R3 diagnostic itself can also be used to predict the metallicity and ionization parameter. Here we would like to define the concept of “best fitting” in high dimensions. Ideally, there exists an optimal photoionization model, which gives consistent predictions over all the line ratios, should the physical properties of the studied H II regions all have narrow distributions about the median values. Still there are scatters around the manifold of the model due to the widths of the distributions of different parameters. Following the previous practices on modeling H II regions, we make the following assumptions.

First, there are two primary parameters, that is to say the metallicity and the ionization parameter, that contribute the largest line ratio variations. Second, at a fixed metallicity and ionization parameter, the variations in other (secondary) parameters result in comparably less variations in line ratios for a statistically large sample of H II regions.

By varying the metallicity and ionization parameter of the photoionization model, we obtain a 2D surface in any given line-ratio space. The variations in secondary parameters appear as small scatters around the model surface in high dimensions, but the best-fit model surface still lies in the center, or the densest part of the data distribution. The location of the “central surface” can then be used to constrain the median values of the secondary parameters. We expect the desired model that best fits the central surface of the data to produce identical relationships between the metallicity and ionization parameter, within the uncertainties, no matter which subset of line ratios is used to derive them (except for those subsets with significant degeneracy). Whether the model predicts the same MI correlation using different line ratios is thus a useful test for its self-consistency.

In the following part of this section, we compare the fitting results from different photoionization models. We used the Bayesian inference to estimate the joint and marginalized probability distribution functions (PDFs) of the logarithmic metallicity, [O/H], and ionization parameter, log(U), for each of the data points, assuming flat priors in both dimensions (see e.g., Blanc et al. 2015, for a general discussion on the Bayesian approach). Basically, for each data point, we calculated the likelihood as $$p(D|M,\theta )=\frac{\text{exp}(-{\chi }^2/2)}{{(2\pi )}^{n/2}\text{det}{(C)}^{1/2}},$$(5)

where D, M, and θ represents the data, the adopted model, and the model parameters, respectively (i.e., [O/H] and log(U)). Furthermore, n is the number of emission line ratios considered, and χ² is given by $${\chi }^2={[X_D-X_M]}^{\text{T}}{C}^{-1}[X_D-X_M],$$(6)

where X_D −X_M is an n dimensional vector describing the logarithmic line-ratio differences between the data point and the model point, and C is the covariance matrix (Hogg et al. 2010). Following Belfiore et al. (2019) who examined the uncertainties reported by DAP using repeated observationsin MaNGA, we inflated all the uncertainties in emission line fluxes by a factor of 1.25. This is equivalent to multiplying all terms in C by 1.5625, which does not have a significant effect on our results. The likelihood was then combined with the flat priors and normalized to give the posterior, p(θ|D, M) (we discuss the influence of nonflat priors in Sect. 5). Finally, we calculated the weighted average metallicity and ionization parameter for each data point using the marginalized posteriors⁵ $$<\log (\text{O/H})>={\displaystyle \sum _{\log {(\text{O/H})}_{\text{min}}}^{\log {(\text{O/H})}_{\text{max}}}p}(\log (\text{O/H})|D,M)\log (\text{O/H}),$$(7) $$<\log (U)>={\displaystyle \sum _{\log {(U)}_{\text{min}}}^{\log {(U)}_{\text{max}}}p}(\log (U)|D,M)\log (U).$$(8)

It should be noted that the metallicity calculated here is the pre-depletion metallicity, that is the metallicity before any depletion onto dust grains occurs. The relation between the pre-depletion metallicity and the post-depletion metallicity depends on the depletion factors adopted by each model. For the JY20 model, this is given by $$12+\log {(\text{O/H})}_{\text{pre}}=12+\log {(\text{O/H})}_{\text{post}}+\mathrm{0.22.}$$(9)

The post-depletion metallicity is the current gas-phase metallicity in the ionized cloud. Thus, one should be cautious when comparing metallicities derived from photoionization models to those derived through other methods.

For each model, we performed three sets of fitting, the first one in the N2O2 versus O3O2 space, the second one in the N2–S2–R3 space, and the last one in the 5D space composed of N2O2, O3O2, N2, S2, and R3. When fitting the data using the N2O2 versus O3O2 method, we applied extinction corrections based on the Balmer decrements, assuming an intrinsic Balmer ratio F_{H α}∕F_Hβ = 2.86. The extinction curve we used is from Fitzpatrick (1999) with R_V = 3.1. We see whether these results agree with each other, and what is the main driver of different MI correlations.

3.1 Comparison of model predictions based on different emission-line ratios

In Fig. 4, we plotted all four SF-ionized models in the P₁ –P₂ diagram as well as their predicted MI correlations. This projection makes the model surfaces roughly edge-on and thus best reveals any discrepancy between the data and model surfaces in the 3D line-ratio space. The leftmost column compares the model surfaces with the H II regions (which are colored in green and yellow). All models were interpolated using the 2D interpolation function GRIDDATA in PYTHON. To better compare the locations of the model surfaces with the data distribution, we show only parts of the models that cover the middle 98% of the data along the P₃ axis, which is perpendicular to the P₁ –P₂ plane in 3D.

By eye, the JY20 model provides the best fit to the densest part of the data distribution in the P₁ –P₂ diagram. For all three combinations of line ratios, the JY20 model predicts positive correlations between 12 + log(O/H) and log(U). The 2D distributions of 12 + log(O/H) and log(U) derived by different sets of line ratios show good consistency overall. The only difference may be that the constraints based on the N2O2–O3O2 method gives a slightly smoother distribution than the other two methods.

In comparison, the L10 model underpredicts P₂ values significantly for P₁ ≲ 0. Since P₂ anticorrelates with metallicity (Ji & Yan 2020), the L10 model fails to describe the majority of the data at relatively low metallicities. As expected, for this model the fitting results from different methods show discrepancies at subsolar metallicities. When usingthe N2O2–O3O2 method, the L10 model predicts no correlation between 12 + log(O/H) and log(U) for spaxels withsubsolar metallicities. Whereas the N2–S2–R3 method gives a clear anticorrelation between 12 + log(O/H) and log(U). When all line ratios are combined, the resulting distributions of 12 + log(O/H) and log(U) closely resemble those predicted by N2, S2, and R3. This could be caused by the fact that N2, S2, and R3 overall have smaller uncertainties than N2O2 and O3O2. However, we note that the position of the model surface relative to the data distribution in 5D could also be relevant.

For the D13 model, it shows better consistency with the data in the N2–S2–R3 space, but it is still offset from the dense part of the reprojected data locus at low metallicities. Using this model, the N2O2–O3O2 method predicts a very broad distribution of log(U). The resulting 12 + log(O/H) and log(U) seem to have no correlation except at very high metallicities. On the other hand, the N2–S2–R3 method again predicts a negative correlation between the two quantities at subsolar metallicities and a positive correlation at supersolar metallicities. When all five line ratios are used, the overall shape of the relation between 12 + log(O/H) and log(U) looks similar to the one obtained using the N2–S2–R3 method, but with a larger scatter in log(U). In summary, the fitting results of the D13 model share similar features with those of the L10 model due to their similar behaviors in the N2–S2–R3 space. They all underestimate these line ratios at (their) subsolar metallicities and show better consistency with the data distribution at supersolar metallicities, which lead to two different correlations between 12 + log(O/H) and log(U) in these two regimes. In regions where these two model surfaces get close to the center of the data distribution in 3D, the derived 12 + log(O/H) and log(U) are clearly positively correlated.

Finally, the B17 model gives a more complicated result as it crosses the dense part of the H II region sample at P₂ ~ 0.2. It underestimates N2, S2, and R3 at very low metallicities, while it overestimates them at higher metallicities. The N2O2–O3O2 predictions show an anticorrelation between 12 + log(O/H) and log(U) until the metallicty reaches the solar value, after which these two quantities seem to start becoming positively correlated again. The predictions given by the N2–S2–R3 method exhibits a similar behavior, despite the turning point of the correlation lying at a subsolar metallicity rather than a solar metallicity. The 5D fitting result is again dominated by N2, S2, and R3.

We see that only the JY20 model predicts a consistent and positive MI correlation no matter which combination of line ratios is chosen. The predictions of other models all have significant dependencies on which a combination of line ratios was used, due to their obvious discrepancies with the data surface in the multidimensional line-ratio space. For these models, the predicted MI relations through the N2–S2–R3 method all have two different trends. That is, with increasing metallicities, the slopes of the relations are at first negative, and then become positive at some metallicities, which seems to be related to where the models start to get close enough to the dense part of the data in 3D.

Figure 5 shows the 1D distributions of 12 + log(O/H) and log(U) predicted by different models and different sets of line ratios. Only the results from the JY20 model show good agreement among different sets of line ratios. Judging from the histograms under the JY20 model, the N2–S2–R3 method seems to yield slightly more high metallicity spaxels, which implies that our linearly extrapolated stellar SED at the highest metallicity is not perfect. Another noticeable difference is that the peak of the log(U) distribution derived by the N2–S2–R3 method is higher than that derived by the N2O2–O3O2 method by ~ 0.1 dex. This difference originates from the fact that the model predicted O2 value is not entirely consistent with the observed O2 value (if we assume that the model predicted N2, S2, and R3 values are all correct). Since in our case O2 is only used in combination with N2 or O3, which have much larger wavelengths, it is also possible that the extinction correction is not correctly applied. We further discuss this issue in Sect. 5.1.

The other models all show bimodal 1D distributions in 12 + log(O/H), consistent with what we have seen in the 2D distributions. Unlike the N2–S2–R3 method, the N2O2–O3O2 method always gives single-peaked metallicity distributions. For log(U), the L10 model, D13 model, and B17 model all predict systematically lower values than those given by the JY20 model. Interestingly, only the B17 model predicts an obvious bimodal distribution in log(U) when using the N2, S2, and R3 lines or all five emission lines. Given the relatively flat stellar mass distribution in our sample, it is reasonable to assume that the true distribution in either 12 + log(O/H) and log(U) should be a single-peaked function; otherwise, one needs to explain why the chemical evolution model would favor a very different metal-enrichment timescale for galaxies with an intermediate metallicity. The double-peaked distributions we see can be qualitatively understood as the following. Even as the model surface significantly deviates from most of the spaxels in high-dimensional line-ratio space, the fitting algorithm is still looking for the closest part of the model surface to the data. Depending on the relative curvature of the model surface to the data distribution, data points that originally correspond to a wider distribution in 12 + log(O/H) (or log(U)) might be preferentially assigned to specific parts of the model grids that are closest to the data surface, causing the double-peaked distribution.

Figure 6 shows how the median absolute distances from the data points to the model surfaces in the N2–S2–R3 space changes as a function of the metallicities predicted by the models. The JY20 model surface lies closest to the majority of the data in 3D space, with the median distances d_3D ≲ 0.05 dex. The D13 model surface is equally close to the high-metallicity data, but fits the low-metallicity data worse. The L10 model surface lies similarly close to the data at 12 + log(O/H) ≈ 9.2, but it moves further away from the data as metallicity gets lower. Finally, the B17 model surface appears closest to the data at 12 + log(O/H) ≈ 8.4, and it starts to deviate from the data as metallicity increases. These results are consistent with what we found in Fig. 4 with the P₁ –P₂ diagram. This again shows that the agreement between the model surface and the central surface of the data distribution in multidimensional line-ratio space is an important metric for model evaluation.

Fig. 4 Different photoionization models viewed in the P1–P2 diagram and their predicted relations between the metallicity and ionization parameter. Leftmost column: Different SF-ionized models and MaNGA MPL-7 data in the P1 –P2 diagram. The models were interpolated and cut so that only the parts that cover the middle 98% of the data along the hidden P3 axis are shown. The spaxels with contributions from SF ionization greater than 90% are colored from yellow to green, while the rest of the spaxels are colored from white to black. Only the former were used to derive the metallicities and ionization parameters. Right columns (a) to (c): derived metallicities and ionization parameters using different SF-ionized models and different combinations of line ratios. In each panel, the green arrow marks the solar oxygen abundance adopted by the corresponding model.

Fig. 5 Histograms showing the number distributions of 12 + log(O/H) and log(U) predicted by different SF-ionized models and different combinations of line ratios.

Fig. 6 Median absolute Euclidean distances in the N2–S2–R3 3D space as a function of fitted metallicities for different SF-ionized models. The metallicities shown were derived based on different models using the N2–S2–R3 method, thus having different ranges. The shaded regions encompass the middle 68% of the data.

3.2 Comparison of model parameters

By comparing four different SF-ionized models, we see that the model which fits the central surface of data in the multidimensional line-ratio space gives the most consistent predictions on the metallicity and ionization parameter. The location of the model surface in the line-ratio space is determined by parameters other than metallicity and ionization parameter (i.e., the secondary parameters). Variations in different secondary parameters move the model surface in different directions by different amounts, as detailed in Ji & Yan (2020). In this section we attempt to find the most relevant secondary parameters that distinguish the JY20 model from other models.

However, before we discuss the effects of the secondary parameters, we note that there are two different photoionization codes used in these models, that is to say CLOUDY and MAPPINGS. According to D’Agostino et al. (2019), the two codes result in almost identical results if the input parameters are the same. However, as pointed out by Law et al. (2021a) and Belfiore et al. (2022), earlier photoionization codes do not have a reliable estimation on the dielectronic recombination rate of S⁺⁺ and use the charge-normalized mean dielectronic recombination rates of C, N, and O as an approximation (Ali et al. 1991). The current version of CLOUDY (v17, Ferland et al. 2017) improves the estimation of this value (Badnell et al. 2015), which results in stronger [S II] lines. Among the four models, only the JY20 model uses the latest photoionization code. Therefore, we should keep in mind that the S2 predicted by other models are underestimated by some amounts, but they could be offset by the effects of other input parameters.

Table 1 lists the input parameters of the four models. We start with the effect of the density structure. The hydrogen densities of the models all lie in the range 10 cm⁻³ < n_H ≤ 100 cm⁻³. In fact, for radiation-bounded clouds with n_H ≲ 100 cm⁻³, the model-predicted line ratios are almost identical, except for some small differences at very high metallicities (Ji & Yan 2020). We also tested the effects of different equations of state on the simulated cloud, and the result shows that a cloud with a constant density produces very similar results as those from an isobaric cloud with the same initial density.

The geometry of the cloud is a more complicated issue. The spherical geometry is often assumed in the dynamical modeling of H II regions, while in the real world it is difficult to have perfect symmetry. Density inhomogeneities in the molecular clouds cause asymmetries and the feedback from the central star clusters could create blister H II regions.In such cases, a plane-parallel geometry might be a more realistic choice (e.g., the Orion nebular, Wen & O’dell 1995). We computed models with these two different geometries while keeping other parameters the same. For the spherical model, we varied the inner radius from 3 to 10 parsecs. We found that these model surfaces cover a similar area in the line-ratio space (regardless of whether the geometry is closed or open for the spherical case) except at very highmetallicities, suggesting that the geometry is not likely to be the dominant factor that causes the differences we see among the four models. For an ionized cloud with a high ionization parameter at the illuminated face, the H II region becomes thick. If the inner radius of the H II region is comparable to the thickness of the ionized layer, then the geometric dilution of the ionizing photons becomes important and the resulting emission line spectrum would be different from that produced by a plane-parallel model with the same ionization parameter. This could explain the difference we see in the models with different inner radii. In Sect. 4, we discuss whether this effect would become significant if we were to change the inner radius according to dynamical models.

The next input parameter in Table 1 is the stellar SED, which determines the relative number of ionizing photons andit is quite important in setting the line ratios of H II regions. The first three models we compare all use stellar SEDs generated by STARBURST99. Both the JY20 model and the D13 model assume a continuous SFH of 4 Myr. Different from the JY20 model, the D13 model assumes a Salpeter IMF (Salpeter 1955). However, we tested models with different IMFs and found that the IMF effect is negligible for the line ratios considered here. The remaining difference is that the D13 model uses an earlier STARBURST99 SED adopted by Dopita et al. (2000), which is harder compared to the newer SEDs generated by later versions of STARBURST99⁶. A harder input SED would strengthen the low-ionization lines such as [N II] and [S II] by increasing the relative size of the partially ionized zone inside the ionized cloud. However, compared to the JY20 model, the D13 model instead shows comparable N2 ratios, but lower S2 ratios. Thus, the effect of the earlier STARBURST99 SED is either not large enough or is offset by the effects of some other parameters (e.g., the underestimated dielectronic recombination rate of S⁺⁺). Meanwhile, the L10 model assumes a continuous SFH of 6 Myr, which should give a very similar SED as that of the JY20 model (see Fig. 7 of Ji & Yan 2020). Unlike the first three models, the B17 model adopts a much harder SED generated by FSPS (Conroy et al. 2009), which contributes to its overpredictions of N2, S2, R3, and equivalently P₁. Such an SED might allow one to describe H II regions undergoing intense star formation, but it is too hard for the majority of the H II regions in our sample. In Fig. 7, we plotted a model with the B17 SED set, but kept the other parameters the same as those in the JY20 model. This model shows relatively large P₁ values as expected.

Finally, the chemical abundance is a fundamental ingredient for setting the line ratios. The eventual chemical abundance in the cloud is determined by the adopted solar abundance, the N/O ratio, and dust depletion factors. The JY20 model and the D13 model both use the solar abundance set in Grevesse et al. (2010), while the L10 model and the B17 model use the solar abundance set in Anders & Grevesse (1989). The abundances of relevant elements, such as nitrogen, oxygen, and sulfur, are more abundant inthe latter set. However, the abundances are further altered by the N/O ratio and dust depletion factors, as we detail below.

The nitrogen is mainly created by the CNO cycles in stars. Depending on whether the carbon and oxygen that get converted to nitrogen are newly created or pre-existed in the clouds that formed the stars, the resulting nitrogen is called to have a primary origin or a secondary origin. The amount of secondary nitrogen is very sensitive to metallicity, making nitrogen-related lines good metallicity tracers (Alloin et al. 1979). Theoretical metallicity calibrations that use nitrogen lines thus strongly rely on the assumed relation between the N/O ratio and metallicity (which we call the nitrogen prescription hereafter). Higher N/O ratios result in higher N2, and slightly lower S2 and R3 due to the thermal balance. Ji & Yan (2020) found that the nitrogen prescription provided by Dopita et al. (2013) fits the MaNGA H II regions best and they used this for their SF-ionized model. Therefore, the nitrogen prescription does not contribute to the difference between the JY20 model and the D13 model. It is unclear what kind of nitrogen prescription is used in the L10 model. Judging from the shape of the L10 model in the P₁ –P₂ diagram, it could have used a nitrogen prescription that has too much primary nitrogen, causing overestimated N2 or equivalently underestimated P₂ at lower metallicities. The B17 model adopts the nitrogen prescription given by Dopita et al. (2000), which yieldsa higher N/O at a given metallicity, as is shown in Fig. 8. Interestingly, due to the effectof dust depletion and higher solar abundances, the actual N/O ratio (at a given O/H) measured in the gas phase is similar to that in the JY20 model. However, if the relative metallicity (e.g., [O/H] or Z/Z_⊙) is used when comparing different models, one should be cautious about the difference in the adopted solar abundances, which we discuss in the next paragraph.

The dust depletion factors modify the gas-phase abundance in H II regions by removing part of the elements that condensed into dust grains. The depleted elements have two effects on the line ratios. First, the intensities of lines emitted by elements that are significantly depleted decrease due to the lower abundances. Second, the removed coolants influence the thermal balance of the cloud and raise the equilibrium temperature, thus strengthening lines emitted by elements that are not significantly depleted. Meanwhile, the increasing amount of dust also heats the cloud up through the photoelectric effect. The JY20 model uses the default depletion set in CLOUDY, which depletes oxygen by 0.22 dex and does not deplete nitrogen at all (Cowie & Songaila 1986; Jenkins 1987). In comparison, the D13 model depletes oxygen less by 0.15 dex and depletes nitrogen more by 0.05 dex. Although the difference in oxygen depletion is substantial, it is not the main cause of the apparent shift of the D13 model surface from the data at low metallicities. For comparison, we computed a model using the JY20 parameters but with the D13 depletion set. This model deviates from the JY20 model in the P₁ –P₂ diagram in a different way, as shown by Fig. 7. The largest separation between this model and the JY20 model is roughly 0.1 ~ 0.2 dex along P₁ (with the JY20 model showing overall larger P₁ values) and it is negligible along P₂, which occurs in high metallicities rather than low metallicities.

For the B17 model, it uses the depletion set provided by Dopita et al. (2000). It depletes oxygen by 0.22 dex, just as JY20. However, it depletes 0.22 dex more nitrogen compared to the JY20 model. Its solar nitrogen abundance is 0.22 dex higher, making the gas-phase N/H the same between the two models at fixed O/H. However, B17’s adopted solar oxygen abundance is 0.24 dex higher and its solar sulfur abundance is 0.09 dex higher. There are two ways to interpret the consequence from these differences: (1) At the same relative pre-depletion metallicity (relative to solar), [O/H], or stellar metallicity, Z/Z_⊙ (Z_⊙ ≈ 0.020 for both STARBURST99 and FSPS), the absolute post-depletion O/H adopted by the B17 model is higher than that of the JY20 model, which also makes the corresponding post-depletion N/O larger according to the right panel of Fig. 8. This makes the predicted R3 and N2 values larger. (2) Alternatively, we could understand the difference as a shifted correspondence between the stellar SEDs and the absolute post-depletion metallicities. At the same post-depletion O/H, the JY20 model and the B17 model show similar post-depletion N/O as shown in the left panel of Fig. 8. However, since the solar O/H is higher in the B17 model, the corresponding [O/H] and Z/Z_⊙ are lower. Since lower metallicity SEDs are in general harder, this makes the B17 model predict larger line ratios. We computed another model using the JY20 parameters, but with the B17 abundances and depletion set. The resulting model indeed predicts larger line ratios and it is shifted to the right in the P₁ –P₂ diagram, as shown in Fig. 7.

We note that current photoionization code does not treat the dust depletion in a totally self-consistent way. The total mass of elements depleted does not necessarily match the total mass of dust assumed, and they do not match in composition either (see e.g., Snow & Witt 1996). In addition, the depletion factors could vary in different locations inside galaxies and could depend on the environment, as noted by Jenkins (2009). It also remains an open question whether or not sulfur depletes in the ISM (Sofia et al. 1994; Jenkins 2009; White & Sofia 2011; Amayo et al. 2021). For a more general discussion of the effect of dust depletion on predicted line ratios, readers can refer to Gunasekera et al. (2022).

In summary, there are three major factors that contribute to the differences among the models. First of all, the change in the dielectronic recombination rate of S⁺⁺ in the updated version of CLOUDY makes the JY20 model predicts larger S2, resulting in some of the differences between this model and the other three models. In addition, the harder stellar SED adopted by the B17 model partly offsets the underestimation of S2, but it overestimates N2 and R3. Last but not least, the gas-phase chemical abundance plays an important role. At a given [O/H] or Z/Z_⊙, the B17 model has overall higher post-depletion abundances due to a combined effect of their adopted nitrogen prescription, dust depletion set, and solar abundance set, which also contributes to its overestimation of the line ratios. All these factors combined make the JY20 model a better representation of an average H II region in our sample. A caveat here is that there is still degeneracy between some of the secondary parameters. For example, a harder stellar SED might be able to cancel out the effect of the underestimation of elemental abundances or overestimation of depletion factors. Independent observational constraints are needed to determine the best values of these parameters. It is possible to include more emission lines and break the degeneracy into even higher dimensions, which is beyond the scope of the current paper. Nevertheless, we have seen that by tuning the model parameters within a reasonable range, we can obtain a best-fit model that lies very close to the majority of H II regions in a multidimensional line-ratio space. The true model and its predictions should not be too far from what we have presented here, as long as the adopted secondary model parameters are not significantly biased.

The best-fit model we show predicts a positive MI correlation extending from subsolar to supersolar metallicities. In the next section, we discuss the physical interpretations of this correlation.

Fig. 7 Comparisons of photoionization models with different input parameters in the P1 –P2 diagram. The cyan model is our fiducial model. The red model has the fiducial abundances and depletion set, but uses the SED adopted by the B17 model. The purple model uses the fiducial SED, but adopts the abundances and depletion set of the B17 model. The blue model uses the fiducial SED, but adopts the abundances and depletion set of the D13 model. The modelsare interpolated and cut so that only the parts that cover the middle 98% of the data along the hidden P3 axis are shown. The density distribution of the H II region sample is plotted in green and yellow, while the rest of the MaNGA sample is plotted in black and white.

Fig. 8 Comparisons of pre- and post-depletion nitrogen prescriptions for the BY17 and the JY20 models. Both the input prescription and the output prescription modified by dust depletion are shown. The JY20 model adopts thesame pre-depletion nitrogen prescription as Dopita et al. (2013) and depletes the elements using the CLOUDY default depletion set. The B17 model follows Dopita et al. (2000) for both the pre-depletion nitrogen prescription and the depletion pattern. While the left panel shows the N/O ratios as functions of the absolute oxygen abundance 12 + log(O/H), the right panel shows the N/O ratios as functions of the relative oxygen abundance [O/H], normalized by the solar O/H adopted by each model. The solar oxygen abundances, 12 + log(O/H)⊙, adopted by Dopita et al. (2000) and Dopita et al. (2013) are 8.93 and 8.69 and indicated by the orange triangle and the cyan star in the plot, respectively.

Fig. 9 Correlation between the ionization parameter and Hα surface brightness (left), and the mass-metallicity relation color coded by the Hα surface brightness (middle) and ionization parameter (right). Left panel: 2D histogram of the logarithmic ionization parameter and logarithmic Hα surface brightness. The density distribution shows the results obtained by the JY20 model using all five line ratios. Thedashed red contours and the dotted-dashed black contours show the distributions obtained by using N2O2–O3O2 ratios and N2–S2–R3 ratios, respectively. The contour levels are in logarithmic scales and range from the 16th percentile to the 84th percentile. Middle panel: spatially resolved mass-metallicity relation of the sample spaxels color coded by themedian Hα surface brightness inside each bin. The black contours represent the density distribution of spaxels obtained by using all five line ratios. Right panel: same as the middle panel, but color coded by the median ionization parameter inside each bin.

4 Physical interpretations

The previous section has shown that according to photoionization models, the metallicity and ionization parameter are correlated. The first theoretical analysis on how these parameters are correlated is given by D06. By studying a dynamical model of a wind-driven bubble that expands over time, D06 found $q\propto {(Z/Z_{\odot })}^{-0.8}$, where q = Uc and Z is the stellar metallicity. Contrary to this theoretical prediction, our best-fit photoionization model results in a positive MI correlation (assuming the stellar metallicity directly scales with the gas-phase metallicity). This positive correlation has also been noted by many authors (e.g., Dopita et al. 2014; Poetrodjojo et al. 2018; Mingozzi et al. 2020). Specifically, D14 studied spatially resolved data of ten luminous infrared galaxies (LIRGs) and found strong MI and SFR-I correlations. Hence it is possible that the ionization parameter is related to metallicity through their common dependency on the SFR. Although the dynamical model of D06 predicts a negative correlation between theionization parameter and metallicity, the regulation from star formation could change the statistical behavior in the data. It is also possible that the dynamical model proposed by D06 is simply not realistic enough to describe most of the observed H II regions, or our photoionization model is oversimplified and introduces a bias.

In this section, we explore the potential mechanisms that lead to the positive MI correlation, including the common dependency of the metallicity and ionization parameter on a third parameter (e.g., the SFR and related parameters), the incompleteness of the wind-driven bubble model, and the oversimplification of the photoionization model.

4.1 Impact of star formation

Star formation plays an important role in regulating the chemical enrichment of galaxies. At fixed stellar masses, galaxies with higher SFRs tend to have lower metallicities, leading to the fundamental relation between stellar mass, gas-phase metallicity, and SFR (Ellison et al. 2008; Mannucci et al. 2010). Whereas overall, the metallicity positively correlates with the SFR since more massive SF galaxies tend to have both higher SFRs and higher metallicities. In contrast, the dependence of the ionization parameter on the SFR is less clear since it is not a direct observable.

D14 explored several possibilities to relate the ionization parameter to the SFR. One plausible scenario is the effect of cluster mass. D06 observed in their dynamical model that the ionization parameter weakly depends on the cluster mass $$U\propto M_{\text{cl}}^{1/5}.$$(10)

According to D06, this relation comes from the dependency of U on the relative number density of the H II region to the ambient medium (see their Eq. (12)). There is evidence that galaxies undergoing more intense star formation can hold a larger number of massive clusters (e.g., Bastian et al. 2008; Powell et al. 2013), thus connecting the ionization parameter to the SFR. However, D14 also found results that contradict observations using this relation. Their data show $U\propto {\Sigma }_{\text{SFR}}^{0.34}$. Assuming there are no other dependencies, we have $M_{\text{cl}}\propto {\Sigma }_{\text{SFR}}^{1.7}$. Combining the relation Σ_SFR ∝ Σ_clM_cl (where Σ_SFR is the SFR surface density and Σ_cl is the surface number density of young clusters) and the star formation law of Kennicutt (1998) ${\Sigma }_{\text{SFR}}\propto {\Sigma }_{\text{g}}^{1.4}$ (where Σ_g is the gas surface density), one gets ${\Sigma }_{\text{cl}}\propto {\Sigma }_{\text{g}}^{-0.98}$. D14 then argued that observationally speaking, the density of clusters actually increases with increasing gas density, and thus relation (10) is not a viable solution.

If the ionization parameter is truly regulated by the SFR through $U\propto {\Sigma }_{\text{SFR}}^{\alpha }$ with Σ_cl being the medium, we must have the power-law index α < 0.2 in order to make Σ_cl positively correlate with Σ_g. The left panel of Fig. 9 shows the spatially resolved SFR-I relation of our sample derived using the J20 model. We use the extinction-corrected logΣ_Hα to represent logΣ_SFR as they differ by only a constant (Kennicutt & Evans 2012). We can see a weak positive correlation between the two quantities. The correlation is weakest when derived using the N2O2–O3O2 method. If we fit a linear relation to the data, the slope is below 0.2, with the N2–S2–R3 method giving the largest value and the N2O2–O3O2 method giving the smallest value. There is, however, large uncertainty in this relation due to the large intrinsic scatter in log(U). It also appears that the relation is not linear, but it flattens at high log Σ_Hα. The Pearson correlation coefficients for the (logarithmic) MI correlation, SFR-I correlation, and SFR-M correlation are 0.70, 0.61, and 0.42, respectively. Therefore, despite the derived slope falling into the plausible range, it is unlikely that the weaker SFR-I and SFR-M correlations are able to give rise to a stronger correlation between the metallicity and ionization parameter, with the SFR as the main driver. To produce a stronger MI correlation, the residuals in SFR-I and the residuals in SFR-M have to be correlated, which would mean whatever factor that determines the residual is a more fundamental parameter than the SFR to setup the MI correlation. This proves that the SFR cannot be the dominant factor for establishing the correlation.

D14 also proposed another potential mechanism that modifies the ionization parameter through star formation. The geometry of a realisticH II region could be highly nonspherical. If the fragmented molecular clouds within an H II region can survive long enough and move close to the central OB stars due to turbulent motions, they would become highly ionized and raise the overall ionization parameter we observe. D14 found that for SF regions with lower cluster masses or higher pressures, it is possible to have a molecular cloud to cross an H II region within the lifetime of the OB stars. Specifically, H II regions older than 1 Myr with M_cl ≲ 10⁴ M_⊙ and n_H ≳ 10² cm³ can have a shorter crossing time compared to the expansion time. However, in MaNGA we found that most of the SF regions have hydrogen densities close to 14 cm⁻³ (Ji et al. 2020), which means the gas pressures are not high enough for the majority of our sample. We also found no correlation between the derived ionization parameter and the density (and pressure) sensitive ratio [S II]λ6716/[S II]λ6731 in our sample. Unfortunately, MaNGA does not have enough resolution power for investigations of turbulent motions on subkiloparsec scales. Thus, we cannot conclude observationally on the dependency of the turbulent motion on the SFR. The theoretical calculation of Joung et al. (2009) does predict an increase in the gas velocity dispersion due to supernova explosions. Regardless, given the weakness of the SFR-I correlation, this solution is questionable.

A related question is the dependence of the various relations we consider here on stellar masses. Since we see the ionization parameter positively correlates with both the metallicity and SFR, whether it contradicts the fundamental relation between the stellar mass, gas-phase metallicity, and SFR is worth investigating. In the middle and right panel of Fig. 9, we plotted the spatially resolved mass-metallicity relations color coded by log Σ_Hα and log(U), respectively. The stellar masses are drawn from the MaNGA PIPE3D Value Added Catalogue (VAC) (Sánchez et al. 2016). The fundamental relation does not seem apparent in the local scale, which has also been noted by several works (e.g., Sánchez et al. 2013, 2017; Barrera-Ballesteros et al. 2017). At a fixed stellar mass surface density, increasing the SFR surface density does not significantly lower the gas-phase metallicity on average; however, there could still be a small effect as can be seen from the color gradient. On the other hand, Σ_SFR appears strongly correlated with ${\Sigma }_{M_{*}}$. Meanwhile, the right panel of Fig. 9 shows that at a fixed stellar mass surface density, there is still a strong positive MI correlation. We further checked the median MI relations at different stellar mass surface density bins and Hα surface brightness bins in Fig. 10. One can see that the median MI relations are nearly identical in different ${\Sigma }_{M_{*}}$ bins. This indicates that the stellar mass surface density does not drive the MI correlation.

If we investigate the MI relations in different Σ_Hα bins instead, we see the MI relations persist in most bins. However, as Σ_{H α} increases, the relation becomes flattened and elevated, and nearly vanishes at the highest Σ_{H α}. This is consistent with what we see in the left panel of Fig. 9: there is a large range of possible U at low Σ_{H α} values, but a much narrower range of U at high Σ_{H α} values. The ionization parameter saturates at high Σ_SFR and becomes independent of metallicity. It is possible that at very high U, the radiation pressure becomes important and prevents U from gettingeven larger (Yeh & Matzner 2012). Another potential explanation is related to the diffuse ionized gas (DIG) surrounding the H II regions. Leaking radiation from the H II regions could ionize the DIG, which exhibits low ionization states as a result of the diluted ionizing flux. If we limit the sample to only the low SFR regions, where the DIG contribution to the emission-line spectra becomes important, we see that the MI correlation also becomes stronger. To explain the trend, however, it requires the molecular clouds with low metallicities to be more leaky compared to the high metallicity clouds, which has not been observed to our best knowledge. Finally, the measurement bias induced by the photoionization models might contribute to this effect as well, which we detail later in Sect. 4.3.

In summary, there indeed exists a positive SFR-I correlation in our sample. However, the coupling between the metallicity and ionization parameter appears stronger than that between the metallicity and SFR, or between the ionization parameter and SFR. In fact, the dependence on the SFR weakens the overall MI correlation, as shown in the right panel of Fig. 10. It is thus questionable whether the main driver of the MI relation is the SFR, while the SFR-I correlation could come from the influence of the cluster mass or the turbulent motions of the molecular clouds inside H II regions. This is related to the details of the dynamical models for H II regions, which is the topic of the next subsection.

Fig. 10 Dependence of the log(U) versus 12 + log(O/H) relation on the stellar mass surface density (left) and Hα surface brightness (right). The ionization parameters and metallicities shown were derived using all five line ratios. The contours represent five density levels equally spaced in the logarithmic space from the 16th percentile to the 84th percentile. The dotted-dashed lines indicate the median relations in different bins.

4.2 Dynamical evolution of HII regions

The structure of an H II region is shaped by the feedback from the central young massive stars. The time evolution of the structure of the H II region is described by the dynamical models, which in principle set the initial conditions for photoionization models. In this section, we discuss dynamical models with different assumptions and whether they lead to any correlation between the metallicity and ionization parameter.

Among the various feedback mechanisms considered in the modeling, stellar winds and photoionization are two main factors for shaping the geometry. Without stellar winds and any other source of mechanical energy, the radial structure of the ionized region is solely determined by photoionization, with the simplest case being the two-phase solution given by Spitzer (1978). This solution describes a central star photoionizing a homogeneous cloud purely composed of hydrogen. During the first phase, the ionization front (IF) moves exponentially with a timescale of τ = 1∕n_eα_B, where n_e is the electron density and α_B is the recombination coefficient under the on-the-spot approximation. Once the IF reaches the Strömgren radius given by $$r_{\text{S}}={\left(\frac{3Q_0}{4\pi n_{\text{e}}^2{\alpha }_{\text{B}}}\right)}^{1/3},$$(11)

the ionization rate equals the recombination rate and the second phase of expansion begins. A shock front is created before the IF and the compressed gas is moving toward the neutral ISM. The location of the IF is given by $$r_{\text{i}}=r_{\text{S}}{\left[1+\frac{7c_{\text{i}}(t-t_{\text{S}})}{4r_{\text{S}}}\right]}^{7/4},$$(12)

where c_i is the sound speed in the ionized medium (typically ~ 10 km s⁻¹) and t_S is the time when the IF reaches r_S. This solution is oversimplified and does not take important factors into account, such as metals and dust. Meanwhile, we need to define the ionization parameter. Since it would be a function of radius, we need to specify a representative location or define an average to calculate it. One commonly adopted choice is the imagined ionization parameter at the Strömgren radius, which can be written as $$U_{\text{S}}=\frac{Q_0}{4\pi r_{\text{S}}^2{n}_{\text{H}}c}\approx {\left(\frac{Q_0{n}_{\text{e}}{\alpha }_{\text{B}}^2}{36\pi c^3}\right)}^{1/3}.$$(13)

Clearly, we have U_S = U_VA∕3, where U_VA is the volume-average ionization parameter inside the Strömgren sphere. We note that this expression does not describe the actual ionization parameter measured at the Strömgren radius, as the number of ionizing photons decreases with increasing depth into the cloud.

Given that ${\alpha }_{\text{B}}\approx 2.56\times {10}^{-13}{T}_4^{-0.83} {\text{cm}}^3{\text{s}}^{-1}$ (Draine 2011), we have $$U_{\text{S}}\propto Q_0^{1/3}{n}_{\text{e}}^{1/3}{T}^{-0.55}.$$(14)

If we now consider adding metals to the simple model, T drops due to cooling. In the meantime, dust abundances increase with increasing metal abundances (Draine 2011), which would heat the gas up due to photoelectric heating. Despite this competing effect, the net result of scaling the metallicity and dust abundances up is a drop in T, which causes an increase in U_S defined above. However, we note that the actual radius within which the recombination happens would be smaller than r_S given in Eq. (11) since metals and dust grains also absorb part of the ionizing photons. Therefore, a better way to calculate U_S when metals and dust grains are included is to simply define r_S at the location where half of the hydrogen is ionized (i.e., n(H+)∕n(H) ~ 0.5). The whole H II region would become more compact, which results in a larger U_S, as shown by Haworth et al. (2015). On the other hand, according to the STARBURST99 models, Q₀ provided by the ionizing stars drops as metallicity increases, provided that all models are normalized to the same SFR.

We can check these dependencies with our measurements. Fig. 10 shows that log(U) roughly increases by 0.2 dex as [O/H] changes from 0 to 0.3 (solar to double-solar value). We generated a series of photoionization models corresponding to spherical H II regions with no inner cavity and with Q₀ ranging from 10^47.5 s⁻¹ to 10⁴⁹   s⁻¹. The other conditions were set to be consistent with our fiducial model. For the model with the largest Q₀, the change in log(U_S) with increasing metallicity is comparable to what we see in the observed data⁷. However, there are two difficulties in explaining the observed MI correlation with this scenario. First, at each metallicity, U_S returned by the model is nearly 6 times the median U in the observed data. Lowering Q₀ can lower U_S, but it wouldflatten the MI relation significantly. Also, the ionizing luminosity of the central stars should decrease with increasing stellar metallicity, which further flattens the MI relation predicted by the model. Second, observations of nearby H II regions show that the ionizing stars are separated from the bulk of the ionized clouds by hot diffuse gas (e.g., Pellegrini et al. 2007, 2011; Güdel et al. 2008). This is because stellar wind feedback from young massive stars creates shocked-wind bubbles within their birth clouds. If we include non-negligible inner cavities (with radii comparable to r_S) in H II region models, the MI relation becomes much flatter. The strength of the stellar wind feedback, however, is also a function of metallicity. This adds more complexity to the dynamical modeling of H II regions, as we detail in the following.

Once the stellar winds are included, the structure of the ISM around the central star cluster changes substantially. Figure 11 shows the structure of a spherical wind-blown cloud. A wind-driven bubble is created and gradually sweeps the outer H II region. Due to their high speed (typically 1500 ~ 2500 km s⁻¹ for O-type stars), stellar winds drive shocks into the H II region and produce a region with a very high temperature and low density (region b). For the dynamical models containing stellar winds, the geometry of the H II region is set by the inner (shock) radius, r_in, and the Strömgren radius. When the stellar winds are important, the relative thickness of the H II region is determined by the ionization parameter at r_in and is generally small if U < 10^−2.5~−2 (Dopita et al. 2006, but the exact value depends on assumptions of the dynamical model). If the shocked-wind radius is small and/or the ionization parameter is high, the H II region becomes thick. In this situation, the geometric dilution of the ionization photons would become important. Measuring the ionization parameter of such an H II region with a photoionization model that has a mismatched geometry would introduce a bias. We discuss a thin H II region first, and investigate the geometric bias for measuring athick H II region in the following section.

Weaver et al. (1977) gave an analytical solution for such a wind-driven bubble. In their model, the H II region is a thin isobaric shell. The position of the shell is given by the adiabatic solution $$r_{\text{shell}}={\left(\frac{250}{308\pi }\right)}^{1/5}{L}_{\text{w}}^{1/5}{\rho }^{-1/5}{t}^{3/5}\approx r_{\text{in}},$$(15)

where $L_{\text{w}}=\frac{1}{2}\dot{M}{v}_{\text{w}}^2$ is the mechanical power of the stellar winds and ρ is the density of the ambient medium, that is the H II region beyond the shock front. The pressure of the H II region is $$P=\frac{7}{{(3850\pi )}^{2/5}}{L}_{\text{w}}^{2/5}{\rho }^{3/5}{t}^{-4/5}.$$(16)

Eliminating t and using ρ ∝ n, we have the following: $$P\propto {\left(\frac{L_{\text{w}}}{n{r}_{\text{in}}^2}\right)}^{2/3}n={\left[4\pi cU\left(\frac{L_{\text{w}}}{Q_0}\right)\right]}^{2/3}n,$$(17)

where we use n to represent the hydrogen density of the H II region. The shell would keep expanding following Eq. (15) until its internal pressure equals the ambient pressure (i.e., the stall condition). After that, the shell moves in a momentum-conserving way and is finally destroyed by turbulence (Oey & Clarke 1997; Dopita et al. 2005). Combining Eq. (17) with P = nkT, we have $$U\propto \left(\frac{Q_0}{L_{\text{w}}}\right)T^{3/2}.$$(18)

This equation is a function of time. According to this equation, increasing the gas-phase metallicity would lower T, thus decreasing U. Meanwhile, increasing the stellar metallicity would generally lower the Q₀∕L_w ratio as the stellar winds and the stellar atmosphere become more opaque. The combined effect is a decrease in U, which gives an anti-correlation between the ionization parameter and metallicity. For a single massive star, the mass loss rate Ṁ roughly scales with the stellar metallicity following Ṁ ∝ Z^0.5~0.85 (Puls et al. 2008). Based on the STARBURST99 models, D06 estimated the anti-correlation to roughly follow U ∝ Z^−0.8. It is noteworthy that the mechanical energy from the stellar winds have to be manually scaled down by roughly 1 dex to reproduce the observed range of U (Dopita et al. 2005). The same discrepancy has been noted by Nazé et al. (2001) and Harper-Clark & Murray (2009) when comparing the observed sizes and expanding velocities of the interstellar bubbles with theoretical predictions. Hence it is likely that there are mechanisms in reality that lower the efficiency of the power output from stars, which could be related to the fact that the stellar winds are clumpy (Evans et al. 2004; Bouret et al. 2005; Fullerton et al. 2006; Puls et al. 2006, 2008). In addition to the inefficiency of stellar winds, Harper-Clark & Murray (2009) suggested that the leakage of the hot shocked gas from H II regions could also lower the radius and increase the ionization parameter by lowering the internal pressure. Whether and how these hyper parameters depend on the metallicity is unknown and requires further investigation.

The above adiabatic solution is based on the assumption that the shocked stellar winds do not have enough time to cool. However, as shown by Mac Low & McCray (1988), the cooling time of the wind bubble is typically short in dense molecular clouds, with t_cool ≲ 10⁴ yr. In such cases, the pressure directly comes from the momentum of the winds, $$P=\frac{{\dot{p}}_{\text{w}}}{4\pi r_{\text{in}}^2}=\frac{{(2\dot{M}{L}_{\text{w}})}^{1/2}}{4\pi r_{\text{in}}^2}.$$(19)

Where ṗ_w is the force imposed by the wind momentum and Ṁ is the mass loss rate of the central star. Equating the above equation with P = nkT, we obtained $$U=\frac{k}{2^{1/2}c}\frac{Q_0}{{(\dot{M}{L}_{\text{w}})}^{1/2}}T.$$(20)

If the metallicity of the cloud is increased, T and Q₀ would decrease, while ${(\dot{M}{L}_{\text{w}})}^{1/2}$ would increase. Therefore, the model still predicts an anti-correlation between the ionization parameter and metallicity, despite being in a different form.

Thus far, we have discussed dynamical models with and without stellar winds. A related question is whether the stellar winds are important in H II regions in general. Geen et al. (2020) studied this problem and conclude that dynamically speaking, stellar winds do not play an important role in most H II regions. However, stellar winds do shape the geometry of the H II regions even when they are dynamically unimportant. In other words, a considerably large shocked wind bubble can exist when the location of the IF is barely affected by the presence of stellar winds.

We have seen that none of these dynamical models are able to explain the positive correlation we found between U and Z. The complicating factor in this problem is that U is not a direct observable and it has to rely on certain assumptions of the geometry of the H II regions. In reality, the geometries of H II regions are complicated and can be far from spherical, and U also varies with locations inside H II regions. It is also unclear how stellar winds shape the shocked bubble when the ambient medium is highly asymmetric and inhomogeneous. The mysterious inefficiency of stellar winds in the wind-driven bubble model could be the key to resolving the problem, which we intend to investigate in future work. Quite surprisingly, despite the complicated picture shown here, the resulting correlation between U and Z is relatively clear and strong. This raises the question about the robustness and reliability of our derivations of these quantities based on photoionization models. In the following section, we discuss this point in detail.

Fig. 11 Schematic plot of a wind-driven bubble model. Region a is the free-wind region filled with hypersonic stellar winds. Region b is the shocked-wind region, which sets the inner boundary for the H II region. Region c is composed of shocked gas. Region d is the H II region. Depending on the ionization rate, region d can be trapped within region c or extend beyond it. Region e is the ambient neutral medium. The relative sizes of the regions are not to scale.

4.3 Geometric bias

When measuring the ionization parameter of a given H II region, one needs to define at which location the ionization parameter is measured. In addition to choosing the imagined U at the Strömgren radius, choosing the U at the inner radius is common and easy to carry out in photoionization modeling. Thus, this definition has been adopted by many photoionization models (e.g., Levesque et al. 2010; Dopita et al. 2013; Byler et al. 2017). However, the definition of the ionization parameter inevitably introduces an uncertainty associated with the actual geometry of the H II region, or more precisely, the relative thickness of the H II region (i.e., r_in ∕r_i).

For a spherical H II region, the ionizing flux decreases more rapidly with radius due to the geometric dilution. Hence, the spherical cloud would have a thinner ionized layer compared to a plane-parallel one with the same inner U. Their ionization structures are also different. Roughly speaking, the average location at which we measure the average line ratios would exhibit a lower U for a spherical model since the ionizing flux dilutes more. In other words, suppose we measure a spherical H II region with a photoionization model with a plane-parallel geometry, the measured ionization parameter would be lower compared to the ionization parameter defined at the inner radius. Similarly, if we use a spherical model with a larger inner radius to measure a spherical H II region with a smaller radius, the resulting U would also be underestimated.

The amount of the bias depends on the relative thickness of the H II region. If r_in ∕r_i is small, the H II region is thick and the effect of geometric dilution becomes important. Both the ionization parameter and metallicity affect the thickness of the H II region. For a higher ionization parameter, there is a larger number of ionizing photons and thus the ionized layer becomes thicker. A higher metallicity, on the other hand, would have two effects. First of all, r_i deceases due to both a lower equilibrium temperature and a reduction in the number of ionizing photons, which reduces the thickness of the H II region. In addition, according to the wind-driven bubble model by Weaver et al. (1977), the inner radius would increase as a result of stronger stellar winds (see Eq. (15)). If the inner U is held fixed, there would be less geometric dilution as the H II region becomes more plane-parallel-like. The combined effect is a reduced underestimation of the inner U for higher metallicities. Therefore, a spurious MI correlation may be induced by the geometric bias, as the ionization parameters of higher metallicity H II regions are less underestimated. The final slope of the MI relation depends on both the slope of the intrinsic relation and the amount of the geometric bias. It might also explain the saturation of the MI relation at a high SFR. If a cluster has a higher SFR and thus holds more stars, both the number of ionizing photons and the power of stellar winds increase, but their ratio remains roughly the same. This keeps the inner U unchanged (see Eq. (18)), but increases the inner radius. If the SFR is high enough, the metallicity dependence of the geometric bias would become negligible as the geometry becomes closer to plane-parallel overall, producing an apparent saturation of U. Still this depends on whether the intrinsic MI relation is flat or not, which is the weak point of this argument as we discuss shortly.

Figure 12 shows an example of the geometric bias. Here we measure the metallicities and ionization parameters of a series of spherical H II region models with the inner radii given by $$r_{\text{in}}=r_0 {(Z/Z_{\odot })}^{0.17}=10 {(Z/Z_{\odot })}^{0.17} pc,$$(21)

with our best-fit plane-parallel model. We used Eq. (15) and assumed the most extreme scaling relation for the mass loss rate, that is, $\dot{M}\propto {(Z/Z_{\odot })}^{0.85}$ (Puls et al. 2008). We chose the normalization to be r₀ = 10 pc, which is typical for a late-stage wind-driven bubble with t_age ≳ 1 Myr (e.g., Dopita et al. 2005; Dale et al. 2014; Geen et al. 2020). We treated the predicted line ratios from the spherical models as observations and added uncertainties typically found in MaNGA data to the line ratios. As expected, for H II regions with low ionization parameters at the inner radii, the measured U is closer to U_in and it barely depends on the metallicity. On the other hand, for H II regions with high U_in, there are positive MI correlations caused by the geometric dilution. Still the bias-induced correlation is insufficient to explain the MI correlation we measured in the MaNGA data, assuming the intrinsic MI relation is flat. Scaling down r₀ would steepen the slopes of the constant U_in lines, but even after reducing r₀ to 3 pc (which is too small for common H II regions with densities of 10 ~ 100 cm⁻³), we still cannot explain the median trend in the MaNGA data with a single constant U_in line. Making the metallicity dependence of the inner radius stronger would also yield a steeper slope, but it would require the power-law index in Eq. (21) to be on the order of 1, which is not evident in any dynamical model. Furthermore, once we consider the inner radius as a function of the metallicity, the intrinsic MI correlation is unlikely to be flat anymore. If we follow D06’s argument and use $U_{\text{in}}\propto {(Z/Z_{\odot })}^{-0.8}$, it is even harder to reconcile with the results from the MaNGA data.

To summarize, the geometric bias could potentially influence the MI trend observed at high ionization parameters. However, it is insufficient to explain the whole MI relation.

One might wonder if we can use a better definition of U to avoid this geometric bias. In Fig. 13, we compare the measured ionization parameters of the same set of spherical models described by Eq. (21) (using our fiducial model) with the ionization parameters computed at their inner radii (U_in), at their Strömgren radii (U_S), and using the volume-averaged values (U_VA), respectively. At nearly all input ionization parameters and metallicities, we have U_S < U_measured ≈ U_VA < U_in. The variation in U_measured (due to metallicity variation) at fixed U_S or fixed U_VA is much smaller than that at fixed U_in. Therefore, the geometric bias would be much weaker if we replace U_in with U_S or U_VA. This simple test shows that our plane-parallel fiducial model is actually roughly measuring the volume-averaged ionization parameter of these spherical models, although we have defined its ionization parameter at the inner surface of the ionized cloud. Even so, the definition of U_in has its own advantage. First of all, it is easy to set U_in as the input parameter for a photoionization model, while U_S and U_VA have to be calculated after the model is computed⁸. In addition, the dynamical models we compare also explicitly define the ionization parameter at the inner radii of H II regions (see Dopita et al. 2005, 2006). For the purpose of comparison, it is more straightforward to use the same definition. Furthermore, our test shows that the bias in the MI relation is not significant even when we consider an extreme scaling relation for the size of the wind-driven bubble in H II regions. Future works with more sophisticated combinations of dynamical models and photoionization models of H II regions will be helpful to provide a clearer picture on which definition is most useful.

We have seen that there is no simple physical picture explaining the entire MI relation in our data. However, we have not investigated the uncertainties associated with the self-consistency of the photoionization model, the method to derive the parameters, and the sample so far. We explore these points in the next section.

Fig. 12 Measured ionization parameters and metallicities for a spherically ionized cloud of which the inner radius varies with the metallicity accoridng to Eq. (21). Each colored line corresponds to H II regions with the same ionization parameter at their inner radii. The black line shows the median trend we measured in the MaNGA data, with the error bars indicating the standard derivations of log(U) in individual metallicity bins. The measurements were performed using N2, S2, R3, N2O2, and O3O2 line ratios.

Fig. 13 Comparisons between the ionization parameters measured by our plane-parallel fiducial model (Umeasured) and the ionization parameters under different definitions for a spherically ionized cloud (UX). The inner radius of the cloud was set to vary with metallicity according to Eq. (21). The ionization parameters of the spherical cloud are defined using the value at the inner radius (rin), the value at the Strömgren radius (rS), and the volume-averaged value, respectively. The dashed red line is a diagonal line for reference.

5 Discussions

As we have seen in Fig. 4, the derived correlation between the ionization parameter and metallicity strongly depends on the choice of photoionization models. The idea behind our derivations is that there exists a best-fit photoionization model surface sitting in the center of the data distribution in the high dimensional line-ratio space. This implicitly assumes that the ionization parameter and metallicity are two primary parameters that determine a 2D manifold embedded in high dimensions, and observed data tend to cluster around this manifold. The secondary parameters, on the other hand, contribute to the scatters around the manifold. Starting from this point, a few questions have yet to be answered.

First, it remains to be checked whether our best-fit model is still good enough if we include more emission lines. We already see that our model outperforms others in a 3D line-ratio space. In higher dimensions, there are more constraints on model parameters. Models that appear to fit the data in lower dimensions do not necessarily work in higher dimensions, which we have already seenin the case of 2D diagnostic diagrams. Second, given the existence of the intrinsic scatters induced by the secondary parameters, it is important to check whether our Bayesian approach is the right choice to capture the true distributions of the primary parameters. In many similar practices, people have invoked priors in their derivations. We need to carefully interpret the effect of priors. Finally, the choice of the sample certainly affects the derivations as it could change the mean values of the secondary parameters. Our best-fit model is based on MaNGA data. We need to check whether the MI correlation is still the same if it is derived from other samples. In what follows, we investigate the above points.

5.1 Offset between the observed O2 and model predictions

In Fig. 5, one can see that the ionization parameters predicted by the N2–S2–R3 method are systematically higher than those predicted by the N2O2–O3O2 method. This implies that the inclusion of the [O II] doublet no longer makes our model optimal in describing the data distribution. To check how much the model surface is offset from the data along the axes containing [O II] in high dimensions, we can define the logarithmic emission-line ratio difference as $$\Delta \log \text{ER}=\log {\text{ER}}_{\text{data}}-\log {\text{ER}}_{\text{model}}[{(\text{O/H})}_{3D},U_{\text{3D}}].$$(22)

For a given data point, ER_data is its extinction-corrected emission-line ratio (involving [O II]), and ${\text{ER}}_{\text{model}}[{(O/H)}_{3D},U_{3D}]$ is the same emission-line ratio predicted by the model with the metallicity and ionization parameter given by the N2–S2–R3 method. The left and middle panels of Fig. 14 show the O3O2 difference as a function of $12+log{(\text{O/H})}_{3D}$ and U_3D. The median difference in O3O2 reaches ~0.1 dex at high metallicities and ionization parameters. In comparison, the median O3O2 (not the difference, but the value itself) in our sample changes roughly by 0.4 dex as 12 + log(O/H) changes from 8.3 to 9.1. There is a clear trend of decreasing O3O2 as U_3D increases, indicating a systematic bias in the model-predicted O3O2. At high ionization parameters and high metallicities, our model tendsto overestimate O3O2, or underestimate O2.

However, there is a caveat for this interpretation. The O3O2 used here has an extra uncertainty from the extinction correction. It is observed that the overall Balmer decrement (and thus the extinction) is positively correlated with the metallicity. Therefore, it is possible that the extinction correction is increasingly overestimated as the metallicity (and thus the ionization parameter) increases. The right panel of Fig. 14 shows how the O3O2 difference changes as the Balmer decrement increases. One can see a similar relation as shown in the middle panel. The extinction correction we used is based on Balmer decrements and a Fitzpatrick (1999) extinction curve with R_V = 3.1 being adopted. If we instead use a Cardelli et al. (1989) extinction curve, the trend still remains. Interestingly, the absolute O3O2 difference would become smaller if we do not apply the extinction correction (but this time the difference is positive), which again implies that the extinction correction is overestimated.

In summary, there could be a systematic offset in both the model predictions concerning O2 and the estimation of the extinction based on Balmer decrements. Although the relative offset in O3O2 is not significant, it is correlated with the predicted ionization parameter and the Balmer decrements. We investigate the model predicted O2 and the potential systematic uncertainties associated with the extinction correction in detail in a future paper.

It is unsurprising that the model already becomes “imperfect” in a 4D space. Although photoionization models work well in predicting a large number of emission lines for individual well-observed H II regions (e.g., Baldwin et al. 1991; Jamet et al. 2005; Stasińska et al. 2013), finding a model that reproduces all strong emission lines for the general population of H II regions is notoriously difficult (see e.g., Sect. 6 of Law et al. 2021a, Sect. 3.2 of Mingozzi et al. 2020, and references therein). As a compromise, when using photoionization models with the Bayesian inference, people usually assume some nonflat priors for metallicities or ionization parameters, and adding terms representing the intrinsic scatters not reflected by the models (e.g., Pérez-Montero 2014; Blanc et al. 2015; Mingozzi et al. 2020). In the next section, we investigate the effects of these methods and their interpretations.

Fig. 14 O3O2data −O3O2model as a function of the gas-phase metallicity (left panel), ionization parameter (middle panel), and Balmer decrement (right panel). The median trend in each panel is shown as the black line, with the standard deviations indicated by the error bars.

5.2 Nonflat priors and intrinsic scatters

The introduction of nonflat priors changes the form of the posteriors. Using the likelihood we introduced in Eq. (5), the posterior is given by $$p(\theta |D,M)=\frac{p(\theta |M)p(D|M,\theta )}{p(D|M)},$$(23)

where p(θ|M) is the prior and p(D|M) is the normalization factor. Priors come in different forms, which usually assume the distribution of the parameter to follow a certain formwithin a given range, or correlate with some observables. Mingozzi et al. (2020), for example, used [S III]λλ9068, 9532/[S II]λλ6716, 6731 to set a prior for the ionization parameter. Photoionization models predict a very tight correlation between [S III]/[S II] and the ionization parameter, with little dependence on the metallicity (Kewley et al. 2019). Mingozzi et al. (2020) adopted the relation given by Diaz et al. (1991) and found that the smaller peak in the metallicity distribution predicted by the D13 model vanishes (see the lower left panel of Fig. 5). Meanwhile, the correlation between the ionization parameter and metallicity seemed to disappear at fixed stellar masses (although the overall correlation still exists). Despite the seeming advantage of using such a prior for “correcting” the model predictions, one caveat needs to be taken into account. Most of the current photoionization models (including ours) fail to predict the observed range of the [S III]/[S II] ratios. Under this circumstance, using the correlations between [S III]/[S II] and U set by photoionization models themselves as priors can introduce a bias. There is also a risk concerning self-consistency if the relation derived from a model is applied to correct the predictions of another model with different model assumptions. In this work, we have seen that the double-peak features in the metallicity distribution are likely a result of the deviation of the model surface from the dense region of the data distribution in the line-ratio space. The effect of priors is to regulate the mappings between themodel parameters and the data positions, which can also be achieved by adjusting the secondary model parameters and make the model surface match the data. We think the latter approach is more physically motivated. Even so, it would be interesting to see how the photoionization models for H II regions can be improved to correctly reproduce observed [S III]/[S II] and update the ionization parameter calibration, which we leave for future investigations.

Another potentially useful treatment in the Bayesian inference is the inclusion of model uncertainties. As we have mentioned, the variations in the secondary model parameters manifest themselves as intrinsic scatters in the data distribution. Since each photoionization model uses a single set of secondary model parameters, including terms of constant uncertainties for the model-predicted line ratios might help improve the fitting (Blanc et al. 2015). Indeed, if we add a systematic uncertainty of 0.1 dex for each line ratio, the agreement between the N2–S2–R3 method and the N2O2–O3O2 method becomes slightly better, and their predicted distributions for the ionization parameter are more consistent. This might be because the 2D method is more susceptible to intrinsic scatters as there are fewer degrees of freedom. Whereas we need to be cautious that the manually added uncertainties could hide the true discrepancy between different line ratios. An alternative way is to include some of the secondary parameters in the modeling as free parameters. In principle, this would reduce the impact of intrinsic scatters. However, more degeneracy would occur as more free parameters are considered. It would be necessary to add more emission lines and go to higher dimensions, which in turn requires models to accurately simulate more lines. No doubt such a self-consistent treatment is complicated and may be computationally expensive, but it might be worth doing if one wants to understand subtle correlations between different parameters.

5.3 Sample selection effect

Our derived positive correlation between the metallicity and ionization parameter is entirely based on the spatially resolved MaNGA H II regions. In this section, we discuss the impact of the sample selection on our result.

MaNGA’s sample galaxies have a relatively flat mass distribution (see Sect. 2), and this explains why a large fraction of our sample spaxels show high metallicities. The data we used are spatially resolved spaxels. For individual spaxels, the measured spectra are light-weighted, while the density distribution of the spaxels are heavily influenced by the on-sky areas of the H II regions in galaxies. This makes the MaNGA data more weighted toward the outer SF regions in galaxies compared to single fiber spectroscopic data. Kewley et al. (2019) suspect that the MI correlation depends on the scales of observations, with the single fiber spectroscopic data yielding anti-correlations as expected by the wind-driven bubble models, while the spatially resolved data give positive or no correlations. To check this possibility, we applied our method to the SF galaxies in the DR7 of SDSS (York et al. 2000; Abazajian et al. 2009). We obtained the emission line measurements using the updated code of Yan et al. (2006) with flux calibrations and zero-point corrections to the equivalent widths for emission lines (Yan 2011, 2018). The result is shown in Fig. 15. Clearly, there is a similar positive trend similar to what we see in the spatially resolved case. We thereby conclude that the discrepant MI correlations seen in previous works are not due to the difference between single-fiber observations and IFS observations or the difference between the centers of galaxies andoutskirts of galaxies. They are more likely a result of different metallicity and ionization parameter calibrators.

Even so, the spatial resolution of MaNGA is not enough to sample individual H II regions. The full width at half maximum (FWHM) of MaNGA’s point spread function (PSF) is roughly $2_{.}^{\prime \prime }5$ (Law et al. 2015), which corresponds to 1 ~ 2 kpc at the typical redshifts of MaNGA galaxies. As suggested by Sanders et al. (2020) and Mannucci et al. (2021), the observed line ratios from H II regions could change considerably if individual H II regions are resolved. There are two effects involved. First, the emission from the diffuse ionized gas (DIG) around H II regions can be blended with that from within H II regions. The DIG emission has enhanced low-ionization lines and can impact the determination of metallicities (Zhang et al. 2017). The DIG impact can largely be removed by applying a cut to the equivalent width of the Hα line, the surface brightness of the Hα line, or the fractional contribution to the total luminosity by young stars (e.g., Cid Fernandes et al. 2011; Sánchez et al. 2014; Zhang et al. 2017; Espinosa-Ponce et al. 2020). The majority of our sample spaxels have a large enough EW(Hα). Thus, the DIG is not likely to be dominant in our case. Second, within a single H II region, the emission line spectrum changes with the aperture size as well. The inner region would be featured by stronger high ionization lines but weaker low ionization lines. How exactly the spectrum changes would depend on the geometry of the H II region. The effective ionization parameter, being sensitive to geometry, would also change with the aperture size. The derived ionization parameter of the entire H II region should thus be related to its internal structure. As we have mentioned, the ionization parameter is not a direct observable most of the time and it is not well defined if the geometry of the H II region is unclear. Therefore, it is vital to have detailed analyses on well-resolved nearby H II regions to understand how well the geometry is reflected by the emission line measurements, which could be addressed by future surveys such as the Affordable Multi-Aperture Spectroscopy Explorer (AMASE, Yan et al. 2020) and SDSS-V/Local Volume Mapper (LVM, Kollmeier et al. 2017).

Last but not least, the redshift of the H II region sample could also impact the result. MaNGA’s sample is mainly composed of low-redshift galaxies (0.01 < z < 0.14). Sanders et al. (2016) and Sanders et al. (2020) found no evidence of evolution in the IM relation for their high redshift sample from the MOSDEF survey (Kriek et al. 2015), but the range of the MI relation they could measure is limited. Star-forming galaxies at higher redshifts form an offset sequence in the line-ratio space compared to the local ones (e.g., Shapley et al. 2005, 2015; Erb et al. 2006; Steidel et al. 2014). They could have higher electron densities (Kewley et al. 2013; Sanders et al. 2016), elevated N/O versus O/H relations (Shapley et al. 2015; Sanders et al. 2016), or harder ionizing spectra (Steidel et al. 2014), which would generally result in larger N2, S2, and R3 line ratios. The density effect on line ratios of H II regions is small compared to the changing SED effect (Ji & Yan 2020). To fit a high redshift sample, one needs to consider the average changes in all secondary parameters and the potential model degeneracy associated with the changes. The accuracy of the model predictions depends on how well these secondary parameters can be constrained. This requires a better understanding of the mappings between different secondary parameters and observed line ratios in H II regions. Precise and self-consistent photoionization modeling as well as carefully chosen visualizations of data distribution in multidimensional line-ratio space would be the key to solving the problem.

Fig. 15 MI correlation derived from the SF galaxies in the DR7 of SDSS. The counters represent the density distribution of the MaNGA data, with five density levels equally spaced in the logarithmic space from the 16th percentile to the 84th percentile.

6 Conclusions

In this work, we study the correlation between the gas-phase metallicity and ionization parameter (MI correlation) in general H II regions and investigate its physical origin. We measured the gas-phase metallicities and ionization parameters for H II regions in MaNGA data released in SDSS DR15 using four different photoionization models (L10, D13, B17, and JY20). We performed the measurements with Bayesian inference and calculated the weighted average values assuming flat priors. A total of five emission line ratios were used for the measurements, including N2, S2, R3, N2O2, and O3O2. We compared the results from different models and different combinations of line ratios. Our conclusions are summarized as the following:

• 1.

  When compared with the data distribution in the 3D line-ratio space spanned by N2, S2, and R3, our updated model (JY20 model) provides a better fit to the location and the shape of the central surface of the data.

• 2.

  The goodness of fit in the 3D space can reflect the consistency of the model predictions over different line ratios. This is confirmed by the measurements of metallicities and ionization parameters using different subsets of line ratios. For each model, we carried out three sets of measurements using two line ratios (N2O2 and O3O2), three line ratios (N2, S2, and R3), and five line ratios (N2O2, O3O2, N2, S2, and R3), respectively. Compared to the other three models considered in this work (L10, D13, and B17), our model provides the most consistent results no matter which set of line ratios is used.

• 3.

  With our best-fit photoionization model, we found a positive MI correlation, which is in contrast to the prediction of a wind-driven bubble model frequently used for describing the dynamical evolution of H II regions.

• 4.

  In orderto solve this discrepancy, we investigated the potential influence from star formation. We found that star formation activities alone are not enough to explain the MI correlation. Our results indicate that the variation in the SFR actually contributes to the scatter in the MI correlation.

• 5.

  We also study the bias brought by the assumed geometry in the photoionization model. Since the relative thickness of a spherical H II region would change with the metallicity due to changes in temperature and stellar wind power, measuring it using a photoionization model with a fixed geometry introduces a bias correlated with the metallicity. We found that this geometric bias can result in a slight positive MI correlation, but it is still not enough to fully explain the MI correlation in the MaNGA data.

As a concluding remark, the discrepancy between the measurements from photoionization models and the prediction from the dynamical model is still an open question. Both a better photoionization model that could self-consistently reproduce more emission line ratios, and a better dynamical model that could accurately describe the impact from stellar wind as well as other feedback mechanisms are vital for future studies. Equally important is the detailed analysis of the structures of the nearby H II regions. A future data release from the Affordable Multi-Aperture Spectroscopy Explorer (AMASE, Yan et al. 2020)and SDSS-V/Local Volume Mapper (LVM, Kollmeier et al. 2017) would be useful for resolving this issue.

Acknowledgements

We thank the anonymous referee, whose thoughtful suggestions improved the clarity of this work. We acknowledgesupport by NSF AST-1715898 and NASA grant 80NSSC20K0436 subaward S000353. RY acknowledges support by the Hong Kong Global STEM Scholar scheme and the Direct Grant of CUHK Faculty of Science. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

All Tables

All Figures

Fig. 1 Density distribution of MaNGA MPL-7 spaxels in the 3D space spanned by N2, S2, and R3. Our sample H II region spaxels are colored from yellow to purple, while the rest of the spaxels are colored from white to black. Two photoionization model surfaces are shown. The cyan model and the red model are our fiducial SF model and AGN model, respectively. We show three viewing angles that lead to different 2D projections. Specifically, panel a and panel b correspond to the [S II]- and [N II]- BPT diagrams. Panel c shows that the two model surfaces are separate in 3D and that they are connected by a continuous mixing sequence.

In the text

Fig. 2 Density distribution of the MaNGA data in the 2D P1–P2 diagram, where the relevant parts of the (interpolated) SF model (cyan grid) and AGN model (red grid) appear edge-on and well separated. The sample H II regions are colored from yellow to purple. Here, we only plotted the parts of the model surfaces that cover the middle 98% of the data along the P3 axis (line of sight), which is perpendicular to the P1 versus P2 plane. This projection corresponds to a line of sight at (θ, ϕ) = (36°, 219°) in the 3D space of (N2, S2, R3), where θ and ϕ are the polar angle and the azimuthal angle, respectively (Ji & Yan 2020).

In the text

	Fig. 3 Redshift and stellar mass distributions of our sample galaxies. The histogram shows a relatively flat distribution in stellar mass. The sample galaxies include both primary and secondary MaNGA galaxies. The former cover spatial regions out to 1.5 Re, while the latter cover spatial regions out to 2.5 Re.
In the text

Fig. 4 Different photoionization models viewed in the P1–P2 diagram and their predicted relations between the metallicity and ionization parameter. Leftmost column: Different SF-ionized models and MaNGA MPL-7 data in the P1 –P2 diagram. The models were interpolated and cut so that only the parts that cover the middle 98% of the data along the hidden P3 axis are shown. The spaxels with contributions from SF ionization greater than 90% are colored from yellow to green, while the rest of the spaxels are colored from white to black. Only the former were used to derive the metallicities and ionization parameters. Right columns (a) to (c): derived metallicities and ionization parameters using different SF-ionized models and different combinations of line ratios. In each panel, the green arrow marks the solar oxygen abundance adopted by the corresponding model.

In the text

	Fig. 5 Histograms showing the number distributions of 12 + log(O/H) and log(U) predicted by different SF-ionized models and different combinations of line ratios.
In the text

	Fig. 6 Median absolute Euclidean distances in the N2–S2–R3 3D space as a function of fitted metallicities for different SF-ionized models. The metallicities shown were derived based on different models using the N2–S2–R3 method, thus having different ranges. The shaded regions encompass the middle 68% of the data.
In the text

Fig. 7 Comparisons of photoionization models with different input parameters in the P1 –P2 diagram. The cyan model is our fiducial model. The red model has the fiducial abundances and depletion set, but uses the SED adopted by the B17 model. The purple model uses the fiducial SED, but adopts the abundances and depletion set of the B17 model. The blue model uses the fiducial SED, but adopts the abundances and depletion set of the D13 model. The modelsare interpolated and cut so that only the parts that cover the middle 98% of the data along the hidden P3 axis are shown. The density distribution of the H II region sample is plotted in green and yellow, while the rest of the MaNGA sample is plotted in black and white.

In the text

Fig. 8 Comparisons of pre- and post-depletion nitrogen prescriptions for the BY17 and the JY20 models. Both the input prescription and the output prescription modified by dust depletion are shown. The JY20 model adopts thesame pre-depletion nitrogen prescription as Dopita et al. (2013) and depletes the elements using the CLOUDY default depletion set. The B17 model follows Dopita et al. (2000) for both the pre-depletion nitrogen prescription and the depletion pattern. While the left panel shows the N/O ratios as functions of the absolute oxygen abundance 12 + log(O/H), the right panel shows the N/O ratios as functions of the relative oxygen abundance [O/H], normalized by the solar O/H adopted by each model. The solar oxygen abundances, 12 + log(O/H)⊙, adopted by Dopita et al. (2000) and Dopita et al. (2013) are 8.93 and 8.69 and indicated by the orange triangle and the cyan star in the plot, respectively.

In the text

Fig. 9 Correlation between the ionization parameter and Hα surface brightness (left), and the mass-metallicity relation color coded by the Hα surface brightness (middle) and ionization parameter (right). Left panel: 2D histogram of the logarithmic ionization parameter and logarithmic Hα surface brightness. The density distribution shows the results obtained by the JY20 model using all five line ratios. Thedashed red contours and the dotted-dashed black contours show the distributions obtained by using N2O2–O3O2 ratios and N2–S2–R3 ratios, respectively. The contour levels are in logarithmic scales and range from the 16th percentile to the 84th percentile. Middle panel: spatially resolved mass-metallicity relation of the sample spaxels color coded by themedian Hα surface brightness inside each bin. The black contours represent the density distribution of spaxels obtained by using all five line ratios. Right panel: same as the middle panel, but color coded by the median ionization parameter inside each bin.

In the text

Fig. 10 Dependence of the log(U) versus 12 + log(O/H) relation on the stellar mass surface density (left) and Hα surface brightness (right). The ionization parameters and metallicities shown were derived using all five line ratios. The contours represent five density levels equally spaced in the logarithmic space from the 16th percentile to the 84th percentile. The dotted-dashed lines indicate the median relations in different bins.

In the text

Fig. 11 Schematic plot of a wind-driven bubble model. Region a is the free-wind region filled with hypersonic stellar winds. Region b is the shocked-wind region, which sets the inner boundary for the H II region. Region c is composed of shocked gas. Region d is the H II region. Depending on the ionization rate, region d can be trapped within region c or extend beyond it. Region e is the ambient neutral medium. The relative sizes of the regions are not to scale.

In the text

Fig. 12 Measured ionization parameters and metallicities for a spherically ionized cloud of which the inner radius varies with the metallicity accoridng to Eq. (21). Each colored line corresponds to H II regions with the same ionization parameter at their inner radii. The black line shows the median trend we measured in the MaNGA data, with the error bars indicating the standard derivations of log(U) in individual metallicity bins. The measurements were performed using N2, S2, R3, N2O2, and O3O2 line ratios.

In the text

Fig. 13 Comparisons between the ionization parameters measured by our plane-parallel fiducial model (Umeasured) and the ionization parameters under different definitions for a spherically ionized cloud (UX). The inner radius of the cloud was set to vary with metallicity according to Eq. (21). The ionization parameters of the spherical cloud are defined using the value at the inner radius (rin), the value at the Strömgren radius (rS), and the volume-averaged value, respectively. The dashed red line is a diagonal line for reference.

In the text

	Fig. 14 O3O2data −O3O2model as a function of the gas-phase metallicity (left panel), ionization parameter (middle panel), and Balmer decrement (right panel). The median trend in each panel is shown as the black line, with the standard deviations indicated by the error bars.
In the text

	Fig. 15 MI correlation derived from the SF galaxies in the DR7 of SDSS. The counters represent the density distribution of the MaNGA data, with five density levels equally spaced in the logarithmic space from the 16th percentile to the 84th percentile.
In the text