Analytic Formulation of 21 cm Signal from Cosmic Dawn: Lyα Fluctuations

The hydrogen-ionizing (ultraviolet—UV) photons emitted from the star within the collapsed structures are absorbed in the immediate vicinity of the sources and carve out H ii regions around them in the IGM. We use excursion set formalism to compute a size distribution of the ionization bubbles by defining self-ionized regions (Furlanetto et al. 2004b). Such regions have enough sources to ionize all gas within them. Their sizes are determined by the ionization efficiency factor, defined as $\zeta =1/{f}_{\mathrm{coll}}={M}_{\mathrm{tot}}/{M}_{\mathrm{coll}}$. Here, f_coll is the fraction of collapsed mass inside the self-ionized region, and ζ is a function of the properties of the sources as well as the surrounding halo (Furlanetto et al. 2004b):

$$\begin{eqnarray}&&\zeta ={f}_{\star }{f}_{\mathrm{esc}}{N}_{\mathrm{ion}}{N}_{\mathrm{rec}}^{-1}.\end{eqnarray} \tag{ 5 }$$

Here, f_⋆ is the fraction of collapsed baryons that are converted into stars, f_esc is the fraction of ionizing photons that escape the source halo, N_ion is the number of UV photons created per stellar baryon, and N_rec is the number of recombinations. We assume ζ to be independent of redshift in this paper, even though it can evolve with time owing to the evolution of quantities used to define it. For higher values of ζ, the reionization is completed at higher redshift (e.g., RS18). These self-ionized regions are larger than the H ii regions of a single source because they are created by highly clustered multiple sources in the early universe (for the ΛCDM model). For our work, we assume the region to be spherical.

Photons of energy E ≫ 13.6 eV (X-rays) escape the H ii region into the surrounding medium. These photons ionize the neutral gas up to 10% and can impart up to 20% of their energy into heating the gas through photoionization and secondary collisional ionization and excitation (e.g., Shull & van Steenberg 1985; Venkatesan et al. 2001). In our study, we have assumed the medium outside the ionized region to be completely neutral (n = 1) because the fraction of ionization due to this process is generally small.

The photoionization cross section of X-ray photons falls as ${\sigma }_{i}{(\nu )={{\sigma }_{i}}_{0}(\nu /{\nu }_{i})}^{-3}$, with ν_i being the ionization threshold frequency of species i. In this work, we only consider two species: neutral hydrogen and neutral helium, with their relative fractions x_i = 12/13 and 1/13 of the baryon number, respectively. This is a valid assumption because the metallicity in the early universe is very low. Because the low-energy X-ray photons are absorbed with higher probability, they contribute to heating the medium immediately surrounding the H ii region, whereas the higher energy photons free-stream through the medium and might get absorbed far away from any source. These photons uniformly heat up the whole IGM to some background temperature T_bg.

We assume the X-ray photon source luminosity to be given by a power law (e.g., Mesinger et al. 2011 and references therein), ${\dot{N}}_{\nu }={\dot{N}}_{t}{(\nu /{\nu }_{\min })}^{-\alpha }$, where ν_min is the lowest frequency of X-ray photons escaping from source halos, and N_heat is the number of X-ray photons emitted per stellar baryon. We assume that f_H = 0.15 is the fraction of energy of emitted photoelectrons that goes into heating the medium (Shull & van Steenberg 1985; Venkatesan et al. 2001). Other than the adiabatic expansion of the universe, we neglect all other cooling processes.

We divide the neutral hydrogen volume into two zones. In the near zone, the heating is dominated by X-ray photons from an individual self-ionized region. In the far zone, the contribution from all of the faraway background sources is taken into account. For more details and explanations, see RS18, where the increase in temperature due to a self-ionized region of radius R_x at a distance R₀ from the center of the ionized region was calculated at redshift ${z}_{c}^{{\prime} }$:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}T^{\prime} & = & \displaystyle \frac{{{hf}}_{{\rm{H}}}\alpha {N}_{\mathrm{heat}}{f}_{\star }{n}_{0}{\nu }_{\min }^{\alpha }}{3{k}_{B}\zeta }\displaystyle \frac{{R}_{x}^{3}}{{R}_{0}^{2}}{\left(1+{z}_{c}^{{\prime} }\right)}^{2}\\ & & {\displaystyle \int }_{t({z}_{\star })}^{t({z}_{0})}{dt}^{\prime} \displaystyle \frac{{f}_{i}(t^{\prime} )}{{f}_{i}(t)}\displaystyle \frac{{\dot{f}}_{\mathrm{coll},{\rm{g}}}}{{f}_{\mathrm{coll},{\rm{g}}}}{\left(\displaystyle \frac{1+z^{\prime} }{1+z}\right)}^{\alpha +1}\\ & & {\displaystyle \int }_{{\nu }_{\min }^{{\prime} }}^{\infty }d\nu ^{\prime} \nu {{\prime} }^{-\alpha -4}{e}^{-\tau ({R}_{0},\nu ^{\prime} )}\displaystyle \sum _{i}(\nu ^{\prime} -{\nu }_{i}){x}_{i}{{\sigma }_{i}}_{0}{{\nu }_{i}}^{3}.\end{array}\end{eqnarray} \tag{ 6 }$$

The primed quantities are calculated at the receiving point (point P), unprimed quantities are at the source (point S), and quantities with 0 subscript are comoving quantities. Here, α is the X-ray spectra power-law index, and f_i are global ionization fractions.

For EoR studies, all of the radiation between Lyα and the Lyman limit is referred to as Lyα radiation emitted from the source, and we shall follow this convention. The Lyα contribution at any point arises from two main factors: Lyα emitted from the sources and Lyα created by X-ray photoelectrons (Venkatesan et al. 2001). The latter is generally small, and we neglect it in this paper.

Lyα photons emitted from the sources escape the surrounding H ii region and redshift until their frequency nearly equals the resonant frequency of one of the Lyman series lines. When a photon redshifts to the frequency corresponding to one of the Lyman series lines, it gets scattered by the neutral hydrogen³ and eventually cascades to a frequency corresponding to Lyα. Given the complicated frequency structure of Lyman series lines, these photons are absorbed at varying distances from the source. Thus, the coupling depends on two factors: the region of influence of the Lyα radiation and the coupling coefficient y_α. Photons between Lyα and Lyβ frequencies can be absorbed by H₂ molecules, but we ignore this effect as the density of H₂ is very low in the IGM.

We define the Ly-n influence region as the distance traveled by the Ly-(n + 1) photons to redshift to Ly-n frequency. If these photons were emitted at z = z_e and absorbed at z = z_a with ν_e = ν_n + 1 and ν_a = ν_n, then the comoving distance traveled by the photon before it is absorbed in an expanding universe is (n ≥ 2)

$$\begin{eqnarray}&&{R}_{\max }(n)\simeq \displaystyle \frac{16040\,\,\mathrm{Mpc}}{{\left(1+{z}_{e}\right)}^{1/2}}\left[{\left(\displaystyle \frac{{n}^{3}(n+2)}{{\left(n+1\right)}^{3}(n-1)}\right)}^{1/2}-1\right].\end{eqnarray} \tag{ 7 }$$

These influence regions become smaller with increasing n. We note that R_max(2) (Lyα influence region) is much larger than the mean distance between ionization bubbles at any redshift. For ζ = 7.5, the values of mean comoving distance between bubbles for redshifts 25, 20, and 15 are 7.85 Mpc, 2.29 Mpc, and 0.96 Mpc, respectively. Therefore, Lyα regions are very large and merge very early. However, this would create homogeneous coupling to H i atoms only if the Lyα coupling coefficient, y_α, is high enough (Equation (1)); y_α is a function of Lyα photon (physical) number density ${n}_{\alpha }^{{\prime} }$ (Field 1958; Chen & Miralda-Escudé 2004):

$$\begin{eqnarray}&&{y}_{\alpha }=\displaystyle \frac{\sqrt{2}{{hc}}^{4}}{36\pi {\nu }_{\alpha }^{2}k}\sqrt{\displaystyle \frac{{m}_{{\rm{H}}}}{2{{kT}}_{k}}}\displaystyle \frac{{A}_{\alpha }}{{A}_{10}}\displaystyle \frac{{\nu }_{10}}{{\nu }_{\alpha }}\displaystyle \frac{{n}_{\alpha }^{{\prime} }}{{T}_{K}}.\end{eqnarray} \tag{ 8 }$$

To calculate the number density of Lyα photons (${n}_{\alpha }^{{\prime} }$) received at any point from ionizing sources, we use the method used to calculate X-ray heating in RS18. Assuming a flat spectrum between Lyα and the Lyman limit, the number density of Lyα photons at a comoving distance R₀ from the source is

$$\begin{eqnarray*}&&{n}_{\alpha ,\star }^{{\prime} }=\displaystyle \frac{{\dot{N}}_{\alpha }}{4\pi {{cR}}_{0}^{2}}\displaystyle \frac{2{\rm{\Delta }}{\nu }_{\alpha }}{{\nu }_{\beta }-{\nu }_{\alpha }}\displaystyle \frac{{\left(1+z^{\prime} \right)}^{3}}{1+z}.\end{eqnarray*}$$

Here, ${\rm{\Delta }}{\nu }_{\alpha }=\sqrt{2{{kT}}_{k}/{m}_{{\rm{H}}}{c}^{2}}{\nu }_{\alpha }$ is the Doppler line width. This factor arises because at the source the photons are emitted with frequencies between ν_β and ν_α, but the only frequencies that are absorbed at redshift $z^{\prime}$ are in the range of Δν_α around ν_α. The Lyα luminosity, ${\dot{N}}_{\alpha }$, can be expressed in terms of the size of the ionization halo, assuming that the Lyα luminosity scales with ionizing luminosity with a factor f_L. Using the balance between ionization and recombination in the ionizing region, we have

$$\begin{eqnarray*}&&{\dot{N}}_{\alpha }={f}_{L}\displaystyle \frac{4\pi }{3}{n}_{0}^{2}{\alpha }_{B}{{CR}}_{x}^{3}{\left(1+z\right)}^{3}.\end{eqnarray*}$$

It would be reasonable to assume that f_L > 1 because Lyα photons escape the halo more easily than ionizing photons. However, for the sake of completeness, we take 0.1 < f_L < 1000 in this paper. Combining all these, we get

$$\begin{eqnarray}&&{y}_{\alpha ,\star }{T}_{k}\simeq \displaystyle \frac{{S}_{\alpha }}{54\pi }\displaystyle \frac{{{hc}}^{2}{A}_{\alpha }{\nu }_{10}{n}_{0}^{2}{\alpha }_{B}}{{\nu }_{\alpha }^{2}{{kA}}_{10}({\nu }_{\beta }-{\nu }_{\alpha })}{f}_{L}C\displaystyle \frac{{R}_{x}^{3}}{{R}_{0}^{2}}{\left(1+z^{\prime} \right)}^{3}{\left(1+z\right)}^{2}.\end{eqnarray} \tag{ 9 }$$

Here, S_α is a correction factor of order unity as defined in Chen & Miralda-Escudé (2004), which depends on the photon spectrum around the Lyα line (Hirata 2006; Pritchard & Loeb 2012). This expression has been derived in a somewhat different manner in Pritchard & Loeb (2012).

We find contributions to y_α from higher-order Lyman transitions (e.g., from Lyγ to Lyβ) in a similar way, by counting the number of Lyn influence regions the point falls within. The effect of higher-order transitions is expected to be subdominant because, for a continuum source, the total number of photons between Lyβ to the Lyman limit is smaller than in the frequency range Lyβ and Lyα. However, these photons can have a substantial impact near an ionizing source because they are absorbed closer to the source given their smaller influence regions. If the distance of a point from the source R₀ is such that ${R}_{\max }(n+1)\lt {R}_{0}\lt {R}_{\max }(n)$, then the point in question will have n Lyα photons in its vicinity rather than one photon if only the transition between Lyα and Lyβ is considered. This means that the Lyα flux from the source center generally falls more rapidly than $1/{r}^{2}$ when this effect is taken into account. In this paper, we only consider Lyman series lines that have influence regions larger than the ionization bubble radius.⁴

The collisional coupling of spin temperature T_S and matter kinetic temperature T_K due to scattering of neutral hydrogen and electrons (Equation (1)) can play an important role at lower redshift (z ≤ 30) too. The coupling coefficient is proportional to the number density of colliding particles:

$$\begin{eqnarray}&&{y}_{c}=\displaystyle \frac{({n}_{{\rm{H}}}{k}_{10}^{H}+{n}_{e}{k}_{10}^{e})}{{A}_{10}}\displaystyle \frac{{T}_{\star }}{{T}_{K}}.\end{eqnarray} \tag{ 10 }$$

For collision rate coefficients, we use the following fits (Zygelman 2005; Pritchard & Loeb 2012):

$$\begin{eqnarray}{k}_{10}^{e}=\left\{\begin{array}{cl}\exp \left(-9.607+0.5\mathrm{log}({T}_{K})\exp \left(-\displaystyle \frac{{\left(\mathrm{log}({T}_{K})\right)}^{4.5}}{1800}\right)\right)\,\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1} & \quad {T}_{K}\leqslant {10}^{4}\,{\rm{K}}\\ {k}_{10,e}({T}_{K}={10}^{4}\,{\rm{K}}) & \quad {T}_{K}\gt {10}^{4}\,{\rm{K}}\end{array}\right.\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{l}{k}_{10}^{H}=\left\{\begin{array}{cl}3.6\times {10}^{-16}{T}_{K}^{3.640}\exp \left(\displaystyle \frac{6.035}{{T}_{K}}\right)\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1} & {T}_{K}\leqslant 10\,{\rm{K}}\\ 3.1\times {10}^{-11}{T}_{K}^{0.357}\exp \left(-\displaystyle \frac{32}{{T}_{K}}\right)\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1} & {T}_{K}\gt 10\,{\rm{K}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 12 }$$

These rates increase with temperature, which means that there is stronger collisional coupling for hotter gas than for cool gas. This effect is important if the gas was colder during the cosmic dawn, due to unknown physics: the prereionization absorption trough might be shallower instead of steeper, in spite of having a larger contrast of matter temperature from the CMB temperature. During EoR, the electron scattering is more effective near the sources where there is partial ionization and high temperature. However, since we assume the ionization fraction to be just the residual fraction outside of ionization bubbles, this effect is negligible in our work.

Our aim in this paper is to analytically model the correlation function of H i brightness temperature fluctuations from the early phase of EoR owing to scales that emerge from ionization, heating, and Lyα coupling inhomogeneities. These inhomogeneities are caused by bubbles of a given size distribution, which evolves with time. These bubbles determine the scales of correlation. The details of our modeling and main assumptions are discussed extensively in RS18. We briefly summarize them here for the sake of completeness, as only a subset of the details are given in this paper. We have assumed a topology where there are isolated, spherical self-ionized bubbles, surrounded by isotropic, smooth T_S profiles that might overlap with one another and smoothly merge with the background. Given the statistical isotropy (we neglect redshift-space distortion) and homogeneity of the process of reionization and heating, our assumption of spherical bubbles and isotropic spin temperature profiles holds even though the individual bubbles might not be spherical.

We compute ionization and spin temperature autocorrelation and ionization–spin temperature cross-correlation. We neglect, as noted above, the cross-correlation of density with ionization and spin temperature inhomogeneities because these contribute negligibly on the scales of interest. We neglect the clustering of self-ionized regions; a self-ionized region already accounts for clustering of ionizing sources at smaller scales. In RS18 we present cases that account for the clustering of self-ionized regions. This does not substantially alter our main results, but it might introduce new scales in the problem that correspond to the correlation scales of bubble centers. We note that this assumption is better for higher redshifts as the mean bubble separation is larger.

At any redshift, using excursion set formalism (Section 2.1) and the matter power spectrum given by the ΛCDM model, we generate a size distribution of self-ionized bubbles. Using Equations (1), (6), (9), and (10), we calculate the spin temperature in shells around these bubbles. The ionization volume fraction and volume fraction due to these shells is, respectively,

$$\begin{eqnarray}&&{f}_{i}=\displaystyle \sum _{{R}_{x}}\displaystyle \frac{4\pi }{3}N({R}_{x}){R}_{x}^{3}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{f}_{{hb}}=\displaystyle \sum _{{R}_{x}}\displaystyle \frac{4\pi }{3}N({R}_{x})({R}_{h}^{3}-{R}_{x}^{3}),\end{eqnarray} \tag{ 21 }$$

where R_x are radii of ionization bubbles, $N({R}_{x})$ is the number density of bubbles with ionization radius R_x, and R_h corresponds to the outer radius of the spin temperature profile for a given R_x (Figure 1). During the initial phase, the ionization bubbles and surrounding temperature profiles are nonoverlapping; we discuss the case of overlapping heating and Lyα coupled regions later in this section.

Zoom In Zoom Out Reset image size

For any bubble of size R_x, we take shells of thickness ${\rm{\Delta }}R({R}_{x},s)$ between radii R_x and R_h, having temperature $s={T}_{\mathrm{CMB}}/{T}_{S}$. A detailed description of notations followed in this paper is given in Table 1. To compute correlations, we assume two random points separated by a distance r, as shown in Figure 1. The formalism used for the computation of correlation functions is described in Appendix D.

Table 1.  Notations

Symbols	Explanation
δ	Overdensity of H i gas
n	Ionization state of H i gas: neutral point n = 1 and ionized point n = 0
s	Temperature state defined as s = TCMB/TS
ψ	Dimensionless brightness temperature: $\psi =n(1+\delta )(1-s)$
ξ	Autocorrelation of overdensity δ: $\xi =\langle {\delta }_{1}{\delta }_{2}\rangle$
ϕ	$\phi =n(1-s)$
μ	Autocorrelation of dimensionless brightness temperature ψ: $\mu =\langle {\psi }_{1}{\psi }_{2}\rangle -\langle \psi {\rangle }^{2}$
fi	Average ionized volume fraction
fn	Average neutral volume fraction
fhb	Total volume fraction due to TS profiles (without correcting for the overlaps)
fh	Average TS profile volume fraction after correcting for the overlaps
fb	Average background volume fraction
Rx	Radius of given ionization bubble
Rh	Outer radius of given TS profile: ${R}_{h}={R}_{h}({R}_{x})$
Rs	Inner radius of the shell with spin temperature ${T}_{S}={T}_{\mathrm{CMB}}/s$ around a given bubble
${\rm{\Delta }}{R}_{s}$	Thickness of the shell with spin temperature ${T}_{S}={T}_{\mathrm{CMB}}/s$ around a given bubble
$N({R}_{x})$	Number density of ionization bubbles of radius Rx

Download table as:  ASCIITypeset image

At low redshift, spin temperature profiles around ionization bubbles are very large and overlap significantly. Thus, in such cases, f_hb can be much larger than 1. In our previous paper (RS18), we had discussed a method to consistently take into account the overlap of heated regions. In this paper, our principal aim is to model the fluctuation of Lyα radiation. Unlike the heated regions that merge after the heating transition, the mean free path of lower Lyman series photons is larger than the mean interbubble distance at all times for z < 30, the starting redshift of our study. However, when the overlapped volume is very large (f_hb ≫ 1), our formalism fails to generate a valid correlation function, as its Fourier transform (Appendix B) could fluctuate and yield negative values at large k. To avoid such unphysical results, we have taken a different approach to modeling overlaps in this paper.

We calculate kinetic temperature and Lyα profiles up to very large distances and start shedding outermost shells until f_hb approaches 1. The energy and Lyα photons from excess shells are uniformly distributed in the neutral universe (background as well as remaining shells). The bubbles are still likely to overlap because of randomness in their positions. To account for that, we use (Appendix C)

$$\begin{eqnarray}&&{f}_{h}\simeq \left(1+\displaystyle \frac{{f}_{i}}{{f}_{{hb}}}\right)(1-{e}^{-{f}_{{hb}}})-{f}_{i}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{f}_{b}=1-{f}_{i}-{f}_{h},\end{eqnarray} \tag{ 23 }$$

where f_h corresponds to the actual volume fraction occupied by heating bubbles (single counting of overlapped region), and f_b is the background fraction. Note that f_h approaches f_hb when the T_S profile volume fraction and ionization fraction are small. However, f_h remains less than unity even if the value of f_hb becomes much larger than unity.

In our analysis, we explore two parameters to model X-ray heating and Lyα coupling:

• 1.  
  N_heat: Number of X-ray photons emitted per stellar baryon. For our study, we assume N_heat in the range 0.1–10.0. For larger values of N_heat, the heating is stronger, with higher values of T_K.
• 2.  
  f_L: Ratio of source luminosity of photons between Lyα and the Lyman limit to the luminosity of UV photons. We take 0.1 < f_L < 1000 in this paper. As the value of f_L increases, the coupling between T_S and T_K is stronger. This increases the absolute value of ΔT_B in both emission and absorption.

In this paper, we have not considered scenarios in which the modeling parameters evolve with time.

We can study a simple model with a single bubble size and a shell with uniform spin temperature around it (flat profile). There are small ionization bubbles embedded in larger heated bubbles (Figure 2). Ignoring density perturbations, there are only three values of $\phi =n(1-{T}_{\mathrm{CMB}}/{T}_{S})$ in this universe: ϕ_i = 0, ${\phi }_{h}=1-{T}_{\mathrm{CMB}}/{T}_{\mathrm{heat}}$, and ${\phi }_{b}=1-{T}_{\mathrm{CMB}}/{T}_{\mathrm{bg}}$. When we include the density perturbations, the correlation function can be written as

$$\begin{eqnarray}&&\begin{array}{l}\mu =(1+\xi )\left({\psi }_{h}^{2}({f}_{h}-{f}_{h}{f}_{b}C(r,{R}_{x},{R}_{h},{R}_{h})\right.\\ \ -\,{f}_{i}{f}_{h}C(r,0,{R}_{x},{R}_{h})-\,{f}_{i}(1-{f}_{i})\\ \ \times \,(C(r,0,{R}_{x},{R}_{x})-C(r,0,{R}_{x},{R}_{h})))\\ \ +\,2{\psi }_{h}{\psi }_{b}{f}_{h}{f}_{b}C(r,{R}_{x},{R}_{h},{R}_{h})\\ \ +\,{\psi }_{b}^{2}{f}_{b}(1-{f}_{h}C(r,{R}_{x},{R}_{h},{R}_{h})\\ \ -\,\left.{f}_{i}C(r,0,{R}_{x},{R}_{h}))\right)-{\left({\psi }_{b}{f}_{b}+{\psi }_{h}{f}_{h}\right)}^{2},\end{array}\end{eqnarray} \tag{ 24 }$$

where R_x is the ionization bubble radius and R_h is the heating bubble radius. Here, $C(x,P,Q,R)$ is a function with value between 0 and 1 (Appendix B). This result was derived and analyzed in detail for various limiting cases in RS18. This simple case does not allow for a negative correlation because, within the bubble, the ψ is positively correlated, and it is not correlated outside the bubble.

Zoom In Zoom Out Reset image size

In the case where Lyα coupling is inhomogeneous, we can construct another simple scenario where there is only one ionized bubble size and the profile around the sources has two shells. The first shell is heated and coupled through Lyα radiation. The second shell is nonheated but still coupled. The region outside the second shell (background region) is neither heated nor coupled (Figure 3). This gives three values of $\phi =n(1-{T}_{\mathrm{CMB}}/{T}_{S})$ when ignoring density fluctuations: ϕ_i = 0 in the ionized regions, ϕ₁ = ϕ > 0 in the first shell, ϕ₂ = −ϕ < 0 in the second shell, and ϕ_b = 0 in the background region. For simplicity, we have assumed that the volume occupied by the two shells is the same, ${f}_{1}={f}_{2}={f}_{h}/2$ and ${\phi }_{1}=-{\phi }_{2}=\phi$. Therefore, $\langle \phi \rangle ={f}_{1}{\phi }_{1}+{f}_{2}{\phi }_{2}=0$.

Zoom In Zoom Out Reset image size

From Figure 3, we can see that the correlation will be positive if both points are either in shell 1 or both in shell 2. If one point is in shell 1 and another in shell 2, then the correlation is negative. The two-point correlation function without density perturbation is

$$\begin{eqnarray*}\begin{array}{rcl}\mu & = & {\phi }_{1}^{2}P({\phi }_{1}\cap {\phi }_{1})+{\phi }_{2}^{2}P({\phi }_{2}\cap {\phi }_{2})+2{\phi }_{1}{\phi }_{2}P({\phi }_{1}\cap {\phi }_{2})\\ & = & {\phi }^{2}\left(P({\phi }_{1})-P({\phi }_{1}\cap {\phi }_{i})-P({\phi }_{1}\cap {\phi }_{b})-P({\phi }_{1}\cap {\phi }_{2})\right.\\ & & +\,P({\phi }_{2})-P({\phi }_{2}\cap {\phi }_{i})-P({\phi }_{2}\cap {\phi }_{b})\\ & & -\,\left.P({\phi }_{2}\cap {\phi }_{1})-2P({\phi }_{1}\cap {\phi }_{2})\right),\end{array}\end{eqnarray*}$$

which can be expanded using the same logic used in RS18 and Appendix D:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\mu }{{\phi }^{2}}={f}_{h}(1-C(r,0,{R}_{x},{R}_{h}))-{f}_{i}(1-{f}_{i})\\ \quad \times \,[C(r,0,{R}_{x},{R}_{x})-C(r,0,{R}_{x},{R}_{h})]-\,(1-{f}_{i})\displaystyle \frac{{f}_{h}}{2}\\ \quad \times \,[C(r,{R}_{s},{R}_{h},{R}_{h})-C(r,{R}_{x},{R}_{s},{R}_{h})]\\ \quad -\,(1-{f}_{i}){f}_{h}[C(r,{R}_{x},{R}_{s},{R}_{s})+C(r,{R}_{s},{R}_{h},{R}_{x})\\ \quad -\,C(r,{R}_{s},{R}_{h},{R}_{s})-C(r,0,{R}_{x},{R}_{h})].\end{array}\end{eqnarray} \tag{ 25 }$$

This result can also be obtained by simplifying the complete model (Equation (41)) for approximations used in this section. At large scales, all of the functions C(., ., ., .) tend to unity. In this case, Equation (25) approaches zero. This is expected because, at large scale ($r\gt 2{R}_{h}$), the two-point correlation should vanish if we ignore density correlation. When r = 0, μ = f_hϕ². This is also as expected because there is no probability of two points being in different shells at r = 0. This expression takes a negative value for a range of values of r, where the two points are more likely to be in different shells than in the same shell.

In Figure 4, we have shown a case for this simple model where the autocorrelation of dimensionless brightness temperature is negative at certain scales. For scales ${R}_{x}+{R}_{s}\lt r\lt {R}_{s}+{R}_{h}$, depending on the values of various radii, two randomly chosen points can have a higher probability of being in different shells than in the same shell, driving the correlation to a negative value. For larger scales, ${R}_{s}+{R}_{h}\lt r\lt 2{R}_{h}$, both points have a finite probability of being in the outer shell and a zero probability of being in different shells of the same bubble, which leads the overall correlation at these scales to be positive. On scales $r\gt 2{R}_{h}$, the correlation function is zero. We note that the possibility of negative correlation at any scale depends on values of R_x, R_s, and R_h and their differences. It is entirely possible to have models where the correlation function remains positive at all scales of interest.

Zoom In Zoom Out Reset image size

For the complete model (discussed in detail in the next section and Appendix D), identifying the scales is more difficult because the ionization bubbles have a range of sizes, and the spin temperature profiles have a number of shells. However, if we ignore ξ, which is positively correlated at all scales of interest, we can still identify cases where the autocorrelation of ϕ goes negative. We show one such case in Figure 5. This constitutes the first important result of our analysis, which underlines the importance of correlation analysis in real space to extract the relevant physics. If the correlation function is negative at certain scales, it will point to a very specific geometry of the heated and Lyα coupled regions surrounding the self-ionized bubbles.

Zoom In Zoom Out Reset image size

In this paper, we have explored the correlation function and power spectrum evolution due to inhomogeneity of spin temperature using two modeling parameters: number of X-ray photons per stellar baryon N_heat and the ratio of luminosity of Lyα photons to ionizing photons of the sources f_L. We have taken ζ ∼ 7.5 (Equation (5)), which is in agreement with the reionization optical depth (τ_reion ≃ 0.055) given by Planck (Planck Collaboration et al. 2018). We first present results for the ΛCDM model in the redshift range 10–30, without any additional cooling mechanism needed to explain the recent EDGES results (Bowman et al. 2018). The results below redshift z ≃ 12 are not entirely reliable because several approximations used in our formalism (including excursion set formalism) become less valid and eventually break down when the ionization volume fraction is large (f_i > 0.1; Furlanetto & Oh 2016; Giri et al. 2018). As noted above, to compute the brightness temperature mean and fluctuations, we have taken a size distribution of ionization bubbles, given by excursion set formalism and calculated spin temperature profiles surrounding these bubbles with a number of shells.

In Figure 6, we show the evolution of global H i brightness temperature as a function of redshift for various combinations of modeling parameters. At z ≃ 30, the global signal starts slightly negative as weak collisional coupling drives the spin temperature toward the matter temperature, which is below the CMB temperature (Equations (1) and (2)). At smaller values of z, the behavior is completely determined by the modeling parameters. For higher values of N_heat, the heating starts earlier and the absorption troughs are shallower. For higher values of f_L, the coupling starts earlier as well as stronger, and the overall strength of the signal is larger at all redshifts before complete coupling is achieved. The redshift of the heating transition (when the average gas temperature T_K = T_CMB and the signal goes from absorption to emission) is only dependent on the heating (N_heat) and is independent of Lyα coupling (f_L). It happens sooner for higher values of N_heat.

Zoom In Zoom Out Reset image size

In Figure 7, we show the evolution of the correlation function at scales r = 0.5–8 Mpc for different values of N_heat and f_L, using Equations (2), (6), (9), and (41). The correlation functions are large at small scales and decrease as the distance between two points increases. At very large scale, they approach Equation (16), where the density inhomogeneity is enhanced by the average of (1 − T_CMB/T_S). On intermediate scales, the structure of the correlation function is determined by the size distribution of bubbles and their surrounding T_S profiles.

Zoom In Zoom Out Reset image size

As the main aim of this paper is to model brightness temperature fluctuations owing to incomplete Lyα coupling, we first discuss it qualitatively. As in the case of incomplete heating, the Lyα coupling can also be separated into near and far zones. In the near zone, the emission from a nearby source dominates the coupling. The strength of this coupling is determined by Equation (9), which shows that there always is a distance from the source at which complete Lyα coupling (${y}_{\alpha ,\star }{T}_{k}\gt {T}_{\mathrm{CMB}}$) can be established. Two conditions have to be met to ensure this creates a new length scale in the correlation function. First, this distance must exceed the size of the ionization region. Because the coupling strength scales as R_x³, larger ionization bubbles meet this requirement more easily. Second, the Lyα flux at this point from all of the background sources must be smaller than the flux from the nearby source, or we are in the limit of homogeneous coupling. The background intensity at any point is proportional to R_max(n) (Equation (7)) multiplied by the number density of sources at any redshift. In the initial phase, the background intensity is not high enough to cause complete Lyα coupling. During this phase, there are regions around individual sources that attain complete coupling, and these regions are surrounded by a background in which only partial coupling has been attained. This creates inhomogeneities in Lyα coupling on the scales of these regions. As the intensity builds in the background owing to the birth of new self-ionized regions and it reaches levels sufficient to cause complete coupling, these inhomogeneities disappear. The nature of these inhomogeneities is determined by f_L, the atomic structure of neutral hydrogen, and the excursion set formalism.

At z ≃ 30, only the collisional coupling is effective. Because it is weak and the Lyα intensity is small, all of the curves shown in the figure start with small values of the correlation function with similar strength. The fluctuations are dominated by density perturbations because the number density of ionizing sources is small and the ionization, heating, and Lyα coupled fractions are tiny. For small f_L, the correlation declines with time as the collisional coupling weakens, owing to the fall of the number density of particles with the expansion of the universe. This situation is only reversed when Lyα coupling becomes efficient. For higher values of f_L, this coupling occurs sooner, leading to an increase in fluctuations with time until complete coupling is achieved. After this period, the fluctuations are determined by heating and ionization inhomogeneities. The position of the first peak in Figure 7 depends strongly on the heating parameter as the peak is determined by the thermal evolution of the background; the correlation only starts decreasing when the background is sufficiently heated. This phase is described in detail in RS18.

The correlation function approaches zero near the heating transition. At large scales, the signal vanishes completely when ${T}_{\mathrm{CMB}}/{T}_{K}$ approaches f_n (Equation (16)). We again see the effect of N_heat on the redshift of the heating transition. Inhomogeneous collisional coupling and the shape of the temperature profiles, which are determined by the spectrum of X-ray photons (α and ν_min, see RS18), can potentially change the redshift and depth of the heating transition by a small amount. We do not study this effect in the current paper. After the heating transition, the effect of inhomogeneous T_S decreases, and the main source of fluctuations is ionization inhomogeneities. Their effect is suppressed if the heating or coupling is not saturated (very small values of N_heat and f_L). In general, larger N_heat creates larger heating profiles and causes correlation at larger scales. However, these profiles also merge sooner and wipe out heating fluctuations at those scales.

The H i signal from the epoch of cosmic dawn and reionization has been extensively studied in the literature using semianalytic methods and large-scale simulations (e.g., Pritchard & Furlanetto 2007; Santos et al. 2008, 2010; Baek et al. 2010; Visbal et al. 2012; Mesinger et al. 2013, 2016; Tashiro & Sugiyama 2013; Pacucci et al. 2014; Fialkov et al. 2015, 2017; Ghara et al. 2015; Ross et al. 2017). Our formalism has been built on geometric arguments, which are intuitively easier to visualize in real space. The entire correlation structure can be written in terms of a single polynomial function: C(., ., ., .) (Appendix B). Taking a Fourier transform with respect to the first argument r of this function yields the power spectrum. In Figure 8, we show the evolution of the power spectrum ${{\rm{\Delta }}}^{2}={k}^{3}P(k)/2{\pi }^{2}$ for a range of k. These figures show an evolutionary trend similar to the correlation functions (Figure 7).

Zoom In Zoom Out Reset image size

While comparing our results with simulations, we have focused on three features: (1) the number of peaks in the power spectrum, (2) the amplitude of Δ²(k) for a range of scales k ∼ 0.1–0.5 Mpc⁻¹, and (3) the difference between the redshift of the heating transition and the redshift of the power spectrum minimum.

Existing results in the literature show that, for k ≃ 0.1–0.5 Mpc⁻¹, there are generally two or three peaks of a power spectrum as a function of redshift (Santos et al. 2008; Baek et al. 2010; Mesinger et al. 2013, 2016; Ghara et al. 2015; Fialkov et al. 2017). At high redshifts, when the Lyα coupling and X-ray heating commence, they create fluctuations of T_S in the medium. If fluctuations in these two fields dominate at widely different times, there will be two distinct peaks at high redshift: one due to coupling inhomogeneities and the other (generally at lower redshift than the former) due to heating inhomogeneities (Pritchard & Furlanetto 2007; Chen & Miralda-Escudé 2008; Ahn et al. 2015). After the heating transition, there is a third, smaller peak at low redshifts, when the power spectrum is dominated by ionization inhomogeneities (e.g., Pritchard & Furlanetto 2007; Fialkov et al. 2014; Ghara et al. 2015). We studied the two peaks that are due to heating and ionization inhomogeneities in RS18.

The third peak owing to the inhomogeneous Lyα coupling is expected at early times for the following reason: as large self-ionizing regions are born, these fluctuations should initially build and then diminish as the contrast between Lyα coupling in the near and far zone reduces, finally disappearing when complete Lyα coupling is established. We see this additional peak in Figure 9 in which the effect of density perturbations is neglected. However, in the complete model, we find a weak peak owing to this effect only at small scales in the power spectrum (k = 2 Mpc⁻¹ in Figure 8). Generally this possible additional peak is masked by density perturbations. The strength of this peak can be understood in terms of the evolution of number density of large self-ionized regions at early times and the influence region of Lyα photons. For z < 30, the number density of the self-ionizing bubble builds exponentially in the excursion set formalism. While this creates inhomogeneities owing to the geometry seen in Figure 3, it also causes a rapid build-up of the background Lyα photons, rapidly destroying the contrast between the near and far zones. At any point, the background flux gets a nearly equal contribution from sources within the (comoving) radius R_max(n) (Equation (7)). This radius is close to 600 comoving Mpc for n = 2 at z ≃ 25. This large influence region contributes to wiping out the contrast in a short span. We note that R_max(n) (for small principal quantum numbers n > 2) determines the length scale at which the physics needs to be captured to study the Lyα generated inhomogeneities. It might be difficult to achieve it using an N-body simulation as the box size is generally smaller than this length scale. In our results for large N_heat and small f_L, we get three peaks at large k (top left panel of Figure 8) because the heating inhomogeneities, and by extension collisional coupling inhomogeneities, dominate before the Lyα inhomogeneities commence.

Zoom In Zoom Out Reset image size

While Δ²(k) agrees with the results of simulations for smaller scales, for $k\simeq 0.1\,{\mathrm{Mpc}}^{-1}$, our results give less power as compared to simulations (bottom right panel of Figure 8). Given the small sizes of ionization bubbles at high redshifts, the only contribution to fluctuations at k ∼ 0.1 Mpc⁻¹ is due to density and spin temperature fluctuations. For the models of heating and Lyα coupling used in this paper, the contribution due to heating and coupling is dominated by faraway sources, which diminishes fluctuations at large scales. Therefore, at these scales, the spin temperature inhomogeneities are negligible, and only the density fluctuations are enhanced by the average contrast of H i spin temperature with CMB temperature. While it is conceivable that the higher power at large scales in the simulation is due to the finite box size (for a discussion, see Zahn et al. 2011; Ghara et al. 2015), which might not allow one to take into account the contribution of faraway sources whose effect tends to homogenize the fluctuations of T_s at large scales, a more detailed comparison with simulations is hard as the parameter range used is generally not the same.

When the average spin temperature equals the CMB temperature during the heating transition, the power spectrum at small k reaches a minimum value (Figure 8). However, for larger k (small r), even during the heating transition there are significant fluctuations due to inhomogeneities of spin temperature and ionization, which delay the minima of the power spectrum. Figure 8 shows that, for k = 2 Mpc⁻¹, the minimum of the power spectrum depends on the value of f_L even though the heating transition is independent of it. In general, the minimum of the power spectrum occurs during or after the global heating transition, depending on the scales and modeling parameters. This is in agreement with simulations that explore the dependence of these inhomogeneities on modeling parameters (Mesinger et al. 2013, 2016; Ghara et al. 2015; Fialkov et al. 2017).

Recent EDGES observations (Bowman et al. 2018) reported a sky-averaged absorption feature of strength ΔT ≃ −500 mK in the frequency range 70–90 MHz, corresponding to a redshift range 15–19 for the redshifted H i line. It can be shown that for a standard recombination and thermal history, the minimum temperature of the gas at z ≃ 17 is T_K ≃ 6 K. It follows from Equations (1) and (2) that the absorption trough should not have been deeper than −180 mK.

One possible explanation of the EDGES result is that there is additional radio background in the redshift range 15 < z < 19; in this case, we can replace T_CMB with the T_CMB + T_radio in Equation (1) in the relevant redshift range (Ewall-Wice et al. 2018; Feng & Holder 2018; Sharma 2018). With this replacement and a suitable choice of T_radio, we can rederive all our results in this paper for compatibility with the EDGES result.

Another plausible explanation invokes the additional cooling of baryons that is due to interaction between dark matter and baryons.⁷ In this case, we can explain the EDGES detection using Equations (1) and (2) if (1) Lyα photons globally couple the spin temperature to matter temperature, that is, T_K y_α ≫ T_CMB, such that T_S = T_K at z ≃ 19; and (2) T_K ≃ 2.5 K.⁸ We consider this theoretical extension to explain the EDGES result in this paper.⁹

We consider a fiducial model that roughly fits the EDGES results. We assume dark matter–baryon interaction of the form described in Barkana (2018), with a cross section ${\sigma }_{1}=5\times {10}^{-24}\,{\mathrm{cm}}^{2}$ and the ratio of dark matter to proton mass ${m}_{{dm}}/{m}_{p}=0.001$. Such an interaction helps cool the baryon gas temperature sufficiently to explain the EDGES result. For f_L=2 and N_heat = 0.08, we get an absorption trough in the global signal similar to the one observed by EDGES data (Figure 6). In Figures 7 and 8, we show the correlation functions and power spectrum for a range of scales, taking into account our fiducial model to replicate the EDGES results.

The main effect of EDGES observations on the expected correlations is to boost the signal by nearly an order of magnitude in the redshift range 15 < z < 19, even as compared to the most optimistic models¹⁰ (low N_heat and high f_L) in the usual case.¹¹ This is entirely owing to a decrease in T_S that boosts s in Equation (3). The correlation function scales roughly as ${(1-s)}^{2}$. Even though this result has been arrived at within the framework of a model involving baryon cooling, this result holds when excess radio background is responsible for the deep absorption feature because s increases by a similar factor in that case too.

This result is also robust to change in other modeling parameters f_L and N_heat. Therefore, an important prediction of the EDGES detection is that the deeper absorption trough in the global signal gives a corresponding increase in the fluctuating component too. Even when the spin temperature field is very uniform, with Lyα-coupled and unheated gas, the low spin temperature would enhance the underlying density inhomogeneities (Equation (17)). As we have shown in Section 4.2, the spin temperature field can be negatively correlated at some scales. In such scenarios, the fluctuating component of the signal could be less than given by Equation (17). We tried to produce such models for parameters needed to explain the EDGES data and found it very difficult to anticorrelate the spin temperature field. Hence we infer that the minimum value of the correlation function that corresponds to the EDGES signal is given by Equation (17), using the temperature derived from the trough in the global signal. If the global signal has a trough of ∼500 mK at z ≃ 17, then assuming complete coupling and no heating (no fluctuations due to these two fields), the autocorrelation function at r = 2 Mpc and r = 4 Mpc should be ∼5100 (mK)² and ∼2500 (mK)², respectively.

Many currently operational (e.g., LOFAR, MWA, PAPER, GMRT) and upcoming radio interferometers (HERA, SKA) have the capability to detect the fluctuating component of the H i signal in the redshift range 8 < z < 25 (for details, see, e.g., Mesinger et al. 2014; Ahn et al. 2015; Koopmans et al. 2015). We discuss the detectability of this signal especially in light of the recent EDGES results. These radio interferometers directly measure visibilities and their correlations, which can be related to the power spectrum of the H i signal (e.g., Bharadwaj & Sethi 2001; Zaldarriaga et al. 2004). This data analysis can readily be extended to the image plane (which is a byproduct of the analysis pipeline; e.g., Patil et al. 2017 for LOFAR). This means real-space correlation functions can also be used for computation of the signal (e.g., Sethi & Haiman 2008). As an approximate rule, one could use r ≃ π/k to shift from Fourier to real space.

SKA1-LOW is expected to detect the H i signal at z ≃ 16 with a signal-to-noise ratio varying from 100 to 10 for 0.1 < k < 0.6 Mpc⁻¹ for a signal strength Δ²(k) ≃ 10² (mK)² (for details, see, e.g., Koopmans et al. 2015). EDGES results predict a signal strength nearly a factor of five larger, which means the signal-to-noise ratio would be significantly higher (Figure 8). The ongoing experiment LOFAR's best upper limit corresponds to (80 mK)² (k ≃ 0.05 Mpc⁻¹) in the redshift range 9.6 < z < 10.6 in 10 hr of integration.¹² LOFAR has the frequency range to probe the redshift range of EDGES detection. If the noise properties at smaller frequencies (≃80 MHz) behave roughly as the one observed at higher frequencies (≃110 MHz), LOFAR might be able to detect the signal in a few hundred hours of integration.

SKA1-LOW (and SKA2) have the capability to detect the EoR signal in the redshift range 20 < z < 25. A signal of 100 (mK)² at z ≃ 25 is detectable by a deep SKA1-LOW survey with signal-to-noise ratio of five for k < 0.1 Mpc⁻¹ (Koopmans et al. 2015). However, we notice that the EDGES result places a significant constraint on the signal in this era. This, as noted above, is because the EDGES results imply a smaller value of f_L, which results in a smaller signal for z > 19 as compared to models with a larger f_L, which give a significantly higher signal (Figure 8).