© The Authors 2024

1. Introduction

Metallicity is a key observable of galaxies, as element abundances and their evolution across cosmic epochs provide unique information on the physical processes driving the evolution of galaxies (see, e.g., Maiolino & Mannucci 2019 for a review). In particular, since all elements heavier than lithium and up to iron and nickel are created in stellar cores, the amount of metals provides a measure of the integrated star-formation history of galaxies. Furthermore, the existence of scaling relations between metallicity and galaxy properties such as the mass-metallicity (e.g., Tremonti et al. 2004) and fundamental metallicity relations (e.g., Mannucci et al. 2010) also indicates that metallicity traces the infall and outflow of gas to and from galaxies. Metallicity can therefore be used to constrain feedback, fueling processes and growth timescales in galaxies (e.g., Mannucci et al. 2009; Cresci et al. 2010; Kewley et al. 2019b; Kumari et al. 2021; Curti et al. 2023 and references therein).

The metallicity of a galaxy and of the surrounding medium can be measured from a range of observables: from the integrated emission of stars, the warm (T_e ∼ 10⁴ K) ionized interstellar medium (ISM), the hot (T_e ∼ 10⁶ K) ionized intracluster (ICM) and circumgalactic medium (CGM), and the absorption features caused by cold and warm ICM or CGM on the spectra of distant quasars and galaxies. In this paper, we focus on the gas-phase metallicity of the ISM warm-ionized gas (T_e ∼ 10⁴ K) mainly located within the HII regions surrounding massive O and B stars but present also in the form of diffused ionized gas (DIG).

The emission line spectrum of photoionized gas can be used to constrain the properties of the ionizing source(s) and provides information on gas physical properties, including density, temperature, and metallicity. In particular, the measurement of metallicity depends on the density, temperature, and ionization status of the ionized gas (e.g., Osterbrock & Ferland 2006; Netzer 2008; Stasińska 2009; Kewley et al. 2019b). Three methods are usually adopted to estimate gas metallicity: two observational ones, with one based on direct electron temperature determinations and the other on metal recombination lines, and a third theoretical one based on photoionization models (see, Section 3 of Maiolino & Mannucci 2019 and references therein for a detailed description and discussion on these methods).

The first observational method, the so-called T_e method, is based on density and temperature estimates from the ratios of collisionally excited forbidden emission lines. It assumes that the complex stratification of the physical properties in the ionized gas volume can be modeled with a few (two to four) homogeneous regions characterized by constant values of electron density (N_e), temperature (T_e), and ionization fractions of any ion X⁺ⁱ. The N_e and T_e of each region can then be determined from one or more line ratios sensitive to density, such as [S II]6730/6716, [O II]3729/3728, and temperature, such as [O III]4363/5007, [O II]7330/3728, [N II]5755/6548, and [S III]6312/9532 (e.g., Osterbrock & Ferland 2006; Berg et al. 2020; Mingozzi et al. 2022 and references therein). The emissivity J_k, l of a specific emission line originating from the k, l transition, whether due to collisional excitation or recombination, can be determined using these N_e and T_e values. Subsequently, the relative abundance of the emitting ion can be obtained from the reddening-corrected observed ratio between the emission line of X⁺ⁱ and a hydrogen line, typically Hβ. Both lines are assumed to be emitted from the same volume ΔV. In formulae,

$$\begin{array}{cc}\hfill L_{k,l}(X^{+i})& \simeq N(X^{+i})N_{\mathrm{e}}{J}_{k,l}(X^{+i};N_{\mathrm{e}},T_{\mathrm{e}})\mathrm{\Delta }V(X^{+i}),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill L(\mathrm{H}\beta )& \simeq N({\mathrm{H}}^{+})N_{\mathrm{e}}{J}_{\mathrm{H}\beta }(N_{\mathrm{e}},T_{\mathrm{e}})\mathrm{\Delta }V(\mathrm{H}\beta ),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{N(X^{+i})}{N(\mathrm{H})}& \simeq \frac{L_{k,l}(X^{+i})}{L(\mathrm{H}\beta )}\frac{J_{\mathrm{H}\beta }(N_{\mathrm{e}},T_{\mathrm{e}})}{J_{k,l}(X^{+i};N_{\mathrm{e}},T_{\mathrm{e}})},\hfill \end{array}$$(3)

where L is a line luminosity, and it is assumed that N(H⁺)≃N(H) and ΔV(X⁺ⁱ)≃ΔV(Hβ). The abundance of an element X can then be determined by adding up the abundances of its ions:

$$\begin{array}{c}\hfill \frac{N(X)}{N(\mathrm{H})}={\displaystyle \sum _i}\frac{N(X^{+i})}{N(\mathrm{H})}.\end{array}$$(4)

Reliable values of T_e are needed because of the strong exponential dependence of J_k, l on temperature. Hereafter, we refer to these abundances as T_e-abundances.

The second observational method is based on metal recombination lines and is similar to the aforementioned approach, with the only fundamental difference being that the line emissivity J_k, l(X⁺ⁱ; N_e, T_e) has only a weak and roughly linear dependence on temperature, whose exact value is therefore much less important than for the T_e method (e.g., Peimbert & Peimbert 2013; Esteban et al. 2014; Toribio San Cipriano 2017).

Although commonly considered observational methods, both approaches necessarily rely on photoionization calculations to estimate ionic temperatures and abundances not directly accessible from the observations (see, e.g., Gutkin et al. 2016; Amayo et al. 2021 and references therein). Moreover, abundances derived from the T_e method are known to disagree with those derived from metal recombination lines. A possible explanation is the existence of temperature fluctuations, which mostly affect the determinations based on collisionally excited lines (e.g., Peimbert 1967; Peimbert & Peimbert 2013; Toribio San Cipriano 2017; Nicholls et al. 2020; Méndez-Delgado et al. 2023b). Also, especially for distant galaxies, the temperature is usually derived for only one of the zones of the HII regions, and locally calibrated relations are used to extrapolate the temperatures of the other zone(s) (e.g., Curti et al. 2017), potentially introducing large uncertainties.

The theoretical method is instead based on a comparison between the observed line ratios and those predicted from photoionization models. In practice, given a set of observed emission lines, a suite of photoionization models (single-cloud models) for different ionizing continua and physical properties of the gas (including metallicity) is created to best match the observed emission lines, thus determining the gas metallicity among the other physical properties (e.g., López-Sánchez et al. 2012; Nicholls et al. 2012; Binette et al. 2012; Pilyugin et al. 2012; Dopita et al. 2013; Blanc et al. 2015; Pilyugin & Grebel 2016; Vale Asari et al. 2016; Pérez-Montero & Amorín 2017; Strom et al. 2018; Mignoli et al. 2019; Papovich et al. 2022). Hereafter, we refer to this method as the “theoretical method” and to the metal abundances derived with it as the “model abundances”.

Overall, all these methods require the measurement of several emission lines to constrain metallicities. The T_e method, in particular, requires the measurement of the so-called auroral lines, which are needed for temperature determination but are usually quite faint (less than a few percent of the brightest lines). Notable examples of auroral lines are [S II]4069, [O III]4363, [O II]7330, [N II]5755, and [S III]6312. Instead, most photoionization models, at least in the extragalactic field, exclusively use the brightest emission lines, which are the only ones usually detected. These are such lines as Hβ; Hα; [O III]5007; [O II]3727, 3729; [N II]6584; [S II]6716, 6730; and [S III]9532, and they are usually referred to as strong lines.

For the sources where it is possible to measure only these strong emission lines, it is customary to use the so-called strong-line method, where metallicity is estimated from one or more strong line ratios whose dependence on metallicity has been calibrated using T_e-based and/or Model-based metallicities from HII regions and galaxies with high S/N spectra (e.g., Kewley & Dopita 2002; Pettini & Pagel 2004; Tremonti et al. 2004; Kewley & Ellison 2008; López-Sánchez et al. 2012; Dopita et al. 2013; Curti et al. 2017; Kewley et al. 2019b; Sanders et al. 2023). These line ratios have clearly defined metallicity dependencies because of existing correlations between the physical parameters of gas clouds, such as the anti-correlation between the ionization parameter and metallicity (see, e.g., Carton et al. 2017; Ji & Yan 2022) or the star-formation rate (e.g., Papovich et al. 2022). These methods, and especially the strong-line one, are used extensively in studies of galaxy evolution, providing evidence for the well-known metallicity scaling relations mentioned above.

However, there are still important open issues in metallicity estimates, as contradicting calibrations of the strong-line method exist in the literature. These disagreements are due to the different reference HII or galaxy samples used by the different authors, the differences in the application of the T_e method (i.e., corrections for temperature fluctuations, ionization corrections), and the different assumptions used for the photoionization models (abundance ratios, ionizing continua, dust properties, etc.). Moreover, the T_e method and theoretical methods based on photoionization models usually provide discrepant results, with the latter suggesting metallicities that can be larger by 0.3 dex or more, leading to significant uncertainties on all scaling relations involving metallicity (e.g., Kewley & Ellison 2008; Stasińska 2010; López-Sánchez et al. 2012; Curti et al. 2020).

In general, it is clear that both approaches aim to reproduce the complexity of ionized clouds under oversimplified assumptions. Indeed, a single HII region is a complex structure made of condensations and filaments of gas with non-trivial gradients of physical properties, possibly exposed to the ionizing photons emitted by different sources, and the emission line spectrum of a star-forming galaxy is the superposition of several hundreds of HII regions, each with different physical conditions (e.g., Charlot & Fall 2000).

Overall, the recent important advances in spectroscopic instrumentation, in terms of both quality and quantity of the data brought forward by integral field spectrographs, have not been matched by similar improvements in the methods used for metallicity measurements. On one hand, the T_e method has not changed in decades, apart from (semi)empirical corrections or improved theoretical ionization corrections (e.g., Méndez-Delgado et al. 2023b. On the other hand, the physics within photoionization models has greatly improved, but the modeling of the observed data is still mostly based on single-cloud models or on simple combinations of multiple clouds. There have been a few attempts to model the complexity of real gaseous nebulae using, for example, Monte Carlo techniques (Ercolano et al. 2003; Jin et al. 2022), but it is at the expense of the accuracy of the physical processes included in codes. Still, the successes and improvements of photoionization models over the past years indicate that they could represent the way forward to improving metallicity estimates.

Here we present HOMERUN (Highly Optimized Multi-cloud Emission-line Ratios Using photo-ionizatioN), a new approach to photoionization modeling that builds on a grid of single-cloud photoionization models. In contrast to other models presented in the literature, it combines single-cloud photoionization models with a weighting function that is determined from observations, not assumed a priori. In Section 2, we describe the current approaches to photoionization modeling, while in Section 3, we present our new HOMERUN model. In Sections 3.6 through 5, we describe and discuss its application to the spectra of samples of HII regions. Finally, in Section 6, we show how our model resolves the apparent discrepancy between the observed metallicities of HII regions and stars in the Milky Way. We summarize our findings and draw our conclusions in Section 7.

2. Photoionization models

There are broadly two kinds of approaches in photoionization modeling of HII regions, and in both cases the publicly available CLOUDY¹ (Ferland 1993; Ferland et al. 2017) and MAPPINGS² (Sutherland & Dopita 1993; Dopita & Sutherland 1996) codes are widely used. One approach is based on comparing single-cloud constant-pressure models with observations (e.g., Dopita et al. 2016), also taking into account the star-formation history (e.g., Charlot & Longhetti 2001; Chevallard et al. 2018). The other is based on the combination of different single-cloud constant-density models through an assumed weighting function. This is, for instance, the case of the so-called locally optimally emitting clouds (LOC; Baldwin et al. 1995; Korista et al. 1998; Nagao et al. 2006b). The former approach is extensively used in the analysis of galaxy spectra and in the determination of metal abundances in their ISM, while the latter approaches are mostly used in the analysis of the emission lines from the broad-line regions of active galactic nuclei although it has also been applied to the spectra of star-forming galaxies (Richardson et al. 2016). These models are successful in reproducing the general emission-line properties of populations of galaxies and active galactic nuclei, but they struggle in reproducing the emission-line spectra of single objects: different line ratios in the same galaxy are often reproduced by models with different input parameters because the models fail to reproduce the observed wide range of excitation levels which is manifested by the existence of strong emission lines from different ionization states, like [O III]λ5007, [O II]λ3727 and [O I]λ6300 (e.g., Martin 1997; Iglesias-Páramo & Muñoz-Tuñón 2002; Stasińska et al. 2015; Ramambason et al. 2020; Lebouteiller & Ramambason 2022). More in detail, a single-cloud model is described by a slab of gas in plane-parallel or spherical geometry, ionized by a central source. A single-cloud model therefore is defined by the spectrum of the ionizing continuum, that can vary widely depending on the model used for the ionizing source (stellar age, stellar and gas-phase metallicity, presence of binary stars, etc.); the hydrogen density N_H at the illuminated face of the cloud; the radial distance r_in of the illuminated face of the cloud from the ionizing source; the ionization parameter $U=Q(\text{H})/(4\pi r_{\text{in}}^{2}c{N}_{\text{H}})$, where Q(H) is the rate of ionizing photons and c the speed of light (sometimes different authors use a different definition of U which is computed, e.g., at a different radius or as a volume average, see for example Charlot & Longhetti 2001; Gutkin et al. 2016; Plat et al. 2019); the abundances of all the elements relative to hydrogen; whether the cloud is at constant pressure or constant density; whether the cloud contains dust, and what are its properties and the fraction of chemical elements depleted into dust; the assumption of stationary equilibrium, that is, no time evolution; a criterion to define the outer edge of the cloud, which could be a limit on any parameter among hydrogen column density, ionization fraction, and electron temperature.

The model is integrated starting from the illuminated face of the cloud, numerically solving the equations of ionization equilibrium for all ions, the energy balance equation, and the radiative transfer equation; all other relevant physical processes are also considered (e.g, Osterbrock & Ferland 2006; Stasińska 2009; Netzer 2008; see volume 3 of Cloudy manual “Hazy” for an exhausting description³). The model predicts, among other things, the fluxes or luminosities of all the considered transitions which can then be compared with observations. The metallicity of the model which provides the best match with the observations is then the estimated gas metallicity.

The ionization level of the gas at a given point in the cloud is set, among other things, by the value of the ionization parameter at that point, that is, computed from Q(H) and N_H at r, where r is the distance from the ionizing source. In constant-pressure cloud models, as r increases, the temperature decreases and the density grows, leading to a radial decrease in U which is faster than that predicted by constant-density models. This wider U range allows for a better match with the observed emission lines from different ionization stages, and is the reason why constant-pressure, single-cloud models are typically more successful than constant-density ones (e.g., Dopita et al. 2002; Stern et al. 2014; Kewley et al. 2019a).

In order to provide a better match to the observed emission lines, combinations of single-cloud models have been considered. One notable example is that of the LOC models mentioned above. Here, a grid of single clouds is computed, each characterized by a single value of N_H and $\mathrm{\Phi }(\text{H})=Q(\text{H})/(4\pi r_{\text{in}}^2)$ (the ionizing photon flux), and the emission lines predicted by the single-cloud models are combined with a weighting function which is usually proportional to $\mathrm{\Phi }{(\text{H})}^{-\alpha }{N}_{\text{H}}^{-\beta }$, with α and β two slopes of order of unity which can be free parameters of the model (Nagao et al. 2006b). As in the case of single-cloud models, the metallicity of the model which provides the best match is then the estimated gas metallicity. A different approach, mostly applied to active galactic nuclei, assumes that the clouds are confined by radiation or magnetic pressure and that their covering factor varies as a power law of distance r from the ionizing source, resulting in N_H and U which vary similarly as a power law of r (e.g., Netzer 2020). LOC models are better suited to reproduce the complexity of a real HII region and its irregular geometrical and ionization structure, but the use of a fixed weighting scheme is still a severe limitation.

3. The HOMERUN model

3.1. The basic assumption: A combination of multiple single-cloud models

The luminosity of an emission line l of ion X⁺ⁱ from either recombination or collisional excitation can be written as

$$\begin{array}{c}\hfill L_l{(X^{+i})}_{\mathrm{obs}}={\displaystyle {\int }_V}{N}_{\mathrm{e}}N(X^{+i})J_l(X^{+i},T_{\mathrm{e}})\mathrm{d}V,\end{array}$$(5)

where V is the volume of the ionized gas, N_e and N(X⁺ⁱ) are the electron and ion densities, respectively, T_e is the electron temperature and J_l(X⁺ⁱ, T_e) is the emissivity of line l. The line emissivity J_l is the same as in Eq. (1) but here we simply indicate the transition as l and we do not explicitly mention the usually weak dependence on the electron density. N_e, N(X⁺ⁱ) and T_e are physical quantities which depend on the position within the volume, yet the luminosity, which is an integrated quantity, does not depend on the detailed density and temperature distribution of the gas but on their weighted average values over the entire gas volume. Therefore, any model that has the correct volume, average density and temperature of the X⁺ⁱ region, will reproduce L_l(X⁺ⁱ)_obs, even if the spatial distribution of the gas physical properties does not match the real one.

In this study we assumed that the ion emission line luminosity defined in Eq. (5) can be sufficiently well approximated by a linear combination of m single-cloud models, each characterized by its own ionization parameter and density

$$\begin{array}{cc}\hfill L_l{(X^{+i})}_{\mathrm{mod}}& ={\displaystyle \sum _{j=1}^m}{w}_j{L}_l{(X^{+i})}_{\mathrm{mod},j}\hfill \end{array}$$(6)

$$\begin{array}{cc}& ={\displaystyle \sum _{j=1}^m}{w}_j{\displaystyle {\int }_{V_j}}{N}_{\mathrm{e}}N(X^{+i})J_l(X^{+i},T_{\mathrm{e}})\mathrm{d}V,\hfill \end{array}$$(7)

where L_l(X⁺ⁱ)_mod, j is a set of “simple” single-cloud photoionization models that can be considered as the basic building blocks of the model. The non-negative weights w_j can either be assumed (like in the LOC models) or can be left as free parameters and determined from fitting the observed lines. The main underlying assumption is that if a model is able to accurately reproduce the luminosity of many different emission lines from a few different ions, then it has the correct combinations of volume, average density and temperature for each ion, regardless of their real spatial distribution. With respect to single-cloud models, we expect this approach to provide a more realistic description of the complexity of the warm ionized ISM.

3.2. The single-cloud models

The single-cloud models used in the following were computed using CLOUDY v22.01 (Ferland et al. 2017) with the following input parameters. The ionizing continuum was provided by BPASS simple stellar population models (v2.3, Byrne et al. 2022; Stanway & Eldridge 2018) including binary stellar evolution, characterized by a Kroupa (2001) initial mass function from 0.1 to 300 M_⊙, [α/Fe] = +0.0, metallicities solar and 0.1 solar, a single burst of star formation and ages (log yr) 6, 6.3, 6.7 7.1 (for solar met.) and 6, 6.4, 6.7, 7.3 (for 0.1 Z_⊙). Abundances Z varied from −2 to +0.4 in steps of 0.2 dex relative to solar (the solar abundance pattern was taken from Grevesse et al. 2010); nitrogen was rescaled as a secondary element with respect to oxygen as [N/O]= − 1.5 + 1.4 * ([O/H]−[O/H]_c) for [O/H]≤[O/H]_crit with [O/H]_c = 8.3 (see, e.g., Fig. 49 of Maiolino & Mannucci 2019), and we also considered cases in which [N/O] varied by ±0.6 dex with respect to this scaling relation; sulphur scaled linearly with oxygen following the solar [S/O] abundances ratio, but we also considered the case in which [S/O] varied by ±0.3 dex. We considered gas with and without dust; dust was included with the “grains ism Z log” CLOUDY command line, where Z is the log of the abundance relative to solar; we then used the command “metals deplete” to take into account the depletion of metals into dust (adopted depletion factors are the default values of Cloudy; see Table 7.8 and references therein of Hazy1 of Cloudy v22.1). A grid of single-cloud models in U and N_H was then computed with log U from −5 to −0.5 in steps of 0.5 dex and log N_H from 0 to 7 in steps of 1.0 dex. Here, U represents the ionization parameter at the inner face of the cloud. All single-cloud models had constant density, a plane-parallel open geometry and were radiation bounded with the following integration stopping criteria: maximum H column density 10²³ cm⁻², minimum temperature 3000 K, minimum fraction of electrons 0.05.

Overall, these choices of parameters required the computation of 2 × 4 × 25 × 3 × 3 × 10 × 8 × 2 = 134, 784 CLOUDY models. Then, to increase the sampling, the model line luminosities relative to Hβ were log-linearly interpolated across metallicity in steps of 0.02 dex.

The choice of constant density was made in order to have the simplest possible physical structure of a single-cloud model and it did not affect the final results. Indeed, we verified that constant pressure models can be reproduced by a combination of a constant density models.

The stellar metallicity, which directly affects the shape of the stellar continuum, was not tied to the metallicity of the gas to allow more freedom in the modeling. We additionally verified that the fit results were not much sensitive to the stellar continuum shape which could not be accurately constrained; this issue is beyond the scope of this paper and will be further investigated elsewhere.

A single-cloud CLOUDY model is characterized by a total Hβ surface brightness. Considering the radial symmetry in a spherical model, it is possible to obtain the volume that the cloud must occupy to reproduce the observed luminosity. Models with low density and low ionization parameters may therefore become too large to be physical or simply to be included in the apertures from which the spectra are obtained. It is easy to show that this implies that these low-density, low-ionization-parameter models have an Hβ surface brightness which is too low compared to the observed one. The CHAOS HII regions spectra described below (Sect. 5) were obtained through apertures (e.g., slits) whose linear size is ∼10″ and the observed Hβ surface brightness is, on average, of the order of ∼10⁻⁴ erg cm⁻² s⁻¹ arcsec⁻². Therefore, in this paper, we considered ∼10⁻⁵erg cm⁻² s⁻¹ arcsec⁻² as the minimum surface brightness that a single-cloud model must have to be retained in the computations. This ensured that single-cloud models with weights larger than 10% of the total account for at least 10% of the observed surface brightness. The effect of this choice was to discard between 15 and 20% of the single-cloud models in the low density, low ionization parameter corners of the U, N_H grids computed for a given abundance Z and ionizing continuum, this fraction increasing with increasing metallicity. In principle, one should adopt an iterative procedure: fit the data – as described below – with all single-cloud models in a grid and find the model weights, compute the surface brightness of the single-cloud models given their weights and discard the models with surface brightness lower than the observed surface brightness multiplied by the model weight, redo the fit and iterate until no single-cloud model is discarded.

3.3. The loss function ℒ

The agreement of a multi-cloud model with the observations was quantified with the loss function ℒ defined below. We considered a grid of single-cloud constant-density, radiation-bound models in ionization parameter U and H-density N_H characterized by the spectrum of the ionizing source, 𝒮_ν, the gas metallicity, parameterized with the oxygen abundance, 𝒜_O, and the relative abundances of metal elements, 𝒵. We identified this grid of models with 𝒢 = 𝒢(𝒮_ν, 𝒜_O, 𝒵). Each single-cloud model of the grid 𝒢 predicts the observed emission line luminosities which we normalised to an arbitrary dimensionless Hβ luminosity of, e.g., 100. There are no limitations to the emission lines that can be considered in the analysis provided that they are computed by CLOUDY. All emission lines were then combined with weights which are found by minimising the loss function (ℒ)

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\left(\frac{L_l{(X^{+i})}_{\mathrm{obs}}-L_l{(X^{+i};w_1,w_2,\cdots ,w_m)}_{\mathrm{mod}}}{\mathrm{\Delta }{L}_l{(X^{+i})}_{\mathrm{obs}}}\right)}^2,\end{array}$$(8)

where n is the number of observed emission lines and m is the numbers of single-cloud models in the grid. ΔL_l(X⁺ⁱ)_obs can be either the error on the observed emission line luminosity or what we arbitrarily considered an acceptable discrepancy between the observed and model luminosities: in the following, we considered ΔL_l(X⁺ⁱ)_obs as the maximum of the two.

The values of the loss function can be approximately interpreted by assuming that the average relative discrepancy between models and observations is ε and that for all emission lines and all ions

$$\begin{array}{cc}& L_{l,\mathrm{mod}}=L_{l,\mathrm{obs}}(1+{\epsilon }_r{\delta }_l),\hfill \end{array}$$(9)

$$\begin{array}{cc}& \mathrm{\Delta }{L}_{l,\mathrm{obs}}={\epsilon }_a{L}_{l,\mathrm{obs}},\hfill \end{array}$$(10)

where δ_l are random numbers extracted from a normal distribution with zero average, and unit standard deviation. We distinguished between “real” relative discrepancy ε_r between the models and the observations, and the acceptable discrepancy ε_a, that is, the one we adopted when computing the loss function. Therefore the loss function is given by

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{n}{\displaystyle \sum _l}{\left(\frac{{\epsilon }_r{\delta }_l}{{\epsilon }_a}\right)}^2.\end{array}$$(11)

For example, if the acceptable discrepancy is equal to the real average discrepancy, the loss function is ℒ ≃ 1.0. Conversely, assuming an acceptable discrepancy ε_a = 0.1, real average discrepancies of ε_r = 0.01, 0.05, 0.1, 0.11, 0.15, 0.2 will correspond to loss function values of ℒ ≃ 0.01, 0.25, 1.0, 1.2, 2.3, 4.0, respectively. This numbers should serve as a guideline to estimate the real average discrepancy for any given value of the loss function.

3.4. Varying the N/O and S/O abundance ratios

Usually, the relative abundances of all elements relative to oxygen are constant and fixed to some predetermined value, like the solar one. However, in our model fitting, we allowed for a scaling of the relative abundances of the elements, especially for N and S, which have the strongest observed lines in the optical. This choice stems from the evidence presented in the literature and is physically expected given the wide range of conditions in which stars and galaxies form and evolve (Steidel et al. 2016; Morisset et al. 2016; Strom et al. 2018; Maiolino & Mannucci 2019; Berg et al. 2020).

In order to allow for different relative abundances other than the solar ones used in the computation of the single-cloud models, we rescaled each set of model emission lines of a given element by a scaling factor which is a free parameter of the fit. The only element to which we did not apply such rescaling was oxygen since some of its forbidden emission lines provide the main cooling mechanism of the ionized zone, and therefore the line luminosities are not trivially proportional to abundance. Considering, for instance, observed lines from oxygen, N, and S, the loss function becomes

$$\begin{array}{cc}\hfill \mathcal{L}& ={\displaystyle \sum _{i=1}^{n_{\mathrm{O},\mathrm{H}}}}{\left(\frac{L_l{(X^{+i})}_{\mathrm{obs}}-L_l{(X^{+i};\mathbf{w})}_{\mathrm{mod}}}{\mathrm{\Delta }{L}_l{(X^{+i})}_{\mathrm{obs}}}\right)}^2\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^{n_{\mathrm{N}}}}{\left(\frac{L_l{(X^{+i})}_{\mathrm{obs}}-\mathcal{N}\ast L_l{(X^{+i};\mathbf{w})}_{\mathrm{mod}}}{\mathrm{\Delta }{L}_l{(X^{+i})}_{\mathrm{obs}}}\right)}^2\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^{n_{\mathrm{S}}}}{\left(\frac{L_l{(X^{+i})}_{\mathrm{obs}}-\mathcal{S}\ast L_l{(X^{+i};\mathbf{w})}_{\mathrm{mod}}}{\mathrm{\Delta }{L}_l{(X^{+i})}_{\mathrm{obs}}}\right)}^2,\hfill \end{array}$$(12)

where w = (w₁, w₂, …, w_m) are the weights of the m single-cloud models, n_O, H indicates that the sum is made on all H and O lines while n_N and n_S indicate summation on N and S lines respectively, and 𝒩, 𝒮 are the scaling factors which are free parameters of the fit.

This procedure assumes that the N and S emission lines scale linearly for small variations of the abundances of these elements. This is not always true, but we verified that if the abundance of N is varied by less than ±0.3 dex with respect to the solar value, and S by less than ±0.15 dex, the line luminosities scale linearly to better than ∼10%.

Therefore, to correct for the non-linearity of the N and S line intensities with the elements’ abundances, we considered only the models where the scaling factors 𝒩 and 𝒮 are smaller than ±0.3 and ±0.15 dex, respectively. With the choices of [N/O] and [S/O] adopted in the computation of the CLOUDY models, this ensures the accuracy of our computations when varying [N/O] and [S/O] within ±0.9 and ±0.45 dex of the solar values, respectively.

3.5. The fitting procedure

In the fitting procedure we considered first a U, N_H grid of single-cloud models computed with a given ionizing continuum and set of element abundances, 𝒢(𝒮_ν, 𝒜_O, 𝒵), and found the multi-cloud model which provided the minimum value of the loss function ℒ_min(𝒮_ν, 𝒜_O, 𝒵). For a given set of scaling factors 𝒩, 𝒮, the best weights (w₁, w₂, …, w_m) were found by solving the non-negative least square (NNLS) problem, in a similar way as in the pPXF python procedure by Cappellari & Emsellem (2004) and by Cappellari (2017). Indeed, this model fitting is similar in principle to the fitting of the stellar continua of galaxies where many hundreds or thousands of stellar continuum templates are combined with non-negative weights. Here, emission-line templates from the single-cloud models are combined with non-negative weights. Since in the single-cloud models the emission lines are normalized such that L(Hβ) = 100, the sum of the weights will be equal to the total observed Hβ luminosity. Finally, in finding the best weights we did not apply any regularization, a standard mathematical procedure to mitigate the high-frequency variations in the solutions (e.g., Cappellari 2017).

We repeated for all U, N_H grids of single-cloud models available, iterating over abundances and ionizing continua and, finally, we found the minimum of all ℒ_min values, thus identifying the combination of 𝒮_ν, 𝒜_O, 𝒵 which provided the best fit to the observed line ratios. A schematic representation of the fitting procedure is presented in Fig. 1.

Fig. 1. Schematic representation of the fitting procedure. For any U, NH grid of single-cloud models, we computed the weighted combination of clouds in the grid that gives the best match with the observed lines. Usually, only a small fraction of the clouds in the grid have non-zero weights. The goodness of the fit was then quantified with a value of the loss function ℒ. The procedure was repeated for all U, NH grids computed with different gas metallicities, and the behavior of ℒ as a function of metallicity was obtained. The metallicity of the observed spectrum is identified by the minimum value of ℒ. This procedure can be repeated with different ionizing continua, thus estimating the effect of varying the ionizing continuum on the metallicity estimate.

3.6. An example of a model fit

We describe the model by presenting an example of its application to the emission line fluxes of an HII region. The test HII region is randomly extracted from the sample presented by Zurita et al. (2021) and it has all the emission lines listed in the bottom right panel of Figure 2 measured with Signal-to-Noise ratio (S/N) larger than 5. The results from the application of the model to the entire Zurita et al. (2021) sample of HII regions will be presented in Section 4.

Fig. 2. Example of model fit to the emission line spectrum of an H II region. Top-left panel: variation of the loss function ℒ as a function of the oxygen abundance, 12 + log(O/H), for the emission line spectrum of the test HII region. In this scale, the solar metallicity is 8.69. The blue vertical dashed line represents the metallicity of the model with the minimum value of ℒ. Bottom-left panel: variation of the N and S abundances as a function of the oxygen one. The dashed lines represent the values adopted for the computation of the photoionization models, while the filled circles correspond to the abundances rescaled as described in the text. Top-right panel: grid of single-cloud models in log(U) and log(NH). The colors represent the weights of each single-cloud model, as indicated by the color bar, when ℒ reaches the minimum value. In this case only four out of 72 models of the grid have non-zero weights. The black dot represents the weighted density and ionization parameter of the single-cloud models. Bottom-right panel: ratio between the model and observed line luminosities. The gray areas represent the acceptable discrepancy of ±10% or, when larger, the relative error on the observed value.

We fit the observed lines with the HOMERUN model using as ionizing continuum a BPASS model with age 10⁶ yr and solar stellar metallicity (Section 3). As mentioned in the previous section, the nitrogen lines ([N II]λ5755, [N II]λ6583) and sulphur lines ([S II]λλ4069,4076, [S II]λλ6716,6731], [S III]λ6312, [S III]λ9531) predicted by the model are rescaled by the factors 𝒩 and 𝒮 which are free parameters of the fit and which take into account that the relative N/O and S/O abundances may be different than those assumed in the photoionization calculations. The same rescaling is not applied to the oxygen lines ([O III]λ5007, [O III]λ4363, [O II]λλ3726,3729, [O II]λλ7320,7330) since they do not depend linearly on abundances. For the grids defined by a specific ionizing continuum and gas metallicity value, we find the weights w_i and 𝒩, 𝒮 factors which minimise the loss function ℒ; finally we select the model with the minimum ℒ across all metallicities considered.

The top left panel of Fig. 2 shows the variation of the loss function as a function of the gas-phase metallicity, that is, after dust depletions have been taken into account; the blue vertical dashed line represents the estimated model metallicity, which trivially corresponds to that of the model minimizing ℒ. For all lines, we have adopted an acceptable discrepancy of 10% (ΔL_l(X⁺ⁱ)_obs/L_l(X⁺ⁱ)_obs = 0.1) when computing ℒ.

The bottom left panel of Figure 2 represents the nitrogen and sulphur metallicities provided as input in the CLOUDY simulations (dashed colored lines) and derived from the fit after correcting the input model values with the 𝒩, 𝒮 factors (filled circles). The values derived from the fit are close to the relative metal abundances assumed in the Cloudy models.

The top right panel of Figure 2 shows the weights of the single-cloud models with different U and N_H which were derived for the best fit: four single clouds out of the 72 in the U, N_H grid have non-zero weights, three of them account for more than 90% of the total Hβ flux. Using the weights, it is possible to compute the average density and ionization parameter of the clouds, which are identified by the black filled circle in top right panel of Figure 2.

While it is not possible to affirm that the single clouds used in the model have a direct correspondence to the real N_H, U distribution, the average density and ionization parameter derived from the fit are representative of the average density and ionization parameter of the HII region.

In the bottom right panel of Figure 2, we compare the model predictions with the observations: the x-axis represents the ratio between model and observed line fluxes for a certain line, which is indicated on the y-axis. The length of the gray bar represents a ±10% discrepancy. In the case of the [O III]4363 line the gray bar reflects the relative error of the observations which is larger than the adopted 10% discrepancy. Finally, the green color indicates a ratio which is within 0.9 and 1.1 (10% discrepancy max), while orange bars (when present) indicate values outside that range.

As discussed in Section 3, a loss function value of 0.11, which is the minimum value in the top left panel, can be interpreted as an average discrepancy of ∼3.3% between observations and models which is visually confirmed in Figure 2.

The model is able to reproduce a large number of emission lines from different ionization stages within 5% discrepancy (with only the [O III]4363 line reproduced with a 10% discrepancy but still consistent with the observational error): this is a definite improvement with respect to previous models. [S III] lines are in very good agreement with observations, at variance with previous work in which [S III] lines are usually so difficult to reproduce that problems with their atomic parameters have been suggested (Mingozzi et al. 2020). The physical parameters of the HII region inferred from the model, such as metallicity 12 + log(O/H) = 8.56, average ionization parameter U = −2.47 and hydrogen density N_H = 1.96, are also realistic and in agreement with our expectations for HII regions. Finally, the plot of the loss function in Figure 2 presents a minimum so deep it requires a logarithmic scale and indicates that the problem is well defined and not affected by degeneracies, as one might naively expect given the freedom in the choice of the weights. In other words, this model has potentially many parameters (the weights of each single CLOUDY model) but it is built in such a way that if the metallicity is not the right one, the lines cannot be well reproduced.

4. The sample of HII regions by Zurita et al.

After presenting the model and its characteristics with a few selected examples, we present its application to larger samples of emission-line sources. In this paper, we focus on HII regions within the Milky Way or local spiral galaxies.

We consider the sample by Zurita et al. (2021), who collected the emission line fluxes of 2831 published HII regions from 51 nearby galaxies to study how the presence of bars affects metallicity gradients. The compiled line fluxes, reddening corrected and normalized to Hβ, include both bright emission lines ([O II]λλ3726,3729, [O III]λλ4959,5007, [N II]λλ6548,6583, Hα, [S II]λλ6717,6731, [S III]λλ9069,9532) and auroral lines ([S II]λλ4068,4076, [O III] λ4363, [N II]λ5755, [S III] λ6312, and [O II]λλ7320, 7330) when available.

For our analysis we selected only the sources where the fluxes of all bright and auroral emission lines are measured with S/N >  5, resulting in a sub sample of 24 HII regions. The choice of S/N threshold for auroral lines aims to reduce uncertainty in metallicity determination. Section 6 will demonstrate how relaxing this constraint can lead to less precise metallicities, the degree of which depends on the accuracy of observed lines. This is crucial when analyzing galaxy spectra with potentially lower signal-to-noise ratios. We will explore this further in a forthcoming paper dedicated to measuring galaxy metallicities (Amiri et al., in preparation).

We performed model fitting with both dusty and non-dusty clouds, that is, using a weigthed combination of single-cloud photoionization models either with or without dust mixed to the gas – and with all the ionizing continua presented in Section 3. Each model set M is then characterized by the ionizing continuum and by whether there is dust mixed with the gas, and consists of U, N_H grids of single-cloud models, each computed for a given metallicity Z. For each model set M, we then find the minimum value of ℒ by varying Z, which we label ℒ_M, min. We then consider ℒ_min as the smallest of all ℒ_M, min values, selected by having scaling factors for N and S abundances within ±0.3 and ±0.15 dex, respectively. Finally, we consider acceptable all the models for which

$$\begin{array}{c}\hfill \mathcal{L}\le {\mathcal{L}}_{min}+0.25\end{array}$$(13)

Following equation (11), the increase on average “discrepancy” between model and observations depends on ε_a and on ε_r of the best model with ℒ = ℒ_min. For instance, given ε_a = 0.1, models with ℒ = ℒ_min + 0.25 have real “discrepancies” ε_r = 0.05, 0.07, 0.11 if the best models have ε_r = 0.01, 0.05, 0.1, respectively. The choice of the ℒ_min + 0.25 threshold is arbitrary and does not affect the best estimate of metallicity determined by the models, as this depends solely on ℒ_min. Increasing the acceptability threshold only widens the range of model metallicities considered consistent with the data.

4.1. Observed versus model lines

We now compare models and observations by comparing the observed and model line fluxes. The left panel of Figure 3 shows the distribution of the ℒ_min values for all the HII regions in the selected subsample. We distinguish the cases where, for each HII region, ℒ_min is found considering all models, only models with dust, and those without. When considering models with or without dust, we find ℒ_min < 1.1 except for three cases each where ℒ_min is ∼2 or larger. However, when finding ℒ_min over all models we always have ℒ_min < 1.1 and the best-fit models are given in 20 out of 24 sources by dusty models, while in the rest by models without dust. Indeed, the left panel of Figure 3 shows that, overall, dusty models provide in general lower ℒ_min values. The distributions of ℒ_min values thus imply that our multi-cloud models are able to reproduce all emission lines with an accuracy better than 10%; this is an extremely good achievement very rarely matched in previous works (see, e.g., Chevallard et al. 2018, Blanc et al. 2015, Morisset et al. 2016, D’Agostino et al. 2019, Mingozzi et al. 2020, Olivier et al. 2022, Fernández et al. 2022). It is important to mention that only a small fraction of the single-cloud models in a 8 × 9 (U, N_H)-grid have non-zero weights: in 21 out of 24 best-fit models, the numbers of clouds with non-zero weights is between four and nine, with only three models being made of 11 or 12 clouds.

Fig. 3. Model fit of the HII regions from the Zurita et al. (2021) sample. Left: distribution of ℒmin values for the subsample of HII regions selected from Zurita et al. (2021). Right: distributions of the ratios between the model and observed line fluxes for all the selected HII regions and all the models with ℒ ≤ ℒmin + 0.25. The numbers in the top-right corner of the panels represent, from top to bottom, the median and standard deviation of the distribution of ℒ values. The vertical dashed lines represent the ±10% discrepancy and show how the large majority of models reproduce all lines within ±10% of the observed values.

The performance of the model with each line is shown in the right panel of Figure 3 where we present the distributions of the ratios between observed and model line fluxes for all the selected HII regions and all the models with ℒ ≤ ℒ_min + 0.25. The HOMERUN model is able to reproduce to within 10% both strong and auroral lines, from high ([O III]), intermediate ([S III]) and low excitation ([O II], [S II], [N II]) lines. In particular, the model can reproduce both [S II] and [S III] lines at variance with what is found in other studies as already mention in Sec. 4.1. In our case [S III] lines are those for which we observe the largest discrepancies which are still below 20% for [S III]λ9531 in all models while only 11% of the models present discrepancies larger than 20% for [S III]λ6312.

4.2. Observed versus model T_e metallicities

We now compare models and observations by comparing the T_e-metallicities obtained by both observed and model line fluxes. We estimated T_e-metallicities using the python package pyneb, which computes the physical conditions and ionic and elemental abundances of photoionized gas by solving the equilibrium equations for an n-level atom (Luridiana et al. 2015). We selected the same sets of atomic parameters used in CLOUDY and we considered a two-zone cloud where the density of both zones is given by N_e(S⁺), the temperature of the low-ionization zone is given by T_e(N⁺) and that of the high-ionization zone by T_e(O⁺³). The choice of T_e(N⁺) for the temperature of the low-ionization zone is justified by the fact that T_e(O⁺) can lead to underestimated T_e-metallicities and is the same adopted by Zurita et al. (2021) (see also Pérez-Montero & Díaz 2003; Méndez-Delgado et al. 2023a). In particular, we adopt the following procedure.

T_e(O⁺³) and N_e(S⁺) were computed simultaneously by combining the [O III] λ5007/λ4363 and [S II] λ6731/λ6716 ratios with the getCrossTemDen procedure of pyneb. The [S II] flux ratio is not an accurate tracer of density below the low-density limit (less than a few ×10 cm⁻³) which characterises most HII regions (e.g., Nagao et al. 2006a); however, the temperature-sensitive line ratios that we used are almost independent of density in this regime.

We used N_e(S⁺) to compute T_e(N⁺) from [N II]λ5755/λ6584, T_e(S⁺) from [S II] λ4069/λ6720, T_e(S⁺²) from [S III]λ6312/λ9531 and T_e(O⁺) from [O II]λ7325/λ3727, but we used T_e(N⁺) to compute N_e(O⁺) from [O II]λ3726/λ3729. The [O⁺/H] and [O⁺²/H] abundance ratios were computed from [O II]λ3727/Hβ using T_e(N⁺), N_e(S⁺) and from [O III] λ5007/Hβ line ratios using T_e(O⁺²), N_e(S⁺); [O/H] was then the sum of the abundances of the two ions, neglecting any contribution from either neutral oxygen or O⁺³, which are typically negligible for HII regions (e.g., Zurita et al. 2021). The [N/O] and [S/O] abundance ratios were similarly computed from [N II]λ6584/[O II]λ3727 and [S II]λ6720/[O II]λ3727 line ratios. Here we assumed that in the regions emitting these lines, N, S and O are mostly single ionized ions.

Errors are estimated with 100 Monte Carlo realizations where the observed line sets are replicated by adding to line fluxes random numbers extracted from a Gaussian distribution with 0 mean and standard deviation equal to the observed relative error (i.e., equal to the assumed discrepancy of ε = 10%). Figure 4 (left) compares the T_e metallicities derived from the observed line fluxes to the ones predicted by the models: considering all acceptable models (ℒ < ℒ_min + 0.25) for all objects, the agreement is good as model and observed metallicities have a median difference of −0.03 dex with a 0.07 dex standard deviation. Similarly, as shown in Figure 4 (right), when considering only the best models (ℒ = ℒ_min) the median difference is −0.03 dex with a 0.06 dex standard deviation. These values are consistent with the < 10% discrepancy between model and observed line flux that we required for the HOMERUN model and provide a different measure of the (dis)agreement between models and observations. Fig. 3 of Chevallard et al. (2018) shows a similar agreement between O/H from Te and (single cloud) models.

Fig. 4. Comparison between the oxygen abundances derived with the Te method using observed and HOMERUN model line fluxes from the selected HII regions from the sample of Zurita et al. (2021). In the left panel, we consider, for all sources, all the models with ℒ < ℒmin + 0.25; the vertical dashed line represents the median of the distribution. In the right panel, we consider only the models with ℒ = ℒmin.

Figure 5 presents the comparison between the metallicities used as input for the HOMERUN model calculations, labeled as “log(O/H)_HR [in models]”, and the T_e metallicities derived with pyneb from the line fluxes predicted by the models, labels as “log(O/H)_HR [T_e method]”. This is not a comparison between models and observations, as it only involves model quantities, but provides an indication on how (non-)accurate the T_e method is in retrieving the metallicity of a system. The T_e metallicities are always smaller than the true values and can be discrepant up to ∼0.5 dex: the median discrepancy and standard deviation of the difference between T_e and model metallicities considering only the best models are −0.09 and 0.07 dex, respectively. These values increase to −0.13 and 0.14 dex when considering all models with ℒ < ℒ_min + 0.25. We believe that this discrepancy is due to the simplified assumptions underlying the T_e method (e.g., 2-3 zones with constant density and temperature), confirming the findings by other authors (see, e.g., Esteban et al. 2014; Cameron et al. 2022), and showing that, with the computational capabilities and the quality and amount of data now available, it is important to adopt more accurate methods for abundance determinations from emission line ratios. A companion paper by Amiri et al. (in preparation) will address in more details the accuracy problems of the T_e method.

Fig. 5. Comparison between the model metallicities inferred from the multi-cloud modeling of the Zurita et al. (2021) sample (x axis) and the Te metallicities derived from the line fluxes predicted by the same models (y models). The red dots represent the best models for each HII region of the sample, while the green dots represent all “acceptable” models, that is, those with ℒ < ℒmin + 0.25. The solid blue line represents the 1:1 relation, while the dotted blue line represents the case where Te metallicities are 0.1 dex lower than model metallicities.

4.3. The effect of varying N/O and S/O relative abundances

Figure 6 outlines the importance of allowing for a scaling factor of the N and S abundances. The left panel shows how in the example model the minimum value of the ℒ function increases by a factor ≃10 when no abundance scaling is taken into account. Overall, when fitting all the HII regions in the Zurita et al. (2021) sample (see Section 4), minimum ℒ values increase by factors up to over 100 (center/middle panel), while, at the same time, metallicities are over- or under-estimated by up to 0.4 dex in case of no abundance scaling (right panel). Clearly, rescaling the N and S abundances with respect to the ones assumed in the photoionization computations is extremely important to improve the comparison of models and observations and the reliability of the metallicity estimates (e.g., Strom et al. 2018).

Fig. 6. Effect of rescaling N/O and S/O relative abundances with respect to the values used in photoionization calculations. Left: variation of the loss function ℒ as a function of metallicity when fitting the emission line spectrum of the test HII region in Figure 1. The blue line with dots represents the HOMERUN model with N and S rescaling (standard model), while the red line with dots represents the same model but without rescaling (i.e., no N,S abundance scaling). We note that the minimum value of the loss function when excluding rescaling (represented by the dotted red line) is around ten times higher than that for the standard model (represented by the dotted blue line). The vertical red dashed line represents the metallicity of the models with the minimum value of the loss function ℒ. Middle: ratios of ℒ minimum values from models with no N,S abundance scaling and standard models when applied to the HII regions of the Zurita et al. (2021) sample. The effect of non-rescaling N and S is to increase the minimum ℒ value by up to a factor 100. Right: variation of abundances of models with no N,S abundance scaling with respect to standard models compared to the abundances of standard models.

4.4. Comparison with constant-pressure single-cloud models

Figure 7 shows the large improvement of our multi-cloud approach with respect to the case of the single-cloud constant pressure models which are used in the literature by many authors (e.g., Blanc et al. 2015; Dopita et al. 2016; Pérez-Montero et al. 2021 and references therein).

Fig. 7. Comparison of HOMERUN with single-cloud constant-pressure models. Left: variation of the loss function ℒ as a function of metallicity for the emission line spectrum of the test HII region. The blue line with dots represents the HOMERUN model presented here (standard model), while the red line with dots represents the case of constant pressure single-cloud models (Constant-P single cloud). The vertical red and blue dashed lines represent the metallicity of the models with the minimum value of the loss function ℒ. Middle: ratios of ℒ minimum values from Constant-P single-cloud and HOMERUN standard models. Right: Variation of abundances of Constant-P single-cloud models with respect to HOMERUN standard models.

The left panel compares the ℒ curves obtained with the single (red) and multi-cloud (blue) approach applied to the test HII region from Zurita et al. (2021) used before. The blue curve is the same as in Figures 2 and 6.

In order to obtain the orange curve, we have obtained grids of constant-pressure, single-cloud models following the same approach as before: for each given ionising continuum and gas metallicity, we have computed a grid of models with varying U and N_H. For each single-cloud model in the grid we have explored the same (U, N_H) grid used in the multi-cloud model to determine the optimal (U, N_H) pair that minimises the loss function, along with the best scalings for N and S abundances, at given gas metallicity and ionising continuum. We have then considered the single-cloud model with the lowest ℒ value. In practice, we have selected the single-cloud, constant-pressure model which provides the minimum value of the loss function for a given continuum, gas abundance, ionization parameter and density. The difference with respect to other models in the literature is that we allow for a variation of N/O and S/O and/or we consider a much wider grid in density (e.g., van Zee et al. 1998; Dopita et al. 2013; Blanc et al. 2015; Strom et al. 2018; Yan 2018).

In the case of this particular example, the minimum ℒ value of our multi-cloud approach is more than a factor ten lower although the estimated gas metallicity is similar. The comparison between the orange curves in Figures 2 and 6 shows that rescaling N and S abundances is crucial for having a low ℒ minimum value. Indeed, the single-cloud constant pressure model has a similar ℒ_min value as the multi-cloud model where the N and S abundances are those assumed in the photoionization calculations.

The center and right panels in the figure compare the minimum ℒ values and the oxygen abundances obtained by applying the two methods to the HII regions of the Zurita et al. (2021) sample: ℒ values obtained with the multi-cloud approach are more than a factor 10 lower, but the discrepancy can increase by two orders of magnitude. On the other hand, discrepancies in oxygen abundances are less severe and are up to ±0.3 − 0.4 dex.

Finally, we note that the difference between the ionization parameter of the single-cloud, constant-pressure models and the average U of the HOMERUN models is on average −0.2 dex with a standard deviation of 0.5 dex. Similarly, the differences in densities are on average 0.2 dex but with a standard deviation of 1.2 dex.

5. The sample of HII regions from the CHAOS project

We now apply the model to a sample of HII regions with more homogeneous data, higher S/N and a larger number of emission lines.

We consider the HII regions of six galaxies from the CHemical Abundances of Spirals (CHAOS) project (Berg et al. 2015; Croxall et al. 2015, 2016; Berg et al. 2020; Rogers et al. 2021, 2022) and we apply the same procedure followed in the previous section. At variance with the compilation by Zurita et al. 2021, the data now include 11 additional lines (see Figure 8 for the list) including [O I]λ6300, a line which has always been difficult to account for in photoionization models (e.g., Dopita et al. 2013). We select the HII regions with the same criterion used for the Zurita et al. (2021) sample, that is, the auroral lines detected with S/N >  5.

Fig. 8. Model fit of the HII regions from the CHAOS sample. Left: distribution of ℒ values for the best or acceptable models of the HII regions from the CHAOS sample. Right: distributions of the contributions to ℒ from each of the lines used for model fitting. The adopted error on the line fluxes is the maximum between the observed error and 10%.

Figures 8 and 9 confirm the results of the previous section: the ℒ_min values of the best models are all lower than 1 (red bars in figure) indicating an average agreement between model and observed line fluxes better than 10%. All lines are well reproduced and in particular, the models are able to reproduce at the same time [O I], [O II] and [O III], as well as [S II] and [S III], which is a significant improvement compared to previous works (see, e.g., Mingozzi et al. 2020 on [S III] and Dopita et al. 2013 on [O I]). We note that the [NeIII]λ3869 line does not provide any constraint to the model since it is the only Neon line, and it is used only to estimate Ne abundance. It is remarkable that when using only the emission lines from the Zurita et al. (2021) sample and [O I]λ6300, the ℒ_min values are all lower than 0.3. Similarly to the case of the Zurita et al. (2021) sample, only a small fraction of the single-cloud models in a 8 × 9 (U, N_H)-grid have non-zero weights: in 23 out of 24 best fit models, the numbers of clouds with non-zero weights is between five and ten, with only one model being made of 12 clouds.

Fig. 9. Comparison between the oxygen abundances derived with the Te method using observed and HOMERUN model line fluxes from the selected HII regions from the CHAOS sample. In the left panel, we consider, for all sources, all the models with ℒ < ℒmin + 0.25; the vertical dashed line represents the median of the distribution. In the right panel, we consider only the models with ℒ = ℒmin.

We have estimated T_e metallicities as in the previous section and found that the agreement between estimates from observed and model line fluxes is similar to what found for the Zurita et al. (2021) sample (Figure 9). Finally, we confirm as well the discrepancy of up to 0.2 dex between the metallicities used in the photoionization model computations and the T_e metallicities derived from the line fluxes predicted by the same models (Figure 10).

Fig. 10. Comparison between the model metallicities inferred from the multi-cloud modeling of the CHAOS sample (x axis) and the Te metallicities derived from the line fluxes predicted by the same models (y models). The red dots represent the best models for each HII region of the sample, while the green dots represent all “acceptable” models, that is those with ℒ < ℒmin + 0.25. The solid blue line represents the 1:1 relation, while the dotted blue line represents the case where Te metallicities are 0.1 dex lower than model metallicities.

Outside of H II regions, the warm (T_e ∼ 10⁴ K), ionized component of the ISM is made of a more diffuse medium, the so-called DIG, particularly evident above and below the galactic plane. Such a medium can significantly affect the integrated spectra of H II regions and star-forming galaxies, depending on the size of the aperture from which spectra are extracted. Mannucci et al. (2021) investigated the contribution of the DIG which is often invoked to explain the differences between the spectra of HII regions in local galaxies and the total or kpc-scale spectra of star-forming galaxies when the flux ratios involving low-ionization emission lines are considered (e.g., Belfiore et al. 2016; Kumari et al. 2019; Vale Asari et al. 2019; Della Bruna et al. 2020; Sanders et al. 2020). Comparing single slit spectra of HII regions from CHAOS with those extracted with different aperture sizes from integral field spectra, they concluded that these differences are mostly ascribed to the small anglular size of the slits, which is not covering the full extent of the HII region, rather than to the DIG contamination. In principle, our model approach, based on combining multiple clouds to reproduce the correct volume of gas occupied by any given ion X⁺ⁱ with the correct T_e and N_e, can naturally take into account both aperture effects and the contributions from regions with different physical properties like the DIG. However, this issue will be addressed in a forthcoming paper.

6. Gas and stellar metallicities in the Milky Way

We present the last application of our method and its potential to solve the observed discrepancy between the metallicities of HII regions and those of the stars in our Galaxy.

The metallicity of the ionized gas is generally in good agreement with that of the stars, at least in external galaxies (see, e.g., Fig. 11 of Bresolin et al. 2016). However, there is a discrepancy between the metallicities of HII regions in our Galaxy, derived from the application of the T_e method to gas emission lines, and the stellar metallicities derived from B-type stars, Classical Cepheids and Open clusters, which are ∼0.2 dex larger (e.g., Fig. 4 of Esteban et al. 2022). Such a difference is not present when considering gas-phase metallicities from recombination lines and is part of the so-called abundance discrepancy problem which, as described in Sec. 1, is quite a common occurrence when comparing T_e and recombination-line-based metallicities. As mentioned by Esteban et al. (2022), such discrepancies might be related to temperature fluctuations which bias auroral lines toward the hotter parts of the HII regions providing higher than average temperatures and consequently lower abundances (Peimbert 1967). Indeed, this discrepancy disappears when correction for temperature fluctuations and ionization fractions are applied (Méndez-Delgado et al. 2022). Still corrections for temperature fluctuations are difficult to estimate and ionization corrections must be derived from photoionization models, thus adding up additional assumptions and uncertainties to the T_e method.

We consider all the Galactic HII regions from table 1 of (Méndez-Delgado et al. 2022) which have at least three of the five auroral lines considered so far detected with S/N >  5 and estimate their metallicities by applying our method to the emission lines fluxes from Arellano-Córdova et al. (2020, 2021), and Zurita et al. (2021). Figure 11 shows the HII region gas-phase metallicities plotted as a function of GAIA-based galactocentric distances from Méndez-Delgado et al. (2022), and compare them with the metallicity gradient of the same HII regions estimated with the T_e method without any ionization correction (Arellano-Córdova et al. 2021), the Cepheid metallicity gradient from Mollá et al. (2019) and the intermediate-age (1–3 Gyr) open clusters metallicity gradient from Magrini et al. (2023). The linear fit considers all HII regions and was performed with the ltsfit⁴ python package which implements a Robust linear regression with scatter in one or two dimensions (Cappellari et al. 2013). Unlike Sections 4 and 5, where we required HII regions to have five auroral lines with S/N >  5, we have reduced the selection threshold to three auroral lines with S/N >  5 to expand our sample size. This change, as the figure demonstrates, does not significantly alter our conclusions but may slightly increase the scatter. The plot clearly shows how our model-based metallicities are, on average, in agreement with the stellar metallicities especially when considering that oxygen is a refractory element and is likely depleted into stars: indeed, it is important to recall that these are gas-phase metallicities which do not take into account the depletion of metals into dust. Oxygen depletion can vary between 10–20% (Psaradaki et al. 2023; Amayo et al. 2021) and the 40% assumed by CLOUDY. Therefore, any abundance discrepancy here can be easily explained with a ≲0.1 dex depletion of oxygen into dust without having to rely on uncertain temperature fluctuations or ionization corrections.

Fig. 11. Comparison between our metallicity estimates of HII regions in the Galaxy, the metallicity gradient from the Te abundances of HII regions (Arellano-Córdova et al. 2021), the Cepheids (Mollá et al. 2019), and intermediate age open clusters (Magrini et al. 2023). The blue line surrounded by the blue shaded area represents the linear fit, with uncertainties, to all our metallicity estimates.

7. Summary and conclusions

We have presented a new approach to modeling emission lines from photoionized gas that allows the observed line ratios from a wide range of ionization to be reproduced with an accuracy to better than 10% and whose main application is the accurate determination of gas metallicity.

Our approach is based on the weighted combination of multiple single-cloud photoionization models computed with CLOUDY, as in previous works. The novelty of our approach is that the weights of the single-cloud models are not parametric nor assumed a priori but are free parameters of the fit. The resulting best-fit model accounts for the varying physical conditions (temperatures, densities, ionization conditions, and gas metallicities) observed along the line of sight of the interstellar gas of HII regions and galaxies.

Since observed line luminosities are integrated quantities over the entire volume of an emitting source, the single-cloud models contributing to a best fit do not necessarily correspond to physical entities but are the basic building blocks of the model. Therefore, if a best-fit model has a single cloud with specific U, N_H, this does not necessarily mean that such a cloud exists but that the total line emission has a contribution from gas with those physical parameters, regardless of where and how it is spatially distributed.

Our model has as many free parameters as the single-cloud models considered. However, the fitting of the observed data is not degenerate, and the loss function (ℒ) curves, computed as a function of gas metallicity, have very deep and well-defined minima. Indeed our approach is similar to the one used for fitting the stellar continua of galaxies, where many hundreds of templates can be combined with weights that are free parameters of the fit (see, e.g., the ppxf code to analyze stellar continuum spectra).

We have shown that a critical point for our modeling is to allow for different abundances for elements, such as N and S, from those used in the photoionization calculations. The N and S emission lines are not the major coolants of the photoionized region, such as O lines, and therefore, their luminosities scale almost linearly with gas abundances. Not accounting for N and S abundance variations can significantly affect the quality of the fit (with minimum ℒ values being worse by even a couple of orders of magnitude) and the estimated gas metallicity.

We have also shown that our approach provides a significant improvement compared to the single-cloud, constant-pressure models commonly used in the literature and that it is able to model a larger set of emission lines. Indeed, several works based on single-cloud photoionization modeling claimed to pin down the metallicity with high accuracy by using just a small set of lines, typically the optical strong lines, and they rarely, if ever, include auroral lines, which are instead successfully reproduced in HOMERUN.

We have applied our model to the HII regions from the samples by Zurita et al. (2021) and from the CHAOS project. For the sample by Zurita et al. (2021), we have considered HII regions where all the following lines were detected with an S/N greater than five: [O II]λλ3726,3729; [O III]λλ4959,5007; [N II]λλ6548,6583; Hα, [S II]λλ6717,6731; [S III]λλ9069,9532; [S II] λλ4068,4076; [O III]λ4363; [N II]λ5755; [S III]λ6312; and [O II]λλ7320, 7330. All the observed emission lines for the 24 selected HII regions were reproduced by the model to better than 10%. In particular, the model was also able to reproduce the [S III]9069,9532 lines, which are well known to be problematic. The T_e-metallicities derived from the observed emission lines and those derived from the model emission lines agree to within 0.05 dex, another indication of the extremely good performance of the models.

However, we showed that T_e-metallicities suffer from uncertainties related to the simplified assumptions underlying the method. The T_e-metallicities computed from the emission lines predicted by the models can be more than 0.1 dex smaller than the real gas metallicities of the models themselves.

The CHAOS sample is more homogeneous in terms of data quality, characterized by a better S/N, and provides the fluxes of many more emission lines. Overall, we have modeled the same emission lines used for the Zurita et al. (2021) sample, but we have also added other lines, doubling their total number.

Even with these increased constraints, the models reproduce the observed emission lines to better than 10%. It is worth noting that the HOMERUN models are able to simultaneously reproduce the emission lines from very different ionization stages, including [O I], [O II] and [O III].

Finally, we have shown that the gas metallicities estimated with our models in HII regions of the Milky Way are in better agreement with the stellar metallicities than the estimates based on the T_e method.

Overall, our method provides a new accurate tool to estimate the metallicity and the physical conditions of the ionized gas. It can be applied in many different science cases and wherever there are emission lines from photoionized gas.

There are many applications of our method, and they will be presented in a series of forthcoming papers. In particular, the second paper of this series (Amiri et al., in preparation) will present the application of this method to the spectra of galaxies, determine their metallicities, and provide a new self-consistent calibration of the strong line method. Future studies will explore the limitations and informational content revealed by incorporating ultraviolet, mid-infrared, and far-infrared spectral lines when modeling spectra of HII regions, star-forming galaxies, and active galactic nuclei.

Acknowledgments

We are grateful to Gary Ferland and his collaborators for the development of CLOUDY and for allowing its free use for the community. This work has been partially financed by the European Union with the Next Generation EU plan, through PRIN-MUR project “PROMETEUS” (202223XPZM). We acknowledge the support of the INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities” and INAF Large Grant 2022 “Dual and binary SMBH in the multi-messenger era”. SC and GV acknowledge support by European Union’s HE ERC Starting Grant No. 101040227 – WINGS.

References

All Figures

Fig. 1. Schematic representation of the fitting procedure. For any U, NH grid of single-cloud models, we computed the weighted combination of clouds in the grid that gives the best match with the observed lines. Usually, only a small fraction of the clouds in the grid have non-zero weights. The goodness of the fit was then quantified with a value of the loss function ℒ. The procedure was repeated for all U, NH grids computed with different gas metallicities, and the behavior of ℒ as a function of metallicity was obtained. The metallicity of the observed spectrum is identified by the minimum value of ℒ. This procedure can be repeated with different ionizing continua, thus estimating the effect of varying the ionizing continuum on the metallicity estimate.

In the text

Fig. 2. Example of model fit to the emission line spectrum of an H II region. Top-left panel: variation of the loss function ℒ as a function of the oxygen abundance, 12 + log(O/H), for the emission line spectrum of the test HII region. In this scale, the solar metallicity is 8.69. The blue vertical dashed line represents the metallicity of the model with the minimum value of ℒ. Bottom-left panel: variation of the N and S abundances as a function of the oxygen one. The dashed lines represent the values adopted for the computation of the photoionization models, while the filled circles correspond to the abundances rescaled as described in the text. Top-right panel: grid of single-cloud models in log(U) and log(NH). The colors represent the weights of each single-cloud model, as indicated by the color bar, when ℒ reaches the minimum value. In this case only four out of 72 models of the grid have non-zero weights. The black dot represents the weighted density and ionization parameter of the single-cloud models. Bottom-right panel: ratio between the model and observed line luminosities. The gray areas represent the acceptable discrepancy of ±10% or, when larger, the relative error on the observed value.

In the text

Fig. 3. Model fit of the HII regions from the Zurita et al. (2021) sample. Left: distribution of ℒmin values for the subsample of HII regions selected from Zurita et al. (2021). Right: distributions of the ratios between the model and observed line fluxes for all the selected HII regions and all the models with ℒ ≤ ℒmin + 0.25. The numbers in the top-right corner of the panels represent, from top to bottom, the median and standard deviation of the distribution of ℒ values. The vertical dashed lines represent the ±10% discrepancy and show how the large majority of models reproduce all lines within ±10% of the observed values.

In the text

Fig. 4. Comparison between the oxygen abundances derived with the Te method using observed and HOMERUN model line fluxes from the selected HII regions from the sample of Zurita et al. (2021). In the left panel, we consider, for all sources, all the models with ℒ < ℒmin + 0.25; the vertical dashed line represents the median of the distribution. In the right panel, we consider only the models with ℒ = ℒmin.

In the text

Fig. 5. Comparison between the model metallicities inferred from the multi-cloud modeling of the Zurita et al. (2021) sample (x axis) and the Te metallicities derived from the line fluxes predicted by the same models (y models). The red dots represent the best models for each HII region of the sample, while the green dots represent all “acceptable” models, that is, those with ℒ < ℒmin + 0.25. The solid blue line represents the 1:1 relation, while the dotted blue line represents the case where Te metallicities are 0.1 dex lower than model metallicities.

In the text

Fig. 6. Effect of rescaling N/O and S/O relative abundances with respect to the values used in photoionization calculations. Left: variation of the loss function ℒ as a function of metallicity when fitting the emission line spectrum of the test HII region in Figure 1. The blue line with dots represents the HOMERUN model with N and S rescaling (standard model), while the red line with dots represents the same model but without rescaling (i.e., no N,S abundance scaling). We note that the minimum value of the loss function when excluding rescaling (represented by the dotted red line) is around ten times higher than that for the standard model (represented by the dotted blue line). The vertical red dashed line represents the metallicity of the models with the minimum value of the loss function ℒ. Middle: ratios of ℒ minimum values from models with no N,S abundance scaling and standard models when applied to the HII regions of the Zurita et al. (2021) sample. The effect of non-rescaling N and S is to increase the minimum ℒ value by up to a factor 100. Right: variation of abundances of models with no N,S abundance scaling with respect to standard models compared to the abundances of standard models.

In the text

Fig. 7. Comparison of HOMERUN with single-cloud constant-pressure models. Left: variation of the loss function ℒ as a function of metallicity for the emission line spectrum of the test HII region. The blue line with dots represents the HOMERUN model presented here (standard model), while the red line with dots represents the case of constant pressure single-cloud models (Constant-P single cloud). The vertical red and blue dashed lines represent the metallicity of the models with the minimum value of the loss function ℒ. Middle: ratios of ℒ minimum values from Constant-P single-cloud and HOMERUN standard models. Right: Variation of abundances of Constant-P single-cloud models with respect to HOMERUN standard models.

In the text

	Fig. 8. Model fit of the HII regions from the CHAOS sample. Left: distribution of ℒ values for the best or acceptable models of the HII regions from the CHAOS sample. Right: distributions of the contributions to ℒ from each of the lines used for model fitting. The adopted error on the line fluxes is the maximum between the observed error and 10%.
In the text

	Fig. 9. Comparison between the oxygen abundances derived with the Te method using observed and HOMERUN model line fluxes from the selected HII regions from the CHAOS sample. In the left panel, we consider, for all sources, all the models with ℒ < ℒmin + 0.25; the vertical dashed line represents the median of the distribution. In the right panel, we consider only the models with ℒ = ℒmin.
In the text

Fig. 10. Comparison between the model metallicities inferred from the multi-cloud modeling of the CHAOS sample (x axis) and the Te metallicities derived from the line fluxes predicted by the same models (y models). The red dots represent the best models for each HII region of the sample, while the green dots represent all “acceptable” models, that is those with ℒ < ℒmin + 0.25. The solid blue line represents the 1:1 relation, while the dotted blue line represents the case where Te metallicities are 0.1 dex lower than model metallicities.

In the text

Fig. 11. Comparison between our metallicity estimates of HII regions in the Galaxy, the metallicity gradient from the Te abundances of HII regions (Arellano-Córdova et al. 2021), the Cepheids (Mollá et al. 2019), and intermediate age open clusters (Magrini et al. 2023). The blue line surrounded by the blue shaded area represents the linear fit, with uncertainties, to all our metallicity estimates.

In the text