© ESO 2019

1. Introduction

Aside from cool molecular or atomic clouds with temperatures of approximately ten to a few hundred K, the warm neutral gas, and the very hot and thin fraction of the interstellar gas with temperatures on the order of one million K (which is heated, for instance, by supernova shocks), the interstellar medium contains a component of ionized material with an intermediate temperature of T ≈ 10 000 K. H II regions around hot stars (dwarf or supergiant stars of spectral class O or B, or central stars of planetary nebulae) are one part of this component. But a significant amount of the mass and the line emission by intermediate-temperature gas is additionally contributed by dilute gas with a mean density of ≲ 0.1 hydrogen ions per cubic centimeter, which can be found – in contrast to H II regions immediately connected to star formation – at heights of up to a few kpc above the galactic plane. This dilute gas is called warm interstellar medium (WIM, often used when referring to gas in the Milky Way) or diffuse ionized gas (DIG, usually referring to the diffuse gas in other galaxies)¹.

The nature of the energy sources that keep the DIG ionized is still not well understood. Photoionization by the radiation of the embedded sources is the dominant type of energy input for classical H II regions around hot stars; observations suggest that it also plays an important role for the ionization of the DIG.

For instance, Haffner et al. (1999) and Zurita et al. (2002) found a correlation between spiral arms and the presence of DIG. The ratio of the hydrogen recombination radiation photons originating from H II regions and the DIG has been determined observationally by Zurita et al. (2000), who found, based on an analysis of the observed emission measures in six nearby galaxies², that the total emission from the DIG is in the range of approximately 30–70% of the emission originating from the denser H II regions. Assuming ionization of the DIG predominantly by photons leaking from H II regions, this would require about half of the ionizing photons to escape from these H II regions.

Some of the spectral features found in the DIG are, however, considerably different from the features of dense H II regions ionized by embedded hot stars. Observations of the warm interstellar medium in our galaxy by Reynolds (1985) showed that the [O III]∕Hα line ratio is considerably smaller in the Galactic WIM than in typical H II regions. Reynolds & Tufte (1995) found only weak emission from the helium recombination line at 5876 Å, which indicates that the helium within the WIM of the Milky Way is mostly neutral. Measurements of the temperature of the gas, such as the comparison of the nebular and auroral lines of N II performed by Reynolds et al. (2001), show temperatures clearly above the value characteristic for the gas in H II regions, which also leads to enhanced emission in the collisionally excited optical lines of singly ionized nitrogen and sulfur. An extensive study of the properties of the WIM in the solar environment of the Milky Way, based on observations using the Wisconsin H α Mapper spectrograph, has been presented by Madsen et al. (2006).

Observations by Greenawalt et al. (1997) indicate a softer ionizing spectrum than in typical H II regions for the DIG of M 31 as well, but this difference appears to be less pronounced than in the Milky Way because the spectra show [O III] emission and the presence of the helium recombination line at 5876 Å, which indicates an appreciable fraction of ionized helium. The dominance of only singly ionized metals in the DIG of the solar environment in the Milky Way cannot, however, be considered a general property of the entire DIG in spiral galaxies. Observations of edge-on galaxies (Otte et al. 2001; Collins & Rand 2001) show that in a part of these objects the [O III]∕Hα ratio rises at large heights above the galactic planes, while large [N II]∕Hα and [S II]∕Hα ratios are also observed. Evidence for a hardening of the ionizing spectra at high galactic altitudes has also been provided by Spitzer Space Telescope observations of three edge-on galaxies in the mid-infrared by Rand et al. (2011), who found a rising Ne III/Ne II ratio for larger distances from the galactic disks. The differences between the observations indicate that the energy sources and the temperature and ionization structure of the DIG may differ considerably from galaxy to galaxy or among different parts of the same galaxy.

Two main explanations for the differences between the spectra of classical H II regions and the DIG have been put forward, which are likely to be complementary: (a) filtering effects modifying the spectral energy distribution (SED) of the ionizing radiation, and (b) the presence of additional sources of ionization and heating. Filtering effects have been studied, for example, by Hoopes & Walterbos (2003) and Wood & Mathis (2004) and are partly able to explain the larger temperatures of the DIG (compared to H II regions) as a result of the radiation hardening effect. We discuss filtering effects in more detail in Sect. 2.3.

An additional contribution of hot low-mass evolved stars (central stars of planetary nebulae or white dwarfs) has been discussed, for example, by Sokolowski & Bland-Hawthorn (1991) and Flores-Fajardo et al. (2011) to explain the fact that in some galaxies the [O III]/[O II] ratio increases for larger distances from the disk. Such objects can occur at larger distances from the plane of the star-forming galactic disk owing to their older ages compared to hot massive main-sequence or supergiant stars.

A possible source for additional heating of the DIG is the dissipation of a part of the (macroscopic) mechanical energy of the gas. There have been several approaches to explain the underlying mechanisms. One of these is shocks that can, for example, be produced by supernova explosions or strong stellar winds from massive OB associations (Collins & Rand 2001). Another possibility for providing the DIG with energy is via magnetic reconnection processes as described by Hoffmann et al. (2012), who analyzed the contributions of photoionization and magnetic reconnection by computing synthetic spectra for purely photoionized or purely thermally (by magnetic reconnection) excited gas and for gas in which both mechanisms are present. These authors found that while the total energy input provided by magnetic reconnection alone may be sufficient to keep the DIG ionized at large heights above the galactic plane, the observed spectra can be explained best by a mixture of both photoionization and thermal excitation.

In this work, however, we concentrate on the ability of hot massive stars close to the midplane of the galactic disk to ionize the DIG. Not much is known about the ionizing fluxes of these objects through direct observation. For this reason we employ realistic SEDs from state-of-the-art, non-LTE³ atmosphericmodels that include the effects of stellar winds and line blocking and blanketing. The SEDs from our group have had a consistently good track record as ionizing sources in models of H II regions (Sellmaier et al. 1996; Giveon et al. 2002; Rubin et al. 2007, 2008, 2016). The use of such SEDs becomes even more important as we model gas and stars at different metallicities, where even fewer detailed direct and indirect observational constraints can be obtained from our nearer environment.

In this paper we systematically explore a plausible subset of the (rather large) parameter space in metallicity, stellar effective temperature, and filtering through surrounding H II regions, to show how observable characteristics of the DIG vary quantitatively with these parameters. As a parameter study, none of the models in our current paper is intended as a detailed description of the DIG in any particular galaxy, and neither is it expected that each particular parameter combination will be realized by nature somewhere in our immediate surroundings. We expect, however, that the large set of theoretical predictions that we present here will serve as a future guide in exploring more detailed quantitative models.

In Sect. 2 we study the properties of the DIG as a function of the temperature- and metallicity-dependent SEDs of the ionizing sources (based on the model stars presented by Weber et al. 2015) and the escape fractions from the H II regions surrounding the ionizing stars, using the spherically symmetric method described by Hoffmann et al. (2012) and Weber et al. (2015). Together, our study comprises emission line ratios for a grid of 60 spherically symmetric DIG models.

A spherically symmetric approach is suitable for systematically examining the differential variations of the emission line spectrum in dependence of a certain input parameter, such as the SED of the ionizing source. There are, however, two drawbacks regarding simulations of the diffuse ionized gas. First, star-forming galaxies usually contain numerous H II regions whose stars may act as sources of ionization for more distant gas. Second, the density distribution of the absorbing gas is usually not spherically symmetric, and thus the attenuation of the radiation field is anisotropic. To account for these effects of the three-dimensional (3D) structure of the interstellar gas, we present in Sect. 3 a set of 3D simulations whose geometry corresponds to a galactic region above a star-formingspiral arm. We show the first synthetic emission line ratio maps of the DIG around a spiral arm with an SED based on a realistic stellar initial mass function (IMF) and compare the results with those obtained when assuming a softer ionizing spectrum corresponding to the SED of a model star with an effective temperature of 35 000 K. We also show how, even for the same mean density and average filling factor, the ionization structure of the DIG is affected by its porosity, which is a function of the size and spatial distribution of the overdense “clumps”.

2. Spherically symmetric simulations of dilute ionized gas surrounding leaky H II regions

In this section we describe a grid of H II region and DIG models in which we vary the effective temperatures T_eff and metallicities Z of the ionizing stars, as well the escape fraction of ionizing photons from the H II regions surrounding the stars into the DIG⁴. Our approach corresponds to that of Hoopes & Walterbos (2003), but our study is based on stellar SEDs from non-LTE atmospheric models taking into account the effects of stellar winds and line blocking in the extreme ultraviolet (EUV), namely the model grid of Weber et al. (2015) at 0.1, 0.4, and 1.0 solar metallicities. We first describe our numerical approach to accurately model the atmospheres of O-type stars that act as sources of ionization. Such detailed modeling is crucial for the determination of the radiation field in the ionizing spectral range above 1 Ryd, which in turn determines the temperature and ionization structure of the gas affected by the stellar radiation. Then we outline the method used to describe the processes in the ionized nebulae. After an explanation of the simulation setup used, we discuss the effects of the filtering on the ionizing radiation field by the gas in the stellar environment and its effect on the emission line spectrum of both the H II regions in the stellar environment and the DIG. Finally we compare our results with observations and discuss the implications of the comparison.

2.1. 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). We summarize the salient points in this section. 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 the following key components: (a) a solution of the hydrodynamics describing velocity and density structures of the outflow based on radiative acceleration by line, continuum, and Thomson absorption and scattering; this hydrodynamic description is an essential aspect of the model because the expansion of the atmosphere alters the emergent flux considerably compared to a hydrostatic atmosphere;(b) the computation 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) a calculation of the radiation field from a detailed radiative transfer solution taking into account the Doppler-shifted line opacities and emissivities along with the continuum radiative transfer⁶; and (d) a computation of the temperature from the requirement of radiative (absorption/emission) and thermal (heating/cooling) balance. An accelerated Lambda iteration procedure⁷ (ALImI; explained in detail by Pauldrach et al. 2014) 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, these effects have not been included in the stellar models used to describe the ionizing sources in this paper, since they affect only high ionization stages such as He III, which are not relevant for the analysis of emission line spectra of H II regions around O-type stars (see Kaschinski & Ercolano 2013).

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, 2016). A more stringent test can be provided by a comparison of the synthetic and observed UV spectra of individual massive stars. This comparison involves hundreds of spectral signatures of multiple ionization stages of various elements and with different ionization thresholds, 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. Numerical modeling of the H II regions

For the computation of spherically symmetric models of gaseous nebulae we have adapted the stellar atmosphere code discussed above (Pauldrach et al. 2001, 2012) to the dilute radiation fields and low gas densities in H II regions (Hoffmann et al. 2012, Weber et al. 2015). The results of this approach describe steady-state H II regions in which ionization and recombination, as well as heating and cooling, are in equilibrium at each radius point.

In the stationary state, the ionization fractions of an element are given by the solution n^∞ of the equation $$G\cdot n^{\infty }=0,$$(1)

where G is a matrix formed from the ionization and recombination rate coefficients⁸, and the components of the vector n^∞ are the number densities of successive ionization stages of the elements considered. Because the rate coefficients in general depend on the radiation field and the temperature (themselves in turn dependent on the ionization ratios), the system is solved iteratively together with the radiation field and temperature until convergence is reached.

In our spherically symmetric models we compute the radiative transfer along a series of rays with different impact parameters to the central source, yielding the radiative intensity I_ν at every point along each ray. At any given radius shell r, each ray intersecting that shell corresponds to a different angle ϑ with respect to the outward normal, and the mean intensity J_ν(r) that enters into the radiative rates (photoionization and photoheating) is then obtained as $J_{\nu }(r)=\frac{1}{2}{{\displaystyle \int }}_{-1}^1{I}_{\nu }(r,\mu )\text{d}\mu$, with μ = cosϑ. The symmetry allows taking into account the diffuse radiation field (ionizing radiation created in the gas itself via direct recombination to the ground level) without incurring additional computational cost, and thus does not necessitate the use of simplifying approximations such as “on the spot” or “outward only”.

The emission spectrum of the gas depends on both its ionization structure and its temperature, since the recombination and collisional excitation rates are temperature dependent. Interpretation of the spectra of ionized gas thus also requires considering the microphysical processes that heat and cool the gas and thereby determine its temperature structure. In H II regions,the main heating source is the ionization of H I and He I, which increases the kinetic energy of the gas by $\Delta E=h\left(\nu -{\nu }_0\right)$ per ionization process; ν is the frequency of the ionizing photon and hν₀ the ionization energy. Cooling of the gas happens via recombination processes, radiative de-excitation of collisionally excited states, and free-free radiation (bremsstrahlung).

In H II regions, a large part of the energy of the absorbed ionizing photons is transformed into line emission in the UV, visible, and infrared wavelength range. These lines are either recombination lines like the Balmer lines of hydrogen, or collisionally excited linessuch as the “nebulium” lines of O III at 4959 and 5007 Å. To compute the emission in these lines, the bound–bound transition rates among the different excitation states within the ionization stage must be calculated. Our numerical approach for computing these rates is described in detail by Hoffmann et al. (2012) and Weber et al. (2015).

2.3. Properties of the filtered stellar radiation

In our first set of simulations we study the emission spectrum of diffuse gas ionized by hot stars surrounded by an H II region of denser gas that attenuates and filters the stellar radiation. For this we use spherically symmetric models and assume a constant density in each of the two gas regions considered, namely n_H = 10 cm⁻³ for the H II region and n_H = 0.1 cm⁻³ for the DIG. Our model grid comprises stellar effective temperatures of 30, 35, 40, and 50 kK and uses the SEDs of the dwarf (i.e., main sequence) star models from Weber et al. (2015) with metallicities of 0.1, 0.4, and 1.0 Z_⊙⁹.

Our procedure involves a two-step process. In the first step, we compute models for H II regions, which we assume to be ionized by a single model star each, and iteratively adjust the radii of these (matter-bounded) H II regions until a given fraction (70, 50, 30, 10, and 5%) of the total hydrogen-ionizing photon flux of the stars is able to escape. In the second step, we use the filtered stellar radiation from these models to illuminate a homogeneous region of much thinner gas representing the DIG, which we take to be radiation-bounded in all cases. For better comparability, we scale the ionizing fluxes escaping from the H II regions such that the number of hydrogen-ionizing photons is 10⁵⁰   s^-1 in all cases; i.e., we compare sources with equal photon-number rates but different spectral shapes. We furthermore assume the metallicity of the DIG to be the same as that of the inner H II regions.

In our models we do not currently consider the effects of dust absorption or scattering. While dust is present in H II regions, Mathis & Wood (2005) find the absorption by silicate-based dust to be weak and – in the frequency range of ionizing radiation – almost wavelength independent. (Thus, it may reduce the total amount of escaping ionizing radiation, but has no large effect on its spectral distribution.) They also call attention to the fact that while carbon-based dust (composed of graphite particles or polycyclic hydrocarbons) has an absorption feature that peaks at 17 eV and could therefore alter the SED of the ionizing radiation field, observations showing that there is almost no C-depletion from the gaseous phase indicate that the carbon-based dust in H II regions has been mostly vaporized. At present we also do not account for the different depletion of individual types of metals from the gaseous phase, as discussed, for instance, by Shields & Kennicutt (1995).

With regard to the amount of light in the optical wavelength range (including the diagnostic lines) emitted by H II regions that is scattered by dust in the DIG, the evidence is still inconclusive. While Seon & Witt (2012) find a large fraction of around 50%, Ascasibar et al. (2016) find a much smaller contribution of only about 1–2%. Thus, while we do not consider such scattering in our models and our computed line ratios therefore represent the “true” line emission by the DIG, the possibility must be kept in mind that the observed line ratios from the DIG may be contaminated to a yet undetermined degree by light not directly originating in the DIG.

In Fig. 1 we plot the SEDs F_⋆(ν) of the ionizing sources and the fluxes F(ν) escaping from the leaky (matter-bounded) H II regions around these sources, as a function of effective temperature, chemical composition, and the total fraction of escaping hydrogen-ionizing photons. Figure 2 shows the spectral shape of the escape fraction f_esc (ν) = F(ν)∕F_⋆(ν), i.e., the “filter functions” of the H II regions.

In the energy range between the ionization edges of hydrogen and (neutral) helium, the escape fraction increases for higher photon energies because the absorption cross section of hydrogen decreases roughly proportional to ν⁻³. However, in matter-bounded H II regions diffuse radiation from recombination processes within the gas is also emitted from the H II region into the surrounding gas. The total emission from the H II region results from the sum of the unabsorbed photons from the star and the diffuse radiation that is able to escape from the gas and is responsible for the increased emission in the energy range directly above the ionization threshold of hydrogen.

At the He I ionization edge (hν₀ = 1.81 Ryd) helium becomes an important additional absorber, and the escape fraction again drops significantly. In the case of the 30 000 K stars, practically no ionizing radiation above the ionization edge of neutral helium can escape from the H II region. This is because the total He-ionizing photon emission rate of these stars is very small (approximately 2.5 dex smaller thanthe hydrogen-ionizing flux; see Weber et al. 2015), and thus the helium within the H II region is sufficient to absorb almost the entire He-ionizing radiation. This latter result is in agreement with earlier work by Hoopes & Walterbos (2003) using LTE stellar atmospheres.

For the 35 000 K stars and escape fractions of 5 and 10% of hydrogen-ionizing photons, the helium-ionizing photons are also nearly completely absorbed in the H II regions, whereas for the larger escape fractions a part of the He I-ionizing radiation can escape into the surrounding DIG. For higher stellar temperatures, the effect of the additional absorption by neutral helium becomes less important because the emission rate of helium-ionizing photons by the star grows in relation to that of hydrogen-ionizing photons. Thus, a smaller fraction of the emitted helium-ionizing photons is consumed to keep the helium in the H II regions ionized.

In the energy range above the ionization edge of He I, the escape fraction rises again as a result of ionization cross sections of both hydrogen and helium decreasing for larger frequencies. In this energy range there is also a noticeable contributionof metal ions to the absorption of the stellar radiation. The most important absorbing metal ions are O II (hν₀ = 2.59 Ryd), Ne II (3.25 Ryd), N III (3.49 Ryd), and C III (3.52 Ryd). The effect is more pronounced for higher metallicities, and for T_eff ≥ 40 000 K only the latter two ions are relevant as absorbers in our model H II regions, while absorption by the former becomes negligible because of their decreasing abundance.

None of our nebular models exhibits a significant amount of escaping He II-ionizing radiation (hν ≥ 4 Ryd). This is a consequence of the even lower stellar photon emission rates in this energy range, which are several orders of magnitude lower than the corresponding H and He I-ionizing fluxes. Even in the case of the 50 kK model star with a metallicity of 0.1 Z_⊙, which is the strongest emitter of He II-ionizing radiation of the stars considered here, the emission rate of He I-ionizing photons exceeds that of He II-ionizing photons by a factor of approximately 1500.

Another factor influencing the frequency-dependence of the escape fraction in real H II regions is the density structure of the filtering gas. While radiation hardening effects such as those described above are pronounced in homogeneous H II regions, the attenuation would become independent of wavelength (though strongly varying with angle) if all of the radiation were completely blocked in some directions by discrete, completely opaque clumps, and no absorbing matter at all intervened in other directions (picket-fence case). This latter scenario (wavelength-independent absorption by opaque clumps) and the former (filtering of the radiation by homogeneous H II regions) represent hypothetical limiting cases. Actual H II regions (and the DIG) show a more complex density structure that has characteristics of both cases¹⁰. We examinethe influence of different density distributions in the context of our 3D simulations in Sect. 3 below.

We approximate the picket-fence case by increasing the radii of the H II regions until they become radiation-bounded, while also using the unfiltered SEDs of the model stars as the ionizing fluxes for the DIG (again scaled to 10⁵⁰ ionizing photons per second).

Fig. 1 Computed stellar SEDs (black lines) modified by homogeneous matter-bounded H II regions (nH = 10 cm-3) for different escape fractions (colored lines as labeled). For the sources we have chosen the main-sequence spectra from the hot star model grid by Weber et al. (2015). The chemical composition of the H II regions is assumed to be the same as that of the corresponding stellar atmospheres. Shown are the Eddington fluxes normalized to the stellar radius.

Fig. 2 Energy-dependent escape fractions F(ν)∕F⋆(ν) of photons above the ionization threshold of hydrogen for the fluxes shown in Fig. 1. In general, the escape fraction rises above each of the strong ionization edges due to the decreasing absorption cross sections for higher photon energies. The He II-ionizing photons are almost completely absorbed within the H II regions for all considered stellar sources; for Teff = 30 000 K almost no He I-ionizing radiation can escape from the H II region either. For higher metallicities, in addition to the prominent ionization edges of H I, He I, and He II, there is also notable absorption by metal ions, especially of O II, Ne II, N III, and C III with ionization energies of 2.59, 3.35, 3.49, and 3.52 Ryd, respectively. The dashed lines indicate the corresponding angle-averaged escape fractions for the case of extremely clumpy H II regions that would allow radiation to escape unhindered in some directions while being completely opaque to all wavelengths in others (“picket-fence” case). The spike in the escape fractions of the upper right panel is a result of the extremely low stellar He II-ionizing photon flux, which is approximately 12 orders of magnitude lower than the corresponding hydrogen-ionizing flux. Thus, even very little emission of photons by the H II region in that energy range is sufficient to exceed the flux emitted by the stellar atmosphere.

2.4. Line ratios in the matter-bounded (leaky) H II regions

In Tables 1 and 2 we show the predicted luminosities of collisionally excited line radiation (relative to H β) both of the matter-bounded H II regions and of the DIG irradiated by the ionizing fluxes escaping from such H II regions.

For the matter-bounded H II regions we find that the ratio of the lines emitted by the singly ionized atoms of N, O, and S relative to the Hβ emission increases for lower escape fractions. The reason for this behavior is that these ions occur predominantly in those parts of the H II regions where the optical depth to the emitting source is large, i.e., near the Strömgren radius of a radiation-bounded region. The increase is small for the H II regions ionized by the 30 000 K model stars,because as a consequence of their relatively low temperature, these stars emit only a very small number of photons capable of creating doubly ionized states of those elements (cf. Table 4 of Weber et al. 2015). Thus, for the 30 000 K stars the singly ionized atoms are abundant even in the inner parts of the H II regions.

Another factor that contributes to enhanced emission from the singly ionized species for lower escape fractions is that the strongest radiation hardening arises for the largest optical depths for ionizing radiation, which occur only for low escape fractions. The hardening results in increased heating rates and gas temperatures, and, consequently, an increased emission of collisionally excited lines in the optical wavelength range. In the case of the metal-poor H II regions with Z = 0.1 Z_⊙ the latter effect is compensated by the enhanced cooling via collisional excitation of neutral hydrogen, which is more abundant in the outer parts of H II regions with low escape fraction¹¹.

The [O III]∕Hβ ratio as a function of the escape fraction strongly depends on the effective temperature of the ionizing sources. For the 30 000 K stars, this ratio is ≤10⁻³ for all metallicities and all escape fractions because of the very low O II-ionizing fluxes of stars at that temperature, which lie more than 4 orders of magnitude below the hydrogen-ionizing fluxes (Weber et al. 2015). For an effective temperature of 35 000 K the [O III]∕Hβ ratio drops for lower escape fractions, as the presence of O III ions is limited to a smaller volume around the source compared to ionized hydrogen responsible for the H β recombination radiation. For the model stars with 40 000 and 50 000 K the [O III]∕Hβ ratio as a function of the escape fraction is metallicity dependent. The ratio decreases for lower escape fractions in the case of a (nebular) metallicity of 0.1 Z_⊙ whereas it increases for the metallicities of 0.4 Z_⊙ and 1.0 Z_⊙. This is caused by the different radius dependence of the temperature in the H II region. The temperature – and consequently the excitation of optical O III lines – (slightly) decreases toward larger radii for Z = 0.1 Z_⊙, which is caused by the cooling by collisionally excited neutral hydrogen, which is more abundant at larger optical depths in the outer parts of the ionized region. In contrast, the temperature increases with radius owing to radiation hardening for the metallicities of 0.4 Z_⊙ and 1.0 Z_⊙. Thus, for the higher metallicities the emission of the collisionally excited lines (relative to H β) of those ions still present near the Strömgren radius is enhanced in comparison to that of matter-bounded H II regions with larger escape fractions.

2.5. Line ratios in the DIG

For the DIG regions we find that the strength of the [N II] and [O II] emission lines increases relative to H β for decreasing escape fractions in those models in which the effective temperature of the source is 30 000 or 35 000 K and the metallicity is 0.4 Z_⊙ or 1.0 Z_⊙. This is caused by the higher temperatures resulting from the radiation hardening effect. In contrast, the [S II]∕Hβ emission line ratio overall tends to decrease for diffuse gas ionized by radiation escaping from the optically thicker H II regions, but the line ratio as a function of the escape fraction is often nonmonotonic and reaches it highest values not in the picket-fence case, but for DIG irradiated by radiation passing through H II regions with escape fractions of 30–50%. This is due to two opposing effects on the abundance of singly ionized sulfur in the ionized volume. Since the ionization energy of singly ionized sulfur (23.3 eV) is slightly lower than that of neutral helium (24.6 eV), the fraction of radiation in the interval between these ionization energies relative to the total hydrogen-ionizing radiation rises as a result of radiation hardening and can lead to a reduction of the number of S II ions in the DIG ionized by the filtered radiation. On the other hand, lower escape fractions in H II regions around the model stars with T_eff ≤ 35 000 K also lead to a more pronounced He I/He II absorption edge, and consequently fewer photons with energies above the corresponding ionization threshold are available to ionize the remaining S II. The latter effect is mainly relevant for the 35 000 K model stars; for the 30 000 K model stars the emission of He I-ionizing radiation is very low, and thus the contribution from that wavelength range to the total S II-ionizing radiation is negligible.

For the radiation sources with effective temperatures of 30 000 K, almost no O II-ionizing radiation can escape from the H II regions, and therefore there is no [O III] line radiation from the DIG. The same is true for the H II regions around the 35 000 K model stars with small escape fractions of 10 or 5%. For larger escape fractions the [O III]∕Hβ ratio in the DIG rises quickly, and the largest ratio is reached for the largest escape fractions.

In contrast, for the higher stellar temperatures of 40 000 and 50 000 K the [O III]∕Hβ ratio increases for decreasing escape fractions, as for these sources O II-ionizing photons can leave the H II regions even for small escape fractions of 10 or 5%. Actually, the fraction of O II-ionizing photons of all hydrogen-ionizing photons that can escape from the H II region is larger than the corresponding fraction of ionizing photons originally emitted by the sources, and the ratio increases for smaller escape fractions due to radiation hardening. For instance, we find an increase of approximately 0.6 dex for model stars with an effective temperature of 40 000 K and an escape fraction of 10%. Additionally, the radiation hardening leads to increased temperatures in the DIG. The radiation hardening also leads to decreasing [N II] and [O II] emission in the DIG around H II regions with small escape fractions for stellar temperatures of T_eff ≥ 40 000 K. In this case the effect of the larger fraction of photons able to strip a further electron from these ions is more important than the temperature increase due to radiation hardening.

The situation is different for S II, where the emission increases for lower escape fractions. The ionization energy of neutral sulfur lies below the ionization threshold of hydrogen, while the energy required for doubly ionized sulfur lies below the ionization threshold of neutral helium. Around hot O stars with strong He I-ionizing fluxes, S II ions therefore exist predominantly in the neutral phase of the medium. At the same time, electrons are required for collisional excitation processes. Most of the S II emission therefore takes place in the transition region between the ionized and the surrounding neutral phase. This region is broadened by radiation hardening as higher photon energies lead to reduced ionization cross-sections, allowing the photons to penetrate deeper into the neutral medium near the Strömgren radius.

In addition to the DIG models in which the SED of the ionizing radiation is modified by dense H II regions in the stellar environment, we also considered the case in which the shape of the ionizing spectrum remains unchanged. This case corresponds to either of the following two scenarios: (a) as a limit case for very high escape fractions, and (b) in situations in which the only absorbers are dense clumps that completely absorb the entire ionizing radiation and coexist with“channels” devoid of gas through which photons can escape unhindered from the star (see above). The results of these simulations are presented in Tables 1 and 2 as the picket fence case. In general, the simulated line ratios in the picket fence case continue the trend of the line ratios for higher escape fractions into the DIG at the respective temperature and metallicity. Even if the stellar radiation from the T_eff = 30 000 K model stars is not filtered, the production of photons above the ionization edge of O II is too low to cause a significant amount of [O III] emission from the DIG. For the models with T_eff = 35 000 K sources, the emission lines from N II and O II are weaker than for the filtered spectra, which is caused by the lack of radiation hardening and the resulting cooler temperatures of the gas. The O III line, in contrast, is stronger than with filtered radiation, where for T_eff = 35 000 K a large fraction of the photons above the He I-ionization edge are consumed within the H II regions.

In the stellar temperature range of T_eff ≥ 40 000 K the trend observed for increasing escape fractions in the “filtering” case (i.e., weaker absorption), namely a reduced emission of the O III line, continues. There is, however, no further significant increase in the emission from singly ionized nitrogen and oxygen.

2.6. Comparison with observations and other theoretical work

As noted in the introduction, observations of the DIG in the Milky Way indicate a softer ionizing spectrum than for typical H II regions. For example, measurements by Madsen et al. (2006) of the H α and He I λ 5876 recombination lines and the collisionally excited lines [N II] λ 6583, [N II] λ 5755, [O III] λ 5007, and [S II] λ 6716 emitted by the warm interstellar medium of the Milky Way (in the direction of the Orion-Eridanus Bubble) show line ratios (in terms of energy) of [O III]∕Hα < 0.15¹². In combination with the weak He I λ 5876 recombination line emission (He I∕Hα ≤ 0.015) and the strong emission of [N II] and [S II] ([N II]∕Hα ~ 0.16–0.43; [S II]∕Hα ~ 0.13–0.29), this shows that there is only a small fraction of ions that need photons above the ionization edge of He I to be formed, indicating that the DIG is ionized by a soft radiation field. A comparison with Tables 1 and 2 shows that such conditions can be attained by ionizing radiation that is predominantly emitted by relatively cool O stars (T_eff ≲ 35 000 K).

For example, for our models using an unaltered spectrum of a 35 000 K star we find an [N II]∕Hβ ratio of 1.3 for a metallicity of 0.4 Z_⊙ and 1.4 for a metallicity of 1.0 Z_⊙. In contrast, the radiation hardening and removal of photons able to ionize He I leads to considerably increased [N II]∕Hβ ratios of up to ~3.4. While these clearly exceed the values in the numerical model DIG by Barnes et al. (2014), who also used the SED of a 35 000 K star as the ionizing flux, such high values have indeed been observed in the diffuse ionized gas associated with the Perseus arm by the WHAM survey (Haffner et al. 1999, 2003, and Madsen et al. 2006; collected data have been presented by Barnes et al. 2014). Regarding the [S II] emission, our simulation results of [S II]∕Hβ ~ 0.6 are considerably lower than the simulation results from Barnes et al. 2014 (who found values on the order of 2.0), but are in better agreement with the observational values both in the DIG associated with the Perseus arm and the interarm DIG. The results for the DIG models using a 30 000 K SED are similar, but the [N II]∕Hβ ratio is less dependent on the escape fractions of the filtering H II regions and remains below the H α ∕Hβ ratio. This is not surprising because in both cases most of the ionized hydrogen volume predominantly contains ionization stages requiring photons below the He I ionization edge to be created. Thus, the results are in agreement with the low [O III]∕Hβ ratios found in the DIG of the Milky Way.

Our results are also in agreement with previous work that connects the ionization structure of the DIG with the radiation from stars at the cooler end of the O star temperature range. For example, Reynolds & Tufte (1995) called attention to earlier observationsby Torres-Peimbert et al. (1974), who found that the fraction of hot massive stars that are no longer embedded in H II regions rises for late-type O stars and B stars, i.e., stars with relatively low temperatures and, consequently, low fluxes above the ionization edge of neutral helium. This would signify higher escape fractions for the radiation from these objects and thus the ionization structure of the DIG would be dominated by the ionizing radiation of these objects. However, Reynolds & Tufte (1995) also point out that such objects contribute only about one fifth of the total stellar hydrogen-ionizing radiation, i.e., most of the radiation is contributed by hotter stars. To ionize the DIG, most of the ionizing radiation from these cooler stars would therefore have to remain unabsorbed by gas in their immediate environment.

For comparison with the DIG, Madsen et al. (2006) also observed “classical” H II regions around hot stars, where they have found strongly varying [O III]∕Hα ratios ranging from 0.025 to 0.45, i.e., the highest ratios found in the sample of H II regions exceed the highest values found in the sample of DIG subregions by a factor of approximately four¹³.

A relatively weak emission of doubly ionized oxygen ([O III]∕Hα ~ 0.20) has also been found in the DIG of M 31 by Greenawalt et al. (1997), although there the discrepancy in the [O III] emission (normalized to the recombination lines) between DIG and H II regions is smaller than in the Milky Way. These line ratios show the best agreement with our model DIG around the 35 000 K stars, but indicate a slightly harder ionizing spectrum than found in the nearby DIG of the Milky Way. We note that, as can be seen from Tables 1 and 2, even a slight variation in the effective temperatures of the ionizing sources or the escape fraction may drastically change the [O III] emission.

In contrast, an ionizing spectrum dominated by stars with effective temperatures of T_eff = 40 000 K or above would cause O III to become the strongest emission line in the optical wavelength range. The radiation hardening and the resulting higher temperatures in the DIG would lead to enhanced [O III] ratios compared to the dense H II regions directly surrounding the ionizing stars, which contradicts the observational results obtained by Greenawalt et al. (1997) and Madsen et al. (2006).

However, the dominance of soft ionizing sources is not consistent with observations that show that at least in some galaxies the [O III]∕Hβ ratio rises in the dilute gas at larger distances from the hot stars near the galactic disk (Otte et al. 2001; Collins & Rand 2001). These observations are incompatible with our models with matter-bounded H II regions around the cooler O stars of T_eff ≲ 35 000 K because in these scenarios the harder radiation above the helium ionization threshold is restricted to the inner parts of the ionized volume of the diffuse gas, as can be seen from the results in Table 1. The situation might be different if not only the H II regions surrounding the stars, but also the DIG were matter-bounded, so that even for a softer SED of the ionizing sources radiation hardening may outweigh the absorption of photons above the He I ionization threshold.

Alternatively, these observational results could be explained by additional sources with high effective temperatures that are not as concentrated toward the galactic disk as massive O stars, for example the evolved stars (central stars of planetary nebulae and white dwarfs) proposed by Sokolowski & Bland-Hawthorn (1991) and Flores-Fajardo et al. (2011). Integral field observations by Lacerda et al. (2018), who find the DIG in galaxies with low current star formation rates to have higher [O III]∕Hβ ratios, support the idea that evolved stars may play an important role in the ionization of the more dilute parts of the DIG. In view of these results it appears unlikely that the rising [O III]∕Hβ ratios are mainly the result of hardening of the radiation of hot massive stars. Our results in Sects. 2.3 and 2.5, which are quantified in Table 2, show that for stellar effective temperatures of T_eff ≳ 40 000 K, He I-ionizing photons are able to leave the inner H II regions even for low escape fractions of the total ionizing radiation. In these cases the effect of radiation hardening in the energy range between the ionization thresholds of He I and He II leads to an enhanced [O III]/[O II] ratio in the DIG. But hot massive stars are extremely rare in early-type galaxies (ellipticals and gas-depleted S0 disk galaxies with little ongoing star formation), while there still exists a population of evolved low-mass hot stars. Therefore an ionization by these objects appears to be more plausible.

3. Simulations of the diffuse ionized gas in 3D

In star-forming galaxies there are typically multiple regions containing hot stars that may contribute to the ionization of the DIG. Simulations of the DIG should therefore be able to account simultaneously for several sources located at different positions and embedded in an inhomogeneous interstellar gas. To realize this aim we applied our 3D radiative transfer code to model the DIG above a spiral arm. The underlying methods of this code are described in detail by Weber et al. (2013), who delineated the geometrical aspects of the ray tracing and the computation of the ionization structure of hydrogen and helium, and Weber et al. (2015), who extended the method with the computation of the ionization structure of metals, the temperature structure of the gas, andthe line emission. As in the spherically symmetric models, we assume the ionized gas to exist in two forms: H II regions in the immediate environment of the stars (unresolved in these simulations) and the dilute and extended DIG.

It is not our aim in this work to exactly reproduce the line ratios of any particular spiral galaxy. Instead, we want to know how varyingthe inhomogeneity in the DIG affects the ionization structure and line emission of the diffuse gas. Our approach is complementary to that of Wood et al. (2004, 2005), who concentrated on simulating the escape of ionizing radiation from H II regions in the stellar environment, whereas we focus on inhomogeneities within the DIG itself. Our approach is similar to that of Wood et al. (2010) and Barnes et al. (2014, 2015), who performed photoionization simulations of the DIG based on a density structure obtained from hydrodynamical simulations, or, in the latter case, a fractal density structure as well, while we choose a simpler description, but study how the clumpiness of the DIG affects its observable properties.

3.1. Method: three-dimensional radiative transfer based on ray tracing

The difficulty in modeling the DIG is related to the inherent lack of symmetry of the problem; the presence of multiple ionization sources coupled with the inhomogeneous gas distribution can create a complex three-dimensional pattern in the radiation field. To compute that pattern numerically a suitable discretization scheme is required. The simplest but still fully general such method is to divide up the simulation volume into a Cartesian grid of cells, where the gas density may vary from cell to cell but the ionization structure and density within each cell is taken to be homogeneous. The accuracy and detail level of the simulation may be then chosen by selecting an appropriate cell size.

Because the spatial extent of the ionized region exceeds stellar radii by typically more than six orders of magnitude, the sources of ionizing radiation can be considered to be point-like, and the radiative transfer can be computed along rays emanating from these points. Limiting the placement of the sources to the centers of cells results in no significant further loss of detail beyond the already effected discretization of the volume into cells, but allows conserving computational resources in calculating the radiation field, because then the intersections of the rays with the cell boundaries relative to the source are the same for all sources and need to be computed only once; the results are stored and then accessed via a lookup table. The rays themselves are distributed evenly in all directions around a source (using the tessellation scheme of Górski et al. 2005), each ray representing a fixed solid angle ΔΩ = 4π∕N_rays and thereby a given fraction L_ray = L∕N_rays of the luminosity L of the source. The geometric dilution of the radiation field is automatically provided for by the decreasing number of rays intersecting each cell with increasing distance from the source.

Each source may have a different luminosity and SED, and the radiative transfer is computed for a number of frequency points (typically 100–1000) sufficient to describe the energy distribution of the ionizing flux. We account for the diffuse radiation field (the ionizing radiation emitted by recombination processes within the gas) by applying the case-B approximation for hydrogen and helium, i.e., the assumption that the photons emitted by recombination directly into the ground state are reabsorbed locally. In the context of our 3D radiative transfer, the accuracy of this approximation is quantitatively discussed by Weber et al. (2013).

As a further step in our simulation procedure the time-dependent ionization fractions n ≡ n(r, t) of all ionization stages i of the elements considered are calculated on the basis of the rate equation system $$\frac{\text{d}}{\text{d}t}{n}_i(t)={\displaystyle \sum _{j\ne i}}{\mathcal{P}}_{ji}{n}_j(t)-{\displaystyle \sum _{j\ne i}}{\mathcal{P}}_{ij}{n}_i(t),$$(2)

which describes the rate of change of the number density of all ionization stages i. The rate coefficients ${\mathcal{P}}_{ij}$ are the sums of the radiative (${\mathcal{R}}_{ij}$) and the collisional (${\mathcal{C}}_{ij}$) transition rates. Because the strengths of emission lines depend on the temperature, we must also account for the heating and cooling processes that determine the temperature structure of the gas.

In the present work we restrict ourselves to the computation of the special case of the stationary state¹⁴ for the temperature and ionization structure (i.e., the case in which the left-hand side of Eq. (2) becomes zero), which is consistentwith our assumption (see below) of continuous ongoing star formation. For otherwise fixed conditions (density structure of the gas, SEDs of the ionizing sources), the stationary equilibrium case would be reached for t → ∞¹⁵.

To determine this equilibrium, it is necessary to compute the ionization and recombination rates and the thermal equilibrium. The photoionization rate coefficients ${\mathcal{R}}_{ij}(m)$ from a lowerstage i to an upper stage j are calculated for a grid cell m as $${\mathcal{R}}_{ij}(m)={\displaystyle {\int }_{{\nu }_{ij}}^{\infty }\frac{4\pi a_{ij}(\nu )}{h\nu }}{J}_{\nu }(m)\text{d}\nu ,$$(3)

where J_ν(m) is the mean intensity in the cell, ν_ij the threshold frequency for the corresponding ionization edge, and a_ij(ν) the frequency-dependent ionization cross section. Similarly, the radiative heating rates, i.e., the gain of thermal energy per unit volume due to ionization of state i are computed as $${\Gamma }_{ij}(m)=n_i(m){\displaystyle {\int }_{{\nu }_{ij}}^{\infty }\frac{4\pi a_{ij}(\nu )}{h\nu }}{J}_{\nu }(m)h(\nu -{\nu }_{ij})\text{d}\nu ,$$(4)

where n_i(m) is the number density of ionization stage i in cell m.

The link between the radiative transfer and the photoionization and heating rates is the mean intensity J_ν (m), which can be computed by making use of the fact that the number of ionization processes within a grid cell must equal the number of ionizing photons absorbed in the grid cell (see Weber et al. 2013). This leads to the expression $$J_{\nu }(m)=\frac{1}{4\pi V(m)}{\displaystyle \sum _{\begin{array}{c}\text{ray}\\ \text{segments}\end{array}}}\frac{L_{\text{ray},\nu }^{\text{inc}}(1-e^{-{\tau }_{\nu }})}{{\chi }_{\nu }(m)},$$(5)

where V (m) is the volume of the cell, $L_{\text{ray},\nu }^{\text{inc}}(m)$ the frequency-dependent incident luminosity of a ray entering the cell (obtained from the radiative transfer solution), τ_ν (m) the optical depth of the ray segment in the cell, and χ_ν(m) the opacity of the material. To obtain the mean intensity, the contributions of all rays crossing the grid cell are added up.

The recombination rate coefficients in grid cell m are ${\mathcal{R}}_{ji}(m)={\alpha }_j(T(m))n_{\text{e}}(m)$, where n_e and T are the electron density and temperaturein the cell and α_j(T) is interpolated from tabulated, temperature-dependent recombination coefficients containing both radiative and dielectronic contributions. For hydrogen and helium we use the tables of Hummer (1994) and Hummer & Storey (1998); for the metal ions we use the fit functions given by Aldrovandi & Pequignot (1973, 1976), Shull & van Steenberg (1982), Nussbaumer & Storey (1983), and Arnaud & Rothenflug (1985), as compiled by D. A. Verner¹⁶.

To compute the temperature, we consider ionization heating and recombination cooling by hydrogen and helium, as well as free–free cooling and cooling via collisional excitation and radiative de-excitation of bound states. The latter includes cooling through the collisionally excited optical lines of N II, N III, O II, O III, O IV, Ne III, S II, and S III and the fine-structure lines of N II, N III, O III, O IV, Ne II, Ne III, S III, and S IV.

The computational effort required to create 3D models of the photoionized gas considerably exceeds the effort for spherically symmetric models. For instance, in our three-dimensional simulations the number of volume elements is approximately 10⁶ and we consider the interaction of the matter of the DIG with 50 sources of ionization, whereas in the spherically symmetric models there is a single source of ionization and the temperature and ionization structure is computed for just 101 radius points. To keep the run time and memory requirement of the 3D simulations at a feasible level, we utilize a simplified modeling of the occupation numbers in the ionized gas, in which we use the “nebular approximation”. That is, we assume that the timescales for all processes that require interaction between two particles, such as ionization (requiring interaction between an atom/ion and a photon) and recombination and collisional excitation (requiring interaction between an atom/ion and an electron) are long compared to the timescale for the relaxation of an excited atom/ion into the ground state. Thus, practically all atoms/ions that undergo interactions with photons or electrons are in their ground state and only transitions from the ground state need to be considered. As a consequence, only the abundances of the different ionization stages have to be computed explicitly, whereas in the spherically symmetric case we explicitly compute the occupation numbers for up to 50 excited levels per ion.

For the recombination processes, which can involve various excited intermediate stages during the recombination from the ground state of the higher to the ground state of the lower ionization stage, the 3D method is based on fit formulae that describe the total recombination rate, i.e., the sum of the recombination rates to all levels.

For the densities considered in this work (n_e ≤ 10 cm^-3) the nebular approximation is a good representation of the actual situation in the gas. This is because, even for forbidden transitions such as the ¹D →³P transitions of N II and O III, the fine-structure transitions between the ³P sublevels of these ions, and the ²D →⁴S transitions of O II and S II, spontaneous radiative relaxation is much more probable than a two-particle interaction in the excited upper state¹⁷.

3.2. Simulations of the DIG around a spiral arm

With our 3D simulations, we study the ionization ratios in the DIG above and below a section of the galactic disk corresponding to a spiral arm. In contrast to our parameter study in the previous section using spherically symmetric models and single-star spectra, in which we examine how the line ratios vary with metallicity and hardness of the ionizing spectrum, the intent of the simulations in this section is to examine how the ionization structure of the DIG varies as a function of small-scale inhomogeneities in the gas. For this purpose we choose a more realistic SED from a stellar population synthesis model¹⁸ as a common ionizing spectrum for the different clumping scenarios we consider.

In each of our simulations, the stellar population that emits the ionizing radiation is represented by 50 sources, whose positions roughly resemble the distribution of luminous H α emitters in a spiral arm of M 81 observed by Greenawalt et al. (1998). The distribution of distances of the sources from the galactic plane follows an exponential law with a scale height of 90 pc (Bahcall & Soneira 1984). To determine the ionizing photon emission rate of the sources, we take the ionizing luminosity predicted by the stellar population synthesis model for the star formation rate of a typical star-forming galaxy, scale this luminosity to the fraction of the galactic disk represented by our simulation volume, and then divide it by the number of sources in our simulation. The resulting hydrogen-ionizing flux emitted by each source is 7.1 × 10⁵⁰ photons persecond, and the emission rate of photons able to ionize He I and He II is 1.1 × 10⁵⁰ and 2.1 × 10⁴⁶ s⁻¹, respectively.

Specifically, these values were obtained as follows. We assume a star-forming disk galaxy with a continuous star formation rate of approximately 1 M_⊙ yr⁻¹ like the Milky Way (Robitaille & Whitney 2010) or M 81 (Gordon et al. 2004) and a scale length of the disk of approximately 2 kpc (Porcel et al. 1998), which for a Kroupa IMF with an upper mass limit of 100 M_⊙ in total emits 1.5 × 10⁵³ hydrogen-ionizing photons per second. The ionizing flux was then adapted to a simulated volume of ${(5\text{kpc})}^3$, assuming that the center of the simulated volume has a distance of 7.5 kpc from the galactic center (similar to the distance between the galactic center and the position of the sun) and the result multiplied by 10, as we assume an enhanced star formation rate in a spiral arm of a star-forming galaxy. This star formation rate was then used as input to the population synthesis code Starburst99 (Leitherer et al. 1999, 2010, 2014; Vazquez et al. 2005), in which we chose the recent Geneva 2012 stellar evolutionary tracks (Ekström et al. 2012, models without rotation) at solar metallicity. The luminosity computed for the simulation volume is then distributed evenly¹⁹ among the 50 sources we assumed in our simulations.

In a first benchmark simulation we assume that the entire (unfiltered) radiation is absorbed by H II regions with a hydrogen number density of 2 cm⁻³. Even with this (for H II regions) comparatively low density, the radii of the ionized regions do not exceed 150 pc. These H II regions are shown in the bottom left panel of Fig. 3, where we compare the visual appearance of both the H II regions (left-hand side) and the DIG (right-hand side) from our 3D model with an observation by Greenawalt et al. (1998). Among our simulations of the DIG we find the best agreement to the observations in a model that combines clumping at various scales (see below).

For the density of 10 cm^-3 assumed in our spherically symmetric computations in Sect. 2, the radii of the H II regions are on the order of the spatial grid cell size of Δx ~ 50 pc used in the 3D simulations in this section²⁰, and thus would not be resolved. For the remainder of this section, we concentrate on the diffuse gas and treat the H II regions around the ionizing sources as unresolved (subgrid) point sources that are not considered explicitly, but whose presence is implicitly accounted for by using the filtered stellar population SED as the ionizing flux for the DIG. Using our procedure of Sect. 2, we compute the filtered spectrum assuming a homogeneous H II region with n_H = 10 cm^-3 that allows a fraction f_esc = 10% of the ionizing flux to escape.

To study the effects of clumping and porosity on the ionization structure of the extraplanar gas, we perform a series of five simulations of the ionization of the DIG for a hydrogen number density that is similar to the values found in the DIG of the Milky Way. For all models, the mean density of the gas in the galactic plane is n_H = 0.1 cm^-3 and decreases exponentially with a scale height of 1 kpc away from the galactic plane (cf. Haffner et al. 2009 and references therein). After filtering by their H II regions, each of the sources provides an ionizing flux of Q_H = 7.1 × 10⁴⁹ s^-2 that can penetrate into the diffuse gas. As there are 50 ionizing sources in the simulated volume and the cross section of that volume is 25 pc², the total ionizing flux into the DIG is approximately 1.4 × 10⁵⁰ s^-1 kpc^-2. This value issimilar to the 1.6 × 10⁵⁰ s^-1 kpc^-2 found by Barnes et al. (2015) to ionize the DIG in their model using a fractal density structure.

We note that this treatment is somewhat simplified because, as shown in Sect. 2.3, the filtering of the ionizing radiation depends on the extent of the H II regions and the included stellar population. This means that the spectrum of the radiation that would escape if each of the stellar sources were individually surrounded by an H II region with a given escape fraction is not the same as the spectrum of the radiation escaping from a H II region surrounding a large association of stars. Compared to the energy distribution of the sources, the emergent spectrum is mainly modified by radiation hardening. Whereas forthe sources themselves the number ratio of photons able to form He II to those able to form H II is Q_He I∕Q_H I = 1: 6.7, this ratio rises to 1: 3.75 for the escaping radiation²¹. While the ratio in the unfiltered case is similar to that found in 40 kK stars, the radiation hardening in the filtered case leads to a ratio similar to that found in 45 kK stars²². We show the SED of the filtered flux from the population synthesis model entering into the DIG in Fig. 5.

For the density fluctuations in the diffuse gas we consider five different scenarios:

• 1.

  In the first simulation, the density structure of the gas is smooth, i.e., apart from the exponential drop of the hydrogen density for increasing heights above the galactic plane there are no inhomogeneities. We note that the smooth density distribution refers only to the diffuse gas itself; we still assume the same filtering of the ionizing radiation by gas in the environment of the sources as in the other simulations.

• 2.

  In the second simulation, the gas has a clumping factor of $f_{\text{cl}}=\langle n_{\text{H}}{}^2\rangle /{\langle n_{\text{H}}\rangle }^2=5$, which is similar to the value found in the observations by Reynolds (1991). We assume that 20% of the volume is filled with gas that has a density of five times the mean density at the corresponding galactic height, whereas the rest of the volume is empty. (Thus, in the galactic plane the density of the clumps is n_H = 0.5 cm^-3.) Here we use the microclumping approach, in which the clumps are assumed to be optically thin and much smaller than the resolution elements of our simulation. This microclumping is simulated numerically by multiplying the rate coefficients for the processes involving ion-electron interactions, namely recombination, collisional excitation, and thermal free-free radiation, by the clumping factor, as for otherwise equal conditions the electron number densities in the filled parts of the gas scale with the clumping factor.

• 3.

  In the third simulation, we replace the microclumping approach by the assumption of larger scale inhomogeneities. This macroclumping is realized by filling only 20% of the grid cells, selected randomly and independently of each other, with gas that has a density of five times the mean density at the corresponding galactic height, leaving the rest of the cells empty (resulting in a typical length scale of the inhomogeneities that roughly corresponds to the cell size of 50 pc). While we are certainly aware that this is not necessarily a realistic characterization of the actual statistical distribution of densities in a nonsmooth DIG, we chose this description for two reasons. First, it is a natural extension of the unresolved (subgrid) microclumping model to a spatially resolved description, as both assume a two-phase medium consisting of swept-up matter and voids. Second, it is the most extreme clumping for a given spatial clumpiness scale. This and the smooth model represent the corner cases between which a more graduated description with intermediate density values would lie.

• 4.

  For comparison, we also performed a simulation in which we increased the (macro)clumping factor to f_cl = 25, i.e., only one out of 25 cells is filled with gas having 25 times the mean density (n_H = 2.5 cm^-3 in the galactic plane), while the rest of the cells remain empty.

• 5.

  Finally, we simulated a combination of micro- and macroclumping, in which both the total amount of gas within each of the grid cells corresponds to the macroclumping case with a clumping factor of f_cl = 5 and the internal structure of the gas within the single cells is characterized by an additional microclumping factor of 5.

We show both face-on and edge-on emission measure maps²³ of the simulation results in Fig. 4 and the ionization fractions of hydrogen as functions of the height above the galactic plane in Fig. 6. In the case of the homogeneous gas, ionized hydrogen exists within entire simulation volume, which reaches up to 2.5 kpc above the galactic plane. The ionization fraction actually rises at larger heights above the galactic plane due to the lower hydrogen densities, which lead both to lower recombination rates and a lower opacity of the gas and therefore larger mean free paths of the ionizing photons.

The increased recombination rates in the denser clumps lead to a smaller ionized volume in both the microclumping and macroclumping case with f_cl = 5. We find that in these cases the ionization fraction drops to approximately 10% (20% of the maximal value we obtain near the galactic plane) at a height of 1.2 kpc with only small differences between both cases. While clumping increases the recombination rates it also increases the fraction of the volume with below-average densities. Thus, channels may form²⁴ that allow the ionizing radiation to reach larger distances from its source. This causes the gas to have an ionization fraction of approximately 25% even at large heights above the galactic plane for the f_cl = 25 case.

The combined micro-and macroclumping case is characterized by the lowest ionization fractions among the simulations. The recombination rates in the filled parts of the volume are equal to the f_cl = 25 case described above, but as a result of the larger number of filled cells there are fewer channels for the radiation to leak out to large distances from the galactic plane.

Figure 6 also visualizes the results of a further parameter study in which we assumed the microclumping scenario and varied the clumping factor between 1 (the unclumped case) and 5 (the microclumping case presented before). The mean densities and ionizing sources are the same as in the previous simulations. In the parameter study the ionization fraction of the gas as a function of the galactic height strongly depends on the clumping factor of the gas. For higher clumping factors, the ionization fraction decreases for larger galactic heights until the gas becomes neutral. In contrast, for low clumping factors the ionization fraction at large galactic heights is larger than within the galactic plane and the ionized region becomes matter-bounded, i.e., ionizing radiation can escape into the intergalactic space. For intermediate clumping factors, the ionization fraction drops until a minimum is reached and then rises again for larger distances from the central plane. The height-dependence of the ionization fraction results from several competing effects. On the one hand, photons have to penetrate the denser, more opaque layers of gas to ionize high-altitude gas, resulting in a decreasing number of available ionizing photons with increasing height. On the other hand, the recombination rate is lower in the thinner gas above the galactic plane, which increases the equilibrium value for the ionization fraction at higher altitudes. Furthermore, the opacity of the thin gas at large galactic height is lower than the opacity near the galactic plane. A comparison between the results for a clumping factor of 3 and for a clumping factor of 4 shows that the transition between clearly matter-bounded and clearly radiation-bounded ionized gas occurs within a relatively small range for the clumping factor. It can be expected that the ionization structure of the gas at large galactic heights is also strongly affected by slight variations of parameters of the gas that influence the recombination rates, such as the mean density or the temperature, or by a small variation of the amount of ionizing radiation. We note that the numerical values shown in the plots represent lower limits for the ionization fractions, as we assume open boundary conditions, thus neglecting photons entering the volume from the sides.

For the simulations with a clumping factor of f_cl = 5 and the simulation with combined micro- and macroclumping, we show line ratio maps of the diagnostically important optical lines [N II], [O II], and [O III] relative to Hβ in the face-on view in Fig. 7 and in the edge-on view in Fig. 8. For the reasons discussed in Sect. 2, higher temperatures due to radiation hardening and a larger fraction of singly ionized nitrogen and oxygen, the [N II]∕Hβ and [O II]∕Hβ line ratios reach their maximal value in the outer regions of the ionized gas.

Barnes et al. (2014) show both the [N II]∕Hα ratios from their simulations and the corresponding observations of the DIG around the Perseus arm and the interarm medium. For simulations that do not assume any heating process other than photoionization, these authors obtain [N II]∕Hα ratios of approximately 0.4. The Hα∕Hβ ratio being approximately 2.9 without strong dependence on the gas density, this implies an [N II]∕Hβ ratio close to unity. Additionally, Barnes et al. (2014) have performed simulations including an additional heating term proportional to the electron density, G_add = n_e × 1.5 × 10⁻²⁷ erg s^-1. For these models they found almost unchanged emission line ratios for the relatively dense gas, but an increase of the [N II]∕Hα ratio to values close to or even above unity for more dilute gas.

Our 3D models do not directly correspond to the models of Barnes et al. (2014), both with respect to the ionizing radiation field and the density distribution of the gas. In particular, the population synthesis model also considers hotter stars than the 35 000 K model star used by Barnes et al. (2014) to describe the radiation field, leading to a harder ionizing spectrum. As a result, in our case the effect of the lower N II (and higher N III) abundances near the galactic plane is more pronounced, and it is not surprising that we obtain stronger variations of the [N II]∕Hβ as a function of the galactic height. We find a strong rise of the [N II]∕Hβ ratio near the edge of the ionized region, where in the microclumping case values of up to approximately 8 are reached. These roughly correspond to the maximal observed values in the environment of the Perseus arm, but are clearly above the range of values found in the interarm medium. In the central parts of the ionized region, [N II]∕Hβ remains below two. For the macroclumping cases there is a less steep increase of the [N II]∕Hβ ratio toward the edge of the ionized volume.

We find a better agreement to the observational measurements of DIG shown in the same paper, especially regarding the rise of that ratio at larger galactic heights, but observations of the DIG of the Milky Way (Reynolds & Tufte 1995) suggest the absence of ionizing photons with energies beyond the ionization edge of neutral helium, which would also imply that there is very little doubly ionized nitrogen.

The [O III]∕Hβ ratio rises from the galactic plane to its maximal value at heights of up to 800 pc. For larger distances from the galactic plane the [O III]∕Hβ ratio decreases again. The outermost parts of the ionized gas do not emit [O III] line radiation. The overall O III ∕Hβ ratios obtained by our 3D simulations are clearly too high to be consistent with the low abundance of ions requiring photons above the ionization threshold of neutral helium to be generated as found, for instance by Reynolds (1985) and Madsen et al. (2006). This could be explained by a lower escape fraction from the galactic plane and a significant contribution of other heat sources, such as magnetic reconnection or shocks (cf. Sect. 1) or by higher escape fractions of ionizing radiation from H II regions around relatively cool ionizing sources (as discussed by Reynolds & Tufte 1995). Furthermore, the contribution of the most massive and therefore hottest stars may be lower than assumed in our simulation. There could be two causes of this reduced contribution: the high-mass end of the mass function could have a lower number of stars than in our population model, or these massive stars could form in larger regions of overdense gas, reducing the fraction of their radiation that is able to escape into the DIG.

In our simulations we also do not find an increase of the [O III]∕Hβ ratio at very large heights above the galactic plane as observed in edge-on galaxies by Otte et al. (2001) and Collins & Rand (2001). To resolve this discrepancy, additional sources of ionization at larger distances from the galactic plane, such as central starsof planetary nebulae or white dwarfs as proposed by Sokolowski & Bland-Hawthorn (1991) and Flores-Fajardo et al. (2011) and supported by the observational results for galaxies with low star formation rates by Lacerda et al. (2018), might need to be considered.

In the simulation with the homogeneous density (not shown in Figs. 7 and 8) the gas is almost fully ionizedand optically thin, and thus the effects of photons entering the volume from the sides would need to be considered. However, because we used open boundary conditions in our simulations such photons are not taken into account, thereby possibly leading to different results for the ionization structure of metals and consequently for the line emission. In the simulation with f_cl = 25, on the other hand, most of the grid cells containing gas are optically thick for the ionizing radiation, which means the ionization structure and properties of the radiation field would in reality vary significantly within the volume of a cell. Because a grid-based approach such as ours assumes that the temperature and ionization structure as well as the radiation field within each cell is constant, this could lead to inaccurate values for the computed emission line ratios.

Fig. 3 Comparison between Hα observations of H II regions and the DIG in the face-on galaxy M 81 by Greenawalt et al. (1998; middle row) and our simulation results (bottom row) of the ionization of the DIG by radiation escaping from the H II regions in the environment of the radiation sources. We find a reasonable qualitative agreement regarding both the denseand relatively intense H II regions around the star-forming regions on the left-hand side of the imageand the DIG on the right-hand side of the image. The simulation results presented in this image were computed for a combination of microclumping and macroclumping with a clumping factor of fcl = 5 each. The brightness values in these images as well as those in Fig. 4 represent an emission measure range between 1 and 750 pc cm-6 (logarithmicscaling). Observational images: ©AAS. Reproduced with permission.

Fig. 4 Edge-on and face-on emission measure maps for the different simulated clumping scenarios (see text). In each case we assumed a mean hydrogen number density of nH (h) = 0.1 cm−3 × e−|h|∕h0, where |h| is the distance from the galactic plane and h0 = 1 kpc is the scale height of the gas. The simulations also share the same distribution and luminosities of the ionizing sources. The edge lengths of the maps are 5 kpc. The macroclumping case with fcl = 5 appears generally very similar to the microclumping case with the same clumping factor, as is expected for sufficiently fine grid resolution. For fcl = 25, the macroclumping becomes very porous at the used grid resolution, opening up channels devoid of gas that allow radiation to escape from the immediate environment of the sources without being absorbed. The combined micro- and macroclumping case with fcl = 5 × 5 is very different from the macroclumping case with fcl = 25, giving an overall appearance resembling that of the cases with fcl = 5, but with theemission in each of the ionized bubbles more concentrated; this is somewhat analogous again to the transition from fcl = 1 to fcl = 5.

Fig. 5 Left panel: SEDs from the population synthesis model, filtered as the stellar atmosphere models in Fig. 1. For our 3D simulations we use the SED of the model with fesc = 10% (magenta) as ionizing flux for the DIG. Right panel: corresponding energy-dependent escape fractions, as in Fig. 2. A comparison with Figs. 1 and 2 shows that the SEDs roughly resemble those using a 40 kK star.

Fig. 6 Left panel: ionization fractions of hydrogen as a function of the height above the galactic plane for the different clumping scenarios presented in Fig. 4. For the homogeneous density distribution, the gas is predominantly ionized throughout the simulation volume and the ionization fraction rises with increasing galactic heights due to the lower gas densities, which results in decreasing recombination rates and opacities. Both the micro- and the macroclumping case with fcl = 5 are characterized by a similar dependence on the height above the galactic plane, reaching approximately 0.5 near the galactic plane and dropping to 0.1 at a height of about 1.2 kpc. In contrast, in the macroclumping scenario with fcl = 25 the ionization fraction varies only slightly throughout the simulation volume. The combination of micro- and macroclumping shows the lowest ionization fraction of all the scenarios owing to the increased recombination rates in the filled cells. While the total clumping factor equals that of the fcl = 25 macroclumping case, there are fewer extended empty regions through which radiation can penetrate to greater distances from the galactic plane, where less ionizing radiation is required to maintain the same fraction ofionized hydrogen due to the lower gas densities. Right panel: as before, but for a parameter study varying the microclumping factor between 1 and 5. The ionization fraction rises at larger galactic heights inthe case of low clumping, whereas it decreases and the gas eventually becomes neutral for high clumping factors. For immediate clumping factors, which for our geometry and density are around fcl ~ 3, the ionization fraction as a function of height behaves in a nonmonotonic way, first dropping to a minimum and then rising again.

Fig. 7 Line ratio maps of several collisionally excited lines relative to Hβ in the face-on view of the simulated spiral arm for different density structures of the gas. Gray indicates emission measures less than 1 pc cm-6.

3.3. Three-dimensional simulation assuming a softer ionizing SED

To test the effects of a softer ionizing radiation on the line emission, we repeated our simulation number 2 (the microclumping case with f_cl = 5) using as ionizing flux an unfiltered 35 000 K dwarf star SED. This choice of SED allows a better comparability with the work of Barnes et al. (2014, 2015) and with earlier work by Reynolds & Tufte (1995), who discuss the possibility that the escape fraction tends to be higher for the radiation of cooler O (and even B) stars. We present the computed line ratio maps for this case in Fig. 9.

Comparing the results from this simulation with those using the filtered population SED, we find a considerably reduced [O III] emission, which is concentrated close to the galactic plane and disappears for large heights above the galactic plane due to the absence of doubly ionized oxygen. Observations that find a rising emission of [O III] for larger distances from the galactic plane can therefore clearly not be explained by a cool ionizing spectrum. By contrast, the emission of [N II] and [O II] rises slightly faster as function of the height above the galactic plane as a result of radiation hardening, the higher abundances of singly ionized nitrogen and oxygen, and the lower abundance of the very efficient coolant O III. However, the [N II]∕Hβ and [O II]∕Hβ ratios do not attain as large values because the softer ionizing spectrum leads to cooler gas temperatures, which in turn leads to increased recombination rates and a smaller ionized volume: while in the model using the hardened spectrum of the filtered population SED we find a mean temperature (weighted by square of the electron density) of approximately 9300 K, the corresponding weighted gas temperature in the model using the spectrum of the 35 000 K star is 7800 K.

3.4. Comparison with an analytical model for the ionization of the DIG

In this section we compute the height of the layer of ionized gas analytically and compare with the results from the numerical solutions presented above. Hereby we assume a microclumped gas whose mean density declines exponentially as a function of the height above the galactic plane and a purely plane-parallel geometry of the radiation field, without discrete sources, but with a homogeneous flow of ionizing photons emerging from the central plane. For simplicity, we further assume that the gas consists of pure hydrogen with a constant temperature (7500 K) in the ionized volume. For computing the ionization structure we use the Strömgren volume approximation that the gas is either fully neutralor fully ionized.

Under these assumptions we can analytically compute the height z_max of the ionized layer of gas. In equilibrium, the emission rate Ṅ of ionizing photons equals the number of recombination processes per unit time. Normalized by surface area of the disk the equilibrium can be expressed as $$\begin{array}{ll}\frac{\dot{N}}{A}\hfill & ={\displaystyle {\int }_{-z_{\text{max}}}^{z_{\text{max}}}{\alpha }_{\text{B}}}{n}_{H\text{\hspace{0.05em}II}}(z)n_{\text{e}}(z)f_{\text{cl}}\text{d}z\hfill \\ \hfill & ={\displaystyle {\int }_{-z_{\text{max}}}^{z_{\text{max}}}{\alpha }_{\text{B}}}{n}_{H\text{\hspace{0.05em}II}}(0)e^{-\left|z\right|/h}{n}_{\text{e}}(0)e^{-\left|z\right|/h}{f}_{\text{cl}}\text{d}z\hfill \\ \hfill & =2{\alpha }_{\text{B}}{n}_{H\text{\hspace{0.05em}II}}^2(0)f_{\text{cl}}{\displaystyle {\int }_0^{z_{\text{max}}}{e}^{-2\left|z\right|/h}}\text{d}z\hfill \\ \hfill & ={\alpha }_{\text{B}}{n}_{H\text{\hspace{0.05em}II}}^2(0)f_{\text{cl}}h\left(1-e^{-2z_{\text{max}}/h}\right),\hfill \end{array}$$(6)

where α_B is the (case-B) recombination coefficient, n_H II the number density of ionized hydrogen and electrons, f_cl the clumping factor, and h the scale height of the gas. Equation (6) can now be solved for z_max, $$z_{\text{max}}=-\frac{h}{2}\ln \left(1-\frac{\dot{N}/A}{n_{H\text{\hspace{0.05em}II}}^2(0)f_{\text{cl}}h}\right),$$(7)

which becomes infinite for $$\dot{N}/A\ge {\alpha }_{\text{B}}{n}_{H\text{\hspace{0.05em}II}}{}^2(0)f_{\text{cl}}h,$$(8)

in which case the ionized gas layer around the disk is no longer radiation-bounded and ionizing radiation can escape into intergalacticspace.

For the parameters used in our numerical simulations, namely an ionizing flux of Ṅ = 3.55 × 10⁵¹ s^-1, a section of the galactic plane of ${(5\text{kpc})}^2$, and a (mean) central hydrogen density of 0.1 cm^-3, the resulting height of the ionizing layer is z_max = 221 pc when assuming microclumping with f_cl = 5 and a recombination coefficient of α_B = 2.8 × 10⁻¹³ cm³ s^-1. The clumping factor below which the ionization becomes matter-bounded for otherwise unchanged conditions is f_cl = 1.73. A comparison of these values with our simulation results shows that for a purely plane-parallel geometry the maximal height of the ionizing gas is lower than in the numerical models that assume discrete sources (see Fig. 6). The clustering of the sources in the spiral arm locally causes stronger fluxes of ionizing radiation than in the plane-parallel model, allowing the ionizing radiation to reach larger heights above the galactic plane. Conversely, in the 3D models there are also regions in the galactic plane devoid of ionizing sources, which is not the case in the simplified analytic model; the gas in these regions remains un-ionized as a result of the lack of nearby ionizing sources and the absorption of the radiation from the more distant sources by the dense gas at low galactic heights.

Fig. 8 As Fig. 7, but an edge-on view. Oxygen and nitrogen are predominantly doubly ionized close to the ionizing sources, and singly ionized only in the outer part of the ionized region. Thus, and because of radiation hardening and therefore higher temperatures, the maximum of the line ratios [N II]∕H β and [O II]∕Hβ occurs at larger heights above the galactic plane. For the cases with a clumping factor of fcl = 5 the [O III]∕Hβ ratio rises up to a height of approximately 700 pc above the galactic plane and drops again at larger distances.

Fig. 9 Line ratio maps from microclumped models (fcl = 5) as in the leftmost columns of Figs. 7 and 8, using the same simulation volume and the same integrated hydrogen-ionizing flux as before, but taking the unfiltered spectrum of a 35 kK dwarf star with solar metallicity as the ionizing SED. The most notable difference to the models using the population SED is the weaker emission of the optical [O III] lines. The emission of the [N II] and [O II] lines is stronger close to the galactic plane, but shows a weaker rise close to the edge of the simulated volume.

4. Summary and outlook

The diffuse ionized gas is the major contributor of warm ionized interstellar matter in disk galaxies, but the mechanisms that keep it ionized and determine its energy and density structure are not yet fully resolved. While a number of different ionization scenarios have been suggested, in this work we concentrated on the ionization of the diffuse interstellar gas by radiation from massive hot stars. Because of their short lifetimes, these objects are associated with star formation regions and are therefore concentrated near the galactic disk. Thus, their ability to ionize the diffuse gas depends on the ability of the emitted radiation to penetrate through the interstellar gas close to the disk.

When the radiation passes through the denser gas surrounding the sources, its attenuation is a function of the frequency-dependent absorption cross sections of the atoms and ions contained within that gas. Thus, the SED of the radiation that penetrates into the DIG can be significantly different from the SEDs of the stellar content of the traversed H II region. The modification of the SEDs depends on the original SED, the escape fraction, and the density distribution of the gas in the environment of the sources, most notably of hydrogen and helium, but also – to a lesser extent and depending on the metallicity of the gas – of carbon, nitrogen, and oxygen ions.

We studied the dilution effects using both a sophisticated radiative transfer scheme that assumes the problem to be spherically symmetric, and a more simplified 3D-approach that, however, dispenses with the requirement of spherical symmetry. For the spherically symmetric case we computed model grids of the line emission from the matter-bounded H II regions and from the diffuse gas ionized by radiation escaping from these H II regions. For these grids we used SEDs from state-of-the-art stellar atmosphere models, similar to those that have already previously shown consistently good results as ionizing sources in numerical models of H II regions that were compared with observations. We vary the effective temperature (30 000 K ≤ T_eff ≤ 50 000 K) and metallicity (0.1 Z_⊙≤ Z ≤ 1.0 Z_⊙) of the model stars, as well as the escape fraction (0.05 ≤ f_esc ≤ 0.7) of the ionizing radiation escaping from the H II regions around these stars.

The range of line ratios predicted by the grid of models encompasses the gamut of observed line ratios and thus supports the idea that the DIG is for the most part ionized by (filtered) radiation of hot stars in the galactic plane. As expected, we find that the attenuation accounts for two – partly opposing – effects that result in differences between the emission line spectra of the H II regions and of the DIG. The first effect is radiation hardening: owing to the preferential absorption of photons near the ionization thresholds, the ionizing photons that escape from the H II regions are on average more energetic than the photons absorbed in the environment of the sources. Thus, photoionization heating leads to higher temperatures in the DIG than in H II regions, which in turn results in stronger emission of collisionally excited optical lines such as N II, O II, O III, or S II.

Second, for stars with effective temperatures T_eff ≲ 35 000 K the volume in which helium is ionized is smaller than the corresponding volume of hydrogen, so that an H II region ionized by such stars (or by a group of stars whose SED is dominated by such stars) can be matter-bounded with respect to hydrogen, but radiation-bounded with respect to helium. In this case the escape fraction of helium-ionizing photons is significantly lower than the escape fraction of hydrogen-ionizing photons, and therefore the line emission from ions that require photon energies above the ionization threshold of He I to be formed – such as O III – is reduced or nonexistent in gas ionized by the filtered radiation.

This allows our model DIG to show similar properties as the DIG in the Milky Way, whose observed spectra indicate that both nitrogen and oxygen are predominantly singly ionized and helium appears to be mostly neutral. Nevertheless, stellar radiation with a low number of He-ionizing photons as the sole source of ionization would not be able to account for an increasing [O III]∕Hβ ratio at large distances from the galactic plane as has been observed in the DIG of some edge-on galaxies.

The DIG itself also acts as a filter for the regions farther away from the ionizing sources. To study the effects of different distributions of inhomogeneities in the DIG, we performed a series of 3D simulations representing a section of the galactic diskaround a spiral arm and the diffuse gas above and below it. The models share a common exponential density distribution of the gas but differ in the degree and type of inhomogeneities. We consider three different scenarios: homogeneous gas, inhomogeneous gas with microclumping (i.e., the inhomogeneities are optically thin and not spatially resolved within the simulation), and inhomogeneous gas with macroclumping, for which we have chosen the inhomogeneities to be of similar scale as the resolution element of the simulation of approximately 50 pc. For the ionizing sources we used an SED computed using the stellar population synthesis code Starburst99 corresponding to the star formation rate in a spiral arm and a Kroupa IMF. In our simulation, the ionizing flux is equally distributed among a set of sources that are positioned such that their spatial distribution roughly corresponds to those in a spiral arm, and additionally modified by homogeneous H II regions with an escape fraction of 0.1.

We find that small-scale inhomogeneities – as in microclumping, or in those macroclumping cases in which the average distance between clumps is smaller than the mean free path of photons in a homogeneous medium – considerably reduce the ability of ionizing photons to escape from the environment of the galactic disk. However, macroclumping with high clumping factors leadsto the formation of channels through which the radiation is able escape from the environment of the galactic disk without being absorbed.

We present synthetic line ratio maps that not only show the run of the line ratios as functions of galactic height but also illustrate the spatial variability of these ratios resulting from the nonregular distribution of the ionizing sources (which obviously depends on their chosen arrangement). We find that [O II]∕Hβ and [N II]∕Hβ increase for larger heights above the galactic plane. However, we find considerably larger [O III]∕Hβ ratios than observed in the Milky Way or in M 31, which may indicate that although the overall number of photons emitted by the stellar population in a star-forming galaxy would be sufficient to ionize the DIG, the filtering characteristics by the stellar environment of the hot massive stars may be different from those we assumed in this simulation.

For comparison purposes and because previous results from other work along with our spherically symmetric results indicate ionization by a softer radiation field, we additionally performed a 3D simulation in which we changed the SED of the ionizing radiation to correspond to the unfiltered spectrum of a 35 kK model star. As expected from our spherically symmetric results, we found a considerably better agreement with the observed line ratios.

To broaden the scope of simulations of the DIG, future work may combine 3D simulations with parameter studies regarding the properties of the stellar populations (metallicity, IMF, and time-dependent effects like starbursts) and consider possible additional energy sources such as magnetic reconnection or hot evolved stars outside the galactic disk. It may also account for the distribution of the ionizing sources in individual galaxies.

Acknowledgements

We would like to thank R. Walterbos and the Astrophysical Journal for granting us permission to reproduce the observational images used in Fig. 3. Furthermore, we would like to thank our referees for their useful suggestions that helped to improve this work. 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

Fig. 1 Computed stellar SEDs (black lines) modified by homogeneous matter-bounded H II regions (nH = 10 cm-3) for different escape fractions (colored lines as labeled). For the sources we have chosen the main-sequence spectra from the hot star model grid by Weber et al. (2015). The chemical composition of the H II regions is assumed to be the same as that of the corresponding stellar atmospheres. Shown are the Eddington fluxes normalized to the stellar radius.

In the text

Fig. 2 Energy-dependent escape fractions F(ν)∕F⋆(ν) of photons above the ionization threshold of hydrogen for the fluxes shown in Fig. 1. In general, the escape fraction rises above each of the strong ionization edges due to the decreasing absorption cross sections for higher photon energies. The He II-ionizing photons are almost completely absorbed within the H II regions for all considered stellar sources; for Teff = 30 000 K almost no He I-ionizing radiation can escape from the H II region either. For higher metallicities, in addition to the prominent ionization edges of H I, He I, and He II, there is also notable absorption by metal ions, especially of O II, Ne II, N III, and C III with ionization energies of 2.59, 3.35, 3.49, and 3.52 Ryd, respectively. The dashed lines indicate the corresponding angle-averaged escape fractions for the case of extremely clumpy H II regions that would allow radiation to escape unhindered in some directions while being completely opaque to all wavelengths in others (“picket-fence” case). The spike in the escape fractions of the upper right panel is a result of the extremely low stellar He II-ionizing photon flux, which is approximately 12 orders of magnitude lower than the corresponding hydrogen-ionizing flux. Thus, even very little emission of photons by the H II region in that energy range is sufficient to exceed the flux emitted by the stellar atmosphere.

In the text

Fig. 3 Comparison between Hα observations of H II regions and the DIG in the face-on galaxy M 81 by Greenawalt et al. (1998; middle row) and our simulation results (bottom row) of the ionization of the DIG by radiation escaping from the H II regions in the environment of the radiation sources. We find a reasonable qualitative agreement regarding both the denseand relatively intense H II regions around the star-forming regions on the left-hand side of the imageand the DIG on the right-hand side of the image. The simulation results presented in this image were computed for a combination of microclumping and macroclumping with a clumping factor of fcl = 5 each. The brightness values in these images as well as those in Fig. 4 represent an emission measure range between 1 and 750 pc cm-6 (logarithmicscaling). Observational images: ©AAS. Reproduced with permission.

In the text

Fig. 4 Edge-on and face-on emission measure maps for the different simulated clumping scenarios (see text). In each case we assumed a mean hydrogen number density of nH (h) = 0.1 cm−3 × e−|h|∕h0, where |h| is the distance from the galactic plane and h0 = 1 kpc is the scale height of the gas. The simulations also share the same distribution and luminosities of the ionizing sources. The edge lengths of the maps are 5 kpc. The macroclumping case with fcl = 5 appears generally very similar to the microclumping case with the same clumping factor, as is expected for sufficiently fine grid resolution. For fcl = 25, the macroclumping becomes very porous at the used grid resolution, opening up channels devoid of gas that allow radiation to escape from the immediate environment of the sources without being absorbed. The combined micro- and macroclumping case with fcl = 5 × 5 is very different from the macroclumping case with fcl = 25, giving an overall appearance resembling that of the cases with fcl = 5, but with theemission in each of the ionized bubbles more concentrated; this is somewhat analogous again to the transition from fcl = 1 to fcl = 5.

In the text

	Fig. 5 Left panel: SEDs from the population synthesis model, filtered as the stellar atmosphere models in Fig. 1. For our 3D simulations we use the SED of the model with fesc = 10% (magenta) as ionizing flux for the DIG. Right panel: corresponding energy-dependent escape fractions, as in Fig. 2. A comparison with Figs. 1 and 2 shows that the SEDs roughly resemble those using a 40 kK star.
In the text

Fig. 6 Left panel: ionization fractions of hydrogen as a function of the height above the galactic plane for the different clumping scenarios presented in Fig. 4. For the homogeneous density distribution, the gas is predominantly ionized throughout the simulation volume and the ionization fraction rises with increasing galactic heights due to the lower gas densities, which results in decreasing recombination rates and opacities. Both the micro- and the macroclumping case with fcl = 5 are characterized by a similar dependence on the height above the galactic plane, reaching approximately 0.5 near the galactic plane and dropping to 0.1 at a height of about 1.2 kpc. In contrast, in the macroclumping scenario with fcl = 25 the ionization fraction varies only slightly throughout the simulation volume. The combination of micro- and macroclumping shows the lowest ionization fraction of all the scenarios owing to the increased recombination rates in the filled cells. While the total clumping factor equals that of the fcl = 25 macroclumping case, there are fewer extended empty regions through which radiation can penetrate to greater distances from the galactic plane, where less ionizing radiation is required to maintain the same fraction ofionized hydrogen due to the lower gas densities. Right panel: as before, but for a parameter study varying the microclumping factor between 1 and 5. The ionization fraction rises at larger galactic heights inthe case of low clumping, whereas it decreases and the gas eventually becomes neutral for high clumping factors. For immediate clumping factors, which for our geometry and density are around fcl ~ 3, the ionization fraction as a function of height behaves in a nonmonotonic way, first dropping to a minimum and then rising again.

In the text

	Fig. 7 Line ratio maps of several collisionally excited lines relative to Hβ in the face-on view of the simulated spiral arm for different density structures of the gas. Gray indicates emission measures less than 1 pc cm-6.
In the text

Fig. 8 As Fig. 7, but an edge-on view. Oxygen and nitrogen are predominantly doubly ionized close to the ionizing sources, and singly ionized only in the outer part of the ionized region. Thus, and because of radiation hardening and therefore higher temperatures, the maximum of the line ratios [N II]∕H β and [O II]∕Hβ occurs at larger heights above the galactic plane. For the cases with a clumping factor of fcl = 5 the [O III]∕Hβ ratio rises up to a height of approximately 700 pc above the galactic plane and drops again at larger distances.

In the text

Fig. 9 Line ratio maps from microclumped models (fcl = 5) as in the leftmost columns of Figs. 7 and 8, using the same simulation volume and the same integrated hydrogen-ionizing flux as before, but taking the unfiltered spectrum of a 35 kK dwarf star with solar metallicity as the ionizing SED. The most notable difference to the models using the population SED is the weaker emission of the optical [O III] lines. The emission of the [N II] and [O II] lines is stronger close to the galactic plane, but shows a weaker rise close to the edge of the simulated volume.

In the text