© ESO, 2013

1. Introduction

Electrons and protons combined to form neutral hydrogen when the young universe cooled (at a redshift of z ≈ 1100), but the continuous absorption troughs blueward of the (redshifted) hydrogen Lyα line, that one would expect to see if the intergalactic medium (IGM) consisted of neutral (recombined) hydrogen gas, were not found in the spectra of distant objects around z ≈ 2.0 (Gunn & Peterson 1965). It turns out that one has to look much deeper to find these troughs, and by means of quasar observations, White et al. (2003) and Kashikawa (2007) were able to show that the troughs do exist for objects at a redshift of z ≈ 6.5, but not for objects at a redshift of z ≈ 5.7. Obviously, a process has occurred that re-ionized the intergalactic gas in the universe, and this process must have been completed by a redshift of roughly z ≈ 6¹. As a result of this process the number density of neutral hydrogen is now very low in the IGM. At which redshifts it was still high is, due to the lack of corresponding observations, more difficult to determine, and one currently has to rely on numerical cosmological simulations to estimate when exactly reionization began. The simulations of Wise & Abel (2008) point for instance to a redshift of roughly z ≈ 30 as starting point of the reionization process and marking the end of the “cosmic dark ages”².

Regarding the nature of the sources that have driven reionization we are still left somewhat in the dark. Although it was immediately recognized that quasars, as the brightest single objects in the UV, should have played an important role in this process, a quantitative investigation of the luminosity function of high-redshift quasars showed that their radiation output was not sufficient to reionize the universe (Fan et al. 2006a)³. The other likely candidate objects to have significantly driven the reionization process are stars, the first generation of which (so-called population III stars) are assumed to have developed in the dark matter halos that had already formed at those redshifts.

There is reason to believe that due to their primordial composition and the corresponding inefficient cooling of the gas during the star-formation process these population III stars formed on a larger fragmentation scale (10²...10³ M_⊙) without subfragmenting further (cf. Abel et al. 2000; Bromm et al. 2002), and therefore the initial mass function (IMF) of the primordial stars differs significantly from the IMF at present days with a preference to massive stars (Bromm et al. 2002 give masses of 100 ≤ M/ M_⊙ ≤ 1000, and temperatures of T_eff ~ 10⁵ K, which to a large extent turn out to be independent of the stellar masses for these short living (τ_MS ~ 2 × 10⁶ yr, cf. Schaerer 2002 objects). Although the exact form of this top-heavy IMF is still under debate⁴, the preference for massive stars leads to huge UV luminosities, and this energy – released at high frequencies – may have been crucial for the evolution of the universe in its early stages.

Thus, the existence and the properties of the first stellar population is closely linked to the reionization history of the universe. But the time period for the appearance of these objects was much shorter than the time required to reionize the IGM, and different stellar populations will have contributed to the reionization at different times. In this context the history of reionization may be divided into three stages. The first one is the “pre-overlap” phase which was characterized by isolated H ii regions located in the interstellar medium (ISM) which expanded steadily into the low density IGM. The second one is the “overlap” phase in which the separate H ii regions started to merge; this process required a strong rise in the amount of ionizing radiation. The last stage is the “post-overlap” phase, characterized by ionizing the low density IGM on large scales (Gnedin 2000).

A more detailed description of this history has been presented by Cen (2003) whose analysis indicated that the universe was reionized in two steps. The first one occurred up to a redshift of z ~ 15 where the ionization was primarily driven by population III stars. As the population III stars returned metals via supernova explosions into their environment, the chemical composition and – as a result – the IMF and spectral energy distributions (SEDs) of later populations of stars changed, giving way to metal-enriched population II stars⁵. This transition from population III to population II stars lead to a decline in the emitted radiation power by a factor of ~10, which in turn lead to a phase, from about z ≈ 13 to z ≈ 6, during which hydrogen partly recombined (n_H ii/n_H ≥ 0.6) (Cen 2003). Thus the universe may have come close to a second cosmological recombination at z ~ 13. To overcome the threatening re-recombination and to restart the reionization process driven by powerful ionizing sources and maintained by the dilution of the intergalactic gas due to the expansion of the universe, a second step (finishing at z ~ 6) was necessary⁶. But how could the presumably much less massive population II and population I stars fully reionize the hydrogen and helium content of the IGM? Of course, there is strong evidence of an increasing star formation rate (SFR) in this epoch (Lineweaver 2001; Barger et al. 2000), but in particular regarding the reionization of helium not only the number of stars is important, but also the characteristics of the SEDs of the stars are much more significant (Wyithe & Loeb 2003). Very massive stars (VMSs) which clearly exceed the assumed upper limit of 100 to 150 M_⊙ of present-day massive stars have pronounced fluxes in the EUV at all metallicities (Pauldrach et al. 2012) and in conjunction with a top-heavy stellar mass distribution may provide such a required SED.

To motivate this idea we note that there are indications of a more top-heavy IMF than the standard one even in present-day galaxies. Harayama et al. (2008) for instance studied the IMF of one of the most massive Galactic star-forming regions – NGC 3603 – to verify whether the IMF really does have a universal shape. What they found does not support this assumption: with Γ = − 0.74 they derived a power-law index for this massive starburst cluster which is considerably less steep than the Salpeter IMF (Γ = − 1.35). Their result thus supports the hypothesis that a top-heavy IMF is not unusual for massive star-forming clusters and starburst galaxies even at solar metallicity. They further argued that a common property among starbursts showing such flat IMFs would be the high stellar density in the core of the cluster, and that the variations of the IMF could be linked to the spatial density of the stellar population (according to their results the core mass density of NGC 3603 is at least 6 × 10⁴ M_⊙ pc^-3).

Moreover, very massive stars are efficient emitters of ionizing photons with approximately 10 times more hydrogen-ionizing photons and 10 000 times more He ii-ionizing photons per unit stellar mass than a population with a Salpeter IMF. Pauldrach et al. (2012) discuss a mechanism to form such objects in present-day stellar clusters. They combined consistent models of expanding atmospheres with stellar evolutionary calculations of massive single stars with regard to the evolution of dense stellar clusters, and investigated the conditions necessary to initiate a runaway collision merger which may lead to the formation of a very massive object (up to several 1000 M_⊙) at the center of the cluster, and the possible formation of intermediate-mass black holes (IMBHs)⁷. Interestingly these theoretical models have been observationally supported by two studies which demonstrate the existence of VMSs in the local Universe. The first one was based on HST and VLT spectroscopy from which Crowther et al. (2010) concluded that the dense cluster R136 located in the 30 Doradus region of the Large Magellanic Cloud hosts several stars whose initial masses were up to ~ 300 M_⊙. The second study concerned the discovery of an optical transient which was classified as a Type Ic supernova (SNF20070406-008). As the observed light curve of this object fits that of a pair-instability supernova with a helium core mass of at least 100 M_⊙, the progenitor of this supernova must have been a VMS (Gal-Yam et al. 2009). Runaway collision mergers clearly exceed the assumed upper limit of 100...150 M_⊙ for the direct formation of present-day massive stars, and since such mergers can also occur in chemically evolved clusters of high core mass density, they may have played an important role at least in the late stages of the reionization history.

This conclusion is relevant in particular for the reionization of He ii to He iii, since this process was not yet completed at a redshift of z ~ 6. (Spectroscopic observations of He ii Lyman α absorption troughs of quasars can be observed even at z = 2.8; Reimers et al. 1997; Kriss et al. 2001; Syphers et al. 2011.) The reionization of He ii was therefore considerably delayed compared to the reionization of H i and He i. Whether the appearance of population II and population I VMSs could have been responsible for the reionization of He ii has not yet been investigated; it is presently an open point whether these very special stars, or, as assumed by Wyithe & Loeb (2003), Gleser et al. (2005), and McQuinn et al. (2009), quasars, or a mixture of these objects have been responsible for the reionization of He ii.

To simulate and understand the reionization process, different, numerical approaches have been developed for a description of the radiative transfer and the evolution of the ionization fronts (e.g. Mellema et al. 2006; Gnedin & Abel 2001; Ciardi et al. 2001; Nakamoto et al. 2001; Razoumov & Cardall 2005; Ritzerveld et al. 2003; Alvarez et al. 2006; Reynolds et al. 2009; Iliev et al. 2009; Trac & Cen 2007; McQuinn et al. 2007; Wise & Abel 2011). Many of these codes are specialized for a particular task and do not attempt to provide a comprehensive description of the time-dependent ionization structure, including the detailed statistical “equilibrium” of all relevant elements. E.g., for the purpose of simulating the radiative transfer in the early universe it is usually sufficient to just include tabulated ionization and recombination rates for hydrogen and helium, which are modeled as simple two- or three-level systems, to neglect metals, and to assume simplified models for SED of the sources – e.g., black-body radiators or power-law spectra described by a few discrete wavelengths specifying the fluxes in the ranges important for the particular simulations.

Our focus in contrast lies in a sophisticated description of the 3d radiative transfer with respect to high spectral resolution, in order to quantify the evolution of all relevant ionization structures accurately, in particular also with the intent of determining possible observational features, such as line strength ratios from different ions. Metal lines, for instance, could be used to reconstruct the history of metal enrichment and the characteristics of the spectra of the ionizing sources (see, for instance, Graziani et al. 2013) that might be used to discriminate between different scenarios and/or particular aspects of these scenarios. Although only H and He are treated in the present paper, our algorithm in principle allows the implementation of multiple levels for each ionization stage of the metals, and we are currently working on the inclusion of metals as described by Hoffmann et al. (2012). We note in this regard that while metal ions have no direct significance as absorbers affecting the radiation field and thus the ionization of hydrogen and helium, metal cooling has a significant influence on the temperature structure of the gas and therefore on the temperature-dependent recombination rates of H and He. Thus, the metallicity dependent temperature has an impact on the ionization structure of the most abundant elements (e.g., Osterbrock & Ferland 2006). But this procedure only makes sense if realistic SEDs of the sources driving the reionization of the universe are supplied – from population III to population I stars at various metallicities, temperatures, and masses. To this aim we use a sophisticated model atmosphere code based on a consistent treatment of the expanding atmospheres of massive and very massive stars (cf. Pauldrach et al. 2001, 2012; note that the winds of hot stars modify the SEDs of the ionizing radiation of the stars dramatically, cf. Pauldrach 1987; Pauldrach et al. 1994). We are therefore able to calculate the required SEDs to be used as input for the 3d radiative transfer simulating the ionization structure of the ISM and IGM.

In the following we will first introduce the theoretical basis of our simulations covering the concept of 3d radiative transfer and the physical mechanisms which influence the ionization structure of the gas surrounding the sources of ionization. We will further present a numerical method to describe the temporal expansion of the ionization fronts and which considers multiple ionization sources (Sect. 2). In Sect. 3 we discuss important details of our numerical approach and present tests comparing our 3d results with those of analytical solutions, a radially symmetric radiative transfer code, and results from a comparison project initiated by Iliev et al. (2006). In Sect. 4 we present first applications of our 3d radiative transfer on showcase simulations of the reionization scenario. Although using simplified initial conditions, we will demonstrate the influence of realistic stellar spectra of massive stars at different metallicities on the ionization structure of the surrounding gas under conditions which correspond to those of the early universe, and we will focus on the different behavior of the expansion of the ionization fronts with respect to a homogeneous gas density and an inhomogeneous cosmological density structure. With regard to the He ii reionization problem we further investigate the influence of different input spectra on the ionization structure of He iii by performing a series of multisource simulations. We interpret and summarize our results along with an outlook in Sect. 5.

2. Three-dimensional radiative transfer as a tool for the description of the epochs of reionization

On larger scales radiation is emitted from point sources and this radiation starts to propagate homogeneously and isotropically in all directions up to a certain distance where material of the environment of the sources is encountered in the form of a density structure, which may also involve inhomogeneities or fluctuations, and radiative processes start to have a disturbing influence on this well-suited behavior. Moreover, the radiation propagating from the different sources will overlap at specific locations and this produces a time and frequency dependent 3-dimensional pattern of radiation. The evolution of such a radiation field and all its variables – intensity, radiative flux, photon flux – is however well defined and at every point in space characterized by the equation of radiative transfer: $\begin{array}{ccc}\left[\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{-}\frac{\mathit{H\nu }}{\mathit{c}}\frac{\mathit{\partial }}{\mathit{\partial \nu }}\mathrm{+}{n̂}\mathrm{·}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\right]{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{r}\mathit{,}{n̂}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}_{\mathit{\nu }}\mathrm{-}{\mathit{\chi }}_{\mathit{\nu }}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{r}\mathit{,}{n̂}\mathit{,t}\mathrm{\right)}\mathit{.}& & \end{array}$(1)This Boltzmann equation (cf. Gnedin & Ostriker 1997) describes the transport of radiative energy via a change of the intensity I_ν caused by the absorption coefficient ${\mathit{\chi }}_{\mathit{\nu }}\mathrm{\equiv }{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}{r}\mathit{,}{n̂}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{\left(}{r}\mathit{,}{n̂}\mathit{,t}\mathrm{\right)}\mathrm{+}\frac{\mathrm{3}\mathit{H}}{\mathit{c}}$ and the emissivity as a function of position r (in comoving coordinates), direction , and time t for every frequency ν. In this equation, c is the speed of light, a ≡ a(t) is the cosmological scale factor, H ≡ H(t) = ȧ(t)/a(t) is the Hubble expansion rate, and κ_ν is the opacity; as κ_ν/(3H/c) is on the order of 10⁴...10⁶ in the relevant frequency range of a neutral gas with a density equal to the mean density of the universe in the considered redshift range (0 ≤ z ≤ 15), we assume χ_ν = κ_ν in the remaining part of this paper.

Neglecting the explicit time derivative term (for a rationale see Sect. 3.2) and the frequency derivative term (justified if the universe does not expand significantly before the corresponding photons are absorbed – cf. Petkova & Springel 2009; Wise & Abel 2011) in Eq. (1), and making use of the identity ${n̂}\mathrm{·}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }\mathrm{=}\frac{\mathrm{d}}{\mathrm{d}\mathit{s}}$, which involves the interpretation of ${n̂}\mathrm{·}\frac{\mathrm{1}}{\mathit{a}}\mathrm{\nabla }$ as a directional derivative along a path element s, the equation of radiative transfer becomes for each ray to be considered $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{\mathrm{d}\mathit{s}}\mathrm{=}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{\left(}{\mathit{S}}_{\mathit{\nu }}\mathrm{\right(}\mathit{s}\mathrm{)}\mathrm{-}{\mathit{I}}_{\mathit{\nu }}\mathrm{(}\mathit{s}\mathrm{\left)}\mathrm{\right)}\mathit{,}& & \end{array}$(2)(where, as usual, S_ν represents the source function, defined as the ratio of emissivity to opacity, S_ν = η_ν/χ_ν). Equation (2) can now be solved formally in an analytical way $\begin{array}{ccc}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\hspace{0.17em}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{+}{\mathrm{\int }}_{{\mathit{s}}_{\mathrm{0}}\mathrm{\left(}\mathit{n}\mathrm{\right)}}^{\mathit{s}}{\mathit{\eta }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\hspace{0.17em}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(3)s₀(n) corresponds to the starting point of each ray n, and ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{{\mathit{s}}_{\mathrm{0}}\mathrm{\left(}\mathit{n}\mathrm{\right)}}^{\mathit{s}}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}$ is the usual dimensionless optical depth variable quantifying the absorption characteristics of the medium⁸. Numerically the equation is solved using a suitable discretization scheme dividing the volume into small cells. Since by definition of this discretization scheme the opacity χ_ν(s) stays constant along the cell-crossing distance l_n(m) of each specific cell m crossed by the considered ray n, the corresponding optical depth τ can be expressed as $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\sum _{{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}\mathit{m}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}}^{\mathit{m}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{\mathit{\chi }}_{\mathit{\nu }}\left({\mathit{m}}^{\mathrm{\prime }}\right){\mathit{l}}_{\mathit{n}}\left({\mathit{m}}^{\mathrm{\prime }}\right)& & \end{array}$(4)(where m(s₀(n)) denotes the cell where the corresponding source of the considered ray n is placed and m(s) denotes the cell which has just been crossed before s is reached) and Eq. (3) thus simplifies to $\begin{array}{ccc}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}}\mathit{.}& & \end{array}$(5)It is now straightforward to calculate from Eq. (5) the energy deposited per time in the cells that are crossed by a ray n over a distance s_n(m) $\begin{array}{ccc}\mathrm{\Delta }\mathit{Ė}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}}\mathrm{)}\mathrm{d}\mathit{\nu }\mathit{.}& & \end{array}$(6)Here $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{=}{\mathit{L}}_{\mathit{\nu }}\mathit{/}{\mathit{N}}_{\mathrm{rays}}$ is the fraction of the (spectral) luminosity emitted by the source into each of the N_rays rays considered per source. (All rays start at sources (cf. Fig. 1), which represent single stars or star clusters, each with a specified SED F_ν(M_n,Z_n,T_eff(n),R_n) and spectral luminosity ${\mathit{L}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathit{n}}^{\mathrm{2}}\hspace{0.17em}{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathit{n}\mathrm{\right)}$, characterized by the source’s mass M_n, metallicity Z_n, effective temperature T_eff(n), and radius R_n (for single stars, the stellar radius; for clusters, an equivalent radius describing the total radiating surface) – cf. Sect. 3.2.)

2.1. 3-dimensional radiative transfer based on ray tracing

Because the equation of radiative transfer has to be solved explicitly for every time step at every point in space, a suitable discretization scheme is required that must be able to deal with possible inhomogeneities in the density distribution of the gas in the considered volume and/or the presence of a variable number of different radiative sources which in general will not possess any kind of symmetry (cf. Figs. 1 and 2). The natural property of light to propagate along straight lines straightforwardly leads to the concept of a ray by ray solution which involves distributing rays isotropically around each source and solving the transfer equation for each of these rays, taking into account the interaction of the gas and the photons along the way⁹.

The ray tracing concept along with a spatial discretization scheme described by a Cartesian coordinate system can naturally account for the following three important points:

• a radiation field that is generated by numerous different pointsources arbitrarily distributed in space and characterized byindividually different SEDs,

• an inhomogeneous density distribution of the medium,

• the temporal evolution of the radiation field and the propagation of the ionization fronts.

Although in our discretization scheme the physical conditions are treated as being spatially constant within each cell, and we consider (for reasons of efficiency) each of the emitting sources to be located in the center of its respective cell, we can account for arbitrary distributions of the gas and the sources with a resolution limited only by the computational resources. Time-dependencies can in principle be described accurately with a suitably fine resolution in space and time. With regard to the spatial resolution we note that it is also essential that the radiation field of each source is itself discretized with a matching angular resolution so that each cell is crossed by a sufficient number of rays in order to describe the effects of the radiation on the gas correctly (cf. Fig. 3 and Sect. 3.1).

2.2. Physical processes affecting the state of ionization

The basis of any approach in constructing detailed radiative models for the ionization structure of the gas surrounding clusters of sources of ionization is a concept that includes the time-dependent statistical equilibrium for all important ions with detailed atomic physics – described by rate equations –, the energy equation, and the radiative transfer equation at all transition frequencies required. As all of the involved equations have to be solved simultaneously, and as the replication of the required physical processes, which due to the energy input of the time varying sources continuously modify the physical properties of the surrounding gas, makes the solution of every realistic approach to a formidable problem, the method is not simple at all. We will in the following therefore give a description of the physics to be treated in some detail.

2.2.1. The rate equations

As an essential step of the procedure the time-dependent occupation numbers n ≡ n(r,t) of all ionization stages i, j of the elements considered have to be calculated. The time-dependent statistical equilibrium $\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{n}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{i}\ne \mathit{j}}{\mathrm{𝒫}}_{\mathit{j,i}}{\mathit{n}}_{\mathit{j}}\mathrm{-}\sum _{\mathit{i}\ne \mathit{j}}{\mathrm{𝒫}}_{\mathit{i,j}}{\mathit{n}}_{\mathit{i}}\mathit{,}& & \end{array}$(7)which describes the temporal derivation of the number density of an ionization stage i, and which contains via the rate coefficient all important radiative ℛ_i,j and collisional transition rates, serves as a basis for this step¹⁰. To solve these systems of differential equations, we define a vector n, which contains the number densities of all ionization stages of the considered element, and a tridiagonal matrix G, which contains in its components g_i,j the rate coefficients¹¹. With these definitions Eq. (7) is rewritten as $\begin{array}{ccc}{E}\mathrm{·}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{n}\mathrm{+}{G}\mathrm{·}{n}\mathrm{=}{0}\mathrm{,}& & \end{array}$(8)where for reasons which will become obvious below the unity matrix E has been applied on the time derivative term on the left hand side of the system of equations. By scaling Eq. (8) with the total number density of the element considered, and by replacing the fraction x_k of the ionization stage k by the condition of particle conservation $\begin{array}{ccc}{\mathit{x}}_{\mathit{k}}\mathrm{=}\mathrm{1}\mathrm{-}\sum _{\mathit{j,j}\ne \mathit{k}}^{\mathit{N}}{\mathit{x}}_{\mathit{j}}\mathit{,}& & \end{array}$(9)where N is the number of ionization stages and x contains in its components the number densities of the ions relative to the total number density of the element, one gets the final system of inhomogeneous equations $\begin{array}{ccc}{E}\prime \mathrm{·}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\hspace{0.17em}{x}\mathrm{+}{{G}}^{\mathrm{\prime }}\mathrm{·}{x}\mathrm{=}{b}\mathit{,}& & \end{array}$(10)where the components of G′ are given by ${\mathit{g}}_{\mathit{i,j}}^{\mathrm{\prime }}\mathrm{=}{\mathit{g}}_{\mathit{i,j}}\mathrm{-}{\mathit{g}}_{\mathit{i,k}}$, and those of b are b_i = − g_i,k, and where all coefficients of the redundant kth row of G′ and b have been replaced by 1, and those of the kth column of E by 0 (with this replacement the unity matrix E becomes E′) – with these numbers inserted the corresponding components represent in total the condition of particle conservation.

Modeling the temporal evolution of the ionization structure.

The solution of the time-dependent rate equations (Eq. (10)) required for modeling the temporal evolution of the ionization structures is not straightforward, because this system forms a set of stiff differential equations, and this means that the occupation numbers can change by several orders of magnitude for a gradual increase of the timescale¹². On the other hand, for constant coefficients ¹³ the rate equations form a set of linear differential equations, and such a system can in principle be solved by an eigenvalue approach. With respect to this approach the general solution for the differential equations is based on the time independent “equilibrium” solution x^∞ of the inhomogeneous system $\begin{array}{ccc}{G}\prime \mathrm{·}{{x}}^{\mathrm{\infty }}\mathrm{=}{b}\mathit{,}& & \end{array}$(11)and the solution of the homogeneous system $\begin{array}{ccc}{E}\prime \mathrm{·}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{G}}^{\mathrm{\prime }}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{0}\mathit{,}& & \end{array}$(12)which result in the composed solution $\begin{array}{ccc}{x}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{x}}^{\mathrm{\infty }}& & \end{array}$(13)of Eq. (10).

While the last component x^∞ of the composed vector is obtained by a direct solution of Eq. (11), the more interesting homogeneous component x_hom(t), which contains in its coefficients just difference values, can not always be determined straightly; its physical behavior however results directly from the structure of Eq. (12) $\begin{array}{ccc}{{x}}_{\mathrm{hom}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{\left(}{x}\mathrm{\right(}{\mathit{t}}_{\mathit{k}}\mathrm{)}\mathrm{-}{{x}}^{\mathrm{\infty }}\mathrm{)}\mathrm{·}{\mathrm{e}}^{\mathrm{-}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathit{k}}\mathrm{)}}\mathit{,}& & \end{array}$(14)where t_k marks the begin of the considered timestep, and δ is the eigenvalue of the system¹⁴.

The numerical solution of the time-dependent behavior of the ionization structures is on basis of the eigenvalue approach stable, even for large time steps (cf. Sect. 3.3). Nevertheless is an appropriate choice of the timestep sizes extremely important, since the evolution of the ionization structures has to be traced with a sufficient temporal accuracy. In the frame of our approach it turned out that the timestep-size can be controled quite well by the following extrapolation method, which limits the maximal relative change of the ionization state within one timestep, $\begin{array}{ccc}\mathrm{\Delta }{\mathit{t}}_{\mathit{k}}\mathrm{=}\mathrm{\Delta }{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\underset{\mathrm{all}\mathrm{cells}\mathit{m}}{\min }\left(\mathit{f}\left|\frac{{\mathit{n}}_{\mathrm{H}\hspace{0.17em}{i}}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}}{{\mathit{n}}_{\mathrm{H}\hspace{0.17em}{i}}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}\mathrm{-}{\mathit{n}}_{\mathrm{H}\hspace{0.17em}{i}}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{2}}\mathrm{\right)}}\right|\right)\mathit{,}& & \end{array}$(15)where Δt_k − 1 is the duration of the previous timestep, and n_H i(m,t_k − 2) and n_H i(m,t_k − 1) are the number densities of neutral hydrogen at the beginning of the previous and the current time steps, respectively¹⁵.

An adequate description of the evolution of the ionization structures furthermore requires the consideration of the cosmological expansion for every timestep via a scaling of the path elements Δs_n(m,t) described by ${n̂}\mathrm{·}\frac{\mathrm{1}}{\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{\nabla }\mathrm{=}\frac{\mathrm{d}}{\mathrm{d}\mathit{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}}$ (cf. Eq. (1)), i.e., the correct cosmological scale factor a(t) has to be known for every timestep of the simulation. For the redshift range of our simulations the scale factor can be approximated quite well by the formula $\begin{array}{ccc}\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}\mathit{,}\mathrm{0}}}{{\mathrm{\Omega }}_{\mathrm{\Lambda }\mathit{,}\mathrm{0}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\sinh }^{\mathrm{2}\mathit{/}\mathrm{3}}\left(\frac{\mathrm{3}}{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}\hspace{0.17em}\mathit{t}\sqrt{{\mathrm{\Omega }}_{\mathrm{\Lambda }\mathit{,}\mathrm{0}}}\right)\mathit{,}& & \end{array}$(16)where H₀ is the Hubble constant at z = 0, and Ω_m,0 and Ω_Λ,0 are the usual relative contributions of matter and dark energy to the total energy of the present universe (cf. Carroll & Ostlie 2006). This is used to recompute all path elements Δs_n(m,t) for each timestep t_k according to $\begin{array}{ccc}\mathrm{\Delta }{\mathit{s}}_{\mathit{n}}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\mathrm{\Delta }{\mathit{s}}_{\mathit{n}}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}\frac{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}}{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}}\mathrm{\forall }\mathit{n}\mathrm{\forall }\mathit{m}\mathit{.}& & \end{array}$(17)The particle densities n(m,t_k) giving rise to the emission and absorption coefficients in the radiative transfer equation are expressed in units of the volumes ΔV(m,t_k) of the grid cells, which scale with the cube of the same factor a(t_k)/a(t_k − 1), have to be rescaled analogously via $\begin{array}{ccc}\mathit{n}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\mathit{n}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}\frac{\mathrm{\Delta }\mathit{V}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}}{\mathrm{\Delta }\mathit{V}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}}\mathrm{=}\mathit{n}\mathrm{\left(}\mathit{m,}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}{\left(\frac{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right)}}{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}}\right)}^{\mathrm{3}}\mathrm{\forall }\mathit{m}\mathit{.}& & \end{array}$(18)

2.2.2. The rate coefficients

In the following paragraphs we describe briefly the ionization and recombination rates we used as coefficients for the rate matrix G. Under nebular conditions, only ionization from the ground state will be important, whereas recombination may occur also to excited levels¹⁶. Therefore the total rates between two ionization stages are given by the sums of the rate coefficients of the specific atomic energy levels l – ground state or excited level of the lower ionization stage i – and u – ground state or excited level of the upper ionization stage i + 1 $\begin{array}{ccc}{\mathit{n}}_{\mathit{i}}{\mathrm{𝒫}}_{\mathit{i,j}}\mathrm{=}\begin{array}{c}\sum \\ \mathit{l}\mathrm{\in }\mathit{i}\mathit{u}\mathrm{\in }\mathit{i}\mathrm{+}\mathrm{1}\end{array}{\mathit{n}}_{\mathit{l}}{\mathit{P}}_{\mathrm{lu}}\mathit{,}{\mathit{n}}_{\mathit{j}}{\mathrm{𝒫}}_{\mathit{j,i}}\mathrm{=}\begin{array}{c}\sum \\ \mathit{u}\mathrm{\in }\mathit{i}\mathrm{+}\mathrm{1}\mathit{l}\mathrm{\in }\mathit{i}\end{array}{\mathit{n}}_{\mathit{u}}{\mathit{P}}_{\mathrm{ul}}\mathit{,}& & \end{array}$(19)where P_lu and P_ul are the rate coefficients connecting these levels, and which contain all important radiative (R_ul, R_lu) and collisional (C_ul, C_lu) transition rates.

Photoionization and computation of the mean intensity.

The photoionization rates are obtained by performing the frequency integral over the mean intensity J_ν weighted by the frequency-dependent ionization cross-section α_lu(ν) $\begin{array}{ccc}{\mathit{R}}_{\mathrm{lu}}\mathrm{=}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{lu}}}^{\mathrm{\infty }}\frac{\mathrm{4}\mathit{\pi }{\mathit{\alpha }}_{\mathrm{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathit{h\nu }}{\mathit{J}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu ,}& & \end{array}$(20)where ν_lu denotes the threshold frequency for the corresponding ionization edge. The computation of the ionization structures thus requires the calculation of the mean intensities J_ν for every grid cell. As these physical quantities cannot simply be calculated as numerical integrals of the radiative intensities over the entire sphere – , here dω represents the solid angle –, as is the case for radiative transfer models which are based on plane-parallel or spherically symmetric geometries, a recipe for the evaluation of the J_ν values, which takes the discrete nature of the rays into account, must be provided. This recipe requires on the one hand knowledge of the function value Ṅ_ν, which describes the number of photons transported per time and frequency by a ray n to the edge of a considered cell m positioned at a distance s_n(m) from the starting point of the ray (cf. Sect. 2) $\begin{array}{ccc}\mathit{Ṅ}\mathit{\nu }\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}}}{\mathit{h\nu }}& & \end{array}$(21)(note that this equation results directly from the considerations of Eqs. (4) and (6)), and the function value ΔṄ_ν(m,n) of the number of photons absorbed in the considered cell from ray n via the ray segment l_n(m) is analogously given by $\begin{array}{ccc}\mathrm{\Delta }\mathit{Ṅ}\mathit{\nu }\mathrm{\left(}\mathit{m,n}\mathrm{\right)}\mathrm{=}\mathit{Ṅ}\mathit{\nu }\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}{\mathit{l}}_{\mathit{n}}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{)}\mathit{.}& & \end{array}$(22)Summing over all rays crossing the considered cell by ray segments finally gives the function value ΔṄ_ν(m) of the total number of photons absorbed in the cell $\begin{array}{ccc}\mathrm{\Delta }\mathit{Ṅ}\mathit{\nu }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\sum _{{{}_{\mathrm{segments}\hspace{0.17em}\mathit{n}\mathrm{\left(}\mathit{m}\mathrm{\right)}}^{\mathrm{ray}}}^{}}\mathit{Ṅ}\mathit{\nu }\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}{\mathit{l}}_{\mathit{n}}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{)}\mathit{.}& & \end{array}$(23)The recipe requires on the other hand also knowledge about the total number of radiative processes occurring per time and frequency unit in a cell. This number, the number of the sum of absorbing particles in a cell – given by the product of the number densities of potential absorbers and their cross-sections, which defines as a sum the opacity χ_ν(m), times the volume V of the cell times the number of available photons per time and frequency unit in the cell –, is represented by the function value¹⁷$\begin{array}{ccc}\mathrm{\Delta }\mathit{Ṙ}\mathit{\nu }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{V}\frac{\mathrm{4}\mathit{\pi }{\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}}{\mathit{h\nu }}\mathrm{·}& & \end{array}$(24)The recipe is thus based on the fact that the number of radiative processes consuming photons in a cell has to be equal to the number of photons absorbed along the ray segments in that cell. That is, ΔṘ_ν(m) = ΔṄ_ν(m) and the consistent value of the mean intensity J_ν is obtained from Eqs. (23) and (24) $\begin{array}{ccc}{\mathit{J}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\sum _{{{}_{\mathrm{segments}\hspace{0.17em}\mathit{n}\mathrm{\left(}\mathit{m}\mathrm{\right)}}^{\mathrm{ray}}}^{}}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathit{\nu }}\end{array}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right(}\mathit{n}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathit{n}}\mathrm{\right(}\mathit{m}\mathrm{\left)}\mathrm{\right)}}\mathrm{[}\hspace{0.17em}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}{\mathit{l}}_{\mathit{n}}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{]}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathit{V}}\mathrm{·}& & \end{array}$(25)

Radiative recombination.

The radiative recombination rate coefficients are computed analogously to the photoionization rates $\begin{array}{ccc}{\mathit{R}}_{\mathrm{ul}}\mathrm{=}{\left(\frac{{\mathit{n}}_{\mathrm{l}}}{{\mathit{n}}_{\mathrm{u}}}\right)}^{\mathrm{\ast }}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{lu}}}^{\mathrm{\infty }}\frac{\mathrm{4}\mathit{\pi }{\mathit{\alpha }}_{\mathrm{lu}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathrm{h}\mathit{\nu }}\left(\frac{\mathrm{2}\mathit{h}{\mathit{\nu }}^{\mathrm{3}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{+}{\mathit{J}}_{\mathit{\nu }}\right){\mathrm{e}}^{\mathrm{-}\mathrm{h}\mathit{\nu }\mathit{/}\mathrm{kT}\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{d}\mathit{\nu ,}& & \end{array}$(26)where k is the Boltzmann constant, T(m) is the local value of the temperature stratification of the gas, and ${{}^{\mathrm{\right(}}\frac{{\mathit{n}}_{\mathrm{l}}}{{\mathit{n}}_{\mathrm{u}}}^{\mathrm{\left)}}}^{\mathrm{\ast }}$ is the Saha-Boltzmann factor – the ratio of the occupation numbers that would be reached in the case of thermodynamic equilibrium¹⁸.

2.2.3. The temperature structure

As the recombination rates (but also emission line intensities and other diagnostic observables) depend directly on the temperature (cf. Eq. (26)), accurate ionization structure modeling requires simultaneously determining the local temperature at every point in the gas. The physical basis for this is the microscopic energy equation which, in principle, states that the energy going into and coming out of the material must be conserved. The standard method for determining the temperature is balancing the energy gains and losses to the electron gas, including all processes that affect the electron temperature: bound-free transitions (ionization and recombination), free-free transitions, and inelastic collisions with ions. Pauldrach et al. (2001) and Hoffmann et al. (2012) give an overview of how the considered heating and cooling rates are computed in our simulations. In contrast to the method described there, however, here we do not restrict our computations to the stationary case, but instead account for temporal changes of the energy content in a grid cell with a time-dependent approach¹⁹ in parallel with the time-dependent computation of the ionization structure (cf. Sect. 2.2.1). At the present stage we consider the ionization of only hydrogen and helium and therefore neglect cooling processes by metals (e.g., collisional cooling by the O iii ion) which are dominant for metallicities characteristic of population II and population I stars. In models which simulate an already partly metal-enriched universe we therefore approximate the temperature from case to case with plausible values.

2.2.4. The spectral energy distribution of the radiation field

The simulation of the evolution of the ionization structures – obtained via an iteration of the radiative transfer and the rate equations – requires a representative set of frequency points at which the radiative transfer equation has to be solved. The total number of frequency points which has to be considered in this regard depends on the atomic physics which determines the distribution of the frequency points over the relevant spectral range²⁰. As only hydrogen and helium have been treated as elements for the simulations presented in this paper, just hundred frequency points had yet to be considered. Table 1 shows the distribution of the required points over the relevant spectral range, and Fig. 4 demonstrates the effect of the wavelength-dependent absorption cross-section on the transition of neutral to ionized hydrogen. The images in the figure show the ionization structure obtained for blackbody sources of different temperatures, but the same number of hydrogen-ionizing photons per second (10⁴⁹)²¹. For a soft ionizing spectrum (e.g. T_eff = 20 000 K) most of the ionizing photons have energies slightly above the Lyman-edge, and because of the high value of the ionization cross section at this frequency, the ionization front is displayed by a thin shell (this is the result one usually obtains if the gas is illuminated by a monochromatic source). In contrast, a hard ionizing spectrum (e.g. T_eff = 150 000 K) produces a mean energy of the photons which is considerably above the Lyman-edge, and since the ionization cross section has lower values at such frequencies, the photons penetrate deeper into the gas – because of their larger mean free paths $\begin{array}{ccc}{\mathit{\lambda }}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{\sum _{\mathit{i}}{\mathit{\chi }}_{\mathit{i}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}\mathit{,}& & \end{array}$(27)(where the total opacity χ_ν has contributions χ_i(ν) from each type of absorbing atom/ion i) – and the “ionization front” is therefore displayed by a thick shell. This well-known effect can be reproduced only by models with a sufficient number of frequency points in the ionization continuum.

3. Implementation and optimization details of the algorithm

The physical equations described in Sect. 2 are the guiding principles for the algorithm used to calculate the radiative transfer in a homogeneous or inhomogeneous gaseous medium. The numerical procedures which realize this algorithm are based on an iterative solution of the time-dependent radiative transfer (cf. Sect. 3.2) and the time dependent statistical description of the microphysical state of the gas, which together characterize the ionization structures in the medium surrounding the clusters of ionizing radiation sources with their characteristic properties and spatial distribution. The numerical procedures in turn are based on discretization schemes, the requisite accuracy of which determine the required lengths of the time steps, the extent of the cubic cells, the number and orientation of the rays, and the number and distribution of the frequency points. Using these iteration and discretization schemes the mean intensities are calculated per frequency interval for every cell by tracing all rays in sequential order, and from these intensities the occupation numbers are determined via a solution of the rate equations for each cell and time step. Thus, the implementation of the iteration cycle is governed by the geometrical aspects of the ray by ray solution and the time-dependencies of the equation schemes. To specify these items we will in the following section describe how the geometrical structure in our simulations is established in detail. In this regard special emphasis will be given to the ray distribution used for the radiative sources, and to the tedious calculation of the lengths of the ray segments within the cells (Sect. 3.1). The accuracy of our method will finally be tested in this section by some benchmark tests (cf. Sect. 3.3).

3.1. Ray tracing in a Cartesian grid

The fundamental geometrical task in grid- and ray-tracing-based 3-dimensional radiative transfer is to distribute the rays correctly. The two main geometrical aspects to meet this task are (a) tracing the segments of each individual ray correctly (cf. Sect. 3.1.1); and (b) implementing the isotropic character of the radiation field of every point source by distributing the rays as evenly as possible (cf. Sect. 3.1.2), furthermore assuring that each cell is traversed by a sufficient number of rays²² (cf. Sect. 3.1.3). With respect to the second concern we have applied and tested two different approaches: a conventional latitude/longitude method where the directional vectors are equally-spaced in polar angle, and for each polar angle equally-spaced in azimuth (cf. Abel et al. 1999), and a method based on dividing the surface of a sphere into elements of equal area (Górski et al. 2005) and using the directions of the centers of these areas as directional vectors for the rays. As these tests have shown that the latter method distributes the rays more evenly and thus leads to less discretization artifacts, we have decided to use this one.

3.1.1. Tracing the segments of a single ray

A source at coordinates (x_S,y_S,z_S) embedded in a Cartesian grid of cubic cells of identical size emits rays whose directions are defined by ϑ, the angle between the directional vector and the “north pole” (positive z axis), and the azimuthal angle ϕ (measured from the x axis). To compute the radiative transfer along the rays we have to know which cells are passed by which rays and how long is the segment of a ray that crosses a cell. To obtain this information we follow each ray from its source to the edge of the simulation volume and determine in sequence which border between grid cells will be crossed next by the ray. This is done via a calculation of the distances r_yz, r_xz, r_xy to the next yz-, xz-, xy-plane crossings in the direction of the ray. (For this purpose it is convenient to choose a coordinate system where the unit is equal to the length of a cube; in the simulation the actual physical length of a grid cell is then given by a simple scaling law.) The values r_yz, r_xz, r_xy thus indicate the distances from the origin of the source to the next point where a grid cell is crossed and the shortest of these distances is the one which defines the ray segment l(m) (cf. Eqs. (4) and (25))²³.

3.1.2. Distribution of the rays

As the radiation field of an ideal point source is isotropic, any discretization of the radiation field should approximate this as close as possible and distribute the rays as uniformly as possible. Furthermore, as the luminosity of the source has to be divided among all rays, it would be highly advantageous if all rays were to correspond to approximately the same solid angle as seen from the source. Thus we have selected a method where the surface of a sphere around the radiation source is split into almost uniformly distributed segments of equal area, so that the directional vectors to the segment centers correspond to equal solid angles and are almost isotropically distributed (“hierarchical equal-area iso-latitude pixelization” method, HEALPix; Górski et al. 2005). These vectors are then used as the directions of the rays²⁴.

3.1.3. Exploitation of symmetries

One of our modifications of the algorithm outlined above for computing the geometry of the rays regards the symmetry of the ray distribution of the individual sources. Because of this symmetry it is sufficient to compute the directional vectors just in one octant and to mirror the calculated geometrical factors at the planes defined by the coordinate axes²⁵. By choosing the size of the octant as large as the entire simulated volume, the same geometrical information can be used for all sources assumed within the volume.

Besides the distribution of the rays the actual number of rays used is essential for the quality of the radiative transfer, because higher ray densities lead to a better representation of isotropy with less discretization artifacts, and in the (theoretical) limit of an infinite number of rays, the radiation field defined by the rays correctly describes the true (continuous) radiation field. Contrariwise, isotropy can be badly violated by using an insufficient number of rays, even if a good distribution algorithm is chosen. The dependence of the accuracy of the radiative transfer on the number of rays is demonstrated in Fig. 5.

For this test case we have chosen a single ionizing source embedded in a homogeneous isothermal gas. Because of symmetry we know that the resulting ionized volume in the gas must be spherical. (As the used geometry is independent of the density structure of the gas, the results regarding the minimum ray density necessary for an accurate description are of course also applicable for nonsymmetric problems as well.) We have found that the deviations from the ideal sphere are maximal along the planes defined by coordinate axes. This is the case for both the HEALPix method as well as the equal polar angle intervals method. For an accurate representation (within the limits of the cartesian cell grid) of the true shape of the ionized volume we need at least 10 rays per cell even in the corners of the simulation volume farthest removed from the source.

3.2. Accounting for the finite speed of light

Although the geometrical properties of the rays are essential for a correct description of the radiative transfer, the evolution of the ionization structures are primarily determined by the ionizing luminosities of the sources; and this means that for conditions where the ratio of photon emission rate to the gas density is very high the radius of the physical ionization front can propagate even at velocities which are close to the speed of light c. For cases where recombinations are extremely rare (e.g., at low gas densities) this behavior can even analytically be described: As essentially every emitted photon ionizes one hydrogen atom, which under these idealized conditions then remains permanently ionized, the volume of the ionized region, V = (4π/3)r³(t) as a function of time, multiplied by the particle number density, n, will be equal to the total number of emitted photons, Ṅ·t. Thus, assuming spherical symmetry $\begin{array}{ccc}\mathit{Ṅ}\mathrm{·}\mathit{t}\mathrm{=}\mathit{V}\mathrm{·}\mathit{n}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{r}}^{\mathrm{3}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}\mathit{n}\mathit{.}& & \end{array}$(28)The time-dependent radius of the ionization front and its temporal derivative therefore are $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{3}\sqrt{\frac{\mathrm{3}\mathit{Ṅ}\mathit{t}}{\mathrm{4}\mathit{\pi n}}}\mathit{,}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{3}\sqrt{\frac{\mathrm{3}\mathit{Ṅ}}{\mathrm{4}\mathit{\pi n}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{t}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{3}\sqrt{\frac{\mathrm{3}\mathit{Ṅ}}{\mathrm{4}\mathit{\pi n}}}{\mathit{t}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{·}& & \end{array}$(29)But, at the beginning of the irradiation this simple model predicts an expansion rate of the ionization front which is faster than the speed of light: $\begin{array}{ccc}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{t}}\mathit{>}\mathit{c}\mathrm{for}\mathit{t}\mathit{<}\frac{\mathrm{1}}{\mathrm{3}}\sqrt{\frac{\mathit{Ṅ}}{\mathrm{4}\mathit{\pi n}{\mathit{c}}^{\mathrm{3}}}}\mathrm{·}& & \end{array}$(30)Since physical ionization fronts can not really propagate faster than the speed of light, we avoid such quasi-super-luminous expansions of the ionization fronts by delaying to the correct time step the computation of the radiative transfer for cells that can not in physical reality be reached by light emitted by the source since its appearance. This method is reasonable if, after being turned on, the luminosities of the sources do not vary considerably during the light travel time through the ionized volume. Figure 6 shows the speed of the ionization front as a function of its radius computed using our method, taking both recombination as well as the above described correction procedure into account.

3.3. Benchmark tests

To test the accuracy of our method we have performed a series of calculations which we compare in this section to an analytical approximation, results from other 3-dimensional radiative transfer codes published in the literature (Iliev et al. 2006), and a detailed spherically symmetric model.

3.3.1. Temporal expansion of a Strömgren sphere in a homogeneous hydrogen gas

A very simple test case is offered by the expansion of a Strömgren sphere in a homogeneous isothermal medium consisting of pure hydrogen, since under these approximations the temporal behavior of the sphere can be described analytically $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{3}\sqrt{\frac{\mathrm{3}\mathit{Ṅ}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{n}{\mathit{\alpha }}_{\mathrm{B}}\mathit{t}}\mathrm{)}}{\mathrm{4}\mathit{\pi }{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{B}}}}\mathit{,}& & \end{array}$(31)where Ṅ is the number of emitted ionizing photons per time unit, a_B is the case-B recombination coefficient and n_H is the number density of hydrogen (Iliev et al. 2006). In Fig. 7 we compare this analytical solution with the result of a series of our simulations produced by a variation of the parameter f which controls the time step size (cf. Eq. (15)). The comparison shows on the one hand that for t → ∞ the asymptotic value r(t) is almost identical with the analytic solution, on the other hand it is also shown that the usage of larger time steps leads to slower expansions of the simulated Strömgren spheres. The explanation of this behavior is based on the fact that while the ionization fractions are updated during a time step the radiative transfer is not, and reflects the (higher) opacities at the beginning of the time step. Therefore, the number of photons reaching the cells just behind the ionization front is underestimated for long time steps. (Still, very large time steps can be used if the aim of the simulation is to compute the equilibrium state of the irradiated gas instead of an accurate description of the temporal evolution of the ionization structure.)

3.3.2. Comparison with other 3-dimensional radiative transfer methods

To compare the ionization structure from our code with that from other 3d radiative transfer codes we have computed models using the parameters specified in the cosmological radiative transfer comparison project (Iliev et al. 2006). Figure 8 shows our result for the ionization state of a pure hydrogen gas enclosed in a box of the size of 6.6 kpc and illuminated by a monochromatic ionizing source in comparison with other simulations, namely the C²-Ray simulation, which also uses the ray-tracing method, the CRASH simulation, which utilizes a Monte Carlo technique, and the OTVET simulation, which uses the variable Eddington tensor method. Although the results of the various methods do not reveal identical shapes, the fact that the structures of the C²-Ray simulation on the left hand side of the figure and ours on the right hand side of the figure are almost indistinguishable is a mutually satisfactory result, since both procedures are based on the same technique. (Other codes based on ray tracing give similar results to these two.)

3.3.3. Comparison of Cartesian and spherically symmetric models

For a test which also includes helium as a gas component in our 3d radiative transfer model we used a spherical symmetric equilibrium model (cf. Hoffmann et al. 2012) as reference. We present two sets of models, one using a hydrogen number density of 10 cm^-3 (a typical value for H ii regions) and one using a hydrogen number density of 10^-3 cm^-3 (a typical value for a cosmological density field at early epochs). (As usual we set the He number fraction to n_He/n_H = 0.1.) For both gas densities we show in Fig. 9 resulting ionization structures in the gas: in the upper panels for a T = 30 000 K blackbody emitting 10⁴⁹ hydrogen-ionizing photons per second, representing a late O-type star, and in the lower panels for a 65 000 K blackbody emitting 8 × 10⁵¹ hydrogen-ionizing photons per second, representing a very massive star (cf. Sect. 4.1).

In those models where we make use of the on-the-spot approximation in the spherically symmetric models, we get exactly the same results – within the resolution of the grid – as for our 3d calculations: not only are the Strömgren radii the same, but also the radius-dependent ionization fractions of H and He are almost identical for both the soft and the hard spectrum source. If we don’t apply the on-the-spot approximation and instead treat the diffuse radiation field correctly in the spherically symmetric models, we find deviations on the order of 20% for the occupation numbers of the neutral stages of hydrogen and helium in the inner regions. The reason for this is that part of the reemitted photons is not reabsorbed locally, as suggested by the on-the-spot approximation, but, due to the low number of available absorbers, at larger distances. However, these small deviations are irrelevant with regard to the positions of the ionization fronts and thus the sizes of the Strömgren spheres.

4. Photoionization models for H ii regions around luminous stars in the early universe

With this new method in hand we will present in the following first exemplary simulations of our 3d radiative transfer applied to the reionization scenario. As different characteristics of the SEDs of the sources of ionization in general have a considerable impact on the power of ionization, the application of oversimplified spectra in radiative transfer codes – e.g. black-body radiators which are based on just a few discrete wavelengths – can lead to incorrect or deceptive results. To quantify such possible systematic errors in the context of hot massive stars, we first will perform a comparison of blackbody and realistic SEDs at different metallicities stressing for the ionizing spectra of the relevant stellar populations used as ionizing sources the importance of metal enrichment. As our focus lies on a high spectral resolution of the 3d radiative transfer describing the evolution of the most important ionization structures of the IGM in certain epochs, we will in this concern not just concentrate on a massive star representing the top of a Salpeter IMF, but also on a very massive star (VMS) resulting from a cluster-collapse as a runaway collision merger (cf. Pauldrach et al. 2012, and references therein).

In the introduction we discussed that the reionization process may have occurred in two distinct steps (Cen 2003), where the first one was primarily driven by population III stars and was concluded at a redshift of z ~ 15, whereas the second one, which completed the reionization with radiation from population II/I stars, finished at a redshift of z ~ 6. With respect to time-dependent simulations of the reionization scenario we will therefore perform in a first step multisource simulations for population III stars in homogeneous and inhomogeneous cosmological gas distributions representing the first stage of the reionization process in order to investigate on the one hand the evolution of the ionization fronts into the IGM and on the other hand the fraction of photons which can escape from the ISM surrounding the ionizing sources and thus build up the ionization fronts in the IGM. In an advanced step we further will combine our predicted SEDs of massive and very massive stars with time-dependent 3d radiative transfer simulations of the IGM in order to describe the second stage of reionization. Although our models are still based on simplified initial conditions, the general requirements that might have allowed the much less massive population II and population I stars to fully reionize the IGM – this approach is basically motivated by the increasing star formation rate in this epoch (Lineweaver 2001) – can already be investigated even at this stage. Regarding the reionization of He ii to He iii in particular, we investigate further the influence and relevance of a top-heavy-IMF connected to very massive stars. Whether the appearance of population II and population I VMSs could have been responsible for the reionization of He ii will therefore be examined as a crucial point.

4.1. Comparison of blackbody and realistic spectral energy distributions at different metallicities as ionizing sources

Motivated by the fact that metallicities which are different from zero can have a strong influence on the SEDs of massive stars and that the ionizing fluxes can at certain frequencies ranges therefore depart decisively from those of blackbody radiators we computed comparison models using synthetic SEDs from atmospheric models instead of blackbodies as ionizing sources. To limit the scope of this challenge we do not attempt to trace the metallicity dependent characteristics of the SEDs of a complete set of stellar populations we rather want to demonstrate the effects of metallicity on the spectra and furthermore the ionized structures by means of a series of striking objects.

A massive star representing the top of a Salpeter IMF.

As striking stellar objects are primarily characterized by their mass, we have chosen a 125 M_⊙ star with an effective temperature of T_eff = 50 000 K and a radius of R = 21.7 R_⊙ as a first example. Figure 10 shows SEDs computed for this object assuming different metallicity values (ranging from 0.001 to 1.0 Z_⊙ – the spectra have been calculated on basis of the method described by Pauldrach et al. 2001, 2012). With respect to the mass each of these models represents obviously a star which can be regarded to mark the top of a metallicity dependent Salpeter IMF. As suspected, Fig. 10 shows clearly that an enhancement of metallicity has a drastic influence on the SEDs. While the hydrogen-ionizing flux (λ < 911 Å) is primarily influenced by P-Cygni spectral lines, shorter wavelengths, as the He ii ionization edge at ~228 Å, are also influenced by non-LTE-effects affecting the bound-free thresholds of certain ionization stages (cf. Pauldrach et al. 2012) and producing in comparison with blackbody radiators considerable differences.

The results of our 3d radiative transfer calculations for the radial density distribution of ionized hydrogen surrounding the 125 M_⊙ stars with different metallicities are shown in Fig. 11. (The parameters of the model correspond to a typical H ii region: n_H = 10 cm^-3, n_He/n_H = 0.1, and T_gas = 10 000 K). Although at the effective temperature of 50 000 K of the ionizing source the difference of the solar metallicity spectrum to the blackbody spectrum leads to a considerable difference of ~4% for the hydrogen Strömgren radius – corresponding to a difference of ~12% in the ionized volume –, the variation of the metallicity of the stars does not have a huge influence on the number of hydrogen-ionizing photons, and therefore a replacement of especially the lowest metallicity spectrum (0.001 Z_⊙) by a blackbody spectrum, which leads to a radial error of ~3% which corresponds to a maximal error of ~9% in the ionized volume, may appear to be legitimate.

However, the similarity of the hydrogen Strömgren spheres presented here results primarily from the fact that the influence of the metals in the gas and therefore the corresponding cooling and heating processes have been neglected. A consistent computation including metals would actually lead to smaller temperatures for larger metallicities and thus result in smaller Strömgren spheres (Osterbrock & Ferland 2006). An even larger influence on the ionization structure of the metals in the gas and thus on the gas temperatures and Strömgren spheres will be caused by the non-LTE line blocking and blanketing²⁶ in the expanding atmospheres of the massive stars which ultimately determine the shapes of the SEDs (cf. Sellmaier et al. 1996; and Pauldrach et al. 2001 – a superficial impression of this behavior can also be obtained by a thorough inspection of Fig. 10). From these arguments it is quite obvious that our Fig. 11 does not yet represent an ultimate result!

Regarding the radial density distribution of He iii differences of the Strömgren radii which originate from different metallicities of the ionizing sources are considerably more pronounced (cf. Fig. 13) – as an example, the difference of the blackbody to the spectrum calculated with solar metallicity leads to a difference in the He iii Strömgren radius by a factor of 30. But these differences are not necessarily a monotonic function of increasing metallicity of the ionizing sources. With respect to our calculations the number of emitted He ii-ionizing photons increases from the model with Z = 10^-3 Z_⊙ to the model with Z = 0.2 Z_⊙, whereas it decreases drastically – by 4 orders of magnitude – for the model with solar metallicity (cf. Fig. 12 and Table 2, which show that there is no strict relationship between the metallicity and the emission of He ii-ionizing photons). These results show clearly that the number of He ii-ionizing photons is not just strongly influenced by the effective temperature and hence the mass of the relevant objects, but also by the evolution of the metallicity. Employing blackbody spectra in general therefore leads for the cases of helium and metals whose ionization energies exceed those of H i and He i to inaccurate ionization structures. As exactly the elements which are affected by this approximation are producing the characteristic and in principle observable emission lines, such a weak-point in the computation of the ionization structures would make as a consequence a comparison to observations extremely disputable (cf. Rubin et al. 1991, 2007; Sellmaier et al. 1996).

A very massive star as a runaway collision merger from a cluster-core collapse.

As a second example of our investigation of the influence of realistic SEDs of massive stars on the ionization structures of the surrounding gas we have chosen an even more impressive stellar object: a 3000 M_⊙ star. Such stars are expected to form as runaway collision mergers in cluster-core collapses, and via feedback effects due to their enormous radiation and strong stellar winds they naturally have a considerable impact on their environment and the local star formation (cf. Sect. 1).

Figures 14 and 16 show the SEDs of three model atmospheres of such a star (L = 10⁸ L_⊙, T_eff = 65 000 K; cf. Pauldrach et al. 2012) computed for different metallicities (Z = 0.05 Z_⊙ and Z = 1.0 Z_⊙) and different hydrogen mass fractions X (70% and 40%). Together with the consistently calculated mass loss rates the ionizing fluxes of the model stars are listed in Table 3. The luminosities of these VMSs exceed those of massive O stars at the top of a Salpeter or Kroupa-IMF by one to two orders of magnitude.

Although the hydrogen-ionizing fluxes of these models are close to those of a corresponding blackbody (Fig. 14 and Table 3), increasing metallicity and/or increasing helium mass fraction leads to a clear decrease of the He ii-ionizing fluxes: while the metal-poor star with a normal helium content has a He ii-ionizing flux 29% larger than that of a corresponding blackbody, for the helium-enriched star with solar metallicity it is 53% smaller (Fig. 16 and Table 3). As a consequence of this behavior the hydrogen Strömgren radii differ from the blackbody Strömgren radius by only 3% (Fig. 15), but the He iii Strömgren radii show differences to the corresponding blackbody value of up to 25% (Fig. 17).

The differences between our very massive stars and a black body are small, however, when compared to the differences between our very massive stars and our normal 125 M_⊙ massive stars. The hydrogen-ionizing fluxes of the 3000 M_⊙ stars are a factor of 40 larger than those of the 125 M_⊙ stars, while the He ii-ionizing fluxes of the 3000 M_⊙ stars exceed those of the 125 M_⊙ stars by a factor of 800 to several 10⁶, depending on the metallicity (Tables 2 and 3). Due to the large luminosities and effective temperatures of the very massive stars the radii of the hydrogen Strömgren spheres thus differ by a factor of 3.5 and those of the He iii Strömgren spheres by a factor of up to 100 (Figs. 15 and 17).

We conclude from this that the transition from normal massive stars to very massive stars could strongly influence the relative volumes of He iii to H ii. This assessment does not change if the fluxes are normalized to the same hydrogen-ionizing photon rate: in this case the He ii-ionizing fluxes of the VMSs are still at least a factor of 20 larger. As this finding may be of importance for a possible scenario of the cosmic He ii-reionization, we will return to this point in Sect. 4.2.2.

4.2. Time-dependent simulations for different characteristics of the reionization scenario

To better understand the reionization process of the universe, the evolution of the ionization fronts expanding in the IGM has to be examined in some detail using a numerical treatment of the radiative transfer and the microphysics of ionization and recombination. For this purpose we present in this section typical results of our 3d radiative transfer simulations for representative values of the densities and star formation rates. In these simulations we use clusters of massive population III stars (for the simulations representing the first stage of reionization) and population II/I stars (for the second stage of reionization) as ionization sources, surrounded by an environment reflecting the intergalactic cosmological gas.

4.2.1. Multisource simulations for population III stars in homogeneous and inhomogeneous gas distributions

As already mentioned in the introduction, the time around z ~ 13.5 marks the transition from population III to population II/I stars and leads to a decline in the emitted radiative power. On the other hand, the star formation rate of population III stars was highest around z ~ 15, shortly before this stellar generation died out. For the first part of the reionization process that was driven by population III stars we have therefore selected this redshift for our calculations.

First stage of the reionization scenario in a homogeneous gas distribution.

Representing the first stage of the reionization history we have calculated models which simulate a volume of (330 kpc)³ in proper coordinates at z = 15 (which is a volume of (5.3 Mpc)³ in comoving coordinates at z = 0) resolved into 101³ grid cells. The temperature of the gas has been kept constant (T = 10 000 K, cf. Cen 2003), and the simulations use randomly distributed clusters of massive population III stars surrounded by a primordial environment characterized – as a rough approximation of the real situation – by a homogeneous intergalactic cosmological gas with a hydrogen number density of 8 × 10^-4 cm^-3 (reflecting the mean baryonic density of the universe at the given redshift, cf. Yoshida et al. 2004). The spectra of the multiple ionizing sources with a representative temperature of 95 000 K have been approximated by blackbodies (at such a high temperature and a negligible metallicity the number of hydrogen-ionizing photons is almost perfectly described by blackbody emitters, cf. Sect. 4.1). Although we have not accounted for structures on smaller scales in this simulation, we have compared the assumed distances of our clusters with the distances between the largest clusters investigated by Gnedin & Bertschinger (1996) and found reasonable agreement. Each of our ten clusters emits 4 × 10⁵³ hydrogen-ionizing photons per second, which corresponds to a star formation rate of 9.5 × 10^-3 M_⊙ yr^-1 per comoving Mpc³ (matching the rate deduced from Cen 2003, and Hernquist & Springel 2003 at a redshift of z = 15) and represents the emission of 4000 population III stars with a mass of 100 M_⊙ and a lifetime of 3 Myr (El Eid et al. 1983; Schaerer 2002)²⁷.

The evolution of the ionization fronts expanding into the IGM in our simulation is shown in Fig. 18. This simulation predicts that the clusters of population III stars as predicted by the theoretical star formation rate can ionize the intergalactic hydrogen on a timescale of just t ~ 50 Myr (corresponding to Δz ~ 1.5 at z = 15). Although cosmic variance could modify the reionization timescales when going to larger box sizes, this result is nevertheless significant, since our chosen volume with regard to the assumed and the observed mean distances of the cluster structures in this epoch is representative for the first stage of the reionization history of the universe. The time around z ~ 15 is the “optimum” period for this reionization phase since in this epoch the star formation rate of population III stars (Cen 2003), and thus the photon emission rate, rises rapidly, just before the changeover to population II stars significantly cuts down on the emission rate of ionizing photons. Thus, our results indicate that the small contributions to the reionization before z ~ 15 are mostly insignificant for ionizing the universe completely, and the exact time of the beginning of the process is not essential.

As the large temperatures of massive population III stars result in considerably hard spectra, it is not surprising that the timescale of the reionization of He ii into He iii is not much different from the timescale of the reionization of hydrogen (each of our ten clusters emits 2 × 10⁵² He ii-ionizing photons per second, which is 1/20 of the amount of hydrogen-ionizing photons per second). Thus, if the first stage of reionization had really been driven, as currently suspected, by a stellar population dominated by massive hot primordial stars, helium would also have been reionized to a large extent into He iii at a redshift of z ~ 13 (cf. Fig. 19). This would not have been the case, however, if the population III stars had a much softer spectrum (cf. Sect. 4.2.2) – as is indicated by some simulations which show a stronger fragmentation even for primordial gas (Clark et al. 2011). Whether helium had also been reionized into He iii during the first stage – by a redshift of z ~ 13 – therefore depends crucially on the emergence of massive population III stars in sufficient numbers. Whether these objects had been the dominant species of the stellar population that illuminated the very early universe can only be discovered by observations. In this regard the James Webb Space telescope as well as the European Extremely Large Telescope (E-ELT), the Giant Magellan Telescope (GMT) and the Thirty Meter Telescope (TMT) are expected to be able to observe clusters of bright population III stars at redshifts of z ≈ 13...15, and thus will help to improve our present knowledge of the stellar content governing the first stage of reionization.

Details of the reionization process in an inhomogeneous cosmological density field.

It is obvious that the assumption of a homogeneous universe as applied to our models above is just a rough approximation of the real situation. As a consequence of this approximation effects which depend on small-scale inhomogeneities of the density structure of the gas and which in turn may considerably affect the reionization process can not be taken into account. Examples of such effects are the interrelation of the density structure and the local star formation rate – mutual back-reactions of the processes involved give rise to the inhomogeneous behavior of the ISM and are responsible for the presence of H ii regions –, and the “porosity” of the medium which determines the escape fraction of ionizing photons that can exit their local environment and reach the intergalactic regions.

In order to study in particular the latter point, and thus the influence of multiple sources on their environment consisting of an inhomogeneous density structure on smaller scales, we have calculated a model which simulates in proper coordinates a volume of (50 kpc)³ resolved into 128³ grid cells. The density structure used for this simulation is based on results of hydrodynamical simulations from an N-body gas dynamics code (cf. Ryu et al. 1993; Iliev et al. 2006). As this density profile had been calculated for cosmological conditions which are consistent with a redshift of z = 9 only, we rescaled it to conditions corresponding to a redshift of z = 15. The sources, which are located at the peaks of the structured density profile (pink dots in the upper left image of Fig. 20), cover 16 smaller clusters with a total emission rate of ~4 × 10⁵³ hydrogen-ionizing photons per second (thus, all of the clusters in this simulation collectively represent one of the clusters in the homogeneous case (cf. Fig. 18)). As before, they are treated as blackbody radiators representing stars with an effective temperature of 95 000 K.

The results of this simulation, for which the gas temperature has again been held fixed at a value of 10 000 K, are presented in Fig. 20 for hydrogen and Fig. 21 for helium. The overall impression of the images reflects the behavior of considerably more complex ionization structures for the inhomogeneous density structure compared to those of the homogeneous models. This impression is certainly true and the irregular form of the ionization fronts and the larger fraction of neutral hydrogen in the higher-density filaments is evident. The simulations also show that the ionization stages He ii and He iii coexist over a considerable part of the volume (cf. Fig. 21).

The most important aspect of the inhomogeneous structure of the ISM is its “porosity”, and the characteristic parameter relevant for cosmic reionization simulations is the fraction of emitted ionizing photons that can escape the local environment and reach the intercluster medium and thus affect the reionization process. This so-called escape fraction is defined as the ratio of the number of photons not absorbed in a specified local reference volume to the total number of emitted photons. To calculate escape fractions for realistic scenarios of density distributions surrounding ionizing sources, we have constructed a number of fractal density structures as proposed by (Wood et al. 2005; see also Elmegreen & Falgarone 1996). These structures are meant to represent the even smaller scales of star formation regions and evolved starburst clusters, as, for example, NGC 3603 and its surrounding thin ionized gas.

Figure 22 shows the H ii and He iii ionization structures and the escape fractions as functions of the photon energy for three different representative models. In each case the central radiation source emits 8.5 × 10⁵² H i-ionizing and 4.0 × 10⁵¹ He ii-ionizing photons per second (corresponding to a blackbody with a temperature of 95 000 K). With a mean hydrogen number density of 3.3 cm^-3 the first model corresponds to an H ii region associated with star-formation clouds in the local universe (a volume of (166 pc)³, somewhat larger than the extension of the Carina nebula with a diameter of approximately 100 pc, cf. Smith et al. (2000), but smaller than NGC 604, cf. Melnick 1980). In this simulation, both hydrogen and helium are almost completely ionized and the majority of both H i-ionizing and He ii-ionizing photons are able to escape from the simulated volume. In the next simulation the mean gas density has been retained, but the volume has been extended to (820 pc)³, which is approximately twice the size of NGC 604. In this case, less than 5% of the H i-ionizing and almost none of the He ii-ionizing photons are able to escape the volume. The last model has a mean hydrogen number density of 0.12 cm^-3 and a volume of (2.5 kpc). This model approximates a typical spiral arm surrounded by the thin diffuse ionized gas (DIG) that in present-day disk galaxies contains a large portion of the interstellar ionized gas (cf. Greenawalt et al. 1998; Hoffmann et al. 2012, and references therein). In this simulation the gas is almost transparent for soft H i-ionizing radiation, but it absorbs approximately 50% of the He ii-ionizing radiation (the escape fraction depends on the photon energy and reaches its minimum for photons whose energy is just above the ionization edge of He ii). Within the energy intervals defined by the ionization edges, the frequency-dependent ionization cross section (cf. Sect. 2.2.2) leads to a radiation hardening, i.e. the escape fraction rises with larger photon energies. We note that the numerical values of the escape fractions can vary strongly even for small changes of the density structure and the properties of the ionizing sources, such that our results are just estimates for the escape fractions.

4.2.2. Simulations for the second stage of reionization using SEDs of massive and very massive stars

In this section we will focus on a high spectral resolution description of the 3d radiative transfer applied to the evolution of the ionization structures in the second stage of the reionization process, which was completed for hydrogen by a redshift of z  ~  6 and for helium by a redshift of z ≈ 3 (cf. Sect. 1). As sources of ionization we consider metal-enriched population II stars, which had already come into being during this epoch. Although these objects are in general significantly cooler and less massive than the population III stars, they are numerous and powerful due to the strongly increased star formation rate during this epoch (Lineweaver 2001; Barger et al. 2000). Furthermore, there are indications that these stellar populations may have had top-heavy-IMFs which additionally generated VMSs (cf. Sect. 1).

With regard to reionization we simulate a volume of (5.3 Mpc)³ in comoving units and assume (as in our simulations for the first stage of reionization) a homogeneous gas density corresponding to the mean baryonic density of the universe at this epoch. Although these approximations are quite restrictive (cf. Sect. 4.2.1) we are convinced that the principal circumstances that would allow the much less massive population II and population I stars to fully reionize at least the hydrogen content of the IGM can be investigated in their essence with this approach.

Scenario based on stars representing the top of a Salpeter IMF.

In this first step we combine the calculated SEDs of massive stars (cf. Fig. 10) with our time-dependent 3d radiative transfer simulations. The ionizing sources in these simulations represent clusters of stars, and their total hydrogen-ionizing photon emission rates are modeled on basis of a cluster with a Salpeter IMF. We have, however, approximated the SED of the cluster by the SEDs of its most luminous constituents, which we take to be 125 M_⊙ stars.

Depending on their evolution, galaxies and even different star forming regions within a single galaxy may have significantly different metallicities of their gas and, as a consequence, also of their young stellar populations²⁸. We have therefore not used the same SED for all sources but instead use that of a 125 M_⊙ star with a metallicity of Z = 0.1 Z_⊙ for 5 of the 10 clusters, and that of a Z = 1.0 Z_⊙ model for the others. Based on our simulations in Sect. 4.2.1 we have chosen an escape fraction of f_esc = 0.1 as a plausible value for stars embedded in H ii regions.

To model the luminosity of the clusters, we have used the population synthesis code Starburst99 (Leitherer et al. 1999, 2010) in continuous star formation mode on the basis of the mean cosmological star formation rate at the respective redshift (cf. Lineweaver 2001; Barger et al. 2000). The total luminosity of our sources (measured by their hydrogen-ionizing photon emission rates) was chosen to match the emission rates predicted by the Starburst99 model²⁹ for the corresponding volume. In the redshift range 10 > z > 6 we assume a star formation rate of 0.02 M_⊙ Mpc^-3 yr^-1, in agreement with the mean value given for that interval by Lineweaver (2001), resulting in an emission rate of 4.7 × 10⁵³ H-ionizing and 1.0 × 10⁵⁰ He ii-ionizing photons per second. For 6 > z > 3 we increased the star formation rate to 0.05 M_⊙ Mpc^-3 yr^-1, leading to an enhanced emission of ionizing photons (H-ionizing flux of 1.2 × 10⁵⁴ photons per second and He ii-ionizing flux of 2.5 × 10⁵⁰ photons per second). As shown in Fig. 23 the radiation of these cluster populations is able to reionize the hydrogen content of the simulation volume within the redshift interval from z ~ 10 to z ~ 5.8.

For stars with solar metallicity the He ii-ionizing flux is reduced by more than four decades compared to stars with 0.1 solar metallicity (cf. Sect. 4.1) while there is no large variation of the H-ionizing flux. Only 1.4% of the He ii atoms are ionized by the stellar populations in this simulation up to the redshift of z ~ 2.8 (Fig. 24) for which observations indicate that the reionization of He was completed (Reimers et al. 1997) – we note that replacing the spectra of the Z = Z_⊙ clusters by SEDs of low-metallicity (e.g. Z = 0.1 Z_⊙) stars would not help much, since this would increase the fraction of He iii only by a factor of less than 2, when compared to the “mixed” population presented in this section. Thus, it appears plausible that stellar populations with Salpeter-like IMFs are the main contributors to the reionization process of hydrogen, but due to the low emission rates of He ii-ionizing photons they cannot contribute significantly to the ionization of He ii.

Scenario based on stars representing runaway collision mergers.

As the reionization of He ii was considerably delayed compared to the reionization of H i, it appears likely that the sources of ionization responsible for the reionization of He ii are different from those which caused the reionization of H and neutral He (cf. Sect. 1). Runaway collision mergers clearly exceed the assumed upper limit of the mass of present-day massive stars, but they are 1 to 2 orders of magnitude more efficient in producing He iii-ionizing photons than normal hot stars based on a standard IMF, and since these stars can appear especially in chemically evolved clusters of high core mass density, thus at the temporal beginning of the star bursts, they may play an important role in these late stages of the reionization history.

To compare the efficacy of collisional mergers to that of “normal” stars with regard to the reionization of He, our second simulation uses the same star formation rate as before (0.02 M_⊙ Mpc^-3 yr^-1 for 10 > z > 6 and 0.05 M_⊙ Mpc^-3 yr^-1 for z < 6), but instead of a Salpeter IMF we assume that a fraction of 1% of the newly formed stellar mass is involved in mergers to become very massive stars³⁰. (This assumption is based on the observationally determined ratio between the central black hole mass and the stellar mass in galaxies with active galactic nuclei at redshifts of z ~ 4 (Targett et al. 2012) as well as for the redshift range 1 < z < 2 (Bennert et al. 2011) which is in the order of 1%³¹). For this simulation we approximate the cluster SED by that of a 3000 M_⊙ star with a metallicity of Z = 0.05 Z_⊙ and, given that these objects are not created directly from gas clouds, but are the product of merging processes of already-formed stars, we assume an escape fraction of f_esc = 0.5 on the basis that these objects are able to remove the gas from their environment more quickly than ordinary massive stars due to their extremely strong winds (Pauldrach et al. 2012). The lifetime of such a very massive star is 2 × 10⁶ yr (Pauldrach et al. 2012), resulting in a mean number of 20 VMSs within the simulated volume at any given time for the redshift range 10 > z > 6 (with a total emission of 1.5 × 10⁵³ H-ionizing photons per second and 1.2  × 10⁵¹ He ii-ionizing per second) and 50 VMSs for the redshift range 6 > z > 2.5 (with a total emission of 3.8 × 10⁵³ H-ionizing and 2.9 × 10⁵¹ He ii-ionizing photons per second).

With regard to hydrogen this scenario is similar to the one using the Salpeter IMF, resulting in a reionization completed at a redshift of z ~ 5.7 (cf. Fig. 25). For helium however, the VMSs are much more successful, producing 60% of the photons required to fully ionize He ii up to a redshift of z ~ 2.8 (Fig. 26). As the figure clearly shows, this scenario based on the possible progenitors of IMBHs and SMBHs offers an alternative or at least complementary picture for the pending explanation of the reionization of He.

5. Summary, conclusions, and outlook

The impact of massive stars on their environment is of major importance not only for the evolution of most galaxies but also for the ionization balance of the IGM. Although rare by number, massive stars dominate the life cycle of gas and dust in star forming regions and are responsible for the chemical enrichment of the ISM and IGM. This is mainly due to the short lifetimes of massive stars, which favors the recycling of heavy elements on short timescales. Furthermore, since massive stars mostly group in young clusters, the combined amount of momentum and energy released by these clusters into the ISM controls the dynamical evolution of the gas, leading to the formation of superbubbles and triggered continuous star formation. The evolutionary behavior of these starbursts is obviously controled by the metallicity, as is the steepness of the IMF, which can be more top-heavy than previously thought even in present-day galaxies. A common property among starbursts showing such flat IMFs is a high stellar density in the core of the clusters, and very recently strong evidence has been found that in these cases nature favors the formation of massive stars with masses >100 M_⊙ (Harayama et al. 2008, Crowther et al. 2010). Moreover, during the evolution of dense stellar clusters runaway collision mergers can occur, leading to the formation of VMSs (objects of several 1000 M_⊙) even in chemically evolved clusters (cf. Belkus et al. 2007; Pauldrach et al. 2012), and therefore runaway collision mergers can play an important role in the late stages of the reionization history of the universe, whose detailed description is still one of the key unknowns in present cosmological research.

It is known meanwhile that the first generations of stars have not been the most important contributors to this process, since there is evidence that the universe was reionized in two steps (Cen 2003). The first one occurred – primarily driven by population III stars – up to a redshift of z ~ 15, and after a second cosmological almost-recombination phase (around Δz ~ 5, see Cen 2003) the reionization was driven by other powerful ionizing sources and completed at z  ~  6. Along with an increasing metallicity and star formation rate in this epoch (Lineweaver 2001; Barger et al. 2000), a top-heavy-IMF and VMSs – which are efficient emitters of ionizing photons with approximately 10 times more hydrogen-ionizing photons and 10 000 times more He ii-ionizing photons per unit stellar mass than stars with a Salpeter IMF (cf. Pauldrach et al. 2012) – are determining factors for this process.

In this paper we have investigated whether such an extended view of the evolution of the universe from its dark ages to its final reionized state leads to a better understanding of the reionization process. For this purpose we have developed a 3-dimensional radiative transfer code (based on a ray-tracing method) in order to describe the evolution of the ionization fronts at different length scales. The main focus of this code is a comprehensive description of the time-dependent ionization structure with regard to a high spectral resolution of the ionizing spectra in order to accurately quantify the evolution of the ionization structures of different elements. Although only H and He are treated in the present paper, our algorithm allows the implementation of multiple levels for each ionization stage of the metals, and we are currently working on the inclusion of all important elements as described by Hoffmann et al. (2012). (As metal cooling has a significant influence on the temperature structure of the gas and therefore on the temperature dependent recombination rates and thus on the ionization structure of the most abundant elements, the temperature structure has to be treated consistently). The geometrical aspects and the evolution of the ionization fronts have been extensively tested by comparisons to analytical models, the results of radially symmetric calculations, and tests presented by the “Cosmological Radiative Transfer Comparison Project” (Iliev et al. 2006). We found good agreement in all of these cases, showing that the important physical mechanisms which control the temporal expansions of the ionization fronts in homogeneous and inhomogeneous gas structures surrounding numerous sources of ionization have been implemented accurately.

For a realistic modeling of ionization fronts not only the structure of the ISM and IGM has to be correctly considered, but also the SEDs of the ionizing sources used as input for the radiative transfer. Thus, based on a sophisticated model atmosphere code which treats the physics of expanding stellar atmospheres consistently, an accurate simulation of realistic SEDs from population III to population I stars at various metallicities, temperatures, and masses has been an important aspect of our work. As different characteristics of the SEDs of the sources generally have a strong impact on the power of ionization, the usual application of oversimplified spectra (such as blackbodies) can lead to incorrect or misleading results with respect to the thickness and the extent of the ionization fronts in radiative transfer codes. We have quantified such possible systematic errors in the context of hot massive stars via a comparison of blackbody and realistic SEDs we have calculated at different metallicities. The results we present for the radial ionization structure of hydrogen and helium surrounding a 125 M_⊙ star with different metallicities show that the variation of the metallicity of the star does not have a huge influence on the number of hydrogen-photons. However, this result does not reflect the influence of the metals in the gas which have been neglected up to now. A consistent computation including metals will actually lead to smaller temperatures for larger metallicities and thus result in smaller Strömgren spheres; then, an even larger influence on the ionization structure of the metals in the gas and thus on the gas temperatures and Strömgren spheres will be caused by non-LTE effects acting in the expanding atmospheres of the stars and determining the shapes of the SEDs, since in this case not only the number of ionizing photons is decisive but also their spectral distribution. Regarding the radial density structure of He iii, the differences of the Strömgren radii which result from different metallicities of the ionizing sources are already considerably more pronounced due to this reason. This result clearly shows that the number of He ii-ionizing photons is not just strongly influenced by the effective temperature and hence the mass of the objects, but also by the evolution of the metallicity. Employing blackbody spectra in general therefore leads to – at least in the case of helium and metals whose ionization energies exceed those of H – oversimplified ionization structures which do not represent the real structures. As exactly those elements which are affected by this approximation produce the observable emission lines, such a weakness in the simulation of the ionization structures would make a comparison to forthcoming observations questionable.

As a central application of our models with high spectral resolution we considered a number of different physical scenarios leading to ionizing sources with corresponding characteristics of the SEDs, which we compared to analogous calculations carried out with blackbodies. In particular we considered massive stars with 125 solar masses, representing the top of a Salpeter IMF, and very massive stars with 3000 solar masses which are expected to form as runaway collision mergers in cluster-core collapses. Due to the large luminosities and effective temperatures of the latter objects the emitted number of H i-ionizing photons is, compared to stars at the top of a Salpeter IMF, larger by a factor of 40. This leads to hydrogen Strömgren radii around VMSs which are a factor of 3.5 larger than those around normal massive stars at the top of a Salpeter IMF. The effect is even more pronounced for He iii where the VMS-Strömgren radii exceed those of the 125 M_⊙ stars by a factor of up to 100. The reason for this behavior is of course the ratio of the ionizing fluxes blueward of the He ii-threshold, which favors very massive stars by a factor of millions. Very massive stars could thus strongly influence the ionization structure of their environment via their distinct tendency of the emission of H and He ionizing photons.

In the context of the first stage of the reionization history we have shown results of representative 3d radiative transfer simulations which are based on clusters of massive population III stars surrounded by a primordial environment reflecting a homogeneous intergalactic cosmological gas at z ~ 15. The calculations showed that our population III stars were able to reionize H as well as He on a timescale of just t ~ 50 Myr (corresponding to Δz ~ 1.5 at z = 15) within our chosen volume of (5.3 Mpc)³. As the assumption of a homogeneous universe is just a rough approximation of the real situation, we have further calculated models based on inhomogeneous density structures, simulating volumes that comprise multiple ionizing sources on smaller scales. The “clumpiness” of the gas (including the star-forming density peaks of the ISM) results in considerably more complex ionization structures than those of models with homogeneous gas distributions, and as a consequence the fraction of ionizing photons which can escape their local environment and reach the intercluster medium can differ significantly. To estimate these escape fractions we have modeled inhomogeneous environments with fractal density structures of different sizes and mean densities ranging from the density of star-forming H ii regions to the one of the diluted interstellar gas. In the simulations the escape fractions for the ionizing photons, and thus their ability to leak into the IGM and contribute to the ionization process of the universe, depend crucially on the size of the density peaks the ionizing sources are embedded in, ranging from less than 5% for a very large and therefore radiation-bounded H ii region, to almost 100% for a smaller, matter-bounded H ii region where the hydrogen content is almost completely ionized and the gas is therefore transparent to radiation. We note that even small variations in the number of absorbers can lead to strong changes of the escape fractions; thus, a precise determination of this value would require knowledge of the statistics of the density distribution in the small-scale environment of the sources (i.e., on the order of parsecs to a few kiloparsecs). Furthermore, for a given density distribution the escape-fractions are a function of the photon energy, and reach their minimum just blueward of the ionization edge of He ii. The SED of the ionizing radiation escaping into the IGM may thus be different from the intrinsic SEDs of the stellar sources.

Although the initial conditions with respect to the structure and the content of the clusters are not yet precisely known, the result that population III stars have been able to reionize H and He within the generally accepted time frame is not spectacular, since it just reflects the current state of knowledge. The more interesting question regards the second stage of reionization – which completed at a redshift of z ~ 6 (for H) resp. at a redshift of z ~ 2.8 (for He) –, since it is by no means clear whether the metal-enriched and thus significantly cooler and less massive population II and population I stars, given a reasonable cluster structure and composition, have been able to fully reionize the hydrogen and helium content of the IGM. To investigate this scenario we have combined our predicted SEDs of massive and very massive stars with our time-dependent 3d radiative transfer simulations of the IGM. (Note that while the assumption of blackbody sources is not too bad an approximation for stars of primordial composition, it is not an adequate description of the real SEDs of later stellar populations. For the study of the second stage of reionization the results of realistic model atmospheres – as presented here – must therefore serve as input spectra.)

Based on plausible numbers for a standard population we have approximated in a first step the SED of a cluster with a Salpeter IMF by the SED of its most luminous members, namely 125 M_⊙ stars. Our calculations show that these objects have certainly been able to reionize H again (for our chosen volume the process required a timescale of Δt ~ 500 Myr, corresponding to the redshift interval between z ~ 10 and z ~ 6), but not He ii (not even up to a redshift of z ~ 2 – note that the observations show that the reionization of He ii was considerably delayed compared to the reionization of H). This latter finding is not unexpected, given that other investigations (Wyithe & Loeb 2003; Gleser et al. 2005; McQuinn et al. 2009) had concluded from similar considerations that the sources of ionization responsible for the reionization of He ii are different from those which caused the reionization of H.

To investigate whether the appearance of very massive, already metal-enriched population II and population I stars could have been responsible for the reionization of He ii we replaced the 125 M_⊙ stars by runaway collision mergers of M = 3000 M_⊙, which clearly exceed the assumed upper mass limit of present-day massive stars and which are up to 2 orders of magnitude more efficient in producing ionizing photons than normal hot stars based on a standard IMF are. The result is striking: In contrast to the case where the stellar populations follow a Salpeter or Kroupa IMF our simulations now show that the universe could have been reionized within approximately 2 Gyr (corresponding to a redshift of z ~ 2.5) in an environment where a fraction of 1% of the mass formed into stars consists of mergers resulting in very massive stars. We note that this percentage is consistent with the mass ratio obtained for IMBHs and SMBHs relative to the stellar mass found in active galaxies and globular clusters like ω Centauri, and that such VMSs are possible progenitors of IMBHs and SMBHs. This result therefore offers an alternative or complementary scenario for the explanation of the reionization of He ii to He iii, and implies that a more realistic simulation of the entire reionization process of the universe must take the complete set of stellar populations into account and not just the primordial generation of stars.

It will therefore be necessary to understand the chemical evolution and the structure formation at the corresponding epochs more precisely and to establish comprehensive grids of realistic SEDs representing the time dependent content of the clusters and the cluster structures. Our future work will further focus on a more sophisticated description of the 3d radiative transfer with respect to high spectral resolution in order to quantify the evolution of all relevant ionization structures accurately. On basis of synthetic emission spectra calculated for the 3d structure of the ionized gas we will try to determine possible observational features (such as line strength ratios from different ions) that may be used to discriminate between the details of different reionization scenarios. Knowledge of such features will be important for comparisons with more detailed observations at high redshifts, which are expected to become available during the next years.

---

Acknowledgments

We thank an anonymous referee for helpful comments which improved the paper and we further wish to thank I. Iliev for granting us access to the data of the Cosmological Radiative Transfer Codes Comparison and N. Gnedin for providing us with the IFrIT program used for some visualizations of our results. The calculations were performed on high-performance computing facilities at the LRZ (Leibniz Rechenzentrum) and the University of Munich which is to be acknowledged. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under grant PA 477/9-1 and PA 477/18-1.

References

All Tables

All Figures