Ionization correction factors and dust depletion patterns in giant H ii regions

To the initial grid of models we applied a few filters in order to end up with a subgrid that is representative of our sample of observed H ii regions.

1. 1.

   We selected models with starburst ages lower than 6 Myr. The reason is that at this age the available ionizing photons from the burst have decreased up to a factor of $\sim$10 (Mollá et al., 2009) which in general makes the H ii region too faint to properly observe the weak emission lines. ²²2Models of stellar populations including binaries stars such as those of Eldridge et al. (2017) actually extend the period of time when massive stars ionize the surrounding gas with respect to the stellar populations used here. However, the age limit in our paper is not to be taken literally since what primarily matters is the strength of the ionizing radiation. Of course, the details of the SED matter also but, since we cover a range of ages, we cover by the same token a range of SEDs as well, although this aspect would require more dedication in the future.

2. 2.

   We selected only those density-bounded models obtained by clipping the radiation-bounded models at about 70% of their H$\beta$ intensity. Models more optically thin than this were removed since they likely do not represent the bulk of observations considered in this kind of study. All the radiation bounded models that satisfy filter (i) are included.

3. 3.

   Previous studies of giant H ii  regions have shown that they follow some gross relations between O/H and N/O (e.g., Pilyugin et al., 2012), U and O/H (e.g., Pérez-Montero et al., 2007) and that they form a thin sequence in the BPT diagram (e.g., McCall et al., 1985). We therefore implemented criteria inspired by those of VA16 to define their ‘fake observational sample’ but slightly modified to better represent our observational sample. The expressions we have used are the following:

   $$\begin{split}&\rm log(N/O)=-1.093,\mbox{ for }Z\leq 7.93,\\ &\rm log(N/O)=1.489Z-12.896,\mbox{ for }Z\geq 7.93,\\ &\rm log(N/O)=-1.693,\mbox{ for }Z\leq 8.25,\\ &\rm log(N/O)=1.489Z-14.900,\mbox{ for }Z\geq 8.25,\end{split}$$

   (1)

   $$\begin{split}&\rm log(U)=5-1.25Z,\\ &\rm log(U)=10-1.25Z,\end{split}$$

   (2)

   where $Z$ = 12 + $\log$(O/H) and,

   $$\begin{split}&y_{\rm low}=[a_{\rm low}\tanh{(b_{\rm low})}-c_{\rm low}]-0.70\\ &y_{\rm up}=[a_{\rm up}\tanh{(b_{\rm up})}-c_{\rm up}]+0.10\\ \end{split}$$

   (3)

   where a${}_{\rm low}=-30.79+1.14(y+0.60)+0.27(y+0.60)^{2}$, $b_{\rm low}=5.74(y+0.60)$, $c_{\rm low}=c_{\rm up}=31.09$, $a_{\rm up}=-30.79+1.14(y-0.30)+0.27(y-0.30)^{2}$, $b_{\rm up}=5.74(y-0.30)$, $y=\log$(N ii $\lambda$6284/ H$\alpha$).

Figure 1 shows the models resulting from applying filters i) and ii) (red crosses) and the models obtained after applying also Equations 1–3 of filter iii) (blue circles). Figure 2 shows the BPT diagram for the photoionization models (in colored circles, where the color is related to the oxygen abundance) and the observational sample (the BCG sample in black circles and the GHR sample in black empty squares). The upper panel shows the initial sample of $\sim$31000 photoionization models. The panel on the bottom shows the final sample of 1887 models obtained after applying all our filters. It can be seen that our final sample of models covers successfully the observed H ii regions, and thus seems adequate for the ICFs calculation.

As a last check, we have also compared the final grid of models with our observational sample using line-ratio diagrams which tell about the ionization structure of different elements. We present in Figure 3 the values of $\log$([S iii] $\lambda$9069/[S ii] $\lambda\lambda$6716+31)³³3We do not use the [S iii] $\lambda$9532 line because it is often prone to contamination by sky lines and, although it is intrinsically stronger that the [S iii] $\lambda$9069 line, it may lead to quite erroneous results. The problem is much less severe for [S iii] $\lambda$9069.as a function of $\log($[O iii] $\lambda$5007/[O ii] $\lambda$3727) (left panel) and the values of $\log$([Ar iv] $\lambda$4740/[Ar iii] $\lambda$7135) as a function of $\log($[O iii] $\lambda$5007/[O ii] $\lambda$3727) (right panel).

For the first panel, we can compare our models only with the GHR sample since [S iii] $\lambda$9069 is outside the observed wavelengths in the BCG sample. We see that, in the observed range of [O iii]/[O ii] values, the models do not cover perfectly well the observed [S iii]/[S ii] ratios. A few observations are well below the loci of the models, and the whole distribution of the observed [S iii]/[S ii] values is slightly shifted downwards with respect to the entire subgrid of models, and covers the zone occupied by medium to low oxygen abundances while we know from VA16 (and also Sect. 5 below) that the GHR sample has medium to high oxygen abundance. This can be due to an inadequate description of the SED of the ionizing radiation field, to incorrect values for sulfur dielectronic recombination coefficients, or to the fact that the density distribution of our models is too simple with respect to reality. The last hypothesis is probably the more likely to be correct, but it is out of the scope of the present paper to explore it here. It has been mentioned in Stasińska et al. (2006) and discussed thoroughly in Ramambason et al. (2020). All in all, the situation revealed in this plot is not too bad, and we can expect that our ICFs for sulfur based on [S iii] and [S ii] lines will be reasonable.

For the second panel, we can compare the photoionization models only for the BCG sample since [Ar iv] $\lambda$4740 is absent in our GHR database. Here, we see that the observations occupy a very narrow strip. At a given value of [O iii] /[O ii], the observed values of [Ar iv] $\lambda$4740/[Ar iii] $\lambda$7135 lie at the upper extreme of the range covered by our models. This is reminiscent of the finding by Stasińska et al. (2015), who argued that this might be due to a too soft SED in the stellar atmosphere models. This then indicates that the ICFs we will propose for Ar could give slightly biased abundances, especially at the highest values of [O iii] /[O ii].

In order to obtain more reliable ICFs we have assigned weights to the selected models according to the number of observations they reproduce in the BPT diagram. Specifically, the weights were computed by comparing the position of each model in the [O iii]/H$\beta$ versus [N ii]/H$\alpha$ BPT diagram with the positions of the observed objects.

As said above, we did not use the [O iii] $\lambda$5007 intensities from the BCG sample due to saturation problems. We computed them as $2.98\times I$[O iii] $\lambda 4959$ in order to compare this sample to the models in the BPT diagram. To compute the weights, we first constructed a uniform mesh on the BPT diagram with $\sim$0.3 dex side quadrants (nb: the final results on the ICF do not change significantly when changing the size of the quadrants). Then, we obtained the weights of the models inside each quadrant, W_i, using the following expression$:$

where n_obs(i) is the number of observed objects inside the quadrant and n_mod(i) is the number of photoionization models inside the quadrant $i$. In such a way, models located in a quadrant of the BPT diagram where there are no observations are given only minimal weight in the computation of the ICFs, while models representative of a high number of observed objects are given large weights.

In Figure 4 we show the BPT diagram of the 1887 selected photoionization models (colored circles), the observational sample formed by the BCG and GHR samples (black circles and squares, respectively), and the mesh constructed to compute the weights. The size of the circles represents the weight of each model. The weights range from 0.01 to $\sim$2.0. Models with $\log$[O iii]/H$\beta$ between $-1.2$ and $0.6$ and $\log$[N ii]/H$\alpha$ between $-1.5$ and $0.6$ have in general the highest weights.

The abundance of a particular element X is generally expressed with respect to hydrogen, X/H:

where $\Sigma_{\rm obs.}{X}^{+i}$ represents the sum of the abundances of the observed ions of $X$.

The analytical expression for the ICF in this case may be computed from a fit of the values of:

from the photoionization models. Here $x$ represents the relative ionic fractions of each ion, weighted by the electron density.

One may be interested in the abundance ratio of an element with respect to O, rather than with respect to H. In that case, it may be better to use an ICF computed with such an aim.

To better illustrate this and to explain the notations adopted throughout the paper, we use two examples. The first example concerns nitrogen. The only ion with emission lines in the optical range is N⁺ and thus, the ICF takes into account the contribution of all the other ions present in the nebula. We derived the ICF as$:$

Therefore, N/O is computed as$:$

and N/H = N/O$\times$O/H.

The second example refers to sulfur. We have derived two ICFs, one to be used when only [S ii] lines are observed:

and another to be used when both [S ii] and [S iii] lines are observed:

Therefore, the total S/O ratio is computed as

in the first case, or:

in the second one. The total abundance of sulfur with respect to hydrogen is S/H = S/O$\times$O/H.

For each element, we found that the value of the ICF most of the times depends on the excitation of the object, and sometimes also slightly on the metallicity. Therefore we express the ICFs as a function of O⁺⁺/(O⁺+ O⁺⁺) ⁴⁴4We explored the use of the metallicity as a second parameter, but it turned out inapplicable., which is a quantity easily observed and computed. For simplicity, we will call it $\omega$.

For each abundance ratio for which we compute an ICF we divided the range of $\omega$ values in 10 bins and computed the weighted median (taking into account the weights of the models) in each bin. Then, we performed a fit using these 10 values. We have fitted a fifth-order polynomial expression of the form: $A+B\omega+C\omega^{2}+D\omega^{3}+E\omega^{4}+F\omega^{5}$ for the logarithm of each ICF. As for the uncertainties associated to the ICFs, we computed the weighted 16 and 84 percentiles in each bin and performed a fit to the logarithm of their values using also a fifth-order polynomial expression. We have called the functions representing these fits as $\varepsilon^{+}$ and $\varepsilon^{-}$, respectively.

The analytical expressions of the ICFs derived for C, N, Ne, S, Cl, and Ar, and their associated uncertainties are listed in Table 1. In the following sections we provide some details about each of the elements studied here and we show the figures with the photoionization models and the derived ICFs. A comparison of our ICFs with some of those proposed previously is available in the Appendix A.

In the optical range, the only C line observed is the C ii $\lambda$4267 recombination line. If using this line to compute the C/O ratio, it is necessary to use O ii recombination lines. However, in H ii regions, the C/O ratio has often been derived using the ultraviolet C iii] $\lambda$1909 line, combined with [O iii] $\lambda$1663 or [O iii] $\lambda$5007.

Figure 5 shows the values of ICF(C⁺⁺/O⁺⁺) as a function of $\omega$ for the photoionization models (colored circles). The weighted 16%, 50%, and 84% percentiles for each bin in $\omega$ are shown with black symbols: upward-pointing triangles, circles and downward-pointing triangles, respectively. The fits are shown with solid and dashed lines. Note that, due to the high dispersion of the models, the ICF is less well defined for H ii regions with $\omega<0.05$ or $\omega>0.95$. In addition some of the models with the highest values of $12+\log$(O/H) are the furthest from the fit which indicates that this ICF is more reliable at medium and low metallicities.

The only nitrogen ion that is observed in the optical range is N⁺ (e.g. [N ii] $\lambda$6584, $\lambda$6548).Figure 6 shows the values computed for ICF(N⁺/O⁺) as a function of $\omega$, with the same layout as Fig.5. The figure shows that our proposed ICF is very similar to the traditional one. The median value of the ICFs from our models is larger than the traditional ICF by less than 0.1 dex in the entire range of $\omega$. But the dispersion is not negligible, especially at the highest values of $\omega$. As in the case of carbon, the uncertainty is larger at high metallicity.

The only Ne ion that is observed in the optical spectrum of giant H ii regions is Ne⁺⁺ (whose most intense line is [Ne iii] $\lambda$3869). Figure 7 shows the values of ICF(Ne⁺⁺/O⁺⁺) as a function of $\omega$. It can be seen that the dispersion of our models is large, especially at low values of $\omega$. This is due to two factors. The ionization potential of Ne⁺, 40.96 eV, is actually much higher than that of O⁺, which is 35.12 eV. This induces a significant difference in the photoionization rates of those two ions, especially when the ionizing energy distribution of the stellar radiation field is relatively mild, i.e. at high metallicities. The other factor is the important difference in the charge-transfer recombination rates of O⁺ and Ne⁺, which implies that the relative ionization structures of these elements in H ii  regions also strongly dependent on the ionization parameter. Therefore, Ne/O ratios are very uncertain for low-excitation objects.

There are two ions that may be observed in optical spectra of giant H ii regions: S⁺ (e.g., [S ii] $\lambda$6716, $\lambda$6731 lines) and S⁺⁺ (e.g., [S iii] $\lambda$6312, $\lambda$9069, $\lambda$9532 lines). We study here two cases: 1) when only S⁺ lines are observed and 2) when both S⁺ and S⁺ lines are observed.

The ions that may be observed in the optical range of giant H ii regions are Cl⁺ (e.g., [Cl ii] $\lambda$9123), Cl⁺⁺ (e.g., [Cl iii] $\lambda$5517, $\lambda$5537) and Cl⁺³ (e.g., [Cl iv] $\lambda$8046). However, in this work we only compute an ICF based on Cl⁺⁺ because this ion is very often the only one available in the spectra of giant H ii regions. Figure 10 shows values of ICF(Cl⁺⁺/O⁺⁺) as a function of $\omega$.

The emission lines of argon that are usually observed are: [Ar iii] $\lambda$7135, $\lambda$7751, [Ar iv] $\lambda$4711, 4740. Based on this, we propose two ICFs for argon, one for the case only Ar⁺⁺ emission lines are observed and another one for the case when also Ar⁺³ lines are observed.

First, we have computed the electron temperatures and densities ($T_{\rm e}$ and $n_{\rm e}$) using the version 1.1.10 of PyNeb (Luridiana et al., 2015; Morisset et al., 2020). We used the effective recombination coefficients provided by Storey & Hummer (1995) for H⁺ and the atomic data from Table 2 for the other ions. These atomic datasets were selected following the works of Stasińska et al. (2013); Delgado-Inglada et al. (2016) and Juan de Dios & Rodríguez (2017).

We assumed that each H ii region is characterized by two values of $T_{\rm e}$  and one of $n_{\rm e}$. The temperatures were derived from the [N ii] $\lambda$5755/$\lambda$6584 ratio and the [O iii] $\lambda$4363/$\lambda$5007 or the [O iii] $\lambda$4363/$\lambda$4959 ratio. In the objects where either $T_{\rm e}$([N ii]) or $T_{\rm e}$([O iii]) could not be determined we used the expression given by Garnett (1992): $T_{\rm e}$([N ii]) = 0.70 $\times$ $T_{\rm e}$([O iii]) + 3000 K. This expression is in fair agreement with modern photoionization models (VA16) and with observations (Arellano-Córdova & Rodríguez, 2020).

The value of $n_{\rm e}$ taken for each nebula is the mean of the densities computed with the [O ii] $\lambda$3726/$\lambda$3729 and [S ii] $\lambda$6731/6716 ratios. In the case of the BCG sample, where no density diagnostic was available, a value of 100 cm^-3 was adopted, as suggested by Stasińska (2004).

Second, we have computed the ionic abundances using only the most intense emission lines of each ion. For ions with an ionization potential lower than 20 eV (O⁺/H⁺, N⁺/H⁺, and S⁺/H⁺) we adopted $T_{\rm e}$([N ii]), while for ions with a greater ionization potential (O⁺⁺/H⁺, S⁺/H⁺, S⁺⁺/H⁺, Ar⁺⁺/H⁺, Ar⁺³/H⁺, and Ne⁺⁺/H⁺) we adopted $T_{\rm e}$([O iii]).

Third, we have computed the total abundances of nitrogen, neon, sulfur and argon with the ICFs derived here. For sulfur and argon, we were able to compute the abundances with the two different ICF expressions proposed in this work and compare the results. We did not compute the chlorine abundances due to the lack of relevant intensities in our observational sample.

Uncertainties were computed through a Monte Carlo simulation. First, we generated a Gaussian distribution of 400 values centered in the observed intensity of each line with a $\sigma$ equal to each line uncertainty. The Gaussian distributions may lead to negative values in intensities. We replaced these negative values with an intensity of 0.0001, a value much lower than the minimum intensity observed in our sample⁵⁵5Since we do not present mean values but percentiles, such a procedure is reasonable.. Using the Monte Carlo experiments of each object, we computed the physical conditions, ionic and total abundances. For objects where we used the Garnett (1992) relation to derive either $T_{\rm e}$([N ii]) or $T_{\rm e}$([O iii]), we additionally generated a Gaussian distribution of 400 Monte Carlo experiments with a $\sigma$ equal to 600 K, centered in the nominal temperature derived with this relation. This $\sigma$ corresponds the dispersion found by the ‘fake observational sample’ of VA16 around this relation (see Fig.A2a of VA16).

For the total abundances, we computed two groups of uncertainties: the first one, including the uncertainties associated to the emission lines and to the dispersion in the $T_{\rm e}$([N ii]) vs $T_{\rm e}$([O iii]) relation and a second one including the uncertainties associated to each ICF. To include the latter uncertainty source, we used the analytical expressions from Table 1 and generated a log-normal distribution of 400 values with a confidence interval delimited by the upper and lower uncertainties associated to each ICF. A Gaussian distribution was not used in order to avoid the possibility of obtaining negative ICF values. The log-normal distribution was generated in a way similar to what was done in Garma-Oehmichen et al. (2020), so that its mean ($\mu$) and standard deviation ($\sigma$) were equal to

where $\varepsilon^{-}$ and $\varepsilon^{+}$ are the 16 and 84 weighted percentiles of the ICFs distribution, respectively. These expressions result in a log-normal distribution with the same $\sigma$ than each ICF distribution, but with a slightly different mean value. From this log-normal distribution, we chose a random value and used it to compute the final element abundance.

Figure 14 shows the abundances of oxygen, computed as O⁺/H⁺ + O⁺⁺/H⁺, as a function of $\omega$. There is a clear trend of increasing degree of ionization as O/H decreases. This trend has already been noted in the past (e.g., McGaugh, 1991). It is due to a softer ionizing SED at higher metallicities due to the effect of the metals on the stellar interiors and atmospheres (Pagel et al., 1979; Maeder, 1990) as well as on the ionization parameter through the increasing strength of the winds (Dopita et al., 2006). Since this behavior is related to metallicity, the same trend is expected for other elements.

The upper row of Fig. 13 shows the values of N/O as a function of $\omega$ (left panel) and of O/H (the other two panels). It can be seen that at the highest metallicities, the uncertainties in N/O are much smaller than those in O/H. As already known (c.f. VA16) they anti-correlate with the uncertainties in O/H due to the opposite temperature dependencies of N/O and O/H. At the lowest metallicities, which correspond to high values of $\omega$, the uncertainties in N/O become substantial due the uncertainty in the ICF, and may reach values of up to $\pm\ 0.2$ dex.

The well-known increase of N/O with O/H at high metallicity remains significant despite the uncertainties due to the ICFs. The linear regression gives a slope of $1.20\pm\ 0.23$ for the GHRs. This increase is due to the fact that nitrogen production is both primary and secondary at these metallicities (Mollá & Gavilán, 2010). Primary nitrogen is synthesized from the oxygen and carbon produced before the CNO cycling and thus, is independent of the initial heavy-element abundances while secondary nitrogen is produced during the CNO cycle and depends on the initial abundance of the heavy elements in the star (Henry et al., 2000). As discussed by Pagel (2009), variations in N/O can also occur due to the time-delayed production of N from low- and intermediate-mass stars or due to a preferential loss of oxygen in the galactic winds produced by core-collapse supernovae. For the BCG sample, i.e. at low metallicities, the slope of the N/O versus O/H relation is $0.26\pm\ 0.08$. Some objects, however, are found at quite high values of N/O. Our study shows that uncertainties in the ICF cannot be the cause of this. Izotov et al. (2006) have argued that theoretical models for massive stars including rotation (Maeder & Meynet, 2005) could perhaps explain high values of N/O at low metallicities if taking into account that the density of the N-rich ejecta is larger than that of the H ii regions. This question, however, deserves a more dedicated study.

Neon, argon and sulfur, as well as oxygen are $\alpha$-elements produced by massive stars and the common view (Henry & Worthey, 1999; Prantzos et al., 2020) is that their abundances should evolve in lockstep inside and among galaxies. Panels d through r of Fig. 13 show the values of Ne/O, S/O and Ar/O as a function of $\omega$ in the left column, and as a function of O/H in the middle and right columns. The meaning of the symbols is exactly the same as described for N/O. For S/O and Ar/O, two different determinations are shown, depending on the ions (and thus, ICF prescription) used as indicated along the axes.

We can readily see that in most panels, the observational points do not gather along horizontal lines. We now discuss these elements separately.

For the BCG sample, the Ne/O ratio shows a slight increase as O/H increases. The slope of the regression line is $0.07\pm\ 0.04$, similar to what was found by Izotov et al. (2006) for their sample. Our GHR sample shows a clear opposite trend, with a slope of $-0.66\pm\ 0.26$. One must note, however, that the uncertainties in Ne/O due to the ICFs are much larger in this metallicity region—due to the fact that the values of $\omega$ are generally small. It is not excluded that for the GHRs of our sample, our ICFs produce an incorrect trend, underestimating the neon abundance especially for objects of low excitation.

For Ar/O we have two possible estimates, one based on Ar⁺⁺ only, the second based on Ar⁺⁺ + Ar⁺³. This second estimate cannot be used for our GHRs sample, where the [Ar iv] $\lambda$4740 line is not measured. We note that the ICF based on Ar⁺⁺ + Ar⁺³ leads to a smaller dispersion in the Ar/O ratios for the BCG sample, despite the fact that it uses an additional line thus involving an additional observational error. This argues in favor the latter ICF being more reliable. The uncertainties in Ar/O for the BCG sample are very small. We note that the slope of the regression line with respect to O/H is consistent with zero with a very small dispersion.

We note also that the average value of Ar/O is slightly larger when using the ICF based on Ar⁺⁺ + Ar⁺³ and closer to the solar Ar/O ratio as given by Asplund et al. (2009) (however given the uncertainty in the solar Ar abundance which, as well as for Ne, is probably larger than quoted, we cannot use the solar value as an argument pro or contra our ICFs). Despite the warning given in Sect. 3.1 it turns out that the ICFs based on Ar⁺⁺ + Ar⁺³ give quite reasonable values.

For S/O we again have two possible estimates, one based on S⁺ alone, the second depending on S⁺+ S⁺⁺. The values obtained in the first case are significantly larger than in the second case (by about 0.6 dex at low metallicities). The uncertainties due to the ICFs are especially large for the first case (typically $\pm$0.2 dex or more at high values of $\omega$), less so for the second case We also note that, for the GHR sample, the first ICF (panel m) gives an increase in S/O as $\omega$ increases while the second ICF (panel p) gives a slight decrease as $\omega$ increases. Obviously at least one of the two estimates of S/O must be wrong. The lines emitted by S⁺ come from the outskirts of the H ii regions, and adding the information from S⁺⁺ obviously allows a better representation of the object. Therefore we do not recommend ICFs based on S⁺ alone. If the only information one has on sulfur in a spectrum comes from [S ii] lines one must be very suspicious on the derived abundances.