WOUTHUYSEN–FIELD COUPLING IN THE 21 cm REGION AROUND HIGH-REDSHIFT SOURCES

The property of the ionized and heated regions around an individual luminous object is dependent on the luminosity (or mass), the spectrum of UV photon emission, and the time evolution of the center object (Cen 2006; Liu et al. 2007). We will not work on a specific model, but consider only the common features. These halos generally consist of three spheres. The most inner region of the halo is the highly ionized Strömgren sphere, or the H ii region, in which the fraction of neutral hydrogen f_H i = n_H i/n_H is no more than 10⁻⁵, where n_H i and n_H are, respectively, the number densities of neutral hydrogen H i and total hydrogen. The temperature of the H ii region is about 10⁴ K. The physical radius of this sphere is in the range of a few to a few tens of kpc. The second region is the 21 cm emission shell, which is just outside the H ii region. The physical size of this shell is similar to that of the H ii sphere. The temperature of hydrogen gas in the emission shell is in the range 10² < T < 10⁴ K due to the heating of UV photons. The fraction f_H i is in the range of 0.1–1. The third region is the 21 cm absorption shell, which is outside the 21 cm emission region. The physical size is on the order of 100 kpc. The temperature of this region is lower than T_cmb = T_CMB(1 + z), where T_CMB is the temperature of CMB today. The timescale of the formation of the halos is about 10⁶ years. The lifetime of the halos is about the same as the lifetime of the first stars.

At the epoch of redshift z ≃ 20, the reionization region consists of isolated patches around first sources. Most Lyα photons in the 21 cm regions should come from the central source and the subsequent re-emission processes. If one can estimate the center object as a Lyα emitter, the emission of Lyα photons in number per unit time would be about dN_Lyα/dt = 10⁵³ s⁻¹. The recombination of H ii and electron in the Strömgren sphere is also a source of Lyα photons. If the physical radius of the Strömgren sphere is ∼10 kpc, the emission intensity is about dN_Lyα/dt = 10⁴⁹ s⁻¹.

Using the parameters of the concordance ΛCDM model, the optical depths of photons with Lyα resonant frequency ν₀ in the 21 cm region are

$$\begin{eqnarray} \tau (\nu _0)&=&n_{\hbox{\scriptsize {H\,{\sc i}}}}\sigma _0R=4.9 \times 10^6 f_{\hbox{\scriptsize {H\,{\sc i}}}}\left(\frac{T}{10^4}\right)^{-1/2}\nonumber\\ &&\times\left(\frac{1+z}{20}\right)^3\left(\frac{\Omega _bh^2}{0.022}\right) \left(\frac{R}{{\rm 10 \;kpc}}\right), \end{eqnarray} \tag{ 1 }$$

where σ₀ is the cross section of the resonant scattering at the frequency ν₀, R is the physical size of the considered sphere, and τ is the distance R in the units of the mean free path of ν₀ photons at redshift z. Equation (1) shows that the optical depth of the 21 cm signal regions with f_H i ⩾ 0.1 is τ(ν₀) ⩾ 10⁶.

It is also useful to define a dimensionless time η = cn_H iσ₀t as

$$\begin{eqnarray} t&=&\eta /cn_{\hbox{\scriptsize {H\,{\sc i}}}}\sigma _0=6.7\times 10^{-3}f^{-1}_{\hbox{\scriptsize {H\,{\sc i}}}}\left(\frac{T}{10^4}\right)^{1/2}\nonumber\\ &&\times\left(\frac{20}{1+z}\right)^3\left(\frac{0.022}{\Omega _b h^2}\right)\eta \;\, \rm yr. \end{eqnarray} \tag{ 2 }$$

This is the time in the units of mean free flight time. For a timescale t ≃ 10⁵ yr, the scale of η is ≃10⁷.

There are, at least, two RF problems with the WF coupling in the 21 cm regions: (1) How to provide enough Lyα photons in such opaque medium? (2) Can the RF in the 21 cm region keep the WF coupling to work well?

If photons always keep the frequency ν₀, their spatial diffusion can be described as a random walk process (Chandrasekhar 1943). The mean number of scattering required for the diffusion over a range R is on the order of τ². Thus, the time for the diffusion over the 21 cm range is η ≃ τ², corresponding to the time t = τ²/cn_H iσ₀ ⩾ 10¹² years. This diffusion mechanism, obviously, is useless for the 21 cm regions.

The timescale of the diffusion would be substantially reduced if the frequency of the Lyα photons can have a small shift from ν₀ to ν = ν₀ ± Δν, because the optical depth τ at the frequency ν₀ ± Δν is significantly less than τ(ν₀). One can estimate the frequency shift Δν by the condition of optical depth τ(ν₀ ± Δν) ≃ 1, which is

$$\begin{equation} \tau (\nu _0)\frac{\phi (x,a)}{\phi (0,a)} \le 1, \end{equation} \tag{ 3 }$$

where ϕ(x, a) is the Voigt function for the resonant line ν₀ profile as (Hummer 1965)

$$\begin{equation} \phi (x,a)=\frac{a}{\pi ^{3/2}}\int ^{\infty }_{-\infty } dy \frac{e^{-y^2}}{(x-y)^2+a^2}, \end{equation} \tag{ 4 }$$

where the dimensionless variable x is defined by x = (ν − ν₀)/Δν_D and Δν_D = ν₀v_T/c is the Doppler broadening of hydrogen gas with thermal velocity $v_T=\sqrt{k_bT/2m_{\rm H}}$. Therefore, x measures the frequency deviation Δν = |ν − ν₀| in the units of the Doppler broadening. The parameter a in Equation (4) is the ratio of the natural to the Doppler broadening. For the Lyα line, a = 2.35 × 10⁻⁴(T/10⁴)^−1/2. In terms of x, the solution of Equation (3) is x ⩾ 3.

An effective mechanism of the frequency shift is given by the diffusion of the photon distribution in the frequency space. Considering this mechanism, the number of scattering required for diffusion over a physical distance R is no longer equal to τ², but roughly on the order of τ (Osterbrock 1962; Harrington 1973; Avery & House 1968; Adams 1972). Since τ ∝ R, the timescale of the spatial diffusion over size R is comparable to R/c. Problem (1) would then be solved.

Problem (2) still remains. If photons with the frequency ν < ν₀ − Δν or ν>ν₀ + Δν take a faster spatial diffusion than the ν₀ photons, how can we keep the WF coupling to work? Without ν₀ photons, one cannot have a local Boltzmann distribution around ν₀. Therefore, the photons with the frequency ν < ν₀ − Δν or ν>ν₀ + Δν should be brought back to the frequency ν₀. We must study whether the frequency space diffusion mechanism can restore ν₀ photons from photons with the frequency ν ≠ ν₀.

Cosmic expansion also leads to a deviation of photon frequency from ν₀ to a lower one ν₀ − Δν, which speeds up the spatial transfer. However, we also need a mechanism to restore Lyα photons from Hubble-redshifted photons. Therefore, in terms of the WF coupling in the 21 cm region, we must study the kinetics of Lyα photons in both physical and frequency spaces.

We first estimate the timescale needed for the frequency shift. We can use the RF equation of a homogeneous and isotropically expanding infinite medium consisting of neutral hydrogen. The equation of the mean intensity J in terms of the photon number is (Hummer & Rybicki 1992; Rybicki & Dell'antonio 1994)

$$\begin{eqnarray} \frac{\partial J(x, \eta)}{\partial \eta }& =& - \! \phi (x,a)J(x, \eta) \nonumber \\ & & +\int \mathcal {R}(x,x^{\prime };a)J(x^{\prime }, \eta)dx^{\prime } + \gamma \frac{\partial J}{\partial x} +S(x, \eta).\nonumber\\ \end{eqnarray} \tag{ 5 }$$

The parameter γ⁻¹ measures the number of scattering during a Hubble time, given by

$$\begin{equation} \gamma ^{-1}= 1.4 \times 10^6 h^{-1}f_{\hbox{\scriptsize {H\,{\sc i}}}}\left(\frac{0.25}{\Omega _M}\right)^{1/2}\left(\frac{\Omega _bh^{2}}{0.022}\right) \left(\frac{1+z}{20}\right)^{3/2}. \end{equation} \tag{ 6 }$$

The re-distribution function $\mathcal {R}(x,x^{\prime })$ gives the probability of a photon absorbed at the frequency x', and re-emitted at the frequency x. It depends on the details of the scattering (Henyey 1941; Hummer 1962, 1969). If we consider coherent scattering without recoil, the re-distribution function with the Voigt profile Equation (4) is

$$\begin{eqnarray} \qquad \quad \mathcal {R}(x,x^{\prime };a)&=& \frac{1}{\pi ^{3/2}}\int ^{\infty }_{|x-x^{\prime }|/2}e^{-u^2} \bigg[ \tan ^{-1}\left(\frac{x_{\min }+u}{a}\right)\nonumber\\ && -\tan ^{-1}\left(\frac{x_{\max }-u}{a}\right)\bigg]du, \end{eqnarray} \tag{ 7 }$$

where x_min = min(x, x') and x_max = max(x, x'). In the case of a = 0, i.e., considering only the Doppler broadening, the re-distribution function is

$$\begin{equation} \mathcal {R}(x,x^{\prime })=\frac{1}{2}e^{2bx^{\prime }+b^2}{\rm erfc}[{\rm max}(|x+b|,|x^{\prime }+b|)], \end{equation} \tag{ 8 }$$

where the parameter b = hν₀/mv_Tc = 2.5 × 10⁻⁴(10⁴/T)^1/2 is due to the recoil of atoms. The re-distribution function of Equation (8) is normalized as $\int _{-\infty }^{\infty } \mathcal {R}(x,x^{\prime })dx^{\prime }=\phi (x,0)\equiv \phi _g(x) =\frac{1}{\sqrt{\pi }}e^{-x^2}$. The numerical algorithm to solve Equation (5) has been given in Roy et al. (2009a, 2009b).

On the right-hand side of Equation (5), the first term is the absorption at the resonant frequency x, the second term is the re-emission of photons with frequency x by scattering, and the third term describes the Hubble redshift of photons. We first solve Equation (5) by dropping the terms of absorption and re-emission. The equation is

$$\begin{equation} \frac{\partial J}{\partial \eta }= \gamma \frac{\partial J}{\partial x}+S(x,\eta). \end{equation} \tag{ 9 }$$

Assuming the source is S = Cϕ_s(x), where ϕ_s(x) is the normalized frequency profile of the source photons, the analytic solution of Equation (9) is (Rybicki & Dell'antonio 1994)

$$\begin{equation} J(x,\eta)=J(x+\gamma \eta,0) + C\gamma ^{-1}\int ^{x+\gamma \eta }_{x}\phi _s(x^{\prime }) dx^{\prime }, \end{equation} \tag{ 10 }$$

where J(x, 0) is the initial flux. Define the mean frequency by

$$\begin{equation} \bar{x}(\eta)\equiv \frac{\int x J(\eta, x)dx}{\int J(\eta, x)dx}. \end{equation} \tag{ 11 }$$

If the initial flux is J(x, 0) = 0, one can show that for any profile ϕ_s(x) we have

$$\begin{equation} \bar{x}(\eta)=- \gamma \eta /2. \end{equation} \tag{ 12 }$$

As expected, the speed of redshift $d\bar{x}/d\eta$ is a constant. For the Hubble expansion, a frequency shift x needs a time η ≃ γ⁻¹x. Thus, in order to have the frequency shift x ⩾ 3, the timescale is η ⩾ 6 × 10⁶, corresponding to t ⩾ 4 × 10⁴ years. This scale seems to be short enough compared with the lifetime of first stars. However, Hubble redshift is less effective comparing with the resonant scattering. This point can be seen in Figure 1, which gives the solution of Equation (5) with the re-distribution function Equation (8). The source is taken to be $S(x)= \phi _g(x) =(1/\sqrt{\pi })e^{-x^2}$. We also take the parameter γ to be 10⁻⁶ [Equation (6)], and ignore the recoil, i.e., b = 0.

Zoom In Zoom Out Reset image size

Figure 1 shows that the diffusion in the frequency space leads to a flat plateau with width |x| ⩽ 3 when η is as small as η ≃ 10⁴, or t ≃ 10² years. This timescale is much less than that of the Hubble redshift. Therefore, the major mechanism to produce photons with frequency |x| ⩾ 3 is given by the resonant scattering.

From Figure 1, we can see that the profile of photons is almost symmetric with respect to x = 0 (or to ν = ν₀) until about the time η = 10⁶. The redshift effect can only be seen from the curve of η ⩾ 10⁶. That is, the resonant scattering impedes cosmic redshift. The impediment is due to the "bounce back" of resonant scattering. Regardless of whether the frequency of the absorbed photons is larger or smaller than ν₀, the mean frequency of the re-emitted photons is always ν₀. Therefore, the resonant scattering will bring back some redshifted photons to the frequency ν₀. Thus the net effect of resonant scattering (absorption and re-emission) is, on average, to bounce redshifted photons back to the resonant frequency ν₀, and to restore the symmetry with respect to x = 0. This bounce back mechanism is the key of restoring ν₀ photons and the WF coupling (see Section 4).

To illuminate the bounce back effect, we calculate the mean frequency $\bar{x}(\eta)$ for the solutions shown in Figure 1. The result is plotted in Figure 2. The solution of pure cosmic redshift Equation (12) is also shown in Figure 2 as the dashed curve. We can see from Figure 2 that the speed of photon frequency shift $d\bar{x}/d\eta$ when considering resonant scattering is in general less than that in the case of pure cosmic redshift. When the time η is large, many photons have been redshifted out of the frequency range of ϕ_g(x). In this case, the "bounce back" effect ceases, and the redshift speed is recovered to $d\bar{x}/d\eta =-\gamma /2$. From Figure 2 one can see once again that the Hubble redshift can produce frequency shift from x = 0 to |x| ≃ 3 only when η ⩾ 10⁷, or t ≃ 10⁵ years. This timescale may not be short enough to match the 21 cm region with a short lifetime. When the timescale of resonant scattering is less than the timescale of Hubble expansion, the frequency shift of the Hubble expansion slows down significantly by the resonant scattering.

Zoom In Zoom Out Reset image size

Considering a photon source located at the central region of a uniformly distributed expanding medium, we can use the RT equation of the specific intensity I(η, r, x, μ) as follows

$$\begin{eqnarray} &&\hspace{-1.8pc}{\partial I\over \partial \eta } + \mu \frac{\partial I}{\partial r}+\frac{(1-\mu ^2)}{r}\frac{\partial I}{\partial \mu } - \gamma \frac{\partial I}{\partial x} = - \phi (x)I \nonumber \\[-4pt] & & + \int \mathcal {R}(x,x^{\prime })I(\eta, r, x^{\prime },\mu ^{\prime })dx^{\prime }d\mu ^{\prime } + S, \end{eqnarray} \tag{ 13 }$$

where μ = cos θ is the direction relative to the radius vector r. The dimensionless coordinate r is rescaled from the physical coordinate r_p as

$$\begin{equation} r_p=2.1\times 10^{-3} f^{-1}_{\hbox{\scriptsize {H\,{\sc i}}}} \left(\frac{T}{10^4}\right)^{1/2} \left(\frac{20}{1+z}\right)^3 \left(\frac{0.022}{\Omega _b h^2}\right) r \;{\rm pc}. \end{equation} \tag{ 14 }$$

With these variables, the propagation of a signal with the speed of light will be described by the equation r = η + const. It would still be reasonable to use the isotropic approximation of the re-distribution (Mihalas et al. 1976).

When the optical depth is large, the Eddington approximation would be proper. It is

$$\begin{equation} I(\eta,r,x,\mu)\simeq J(\eta, r, x) + 3\mu F(\eta,r,x), \end{equation} \tag{ 15 }$$

where $J(\eta,r,x)=\frac{1}{2}\int _{-1}^{+1}I(\eta,r,x,\mu)d\mu$ is the angularly averaged specific intensity and $F(\eta,r,x)=\frac{1}{2}\int _{-1}^{+1}\mu I(\eta,r,x,\mu)d\mu$ is the flux. Defining j = r²J and f = r²F, Equation (13) yields the equations of j and f as

$$\begin{eqnarray} {\partial j\over \partial \eta } + \frac{\partial f}{\partial r} & = & -\! \phi (x)j + \int \mathcal {R}(x,x^{\prime })j dx^{\prime }+ \gamma \frac{\partial j}{\partial x}+ r^2S,\nonumber\\ \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \frac{\partial f}{\partial \eta } + \frac{1}{3} \frac{\partial j}{\partial r} - \frac{2}{3}\frac{j}{r} & = & -\! \phi (x)f. \end{eqnarray} \tag{ 17 }$$

The numerical algorithm for Equations (16) and (17) is given in the Appendix.

Since photons with frequency shifted away from ν₀ would be optically thin, Equations (16) and (17) will no longer be a good approximation when the frequency of photons is shifted to the optical thin case. From Figures 1 and 2, one can see that within η ⩽ 10⁶, most photons are still trapped in the frequency |x| ⩽ 3, for which the optical depth is larger than 1. Therefore, the Eddington approximation would be proper, at least, until η is as large as about 10⁶.

We first solve the Equations (16) and (17) by dropping all terms on the transfer in frequency. It has been shown that the source term can be replaced by a boundary condition of f and j at r = r₀ (Qiu et al. 2006). The equations then become

$$\begin{eqnarray} {\partial j\over \partial \eta } + \frac{\partial f}{\partial r} & = & 0, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} {\partial f\over \partial \eta } + \frac{1}{3} \frac{\partial j}{\partial r} - \frac{2}{3}\frac{j}{r}& = & -\! f. \end{eqnarray} \tag{ 19 }$$

We take the boundary condition at r₀ = 1 to be

$$\begin{equation} j(\eta, 1)=3S_0, \qquad f(\eta, 1)=S_0=10^6. \end{equation} \tag{ 20 }$$

Since Equations (18) and (19) are linear, the parameter S₀ is not important if we are only interested in the shape of j and f as a function of η and r. The initial condition is taken to be

$$\begin{equation} j(0,r)=f(0,r)=0. \end{equation} \tag{ 21 }$$

The solution of f(η, r) is presented in Figure 3. The dashed line is f = S₀ = 10⁶, which is the time independently exact solution of Equation (18). It is a horizontal straight line because the flux f = r²F is r-independent and satisfies the conservation of the photon number. Figure 3 shows that the spatial size r of f(η, r) roughly satisfies $r\propto \sqrt{\eta }$. Therefore, without resonant scattering, the diffusion basically is a random walk process as discussed in Section 2.2.

Zoom In Zoom Out Reset image size

We now turn to the solution of Equations (16) and (17) with resonant scattering and Hubble redshift. The relevant parameters are given by b = 0 and γ = 10⁻⁵. We use the boundary condition at r = 0

$$\begin{equation} j(\eta, 0, x)=0, \qquad f(\eta, 0, x)=S_0 \phi _s(x), \end{equation} \tag{ 22 }$$

where the frequency profile is taken to be the Gaussian profile, i.e., ϕ_s(x) = ϕ_g(x). The parameter S₀ is still taken to be 10⁶. The initial condition is similar to Equation (21), i.e., j(0, r, x) = f(0, r, x) = 0. The solutions of f(η, r, x) are plotted in Figure 4.

Zoom In Zoom Out Reset image size

All the solutions of Figure 4 show two remarkable peaks at x ≃ ±(2 − 3) for all radius r. The amplitude of f(η, r, x) at the peaks is higher than that at x = 0 by a factor of 10 to 10². It shows that the flux is dominated by photons with frequency x ≃ ±(2 − 3). That is, the spatial transfer is carried out by photons of x ≃ ±(2 − 3). The amplitude of the flux at the saturated peaks is basically r-independent, and the saturated values of ∫f(η, r, x)dx are r-independent. This is consistent with the conservation of the photon number.

More importantly, the amplitude of the peaks approaches its saturation at about η ≃ 200 for r = 50, η ≃ 400 for r = 100, and η ≃ 2000 for r = 500. That is, the time η needed for the spatial transfer over size r is roughly proportional to r. This is very different from the random walk relation η ∝ r², or the results shown in Figure 3.

It should be emphasized that the photons at the two peaks are not only those with x ≃ ±(2 − 3) directly from the source S₀ϕ_s(x), but also include the photons with frequency shift from x ≃ 0 to x ≃ ±(2 − 3). This point can be shown by a source of the profile with smaller width, say, $\phi _s(x)=(1/\sqrt{\pi /2})e^{-2x^2}$. In Figure 5, we plot the results. For the source of $\phi _s(x)=(1/\sqrt{\pi /2})e^{-2x^2}$, the number of photons with x ≃ ±(2 − 3) is much less than that of the source $(1/\sqrt{\pi })e^{-x^2}$. We can see that all the features in Figure 5 are the same as those in Figure 4. The two peaks are still located at x ≃ ±(2 − 3). The timescale for approaching saturation in Figure 5 is also the same as that in Figure 4.

Zoom In Zoom Out Reset image size

Therefore, photons at the peaks of x ≃ ±(2 − 3) should come from the frequency shift from x = 0 to x ≃ ±(2 − 3) by resonant scattering. That is, resonant scattering provides a shortcut of the spatial transfer: first to shift x = 0 photons to x ≃ ±(2 − 3) photons, and then to speed up spatial transfer. This mechanism is the same as that for the escape of resonant photons from opaque clouds (Osterbrock 1962; Harrington 1973; Avery & House 1968; Adams 1972). It allows the Lyα photons to be able to transfer over to the 21 cm regions in time η proportional to its size r.

We have mentioned in Section 4.2 that the flux f(η, r, x) lacks ν₀ (or x = 0) photons even when f approaches its saturation. This is, obviously, not good for the WF coupling. However, we found that the situation may not be so if we consider the solution of the mean specific intensity j. Figure 6 presents the solution j(η, r, x) of Equations (16) and (17) at r = 50, 100, and 500. The parameters b = 0 and γ = 10⁻⁵ are the same as the solutions in Figure 4. The results are given in Figure 6.

Zoom In Zoom Out Reset image size

When the time η is small, j(η, r, x) is similar to f(η, r, x), having a valley around x = 0 and two peaks at x ≃ ±(2 − 3). However, different from the flux f(η, r, x), the amplitude of j(η, r, x) around x = 0 quickly increases. In the saturated state it is about the same as the peaks. That is, although the flux always has a valley around x = 0, the ν₀ photons are quickly restored in the mean intensity. The timescale of approaching its saturated state is also proportional to r, not r². Therefore, the restoration of resonant photons is due to the resonant scattering "bounce back" (Section 3.3), which pushes photons with frequency x ≃ ±(2 − 3) back to x ≃ 0.

The shape of j(η, r, x) around x = 0 is a flat plateau, which is similar to that in Figure 1. As have been shown in Paper I, the flat plateau of j(η, r, x) at b = 0 will become the local Boltzmann distribution if the recoil is considered. We can expect that the flat plateau of Figure 6 will also show a local Boltzmann distribution if b ≠ 0. We calculate the solutions j(η, r, x) and f(η, r, x) of Equations (16) and (17) at r = 500 with the same parameters as those in Figure 6, but with b = 0.3. We use a large b, because it is easier to see the slope of the local Boltzmann distribution. The results are given in Figure 7.

Zoom In Zoom Out Reset image size

Figure 7 clearly shows a local Boltzmann distribution within the range |x| ⩽ 2 as

$$\begin{equation} j(\eta,r,x)\simeq j(\eta,r,0)e^{-2bx}= j(\eta,r,0)e^{-h(\nu -\nu _0)/k_BT}. \end{equation} \tag{ 23 }$$

We find the slope to be ln j(η, 0) − ln j(η, 1) = 0.59, which is well consistent with 2b = 0.6. Figure 7 shows that at r = 500, the onset of the WF coupling can occur as early as η = 900, but the amplitude of the local Boltzmann distribution at that time is much lower than its saturated value by a factor of 10². The amplitude of the local Boltzmann distribution substantially increases with time. Therefore, the "bounce back" mechanism keeps the WF coupling to work with a timescale η larger than a few hundreds, i.e., a few hundred collisions. This result is the same as that in Paper I.

The solution of the flux f(η, r, x) in Figure 7 is not very different from that in Figure 4. The only difference between the two figures is that the former is asymmetric with respect to x = 0, i.e., the peak at x < 0 is stronger than that of x > 0, while the latter is symmetric. This is simply due to the recoil of b ≠ 0 leading more photons to move to x < 0. Neither flat plateau nor local Boltzmann distribution is shown in the flux f(η, r, x). There is always a valley around ν₀ even when f(η, r, x) is in its saturated state. It once again indicates that the flux is dominated by photons of x ≃ ±(2 − 3).

The two components f and j of the Eddington approximation are effective in revealing the functions of the resonant scattering. The flux f(η, r, x) describes the photons in transit. It shows that the spatial transfer within opaque media is mainly via photons with frequency shifted to x ≃ ±(2 − 3), which are easy for escaping. The mean intensity j(η, r, x) describes the restoration of x = 0 photons and the onset of the WF coupling. Therefore, f(η, r, x) shows a deep valley around x = 0 (or ν₀), while j(η, r, x) shows a plateau.

In Figure 8, we present the solutions of mean intensity j(η, r, x) and flux f(η, r, x) given by Equations (16) and (17) with the Voigt profile Equation (4) and re-distribution function Equation (7). As expected, Figure 8 has the Lorentz wings. However, in the center part |x| ⩽ 3, Figure 8 shows the same features as Figures 4 and 6. That is, the mechanism of "escape via shortcut" plus "bounce back," which mainly relies on photons with |x| < 3, still works well. The Lorentz wing only leads to long tails in the profiles of j and f in the range |x|>3, and the wings have very low amplitudes.

Zoom In Zoom Out Reset image size

Hubble redshift is another mechanism to produce photons with frequency ∼ν₀, which also have the problems of the spatial transfer and the WF coupling. Since the 21 cm region basically is optically thin for photons with x > 3, we first model the cosmic redshift by a photon source with the profile ϕ_s(x) = ϕ_g(x − 3). The solutions of f(η, r, x) and j(η, r, x) with the same parameters as in Figure 7 are shown in Figure 9.

Zoom In Zoom Out Reset image size

Although the injected photon has a peak at x = 3, the flux f(η, r, x) in Figure 9 still shows two peaks at x ≃ ±(2 − 3). The peak at x ≃ (2 − 3) is higher than that at −x ≃ (2 − 3). It is simply because the photon source is at x = 3. The peak at x < 0 is also much higher than f(η, r, x) at x = 0. That is, even though the source photon is at x = 3, the spatial transfer is still dominated by photons of both x ≃ (2 − 3) and x ≃ −(2 − 3). The mean intensity j(η, r, x) shows once again a perfect local Boltzmann distribution with slope 2b in the range −2 < x < 2 (Equation (23)). Any photons injected by the redshift into the 21 cm region will quickly join the WF coupling by the bounce back mechanism.

We now model the source with a continuous frequency spectrum: ϕ_s(x) = S₀ within |x| ⩽ 10, and ϕ_s(x) = 0 at |x|>10. The solutions of the mean intensity j(η, r, x) and flux f(η, r, x) of Equations (16) and (17) are shown in Figure 10, in which the other parameters are the same as those in Figure 8. We can see that in the center part |x| ⩽ 3, j(η, r, x) and f(η, r, x) have the same features as in Figure 8. This result is expected, as the center part has contributions from the photons redshifted to |x| ⩽ 3 from x > 3. Figure 9 shows that the source of ϕ_g(x − 3) does not change the mechanism of "escape via shortcut" plus "bounce back," and therefore, the source with a continuous frequency spectrum will keep these features as well.

Zoom In Zoom Out Reset image size

Although the cosmic expansion is considered in all the above-mentioned solutions, the effect of cosmic expansion seems to be negligible. It is because the optical depth is large and the parameter γ is very small. The number of resonant scattering within the Hubble time is very large. The cosmic expansion is too small in one free flight time. The "bounce back" is dominant.

When f_H i is smaller, optical depth of the ν₀ photons is smaller, and γ is larger, the effect of the Hubble redshift would appear. Figure 11 presents the solution f(η, r, x) and j(η, r, x) of Equations (16) and (17) with the same parameters as those in Figures 4 and 6, except that the parameter γ = 10⁻³ is large. Both j(η, r, x) and f(η, r, x) show the same features as in Figures 4 and 6. The Hubble redshift makes the profile to be asymmetric with respect to x = 0. The red wing is stronger than the blue wing. j(η, r, x) still shows a flat plateau in the range −2 < x < 2. Therefore, the WF coupling will work when γ = 10⁻³, corresponding to f_H i ≃ 10⁻³.

Zoom In Zoom Out Reset image size

The spatial derivative terms in Equations (16) and (17) are approximated by a fifth-order finite difference WENO scheme.

We first give the WENO reconstruction procedure in approximating $\frac{\partial j}{\partial x}$,

$$\begin{equation} \frac{\partial j(\eta ^n,r_i,x_j)}{\partial x} \approx \frac{1}{\Delta x} (\hat{h}_{j+1/2}-\hat{h}_{j-1/2}), \end{equation} \tag{ A1 }$$

with fixed η = ηⁿ and r = r_i. The numerical flux $\hat{h}_{j+1/2}$ is obtained by the fifth-order WENO approximation in an upwind fashion, because the wind direction is fixed (negative). Denote

$$\begin{equation} h_j = j(\eta ^n,r_i,x_j), \qquad j = -2, -1,\ldots,N+3 \end{equation} \tag{ A2 }$$

with fixed n and i. The numerical flux from the WENO procedure is obtained by

$$\begin{equation} \hat{h}_{j+1/2}=\omega _1\hat{h}_{j+1/2}^{(1)}+\omega _2\hat{h}_{j+1/2}^{(2)} + \omega _3\hat{h}_{j+1/2}^{(3)}, \end{equation} \tag{ A3 }$$

where $\hat{h}_{j+1/2}^{(m)}$ are the three third-order fluxes on three different stencils given by

$$\begin{eqnarray*} \hat{h}_{j+1/2}^{(1)} &=& -\!\frac{1}{6}h_{j-1}+\frac{5}{6}h_{j}+\frac{1}{3}h_{j+1},\\ \hat{h}_{j+1/2}^{(2)} &=& \frac{1}{3}h_{j}+\frac{5}{6}h_{j+1}-\frac{1}{6}h_{j+2},\\ \hat{h}_{j+1/2}^{(3)} &=& \frac{11}{6}h_{j+1}-\frac{7}{6}h_{j+2}+\frac{1}{3}h_{j+3}, \end{eqnarray*}$$

and the nonlinear weights ω_m are given by

$$\begin{equation} \omega _m = \frac{\check{\omega }_m}{\sum _{l=1}^3\check{\omega }_l}, \qquad \check{\omega }_l = \frac{\gamma _l}{(\epsilon +\beta _l)^2}, \end{equation} \tag{ A4 }$$

where is a parameter to avoid the denominator to become zero and is taken as = 10⁻⁸. The linear weights γ_l are given by

$$\begin{equation} \gamma _{1} = \frac{3}{10},\qquad \gamma _{2} = \frac{3}{5}, \qquad \gamma _{3} = \frac{1}{10}, \end{equation} \tag{ A5 }$$

and the smoothness indicators β_l are given by

$$\begin{eqnarray*} \beta _1 &=& \frac{13}{12}(h_{j-1}-2h_{j}+h_{j+1})^2 +\frac{1}{4}(h_{j-1}-4h_{j}+3h_{j+1})^2,\\ \beta _2 &=& \frac{13}{12}(h_{j}-2h_{j+1}+h_{j+2})^2 +\frac{1}{4}(h_{j}-h_{j+2})^2,\\ \beta _3 &=& \frac{13}{12}(h_{j+1}-2h_{j+2}+h_{j+3})^2 +\frac{1}{4}(3h_{j+1}-4h_{j+2}+h_{j+3})^2. \end{eqnarray*}$$

To approximate the r-derivatives in the system of Equations (16) and (17), we need to perform the WENO procedure based on a characteristic decomposition. We write the left-hand side of Equations (16) and (17) as

$$\begin{equation} {\bf u}_t + A {\bf u}_r, \end{equation} \tag{ A6 }$$

where u = (j, f)^T and

$$\begin{eqnarray*} A =\left(\begin{array}{@{}cc@{}} 0 & 1 \\ \frac{1}{3} & 0\end{array} \right) \end{eqnarray*}$$

is a constant matrix. To perform the characteristic decomposition, we first compute the eigenvalues, the right eigenvectors, and the left eigenvectors of A and denote them by, Λ, R and R⁻¹. We then project u to the local characteristic fields v with v = R⁻¹u. Now u_t + Au_r of the original system is decoupled as two independent equations as v_t + Λv_r. We approximate the derivative v_r component by component, each with the correct upwind direction, with the WENO reconstruction procedure similar to the procedure described above for $\frac{\partial j}{\partial x}$. In the end, we transform v_r back to the physical space by u_r = Rv_r. We refer the readers to Cockburn et al. 1998 for more implementation details.

We apply the rectangular rule to evaluate ∫R(x, x')jdx. The rectangular rule is known to have spectral accuracy for smooth integrated function R(x, x')j(x') with compact support. In Equations (16) and (17), it is known that

$$\begin{equation} \int R(x, x^{\prime })dx^{\prime } = \phi (x), \end{equation} \tag{ A7 }$$

which is crucial for photon conservation. However,

$$\begin{eqnarray*} &&\sum _{j} R(x, x_j) \Delta x = \phi (x) \end{eqnarray*}$$

may not be true in general. To numerically preserve the photon conservation, we first compute

$$\begin{equation} {\rm ratio}(x, \Delta x) = \frac{\phi (x)}{ \sum _j R(x, x_j) \Delta x}, \end{equation} \tag{ A8 }$$

then we normalize the collision term by approximating ∫R(x, x')j(t, r, x')dx' with

$$\begin{equation} {\rm ratio}(x, \Delta x) \, \sum _j R(x, x_j) j(\eta,r,x_j) \Delta x \end{equation} \tag{ A9 }$$

with fixed η and r.

The source term given in Equation (16) is implemented as a boundary condition on f(η, r = r₀, x),

$$\begin{eqnarray*} &&f(\eta,r=r_0,x) = s_0\phi _s(x). \end{eqnarray*}$$

For the intensity j, a reflective boundary condition is used at r = r₀. At the boundary of r = r_max, x = x_left, and x = x_right, we use zero boundary conditions for both j and f, because of the way we choose r_max, x_left, and x_right.

The time derivatives $\frac{\partial j}{\partial \eta }$ and $\frac{\partial f}{\partial \eta }$ are approximated by the third-order TVD Runge–Kutta time discretization (Shu & Osher 1988). For systems of ODEs u_t = L(u), the third-order Runge–Kutta method is

$$\begin{eqnarray*} u^{(1)} &=& u^n + \Delta t L(u^n,t^n),\\ u^{(2)} &=& \frac{3}{4}u^n + \frac{1}{4}(u^{(1)}+\Delta t L(u^{(1)},t^n + \Delta t)),\\ u^{n+1} &=& \frac{1}{3}u^n + \frac{2}{3}\left(u^{(2)}+\Delta t L\left(u^{(2)},t^n + \frac{1}{2} \Delta t\right)\right). \end{eqnarray*}$$

From Equation (16) we have

$$\begin{equation} \frac{\partial }{\partial \eta } \int j dx + \frac{\partial }{\partial r} \int f dx = 0. \end{equation} \tag{ A10 }$$

Therefore, for a time-independent solution we have ∫f(r, x)dx = const. It yields

$$\begin{equation} \int f(r,x)dx = \int f(0,x)dx= S_0. \end{equation} \tag{ A11 }$$

That is, the flux at saturated states is r-independent. This is the conservation of the number of photons. It can be used to test the algorithm.