© ESO, 2015

1. Introduction

In the present phase of the universe the number and mass fraction of hot massive stars is small compared to the total stellar population; however, the impact of these objects on their environment is decisive for the evolution of the host galaxies. For instance, metals are produced in the cores of massive stars, and in the later evolutionary stages of these objects these metals are distributed via stellar winds in the late phases and supernova explosions into the surrounding interstellar gas. This mechanism also influences the subsequent star formation in these environments (cf. Woosley & Weaver 1995; François et al. 2004; Maio et al. 2010; Hirschmann et al. 2013; Sandford et al. 1982; Oey & Massey 1995; Bisbas et al. 2011). As massive stars have short lifetimes of only a few million years, their locations and their chemical compositions are still closely correlated with the environment where they formed. They can thus be used as tracers for the chemical states and metallicity gradients of galaxies. Furthermore, in disk galaxies hot massive stars are the most important emitters of ionizing radiation. They act as the primary energy supply sources for H ii regions and they are involved in the process of energy maintenance required for the dilute, but extended, diffuse ionized gas (cf. Haffner et al. 1999; Rossa & Dettmar 2000; Hoffmann et al. 2012). Thus, not least with respect to their luminosities of up to L ~ 10⁸L_⊙ (Pauldrach et al. 2012), hot massive stars are important players in the evolution of massive star clusters, especially in starburst clusters and galaxies. The significance of these starbursts is not only that they heat the intergalactic medium and enrich it with metals, but that massive stars themselves are also useful diagnostic tools. Stellar wind lines can presently be identified in the spectra of individual O-supergiants up to distances of 20 Mpc. Because of their high star formation rates and the resulting large number of young stars, starburst galaxies show the distinctive spectral signatures of hot massive stars in the integrated spectra of starburst galaxies even at redshifts up to z ~ 4 (Steidel et al. 1996; Jones et al. 2013). These spectral features can provide important information about the chemical composition, the stellar populations, and thus the galactic evolution even at extragalactic distances. That this is feasible has already been demonstrated by Pettini et al. (2000), who showed that the O-star wind line features can indeed be used to constrain star formation processes and the metallicity. By comparing theoretical population synthesis models with observational results, it is possible to reconstruct the physical properties and the recent evolution of these objects (Leitherer & Heckman 1995; Leitherer et al. 1999, 2010).

We are thus close to the point of making use of complete and completely independent quantitative spectroscopic studies of the most luminous stellar objects. To realize this objective a diagnostic tool for determining the physical properties of hot stars via quantitative UV spectroscopy is required. The status of the ongoing work to construct such an advanced diagnostic tool that includes an assessment of the accuracy of the determination of the parameters involved has recently been described by Pauldrach et al. (2012). They have shown that the atmospheric models developed for massive stars are already realistic with regard to quantitative spectral UV analysis calculated along with consistent dynamics, which allows the stellar parameters to be determined by comparing an observed UV spectrum to a set of suitable synthetic spectra. The astronomical perspectives are enormous, not only for applications of the diagnostic techniques to massive O-type stars, but also extended to the role of massive stars as tracers of the chemical composition and the population of starbursting galaxies at high redshift. In all of these applications a close connection between observations and accurate theoretical modeling of the stellar spectra is required, and this is also the case for the physical properties and observational features of H ii regions, which are closely connected with the spectra and luminosities of their ionizing sources: hot massive stars (cf. Rubin et al. 1991; Sellmaier et al. 1996; Giveon et al. 2002; Pauldrach 2003; Sternberg et al. 2003).

The analysis of H ii regions, which is primarily based on emission line diagnostics, is a powerful method with which to gain information about the chemical properties and evolution of galaxies. Because the strength of the emission lines depends on the properties of the gas (density, chemical composition) and the properties of the sources of ionization (luminosity, spectral energy distribution), quantitative analyses of these lines are used to determine abundances in H ii regions in galaxies of different metallicities and to draw conclusions about the spectral energy distributions (SEDs) of the irradiating stellar fluxes and thus about the upper mass range of the stellar content of the clusters. High-quality far-infrared spectra of extragalactic H ii regions taken with the Spitzer Space Telescope are discussed and interpreted, for example, by Rubin et al. (2007, 2008), work is also in progress to determine abundance gradients within single galaxies, e.g., by Márquez et al. (2002), Stanghellini et al. (2010).

One of the most important ingredients of galaxies, clusters, and structures to be considered and to be determined in this context is obviously the metallicity. Not only are the properties of stars and nebulae influenced to a large extent by the metal abundances, but the metals are in turn also provided locally and their distribution is modified by the physical behavior of the stellar content. Thus, it is not surprising that the metallicities found in the interstellar medium differ considerably even within groups of galaxies. The metallicity found in the Large Magellanic Cloud, for example, has a value of Z_LMC ≈ 0.4 Z_⊙, whereas the metallicity in the Small Magellanic Cloud has been determined to Z_SMC ≈ 0.15 Z_⊙ (Dufour 1984). Both values are significantly smaller than the metallicity in the Milky Way which seems to correspond roughly to the solar value. However, the term “solar abundance” itself is somewhat vague, since studies in the past decade have surprisingly yielded significantly lower number fractions of the most abundant metals such as C, N, O, and Ne (Asplund et al. 2009) than the previously determined values (e.g., Grevesse & Sauval 1998). Other examples of large differences to the galactic metallicity in the local universe are M 83 (Z ≈ 2 Z_⊙, Bresolin & Kennicutt 2002) and blue low surface-brightness galaxies (Z ≈ 0.1 Z_⊙, Roennback & Bergvall 1995). Beyond that, metallicities vary not only from galaxy to galaxy, they also vary within the disks of certain galaxies. Based on observations of H ii regions at increasing distances from the center of our Galaxy, Rudolph et al. (2006) have found a decline of the oxygen abundance, for instance.

In order to make further progress in the diagnostics of the emission line spectra of H ii regions one obviously has to account for the metallicity-dependence of the calculated SEDs of massive stars. One of the most important obstacles in this regard is the fact that massive stars show direct spectroscopic evidence of winds throughout their lifetime, and these winds modify the ionizing radiation of the stars considerably (cf. Pauldrach 1987) and contribute significantly to the state and energetics of the atmospheric structures in a metallicity-dependent way. Modeling hot star atmospheres is complicated by the fact that the outflow dominates the physics of the atmospheres, in particular regarding the density stratification and the radiative transfer, which are substantially modified through the presence of the macroscopic velocity field. In the frame of a consistent treatment of the hydrodynamics, the hydrodynamics influence the non-LTE model¹, and are in turn controlled by the line force determined by the occupation numbers and the radiative transfer of the non-LTE model (cf. Pauldrach 1987, 1990, 1994, 2001; Pauldrach et al. 2012).

One of the objectives of this paper is to present a grid of advanced stellar wind spectra for O-type dwarfs and supergiants at different metallicities (Sect. 2) computed using hydrodynamic atmospheric models that include a full treatment of non-LTE line blocking and blanketing² and the radiative force. In the second part of this paper the influence of the computed SEDs on the properties of the irradiated interstellar gas is differentially investigated. In Sect. 3 we apply our computed SEDs to a series of simulations of sample H ii regions in order to investigate the dependence of the temperature and ionization structures on the ionizing spectra and the metallicity of the gas. At first we restrict the simulations to H ii regions which are illuminated by a single star and which consist of 1D (i.e., spherically symmetric) homogeneous density structures. These restrictions are dropped in the second part of this section, where the effects of multiple radiative sources and inhomogeneous density structures on H ii regions are investigated and discussed using our recently developed 3D radiative transfer models. In Sect. 4 we will summarize our results along with our conclusions.

2. Stellar wind models for O-type dwarfs and supergiants at different metallicities

In this section we apply our method for modeling the expanding atmospheres of hot stars to a basic grid of massive O-stars. The objective of these calculations is to present ionizing fluxes and SEDs for massive dwarfs and supergiants at different metallicities that can be used for the quantitative analysis of emission line spectra of H ii regions.

2.1. The general concept for calculating synthetic spectra and SEDs of massive stars

Our approach to modeling the expanding atmospheres of hot, massive stars has been described in detail in a series of previous papers (Pauldrach 1987; Pauldrach et al. 1990, 1993, 1994, 1998, 2001, 2004, 2012; Taresch et al. 1997; Haser et al. 1998), and we summarize the salient points here. Our method is based on the concept of homogeneous, stationary, and spherically symmetric radiation-driven atmospheres. Although this is an approximation to some extent, it is sufficient to reproduce all important characteristics of the expanding atmospheres in some detail.

A complete model atmosphere calculation consists of (a) a solution of the hydrodynamics describing velocity and density of the outflow, based on radiative acceleration by Thomson, continuum, and line absorption and scattering (an essential aspect of the model, because the expansion of the atmosphere alters the emergent flux considerably compared to a hydrostatic atmosphere); (b) determination of the occupation numbers from a solution of the rate equations containing all important radiative and collisional processes, using sophisticated model atoms and corresponding line lists³; (c) calculation of the radiation field from a detailed radiative transfer solution taking into account not only continuum, but also Doppler-shifted line opacities and emissivities⁴; and (d) computation of the temperature from the requirement of radiative (absorption/emission) and thermal (heating/cooling) balance. An accelerated Lambda iteration (ALI) procedure⁵ is used to achieve consistency of occupation numbers, radiative transfer, and temperature. If required, an updated radiative acceleration can be computed from the converged model, and the process repeated.

In addition, secondary effects such as the production of EUV and X-ray radiation in the cooling zones of shocks embedded in the wind and arising from the nonstationary, unstable behavior of radiation-driven winds can, together with K-shell absorption, be optionally considered (based on a parametrization of the shock jump velocity; cf. Pauldrach et al. 1994, 2001). However, they have not been included in the models presented here, since they affect only high ionization stages like O vi which are not relevant for the analysis of emission line spectra of H ii regions.

Of course, it needs to be clarified whether the SEDs calculated by our method are realistic enough to be used in diagnostic modeling of H ii regions. Although the radiation in the ionizing spectral range cannot be directly observed, the predicted SEDs can be verified indirectly by a comparison of observed emission line strengths and those calculated by nebular models (cf. Sellmaier et al. 1996; Giveon et al. 2002; Rubin et al. 2007, 2008). A more stringent test can be provided by a comparison of the synthetic and observed UV spectra of individual massive stars, which involves hundreds of spectral signatures of various ionization stages with different ionization thresholds, and covering a large frequency range: because almost all of the ionization thresholds lie in the spectral range shortward of the hydrogen Lyman edge (cf. Pauldrach et al. 2012), and the ionization balances of all elements depend sensitively on the ionizing radiation throughout the entire wind, the ionization balance can be traced reliably through the strength and structure of the wind lines formed throughout the atmosphere. In this way a successful comparison of observed and synthetic UV spectra (Pauldrach et al. 1994, 2001, 2004, 2012) ascertains the quality of the ionization balance and thus of the SEDs.

2.2. Synthetic spectra and SEDs from a model grid of massive stars

Our model grid comprises massive stars with effective temperatures ranging from 30 000 to 55 000 K and luminosities from 10⁵L_⊙ to 2.2 × 10⁶L_⊙ (Table 1). The model parameters correspond to those used by Pauldrach et al. (2001) and were chosen in accordance with the range of values deduced from observations.

We compute stellar models for metallicities of 0.1, 0.4, 1.0, and 2.0 Z_⊙ and compare the results for two different data sets for the solar metallicities, the abundances given by Grevesse & Sauval (1998) that had been used by Pauldrach et al. (2001) and the updated values published by Asplund et al. (2009). The latter values were determined using the comparison of the solar spectrum with 3D radiative transfer simulations of the solar photosphere and differ from the former determinations by up to 38% percent for the most abundant metals (see Table 2).

The influence of the metallicity on the winds of hot stars.

As the winds of hot stars are primarily driven by metal lines, the chemical composition of the stellar atmosphere is a decisive factor controlling the strength of the winds. Although larger metallicities lead in general to stronger winds (Kudritzki et al. 1987; Pauldrach 1987; Pauldrach et al. 2012), not all elements act on the winds in the same way because of their different line strength distributions (cf. Pauldrach 1987). For example, an increase in the abundances of the iron group elements (Fe, Ni) will increase the mass-loss rate and correspondingly decrease the terminal velocity, whereas an increase in the abundances of some lighter elements (C, N, O, S, Ne, Ar) will increase the terminal velocity and correspondingly decrease the mass-loss rate. Thus, different wind parameters may be obtained if the abundance ratios are changed, even if the total metallicity is kept constant⁶.

We will, however, not consider such second-order effects in the present work. Instead, our computed mass-loss rates shown in Table 1 are based on the values of the mass-loss rates presented by Pauldrach et al. (2001), scaled by the corresponding factors derived from the metallicity dependence of the wind strengths exhibited by self-consistent hydrodynamical radiation-driven wind calculations (details of this procedure are described by Pauldrach et al. 2012). We have also applied our hydrodynamical radiation-driven wind calculations to a grid of very massive stars (VMS) with initial masses between 150 M_⊙ and 3000 M_⊙, effective temperatures T_eff between 40 000 K and 65 000 K, and metallicities Z of 0.05 Z_⊙ and 1 Z_⊙ (Table 3; the stellar parameters are based on theoretical evolutionary models as described by Belkus et al. 2007; see also Pauldrach et al. 2012).

Synthetic spectra and SEDs obtained from the O-star grid.

As a primary result of our computations, Figs.A.1−A.4 show the emergent ionizing fluxes together with the corresponding shapes of the continuum for the models of our grid. As can be verified from these figures the influence of the line opacities, i.e., the difference between the hypothetical continuum and the real emergent flux, increases from objects with cooler effective temperatures to those with hotter effective temperatures and from dwarfs to supergiants. This result is not surprising since it is a consequence of the behavior of the most important dynamical parameter, the mass-loss rate Ṁ, which is directly coupled to the stellar luminosity. Thus, the mass-loss rate increases with the effective temperature and the luminosity class (Table 1), and the optical depths of the spectral lines, which increase accordingly, produce, along with the increasing mass-loss rate, a more pronounced line blocking effect in the wind part of the atmosphere. As can be seen from Figs.A.1−A.4, this behavior saturates for objects with effective temperatures larger than T_eff = 45 000 K, since in these cases higher main ionization stages are encountered (e.g., Fe v and Fe vi which have a smaller number of bound-bound transitions (cf. Pauldrach 1987). The effect of line blocking is thus strongest for supergiants of intermediate T_eff.

The figures show that the SEDs depend sensitively on the stellar parameters (especially the effective temperature), and also on the metallicity. This influence on the spectra is not only due to the direct line-blocking effect caused by the metals, but also indirectly due to changes in the hydrodynamic structure that occur as a consequence of the influence of the metallicity on the radiative line acceleration.

In Table 4 we list for each model the number of photons emitted per second capable of ionizing H, He, He⁺, O⁺, Ne⁺, and S⁺. The hydrogen ionizing flux determines (along with the density structure and temperature of the gas) the extent of the ionized volume while helium is an important absorber for the hard ionizing radiation. The ionization products of considered ionization stages of metals are effective line radiation emitters in gaseous nebulae and therefore play an important role in the corresponding line diagnostics. Their line emission also significantly contributes the energy balance of the ionized gas. The corresponding ionizing fluxes of stellar sources are therefore decisive for the properties of the gas in their environment.

Although the number of hydrogen-ionizing photons depends only weakly on metallicity in the stellar temperature range of our grid, significant differences of up to several orders of magnitude are found for the ionization stages of He⁺ and Ne⁺. These pronounced anti-correlations between the ionizing fluxes and the stellar metallicity indicate strongly that the influence of stellar metallicity on nebular line strength ratios is of similar importance to that of variations in the nebular metallicity⁷. However, we also note that the relative influence of the stellar metallicity on the number of emitted ionizing photons decreases for larger effective temperatures and hence harder ionizing spectra.

Very massive stars.

As a first indication to what extent very massive stars can be identified by their ionizing influence on the environment, we list in Table 5 the photon emission rates of our very massive star models in the different ionization continua. The SEDs of the hottest of these stars with an effective temperature of 65 000 K and masses of 1000 M_⊙ and 3000 M_⊙ are characterized by ratios of He ii-ionizing to H-ionizing fluxes that exceed those of even the hottest “normal” O-type model stars (the 55 000 K dwarfs) by a factor of 5 to 10, while for the cooler very massive star models the ratios are similar to those of the normal O-star models⁸. We will simulate the effects of the very massive star model SEDs on their surrounding H ii regions in Sect. 3.2.2. Plots of the SEDs of the VMS are shown in Figs. A.5 and A.6. The corresponding data sets can be copied from the links provided there.

3. Application of the calculated SEDs to one- and three-dimensional simulations of H ii regions: influence on the emission line intensities

The ionized gas of H ii regions reemits most of the extreme ultraviolet (EUV) energy output of the star in a limited number of lines in the UV, the optical, and the IR. The comparatively high intensities in these lines allows observing the ionized gas even if the central source is visually much fainter, making it possible to determine the temperatures and metallicities of the ionized gas even of extragalactic star-forming regions (cf. Zaritsky et al. 1994; Moy et al. 2001; Rubin et al. 2008, 2010, 2012; Pilyugin et al. 2013). Additionally, the emission lines of H ii regions can be measured using narrow-band filters allowing for observations which simultaneously provide information about the emission spectrum and the spatial structure of an object. This method has been applied to examine the metallicities of H ii regions as a function of the position within their host galaxies (cf. Cedrés et al. 2012) or to analyze the substructure of the gas in a single H ii region (cf. Heydari-Malayeri et al. 2001).

In this section we investigate the influence of the computed stellar SEDs on the properties of H ii regions with homogeneous and with inhomogeneous density structures. Via simulated sample H ii regions we examine how the temperature and ionization structures in the gas depend on the ionizing spectra and the stellar and nebular metallicity. The question remains, however, to what extent these models may represent real-world H ii regions. The metallicity, the complex geometric structures of H ii regions, and the clumpiness of the medium are usually quantitatively unknown, but may affect the analysis. If discrepancies are encountered, it is often difficult to decide which of the assumptions is responsible for the disagreement found. The computed stellar SEDs themselves can be tested, if the atmospheric models are sufficiently consistent, by their own predicted UV spectra, since the spectral lines and the corresponding ionization fractions in the stellar winds are influenced by the EUV radiation field via the same atomic processes (though at different temperatures and densities) as those that shape the nebular spectra under influence of the emergent stellar flux. In a series of papers we have shown that our consistently calculated synthetic spectra are in sufficiently good agreement with the observed UV spectra of O stars (cf. Pauldrach et al. 1994, 2001, 2004, 2012).

In this work we will therefore study the influence of these effects separately. We begin in Sect. 3.1 with simple homogeneous models with spherical geometry to explore the influence of the stellar SEDs, covering the range of effective temperatures and metallicities given by the grid of stellar models presented in Sect. 2.2.

Next, we will analyze how the assumption of “perfect” spherically symmetric nebulae influences the emission line diagnostics by comparing the results from homogeneous, spherical models to those from 3D models using a fractally inhomogeneous density structure (such as described by Elmegreen & Falgarone 1996; and Wood et al. 2005). The simulations of these structures will be performed using the 3D radiative transfer code described by Weber et al. (2015). This code has been extended for the present paper such that it can account for the ionization structure and the line emission of the most abundant metal ions and the resulting influence of these ions on the energy balance of the H ii regions – for the first time with full time-dependence (see Sect. 3.2.1).

Third, we will compare the gas temperatures and the related line emission of evolving H ii regions with the respective values of steady-state regions. The motivation is that the description of H ii regions as steady-state systems is an approximation which neglects the evolution of star-forming regions (cf. Preibisch & Zinnecker 2007; Murray 2011) and the short lifetimes of the massive stars (cf. Langer et al. 1994) that act as sources of ionization.

Real-world star-forming regions such as the Orion nebula (cf. Muench et al. 2008) and the η Carinae region (cf. Smith 2006) contain not just a single hot star as ionizing source, but several. In the last step of this work we will therefore replace the single-star ionizing sources in the models by clusters of stars. We compute the ionization structure and total emergent fluxes of a single H ii region illuminated by a dense cluster and compare the results with those from a simulation where larger distances between the ionizing sources lead to the formation of partly separate H ii regions. Additionally we will investigate the possibility of finding very massive stars (VMS) in star-forming regions by means of the line emission from the gas ionized by the associations.

As an essential ingredient in the simulations we focus here on the accurate description of the time-dependent ionization structure of the metals. Schmidt-Voigt & Koeppen (1987) had presented a computationally efficient approach where the ionization structure of a given element is determined by interpolating between the initial conditions and the stationary case using a single eigenvalue of the system of rate equations (see below), but this approach is not accurate if more than two consecutive ionization stages of an element have to be considered, because the interpolation does not account for the fact that the transition between two noncontiguous ionization stages requires the creation of ions of the intermediate stages. In the method by Graziani et al. (2013) the radiative transfer and the computation of the occupation numbers are performed consistently for hydrogen and helium. For the occupation numbers of the metal ions and the temperature of the gas pre-computed results from the Cloudy code (Ferland et al. 2013), which describes the stationary states of H ii regions using spherical symmetry, are used. These results are selected from a database such that they match the 3D results for the occupation numbers of H and He ions, and the radiation field.

Our approach on the other hand is to apply a numerical method which treats the ionization stages of all the metals we account for in the same consistent way as hydrogen and helium. The occupation numbers n ≡ n(r,t) of all ionization stages i,j of the elements considered are calculated using the equation of the time-dependent statistical “equilibrium” (Pauli Master Equation, Pauli 1928) $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\ne \mathit{j}}{\mathrm{𝒫}}_{\mathit{j,i}}{\mathit{n}}_{\mathit{j}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}\sum _{\mathit{i}\ne \mathit{j}}{\mathrm{𝒫}}_{\mathit{i,j}}{\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(1)which describes the temporal change of the number density of all ionization stages i, and contains in the rate coefficients all important radiative (ℛ_i,j) and collisional () transition rates.

For the solution of these systems of differential equations we chose an approach that combines integrating the condition of particle conservation within the rate matrices (as described by Mihalas 1978) with providing a robust solution for temporal evolution of the system. To realize this aim we define, following the notation by Marten (1993), a vector x, which contains in its components the number fractions of all considered N ionization stages relative to the total number density of the element and a matrix G′, where the components are defined as ${\mathit{g}}_{\mathit{i,i}}^{\mathrm{\prime }}\mathrm{=}\mathrm{\left(}{\sum }_{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,j}\ne \mathit{k}}^{\mathit{N}}{\mathrm{𝒫}}_{\mathit{i,j}}\mathrm{\right)}\mathrm{+}{\mathrm{𝒫}}_{\mathit{i,k}}$ and ${\mathit{g}}_{\mathit{i,j}\ne \mathit{i}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathrm{𝒫}}_{\mathit{j,i}}\mathrm{+}{\mathrm{𝒫}}_{\mathit{k,i}}$. The fraction x_k(t) of the ionization stage k is replaced by the condition of particle conservation: ${\mathit{x}}_{\mathit{k}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,j}\ne \mathit{k}}^{\mathit{N}}{\mathit{x}}_{\mathit{j}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(2)The resulting inhomogeneous system of differential equations is ${E}\prime \mathrm{·}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\hspace{0.17em}{x}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{G}}^{{\prime }}\mathrm{·}{x}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{b}\mathit{,}$(3)where the components of b are and where all coefficients of the redundant kth row of G′ and b have been replaced by 1, and those of the kth column of E by 0 (with this replacement the unity matrix E becomes E′) – with these numbers inserted the corresponding components represent in total the condition of particle conservation.

For the spherically symmetric models of Sect. 3.1 we focus on the stationary case (where the time-derivative in Eq. (3) disappears) in order to make the results for the H ii regions models using our new stellar SEDs as ionizing sources comparable to the results of other simulations. In Sect. 3.2 we present our results for the time-dependent ionization structure of metals in H ii regions in the context of the description of our 3D radiative transfer code.

3.1. Spherically symmetric models of H II regions

The standard procedure for simulations of emission line spectra of H ii regions is still mostly founded on spherically symmetric models (cf. Stasińska & Leitherer 1996; Hoffmann et al. 2012; and Ferland et al. 2013). They remain useful as a comparative tool because the wavelength region of the ionizing sources blueward of the Lyman edge cannot be observed directly, but the SED of the ionizing flux can nevertheless be studied with such models by means of their influence on the emission line spectra of gaseous nebulae. Such tests are, however, not altogether without uncertainties owing to the additional dependence of the emission line strengths on the chemical composition of the gas in the H ii regions.

Below we outline our numerical approach to investigating these dependencies quantitatively for spherically symmetric model H ii regions. The results obtained with this method are then discussed for a grid of models with different metallicities, using stellar SEDs computed for different temperatures and metal abundances (Sect. 2.2).

3.1.1. The numerical approach applied for the computation of the spherical H II region models

The basic equations describing the temperature and ionization structure of H ii regions are similar to those used to describe non-LTE stellar atmospheres in statistical equilibrium. For the computation of the spherically symmetric models of gaseous nebulae we therefore use a modified version of the WM-basic stellar atmosphere code (Sect. 2.1), which has been adapted to the dilute radiation fields and low gas densities of H ii regions (cf. Hoffmann et al. 2012). This approach yields descriptions of steady-state H ii regions, in which ionization and recombination, as well as heating and cooling, are in equilibrium at every radius point. In such a stationary state Eq. (3) simplifies to ${G}\prime \mathrm{·}{{x}}^{{\infty }}\mathrm{=}{b}\mathit{.}$(4)This equation must be be solved iteratively until it converges to the final value for the stationary state x^∞, because the rate coefficients that define the entries of G′ themselves depend on the occupation numbers x: on the one hand, the recombination and collisional ionization rate coefficients in a gas are proportional to the electron density, which in turn mainly depends on the ionization structure of the most abundant elements hydrogen and helium; on the other hand, the radiative ionization rates are determined by the mean intensity J_ν, which is influenced by the ionization-dependent opacity of the matter between the radiation sources and the considered point in the simulation volume, and by the emissivity of the surrounding gas.

In our approach J_ν is computed using a radiative transfer procedure that is performed along a number of parallel rays intersecting radius shells around the source at different angles, describing the radius- and direction-dependent intensities I_ν. The mean intensity at a given radius is computed by evaluating the integral ${\mathit{J}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{{}^{\mathrm{\int }}}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\mu }$, where μ is the cosine of the angle between a ray and the outward normal at a given radius. The diffuse radiation field created by recombination processes and electron-scattering is treated correctly, avoiding approximations regarding the propagation of photons such as “case B” or “outward-only”.

Apart from the ionization structure, the emission spectrum of an H ii region primarily depends on the temperature of the gas because the recombination and collisional excitation rates are functions of the temperature. The interpretation of observations of H ii regions thus requires an accurate understanding of the microphysical processes that lead to gains and losses of the thermal energy of the gas, which in turn determines the temperature structure. The processes regarded for the computation of the energy balance in H ii regions are heating by photoionization and cooling by radiative recombination, as well as free-free and collisional bound-bound processes.

The low density of the interstellar gas (compared to the gas in stellar atmospheres) leads to small collisional de-excitation rates. Thus radiative transitions of collisionally excited lines are important or even the dominant cooling processes in H ii regions⁹. We extend the modeling described by Hoffmann et al. (2012) to include the cooling rates connected to the forbidden radiative de-excitation processes of collisionally excited substates of the ground levels of C ii, N ii, N iii, O iii, O iv, Ne iii, S iii, and S iv. The cooling by fine-structure transitions is computed by multiplying the collisional transition rates into the excited states with the probability of the corresponding radiative de-excitation processes and the energy of the photons emitted during the relaxation back into the ground state¹⁰.

3.1.2. Dependence of the properties of H ii regions on the ionizing sources

In our systematic series of simulations we examine the ionization and temperature structures and the resulting emission spectra of the gas in homogeneous H ii regions irradiated by single stars. In this series we consider the temperature range of O stars (30 000 K to 55 000 K) and we use the same composition for both nebular and stellar matter¹¹, covering the metallicity range of star-forming regions in the present-day universe. In a second series, we keep the gas at solar metallicity in order to analyze the influence of the metallicity-dependent stellar SEDs (using our 40 000 K dwarf stars) on the temperature and ionization structure of the H ii regions independently from the effects of the metallicity of the gas of the H ii regions.

The dependence of the ionization structure of H  II regions on the metallicity and the stellar SEDs .

In a steady-state H ii region, the number of recombinations in the ionized volume (the Strömgren sphere) must equal the stellar emission rate of ionizing photons. Thus, the sizes of the Strömgren spheres depend on the SEDs of the ionizing stars and on the recombination rates of the ions in the H ii regions.

In the presented grid the radii r_{H ii} of the hydrogen Strömgren spheres grow significantly for lower metallicities, but otherwise equal stellar parameters R, log (g), and T_eff (see Figs. A.7 and A.8). For instance, the hydrogen Strömgren radius for a metallicity of 0.1 Z_⊙ is approximately 50% larger than the Strömgren radius for 2.0 Z_⊙ in the case of the 40 000 K dwarf stars. As shown in Table 4 the hydrogen-ionizing fluxes of O stars are – for otherwise equal stellar parameters – almost independent of metallicity. The different Strömgren radii are therefore primarily a consequence of the different recombination rates, which increase for lower temperatures of the ionized gas as they occur for higher metallicities. This relation results because radiative decays of collisionally excited states of metal ions are the dominant cooling processes in the ionized gas¹². Cooler temperatures in turn increase the recombination coefficients of H ii and other ions (cf. Osterbrock & Ferland 2006) and hence reduce the ionized volumes.

The size of the He ii Strömgren sphere is of particular importance because it marks the boundary where all photons with energies above the He i ionization threshold have been used up, and consequently no significant amounts of metal ions requiring ionization energies above that of He i will be found. Among these ions are N iii, O iii, Ne iii, and S iv. The volumes where helium is ionized are considerably smaller than the hydrogen Strömgren spheres for stellar effective temperatures of T_eff = 30 000 K. The radii of the hydrogen Strömgren spheres exceed the radii of the helium Strömgren spheres by a factor of approximately 3 in the case of the dwarf stars. In contrast, for the supergiants the corresponding factor is approximately 1.5. This results from a lower ratio of He-ionizing to H-ionizing photons for the dwarf models compared to the supergiant models (Table 4). Helium is singly ionized up to the Strömgren radius of hydrogen in H ii regions where the effective temperature of the stellar sources is ≥40 000 K.

The radius range in which helium is fully ionized, i.e., where He iii is the most abundant ionization stage of helium, is small in comparison to the hydrogen Strömgren sphere for dwarf and supergiant O stars. These He iii regions are too small to be resolved in those of our simulations in which the ionizing source is one of the 30 000 K or 40 000 K model stars (with the exception of the supergiant model with 0.1 Z_⊙ where He iii is the most abundant ionization stage for ≈0.05 r_{H ii}). Only the dwarf and supergiant stars with effective temperatures of 50 000 K have appreciable He iii Strömgren spheres that reach ≈0.05 r_{H ii} (for 2.0 Z_⊙) to ≈0.10 r_{H ii} (for 0.1 Z_⊙).

Although the number densities of metal ions are much smaller than the number densities of hydrogen and helium ions, metals have a large impact on the energy balance of the ionized gas, and interpretation of metal line ratios as markers for the galactic evolution (cf. Balser et al. 2011 and the references therein) requires knowledge of the relation between the ionization fractions of metals and the SEDs of the ionizing sources. For example, the O iii/ O ii ratio results from the O ii ionizing flux, which in turn depends not only on the effective temperature of the ionizing sources, but also on their metallicity and their atmospheric density stratification. In the H ii regions around the 30 000 K dwarf stars, O iii is the most abundant ionization stage of oxygen just within less than the innermost ≈0.05 r_{H ii} for all metallicities. The O ii-ionizing fluxes of the supergiants at 30 000 K exceed the fluxes of the dwarf stars of the same effective temperature by approximately 3 dex, which results in a more extended volume in which O iii is the dominant ionization stage of oxygen. The extension of this volume strongly depends on metallicity. Its radius r_{O iii} is ≈0.4 r_{H ii} for 0.1 Z_⊙, but drops to ≈0.1 r_{H ii} for a metallicity of 2.0 Z_⊙.

The O ii-ionizing fluxes of 40 000 K supergiants, which differ by ≈1.3 dex between the model with 0.1 Z_⊙ and 2.0 Z_⊙, show a stronger metallicity dependence than the O ii-ionizing fluxes of the dwarf stars with the same effective temperature, which differ by ≈0.6 dex. Consequently, the radii of the O iii dominated parts of the H ii regions around the supergiants vary more (r_{O iii} ≈ 0.98 r_{H ii} for 0.1 Z_⊙, r_{O iii} ≈ 0.55 r_{H ii} for 2.0 Z_⊙) than the radii around the dwarf stars (r_{O iii} ≈ 0.97 r_{H ii} for 0.1 Z_⊙, r_{O iii} ≈ 0.87 r_{H ii} for 2.0 Z_⊙). The ion O iii is the most abundant ionization stage within the entire Strömgren spheres for dwarf and supergiant stars of all metallicities for an effective temperature of 50 000 K.

Like the O iii fraction relative to the total amount of oxygen, the fraction of S iv decreases for higher metallicities as can be expected in view of the almost identical ionization energies of S iii (2.56 Ryd) and O ii (2.58 Ryd). Still, the relative fraction of S iv is considerably smaller than the relative fraction of O iii within the same H ii region. The reason is the smaller ionization cross section from the ground-state of S iii (0.36 × 10^-18 cm² at the ionization edge) compared to that of O ii (10.4 × 10^-18 cm² at the ionization edge).

The above result that low stellar metallicities lead to harder ionizing spectra and thus to a larger fraction of high ionization stages is in agreement with spectroscopic observations. For example, Rubin et al. (2007, 2008) found rising ⟨ Ne^{2 +} ⟩ / ⟨ Ne⁺ ⟩ and ⟨ S^{3 +} ⟩ / ⟨ S^{2 +} ⟩ ratios for increasing distance from the galactic centers of M 83 and in M 33, based on mid-IR observations with the Spitzer Space telescope. This relation is likely to be connected with the lower metallicities in the outer parts of the galaxies (cf. Rubin et al. 2007). Further observations of the metal-poor galaxy NGC 6822 (Rubin et al. 2012) have found larger fractions of the higher ionization stages than in M 83 (supersolar metallicity) or in M 33 (roughly solar metallicity). These results might, however, additionally be influenced by other factors, such as different stellar mass functions or different effective temperatures of the ionizing stars as a function of the chemical composition of the star-forming gas.

In Fig. 1 we compare the simulated ⟨ Ne^{2 +} ⟩ / ⟨ Ne⁺ ⟩ and ⟨ S^{3 +} ⟩ / ⟨ S^{2 +} ⟩ ratios of our model H ii regions with the corresponding ion ratios determined from the observations described by Rubin et al. (2007, 2008)¹³. The figure also shows that the metallicity-dependence of the ionization structure decreases for larger effective temperatures of the ionizing stars as can be expected from the ionizing fluxes shown in Table 4.

Emission line spectra of H  II regions at different metallicity.

Comparing observed emission line spectra to synthetic nebular models is the most important approach for obtaining information about the temperature, the density, and the ionization structure of H ii regions. To investigate the variations in the emission line ratios, we have computed the fluxes of the collisionally excited optical lines [N ii] 6584 Å + 6548 Å, [O ii] 3726 Å + 3729 Å, [O iii] 5007 Å + 4959 Å, and [S ii] 6716 Å + 6731 Å for each of our H ii region models. In Table 6 the results are shown relative to the corresponding Hβ emission.

In most cases the strongest line emission in the optical range is reached for a metallicity of 0.4 Z_⊙. This non-monotonic behavior is the result of two opposing metallicity-dependent effects. On the one hand, a higher metallicity increases the number density of ions that are potential line emitters. On the other hand, the lower temperature of the gas decreases the probability of a collisional excitation process, which finally leads to the emission of optical line radiation. The effect of the lower temperatures dominates for 2.0 Z_⊙ where in most cases the fluxes for this metallicity are the weakest among the considered metallicities¹⁴. This is explained by the equation for the rate coefficients for the collisional excitation, $\begin{array}{ccc}{\mathit{C}}_{\mathit{lu}}\mathrm{=}{\mathit{n}}_{\mathrm{e}}{\left(\frac{\mathrm{2}\mathit{\pi }}{\mathit{kT}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{ħ}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{e}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\frac{{\mathrm{\Omega }}_{\mathit{lu}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{g}}_{\mathrm{l}}}{\mathrm{e}}^{\mathrm{-}\mathit{h}{\mathit{\nu }}_{\mathit{lu}}\mathit{/}\mathit{kT}}& & \end{array}$(5)(Mihalas 1978; n_e is the number density of electrons, m_e the electron mass, Ω_lu(T) the velocity-averaged collision strength, a slowly varying function of temperature, hν_lu the transition energy, and g_l the statistical weight of the lower level), which leads to qualitatively different temperature dependences for collisionally excited lines in different wavelength ranges. The exponential term in Eq. (5) rises strongly with increasing temperature if the energy difference between the levels is large compared to kT, as is the case for lines in the visible and ultraviolet range. For the lines in the mid-infrared to far-infrared range, the energy difference between the upper and the lower level is significantly smaller than kT. Thus, the exponential term is close to unity in the entire temperature range found in H ii regions and the rate coefficients are roughly proportional to the inverse square root of the temperature. This implies that the cooling by infrared lines becomes more effective for decreasing temperatures, i.e., there is a positive feedback between the lower temperatures and the infrared line cooling. As a result, a small increase in the metallicity has a strong influence on the temperature structure if it causes the infrared transitions mentioned above to become the dominant cooling processes.

The influence of metallicity-dependent SEDs on the temperature and ionization structure of gas with fixed metallicity.

In addition to the models where stellar and gas metallicity are equal, we have computed models where the metallicity of the gas is fixed at 1.0 Z_⊙, but the metallicity of the stellar sources is varied in order to analyze the effects of stellar SEDs independently of the effects of the chemical composition of the gas. The simulations have been performed using the four different 40 000 K dwarf star models.

The sizes of the hydrogen Strömgren spheres show variations of approximately 5%, which are mainly caused by lower recombination rates in the gas around the lower-metallicity stars where the temperatures are higher than in the gas surrounding the model stars with higher metallicities¹⁵. The differences are more pronounced for the ionization structure of oxygen and sulfur, especially for the ratios O ii/O iii and S iii/S iv. For the model with a metallicity of 2.0 Z_⊙, O iii is the most abundant ion of oxygen within the inner ≈0.85 r_{H ii}, while this is the case for ≈0.97 r_{H ii} at a metallicity of 0.1 Z_⊙. For S iv the corresponding values are 0.28 r_{H ii} and 0.59 r_{H ii}, respectively.

The emission of the [O iii] 5007 Å + 4959 Å lines increases with decreasing stellar metallicities as a result of the higher temperatures and the larger fraction of O iii ions. At a metallicity of 0.1 Z_⊙ the emission of these lines is twice as strong as at 2.0 Z_⊙ (last part of Table 6). In contrast, the emission of lines associated with singly ionized atoms, like [N ii] 6584 Å + 6548 Å and [O ii] 3726 Å + 3729 Å, decreases for lower stellar metallicities because the number of singly ionized ions of these elements is reduced as the shells in which they are the dominant ionization stage become thinner.

3.2. Three-dimensional, time-dependent simulations of the ionization structures of the metals in H II regions

We will now abandon the assumption of a homogeneous gas and a spherically symmetric ionized bubble and discuss the effects of inhomogeneities of the gas density. This is motivated by the fact that H ii regions are obviously not homogeneous spheres, but are characterized by a more complex structure. The inhomogeneous density structure of H ii regions has been treated, for example, by Wood et al. (2004, 2005), who focused on the ability of photons to escape from a “porous” H ii region into the diffuse component of the ionized gas in a galaxy; Wood et al. (2013), who studied the consequences of inhomogeneous density structures for the determination of metallicities in H ii nebulae; and by Dale & Bonnell (2011), Walch et al. (2012), and Dale et al. (2013), who studied the interaction between stellar radiation fields and the density structure of ionized gas. Apart from the metallicity it is therefore the complex geometric structure and the embedded clumps that are of importance for the nebular models and a corresponding comparison of the calculated emission line strengths with the observed spectral features.

To include such effects in our computations we have developed a 3D radiative transfer code based on the approach by Weber et al. (2015), but considerably extended to account for the time-dependent ionization structure of the metals within the gas. In the description of our method we will first introduce our ray-tracing approach to 3D radiative transfer. Subsequently, we will present a solution of the time-dependent rate equations and treat the influence of the metal ions on the evolution of the temperature structure. Finally, we will describe the generation of synthetic narrowband images, which can link the theoretical models with observations.

This numerical approach will be applied to inhomogeneous, fractally structured gas with various metallicities and sources of ionization. The results of these computations results will be compared with the results obtained for a homogeneous gas distribution. In addition, we will account for the temporal evolution of H ii regions and study the effects of the distribution of the sources within clusters of hot stars on the emergent emission line flux of the surrounding ionized gas.

3.2.1. Three-dimensional radiative transfer based on ray tracing

The propagation of light along straight lines leads directly to the concept of a ray-by-ray solution, where the luminosity of each source is distributed among a set of rays originating from the source. The main aspects of the ray-by-ray solution are the isotropic distribution of rays around each source and the solution of the radiative transfer equation along each of these rays.

The sources are characterized by a specified spectral energy distribution F_ν and luminosity ${\mathit{L}}_{\mathit{\nu }}\mathrm{\left(}\mathit{l}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathit{l}}^{\mathrm{2}}\hspace{0.17em}{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathit{l}\mathrm{\right)}$, where R_l is the radius of the lth source. The frequency-dependent luminosity L_ν of each source is distributed evenly¹⁶ among N_rays rays, such that each ray is associated with a luminosity of $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{=}{\mathit{L}}_{\mathit{\nu }}\mathit{/}{\mathit{N}}_{\mathrm{rays}}$. The rays are then traced from the source(s) until the border of the simulated volume is reached.

The integrated form of the radiative transfer equation along a ray is $\begin{array}{ccc}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\hspace{0.17em}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{+}{\mathrm{\int }}_{{\mathit{s}}_{\mathrm{0}}\mathrm{\left(}\mathit{n}\mathrm{\right)}}^{\mathit{s}}{\mathit{\eta }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\hspace{0.17em}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\mathit{,}& & \end{array}$(6)where I_ν(s) is the intensity at the position s, s₀(n) corresponds to the starting point of each ray n, η_ν(s′) is the emissivity of the gas¹⁷ at the position s′, and ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{{\mathit{s}}_{\mathrm{0}}\mathrm{\left(}\mathit{n}\mathrm{\right)}}^{\mathit{s}}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}$ is the optical depth with respect to the source.

The simulated volume is discretized into a Cartesian grid of cells and each of the sources is located in the center of one of these cells. The energy deposited per time in the cells crossed by a ray n along the distance s_n(m) between the source and the entry point into a cell m is then calculated as $\mathrm{\Delta }\mathit{Ė}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\nu }$(7)(for details see Weber et al. 2015).

Temporal evolution of the ionization structure of metals.

Modeling the evolution of the ionization structures is not straightforward. The rate equations form a stiff set of differential equations as the timescales for the various ionization and recombination processes in the cell may differ by several orders of magnitude. Furthermore, the timescales differ between various cells depending on whether they are passed by an ionization front and how distant they are from the sources of ionization. However, the eigenvalue approach presented here provides a novel method of solving the rate equations in a stable and efficient way.

The solution ${x}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{x}}^{{\infty }}\mathrm{+}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$(8)of the inhomogeneous system of differential equations defined by Eq. (3) is composed of the equilibrium solution x^∞ of Eq. (4) and the solution of the corresponding homogeneous system ${E}\prime \mathrm{·}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{G}}^{{\prime }}\mathrm{·}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{0}\mathrm{,}$(9)

which is written in components as $\left(\begin{array}{c}\\ \frac{\mathrm{d}{\mathit{x}}_{\mathrm{hom}\mathit{,}\mathrm{1}}}{\mathrm{d}\mathit{t}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathrm{d}{\mathit{x}}_{\mathrm{hom}\mathit{,N}\mathrm{-}\mathrm{1}}}{\mathrm{d}\mathit{t}}\\ \mathrm{0}\end{array}\right)\mathrm{+}\left(\begin{array}{c}\\ {\mathit{g}}_{\mathrm{1}\mathit{,}\mathrm{1}}^{\mathrm{\prime }}& \mathit{...}& {\mathit{g}}_{\mathrm{1}\mathit{,N}\mathrm{-}\mathrm{1}}^{\mathrm{\prime }}& \mathrm{0}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{g}}_{\mathit{N}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{1}}^{\mathrm{\prime }}& \mathit{...}& {\mathit{g}}_{\mathit{N}\mathrm{-}\mathrm{1}\mathit{,N}\mathrm{-}\mathrm{1}}^{\mathrm{\prime }}& \mathrm{0}\\ \mathrm{1}& \mathit{...}& \mathrm{1}& \mathrm{1}\end{array}\right)\mathrm{·}\left(\begin{array}{c}\\ {\mathit{x}}_{\mathrm{hom}\mathit{,}\mathrm{1}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{x}}_{\mathrm{hom}\mathit{,N}\mathrm{-}\mathrm{1}}\\ {\mathit{x}}_{\mathrm{hom}\mathit{,N}}\end{array}\right)\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{0}\\ \mathrm{0}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathrm{0}\end{array}\right)\mathit{.}$(10)In Eq. (10) the condition of particle conservation is as an example inserted for k = N. The particle conservation implies that the sum of all components of a solution of the homogeneous system of differential equations is 0. For the solution x_hom(t) of Eq. (10), the kth component is thus computed as ${\mathit{x}}_{\mathrm{hom}\mathit{,k}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,i}\ne \mathit{k}}^{\mathit{N}}{\mathit{x}}_{\mathrm{hom}\mathit{,i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(11)The general structure of the solution for the components of the homogeneous equation different from k is ${\mathit{x}}_{\mathrm{hom}\mathit{,i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{-}\mathrm{1}}{\mathit{v}}_{\mathit{i,j}}{\mathrm{e}}^{\mathrm{-}{\mathit{\delta }}_{\mathit{j}}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}}{\mathit{\gamma }}_{\mathit{j}}\mathit{,}$(12)where v_i,j is the ith component of the jth eigenvector of the system of differential equations obtained from Eq. (10) by removing the kth line of x and the kth line and kth column of G^′; δ_j are the corresponding eigenvalues. The values of γ_j are defined by the occupation numbers at t₀, the beginning of the respective timestep.

Equations (11) and (12) are combined in matrix-vector notation using a matrix V composed of the eigenvectors v_j, but additionally contains in its kth row and jth column the negative of the scalar product 1·v_j (where the vector 1 contains the entry 1 in each of its components): ${{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{V}\mathrm{·}{\delta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{\gamma }\mathrm{=}\left(\begin{array}{c}\\ {{v}}_{\mathrm{1}}& {{v}}_{\mathrm{2}}& \mathit{...}\hspace{0.17em}& {{v}}_{\mathit{N}\mathrm{-}\mathrm{1}}\\ \mathrm{-}{1}\mathrm{·}{{v}}_{\mathrm{1}}& \mathrm{-}{1}\mathrm{·}{{v}}_{\mathrm{2}}& \mathit{...}\hspace{0.17em}& \mathrm{-}{1}\mathrm{·}{{v}}_{\mathit{N}\mathrm{-}\mathrm{1}}\end{array}\right)\mathrm{·}{\delta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{\gamma }$(13)The diagonal matrix δ(t) has entries δ_j,j = e^{− δ_j(t−t₀)}, and γ is a vector that contains in its jth row the coefficient γ_j (in Eq. (13) we set k = N again).

The vector γ is determined by inserting the occupation numbers for t = t₀ which are known from the previous timestep. In this case δ becomes the identity matrix and from Eq. (12) follows $\begin{array}{c}{˜}\\ {V}\end{array}\mathrm{·}{\gamma }\mathrm{=}{\begin{array}{c}{˜}\\ {x}\end{array}}_{\mathrm{hom}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\begin{array}{c}{˜}\\ {x}\end{array}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{-}\stackrel{}{{{{x}}^{{\infty }}}_{{˜}}}\mathit{,}$(14)where $\begin{array}{c}{{˜}}_{}\\ {V}\end{array}$ corresponds to the matrix V without the kth row (the kth row is also removed for ${\begin{array}{c}{{˜}}_{}\\ {x}\end{array}}_{\mathrm{hom}}$ and $\begin{array}{c}{{˜}}_{}\\ {{x}}^{{\infty }}\end{array}$). Equation (14) is solved for γ by ${\gamma }\mathrm{=}{\begin{array}{c}{˜}\\ {V}\end{array}}^{-1}\left({\begin{array}{c}{˜}\\ {x}\end{array}}_{\mathrm{hom}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{-}\stackrel{}{{{{x}}^{{\infty }}}_{{˜}}}\right)\mathit{.}$(15)With Eqs. (15) and (13), the values for x_hom(t)${{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{V}\mathrm{·}{\delta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{\begin{array}{c}{˜}\\ {V}\end{array}}^{-1}\left(\begin{array}{c}{˜}\\ {x}\end{array}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{-}\stackrel{}{{{{x}}^{{\infty }}}_{{˜}}}\right)$(16)are obtained. The result for x(t) is the initial value for the next timestep and after the recomputation of the rate coefficients the matrix G′ is updated. We currently consider the ionization stages I and II of hydrogen; I to III of helium; and I to IV of carbon, nitrogen, oxygen, neon, and sulfur, but in principle the method can be extended to any number of ionization stages and any set of elements required to describe the respective problem. The presented method has proven to be numerically very stable. Particle conservation is preserved for each element with a relative deviation of less than 10^-8.

In the next paragraphs we describe the computation of the rate coefficients for the recombination and ionization processes which via the rate matrix G′ determine the temporal evolution of the ionization fractions in x(t).

Recombination rates.

The relevant recombination processes in the interstellar gas are radiative and dielectronic recombination. In radiative recombination processes a free electron is captured by an ion and its kinetic and potential energy (relative to the bound state immediately after the recombination) is converted into the energy of the emitted photon. In dielectronic processes the energy of the captured electrons excites another electron of the ion, resulting in a doubly excited intermediate state. Dielectronic recombinations can only occur for discrete electron energies because of the additional bound-bound process. The recombination rate coefficient of an ion j in a grid cell m filled with gas at a temperature of T(m) is given by ℛ_j,j−1(m) = α_j(T(m))·n_e(m), where n_e(m) is the number density of electrons and α_j(T(m)) is the total recombination coefficient, composed of a radiative contribution α_r,j(T(m)) and a dielectronic contribution α_d,j(T(m)): ${\mathit{\alpha }}_{\mathit{j}}\mathrm{\left(}\mathit{T}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}{\mathit{\alpha }}_{\mathit{r,j}}\mathrm{\left(}\mathit{T}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{+}{\mathit{\alpha }}_{\mathit{d,j}}\mathrm{\left(}\mathit{T}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathit{.}$(17)For the temperature-dependent radiative recombination coefficients of metal ions of type j in cell m, we use the approximate formula ${\mathit{\alpha }}_{\mathit{r,j}}\left(\mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\mathrm{=}{\mathit{A}}_{\mathit{j}}{\mathit{T}}^{\mathrm{\prime }\mathrm{-}{\mathit{X}}_{\mathit{j}}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,}$(18)where T′(m) = T(m)/10 000 K is the temperature of the gas in units of 10 000 K. The fit parameters A_j and X_j are taken from Aldrovandi & Pequignot (1973; 1976), Shull & van Steenberg (1982), and Arnaud & Rothenflug (1985).

The dielectronic recombination rates for low temperatures, where kT is smaller than the energy needed to create the excited intermediate states, as is typically the case in H ii regions, are described by the fit formula (Nussbaumer & Storey 1983)¹⁸$\begin{array}{ccc}{\mathit{\alpha }}_{\mathit{dl,j}}\left(\mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)& \mathrm{=}& {\mathrm{10}}^{-12}\mathrm{\times }\left(\frac{{\mathit{a}}_{\mathit{j}}}{{\mathit{T}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\mathit{b}}_{\mathit{j}}\mathrm{+}{\mathit{c}}_{\mathit{j}}{\mathit{T}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{+}{\mathit{d}}_{\mathit{j}}{\mathit{T}}^{\mathrm{\prime }\mathrm{2}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\\ & & \end{array}$(19)Here a_j, b_j, c_j, d_j, and f_j are fit parameters which depend on the type of the ion.

An additional contribution to the dielectronic recombination rate that becomes relevant for higher temperatures where kT is in the order of the excitation energy has been described by (Burgess 1964) and is approximated by ${\mathit{\alpha }}_{\mathit{dh,j}}\left(\mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\mathrm{=}{\mathit{B}}_{\mathit{j}}{\left(\frac{\mathit{T}}{\mathrm{K}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\left(}\mathit{m}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{T}}_{\mathrm{0}\mathit{,j}}\mathit{/}\mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{C}}_{\mathit{j}}{\mathrm{e}}^{\mathrm{-}{\mathit{T}}_{\mathrm{1}\mathit{,j}}\mathit{/}\mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{)}\mathit{.}$(20)The values for the fit parameters B_j, C_j, T_0,j, and T_1,j were obtained from from the same sources as the radiative recombination coefficients. The total dielectronic recombination rate α_d,j is obtained by adding up the low-temperature and the high-temperature contributions (Storey 1983). The fit functions for the recombination rates already consider recombination to all levels. The recombination rates for hydrogen and helium are computed by interpolating the tables provided by Hummer (1994) and Hummer & Storey (1998).

Photoionization rates and computation of the mean intensity.

The photoionization rates R_lu(m) from a level l in one ionization stage to a level u in the next ionization stage are calculated as ${\mathit{R}}_{\mathit{lu}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathit{lu}}}^{\mathrm{\infty }}\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathit{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathit{h\nu }}{\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\nu ,}$(21)where J_ν(m) is the mean intensity, ν_lu the threshold frequency for the considered ionization process, and a_lu(ν) the frequency-dependent ionization cross section¹⁹.

To determine the average intensity J_ν needed to compute the radiative rates, we proceed as follows. From our discussion of the radiative transfer we already know the total energy absorbed (per unit time) by a cell, namely that given by Eq. (7). The number of photons absorbed per unit time in the cell in a particular transition is given by the same integral, with the integrand divided by the energy hν of a photon and weighted by the relative contribution of that transition to the total opacity, $\begin{array}{ccc}\sum _{{{}_{\mathrm{segments}}^{\mathrm{ray}}}^{}}\mathrm{\int }\frac{\stackrel{}{{{\mathit{L}}_{\mathit{\nu }}^{\mathrm{\left(}\mathrm{inc}\mathit{.}\mathrm{\right)}}}_{\mathrm{˜}}}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}}\right)}{\mathit{h\nu }}\frac{{\mathit{\chi }}_{\mathit{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{{\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{tot}}}\hspace{0.17em}\mathrm{d}\mathit{\nu }& \mathrm{=}& \mathit{V}{\mathit{n}}_{\mathit{l}}{\mathit{R}}_{\mathit{lu}}\\ & \mathrm{=}& \mathit{V}{\mathit{n}}_{\mathit{l}}\mathrm{\int }\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathit{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathit{h\nu }}{\mathit{J}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu ,}\end{array}$(24)where $\stackrel{}{{{\mathit{L}}_{\mathit{\nu }}^{\mathrm{\left(}\mathrm{inc}\mathit{.}\mathrm{\right)}}}_{\mathrm{˜}}}$ is again (for every ray passing through that cell) the luminosity of the ray incident on the cell, τ_ν is the total optical depth of the cell along that ray, and V is the volume of the cell. Since χ_lu(ν) is simply n_la_lu(ν), we see that the expression for J_ν that ensures consistency between radiative transfer and rate equations is ${\mathit{J}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi V}}\sum _{{{\mathrm{segments}}_{\mathrm{ray}}^{}}^{}}\frac{\stackrel{}{{{\mathit{L}}_{\mathit{\nu }}^{\mathrm{\left(}\mathrm{inc}\mathit{.}\mathrm{\right)}}}_{\mathrm{˜}}}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}}\right)}{{\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{tot}}}\mathit{,}$(25)which is (as it must be) independent of the actual transition considered in the discussion above, and can be used to compute the photoionization rates of all elements and ionization stages (for more details see Weber et al. 2015).

The temperature structure of evolving H ii regions.

The assumption that in an H ii region the heating and the cooling rates are in balance (as in the time-independent spherically symmetric approach presented in Sect. 3.1) and the temperature of the gas consequently remains constant is only valid for a steady-state H ii region. However, in evolving H ii regions the heating and cooling rates are in general not equal and the variation of the thermal energy content within the considered volume elements of the simulated gas has to be taken into account explicitly. Our simulations account for the photoionization of hydrogen and helium as heating processes. The photoionization heating rate (per volume) from a state l into the state u is ${\mathrm{\Gamma }}_{\mathit{lu}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathit{l}}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathit{lu}}}^{\mathrm{\infty }}\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathit{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{h}\mathrm{(}\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathit{lu}}\mathrm{)}}{\mathit{h\nu }}{\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\nu ,}$(26)where n_l is the occupation number of the lower state l and hν₀ the ionization energy. Again we assume the nebular approximation that all ionization processes take place from the ground state of an ion.

The considered cooling processes are the radiative recombination of hydrogen (Hummer 1994) and helium (Hummer 1994; Hummer & Storey 1998), as well as the radiative decay of collisionally excited ions²⁰. Furthermore, the simulations account for the cooling by free-free radiation (Osterbrock & Ferland 2006).

The heating and cooling rates are used to compute the change of the thermal energy content of the gas within a grid cell m. The temperature T(m) of the gas in a cell m is computed from the total content of thermal energy E_therm(m) in that cell by $\begin{array}{ccc}& & \frac{\mathrm{3}}{\mathrm{2}}{\mathit{N}}_{\mathrm{part}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{kT}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{V}{\mathit{n}}_{\mathrm{part}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{kT}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{therm}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,}\\ & & \mathit{T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}{\mathit{E}}_{\mathrm{therm}}\mathrm{\left(}\mathit{m}\mathrm{\right)}}{\mathrm{3}\mathit{V}{\mathit{n}}_{\mathrm{part}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{k}}\mathrm{·}\end{array}$(27)Here N_part(m) is the number of gas particles (electrons, atoms, and ions) within a grid cell, n_part(m) the number density of the particles, and V(m) the volume of the grid cell. After each timestep Δt, the thermal energy ${\mathit{E}}_{\mathrm{therm}}^{\mathrm{new}}\mathrm{\left(}\mathit{m}\mathrm{\right)}$ is recomputed as ${\mathit{E}}_{\mathrm{therm}}^{\mathrm{new}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{therm}}^{\mathrm{old}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{+}\left(\mathrm{\Gamma }\mathrm{\left(}{\mathit{J}}_{\mathit{\nu }}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{-}\mathrm{\Lambda }\left({\mathit{n}}_{\mathrm{e}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\right)\mathrm{·}\mathit{V}\mathrm{·}\mathrm{\Delta }\mathit{t,}$(28)where $\mathrm{\Gamma }\left({\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\mathrm{=}\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathrm{\Gamma }}_{\mathit{i}}\left({\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)$(29)and $\mathrm{\Lambda }\left({\mathit{n}}_{\mathrm{e}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{ff}}\left({\mathit{n}}_{\mathrm{e}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)\mathrm{+}\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathrm{\Lambda }}_{\mathit{i}}\left({\mathit{n}}_{\mathrm{e}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,T}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right)$(30)are the total heating and cooling rates per volume unit. The summations are carried out over all ionization stages i, where n_i are the number densities of the ions. The heating and cooling rate coefficients, Γ_i(J_ν(m)) and Λ_i(n_e(m),T(m)), depend in turn on the radiation field described by the mean intensity J_ν, on the electron density n_e, and the temperature T within a given grid cell. The free-free cooling rate is denoted by Λ_ff, Δt is the length of a timestep, and ${\mathit{E}}_{\mathrm{therm}}^{\mathrm{old}}\mathrm{\left(}\mathit{m}\mathrm{\right)}$ is the thermal energy content of the cell m before the timestep. The rates Γ_i and Λ_i already include the contributions of the different heating and cooling processes connected to an ionization stage i.

Computing synthetic images.

Images of gaseous nebulae show 2D projections of the 3D emission pattern of these objects. To link the results of our simulations with possible observations we create synthetic images that can be compared with narrow-band images taken in the wavelength range around diagnostically important emission lines.

The dilute gas found in H ii regions is almost transparent for the radiation of lines in the visible or infrared part of the spectrum if they are either emitted during the transition to a non-ground-level state (e.g., the lines of the Balmer series of hydrogen) or by a forbidden transition into a ground-level state (e.g., the forbidden line of O ii at 3729 Å). If the absorption and scattering terms are negligible, the intensity of the radiation at a given wavelength is given by ${\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathrm{\int }}_{{\mathit{s}}_{\mathit{c}}}^{{\mathit{s}}_{\mathit{d}}}{\mathit{\eta }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}$(31)Here η_ν is the emissivity of the medium and the integration over s along the line of sight is carried out between the point closest to the observer, s_c, and the point most distant from the observer, s_d. We consider the frequency-integrated intensities $\overline{\mathit{I}}$ and emissivities²¹$\overline{\mathit{\eta }}$ of the examined emission lines because the wavelength resolution in the 3D radiative transfer is not sufficient to resolve the line profiles. The number of emitted photons dN per solid angle dΩ, detector surface dA, and time dt therefore is

$\frac{{\mathrm{d}}^{\mathrm{3}}\mathit{N}}{\mathrm{d\Omega }\hspace{0.17em}\mathrm{d}\mathit{A}\hspace{0.17em}\mathrm{d}\mathit{t}}\mathrm{\approx }\frac{\overline{\mathit{I}}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\mathrm{=}{\mathrm{\int }}_{{\mathit{s}}_{\mathrm{c}}}^{{\mathit{s}}_{\mathrm{d}}}\frac{\overline{\mathit{\eta }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\hspace{0.17em}\mathrm{d}\mathit{s}$(32)for each of the considered lines, where ν₀ is the frequency of the center of the line. In our discretization scheme the integration along the line of sight is replaced by a summation of the emissivities multiplied by the lengths Δs(m) of the ray segments through the cells.

The small-angle approximation that the rays connecting the observer and all parts of the emitting regions are parallel can be used if the emission regions are small compared to their distance from the observer (as we assume in our simulations). In our case the summation is carried out along one of the coordinate axes (Fig. 2).

3.2.2. Applications of the three-dimensional approach for the simulation of H ii regions around hot stars

The presented 3D approach is now applied to examine various aspects of the interaction between hot stars and inhomogeneous H ii regions which have been outlined in the introduction of Sect. 3.

Comparison between homogeneous and inhomogeneous H  II regions.

First we use our procedure to simulate the interaction of the radiation field of hot stars and the inhomogeneous H ii regions surrounding these stars. The inhomogeneous H ii models are based on a fractal density distribution similar to the distributions that have been used to describe the interstellar gas by Elmegreen & Falgarone (1996) and by Wood et al. (2005). The results are then compared with the ionization structures and emission properties of homogeneous H ii regions with the same mean densities, metallicities, and sources of ionization.

In the inhomogeneous models, the mean number density of hydrogen is set to 9 cm^-3 for the fractal structures in the simulations. Additionally, the simulated gases contains a homogeneous fraction of 1 cm^-3. The total number densities n_H of hydrogen atoms vary between 1 cm^-3 (where the contribution of the fractal density field is zero) and 93 cm^-3. The density distribution of the gas is also characterized by the clumping factor ${\mathit{f}}_{\mathrm{cl}}\mathrm{=}\mathrm{⟨}{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathrm{⟨}{\mathit{n}}_{\mathrm{H}}{\mathrm{⟩}}^{\mathrm{2}}$ of ≈1.75 and the standard deviation of the hydrogen number density $\mathrm{(}\mathrm{⟨}{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}{\mathit{n}}_{\mathrm{H}}{\mathrm{⟩}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{8.7}\mathrm{c}{\mathrm{m}}^{-3}$. In Fig. 3 we show the number density of hydrogen for the mid-plane of the simulated volume, which in the presented simulations contains the source or sources of ionization. As the recombination and collisional excitation rates are proportional to the product of the number densities of the corresponding ions and the number density of electrons, they are approximately proportional to the square of the hydrogen number density in the ionized volumes. In the simulated inhomogeneous volume of gas the line emission per volume unit of the gas thus differs by more than 3 dex for a given temperature and number ratio of the ions. For this density structure we perform simulations of H ii regions around 40 000 K dwarf stars with metallicities of 0.1 Z_⊙ and 1.0 Z_⊙, where the relative metallicity in the gas of the H ii regions is the same as in the stellar atmospheres and the metallicity is assumed to be constant within the simulated volume.

The ionization fractions within the simulated volume, the total emission of some diagnostically important optical lines, and the maximum intensities for these lines in the synthetic images of the simulated H ii regions are compared in Table 7 . We find that the ratio between the number fractions of ionized hydrogen in the homogeneous and in the inhomogeneous models is approximately 1.6, which roughly agrees with the clumping factor. This behavior is a result of the proportionality of the recombination rates and the square of the gas density. By definition the mean square of the hydrogen density in the inhomogeneous case exceeds the mean square of the hydrogen density in the homogeneous case by the clumping factor. Consequently, a smaller number fraction of ionized hydrogen results for identical sources of ionization. Likewise, the clumping increases the neutral fraction for the other considered elements.

Figure 4 shows the projected images of the resulting photon emission rates $\overline{\mathit{I}}\mathit{/}\mathit{h}{\mathit{\nu }}_{\mathrm{0}}$ of the Hα line and of the collisionally excited lines of [N ii] 6584 Å + 6548 Å, [O ii] 3726 Å + 3729 Å, and [O iii] 5007 Å + 4959 Å. In the homogeneous case the density of the electrons and hydrogen ions varies only marginally over most of the ionized volume. The emission of the Hα line, which in the temperature range found in typical H ii regions is primarily a recombination line, consequently shows little structure apart from a decrease from the center to the border of the Strömgren sphere, which is a geometrical effect – the length of a path through a sphere decreases for larger impact parameters with respect to the source. In contrast, N ii and O ii primarily exist within the outer regions of the hydrogen Strömgren sphere. In the projected image of the line emission, this leads to a ring-like distribution of the projected emission line strengths. Figure 4 shows that for the inhomogeneous density structures the emission is more peaked around those lines of sight that intersect the density maxima (“clumps”).

How regions with different intensities in Fig. 4 contribute to the total flux is shown in the histograms in Fig. 5, using the Hα line emission as an example. The histogram covers the range for the values of $\overline{\mathit{I}}\mathit{/}\mathit{h}{\mathit{\nu }}_{\mathrm{0}}$ which appear in our synthetic images with 50 equidistant bins from 0 to 9 × 10⁸ photons s^-1 sr^-1 cm^-2. For each bin we add up the contribution to the total flux of all pixels whose intensities lie within the corresponding intensity intervals. In the homogeneous models the highest Hα intensities are reached for the lines of sight around the center of the ionized volume, while for the inhomogeneous models there is a broader distribution of the intensities. The histogram shows that for the inhomogeneous models the largest intensities exceed the largest intensities for the homogeneous case by a factor of about 4 to 5. Table 7 shows that similar factors also occur for the maximum intensities of the considered collisionally excited lines. Moreover, the contribution of high intensities for 1.0 Z_⊙ is larger than for 0.1 Z_⊙ in the homogeneous and in the inhomogeneous case. The reason is that the recombination coefficient is larger and the ionized volume smaller owing to the lower temperatures for the larger metallicity (Sect. 3.1.2), which leads to a reduced ionized volume. The total Hα emission is, however, nearly proportional to the number of ionization processes of hydrogen, which are almost the same in both models. As the same emission occurs in a smaller volume, larger intensities result. (The small differences between the hydrogen ionizing fluxes of the D-40 models at 0.1 Z_⊙ and 1.0 Z_⊙ are of minor importance.)

Temporal behavior of the ionization and temperature structure in model H  II regions.

We next investigate the time-dependence of the ionization and the temperature structure of the gas during the expansion of H ii regions into previously neutral interstellar gas with a fractal density structure. These simulations were performed for the 40 000 K dwarf stars with metallicities of 0.1 Z_⊙ and 1.0 Z_⊙ (the metallicities within the gas of the H ii regions match the metallicities in the stellar atmospheres). In these simulations the total mass density of the gas in each cell is assumed to remain constant²².

The results are presented in Figs. 6 and 7. In both cases the temperature increases quickly during the ionization process of previously neutral gas crossed by the ionization front. The heating rate is proportional to the number of ionization processes for a given ionizing spectrum. Therefore, the heating rate is increased until the gas is ionized to a large extent, such that it becomes optically thin for the ionizing radiation. Unlike the heating rates, the relevant cooling rates are approximately proportional to the square of the electron density, i.e., they continuously increase during the ionization process. Consequently, the cooling rates reach their maximum later than the heating rates and the rise in the temperature of the gas during the ionization process is followed by a decrease until the equilibrium temperature is reached. In the inner parts of the ionized region where the ionization timescale and consequently the heating timescale are considerably shorter than the recombination and cooling timescale, the temperature maximum reached depends on the ionizing SED rather than on the metallicity-dependent efficiency of the cooling within the gas. This overshooting effect of the temperature is more pronounced for the simulation of the H ii region with solar abundance presented in Fig. 7 than for the metal-poor region in Fig. 6 because of the lower equilibrium temperature for the larger metallicity. When the ionized volume is close to its equilibrium value and the outer parts of the equilibrium Strömgren volume are ionized, the ionization process is slower than for the inner regions of the H ii region because of the more dilute radiation field and the timescale for the heating process becomes comparable to the cooling timescale. As a result, the value of the temperature peak relative to the equilibrium value is lower than in the inner regions of the H ii region.

The overshooting of the temperatures affects the emission spectrum of the expanding nebula such that the intensities of collisionally excited lines like the optical [O ii] and [O iii] lines decrease when the gas cools down. The effect is more pronounced for the gas with 1.0 Z_⊙ (Fig. 7) than for the gas with 0.1 Z_⊙ (Fig. 6).

A comparison of the temperature and the density of the ionized gas in the first two rows of Fig. 7 shows that the cooling process is faster for denser gas. For given temperature and ionization fractions, the thermal energy content increases linearly with the density of the gas. On the contrary, the cooling rates are approximately proportional to the square of the density and the cooling timescale is therefore roughly proportional to the inverse of the gas density.

Simulations of interstellar gas irradiated by clusters of hot stars.

Finally, we want to see how these results change if we replace the single star in the previous models by a small cluster of stars and study the effects of a different spatial distribution of the stellar sources. The aim of this comparison is to get a rough idea of how the spatial distribution influences the total nebular emission, and we therefore simply use single stars as ionizing sources, and do not attempt to model a cluster SED via population synthesis.

We first assume that these stars are in close proximity to each other (which we realize by placing them in the same grid cell) and are formed instantaneously at the same time. For the composition of the cluster we assume three 30 000 K stars and two 35 000 K stars (models D-30 and D-35), which together have approximately the same quantity (89%) of the hydrogen-ionizing photon flux as our 40 000 K model dwarf star at solar metallicity, although the total luminosity of the cluster exceeds that of the 40 000 K model star by a factor of 2.8. We assume solar metallicity for both stars and gas.

The softer SED of the ionizing radiation of the cluster results in a considerably smaller fraction of N iii, O iii, Ne iii, and S iv as a comparison of Tables 7 and 8 shows, although the similar hydrogen-ionizing fluxes (and the identical metallicity which results in comparable temperatures and consequently recombination rates) lead to almost identical hydrogen Strömgren volumina. In the temporal evolution there is an overshooting of the temperature of the gas passed by the ionization fronts like in our simulations containing single stars. The evolution of the ionized gas around this cluster is shown in Fig. 8.

The results for the ionization fractions and the line emission do not change significantly (i.e., by less than 15%) if the distances between the stars are increased even up to the typical sizes of the Strömgren spheres, such that no single H ii region is formed, but the ionized volumes around the stars merge only partially (see Fig. 9).

As is expected from the different SEDs, the ionization structure of oxygen around the different types of stellar sources varies significantly among the single H ii regions. In the regions around the 35 000 K stars O ii is the most abundant ionization stage of oxygen up to approximately half the distance between the sources and the border of the hydrogen Strömgren volume. Around the 30 000 K stars, the volumes in which O iii is more abundant than O ii are too small to be spatially resolved within our simulation. This is in agreement with the ionizing fluxes presented in Table 4 and the results of the spherically symmetric H ii region models in Sect. 3.1.2.

Spectroscopic signature of H ii regions around very massive stars.

The direct detection of very massive stars (Sect. 2.2) is a problematic procedure for two reasons. First, they are rare objects such that even the closest of these objects are likely to exist only on extragalactic distance scales. Second, their occurrence is correlated with massive stellar clusters where they might form by stellar mergers (cf. Pauldrach et al. 2012), such that observations of distant objects are likely to suffer from crowding effects. In this section we therefore examine the possibility of identifying VMS by means of the emission spectrum of the surrounding H ii regions. We simulate H ii regions around the 45 000 K and 65 000 K stars with M_∗ = 3000 M_⊙, which represent the stars with the highest masses in our VMS grid, and the 40 000 K and 50 000 K stars with M_∗ = 150 M_⊙, which are close to the currently known upper mass limit of O stars. The luminosity has been chosen to correspond to a stellar content of 3000 M_⊙, i.e., we simulate single 3000 M_⊙ stars and clusters of twenty 150 M_⊙ stars. The size of the simulation volume is (400 pc)³. In the simulations all stars are located in the center this volume. In Fig. A.11 we show the emission line intensities for [O ii] 3726 Å + 3729 Å, [O iii] 5007 Å + 4959 Å, [O iii] 88.3 μm, and [O iv] 25.9 μm.

The volume containing singly ionized helium is considerably smaller than the hydrogen Strömgren volume around the 3000 M_⊙ model star with T_eff = 45 000 K and solar metallicity because of the strong absorption within the stellar wind of this object (see Sect. 2.2). Outside the He i/He ii Strömgren sphere no ionization stage of oxygen above O ii is reached. Consequently, the [O ii] 3726 Å + 3729 Å line emission exceeds the emission of the [O iii] 5007 Å + 4959 Å lines, which is not the case for the other dwarf or supergiant models with T ≥ 40 000 K. Both O ii and O iii exist in the part of the ionized volume where the predominant ionization stage of helium is He ii, but the fraction of O iii increases for higher effective temperatures of the ionizing stars. The O iv emission line at 25.9 μm and the He ii recombination line at λ = 4696 ^˚A are emitted by regions where helium is ionized to He iii, i.e., around objects with large effective temperatures such as the 65 000 K model VMS. The latter lines have been detected in starburst galaxies (Lutz et al. 1998). Here hot (T> 50 000 K) very massive stars created by mergers in dense clusters may provide a complementary explanation to shocks (Lutz et al. 1998) or the photoionization by Wolf-Rayet stars (Schaerer & Stasińska 1999). The O iv line is, however, strongly metallicity-dependent as for lower metallicity there is a smaller number of potential emitters and the collisional excitation rate coefficient is reduced for larger temperatures of the metal-poor gas (Sect. 3.1.2). Therefore, the emission of the [O iv] by gas around very hot and very massive stars is more likely to be observed in an already chemically enriched environment.

4. Summary and conclusions

In this work we have presented state-of-the-art SEDs of hot, massive stars as a contribution to the study of the line emission from the ionized gas around stars and stellar clusters. Such H ii regions are essential tracers for the state and evolution of galaxies because they are closely connected to star formation processes and their emission line spectra provide information about the chemical composition, the structure, and the dynamics of the interstellar gas in galaxies, but their spectra depend to a large extent on the SEDs of the ionizing stellar sources.

Our stellar model grid comprises O-type dwarf and supergiant stars with temperatures from 30 000 K to 55 000 K and metallicities between 0.1 Z_⊙ and 2.0 Z_⊙, representing the metallicity range of star-forming regions in the present-day universe. A comparison of the ionizing fluxes for different metallicities shows that even for the same effective temperature there are considerable differences (up to several orders of magnitude) in the emission rates of photons in the energy ranges required to produce key ionic species (such as O iii or Ne iii, which have important diagnostic lines). In most cases the fluxes in these energy ranges decrease with increasing stellar metallicities because of the enhanced line blocking in the atmospheres of the stars.

Additionally, we have computed the ionizing fluxes of very massive stars that possibly form through collisional mergers in dense massive clusters and may considerably exceed the mass limit of O-type stars forming directly from molecular gas clouds. For these objects we have considered masses from 150 M_⊙ to 3000 M_⊙ and temperatures from 40 000 K to 65 000 K.

We have used the SEDs in spherically symmetric photoionization models to gauge the effect of stellar and gas metallicity on the ionization structure of the H ii regions. As expected from a direct assessment of the SEDs, for a given stellar effective temperature the integrated number fractions of the different ionization stages of the metals in the gas strongly depend on the metallicity-dependent SEDs for those species whose ionization energies lie above the ionization edge of He i. Even at the same metallicity and temperature, the results differ significantly between dwarfs and supergiants owing to the influence of the much stronger stellar winds of supergiants on the SEDs.

To judge the effects of an inhomogeneous distribution of the gas on the nebular line emission we have used our 3D radiative transfer code, augmented to account for the metals C, N, O, Ne, and S, to compute integrated line-of-sight intensities and synthetic narrowband images. Although an inhomogeneous density structure leads to higher peaks in the distribution of the observable line intensities for a given mean density of the gas, we find no relative changes above 20% in the total emission per line.

With regard to the temporal evolution of the ionization and temperature structure of the H ii regions we find an overshooting of the gas temperatures significantly above the equilibrium value during the transition of the gas from initially neutral to ionized. However, because of the short timescales, and because in reality the ionizing sources themselves evolve, this is unlikely to be of practical significance for the analysis of H ii regions with current models.

Finally, we have studied the effects of clustering of stars in two limiting cases, a dense cluster where the distances between the stars become irrelevant compared to the extent of the combined H ii region, and a loose association where the distances between the stars are of the same order as the Strömgren radii, so that the individual H ii regions only partly overlap. We find that the differences between the two models in the number fractions of the ionization stages and the total emission rates per line remain moderate – the total emitted fluxes of the considered emission lines differ by less than 10%.

Although spherically symmetric modeling of H ii regions is still an important tool for determining the global properties of the ionized gas and for performing systematic examination of large sets of different input parameters, 3D radiative transfer models accounting for inhomogeneous density structures, several sources of ionization, and a time-dependent evolution of the ionized medium are able to provide additional information for studying H ii regions in more detail. A comparison of the results from these models with narrow-band or spectrally resolved observations allows detailed information to be obtained about the internal structure of H ii regions and about its relation to the ionizing stars.

Online material

Appendix A: Figures and model data

---

Acknowledgments

We thank an anonymous referee for helpful comments which improved the paper. We also wish to thank K. Butler for helpful discussions concerning the atomic data used for our computations and B. Ercolano for useful comments. We would like to thank N. Gnedin for providing the IFrIT program used for some visualizations of our results. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under grants PA 477/18-1 and PA 477/19-1.

References

All Tables

All Figures