© ESO, 2012

1. Introduction

The emission and absorbtion lines of star-forming regions can inform us about the physical conditions of the gaseous media, such as abundances, temperature, and ionization degree. They also provide us with estimates of the ages, masses, and composition of the stellar population and the properties of the ionizing stellar clusters (see e.g. Hägele et al. 2011, hereafter H11; and references therein). The electron temperature is a property crucial for interpreting emission line intensities. When we investigated the [S iii] temperatures in the literature, the study of Hägele et al. (2006, hereafter H06), who compared the temperatures inferred from collisionally excited lines of [S iii] with those of [O iii], caught our attention. Reproducing these with photoionization models was set by some of us (RM and LB) as a potentially instructive exploratory exercise that turned out more challenging as we proceeded with our analysis. Our paper reports on this exploration and is a work in progress, mainly because the uncertainties in the data are substantial, as discussed in Sect. 2. In essence, in nebulae where measurement errors are the smallest, we find that the T_O   iii temperatures exceed the values derived for T_S   iii by up to  ~3000 K, a trend that we were not able to reproduce using simple photoionization calculations. We consider it remarkable (yet with insight predictable) that the relationship between these two temperatures is quite insensitive to strong changes in the input parameters of photoionization models. We subsequently explored models with (1) metallicity inhomogeneities; (2) temperature inhomogeneities; (3) non-Maxwell-Boltzmann energy distributions (hereafter a κ-distribution), and completed our investigation with models that (4) combine photoionization with shock heating. We hope that the current work may motivate observers to undertake observational studies that are able to substantially reduce errors in the measurements of the nebular and auroral [S iii] lines.

Considering the extra difficulties in measuring the [S iii] lines in the infrared and the auroral counterpart at 6312 Å, what alternative species can be used to infer the electron temperature from an auroral-to-nebular line ratio? There are the two well-known [N ii] or [O ii] temperature diagnostics, but they carry their own set of uncertainties, which renders an interpretation of these diagnostics ambiguous and more model dependent. For instance, the electron temperatures deduced from the [N ii] 5755   Å/(6548   Å + 6583   Å) ratio are frequently higher than those of the corresponding [O iii] lines (e.g. Tsamis et al. 2003). According to Stasinska (1980) and Copetti (2006), the increase in temperature outward as a function of nebular radius is likely caused by a combination of hardening of the radiation field with increasing optical depth, and stronger cooling from fine-structure lines of [O iii] in the inner parts of the nebulae where O⁺² dominates (see Stasinska 1980). This expectation is not fully confirmed, however, due to the significant contribution of recombination to the excitation of the auroral [N ii]λ5755 line (Rubin 1986; Tsamis et al. 2003) coupled with the potential presence of high-density inclusions in nebulae (Viegas & Clegg 1994), all of which indicates that [N ii] fails to qualify as a suitable alternative. The [O ii] nebular-to-auroral line ratio is even more susceptible to contribution from recombination and to the potential presence of high-density inclusions in nebulae (Viegas & Clegg 1994), which may cause deceptively high temperatures to be derived from the [O ii] 3727   Å/(7320   Å    +    7330   Å) ratio. Furthermore, this line ratio is more sensitive to variations in the extinction curve (R_V) or to second-order effects from an inhomogeneous dust distribution, which standard reddening corrections cannot correct for. In comparison to N⁺ and O⁺, the specie S⁺² presents the clear advantage of strongly overlapping O⁺² in the higher excitation region of a nebula. An appealing alternative might be the [Ar iii] auroral-to-nebular line ratio 5192   Å/(7136   Å    +    7751   Å) (Keenan et al. 1988; Copetti 2006), but there are very little data available for these lines owing to their intrinsic weakness.

Studying physical properties of emission line objects at different redshifts with accurate and reliable methods is a fundamental issue when one aims to detect possible evolutionary effects such as systematic differences in the chemical composition of these objects. We discuss in Sect. 4 that the observation that the T_O   iii temperatures are higher than T_S   iii may be connected to gaseous nebula phenomena, which are not yet fully understood. One of these phenomena is the abundance discrepancy factor (ADF), which is the ratio of the abundance derived for a heavy-element ion from its optical recombination lines (ORLs) to the abundance derived for the same ion from its ultraviolet, optical, or infrared collisionally excited lines (CELs; Tsamis 2004). In Hii regions, the ADF is about 2 and has been the subject of various interpretations (Peimbert et al. 1995; Esteban et al. 2002; Tsamis & Péquignot 2005; Pérez-Montero & Díaz 2007; Tsamis et al. 2003, 2011; Nicholls et al. 2012; López-Sánchez et al. 2012, and references therein).

For the present study we carried out several models using the multipurpose photoionization-shock code mappings ie that allows us to explore five options: (1) simple photoionization models where we varied the main input parameters; (2) models that consider the effect of temperature inhomogeneities¹; (3) a combination of models of widely different metallicities; (4) models that consider κ-distributions of electron energies (non-Maxwell-Boltzmann energy distributions); and (5) models that combine shock excitation with photoionization. The article is structured as follows: Sect. 2 contains the data analysis, our modeling strategy is described in Sect. 3, and the comparison of models with the data is discussed in Sect. 4. Brief conclusions are presented in Sect. 5.

Fig. 1 Relation between TS   iii and TO   iii temperatures inferred from the [S iii] and [O iii] auroral-to-nebular line ratios for the objects studied by PM06, including the three Hii galaxies later observed by H06. The light gray dashed line depicts the locus of equal temperatures (TS   iii = TO   iii), while the solid line indicates the error-weighted least-squares fit to the data (see Sect. 2).

2. The data base

The analysis carried out in this paper relies on two temperature-sensitive line ratios that involve auroral and nebular line transitions of the same ion. These are [O iii] 4363 Å/(4959 Å + 5007 Å) and [S iii] 6312 Å/(9069 Å + 9532 Å). The inferred temperatures assuming the low-density limit hereafter are labeled T_O   iii and T_S   iii. Our sample corresponds to a compilation of data (99 objects) initially published by Pérez-Montero et al. (2006: hereafter PM06) and subsequently augmented by three Hii galaxies observed by H06. In total, the sample consists of the spectra of 42 Hii galaxies, 34 giant extragalactic Hii regions (GEHRs), 14 Galactic Hii regions, and 12 Hii regions from the Magellanic Clouds. Figure 1 shows the temperatures inferred from a total of 102 pairs of [O iii] and [S iii] line ratio measurements as well as the associated error bars. For clarity’s sake, the error bars are not drawn in the figures. Hereafter, we will use the terms Hiiregion or nebula as generic expressions representative of objects from the sample.

As is evident from the distribution, a correlation between the two temperatures is observed in Fig. 1. An error-weighted least-squares fit of the data results in the following linear regression: (1)with a correlation coefficient of 0.85. The light gray dashed line depicts what would be the locus of equal temperature values in this plot and in subsequent figures. The correlation is markedly steeper than that inferred by Kennicutt et al. (2003) of T_S   iii = 0.83 × T_O   iii − 1700   K from a study of 20 Hii regions in the spiral galaxy M101 (these objects are included in our much larger sample). We note a systematic inequality between the two temperatures in Fig. 1. In particular, lower temperature objects T_S   iii < 14   000 K tend to lie below the equal temperature dashed line. Towards the high temperature end (T_S   iii > 14   000   K), however, the error bars and the overall dispersion of the different measurements become quite large. Due to this and to the absence of any clear trend in the data, our analysis will disregard the high temperature upper right quadrant in Fig. 1.

Given the distribution of the error bars and the increasing dispersion of the points toward higher temperatures, it is doubtful that a single linear regression provides an appropriate description of the current data set. We suggest that a bimodal distribution might be a better description of the data: 1) the lower left quadrant, which shows measurements that are correlated and possess smaller error bars; and 2) the upper right quadrant, where no clear trend is apparent and where the error bars are noticeably larger. Our study and the modeling presented below focus on the temperature behavior within this lower left quadrant and models the temperature departure from the dashed line in Fig. 1, that is, the domain T_O   iii > T_S   iii below T_O   iii < 14   000 K. Even if the error bars had been underestimated by PM06, a proportionally higher noise amplitude would not erase the general trend that the [O iii] temperatures lie above those of [S iii]. Could this trend be the result of a bias affecting the weaker lines? Rola & Pelat (1994) established that the “measured” intensities in poor signal-to-noise (S/N) line measurements are strongly biased toward overestimating their intrinsic values. The errors reported by PM06 for the vast majority of objects are markedly larger for T_S   iii than for T_O   iii . Because auroral lines are much weaker and the error in subtracting the underlying continuum larger, the errors pertaining to the weaker auroral lines are larger. The statistical bias studied by Rola & Pelat (1994) would in our case lead to an overestimate of [S iii]λ6312, hence of T_S   iii , which is the opposite of the observed trend. Telluric absorption of the infrared [S iii] 9069 Å and 9532 Å lines or contamination of the 6312 Å line by atmospheric emission² from nearby [O i] 6300 Å would both contribute to enhancing the inferred T_S   iii values, in contrast to the observed trend. In our view, the spread in temperatures in the lower left quadrant is genuine, but we cannot rule out that unforeseen systematic errors are contributing to it.

Improving the accuracy of temperature measurements will be required before definitive conclusions can be drawn. The infrared [S iii] lines are intrinsically more difficult to measure. The risk of systematic errors caused by atmospheric telluric emission, which increases toward longer wavelengths near 9500 Å and beyond, remains a concern. As databases will grow, a significant improvement in the quality of the data should be the first priority. For instance, it might become preferable to discard data for which the 9532 Å/9069 Å ratio does not agree with the theoretical value of 2.46 (from Podobedova et al. 2009) within a reasonable S/N margin.

3. Model calculations

An abundance analysis of the complete sample has already been carried out by H06 and PM06, in which the authors took into account the measured nebular temperatures as well as the intensities of the forbidden lines with respect to the recombination lines. Such an abundance analysis will not be repeated in this paper. Rather, our aim will be to investigate whether one can model [O iii] temperatures that are higher than those derived from [S iii], as is observed in the lower left quadrant of Fig. 1. It follows that our study will focus on nebulae that tend to lie on the high-metallicity end of the data points. As an exercise, we can estimate the range of abundances that the objects in this quadrant cover. To do so, we can use the abundance calibration presented in Appendix A. It was derived from the data used by H06 and consists of a least-squares fit of the observed [O iii] temperatures as a function of the deduced O/H abundance ratio. When adopting the value of T_S   iii = T_O   iii = 14   000 K as a marker that defines the lower left quadrant, this corresponds to Z ≈ 0.15   Z_⊙ (see Appendix A). If we select instead the temperature at which the linear regression crosses the dashed line, which is at T_O   iii = 12   400 K, the inferred metallicity is Z ≈ 0.25   Z_⊙. In summary, the nebulae whose temperatures we attempt to model have metallicities Z ≳ 0.2 Z_⊙.

In standard (steady-state) photoionization models, we disregard hydrodynamical effects and assume that equilibrium ionization and equilibrium temperature are reached at every point in the nebula. In that case, the temperature is determined locally by the condition that the energy is gained by photoionization and is balanced by the energy loss through line and continuum emission. As one increases the gas metallicity, the equilibrium temperature decreases. A systematic exploration of the [O iii] and [S iii] temperature behavior with respect to one another was initially carried out as we varied in turn the metallicity Z, the ionization parameter (U), and the spectral energy distribution (SED). It turned out to be a daunting challenge to have T_O   iii significantly exceed T_S   iii . With standard photoionization, we only obtained promising results when we considered combining ionized regions of greatly varying metallicities (Sect. 3.3). As an alternative to metallicity inhomogeneities, we subsequently initiated the exploration of three “extraneous” factors: (1) the possible contribution of shock excitation; (2) the presence of temperature inhomogeneities as introduced by Peimbert (1967; cf. Peña-Guerrero et al. 2012a, and references therein); and (3) the effects of a departure from a Maxwell-Boltzman distribution for the electron energies (Nicholls et al. 2012: hereafter NDS12).

3.1. Photoionization calculations with mappings i

To compute nebular spectra, we used the multipurpose shock/photoionization code mappings ic described by Ferruit et al. (1997; see also Binette et al. 1985), in which the ionization and thermal structure of the nebula are computed in a standard fashion assuming local ionization and thermal equilibria throughout the nebula. The calculation proceeds at a progression of depths in the nebula until the entire ionization radiation has been absorbed. The updated version of the current code (mappings ie) now includes the option of considering a non-Maxwell-Boltzmann distribution of the electron energies when calculating excitation rates, as described in Sect. 3.4.

3.1.1. Geometry

We adopted a simple slab geometry for the ionized gas, since for the lines of interest the integrated line spectrum does not appreciably differ from that calculated with spherical geometry (Shields 1974). This is the geometry preferred for the Orion nebula that Baldwin et al. (1991) modeled in detail. Hii regions are very chaotic in appearance, showing evidence of strong density inhomogeneities, with much of the emission coming from condensations or ionized “skins” of condensations. Modeling this level of complexity is beyond the current means of photoionization models and is undesirable, unless the observations were sufficiently detailed to constrain such models. Finally, the definition of the ionization parameter in the spherical case depends crucially on the Strömgren radius size (R_s), which depends on the volume filling factor of the gas and the radius R_g of the empty gas volume near the ionizing stars, both of which are poorly constrained if not downright ill-defined.

We adopted the usual definition of the ionization parameter: (2)were c is the speed of light, h the Planck constant, and R_g the distance of the slab from the ionizing star, in other words, U is the ratio between the density of ionizing photons impinging on the slab (ϕ_H/c) and the slab H density n_H . With a slab geometry, we assume that the geometrical thickness R_s − R_g of the ionization-bounded slab is insignificant with respect to R_g.

In our grid of photoionization calculations, we varied U over a wide range from 0.001 to 0.046, although in practice a more limited range may suffice (e.g., Evans & Dopita 1985). Campbell (1988) found that most of the nebulae in her sample lie in the range  − 1.9 ≳ lgU ≳  − 2.5. All photoionization calculations were ionization-bounded and we assumed the low-density regime, by adopting n_H = 10    cm^-3. Some emission lines of low critical density, such as [S ii]λλ6716, 31 or [O ii]λλ3727, 29, would be affected by collisional deexcitation in those few cases where the nebular densities reach values as high as  ~500  cm^-3. Nevertheless, since the overall energy budget remains set by the majority of other lines with much higher critical density, the conclusions we reached about the T_O   iii and T_S   iii temperatures do not depend on the density for the overwhelming majority of Hii regions.

3.1.2. Reference set of abundances

We express the gas metallicity Z with respect to the solar abundance set of Asplund et al. (2006), where He/H = 0.085 and the other elements are C,N,O,Ne,Mg,Si,S,Ar,Ca,Fe:   H = Z × (250,60,460,69,34,32,14,1.5,2,28) × 10^-6   Z_⊙, where Z is the metallicity in solar units. We computed nebular models of varying metallicities covering a range from 1% solar (Z = 0.01   Z_⊙) to more than twice solar (Z = 2.5   Z_⊙).

3.1.3. Stellar spectral energy distributions

We explored different stellar SEDs for the ionization source of the Hii regions. This did not turn out to be a sensitive parameter. In the models reported below, we compare photoionization models with an SED from a zero-age instantaneous stellar burst assuming a Chabrier (2003) initial mass function with M_sup  = 100 M_⊙, and solar metallicity, obtained from the stellar population synthesis models by Bruzual & Charlot (2003). The stellar atmospheres used to derive the stellar ionizing SED are taken from the Kurucz spectra library (distributed by Lejeune et al. 1997) for stars with T_eff   < 50 000   K, and the Rauch (1997) NLTE model atmospheres for the higher temperature stars (Magris et al. 2003, hereafter M03). We also explored the alternative case of a single stellar atmosphere of temperature T_eff, using the LTE atmosphere models of Hummer & Mihalas (1970, hereafter HM70).

3.2. Models that consider temperature inhomogeneities

It has been proposed that the differences between the temperatures from ORLs and CELs could be due to temperature inhomogeneities (e.g. Peimbert 1967). To evaluate how these inhomogeneities affected the line ratios, we adopted the temperature inhomogeneity characterization described in Binette & Luridiana (2000, hereafter BL00), who used different temperatures for the CELs in accordance with the energy separation of the atomic levels involved for each line emission calculation. BL00 did not assess the physical driver of these inhomogeneities. These authors simply assumed their existence and adopted a simple scheme in which the inhomogeneities corresponded to hot-spots as a result of an unknown heating process, distinct from the photoelectric heating. These authors considered the global effects of temperature inhomogeneities on the emission lines, based on the statistical approximation introduced by Peimbert (1967), which quantifies the brightness increase or decrease of each line in terms of the departure of the temperature from the average temperature T_o. For each transition, Peimbert used an expansion in Taylor series of the temperature that takes the functional dependence with temperature of each line emissivity into account. Following Peimbert (1967), the mean nebular temperature T_o , is defined as (3)for a homogeneous metallicity nebula characterized by small temperature inhomogeneities, n_e is the electron density, T the electron temperature, and V the volume over which the integration is carried out. The rms amplitude t of the temperature inhomogeneities is defined as (4)We simplified the expression presented by Peimbert (1967), in which t² depends on the density of the ionic species considered (n_i) while in the above equations we implicitly consider only ionized H (by setting n_i = n_H ⁺  = n_e).

In photoionization calculations, one determines the local equilibrium temperature at every point in the nebula, T_eq , that satisfies the condition that the cooling by radiative processes equals the heating due to the photoelectric effect. Assuming the existence of temperature inhomogeneities and adopting the scheme of BL00 that describes these as hot-spots, we can define T_eq as the temperature floor above which putative temperature inhomogeneities take place. These inhomogeneities are taken into account only in the statistical sense, that is, through the use of appropriately “biased” temperatures that depend on the mean squared amplitude of the inhomogeneities and on the functional dependence of the (de)excitation rate on temperature. In this context, t² also represents the amount of external energy that the hot spots contribute to the nebular energy budget (BL00). Computing the effects of the hot-spots on the line intensities is a two-step process: first T_o is derived, and second the biased temperatures corresponding to each atomic transition are calculated. For the hot-spot scenario that BL00 simulated, the following expression for T_o provided a satisfying parametric fit (5)A value of γ of  ≃  − 15 was found to provide an adequate parametrization of the behavior of T_o .

To compute individual recombination line intensities in the presence of small inhomogeneities (which are proportional to T^α), the code uses a biased mean temperature T_rec given by (6)Because temperature inhomogeneities have a lesser impact on recombination lines than on collisionally excited lines, we adopted the simplification of considering a single constant value of the power α =  −0.83 for all recombination processes.

For collisionally excited lines, when evaluating the excitation rates ( ∝ T^β_ij   exp [ − ΔE_ij/kT] ) and deexcitation ( ∝ T^β_ji) of a given multi-level ion, each rate ij (population) or ji (depopulation) is calculated using T_o (instead of T_eq ), and then multiplied by the appropriate correction factor, either (7)in the case of excitation, or (8)in the case of deexcitation. k_B is the standard Boltzmann constant and β_ji = β_ij. These factors, adapted from the work of Peimbert et al. (1995), were applied to all collisionally excited transitions calculated by mappings ie. It typically results in a significant enhancement of collisional rates by an amount that strongly depends on the ratio ΔE_ij/kT_o for the selected transition ij.

3.3. Models that consider abundance inhomogeneities

We explored higher complexity photoionization calculations in the form of small-scale metallicity inhomogeneities within the ionized gas. One could imagine, for instance, enriched condensations embedded in a more diffuse gas of lower metallicity. A justification for such a scenario was provided by Tenorio-Tagle (1996). It is based on the fact that the metal-rich ejecta from type II supernovae (SNe) ought to follow a long excursion into the galactic environment before they are able to mix with the ISM. In this scenario, the violently ejected matter generates large-scale superbubbles, which eventually cool non-uniformly, creating in the process multiple re-pressurizing shocks. This leads to the formation of small, dense cloudlets of metal-rich gas, which eventually fall back toward the disk of the galaxy because of gravity. An efficient metallicity mixing with the surrounding unenriched ISM will only take place when these “metallic droplets” (or cloudlets) are caught within the nebulae associated to star-forming regions. Stasinska et al. (2007) quantitatively investigated the hypothesis that the photoionization of these droplets can explain the ADFs observed in Hii regions. The authors derived bounds of 10¹³–10¹⁵ cm on the droplet sizes to (1) explain that they have not yet been spatially resolved and (2) to ensure that they survive long enough their destruction by diffusion.

The possibility of metallicity inhomogeneities in Hii regions has been studied on two scales: (1) on a scale in which they have not yet been spatially resolved³ (Tsamis et al. 2003); and (2) on a scale associated with proplyds observed in Hii regions (Tsamis et al. 2011). An example of case (1) is given by the work of Tsamis & Péquignot (2005), who modeled the line spectrum of the central region of 30 Doradus in the LMC in detail. Their dual metallicity model consists of gas of approximately normal LMC composition (0.36 Z_⊙) co-extensive with a hydrogen-poor component, in which most elements are enhanced by  ≈ 0.9 dex. The relative weight of the hydrogen-poor component is relatively low, at an 8% level of the integrated Hβ flux. In addition to successfully fitting the ORLs from CNO elements (in addition to the UV/optical CELs), the model of Tsamis & Péquignot leads to greatly improved predictions for several other optical and IR lines as well. For planetary nebulae (PNe), the idea has been studied earlier (Liu et al. 2000; Ercolano et al. 2003; Tsamis et al. 2004; Bohigas 2009; Yuan et al. 2011), although the enriched gas component has a totally different origin than that envisaged by Tenorio-Tagle (1996).

To study abundance inhomogeneities, we adopted a dual metallicity modeling approach in which we combined pairs of isobaric photoionization models that simply differ in metallicity. Although we did not consider all physical implications of embedded condensations, such as an intermingled diffuse field, thermal conduction, or a change in ionization parameter, our calculations suffice in evaluating how abundance inhomogeneities can affect the [S iii] and [O iii] temperatures. The technique we employed is the following: for a given pair of photoionization models, A and B, we multiplied all line fluxes of model A by a weight factor w_i and added them to the corresponding fluxes from model B after multiplying these by the complementary weight 1 − w_i. Subsequently, we calculated the line ratios of the resulting “averaged” model and derived physical quantities of interest, such as the “apparent” T_S   iii and T_O   iii temperatures. We built sequences of these models in which only w_i is varying along the sequence. We explored the case where models A and B differ only in their metallicity and the case where they differ in both U and Z. Since the results turned out to be qualitatively similar in our temperature diagram (provided the contrast in U was  ≲ 10), we will only show cases where both models A and B share the same value in U . Our sequences can be defined by three parameters: the ionization parameter U_A , the metallicity Z_A of the (lower abundance) model A, and the metallicity contrast between models A and B given by θ = lg(Z_B/Z_A).

3.4. Models with a non-Maxwell-Boltzmann energy distribution

Recently, NDS12 explored the possibility that electrons in Hii regions and PNe depart from a Maxwell-Boltzmann equilibrium energy distribution. To represent the electron energy distributions, these authors adopted the parametrization given by a κ-distribution, which is a generalized Lorentzian distribution initially introduced by Vasyliunas (1968). It is widely used in the studies of solar system plasmas (Livadiotis & McComas 2011; Livadiotis et al. 2011) where evidence of suprathermal electron energies abounds. NDS12 found that a κ-distribution of electron energies resolved many of the difficultie encountered when attempting to reproduce the temperatures observed in nebulae. They concluded that mild departures from a Maxwell-Boltzmann distribution is sufficient for reproducing the spread in temperature values found in diagnostic diagrams such as T_rec vs. ⟨T_OII + T_NII⟩ or T_SII vs. T_OII in Hii regions or T_rec vs. T_O   iii in PNe. The κ-distribution is a function of temperature and κ. In the limit as κ → ∞, the distribution becomes a Maxwell-Boltzmann distribution. From examining data from Hii regions and PNe, NDS12 found that κ ≳ 10 is sufficient to account for most objects they compared their calculations with. These authors proposed various mechanisms that would result in a suprathermal distribution, such as (1) magnetic reconnection followed by the migration of high-energy electrons along field lines, and by the development of inertial Alfvén waves; (2) the ejection of high-energy electrons from the photoionization process itself (when the ultraviolet source is very hard, as in PNe or AGN); (3) or from photoionization of dust (Dopita & Sutherland 2000); and (4) by X-ray ionization, resulting in highly energetic (~1 keV) inner-shell (Auger process) electrons (e.g., Aldrovandi & Gruenwald 1985).

To implement this new option in mappings ie, we proceeded as follows. In the case of collisionally excited lines of a given multi-level ion or atom (including Hi and He ii), each excitation rate ij is calculated first, using the collision strength Ω_ij, the excitation energy E_ij, and the local equilibrium temperature T = T_eq, as is customary (Osterbrock 1989). Next, each rate is multiplied by the corresponding correction factor (9)which is an expression taken from the work of NDS12 that corresponds to the ratio of the rate computed assuming a κ-distribution over the rate that would be derived assuming a Maxwell-Boltzmann distribution instead. To prevent numerical overflows, the natural logarithm of Eq. (9) is evaluated first. To compute the function Γ(z), we adopted the expression (Nemes 2010) (10)To derive the correction factors for the deexcitation rates ji, we used the above expression, but evaluated using ΔE_ji = 0. This coefficient’s value () is unique for a given combination of κ and T. It can be applied to threshold-free processes when their cross section approximately behaves as  ∝ E^-1, as is the case with recombination rates. However, rather than applying correction factors to the recombination rates, we opted to use an effective recombination temperature T_rec that would result in the same rate coefficients as those obtained by applying the correction factor cf₀. Since recombination rates approximately vary as  ∝ T ^− 1/2, it follows that the effective temperature to be used when computing recombination rates is given by T_rec ≃ T/(cf₀)².

We also implemented the enhancement of collisional ionization that can result from κ-distributions. There are two sets of correction factors: one applies to the excitation-autoionization process (E-A), and the other to direct collisional ionization process (D-I). In the case of E-A, we found from the work of Owocki & Scudder (1983) that Eq. (9) can be used to derive the correction factors when the excitation energy is replaced by the relevant autoionizing energies. These coefficients are then applied to the E-A rates. In the case of D-I, mappings ie uses the method proposed by Arnaud & Rothenflug (1985) to derive the ionizing rates, which is based on a four-parameter functional fit of the D-I cross-sections σ(E_ij) (their Eq. (1)), for which we lack an analytical prescription that would allow us to compute the correction factors when considering a κ-distribution. We therefore opted for a numerical integration. For each electron shell j of ion i we numerically integrated σ(E_ij) assuming the selected κ-distribution for the electron energies. We divided the result by the same numerical integration, assuming instead a Maxwell-Boltzmann distribution. The resulting ratio defines the correction factor , which was then applied to the corresponding D-I rate computed with the Arnaud & Rothenflug (1985) expression. For the calculations presented here, the enhancement of collisional ionization caused by a κ-distribution was found to have no impact on the ionization structure of Hii regions.

3.5. Photoionization models with shock heating

To explore the possible impact of shocks in Hii regions, we used mappings ie to calculate steady-state plane-parallel shock models, assuming a zero preshock magnetic field. This shock wave, by construction, is set to propagate in a nebula that is already photoionized by hot stars. The nebular gas constitutes the preshock zone with a temperature T_pre that is evaluated assuming photoionization equilibrium. We did not include the surrounding unshocked nebular emission in our calculations. Rather than assuming an arbitrary opacity ahead of the shock front, we adopted an unabsorbed ionizing SED. This means that we assumed that the shock wave occurs close to the inner nebula boundary, if homogeneous, or close to the inner boundary of the brightest filaments otherwise. Throughout the calculations, we assumed the low-density limit regime by adopting a density n_pre = 1.0  cm^-3, which we term the preshock density. The basic input parameter of these models will be the postshock temperature T_post . Using the values of n_pre , T_pre and T_post , we computed the Rankine-Hugoniot conditions to determine the shock velocity V_s and the postshock density n_post . Both quantities appear in Table 1. This density n_post is about four times that of n_pre when V_s corresponds to a high Mach number. Downstream of the shock front, as the postshock gas progressively cools, its density increases since the gas pressure is conserved throughout the postshock region. We adopted the previous zero-age instantaneous starburst model of M03 (Sect. 3.1.3) as stellar SED that is responsible for ionizing and heating the preshock zone. The ionization parameter is defined by Eq. (2) as before, except that n_pre substitutes n_H .

One caveat of our models is, however, that the putative shock wave is set to travel “inward”, toward the front face of the nebula, while one would favor the opposite, that is, we expect the shock wave to enter the nebula from the front face (that is, from the side exposed to the incoming stellar UV and stellar wind) and to propagate outward. The reason we are limited to the “inward” direction for the shock propagation has to do with the way the opacity is integrated, which is subsequently used to determine the radiative transfer of the ionizing flux. In mappings ie, this integration starts at the shock front and proceeds downstream in parallel with the temporal evolution of the shocked gas. The opacity is integrated along the flow. This way, the radiative transfer of the radiation emitted within the hot layers of the shock as well as that from the hot stars can be resolved. However, the opposite “outward” shock propagation might be preferable, for instance in scenarios where shocks are generated by a stellar wind or other stellar ejection mechanisms. In these scenarios where the ionizing photons enter downstream, from the backside of the shock, one needs to know from the start how the opacity behaves across the entire shocked layer. This requires an iterative method to ensure a proper radiative transfer determination. Aldrovandi-Viegas & Contini (1985) developed such an iterative procedure for their code SUMA (Viegas & Contini 1997) and applied it to model the narrow line region (NLR) of active galactic nuclei (AGN). This option, however, is not available with the code mappings ie. To circumvent this limitation, we computed models for which the external ionizing flux is sufficiently intense so that the shocked layer’s depth remains a small fraction of the total ionized depth of the nebula in absence of shocks. When this condition is satisfied, the shocked layer can be expected to show a similar structure, whether it is propagating inward or outward, with respect to the direction of the dominant radiation field.

To determine the preionization conditions of the gas before it enters the shock, we assumed equilibrium ionization and equilibrium temperature. This is a straightforward assumption, since the preshock gas is immersed in a radiation field so intense that it is already highly ionized and at or near equilibrium temperature. The ion and electron temperatures were set to be equal in all shock models. Starting at the shock front, the adiabatic cooling (and recombination) of the plasma is followed only until such time as cooling by radiative processes has sufficiently decreased, so that it becomes equal again to the heating due to photoionization. Hence, the line intensities that we computed strictly apply to the shock-heated layers (that is, the out-of-equilibrium layers beyond which steady-state photoionization eventually settles again, but this zone is not computed). We recall that the models do not include line emission from the photoionized preshock region either.

Fig. 2 Photoionization calculations with an SED corresponding to a zero-age instantaneous starburst model from M03 with Msup = 100    M⊙. A sequence of models of increasing metallicity Z is represented by the magenta solid line, for which U is constant at 0.01. The metallicities covered a range from Z = 0.01   Z⊙ to Z = 2.5   Z⊙, where Z increases by 0.2 dex from model to model. Three sequences of models are superimposed in which U is increased (rather than Z), covering the range 0.001 to 0.046. Each sequence has a different metallicity, as follows: Z = 0.01   Z⊙ (cyan dotted line), Z = 0.25   Z⊙ (light green dotted line) and Z = 1.0   Z⊙ (blue dotted line). In this and subsequent figures, a light gray dashed line represents the locus of equal temperatures.

To facilitate comparisons between photoionized plasma and steady-state shock waves, we will make use of column densities. For the assumed plane-parallel geometry, the line emission measures are what models integrate to derive line luminosities. However, our interest in referring to integrated column densities is that the internal structure (i.e., the behavior with depth) of the ionization, the temperature, or the cumulative line ratios are the same when comparing homologous models and when the quantities of interest in these models are plotted as a function of column depth. In the familiar photoionization case, models with different densities are homologous, provided the ionization parameter U is the same, at least in the regime of sufficiently low densities where collisional deexcitation does not affect the thermal structure. In the case of shocks, constant V_s models (in the low-density regime) are homologous, and their structure remains unchanged as a function of column density, independently of the postshock density n_post . The reason is that although the internal energy of the shocked gas after the shock discontinuity scales as n_post , the local cooling rate scales as the density squared, which means that the cooling will proceed on a timescale that evolves as n_H^-1, that is, decreasing with increasing postshock density. This keeps the behavior of the temperature or ionization invariant as a function of depth, when expressed as a column density. In the case of combined shock-photoionization models, models will be homologous if both U and V_s are kept constant when changing the input density. In the case of pure photoionization, the integrated column density of ionized H for an ionization-bounded model is given⁴ by , where α_B is the case-B recombination rate coefficient of hydrogen. For the ionization parameter value of U  = 0.01 that we assumed in Sect. 4.5 for the combined shock-photoionization models, is 8.2    ×    10²⁰  cm^-2 in the pure photoionization case with Z = 1.0   Z_⊙ (in absence of shocks). In the case of a steady-state shock wave, increases monotonically with shock velocity⁵V_s without external photoionization. It is a manifestation of the longer elapsed time required for the shocked plasma to recombine and cool from T_post to  ~100 K. In our calculations, for which the shock wave was immersed in a photoionized plasma, when integrated over the out-of-equilibrium layers only, is as low as 2.5    ×    10¹⁹  cm^-2 for a shock velocity of 56 km   s^-1. Hence, the aforementioned condition of considering only shocked layers that are geometrically thin with respect to an ionization-bounded photoionized layer is clearly satisfied with the velocities considered in this work ( ≤ 65   km   s^-1).

4. Model results and discussion

In figures without the error bars, such as Fig. 2, the span in [O iii] temperatures to the right of the dashed line in the lower left quadrant stands out more clearly. At the low-temperature end, T_O   iii can exceed T_S   iii by as much as  ~3000 K. We now investigate which types of models or parameter values would cause the [O iii] temperatures to be appreciably higher than those from [S iii]. We begin with homogeneous photoionization calculations.

Fig. 3 Comparison of two sequences in metallicity that differ by the shape of their ionizing SED: a single temperature star of either Teff = 55   000 K (light green line) or Teff = 40   000 K (blue dotted line). The previous zero-age instantaneous starburst SED of M03 with Msup = 100    M⊙ appears underneath (magenta line). Each sequence comprises 13 models with the metallicity increasing by a factor 0.2 dex from model to model, starting at Z = 0.01   Z⊙ and ending at 2.5 Z⊙. The ionization parameter is U = 0.01 in all calculations.

4.1. The pure homogeneous photoionization case

When the plasma temperature is set by the condition that equilibrium ionization and equilibrium temperature apply locally everywhere within the nebula, and that there are no inhomogeneities in either the local densities or the metallicities, one finds that varying the input parameters gives little leverage on the behavior of both temperatures. This is illustrated by the models presented in Fig. 2. The magenta solid line shows a sequence of 13 models (all with U = 0.01 and the zero-age instantaneous starburst SED of M03), in which only the metallicity increases along the sequence, by 0.2 dex from model to model, covering the range Z = 0.01 to 2.5 Z_⊙. Higher metallicity models result in cooler nebulae, while the lowest Z models occupy the upper right corner. What we infer from this plot is that below 14 000 K, there is no significant departure of the models from the equal temperature (dashed) line. Only beyond 16 000 K does a modest departure from the equal temperature line begin to occur. To complete the picture, we show three sequences, of six models each, starting at U = 0.001 and ending at 0.046, along which the ionization parameter increases by 0.33 dex from model to model. These three sequences are characterized by different metallicities Z of 0.01 Z_⊙ (cyan dotted line), 0.25 Z_⊙ (light green dotted line) and 1.0 Z_⊙ (blue dotted line), respectively. They illustrate quantitatively that U is not a significant parameter, except at very low metallicities. In this case, the value of U matters because of the cooling provided by collisional excitation of Hi, which plays an increasingly important role in establishing the equilibrium temperature when metal cooling becomes progressively inoperative beyond Z ≲ 0.05 Z_⊙. Increasing U reduces the neutral fraction of hydrogen H⁰/H, which diminishes the cooling that Hi collisional excitation can provide, resulting in a hotter nebula.

In Fig. 3, we compare the case of using an SED consisting of a single star atmosphere model (HM70) with T_eff = 55   000 K (light green solid line) with the case of using the previous M03 SED (magenta line). Both sequences are metallicity sequences covering the interval Z = 0.01 to 2.5 Z_⊙. We also explored stellar SEDs of lower T_eff or SEDs from M03 of lower M_sup and/or of an evolved population. To prevent cluttering in the figure with too many overlapping lines, these models are not plotted, except for the T_eff = 40   000 K HM70 stellar atmosphere SED (blue dotted line). Those models showed that reducing either M_sup or T_eff does not cause any improvement in reproducing the spread in the observed T_S   iii and T_O   iii temperatures of the lower left quadrant. The SEDs that we adopted for the calculations shown in Fig. 3 illustrate that a change in the continuum hardness is not effective in resolving the T_S   iii –T_O   iii temperature issue.

Fig. 4 Sequences of models in which t2 is increased from 0.004 to 0.16 by successive factors of 1.58 ( = 0.2 dex). All calculations were carried out with U   = 0.01 and the zero-age instantaneous starburst model SED with Msup = 100      M⊙ from M03. Each sequence corresponds to a different metallicity Z, which is indicated in the line color legend in Z⊙ units.

4.2. Results of models with temperature inhomogeneities

Because mappings ie offers the possibility of simulating the effects of temperature inhomogeneities following the method described in Sect. 3.2, we explored various runs of models in which the mean squared amplitude of these inhomogeneities (t²) is progressively increased⁶. The t² values that we explored range from 0.004 to 0.16, fully encompassing that found in Hii regions. For instance, in the sample of 28 nebulae used by Peña-Guerrero et al. (2012b, hereafter PG12) to recalibrate the Pagel method, 70% of the objects have t² values within the interval 0.02 ≤ t²   ≤ 0.04, one object has a lower value, while the rest (25%) lie within the range 0.06–0.12. The results of our calculations are shown in Fig. 4. All calculation were carried out with U = 0.01 and the starburst SED from M03. Each t² -sequence comprises nine models, with t² increasing by a factor 0.2 dex (=1.58) from model to model. The different t²-sequences differ only in their metallicity, whose value appears in Z_⊙ units in the legend of Fig. 4. We found that the modifications to the procedure for calculating emissivities (with respect to pure photoionization) suffice to explain the growing divergence of low-metallicity sequences from the equal temperature line as t² is increased. Changes in ionization structure due to calculating the recombination rates at T_rec (see Eq. (6)), rather than at the equilibrium temperature, are minor.

Figure 4 shows that, for models within the lower left quadrant, increasing t² increases the temperatures T_S   iii and T_O   iii by comparable amounts. As a result, the t² -sequences do not cover the spread in the data observed to the right of the equal temperature line in the lower left quadrant. Only within the upper right quadrant do the calculated temperatures diverge from the dashed line as t² is increased. Hence, temperature inhomogeneities cannot account by themselves for the excess in T_O   iii temperatures observed in the lower left quadrant.

Fig. 5 Sequences of models in which the relative weight wi of two photoionization calculations is varied, as described in Sect. 3.3. For each sequence, the metallicity of its model A in solar units appears in the line color legend followed by the contrast θ in dex. An open circle indicates the position where the relative weight applied to model A and B is equal (wi = 0.5). In all the above models, the ionization parameter is U = 0.01 and the SED corresponds to the zero-age instantaneous starburst model with Msup = 100      M⊙.

4.3. Results of models with metallicity inhomogeneities

We now turn to the dual metallicity modeling approach described in Sect. 3.3, which has the advantage of not requiring an external contribution to the energy balance of the nebula. We report on sequences of models that have identical values in ionization parameter U_A = U_B = 0.01. Reducing this parameter to 0.001 did not affect the results in a noticeable way, except in cases where Z_A was  ≪ 0.2   Z_⊙. We recall that along any sequence, only the relative weight w_i applied to models A and B varies, within the interval 0 ≤ w_i ≤ 1. In Fig. 5 we present examples of dual metallicity sequences. The θ = 0.7 dex sequence (cyan line) corresponds to a weighted average of a model A with Z_A = 0.2   Z_⊙ with a model B with Z_B = 1.0   Z_⊙. The blue line sequence and the light green line sequence have the same Z_A, but a higher abundance contrast θ of 0.8 and 0.9 dex, respectively. The maroon dashed line sequence was calculated with a higher value Z_A = 0.32   Z_⊙ and a metallicity contrast of 0.7 dex. We infer from Fig. 5 that sequences with a sufficiently high contrast θ can cover the observed domain of points to the right of the equal temperature (dashed) line. The locus of the various models also suggests that a metallicity contrast extending up to  ≲ 1 dex is required to account for the nebulae with the largest temperature differences of T_O   iii − T_S   iii. A possible drawback might be that the modeling relies on some fine-tuning of w_i to reach the observed T_O   iii value. To give an idea of how T_S   iii and T_O   iii varies with w_i, we use an open circle to indicate the position along each sequence where both models A and B have equal weight values of 0.5. This corresponds to the case where models A and B reprocess the same amount of ionizing photons (and emit about the same flux in recombination lines).

One interesting feature of models that consider abundance inhomogeneities is that the divergence from the equal temperature line becomes markedly smaller toward the high-temperature end. The reason is that beyond Z ≲ 0.1 Z_⊙, metals play a progressively vanishing role in affecting the plasma temperature. This behavior is illustrated by the seven sequences that are presented in Fig. 6, which extend over a much wider metallicity range compared to the previous figure. The abundance contrast θ is fixed at 1.0 dex and the metallicity of the model A of each sequence appears in Z_⊙ units in the figure legend. Data of higher S/N would be essential to ascertain how nebulae behave in the low-metallicity (high-temperature) quadrant of the diagram.

Fig. 6 Seven sequences of dual metallicity models covering a wider range in metallicities. The metallicity contrast between ZB and ZA is fixed at θ = 1.0 dex. Each sequence corresponds to a different metallicity ZA value, which is indicated in the line color legend in Z⊙ units. An open circle along each sequence indicates the position where the relative weight applied to each model is equal (wi = 0.5). In all the above sequences, the ionization parameter is U = 0.01, and the SED corresponds to the zero-age instantaneous starburst model with Msup = 100      M⊙.

Fig. 7 Sequences of models in which the parameter κ has the values of 80, 40, 20, 10, and 5. The larger κ, the closer the model to standard models with a Maxwell-Boltzmann electron energy distribution. These κ sequences differ in their metallicity Z (expressed in Z⊙ units in the line color legend). The dashed magenta line represents a metallicity sequence of standard photoionization models (κ = ∞) (same sequence as in Fig. 2). The light gray arrow illustrates the increase in the TS   iii and TO   iii temperatures when enhanced Hi excitation is left out from the calculations of the model with Z = 0.01 Z⊙ and κ = 5. In all calculations, we adopted U = 0.01 and an SED consisting of a zero-age instantaneous starburst model from M03 with Msup = 100    M⊙.

4.4. Results of models with a κ-distribution

After implementing the effects of a κ-distribution on the line intensities and recombination processes as well as on the equilibrium temperature of a photoionized plasma, we proceeded to compute sequences of models where κ has the values of 80, 40, 20, 10, and 5 (which encompasses the range 6 − 20 that NDS12 explored). In Fig. 7, we present seven such sequences that differ only in their gas metallicity. The gas metallicity of each sequence appears in the figure legend, expressed in Z_⊙ units. The dashed magenta line represents a metallicity sequence with a pure Maxwell-Boltzmann distribution (the same Z-sequence that was plotted in Fig. 2). All calculations assumed U = 0.01. Our results indicate that κ-distributions are able to reproduce the spread in the observed values of T_O   iii vs. T_S   iii inside the lower left quadrant and therefore constitute an interesting alternative to the other options studied in this paper.

With κ-distributions, the higher the excitation energy of a level, the larger the increase in the corresponding emission line intensity. This explains why [O iii]λ4363 grows considerably more than [S iii]λ6312. At a fixed κ value, the ratio between the respective T_S   iii and T_O   iii temperature increases is constant. Hence, independent of metallicities, models with the same κ lie at about the same distance from the equal temperature dashed line in Fig. 7. This is also the case for the other temperature diagrams presented by NDS12. Because mappings ie includes the enhancement in Hi excitation cooling due to the high-energy tail of κ-distributions, T_O   iii and T_S   iii are seen to increase less and less in Fig. 7 when decreasing metallicity at a fixed κ value. A measure of this effect is represented by the arrow in Fig. 7, which represents the increase in the T_S   iii and T_O   iii temperatures when enhanced Hi excitation is left out from the calculations of the model with Z = 0.01 Z_⊙ and κ = 5. Since collisional excitation of Hi is a strong function of the neutral fraction of hydrogen H⁰/H (Luridiana et al. 2003), the importance of the effect described above at low metallicities is governed by the ionization parameter⁷. If we were to increase U , it would result in a lower H⁰/H fraction and this would reduce the importance of Hi excitation cooling present at very low metallicities and therefore result in a higher nebular temperature.

Fig. 8 Six sequences of combined shock-photoionization models along which the shock velocity Vs is increased in locked steps. The preshock conditions are set by photoionization equilibrium, using the parameter U = 0.01 and the M03 SED, and assuming thermal and ionization equilibria (see Sect. 3.5). Each sequence is characterized by a different metallicity Z value, which is expressed in Z⊙ units in the line color legend. For each model, the preshock temperature, the postshock density, the postshock temperature, and the resulting shock velocity are given in Table 1. The first two models of the Z = 0.5 Z⊙ and 1.0 Z⊙ sequences overlap too closely to be distinguished.

4.5. Results from photoionization models with shock heating

A set of calculations was carried out in which we computed the line intensities of a plane-parallel shock wave that propagates within a photoionized gas layer. To simplify the exercise, we adopted the zero-age instantaneous starburst SED and the ionization parameter value of U = 0.01 for the preshock photoionized gas in all calculations. More information about the models is presented in Sect. 3.5.

In Fig. 8 we present six sequences of combined shock-photoionization calculations, which differ by their metallicity Z (see the adopted values expressed in Z_⊙ units in the figure legend). Preshock temperatures (T_pre) are given for each model in Table 1. Along each sequence, it is the post-shock temperature (hence, the shock velocity) that is progressively increased from model to model. The exact values of the postshock temperature, T_post , and the associated shock velocities, V_s , are listed for each model in Table 1 (along with n_post ). The associated shock velocities typically cover the interval 20 to 60 km   s^-1. This velocity range is very similar to what is obtained in detailed kinematical studies of GEHRs and circumnuclear star-forming regions (see e.g. Rozas et al. 2006a,b; Firpo et al. 2011; Hägele et al. 2010). The results depicted in Fig. 8 turn out to be quite encouraging. In the lower left quadrant, the models are surprisingly successful in covering the spread in object positions to the right of the equal temperature (dashed) line. As the value of the sequences’ metallicity decreases, the calculated temperatures move closer to the dashed line. A discussion on the behavior of T_O   iii in combined shock-photoionization models can be found in Appendix B.

Despite the relative success of our calculations, it should be noted that they represent an oversimplification of a very complex problem. A subsequent analysis with 3D hydrodynamical simulations is necessary. As stated earlier, the line intensities and the temperatures calculated in the models shown in Fig. 8 correspond to the zone that is shock-heated. That is, they include line emission from the shock discontinuity down to the point where cooling and heating (by the UV flux) become equal again, at which point the calculations stop⁸. As indicated in Sect. 3.5, the total column density of an isobaric ionization-bounded photoionization model with U = 0.01 is 32 times that of the shocked layer of the 56 km   s^-1 model from the Z = 1.0   Z_⊙ sequence (last model of sequence No. 4 in Table 1). This is a concern, since it appears to suggest that one might need 32 successive shocks propagating across the whole nebula for the resulting emission lines to reflect the presence of shocks at the level shown in Fig. 8, a rather implausible proposition. Instead of moving inward, we assumed a possibly more realistic scenario in which the shocks are moving outward as a result of a stellar wind generated by the ionizing stars. In this case, the following scenario emerges. A high Mach number shock results in a compression factor of  ~4 of the pressure (across the shock front) and, from that point on, the density increases because the pressure is conserved across the adiabatic flow, as the gas cools and the temperature decreases. At the point downstream where photoheating balances cooling again, the density has reached a value of about 4T_post/T_pre, where T_pre is an estimate of the local equilibrium temperature. This amounts to a compression factor of  ~28. Hence the effective ionization factor near the back-photoionized zone is  ~3.5  ×  10^-4, and even if this layer was kept photoionized, its [O iii]λ5007 emission contribution could not be significant. As the shock wave progresses outward within the nebula, the accumulating denser gas shell behind the shock will eventually quench the photon flux responsible for photoionizing the preshock zone ahead. The latter would quickly evolve into a low-excitation fossil nebula (ahead of the shock) without significant emission in [O iii]λ5007. If we consider the back of the outward-moving shock (facing the incoming ionizing radiation), we may expect the accumulating dense shell to break up into filaments as a result of instabilities.

Another complication that our simple plane-parallel calculations does not consider is that Hii regions are characterized by a low volume filling factor (Osterbrock 1989). The fine spatial structure of Hii regions has so far not been fully resolved. Although progress has been significant, we have not yet reached the microstructure scale of nebular gas (e.g. O’Dell et al. 2003; Mesa-Delgado et al. 2008; Mesa-Delgado & Esteban 2010). It has been proposed that nebulae consist of gas condensations or filaments close to being ionization-bounded (e.g. Tsamis & Péquignot 2005). The shock scenario described above can in principle be transposed to the case of dense condensations. 3D hydrodynamical models, such as the one developed by Raga et al. (2008), would be required to simulate such a case and to determine the temperature of the diffuse gas in which the large-scale shock surrounding the condensations (or the filaments) propagates.

To summarize, despite the caveats just discussed, we consider that shock-photoionization models have the potential of reproducing the observed spread in temperatures observed below the equal temperature line. These models show a tendency quite similar to that of abundance inhomogeneity sequences presented in Sect. 4.3, since T_O   iii approaches T_S   iii as the gas metallicities decrease (toward the upper right quadrant). This behavior is the opposite of that observed with temperature inhomogeneities as modeled with t² parametrization (Sect. 4.2).

Fig. 9 Recombination temperature as a function of TO   iii for various models. Along each sequence, only the metallicity varies. Solid and dashed lines are used to represent . A dotted line of the same color depicts , but it is plotted only for cases where differs from by more than 3%. The ionization parameter is always U = 0.01. In the case of the 50 km   s-1 combined photoionization-shock model sequence (light green), Z assumes the values of 0.01 Z⊙, 0.10 Z⊙, 0.25 Z⊙  0.49 Z⊙  1.0 Z⊙, and 1.6 Z⊙ (the higher the metallicity, lower is Trec ). All other sequences cover the same metallicity range of 0.01 Z⊙ to 2.5 Z⊙, with a fixed increment of 0.2 dex between successive models. Three sequences represent κ-distributions (20, 10, or 5, see legend) and two sequences represent the case of temperature inhomogeneities (t2 = 0.05 or 0.10, see legend). The homogeneous pure photoionization sequence is represented by the magenta line (same as in Fig. 2).

4.6. Recombination temperatures and the Balmer decrement

4.6.1. Behavior of the recombination temperatures

The abundance discrepancy factor can be defined as the ratio, or the logarithmic difference, between abundances of a same ion derived from ORLs and CELs. The ADF is a result of the large difference in the dependence of CELs and ORLs on temperature. To evaluate how the different models presented so far can account for the ADFs reported in the literature, we computed the recombination temperature characterizing each model. To reflect the fact that recombination line emissivities decrease with increasing temperature, we averaged the quantity for each selected ion X⁺ⁱ using a weight across the nebular structure that is proportional to the product of this ion density with the electron density. We therefore define the average recombination temperature of species X⁺ⁱ as . We derived in this way the three recombination temperatures , and , corresponding to the ions H⁺, He⁺ and O⁺⁺, respectively. We adopted α =  −0.83 (Eq. (6)) for H, He and  − 0.87 for O. Since in any model and do not differ by more than 30 K, the temperature will not be discussed further.

To characterize the behavior of as a function of T_O   iii , we present various sequences of models (solid lines) in Figs. 9 and 10. The metallicity is the only parameter that varies along any of the sequences. A dotted line is used when plotting , but it is only shown for cases where one model of the sequence has differing from by more than 3%. The behavior of in the case of κ-distribution (with κ = 20, 10 and 5) is shown in Fig. 9 together with the case of t² models (with t² = 0.05 or 0.10). The Z sequence for the combined shock-photoionization case with V_s ≃ 50   km   s^-1 is overlaid as well. Interestingly, the t² sequences show a behavior similar to those of κ-distributions. Since t² = 0.10 is close to the highest value reported by PG12, and since the corresponding sequence (purple line) overlaps the κ = 10 sequence, we infer that values as low as κ ≃ 5 should not be required to reproduce even the largest ADFs observed in some Hii regions. NDS12 came to the same conclusion using different temperature diagnostics.

Fig. 10 Recombination temperatures as a function of TO   iii for models that consider metallicity inhomogeneities as discussed in Sect. 4.3. Along each sequence, only the metallicity varies (see legend for values of ZA and θ). In no model does differs from by more than 3%. All four sequences from Fig. 5 with varying contrast θ are included as well as the seven sequences from Fig. 6 with θ = 1.0 (color–coding and nomenclature remain the same, see legend). An open circle indicates the position where wi = 0.5. The ionization parameter is U = 0.01 throughout.

All models representing metallicity inhomogeneities previously shown in Figs. 5 and 6 are displayed in Fig. 10 (see the figure line color legend for values of Z_A and θ). Wherever substantially exceeds T_O   iii , the models imply a large ADF, which happens only with the higher metallicity sequences. A pertinent result from Fig. 10 is that for Z_A ≲ 0.05 sequences, the models diverge very little from the equal temperature gray dashed line. Moreover, the upward curvature⁹ above the equal temperature line (sequences Nos. 8–11 in Fig. 10) is opposite to that of the higher Z sequences, since in sequences Nos. 8–11 lie above T_O   iii . These models become unsuitable to account for any ADF among metal-poor nebulae. If we increase the contrast θ much beyond the assumed value of 0.1 dex, we find that it is possible to reverse the sign of the curvature. For instance, for the navy sequence with Z_A ≲ 0.01, the metallicity contrast θ had to be increased to 2.4  dex  for to become again lower than T_O   iii . Even though obtaining  is always possible with sufficiently high θ values, targeting metal-poor nebulae comes with the proviso that only low values of w_i are allowed. It appears to us that such a fine-tuning in the parameter space is an implausible requirement. Unlike the alternate mechanisms presented in previous Fig. 9, metallicity inhomogeneities cannot account for sizable ADFs in nebulae that have T_O   iii ≥ 15   000 K, that is, in Hii regions that are metal-poor. Yet, in their recalibration of the Pagel method, PG12 reported five Hii regions beyond  > 15   000 K. These nebulae are characterized by an equivalent t² in the range 0.02–0.12, indicative of a temperature markedly lower than T_O   iii .

Fig. 11 Calculations of the Balmer decrement as a function of metallicity for various models. The thick light gray line is a metallicity sequence representing the homogeneous pure photoionization case with U = 0.01 (same sequence as in Fig. 2). The (short vertical) black line shows a U-sequence at a fixed metallicity Z = 0.01   Z⊙ (see cyan dotted line in Fig. 2). Three metallicity sequences representing the κ-distribution case are shown. The corresponding κ value appears in the line color legend. The light green line represents a metallicity sequence for a 50 km   s-1 shock combined with photoionization. Two sequences (labeled 7 and 8 in the legend) are shown that represent models with metallicity inhomogeneities. The values of ZA (in Z⊙ units) and θ (in dex) appear in the line color legend. An open circle indicates the position where wi = 0.5. Yellow dashed lines represent the low-density recombination case-B Balmer decrement at three different temperatures (Osterbrock 1989).

4.6.2. Enhanced Balmer decrement and dereddening

The reddening correction that is routinely applied to observed line intensities can be a concern, even if the value of the Balmer decrement adopted for determining A_V is apparently consistent with the subsequently derived [O iii] temperature. In particular, an excessive dereddening of the auroral and nebular line intensities would result in overestimates of the T_O   iii and other CEL temperatures when these auroral lines occur at shorter wavelengths than the corresponding nebular lines¹⁰. This inconsistency can have consequences for the case of κ-distributions or large temperature inhomogeneities, since the intrinsic Balmer decrement is not set by T_O   iii , but rather arises from a plasma at a considerably lower temperature, T_rec . For instance, the intrinsic Balmer decrement can approach values  ≃3 when the abundances are above solar and T_rec ≪ T_O   iii. A better understanding of temperatures in Hii regions may require that we relate the process of dereddening with temperature determination. To illustrate our concern, we carried out a comparison of the Balmer decrement from models that reproduce the spread in T_O   iii . The results are shown¹¹ in Fig. 11 where the calculated Balmer decrement is plotted as a function of nebular metallicity. Two metallicity inhomogeneity sequences with θ = 1   dex are plotted as a function of average metallicity defined as . These correspond to a model A with Z_A = 0.04 Z_⊙ (purple line) and Z_A = 0.16 Z_⊙ (dark gray line), respectively. We note that they do not depart appreciably from the homogeneous photoionization case. The combined photoionization-shock models show evidence at low metallicities of a slightly higher decrement. For κ-distributions, the lower the adopted κ value, the higher the decrement toward low metallicities.

We recall that the lower left quadrant in our temperature diagram contains objects whose metallicities are likely to be higher than 0.2 Z_⊙ (see Sect. 3). Yet, for our data set, the objects with the lowest T_S   iii values can have corresponding T_O   iii temperatures that exceed 10 000 K. At metallicities  ~1–2 Z_⊙, which κ-distribution models favor to reproduce the strongest excesses in T_O   iii (red line in Fig. 7), the Balmer decrement predicted by these models in Fig. 11 is nearing or even exceeding 3. Hence, when comparing data with models that predict different behaviors of CEL temperatures, it would make sense to first apply a reddening correction consistent with the Balmer decrement that these same models predict.

5. Conclusions

An analysis of 102 pairs of [O iii] and [S iii] temperatures of Hii galaxies, giant extragalactic Hii regions, Galactic Hii regions, and Hii regions from the Magellanic Clouds was carried out. The [O iii] temperatures appear to be higher than the corresponding values derived from [S iii] in objects with T_S   iii lower than 14 000 K. This trend is present in objects with metallicities Z ≳ 0.2. For the coolest nebulae (the highest metallicities) the [O iii] temperature excess can reach +3000 K. The results from grids of models that attempt to account for this excess can be summarized as follows:

• 1.

  Simple photoionization calculations result in essentially equal [O iii] and [S iii] temperatures in the range of interest (T_S   iii < 14 000   K). Varying either the ionization parameter U, the metallicity Z, or the hardness of the ionizing SED over a wide range does not alter this result.

• 2.

  Including temperature inhomogeneities of large mean squared amplitude (t²) does not result in values of T_O   iii higher than T_S   iii in the quadrant of reliable data (T_S   iii < 14 000   K).

• 3.

  Metallicity inhomogeneities can successfully reproduce the observed excess in T_O   iii temperatures. By combining two photoionization models of widely differing metallicities (for instance a model with Z = 0.2   Z_⊙ with another with Z = 1.26   Z_⊙), most of the observed spread in temperatures can be reproduced.

• 4.

  Adopting a κ-distribution for the electron energy distribution results in models that successfully reproduce the observed excess in T_O   iii temperatures. Values of κ in the range  ≳ 5 to 20 are favored.

• 5.

  Shock waves that propagate in the photoionized gas can also account for the observed excess in T_O   iii . The 1D models that we computed are simplistic and cannot be considered a fully autoconsistent description of shock waves propagating in a chaotic medium. 3D hydrodynamical calculations would be required that include proper radiative transfer in the context of outward-propagating shocks.

• 6.

  Both the shocked nebula model and the photoionization model with abundance inhomogeneities share the property of a diminishing excess in T_O   iii with decreasing nebular metallicity, while the κ-distributions result in a near constant T_O   iii –T_S   iii offset when κ is set to a fixed value. Through improved observations, it would be important to determine the behavior of the temperatures in the top right quadrant of the T_S   iii vs. T_O   iii diagram.

• 7.

  We derived the recombination temperatures for the various models. Models that consider κ-distributions or temperature inhomogeneities t² result in recombination temperatures that are substantially lower than T_O   iii over the whole metallicity range. On the other hand, models with metallicity inhomogeneities have higher than T_O   iii in the upper right quadrant where T_O   iii ≳ 15   000 K. Metallicity inhomogeneities are therefore only applicable to metallicity-enriched nebulae and could not account for any sizable ADF in metal-poor nebulae (Z < 0.2 Z_⊙).

Acknowledgments

L.B. acknowledges support from CONACyT grant CB-128556. R.M. received financial support from project IN105511 from PAPIIT (UNAM). We are thankful to Jana Benda Klouda for providing us with the initial translation in English. We thank Diethild Starkmeth for proofreading the article. We acknowledge that many suggestions from an anonymous referee contributed to improve the usefulness and scope of the final version.

References

Appendix A: Sample abundance calibration

Figure A.1 shows the relationship between the electron temperatures inferred from the [O iii] lines and the total oxygen abundances 12 + lg(O/H). This sample was built from the literature and is a collection of objects for which measurements of the nebular and auroral emission lines of [O ii] and [O iii] were available. This sample comprises objects from Pérez-Montero & Díaz (2005) for which the determination of T_O   iii and 12 + lg(O/H) was possible, using the direct method (see Hägele et al. 2008; H08), as well as the Hii galaxies from H06, H08, and H11. The solid line corresponds to a least-squares fitting of the data, taking their individual errors into account:

where t_e([O iii]) is the [O iii] electron temperature in units of 10⁴ K. The O/H determination above is based on CELs only. It likely underestimates the true O/H abundance ratio by a factor  ~2, which would likely have been inferred if the ADF had been determined from ORLs and had been taken into account in the above regression. If the electron energies exhibit a κ non-equilibrium distribution, this is expected to increase the oxygen abundance, as shown by NDS12. Likewise with models that consider temperature or abundance inhomogeneities.

Fig. A.1 Relation between te([O iii]) and total oxygen abundances (12 + lg(O/H)) for Hii galaxies from H06, H08 and H11 (solid circles), and Hii galaxies (open circles), GEHRs (upward triangles) and diffuse Hii regions in the Galaxy and the Magellanic Clouds (squares) from Pérez-Montero & Díaz (2005). The solid line is the actual fit to the data. The temperatures are in units of 104 K.

Appendix B: Behavior of combined shock-photoionization models when V_s is increased

In combined shock-photoionization models the span in temperature and in ionization stages of metals covers a much wider range than with photoionization alone. Furthermore, the temporal evolution of the plasma that characterizes shocks is another significant difference. To analyze why the T_O   iii temperatures increase markedly with V_s , we compare two models from the cyan sequence with Z = 1.0   Z_⊙ in Table 1. When one derives volume-averaged quantities that are weighted by the density square of the ion O⁺² or S⁺², it is found that the 45.3 km   s^-1 and 56.1 km   s^-1 models are characterized by a volume-averaged ionic temperature ⟨T_O++⟩ of 8400 K and 9700 K, respectively, while for ⟨T_S++⟩ the corresponding values are 8070 K and 7980 K, respectively. Hence, not much separates the two models as far as average properties are concerned. In the 45.3 km   s^-1 model, the mean doubly ionized fraction for oxygen ⟨n_O++/n_O⟩ is 0.164 while the equivalent fraction for sulphur ⟨n_S++/n_S⟩ is 0.889. In the 56.1 km   s^-1 model, the mean doubly ionized fractions are 0.087 and 0.844, respectively. Photoionization is responsible for the small reverse effect of a reduction in ⟨n_O++/n_O⟩ because the mean density ⟨n_H + ⟩ actually increases from 21.1  cm^-3 to 32.1  cm^-3, (between the 45.3 km   s^-1 and the 56.1 km   s^-1 model, respectively), which has the effect of reducing the effective ionization parameter. We can conclude that combined photoionization-shock models in this velocity regime do not drastically alter the mean degree of ionization of oxygen or sulphur, or their mean temperatures where the doubly ionized ions reside. Indeed, both shock models start with the same ionization fraction at the shock front, of 0.93 for n_O++/n_O and 0.49 for n_S++/n_S, but they possess quite different postshock temperatures (35 000 K and 50 400 K, respectively).

To summarize, there are two main factors that account for the behavior of T_O   iii in Fig. 8: (1) ⟨n_O++/n_O⟩ ≪ ⟨n_S +  + /n_S⟩, which causes the much hotter region near the shock front to contribute proportionally more to the integrated line spectrum of the [O iii] lines than it does for the [S iii] lines, and (2) the strong dependance of the CELs on the factor  ∝ T^-0.5exp( − ΔE/kT) is responsible for favoring the copious emission of lines that have the highest energy transitions ΔE at the head of the shock where T ≲ T_post, as is the case for [O iii]λ4363. This is the main factor that accounts for the high 4363   Å/5007   Å line ratio, and hence for the higher values of T_O   iii . [S iii]λ6312 on the other hand does not beneficiate much from either factor.

The decrease of both T_O   iii and T_S   iii at low velocities in sequences of very low-metallicity gas (light green line and red line sequences in Fig. 8) reflects the increase in the mean density as the postshock temperature is successively increased, which causes a reduction of the effective ionization parameter and hence an increase in the neutral fraction H⁰/H, which favors an increase in the cooling due to collisional excitation of Hi. This increased cooling in turn decreases ⟨T_O++⟩ and ⟨T_S++⟩. A similar mechanism is at work at very low metallicities in U -sequences of homogeneous photoionization models (e.g. cyan dotted line in Fig. 2), as discussed in Sect. 4.1.

All Tables

All Figures

Fig. 1 Relation between TS   iii and TO   iii temperatures inferred from the [S iii] and [O iii] auroral-to-nebular line ratios for the objects studied by PM06, including the three Hii galaxies later observed by H06. The light gray dashed line depicts the locus of equal temperatures (TS   iii = TO   iii), while the solid line indicates the error-weighted least-squares fit to the data (see Sect. 2).

In the text

Fig. 2 Photoionization calculations with an SED corresponding to a zero-age instantaneous starburst model from M03 with Msup = 100    M⊙. A sequence of models of increasing metallicity Z is represented by the magenta solid line, for which U is constant at 0.01. The metallicities covered a range from Z = 0.01   Z⊙ to Z = 2.5   Z⊙, where Z increases by 0.2 dex from model to model. Three sequences of models are superimposed in which U is increased (rather than Z), covering the range 0.001 to 0.046. Each sequence has a different metallicity, as follows: Z = 0.01   Z⊙ (cyan dotted line), Z = 0.25   Z⊙ (light green dotted line) and Z = 1.0   Z⊙ (blue dotted line). In this and subsequent figures, a light gray dashed line represents the locus of equal temperatures.

In the text

Fig. 3 Comparison of two sequences in metallicity that differ by the shape of their ionizing SED: a single temperature star of either Teff = 55   000 K (light green line) or Teff = 40   000 K (blue dotted line). The previous zero-age instantaneous starburst SED of M03 with Msup = 100    M⊙ appears underneath (magenta line). Each sequence comprises 13 models with the metallicity increasing by a factor 0.2 dex from model to model, starting at Z = 0.01   Z⊙ and ending at 2.5 Z⊙. The ionization parameter is U = 0.01 in all calculations.

In the text

	Fig. 4 Sequences of models in which t2 is increased from 0.004 to 0.16 by successive factors of 1.58 ( = 0.2 dex). All calculations were carried out with U   = 0.01 and the zero-age instantaneous starburst model SED with Msup = 100      M⊙ from M03. Each sequence corresponds to a different metallicity Z, which is indicated in the line color legend in Z⊙ units.
In the text

Fig. 5 Sequences of models in which the relative weight wi of two photoionization calculations is varied, as described in Sect. 3.3. For each sequence, the metallicity of its model A in solar units appears in the line color legend followed by the contrast θ in dex. An open circle indicates the position where the relative weight applied to model A and B is equal (wi = 0.5). In all the above models, the ionization parameter is U = 0.01 and the SED corresponds to the zero-age instantaneous starburst model with Msup = 100      M⊙.

In the text

Fig. 6 Seven sequences of dual metallicity models covering a wider range in metallicities. The metallicity contrast between ZB and ZA is fixed at θ = 1.0 dex. Each sequence corresponds to a different metallicity ZA value, which is indicated in the line color legend in Z⊙ units. An open circle along each sequence indicates the position where the relative weight applied to each model is equal (wi = 0.5). In all the above sequences, the ionization parameter is U = 0.01, and the SED corresponds to the zero-age instantaneous starburst model with Msup = 100      M⊙.

In the text

Fig. 7 Sequences of models in which the parameter κ has the values of 80, 40, 20, 10, and 5. The larger κ, the closer the model to standard models with a Maxwell-Boltzmann electron energy distribution. These κ sequences differ in their metallicity Z (expressed in Z⊙ units in the line color legend). The dashed magenta line represents a metallicity sequence of standard photoionization models (κ = ∞) (same sequence as in Fig. 2). The light gray arrow illustrates the increase in the TS   iii and TO   iii temperatures when enhanced Hi excitation is left out from the calculations of the model with Z = 0.01 Z⊙ and κ = 5. In all calculations, we adopted U = 0.01 and an SED consisting of a zero-age instantaneous starburst model from M03 with Msup = 100    M⊙.

In the text

Fig. 8 Six sequences of combined shock-photoionization models along which the shock velocity Vs is increased in locked steps. The preshock conditions are set by photoionization equilibrium, using the parameter U = 0.01 and the M03 SED, and assuming thermal and ionization equilibria (see Sect. 3.5). Each sequence is characterized by a different metallicity Z value, which is expressed in Z⊙ units in the line color legend. For each model, the preshock temperature, the postshock density, the postshock temperature, and the resulting shock velocity are given in Table 1. The first two models of the Z = 0.5 Z⊙ and 1.0 Z⊙ sequences overlap too closely to be distinguished.

In the text

Fig. 9 Recombination temperature as a function of TO   iii for various models. Along each sequence, only the metallicity varies. Solid and dashed lines are used to represent . A dotted line of the same color depicts , but it is plotted only for cases where differs from by more than 3%. The ionization parameter is always U = 0.01. In the case of the 50 km   s-1 combined photoionization-shock model sequence (light green), Z assumes the values of 0.01 Z⊙, 0.10 Z⊙, 0.25 Z⊙  0.49 Z⊙  1.0 Z⊙, and 1.6 Z⊙ (the higher the metallicity, lower is Trec ). All other sequences cover the same metallicity range of 0.01 Z⊙ to 2.5 Z⊙, with a fixed increment of 0.2 dex between successive models. Three sequences represent κ-distributions (20, 10, or 5, see legend) and two sequences represent the case of temperature inhomogeneities (t2 = 0.05 or 0.10, see legend). The homogeneous pure photoionization sequence is represented by the magenta line (same as in Fig. 2).

In the text

Fig. 10 Recombination temperatures as a function of TO   iii for models that consider metallicity inhomogeneities as discussed in Sect. 4.3. Along each sequence, only the metallicity varies (see legend for values of ZA and θ). In no model does differs from by more than 3%. All four sequences from Fig. 5 with varying contrast θ are included as well as the seven sequences from Fig. 6 with θ = 1.0 (color–coding and nomenclature remain the same, see legend). An open circle indicates the position where wi = 0.5. The ionization parameter is U = 0.01 throughout.

In the text

Fig. 11 Calculations of the Balmer decrement as a function of metallicity for various models. The thick light gray line is a metallicity sequence representing the homogeneous pure photoionization case with U = 0.01 (same sequence as in Fig. 2). The (short vertical) black line shows a U-sequence at a fixed metallicity Z = 0.01   Z⊙ (see cyan dotted line in Fig. 2). Three metallicity sequences representing the κ-distribution case are shown. The corresponding κ value appears in the line color legend. The light green line represents a metallicity sequence for a 50 km   s-1 shock combined with photoionization. Two sequences (labeled 7 and 8 in the legend) are shown that represent models with metallicity inhomogeneities. The values of ZA (in Z⊙ units) and θ (in dex) appear in the line color legend. An open circle indicates the position where wi = 0.5. Yellow dashed lines represent the low-density recombination case-B Balmer decrement at three different temperatures (Osterbrock 1989).

In the text

	Fig. A.1 Relation between te([O iii]) and total oxygen abundances (12 + lg(O/H)) for Hii galaxies from H06, H08 and H11 (solid circles), and Hii galaxies (open circles), GEHRs (upward triangles) and diffuse Hii regions in the Galaxy and the Magellanic Clouds (squares) from Pérez-Montero & Díaz (2005). The solid line is the actual fit to the data. The temperatures are in units of 104 K.
In the text