AN IMPROVED TREATMENT OF COSMOLOGICAL INTERGALACTIC MEDIUM EVOLUTION

The evolution of galaxies is intertwined with that of the intergalactic medium (IGM). Mechanical heating of IGM by active galactic nuclei (Bower et al. 2006; Croton et al. 2006) and radiative heating by X-rays produced in supernovae (White & Rees 1978; Dekel & Silk 1986; Cole 1991; Lacey & Silk 1991; White & Frenk 1991; Oh & Haiman 2003) together with ionizing photons emitted by young stars (Ikeuchi 1986; Rees 1986; Miralda-Escudé & Ostriker 1990, 1992; Shapiro et al. 1990; Efstathiou 1992) modify the temperature and ionization state of the IGM, which in turn alters subsequent galaxy formation.

The physics involved in the coupled evolution of IGM and luminous sources is so complex and covers such a wide range of scales that its treatment involves important approximations. In fact, most studies focusing on galaxy formation adopt an IGM with fixed ad hoc properties. Only studies of reionization follow the IGM evolution in more or less detail.

IGM evolution is described by a couple of differential equations for its ionization state and temperature with some source functions provided by a galaxy model. It is in this latter part where most approximations and simplifying assumptions are made, depending on the particular approach followed, namely hydrodynamic simulations (Quinn et al. 1996; Navarro & Steinmetz 1997; Weinberg et al. 1997; Ciardi et al. 2000; Wyithe & Loeb 2003; Iliev et al. 2007; Okamoto et al. 2008; Trac et al. 2008; Battaglia et al. 2013; Sobacchi & Mesinger 2013a), numerical and seminumerical simulations (Zhang et al. 2007; Faucher-Giguère et al. 2009; Zahn et al. 2011; Sobacchi & Mesinger 2013b), pure analytic models (Haiman et al. 1996; Thoul & Weibnberg 1996; Dijkstra et al. 2004; Furlanetto et al. 2004; Alvarez et al. 2012; Kaurov & Gnedin 2013), and semianalytic models (Babul & Rees 1992; Efstathiou 1992; Shapiro et al. 1994; Mesinger & Dijkstra 2008; Font et al. 2011; Wyithe & Loeb 2013), each with its pros and cons.

The treatment of the IGM itself, a composite inhomogeneous multiphase medium, is not fully accurate either. In principle, the problem is less severe for hydrodynamic simulations than for (semi)numerical and (semi)analytic models because these equations apply locally, so one must not worry about the spatially fluctuating properties of IGM. However, current simulations do not resolve the different ionized phases (Finlator et al. 2012).

A usual procedure (e.g., Shapiro et al. 1994; Wyithe & Loeb 2003; Benson et al. 2006; Zhang et al. 2007) is to consider the IGM as having a simple hydrogenic composition and constant, uniform temperature, equal to the characteristic temperature of photoionized hydrogenic gas (∼10⁴ K), and to focus on the evolution of the ionization state through the simple equation derived by Shapiro & Giroux (1987). However, the IGM temperature is crucial not only for estimating the minimum galaxy mass but also for computing the recombination coefficients, so such an approximation also affects the ionization state of the IGM.

Hui & Gnedin (1997) derived the first coupled equations for the ionization state and temperature of the IGM taking into account the dependence of the latter on hydrogen and helium abundances and local density of the gas (Miralda-Escudé & Rees 1994). However, these equations only held for the cooling phase after ionization, and Haiman & Holder (2003) and Hui & Haiman (2003) extended them to include the ionization period.

However, the IGM is also multiphasic (Miralda-Escudé et al. 2000 and references therein): the neutral, singly, and doubly ionized regions are separated. Choudhury & Ferrara (2005) derived the equations for IGM evolution since the dark ages taking into account the full composite, inhomogeneous, and multiphase nature of IGM. However, instead of taking the average recombination coefficients in each (ionized) phase they use the value these coefficients would take for the average (approximately mass-weighted) IGM temperature. On the other hand, they ignored the mass exchanges between halos and IGM, although about 90% of the initial diffuse gas ends up locked into halos, and the current IGM metallicity shows that halos also eject substantial amounts of gas into the medium.

These mass exchanges affect the volume filling factors of the various ionized species as well as the mean particle kinetic energy, so they must be taken into account. In principle, this would introduce one explicit differential equation for each of the varying comoving densities. However, taking into account their trivial form (i.e., the variation in each quantity is equal to the corresponding source function), these variations can be directly included in the usual master equations. Note that the IGM metallicity determining its mean molecular weight also changes. However, the mass fraction in metals in the IGM is so small (∼10⁻² Z_⊙ at z ∼ 5; e.g., Simcoe et al. 2011 and D'Odorico et al. 2013) that these variations have a negligible effect.

Lastly, the source functions in the IGM master equations were calculated by averaging the feedback of luminous objects over ionized regions, assuming an evolving universal halo mass function (MF). Yet, as the mass of halos able to trap gas and to form stars depends on the temperature and ionization state of the surrounding IGM, the halo MF itself depends on the environment. That is, the MF of halos lying in ionized or neutral regions differs. This bias, hereafter referred to as the ionization-bias to distinguish it from the well-known mass-bias (e.g., Tinker et al. 2010 and references therein)¹ must thus be corrected for.

The aim of the present paper is to improve the analytical treatment of IGM evolution by deriving new, more accurate, master equations for its ionization state and temperature, and by estimating the halo ionization-bias necessary to properly compute the source functions in these equations. Such an improved treatment of IGM can be incorporated into any given (semi)numerical or (semi)analytic model of galaxy formation such as the one developed by Manrique et al. (2015). The IGM properties shown throughout the paper to illustrate the effects of the new treatment have been obtained from that model.

In Sections 2 and 3, we derive the new equations for the IGM ionization state and temperature, respectively. In Section 4, we derive the halo MFs that result in neutral and ionized environments. Our results are discussed and summarized in Section 5.

The structure of IGM is determined by the ionizing radiation from luminous sources. UV photons with a short mean free path ionize small regions around these sources. Their less energetic fraction gives rise to singly ionized hydrogen and helium bubbles, while the less abundant, more energetic fraction gives rise to doubly ionized helium subbubbles. Bubbles and subbubbles grow and progressively overlap or retract and fragment, depending on the intensity of the ionizing flux is. In any case, the neutral, singly and doubly ionized phases are kept well separated at any time.

As mentioned, IGM is not only multiphasic but also inhomogeneous. All IGM properties, such as temperature, baryon density or H i number density, are random fields characterized by their respective probability distribution functions (PDFs). We are here interested in the time evolution of the IGM properties averaged over different regions. When these averages refer to the neutral, singly, and doubly ionized phases, they will be denoted by angular brackets with subscripts I, II, and III, respectively; when they refer to regions encompassing one particular chemical species, such as H ii (i.e., all ionized regions), the subscript will explicitly indicate that chemical species; and when the average is over the entire IGM, there will be no subscript. Averages of the product of several (either correlated or uncorrelated) quantities are for their joint PDF, so they will differ in general from the product of the averages of the individual quantities.

The local comoving density of H ii ions, ${{n}_{{\rm HII}}}$, at the cosmic time t satisfies the balance equation

$$\begin{eqnarray}&&\frac{d{{n}_{{\rm HII}}}}{dt}={{\dot{N}}_{{\rm HII}}}-\frac{{{\alpha }_{{\rm HI}}}(T)}{{{a}^{3}}}{{n}_{{\rm HII}}}\;{{n}_{{\rm e}}},\end{eqnarray} \tag{ 1 }$$

where ${{n}_{{\rm e}}}$ is the comoving density of free electrons, ${{\dot{N}}_{{\rm HII}}}$ is the local metagalactic emissivity of H i-ionizing photons due to luminous sources and recombinations, including redshifted photons emitted and not absorbed at higher z's, and the second term on the right is the recombination rate density to H i. Note that the temperature-dependent recombination coefficient for optically thin regions, ${{\alpha }_{{\rm HI}}}(T)$ (see, e.g., Faucher-Giguère et al. 2009; Meiksin 2009), is divided by the cube of the cosmic scale factor a(t) so as to express it in comoving units.

Taking the average of Equation (1) over the whole IGM, with the average of the second term on the right decomposed in the sum of the averages over the different phases I, II, and III, duly weighted by their respective volume filling factors, ${{Q}_{{\rm I}}}=1-{{Q}_{{\rm HII}}}$, ${{Q}_{{\rm II}}}={{Q}_{{\rm HII}}}-{{Q}_{{\rm HeIII}}}$, and ${{Q}_{{\rm III}}}={{Q}_{{\rm HeIII}}}$, with ${{Q}_{{\rm HII}}}$ and ${{Q}_{{\rm HeIII}}}$ standing for the H ii and He iii volume filling factors, respectively defined as $\langle {{n}_{{\rm HII}}}\rangle /\langle {{n}_{{\rm H}}}\rangle$ and $\langle {{n}_{{\rm HeIII}}}\rangle /\langle {{n}_{{\rm He}}}\rangle$, gives rise to the rigorous equation

$$\begin{eqnarray}&&\frac{d\langle {{n}_{{\rm HII}}}\rangle }{dt}=\langle {{\dot{N}}_{{\rm HII}}}\rangle -\mathop{\sum }\limits_{i={\rm II}}^{{\rm III}}{{\left\langle \frac{{{\alpha }_{{\rm HI}}}(T)}{{{a}^{3}}}{{n}_{{\rm HII}}}\;{{n}_{{\rm e}}} \right\rangle }_{i}}{{Q}_{i}}.\end{eqnarray} \tag{ 2 }$$

Approximating ${{\alpha }_{{\rm HI}}}(T)$ in ionized regions by a uniform value corresponding to the characteristic temperature ${{T}_{{\rm typ}}}$ of photoionized hydrogenic gas (∼10⁴ K), and dividing by the approximately constant value (ignoring inflows and outflows) of the mean comoving hydrogen density, $\langle {{n}_{{\rm H}}}\rangle$, we arrive at the following simple equation for the H ii volume filling factor ${{Q}_{{\rm HII}}}$ (Shapiro & Giroux 1987),

$$\begin{eqnarray}&&\frac{d{{Q}_{{\rm HII}}}}{dt}=\frac{\langle {{{\dot{N}}}_{{\rm HII}}}\rangle }{\langle {{n}_{{\rm H}}}\rangle }-\frac{{{\alpha }_{{\rm HI}}}({{T}_{{\rm typ}}})}{{{a}^{3}}}{{C}_{{\rm HII}}}\;{{\langle {{n}_{{\rm e}}}\rangle }_{{\rm HII}}}\;{{Q}_{{\rm HII}}},\end{eqnarray} \tag{ 3 }$$

where ${{C}_{{\rm HII}}}\equiv \langle n_{{\rm HII}}^{2}\rangle _{{\rm HII}}^{{}}/\langle {{n}_{{\rm HII}}}\rangle _{{\rm HII}}^{2}$ is the so-called clumping factor. To write Equation (3), we have made two approximations: $\langle {{n}_{{\rm e}}}\rangle _{{\rm HII}}^{{}}\approx \langle {{n}_{{\rm HII}}}\rangle _{{\rm HII}}^{{}}$ and $\langle {{n}_{{\rm H}}}\rangle _{{\rm HII}}^{{}}\approx \langle {{n}_{{\rm H}}}\rangle$. The former presumes hydrogenic composition, and the latter presumes that ionized regions have the same average properties as the whole IGM.

However, $\langle {{n}_{{\rm H}}}\rangle$ is not constant, but evolves due to inflows and outflows into and from halos. In addition, the IGM is not strictly hydrogenic, as its temperature varies both in space and time. Lastly, there should be, as mentioned earlier, some halo ionization-bias, so the average IGM properties in ionized regions should differ in general from the global average properties. We should thus try to do better.

Let us come back to the rigorous Equation (2). Neglecting metals, we have ${{n}_{{\rm e}}}={{n}_{{\rm HII}}}(1+Y/4X)+{{n}_{{\rm HeIII}}}$, where the comoving density of He iii ions, ${{n}_{{\rm HeIII}}}$, takes the approximate form $f(X,Y,\gamma ){{n}_{{\rm HII}}}$, with f equal to a universal function of the hydrogen and helium mass fractions, X and Y, respectively, and the typical spectral index γ of ionizing sources. Thus, the average in the summation on the right of Equation (2) splits into a sum of two products of the form: average of a function of T times average of $n_{{\rm HII}}^{2}$. This is possible thanks to the fact that there is essentially no correlation between T and ${{n}_{{\rm HII}}}$. The reason for this is that, in ionized regions, ${{n}_{{\rm HII}}}$ is essentially equal to ${{n}_{{\rm H}}}=X\;n$, where $n$ is the baryon density. Furthermore, the only terms in Equation (7) for the evolution of the IGM temperature coupling n and T are the second and fifth ones giving the heating/cooling by adiabatic compression/expansion of the fluid element, and the heating/cooling by the loss/gain of baryons due to inflows/outflows, respectively, which are less than the first term giving the cosmic adiabatic cooling, and much less than the third and fourth terms including the stochastic effects of nearby luminous sources. Under these justified approximations, Equation (2) becomes

$$\begin{eqnarray}&&\frac{d\langle {{n}_{{\rm HII}}}\rangle }{dt}=\langle {{\dot{N}}_{{\rm HII}}}\rangle -{{\left\langle \frac{{{\alpha }_{{\rm HI}}}(T)}{X{{\mu }^{{\rm e}}}{{a}^{3}}} \right\rangle }_{{\rm HII}}}{{C}_{{\rm HII}}}\;{{\langle {{n}_{{\rm H}}}\rangle }^{2}}\;{{Q}_{{\rm HII}}},\end{eqnarray} \tag{ 4 }$$

where ${{\mu }^{{\rm e}}}$ is the electronic contribution to the mean molecular weight. Then, dividing Equation (4) by $\langle {{n}_{{\rm H}}}\rangle$, we arrive at the new equation

$$\begin{eqnarray}\begin{array}{rcl} \frac{d{{Q}_{{\rm HII}}}}{dt} & = & \frac{\langle {{{\dot{N}}}_{{\rm HII}}}\rangle }{\langle {{n}_{{\rm H}}}\rangle }-\left[ {{\left\langle \frac{{{\alpha }_{{\rm HI}}}(T)}{{{\mu }^{{\rm e}}}{{a}^{3}}} \right\rangle }_{{\rm HII}}}{{C}_{{\rm HII}}}\;\langle n\rangle \right. \\ {} & {} & +\left. \frac{d\;{\rm ln} \langle {{n}_{{\rm H}}}\rangle }{dt} \right]{{Q}_{{\rm HII}}}. \\ \end{array}\end{eqnarray} \tag{ 5 }$$

Moreover, taking the Taylor expansion around the average temperature in phase i, ${{\langle T\rangle }_{i}}$, of the function of temperature f(T) given by the first term in claudators on the right-hand side of Equation (5), we find that the average over ionized regions, i = II + III, of f(T) is well-approximated by $f({{\langle T\rangle }_{i}})+{{({{d}^{2}}\;f/d{{T}^{2}})}_{{{\langle T\rangle }_{i}}}}\sigma {{_{T}^{2}}_{i}}/2$, where ${{\sigma }_{T}}_{i}$ is the dispersion in temperatures around the mean.

Besides being better justified than expression (3), expression (5) is also more accurate for the following reasons: (1) instead of taking the recombination rate density at a fixed typical temperature divided by ${{\langle {{\mu }^{{\rm e}}}\rangle }_{{\rm HII}}}$, it uses the z-dependent average of ${{\alpha }_{{\rm HI}}}(T)/{{\mu }^{{\rm e}}}$ in the H ii region, and (2) the last term on the right accounts for the changing comoving hydrogen density due to inflows/outflows. In Figure 1, we compare ${{\langle {{\alpha }_{{\rm HI}}}(T)/{{\mu }^{{\rm e}}}\rangle }_{{\rm HII}}}$ for the z-dependent temperature shown in Figure 2 to the uniform constant value ${{\alpha }_{{\rm HI}}}({{T}_{{\rm typ}}})/{{\langle {{\mu }^{{\rm e}}}\rangle }_{{\rm HII}}}$ with ${{T}_{{\rm typ}}}={{10}^{4}}$ K appearing in Equation (3). As can be seen, the difference is noticeable, particularly around the redshifts z = 10.3 and 5.5 of complete ionization in the particular galaxy model with double reionization considered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When ${{Q}_{{\rm HII}}}$ reaches the value of one and $\langle {{\dot{N}}_{{\rm HII}}}\rangle$ is sufficient to balance recombinations, a period of ionization equilibrium begins in which ${{Q}_{{\rm HII}}}$ stays equal to one. However, if $\langle {{\dot{N}}_{{\rm HII}}}\rangle$ becomes insufficient to keep ionized regions growing (or stable), a recombination will begin. The constant and decreasing values of ${{Q}_{{\rm HII}}}$ in those two regimes are also governed by Equation (5), in the former case with $\langle {{\dot{N}}_{{\rm HII}}}\rangle$ replaced by the equilibrium value, with the leftover metagalactic emissivity eventually used, duly redshifted, to ionize more hydrogen atoms at lower z's.

A similar derivation leads to the homologous equation for the He iii volume filling factor,

$$\begin{eqnarray}\begin{array}{rcl} \frac{d{{Q}_{{\rm HeIII}}}}{dt} & = & \frac{\langle {{{\dot{N}}}_{{\rm HeIII}}}\rangle }{\langle {{n}_{{\rm He}}}\rangle }-\left[ {{\left\langle \frac{{{\alpha }_{{\rm HeII}}}(T)}{{{\mu }^{{\rm e}}}{{a}^{3}}} \right\rangle }_{{\rm III}}} \right. \\ {} & {} & \times {{C}_{{\rm HII}}}\;\langle n\rangle +\left. \frac{d\;{\rm ln} \langle {{n}_{{\rm He}}}\rangle }{dt} \right]{{Q}_{{\rm HeIII}}}. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

Again, if at any point a period of He iii ionization equilibrium or recombination takes place, then ${{Q}_{{\rm HeIII}}}$ stays equal to one or begins to diminish, respectively, according to the same Equation (6).

Photo-ionization leads to photo-heating of the different IGM phases. Other heating mechanisms acting on the IGM are Compton heating by X-rays and by cosmic microwave background (CMB) photons at very high-z (after decoupling of baryons from radiation at z ∼ 150). Such heating is partially balanced by the cooling due to recombinations and deexcitations, cosmic expansion, Comptonization from CMB photons at low z, and collisional cooling (significant only in very hot neutral regions, if any). In addition, density fluctuations suffer gravitational contraction/expansion causing extra heating/cooling. These are the main mechanisms causing the thermal evolution of the IGM. Below we mention (in italics) a few additional mechanisms that are included in the present more accurate treatment (see also Hui & Gnedin 1997 for other possible heating and ionizing mechanisms, due to decaying or annihilating dark matter, not included herein).

The local temperature of the IGM evolves according to the differential equation (e.g., Choudhury & Ferrara 2005)

$$\begin{eqnarray}\begin{array}{rcl} \frac{dT}{d\;{\rm ln} (1+z)} & = & T\left[ 2+\frac{2}{3}\frac{d\;{\rm ln} (n/{{\langle n\rangle }_{i}})}{d\;{\rm ln} (1+z)} \right. \\ {} & {} & +\frac{d\;{\rm ln} \mu }{d\;{\rm ln} (1+z)}+\frac{d\;{\rm ln} }{d\;{\rm ln} (1+z)} \\ {} & {} & -\left. \frac{d\;{\rm ln} n}{d\;{\rm ln} (1+z)} \right]. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

The first term in claudators on the right, equal to 2, gives the cosmological adiabatic cooling of the gas element; the second term gives its adiabatic heating/cooling by gravitational compression/expansion for the baryon density $n$ around the mean value ${{\langle n\rangle }_{i}}$ in region i, taking into account that most diffuse IGM is in a linear or moderately nonlinear regime; the third term gives the cooling due to the increase in mean molecular weight, μ, caused by ionization and outflows from halos; the fourth term gives the Compton cooling from CMB photons, and the gain/loss of energy density, ε, due to photo-ionization/recombination, Compton heating from X-rays, the achievement of energy equipartition by newly ionized/recombined fraction of gas (the different phases have distinct temperatures in general) plus mechanical heating accompanying outflows from halos; and the fifth term gives the cooling/heating by the gain/loss of baryon density, $n$, due to outflows/inflows (this changes the average specific energy of the IGM). As outflows take place from halos harboring luminous sources, we assume that they only affect ionized regions.

Multiplying Equation (7) by μ n, and taking the average over each specific phase under the approximation, for the reasons mentioned in Section 2, that μ, ε, n, and T do not correlate with each other, we arrive at

$$\begin{eqnarray}&&\frac{d\;{\rm ln} {{\langle T\rangle }_{i}}}{d\;{\rm ln} (1+z)}=2+\frac{d\;{\rm ln} ({{\langle \mu \rangle }_{i}}{{\langle \rangle }_{i}}/{{\langle n\rangle }_{i}})}{d\;{\rm ln} (1+z)},\end{eqnarray} \tag{ 8 }$$

with i = I, II, or III. Note that, in neutral regions (i = I), there are no stochastic effects of luminous sources: ε does not change either through photo-ionization or by X-rays, μ is kept strictly equal to the primordial value, and there is only a small change in n due to inflows. Consequently, a strong correlation is foreseen between the quantities ε and n and temperature. Yet, we still ignore such a correlation for simplicity. This approximation is only necessary during the initial period of increasing ionization; in recombination periods, the gas properties in the new neutral phase remain uncorrelated as they have suffered important stochastic feedback effects from luminous objects over the previous ionized phase.

And what about the temperature dispersion around the mean in the different IGM phases, also required in Equations (5) and (6)? To calculate $\sigma {{_{T}^{2}}_{i}}={{\langle {{T}^{2}}\rangle }_{i}}-\langle T\rangle _{i}^{2}$ we need to consider the relation

$$\begin{eqnarray}&&\frac{\frac{1}{2}d{{T}^{2}}}{d\;{\rm ln} (1+z)}={{T}^{2}}\left[ 2+\frac{2}{3}\frac{d\;{\rm ln} (n/{{\langle n\rangle }_{i}})}{d\;{\rm ln} (1+z)}+\frac{d\;{\rm ln} (\mu /n)}{d\;{\rm ln} (1+z)} \right]\end{eqnarray} \tag{ 9 }$$

following from Equation (7). The same steps above lead to

$$\begin{eqnarray}&&\frac{d\;{\rm ln} \langle {{T}^{2}}\rangle _{i}^{1/2}}{d\;{\rm ln} (1+z)}=2+\frac{d\;{\rm ln} ({{\langle \mu \rangle }_{i}}{{\langle \rangle }_{i}}/{{\langle n\rangle }_{i}})}{d\;{\rm ln} (1+z)}=\frac{d\;{\rm ln} {{\langle T\rangle }_{i}}}{d\;{\rm ln} (1+z)}.\end{eqnarray} \tag{ 10 }$$

The initial conditions for Equation (8) are ${{\langle T\rangle }_{{\rm I}}}({{z}_{{\rm ini}}})={{T}_{{\rm CMB}}}({{z}_{{\rm ini}}})$ and ${{\langle T\rangle }_{{\rm II}}}({{z}_{{\rm ini}}})={{\langle \mu \rangle }_{{\rm II}}}{{\langle T\rangle }_{{\rm I}}}({{z}_{{\rm ini}}})/{{\langle \mu \rangle }_{{\rm I}}}$, where ${{z}_{{\rm ini}}}$ is the redshift at which the IGM temperature begins to deviate from the temperature of CMB photons,² satisfying $1+{{z}_{{\rm ini}}}=100{{({{{\Omega }}_{{\rm b}}}{{h}^{2}}/0.0125)}^{2/5}}$, where ${{{\Omega }}_{{\rm b}}}$ and h are the baryonic density parameter and the Hubble parameter scaled to 100 km s⁻¹ Mpc⁻¹. Similarly, the initial condition for Equation (10) is ${{\langle {{T}^{2}}\rangle }_{{\rm II}}}({{z}_{{\rm ini}}})=\langle \mu \rangle _{{\rm II}}^{2}{{\langle {{T}^{2}}\rangle }_{{\rm I}}}({{z}_{{\rm ini}}})/\langle \mu \rangle _{{\rm I}}^{2}$, where ${{\langle {{T}^{2}}\rangle }_{{\rm I}}}({{z}_{{\rm ini}}})$ is equal to $\sigma {{_{T}^{2}}_{{\rm CMB}}}({{R}_{{\rm J}}},{{z}_{{\rm ini}}})+T_{{\rm CMB}}^{2}({{z}_{{\rm ini}}})$, where $\sigma {{_{T}^{2}}_{{\rm CMB}}}({{z}_{{\rm ini}}})=2\;T_{{\rm CMB}}^{2}({{z}_{{\rm ini}}}){{\sigma }_{0}}{{({{R}_{{\rm J}}},{{z}_{{\rm ini}}})}^{2}}/3$ is the CMB temperature variance at the Jeans scale at recombination, ${{R}_{{\rm J}}}$, evolved to ${{z}_{{\rm ini}}}$, with σ₀(R, z) the zero-order spectral moment at the scale R and redshift z. The reason for the filtering at the scale ${{R}_{{\rm J}}}$ is that, at smaller scales, there were no temperature fluctuations at recombination, and the uniform temperature on those scales only suffered cosmological adiabatic cooling and the effects of luminous sources, uncorrelated with T. If there is a period of increasing recombination, the initial mean temperature and variance in the recombined region are equal (except for different mean molecular weights) to those in the ionized phase giving it rise. We have checked that $\sigma _{T}^{2}({{\langle T\rangle }_{i}})$ is always much less than $\langle T\rangle _{i}^{2}$, meaning that the second order Taylor expansion around ${{\langle T\rangle }_{i}}$ of any arbitrary function f of temperature is really close to the value f(T).

In Figure 2, we show the temperature evolution that results from the present improved treatment of IGM evolution for the source functions, ${{\dot{N}}_{{\rm HI}}}$, ${{\dot{N}}_{{\rm HeIII}}}$, $d\;{\rm ln} {{\langle \rangle }_{i}}/d\;{\rm ln} (1+z)$ and $d\;{\rm ln} {{\langle n\rangle }_{i}}d\;{\rm ln} (1+z)$, provided by the same galaxy model as in Figure 1.

The chance that halos with a given mass M at t will trap gas, and that the trapped gas will cool either through molecular bands or atomic lines and form metal-poor or metal-rich stars, respectively, depends on the temperature and ionization state of the IGM in which the halos are embedded. Consequently, the halo MF itself must vary between neutral and ionized environments. Note that, given the homogeneity of the universe, these probabilities are not a function of a specific point. In particular, the probability that a given arbitrary point lies in a ionized or neutral region is uniform and equal to ${{Q}_{{\rm HII}}}(t)$ and $1-{{Q}_{{\rm HII}}}(t)$, respectively.

To calculate the probability that a halo with mass M is located in an ionized region at the cosmic time t, P_M(H ii, t), we will first consider the conditional probabilities ${{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm II},{{t}_{{\rm f}}})$ and ${{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm I},{{t}_{{\rm f}}})$ that the halo is in an ionized region at t given that it was either in an ionized or neutral region, respectively, at its formation at ${{t}_{{\rm f}}}$. The former of these two quantities is simply

$$\begin{eqnarray}&&{{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm II},{{t}_{{\rm f}}})=1-{{P}_{{\rm HI}}}(t,{{t}_{{\rm f}}}),\end{eqnarray} \tag{ 11 }$$

where ${{P}_{{\rm HI}}}(t,{{t}_{{\rm f}}})$ is the probability that the halo environment recombines between ${{t}_{{\rm f}}}$ and t because of the absence of nearby sufficiently luminous sources. The latter is given by

$$\begin{eqnarray}&&{{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm I},{{t}_{{\rm f}}})=P_{M}^{\star }(t,{{t}_{{\rm f}}})+{{P}_{{\rm HII}}}(t,{{t}_{{\rm f}}}),\end{eqnarray} \tag{ 12 }$$

where $P_{M}^{\star }(t,{{t}_{{\rm f}}})$ is the probability that star formation begins to take place in a halo with M lying in a neutral environment between ${{t}_{{\rm f}}}$ and t (we say "begins" because newborn stars soon ionize the medium around the halo), and ${{P}_{{\rm HII}}}(t,{{t}_{{\rm f}}})$ is the probability that the halo environment will become ionized in the same period of time because of the presence of nearby external ionizing sources.

To derive Equations (11) and (12) we have assumed that the probabilities ${{P}_{{\rm HI}}}(t,{{t}_{{\rm f}}})$ and ${{P}_{{\rm HII}}}(t,{{t}_{{\rm f}}})$ are independent of halo mass. This may not be the case if there is some correlation between the halo mass- and ionization-biases. However, in terms of the effect of density on the ionization state of a region the tendency for halos harboring more powerful ionizing sources to lie in higher-density regions contrasts with that for ionized bubbles to stretch more rapidly in lower-density regions, so they tend to balance one another. Therefore, even though the importance of this correlation is hard to assess without performing accurate hydrodynamic simulations with ionizing radiative transfer, we do not expect it to be too marked. In other words, the present treatment should be reasonably approximate.

The total probability of finding a halo ionized at t can be expressed in terms of the above conditional probabilities and ${{P}_{M}}({\rm H}\;{\rm II},t)$ upon formation,

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{M}}({\rm H}\;{\rm II},t) & = & {{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm II},{{t}_{{\rm f}}}){{P}_{M}}({\rm H}\;{\rm II},{{t}_{{\rm f}}}) \\ {} & {} & +{{P}_{M}}({\rm H}\;{\rm II},t|{\rm H}\;{\rm I},{{t}_{{\rm f}}})\left[ 1-{{P}_{M}}({\rm H}\;{\rm II},{{t}_{{\rm f}}}) \right]. \\ \end{array}\end{eqnarray} \tag{ 13 }$$

Substituting the conditional probabilities on the right of Equation (13) by Expressions (11) and (12), setting $t={{t}_{{\rm f}}}+{\Delta }t$, and taking the limit of small Δt, Equation (13) leads to the following differential equation governing the evolution of ${{P}_{M}}({\rm H}\;{\rm II},t)$

$$\begin{eqnarray}&&\frac{d{{P}_{M}}({\rm H}\;{\rm II},t)}{dt}=\frac{d{{Q}_{{\rm HII}}}(t)}{dt}+\frac{dP_{M}^{\star }(t,{{t}_{{\rm f}}})}{dt}\left[ 1-{{Q}_{{\rm HII}}}(t) \right].\end{eqnarray} \tag{ 14 }$$

To derive Equation (14) we have taken ${{P}_{{\rm HI}}}(t,{{t}_{{\rm f}}})=0$ and ${{P}_{{\rm HII}}}(t,{{t}_{{\rm f}}})=[{{Q}_{{\rm HII}}}(t)-{{Q}_{{\rm HII}}}({{t}_{{\rm f}}})]/[1-{{Q}_{{\rm HII}}}({{t}_{{\rm f}}})]$ in periods of increasing ionization, and ${{P}_{{\rm HI}}}(t,{{t}_{{\rm f}}})=[{{Q}_{{\rm HII}}}({{t}_{{\rm f}}})-{{Q}_{{\rm HII}}}(t)]/{{Q}_{{\rm HII}}}({{t}_{{\rm f}}})$ and ${{P}_{{\rm HII}}}(t,{{t}_{{\rm f}}})=0$, in periods of increasing recombination. Interestingly, in both cases one is led to the same differential Equation (14),whose solution for the initial condition ${{P}_{M}}({\rm H}\;{\rm II},0)=0$ yields the desired probability ${{P}_{M}}({\rm H}\;{\rm II},t)$ of finding a halo with M in a ionized region at t, its complementary value giving the probability ${{P}_{M}}({\rm H}\;{\rm I},t)$ of finding it in a neutral region.

The probability $P_{M}^{\star }(t,{{t}_{{\rm f}}})$ in Equation (14) is hard to estimate analytically because it depends on the number fraction of H₂ molecules, ${{f}_{{{{\rm H}}_{2}}}}$, at the center of halos with M, whose PDF cannot be established without making an appeal to the whole halo aggregation history. Thus, this function must be drawn from a full treatment of galaxy and IGM evolution. In Figure 3, we plot this function obtained from the same galaxy model as in previous Figures. The halo MFs in ionized and neutral regions resulting from a global MF of the Sheth & Tormen (2002) form at two different redshifts are plotted in Figure 4. As can be seen, the higher the redshift, the more marked the effect,which is only visible, of course, before full ionization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the present paper, we have derived an improved version of the master equations for the evolution of IGM ionization state and temperature, accounting for the composite, inhomogeneous, multiphase nature of this medium. Besides all the usual effects, the new version includes collisional cooling in hot neutral regions (necessary to deal with recombination periods as found in double reionization), mass exchanges between halos and IGM, and the achievement equipartition for newly ionized/recombined gas. In addition, we have derived the probability that a halo with a given mass M at z is located in an ionized or neutral environment, which is needed to accurately compute the source functions required in the IGM master equations.

To check the performance of this improved treatment of IGM we coupled it to the galaxy model by Manrique et al. (2015) for realistic values of the parameters leading to double reionization (E. Salvador-Solé & A. Manrique 2015, in progress). The main results were as follows.

• 1.  
  The average temperatures in the three IGM phases show marked variations over the different ionization/recombination periods. This harbors relevant information on the epoch of reionization. The usual treatment dealing with the average temperature over the whole IGM (or at mean IGM density, T₀) loses this information.
• 2.  
  The inclusion of collisional cooling is mandatory to recover the sudden decrement in the average temperature of neutral regions after first ionization in double reionization (see Figure 2). In the only work to date, by Choudhury & Ferrara (2005), dealing with the evolution of the average temperature in the different IGM phases, neutral regions cooled adiabatically after decoupling.
• 3.  
  The average temperatures of singly and doubly ionized regions show a maximum similar to that found by Choudhury & Ferrara (2005; see panel f of their Figure 1). However, our temperatures also show a minimum, due to the recombination after first ionization. More importantly, the average temperature in doubly ionized regions is always higher than in singly ionized ones, while this was surprisingly not the case in Choudhury & Ferrara's solution.
• 4.  
  Although the average temperature in singly ionized regions is not as high as that reported by Choudhury and Ferrara, it is still notably higher (by a factor of ∼3) than the value of 10⁴ K often adopted in reionization studies (Shapiro et al. 1994; Wyithe & Loeb 2003; Benson et al. 2006; Zhang et al. 2007).
• 5.  
  This difference translates into the average recombination coefficients. The values we find are substantially smaller (by a factor of ∼4) than found for the temperature of 10¹⁰ K, and somewhat greater (by a factor of ∼2) than the minimum value at the average temperature reached in Choudhury & Ferrara's solution.
• 6.  
  This affects the evolution of the volume filling factors of ionized hydrogen and helium for identical source functions (identical galaxy models). However, this makes a small difference compared to that arising from the galaxy models used, which may lead, for instance, to single or double reionization.
• 7.  
  We have computed the halo ionization-bias in the calculation of the source functions appearing in the IGM master equations. The ratios between the halo MF in ionized and all environments found for low mass newly star-forming halos and for the rest are, respectively, equal to ∼3 × 10⁻⁴ (0.1) and ∼0.3 (0.4) at z = 30 (z = 15).

This improved treatment of IGM can be easily implemented in any model of galaxy and IGM evolution. This is particularly advisable for accurate models of galaxy formation or reionization when contrasting them with current observations (e.g., E. Salvador-Solé & A. Manrique 2015, in progress) or future ones (e.g., 21 cm line experiments).

This work was supported by the Spanish DGES grant AYA2012-39168-C03-02, and the Catalan DIUE grant 2009SGR00217.