Ly$\alpha$ Radiative Transfer: Monte-Carlo Simulation of the Wouthuysen-Field Effect

In the code, we generate photon packets and track each of them in physical and frequency spaces. The Monte-Carlo method for the scattering of Ly$\alpha$ photons by neutral hydrogen atoms proceeds through the following steps. First, we generate a position $\mathbf{r}$, a frequency $\nu$, and a propagation direction $\mathbf{k}_{i}$ of the photon packet. The initial location and frequency of the photon are drawn from the source’s spatial distribution function and an intrinsic frequency distribution under consideration, respectively. The initial propagation direction is generated to be isotropic. Second, the optical depth $\tau$ that the photon will travel through before it interacts with a hydrogen atom is randomly drawn from the probability density function $P(\tau)=e^{-\tau}$ as follows:

$$\tau=-\ln\xi,$$

(1)

where $\xi$ is drawn from a uniform random distribution. The physical distance is then calculated by inverting the integral of optical depth over the photon pathlength, defined in Equation (2). Third, we generate the velocity vector of the atom that scatters the photon, as described later. Fourth, we draw the scattering angle of the photon from an appropriate scattering phase function and calculate the new propagation direction $\mathbf{k}_{f}$. The velocity of the scattering atom and the old and new direction vectors of the photon determine a new frequency following the energy and momentum conservation laws. The process of propagation and scattering is repeated until the photon escapes the system. For the present study, we used a number of photon packets of $10^{5}$ up to $10^{10}$.

The optical depth $\tau_{\nu}(s)$ of a photon with frequency $\nu$ along a path length $s$ is given by

$$\tau_{\nu}(s)=\int_{0}^{s}\int_{-\infty}^{\infty}n(v_{\parallel})\sigma_{\nu}dv_{\parallel}dl,$$

(2)

where $n(v_{\parallel})$ is the number density of neutral hydrogen atom with velocity component $v_{\parallel}$ parallel to the photon’s traveling path and $\sigma_{\nu}$ is the scattering cross-section. It is convenient to express frequency in terms of the relative frequency, $x\equiv(\nu-\nu_{\alpha})/\Delta\nu_{{\rm D}}$, where $\nu_{\alpha}=2.466\times 10^{15}$ Hz is the central frequency of Ly$\alpha$, $\Delta\nu_{{\rm D}}=\nu_{\alpha}(v_{{\rm th}}/c)=1.057\times 10^{11}(T_{{\rm K}}/10^{4}{\rm K})^{1/2}$ Hz is the thermal, Doppler frequency width, and $v_{{\rm th}}=(2k_{{\rm B}}T_{{\rm K}}/m_{{\rm H}})^{1/2}=12.85(T_{{\rm K}}/10^{4}{\rm K})^{1/2}$ km s^-1 is the thermal velocity dispersion, multiplied by $\sqrt{2}$, of atomic hydrogen gas with a kinetic temperature $T_{{\rm K}}$. Here, $k_{{\rm B}}$ is the Boltzmann constant and $m_{{\rm H}}$ is the hydrogen atomic mass. For the ray tracing of photon packets in evaluating the optical depth integral of Equation (2), we use the fast voxel traversal algorithm (Amanatides & Woo, 1987), which is used in our dust radiative transfer code MoCafe (Seon & Witt, 2012; Seon, 2015; Seon & Draine, 2016).

The scattering cross-section of a Ly$\alpha$ photon in the rest frame of a hydrogen atom is given by

$$\sigma_{\nu}^{{\rm rest}}=f_{\alpha}\frac{\pi e^{2}}{m_{e}c}\phi_{\nu}^{{\rm rest}}=f_{\alpha}\frac{\pi e^{2}}{m_{e}c}\frac{\Gamma/4\pi^{2}}{(\nu-\nu_{\alpha})^{2}+(\Gamma/4\pi)^{2}},$$

(3)

where $f_{\alpha}=0.4162$ is the oscillator strength of Ly$\alpha$, and $\Gamma(=A_{\alpha})$ is the damping constant. The resonance profile¹¹1The normalized profile of the resonance cross-section as a function of frequency is often referred to as the line profile. In this paper, we call it “the resonance profile” to distinguish it from the line profile arising as a result of the RT process. $\phi_{\nu}^{{\rm rest}}$ is normalized such that $\int_{0}^{\infty}\phi_{\nu}^{{\rm rest}}d\nu=1$. The Einstein A coefficient of the Ly$\alpha$ transition between the $n=2$ and $n=1$ levels of a hydrogen atom is $A_{\alpha}=6.265\times 10^{8}$ s^-1. Integrating over the Maxwellian velocity distribution of hydrogen atoms, the scattering cross-section in the laboratory frame is found to be

$$\sigma_{\nu}=\frac{\chi_{0}}{\Delta\nu_{{\rm D}}}\frac{H(x,a)}{\sqrt{\pi}}=\frac{\chi_{0}}{\Delta\nu_{{\rm D}}}\phi_{x},$$

(4)

where $\chi_{0}\equiv f_{\alpha}\pi e^{2}/m_{e}c=\left(3c^{2}/8\pi\nu_{\alpha}^{2}\right)A_{\alpha}=1.105\times 10^{-2}$ cm² Hz, $\phi_{x}\equiv H(x,a)/\sqrt{\pi}$, and

$$H(x,a)=\frac{a}{\pi}\int_{-\infty}^{\infty}\frac{e^{-u^{2}}}{(x-u)^{2}+a^{2}}du$$

(5)

is the Voigt-Hjerting function. Here, $a=\Gamma/(4\pi\Delta\nu_{{\rm D}})=4.717\times 10^{-4}\ (T_{{\rm K}}/10^{4}{\rm K})^{-1/2}$ is the natural linewidth relative to the thermal frequency width. The atomic velocity v is normalized by the thermal speed, i.e., $\mathbf{u}=\textbf{v}/v_{{\rm th}}$. We also note that the resonance profile $\phi_{x}$ is normalized, e.g., $\int_{-\infty}^{\infty}\phi_{x}dx=1$. In this paper, the monochromatic optical thickness $\tau_{0}\equiv(\chi_{0}/\Delta\nu_{{\rm D}})N_{{\rm HI}}\phi_{x}(0)\simeq(\chi_{0}/\sqrt{\pi}\Delta\nu_{{\rm D}})N_{{\rm HI}}$ measured at the line center is used to represent the total optical thickness of media. Here, $N_{{\rm HI}}$ is the column density of neutral hydrogen atoms. The optical depth has a temperature dependence of

$\displaystyle\tau_{0}$

$\displaystyle\simeq$

$\displaystyle 5.90\times 10^{6}\left(N_{{\rm HI}}/10^{20}\,{\rm cm}^{-2}\right)\left(T_{{\rm K}}/10^{4}\,K\right)^{-1/2}.$

(6)

The monochromatic optical thickness ($\tau_{0}$) is related to the optical thickness integrated over the relative frequency ($\tau_{*}\equiv\int_{-\infty}^{\infty}\tau_{x}dx$), which is used in Neufeld (1990), by

$$\tau_{0}\simeq\frac{\tau_{*}}{\sqrt{\pi}}\left(1-\frac{2a}{\sqrt{\pi}}+a^{2}-\cdots\right).$$

(7)

We need an accurate yet fast algorithm to calculate the Voigt function because the function is repeatedly estimated in simulations. An approximate formula suggested by N. Gnedin, as described in Tasitsiomi (2006) and Laursen et al. (2009), is simple and thus might be appropriate in calculating the Ly$\alpha$ spectrum emerging from the hydrogen gas. However, we found that the approximation yields a systematic, spurious wavy feature at the line center of the Ly$\alpha$ spectrum calculated inside the medium. We, instead, adopt a slightly modified version of the function “voigt_king” in the VPFIT program (Carswell & Webb, 2014). This routine was found to be accurate and as fast as or even quicker than that given in Tasitsiomi (2006).

The velocity of the scattering atom is decomposed of the parallel ($u_{\parallel}$) and perpendicular ($\mathbf{u}_{\perp}$) components with respect to the photon’s propagation direction $\mathbf{k}$. The two (mutually orthogonal) perpendicular components $\mathbf{u}_{\perp}$ are independent of the photon frequency and thus are drawn from the following two-dimensional (2D) Gaussian distribution with zero mean and standard deviation $1/\sqrt{2}$:

$$f(\mathbf{u}_{\perp})=\frac{1}{\sqrt{\pi}}e^{-|\mathbf{u}_{\perp}|^{2}}.$$

(8)

The two random velocity components are then generated as

	$\displaystyle u_{\perp,1}$	$\displaystyle=$	$\displaystyle(-\ln\xi_{1})^{1/2}\cos(2\pi\xi_{2})$
	$\displaystyle u_{\perp,2}$	$\displaystyle=$	$\displaystyle(-\ln\xi_{1})^{1/2}\sin(2\pi\xi_{2}),$		(9)

where $\xi_{1}$ and $\xi_{2}$ are two independent random variates between 0 and 1. When a photon is in the core of the Voigt profile, it is trapped by resonance scattering in a small region, where it started, and thus cannot travel long distance until it is scattered off to the wing of the profile. Therefore, to speed up the calculation, one can use an acceleration scheme, known as “core-skipping,” which was first suggested by Ahn et al. (2002). The scheme artificially pushes the photon into the wing by using a truncated Gaussian distribution. This is achieved by generating the perpendicular components as

	$\displaystyle u_{\perp,1}$	$\displaystyle=$	$\displaystyle(x_{{\rm crit}}^{2}-\ln\xi_{1})^{1/2}\cos(2\pi\xi_{2})$
	$\displaystyle u_{\perp,2}$	$\displaystyle=$	$\displaystyle(x_{{\rm crit}}^{2}-\ln\xi_{1})^{1/2}\sin(2\pi\xi_{2}),$		(10)

in which the critical frequency $x_{{\rm crit}}$ can be chosen as described in Dijkstra et al. (2006), Laursen et al. (2009), or Smith et al. (2017). This acceleration scheme is also implemented in LaRT. However, the scheme is not used for the present study because it gives rise to a serious underestimation of the number of scatterings. As we shall show later, measuring the number of scatterings is the most crucial point to study the WF effect. We also attempted to apply various bias techniques, such as the composite bias method (Baes et al., 2016) or others (Juvela, 2005), which have been developed for dust and nuclear physics, to improve the computation speed. However, none of the methods was found to be appropriate for the Ly$\alpha$ RT. This failure is because the number of scatterings in the Ly$\alpha$ RT is exceptionally high compared to that for dust.

The parallel component $u_{\parallel}$ of the scattering atom depends on the incident photon frequency and thus should be drawn from the following distribution function:

$$f(u_{\parallel}|x,a)=\frac{a}{\pi H(x,a)}\frac{e^{-u_{\parallel}^{2}}}{a^{2}+(x-u_{\parallel})^{2}}.$$

(11)

Note that this probability distribution function depends on not only the photon frequency but also the gas temperature through $a$. The most commonly adopted algorithm is based on the acceptance/rejection method of Zheng & Miralda-Escudé (2002a). Detailed descriptions are given in Semelin et al. (2007) and Laursen et al. (2009). The algorithm is elaborated by Smith et al. (2017). For LaRT, we developed a novel method to generate random $u_{\parallel}$ using the ratio-of-uniforms method (Kinderman & Monahan, 1977), as described in Appendix 28. Our algorithm is efficient in the sense that fewer random numbers are rejected, especially for high $x$ values and high temperatures $T_{{\rm K}}\gtrsim 10^{3}$ K.

After the velocity components of the scattering atom is selected, the scattering direction is generated assuming the Rayleigh phase function $P(\theta)\propto 1+\cos^{2}\theta$, where $\theta$ is the angle between the incident direction $\mathbf{k}_{i}$ and the outgoing direction $\mathbf{k}_{f}$. The random, scattering angle $\theta$ is determined as described in Seon (2006):

	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle\left[2(2\xi-1)+\sqrt{4(2\xi-1)^{2}+1}\right]^{1/3}$		(12)
			$\displaystyle-\left[2(2\xi-1)+\sqrt{4(2\xi-1)^{2}+1}\right]^{-1/3},$

where $\mu=\cos\theta$ and $\xi$ is a uniform random variate. The azimuthal angle $\phi$ is generated by $\phi=2\pi\xi^{\prime}$ using another random variate $\xi^{\prime}$. More generally, the scattering phase function is known to have a form of $P(\theta)\propto 1+(R/Q)\cos^{2}\theta$, in which $R/Q$ is the degree of polarization for 90^∘ scattering and depends on the incident frequency (Stenflo, 1980). We will ignore the dependence of the scattering angle distribution on the photon frequency, which is not significant for the present purpose. The algorithm to draw the scattering angle $\theta$ for this phase function will be described in a forthcoming paper, in which the Ly$\alpha$ polarization is investigated.

We assumed, unless otherwise stated, that photons are injected at a monochromatic frequency of $x=0$. If the initial photon packets are supposed to have the Voigt profile²²2The emission line profile is often regarded to be a Gaussian function. However, the ‘intrinsic’ profile of an emission line is a convolution of Lorentzian and Gaussian functions, representing the natural line profile from an undriven harmonically bound atom and the effect due to the gas thermal motion, respectively (Rybicki & Lightman, 1985; Hubeny & Mihalas, 2014). The ‘intrinsic’ absorption profile has the same shape as the ‘intrinsic’ emission profile by the principle of detailed balance (Hubeny & Mihalas, 2014)., its frequency is obtained by adding a random number for a Lorentzian distribution and a random number for a Gaussian distribution:

$$x=a\tan\left[\pi\left(\xi-\frac{1}{2}\right)\right]+\frac{\xi_{{\rm G}}}{\sqrt{2}},$$

(13)

where $\xi$ is a uniform random number between 0 and 1, and $\xi_{{\rm G}}$ is a random number drawn from a Gaussian distribution with zero mean and unit variance. Here, we note that $\xi_{{\rm G}}/\sqrt{2}$ can be obtained, for instance, using Equation (9).

In the code, each cell has a bulk velocity ($\mathbf{u}^{{\rm bulk}}=\textbf{v}^{{\rm bulk}}/v_{{\rm th}}$), which is expressed in units of the thermal speed that depends on the gas temperature. When a photon packet traverses from cell to cell (for instance, $i$ to $j$), the gas temperature $T_{{\rm K}}$ and the bulk velocity $\mathbf{u}^{{\rm bulk}}$ may vary. The relative frequency should then be changed to

$$x_{j}=\left(x_{i}+\mathbf{k}\cdot\mathbf{u}_{i}^{{\rm bulk}}\right)\left(\Delta\nu_{{\rm D}}^{i}/\Delta\nu_{{\rm D}}^{j}\right)^{1/2}-\mathbf{k}\cdot\mathbf{u}_{j}^{{\rm bulk}},$$

(14)

as the photon crosses the cell $i$ to $j$, where $\mathbf{k}$ is the photon’s propagation vector with unit length.

After a scattering event, the incident photon frequency $x_{i}$ is changed to the outgoing frequency $x_{f}$ given by

$$x_{f}=x_{i}-u_{\parallel}+\mathbf{k}_{f}\cdot\mathbf{u}_{{\rm atom}}-g_{*}\left(1-\mathbf{k}_{f}\cdot\mathbf{k}_{i}\right),$$

(15)

where $g_{*}=(h\nu_{\alpha}^{2}/m_{{\rm H}}c^{2})/\Delta\nu_{{\rm D}}=2.536\times 10^{-4}(T_{{\rm K}}/10^{4}{\rm K})^{-1/2}\simeq 0.54a$ is the recoil parameter. Here, $h$ is the Planck constant. The recoil parameter might be negligible for applications in which only the spectra emergent from galaxies are of interest. However, as we will show, it is this small effect that causes the WF effect.

In addition to being scattered by neutral hydrogen atoms, Ly$\alpha$ photons also interact with dust grains; they can either be scattered or destroyed by dust. When the interaction with dust is included, the total optical depth along a pathlength $s$ is the sum of optical depths due to neutral hydrogen and dust, $\tau=\tau_{{\rm H}}+\tau_{{\rm dust}}$. The absorption optical depth $\tau_{{\rm abs}}$ by dust is defined to be $\tau_{{\rm abs}}=(1-\eta)\tau_{{\rm dust}}$ for the dust extinction optical depth $\tau_{{\rm dust}}$ and the scattering albedo $\eta$.

The probability of being interacted with dust is given by

$$\mathcal{P}_{{\rm dust}}=\frac{n_{{\rm dust}}\sigma_{{\rm dust}}}{n_{{\rm H}}\sigma_{{\rm H}}+n_{{\rm dust}}\sigma_{{\rm dust}}}.$$

(16)

If a random number $\xi$ is generated to be less than $\mathcal{P}_{{\rm dust}}$, then the photon interacts with a dust grain; otherwise, it is scattered by a hydrogen atom.

To take into account the interaction with dust grains, we implemented two different methods in LaRT. The first method is to use the weight of photon packets, where the weight is initially set to unity ($w=1$) and reduced by a factor of albedo ($w^{\prime}=\eta w$) as the photon packet undergoes dust scattering. The second method is to compare a uniform random number with the dust albedo and then assume the photon packet will be absorbed and destroyed by a dust grain if the random number is larger than the albedo.

When the photon is scattered by dust, the scattering angle $\theta$ is found using the inversion method to follow the Henyey-Greenstein phase function for a given asymmetry factor $g$ (e.g., Witt, 1977):

$$\mu=\frac{1+g^{2}}{2g}-\frac{1}{2g}\left(\frac{1-g^{2}}{1-g+2g\xi}\right)^{2},$$

(17)

where $\xi$ is drawn from a uniform random distribution.

We adopt the Milky-Way dust model of Weingartner & Draine (2001) and Draine (2003). The dust extinction cross-section per hydrogen atom is $n_{{\rm dust}}\sigma_{{\rm dust}}/n_{{\rm H}}=1.61\times 10^{-21}$ cm²/H at Ly$\alpha$. The asymmetry factor and albedo at Ly$\alpha$ are given by $g=0.676$ and $\eta=0.325$, respectively. The dust parameters $(n_{{\rm dust}}\sigma_{{\rm dust}}/n_{{\rm H}},\ \eta,\ g)$ for the other dust types are ($5.32\times 10^{-22}$ cm²/H, 0.255, 0.648) for the Large Magellanic Cloud, and ($3.76\times 10^{-22}$ cm²/H, 0.330, 0.591) for the Small Magellanic Cloud.

The hyperfine structure of the $n=2$ state should be considered in order to derive the formula for the spin temperature due to the indirect Ly$\alpha$ pumping of the hyperfine levels in the ground state (the lower level “0” and the upper level “1” in Figure 33). However, one can ignore the hyperfine splitting when calculating only the Ly$\alpha$ spectrum through the RT process because the Ly$\alpha$ line profile is much broader than the splitting. The fine structure of the $n=2$ state also can be ignored in most of the cases presented in this paper. At low gas temperatures and low optical thicknesses, however, the fine-structure splitting of the $n=2$ state must be taken into account; in these circumstances, the line width is comparable to or smaller than the frequency gap between the fine-structure levels of the excited state. The level splitting of 10.8 GHz in the excited state is equivalent to a Doppler shift of $1.34$ km s^-1; the Doppler width ($\Delta\nu_{{\rm D}}\lesssim 1.28$ km s^-1) for the gas of $T\lesssim 10^{2}$ K is thus smaller than the fine-structure splitting.

We use the term $\nicefrac{{1}}{{2}}$ to denote the transition $1\,{}^{2}{\rm S}_{1/2}\leftrightarrow 2\,{}^{2}{\rm P}_{1/2}^{{\rm o}}$ and $\nicefrac{{3}}{{2}}$ the transition $1\,{}^{2}{\rm S}_{1/2}\leftrightarrow 2\,{}^{2}{\rm P}_{3/2}^{{\rm o}}$. Then, the cross-section is given by

$\displaystyle\sigma_{\nu}^{{\rm rest}}$

$\displaystyle=$

$\displaystyle\chi_{0}\left[\frac{(1/3)\Gamma/4\pi^{2}}{(\nu-\nu_{\nicefrac{{1}}{{2}}})^{2}+(\Gamma/4\pi)^{2}}+\frac{(2/3)\Gamma/4\pi^{2}}{(\nu-\nu_{\nicefrac{{3}}{{2}}})^{2}+(\Gamma/4\pi)^{2}}\right]$

(18)

in the rest frame of a hydrogen atom (e.g., Ahn & Lee, 2015). The cross-section is a combination of two Lorentzian functions with weights of 1:2. In a gas medium having a Maxwellian velocity distribution at temperature $T$, the scattering cross-section is the sum of two Voigt functions:

$$\sigma_{\nu}=\frac{\chi_{0}}{\Delta\nu_{{\rm D}}}\left[\frac{1}{3}\phi(x_{\nicefrac{{1}}{{2}}})+\frac{2}{3}\phi(x_{\nicefrac{{3}}{{2}}})\right],$$

(19)

where $x_{\nicefrac{{1}}{{2}}}\equiv(\nu-\nu_{\nicefrac{{1}}{{2}}})/\Delta\nu_{{\rm D}}$ and $x_{\nicefrac{{3}}{{2}}}\equiv(\nu-\nu_{\nicefrac{{3}}{{2}}})/\Delta\nu_{{\rm D}}$. We also define $x\equiv\left[\nu-(\nu_{\nicefrac{{1}}{{2}}}+\nu_{\nicefrac{{3}}{{2}}})/2\right]/\Delta\nu_{{\rm D}}=(x_{\nicefrac{{1}}{{2}}}+x_{\nicefrac{{3}}{{2}}})/2$. The parallel component $u_{\parallel}$ of scattering atom is drawn from the following distribution function:

$\displaystyle f_{{\rm FS}}(u_{\parallel}|x)$

$\displaystyle=$

$\displaystyle\mathcal{P}_{\nicefrac{{1}}{{2}}}f(u_{\parallel}|x_{\nicefrac{{1}}{{2}}})+(1-\mathcal{P}_{\nicefrac{{1}}{{2}}})f(u_{\parallel}|x_{\nicefrac{{3}}{{2}}}),$

(20)

where

$$\mathcal{P}_{\nicefrac{{1}}{{2}}}\equiv\frac{H(x_{\nicefrac{{1}}{{2}}},a)}{H(x_{\nicefrac{{1}}{{2}}},a)+2H(x_{\nicefrac{{3}}{{2}}},a)}.$$

(21)

The random variates for this distribution are obtained using the composition method. If a uniform random number $\xi$ is selected to be smaller than the probability $\mathcal{P}_{\nicefrac{{1}}{{2}}}$, the photon is scattered by the $\nicefrac{{1}}{{2}}$ transition; otherwise, by the $\nicefrac{{3}}{{2}}$ transition. Then, either $f(u_{\parallel}|x_{\nicefrac{{1}}{{2}}})$ or $f(u_{\parallel}|x_{\nicefrac{{3}}{{2}}})$, depending on the transition, determines the parallel velocity component. The initial photon frequency is chosen to be either $\nu=\nu_{\nicefrac{{1}}{{2}}}$ or $\nu=\nu_{\nicefrac{{3}}{{2}}}$ with a probability ratio of 1:2. We note, in an optically thick medium, that the line ratio becomes equal to 1:1 as a result of the resonance-line scattering (e.g., Long et al., 1992).

In this paper, we present the results for $T_{{\rm K}}=10$ K and $10^{4}$ K, unless otherwise stated. We consider the fine-structure splitting only when calculating the Ly$\alpha$ spectrum inside the medium of $T_{{\rm K}}=10$ K and occasionally of $T_{{\rm K}}=10^{2}$ K in Section 4; the fine-structure splitting is negligible at $T_{{\rm K}}\gtrsim 10^{3}$ K. In Sections 3 and 5, the splitting effect is not taken into account. A more complete description of how to deal with the fine-structure splitting in Ly$\alpha$ RT will be given in a separate paper.

Field (1958) found that the 21-cm spin temperature is a weighted mean of the three temperatures $T_{{\rm R}}$, $T_{{\rm K}}$, and $T_{\alpha}$, as follows:

$\displaystyle T_{{\rm s}}$

$\displaystyle=$

$\displaystyle\frac{T_{{\rm R}}+y_{{\rm c}}T_{{\rm K}}+y_{\alpha}T_{\alpha}}{1+y_{{\rm c}}+y_{\alpha}},$

(22)

where $T_{{\rm R}}$ is the brightness temperature of the background radio radiation at $21$ cm (e.g., $T_{{\rm R}}\sim 2.7$ K for the CMB), and $T_{{\rm K}}$ is the kinetic temperature of hydrogen gas. The color temperature $T_{\alpha}$ of Ly$\alpha$ is defined by the exponential slope of the mean intensity spectrum at the line center, as follows:

$$J_{\nu}=J_{\nu}(\nu_{\alpha})\exp\left[-\frac{h\left(\nu-\nu_{\alpha}\right)}{k_{{\rm B}}T_{\alpha}}\right].$$

(23)

The weighting factors in Equation (22) are

$\displaystyle y_{{\rm c}}$

$\displaystyle=$

$\displaystyle\frac{T_{*}}{T_{{\rm K}}}\frac{P_{10}^{{\rm c}}}{A_{10}}\ \ {\rm and}\ \ y_{\alpha}=\frac{4}{27}\frac{T_{*}}{T_{\alpha}}\frac{P_{\alpha}}{A_{10}},$

(24)

where $T_{*}\equiv h\nu_{10}/k_{{\rm B}}=0.0681$ K, $P_{10}^{{\rm c}}$ is the downward transition rate (per sec) from the hyperfine level “1” to “0” due to particle collisions (see Figure 33 for the level definition), and $P_{\alpha}$ is the scattering rate (the number of Ly$\alpha$ scatterings per unit time for a single hydrogen atom). Here, $A_{10}$ is the Einstein coefficient for the spontaneous emission of 21-cm line. In Appendix 33, we describe the above equations in more detail.

From the above Equations (22)-(24), we can identify two necessary conditions for the WF effect. First, the central part of the Ly$\alpha$ spectrum in the medium should be the exponential shape of Equation (23) with a color temperature that is equal to the kinetic temperature (e.g., $T_{\alpha}=T_{{\rm K}}$). Second, the weighting factor for the Ly$\alpha$ pumping should be much higher than that of the background radiation (e.g., $y_{\alpha}\gg 1$). We note that the first condition requires individual Ly$\alpha$ photons to experience a sufficient number of resonance-line scatterings to thermalize the spectrum via atomic recoil. On the other hand, the second condition is primarily a matter of the amount of Ly$\alpha$ photons that are supplied to the medium. The amount of dust grains and the kinematics of gas are also critical factors that affect both requirements.

To study the condition for the first requirement of the WF effect, we first need to examine the spectral shape of Ly$\alpha$ in the medium. To calculate the spectrum of the Ly$\alpha$ radiation field in a volume element of the medium, we use the method that is described in Lucy (1999) and has been extensively adopted to estimate the absorbed energy of stellar radiation by dust grains and the re-emitted infrared emission in dust RT simulations (e.g., Misselt et al., 2001; Baes et al., 2011).

Let $\delta\ell_{i}$ denote the pathlength between successive events and $\Delta t$ the total duration of the “sequential” simulation. Here, events indicate not only scatterings but also crossings across boundaries between volume elements. The energy carried by a single photon packet is denoted by $\epsilon_{*}=L_{*}\Delta t/N_{{\rm packet}}$, where $L_{*}$ and $N_{{\rm packet}}$ are the total luminosity and the total number of photon packets used in the simulation, respectively. The segment of the packet’s trajectory between two successive events contributes $\epsilon_{*}\delta t_{i}/\Delta t$ to the time-averaged energy content of a volume element, where $\delta t_{i}=\delta\ell_{i}/c$ is the traveling time between events. Here, we note that $\Delta t=\sum_{k=1}^{N_{{\rm packet}}}\sum_{i}\delta t_{i}$. The mean intensity $J_{\nu}$ of the radiation field, expressed in the unit of photon number, in the frequency interval $(\nu,\ \nu+\Delta\nu)$ can be expressed in terms of the energy density $u_{\nu}$ by $h\nu_{\alpha}J_{\nu}=(c/4\pi)u_{\nu}$.³³3In this paper, $J_{\nu}$ is used to denote the mean intensity in terms of photon number. The mean intensity expressed in terms of energy is denoted by $I_{\nu}$. They are related by $I_{\nu}=h\nu_{\alpha}J_{\nu}$. The mean intensity is then given by summing all energy segments that contribute to the volume element:

	$\displaystyle h\nu_{\alpha}J_{\nu}\Delta\nu$	$\displaystyle=$	$\displaystyle\frac{c}{4\pi}\frac{1}{\Delta V}\sum_{i}\epsilon_{*}\frac{\delta t_{i}}{\Delta t},$
	$\displaystyle J_{\nu}$	$\displaystyle=$	$\displaystyle\frac{1}{4\pi\Delta V\Delta\nu}\frac{L_{*}/h\nu_{\alpha}}{N_{{\rm packet}}}\sum_{i}\delta\ell_{i}.$		(25)

Thus, the spectrum of the radiation field $J_{\nu}$ can be obtained by adding the pathlengths of the packet’s trajectory between events. In the case of an infinite, plane-parallel slab, the luminosity should be replaced by $L_{*}=\mathcal{L}_{*}\Delta A$, where $\mathcal{L_{*}}$ is the photon energy per unit area per unit time, and $\Delta A$ is the area of the simulation box in the $XY$ plane. The volume element is $\Delta V=H\Delta A$ for the height $H$ of the slab. We also note that all intensity spectra $J_{x}\equiv J_{\nu}\Delta\nu_{{\rm D}}$ in this paper (including the emergent spectra) were angle-averaged and normalized for a unit production-rate of photons; thus, they should be multiplied by the Ly$\alpha$ production rate ($L_{*}/h\nu_{\alpha}$ or $\mathcal{L}_{*}/h\nu_{\alpha}$).

In Appendixes B and C, we derive the analytic formulae (Equations (B5) and (C3)) for the mean intensity at an arbitrary position in a static, infinite slab and sphere, respectively. The formulae are used to validate the Monte-Carlo simulation. Neufeld (1990) and Dijkstra et al. (2006) provided series solutions for the mean intensity in the slab and sphere, respectively. Our formulae are based on the series solutions, but we offer the solutions in more compact, analytic forms.

We also calculate the scattering rate $P_{\alpha}$, which is defined as the number of scatterings per unit time that a single atom undergoes. By the definition of the cross-section, the scattering rate can be calculated by integrating the scattering cross-section multiplied by the mean intensity:

$$P_{\alpha}=4\pi\int J_{\nu}\sigma_{\nu}d\nu.$$

(26)

We note that, at the central part of the Ly$\alpha$ line, the spectrum can be approximated to be constant as the cross-section peaks sharply. We then obtain

	$\displaystyle P_{\alpha}$	$\displaystyle\simeq$	$\displaystyle 4\pi J_{\nu}(\nu_{\alpha})\int\sigma_{\nu}d\nu=4\pi\chi_{0}J_{\nu}(\nu_{\alpha})$		(27)
		$\displaystyle=$	$\displaystyle 4\pi\frac{\chi_{0}}{\Delta\nu_{{\rm D}}}J_{x}(0)$
		$\displaystyle=$	$\displaystyle 1.315\times 10^{-12}\left(T_{{\rm K}}/10^{4}\,{\rm K}\right)^{-1/2}J_{x}(0)$

In Section 4, we show that the radiation field spectrum is indeed constant if the atomic recoil effect is ignored or exponential if the recoil effect is considered, over a wide range of frequency. The exponential function can be approximated to a linear function at the central portion of the spectrum and thus the integral in Equation (27) should give rise to precisely the same result as for a constant spectrum. In Appendixes B and C, we provide useful formulae (Equations (B) and (C)) for the mean intensity at the line center $J_{x}(0)$ in an any position of a slab and sphere to calculate the scattering rate $P_{\alpha}$.

In the Monte-Carlo RT, it is more natural to calculate the scattering rate by directly counting the number of scattering events in each volume element and normalize it by the number of atoms. The total number of photons injected during $\Delta t$ is $L_{*}\Delta t/h\nu_{\alpha}$. The scattering rate can then be obtained by

$\displaystyle P_{\alpha}$

$\displaystyle=$

$\displaystyle\frac{1}{n_{{\rm H}}\Delta V}\frac{L_{*}/h\nu_{\alpha}}{N_{{\rm packet}}}N_{{\rm scatt}}^{{\rm cell}},$

(28)

where $N_{{\rm scatt}}^{{\rm cell}}$ is the number of scattering events that occur in the volume element $\Delta V$. Again, for the infinite slab geometry, the luminosity and volume element should be replaced by $L_{*}=\mathcal{L}_{*}\Delta A$ and $\Delta V=H\Delta A$, respectively.

The third estimator of $P_{\alpha}$ can be obtained by combing Equations (25) and (26), as follows:

$$P_{\alpha}=\frac{1}{n_{{\rm H}}\Delta V}\frac{L_{*}/h\nu_{\alpha}}{N_{{\rm packet}}}\sum_{i,\,\nu}\delta\tau_{i,\nu},$$

(29)

where $\delta\tau_{i,\nu}=n_{{\rm H}}\sigma_{\nu}\delta l_{i}$ is the optical depth segments between events. Here, the optical depths include only those estimated for hydrogen atoms, but not due to dust extinction.

Directly counting the number of scattering events is much simpler than measuring the radiation field spectrum at the line center. However, the method of using the radiation field spectrum is advantageous to derive approximate formulae for $P_{\alpha}$ from analytic solutions of the RT equation and to check the self-consistency of the Monte-Carlo simulation. We also note that Equation (29) is the same as Equation (28), except that the number of scattering events is replaced with the summation of optical depths. The probability of scattering events in an optical depth interval of $\delta\tau$ is $\mathcal{P}=1-e^{-\delta\tau}\simeq\delta\tau$; therefore, Equations (28) and (29) are in fact equivalent. An advantage of using Equation (29) is that the estimator gives a smooth, non-zero $P_{\alpha}$ even in the cells in which the density of neutral hydrogen is very low and thus no scattering occurs ($N_{{\rm scatt}}^{{\rm cell}}=0$), in a practical sense, because of a limited number of photon packets. On the other hand, the estimator in Equation (28) yielded $N_{{\rm scatt}}^{{\rm cell}}=0$ in more than 50% of cells, corresponding to hot and fully ionized gas, when it was applied to the simulation snapshots in Section 5.4. The results in Section 5 were thus obtained using Equation (29). However, the above three methods of calculating $P_{\alpha}$ gave no significant differences in our simulations.

The average number of scatterings per a photon before escape can be obtained by integrating the scattering rate for a single photon over the whole volume, after multiplying by the neutral hydrogen density:

$$\left\langle N_{{\rm scatt}}\right\rangle=\int P_{\alpha}n_{{\rm H}}dV.$$

(30)

Using this equation and Equation (27), we can also derive the analytic formulae for the mean number of scatterings in the cases of a static slab and sphere, as described in Appendixes B and C, respectively. The resulting formulae for the mean number of scatterings are

	$\displaystyle\left\langle N_{{\rm scatt}}^{{\rm slab}}\right\rangle$	$\displaystyle=$	$\displaystyle 1.612\tau_{0},$		(31)
	$\displaystyle\left\langle N_{{\rm scatt}}^{{\rm sphere}}\right\rangle$	$\displaystyle=$	$\displaystyle 0.9579\tau_{0}$		(32)

for a slab and sphere, respectively. The equation for a slab is the same as the first derived by Harrington (1973); the equation for a sphere is newly-derived in Appendix C.

As a test of our code, we performed the RT calculation for the cases of a static, plane-parallel slab, and sphere, for which Neufeld (1990) and Dijkstra et al. (2006), respectively, derived analytic approximations of the emergent mean intensity at the limit of large optical thickness. We also present new analytic formulae, which are slightly different from theirs, in Appendixes B and C. For these tests, we ignored the recoil effect unless otherwise stated; in other words, the recoil parameter $g_{*}$ in Equation (15) was assumed to be zero.

We first consider the infinite slab case in which the Ly$\alpha$ source is located at the mid-plane and injects monochromatic photons at frequency $x=0$ ($\nu=\nu_{\alpha}$). Figure 1 shows our simulation results for the slab with a temperature of (a) $T_{{\rm K}}=10$ K and (b) $T_{{\rm K}}=10^{4}$ K together with the approximate, analytic solution of Neufeld (1990). The optical thickness at line center ($\tau_{0}$) varies from $10^{5}$ to $10^{8}$, as indicated in the figure. Our results agree with the analytic solution, in general. However, as the optical thickness decreases and the temperature increases, the analytic solution begins to deviate from the simulation results. This discrepancy is because the analytic solution was derived in the limit of optically thick media ($2a\tau_{0}=9.4(T_{{\rm K}}/10^{4}{\rm K})^{-1/2}(\tau_{0}/10^{4})\gtrsim 10^{3}$).

As the second test, Figure 2 shows the simulation results for the case of a sphere in which photons are injected at frequency $x=0$ from the center of the sphere. The temperature of the medium is assumed to be (a) $T_{{\rm K}}=10$ K and (b) $T_{{\rm K}}=10^{4}$ K, and the optical depth at line center varies from $10^{5}$ to $10^{8}$, as for the slab geometry. In the figure, Equation (C3) given in Appendix C is compared. The analytic curves are in good agreement with the simulation results, except for $T_{{\rm K}}=10^{4}$ K and $\tau_{0}\lesssim 10^{6}$.⁴⁴4The analytic solutions (Equations (B5) and (C3)) for slab and sphere geometries were derived by assuming the line wing approximation for the Voigt function, i.e., $H(x,a)\approx a/(\pi^{1/2}x^{2})$. They, thus, do not well reproduce the broad and deep U-shape feature at $x\approx 0$, which is caused by the core scattering, as shown for $T_{{\rm K}}=10^{4}$ K and $\tau_{0}=10^{5}$ in Figure 1. If we use the exact Voigt function in the equations, the U-shape is well reproduced. However, in that case, the resulting line profile gives rise to an underestimation of the total flux, for the models with $\tau_{0}\lesssim 10^{6}$, integrated over frequency.

To test the code for the dynamic case, we examine the Hubble-like flow model. In the model, we consider an isothermal, homogeneous sphere, which is isotropically expanding or contracting. The bulk velocity of a fluid element at a distance $r$ from the center is assumed to be

$$\textbf{v}_{{\rm bulk}}(\mathbf{r})=\left(\frac{v_{{\rm max}}}{R}\right)\mathbf{r},$$

(33)

where $R$ is the radius of the medium and $v_{{\rm max}}$ is the maximum velocity at the edge of the sphere ($r=R$). There is no analytical solution of the emergent spectrum for a moving medium with non-zero temperature.⁵⁵5The analytic solution presented in Loeb & Rybicki (1999) is for the zero temperature. In other words, no thermal broadening was taken into account. We thus use the same parameters as the models of Laursen et al. (2009), Yajima et al. (2012), and Smith et al. (2017) and compare our simulation results with theirs. In the models, the gas medium is assumed to have a temperature of $T_{{\rm K}}=10^{4}$ K and a column density of $N_{{\rm HI}}=2\times 10^{20}$ cm^-2 (corresponding to $\tau_{0}\simeq 1.2\times 10^{7}$). The sphere of gas expands isotropically with $v_{{\rm max}}=20$, 200, or 2000 km s^-1. Figure 3 shows excellent agreement of our results with those of Laursen et al. (2009), Yajima et al. (2012), and Smith et al. (2017). We note that, as $v_{{\rm max}}$ increases, the blue peak is suppressed and disappears entirely at $v_{{\rm max}}\sim 200$ km s^-1. The peak location is also pushed further from the line center. These trends are because photons are Doppler shifted out of the line center due to the velocity gradient. However, at the extreme velocity of $v_{{\rm max}}=2000$ km s^-1, the peak approaches the line center again because the extreme velocity gradient allows photons to escape more easily. If a collapsing sphere ($v_{{\rm max}}<0$) is considered, then the red wing is suppressed, and the blue wing is enhanced, as opposed to the outflowing case.

The average number of scatterings of Ly$\alpha$ experienced before escaping from the gas cloud is also an interesting quantity to examine. Figures 4(a) and (b) show the average number of scatterings for the infinite slab and sphere, respectively, as a function of the central optical thickness $\tau_{0}$. We calculated the number of scatterings in the cases of the gas temperature $T_{{\rm K}}=10$ K (blue circles) and $T_{{\rm K}}=10^{4}$ K (red dots) for a wide range of optical thickness $\tau_{0}=1-10^{9}$. The results cover the widest range reported to date. Overplotted as black dashed lines are the theoretical prediction (a) for the slab geometry (Equation (31)) and (b) for the spherical geometry (Equation (32)).

It is evident that the mean number of scatterings in the slab is higher than that in the sphere, because the slab geometry is infinitely extended along the XY plane and thus photons can escape only in the direction perpendicular to the XY plane.We also note that, in the optical thickness range of $10^{2}\lesssim\tau_{0}\lesssim 10^{6}$, Ly$\alpha$ photons in the lower-temperature medium undergo less scattering than in the higher-temperature medium. This is not caused by poor statistics; we increased the number of photon packets from $10^{5}$ up to $10^{9}$ and found no significant differences in the results.

Note that the mean numbers of scatterings for $\tau_{0}=10^{4}$ and $10^{5}$ in the bottom left panel of Figure 1 in Dijkstra et al. (2006) are lower than our results as well as the analytic approximation. This may be due to Dijkstra et al. (2006) having used the weighting scheme developed by Avery & House (1968), as described in Appendix B1 of Dijkstra et al. (2006). The distribution of the number of scatterings has a long tail to a large number of scatterings. The weighting scheme gives less weight to the photons that experience a large number of scatterings and thus underestimates the average number of scatterings.

So far, dust is neglected in the tests. The presence of dust grains in the medium suppresses the Ly$\alpha$ intensity escaping from the medium. Figure 5 shows the emergent spectra ($J_{{\rm esc}}$) from a dusty, homogeneous sphere of which the gas-to-dust ratio is that of the Milky Way ($\sim 100$ by mass). The spectra were obtained for four combinations of the gas temperature ($T_{{\rm K}}=10$ and $10^{4}$ K) and the optical thickness ($\tau_{0}=10^{5}$ and $10^{6}$). As shown in the figure, the overall shape of the emergent Ly$\alpha$ spectrum is not significantly altered by dust. In the figure, we also show the spectra of photons that are absorbed by dust ($J_{{\rm abs}})$. The absorbed spectra are not observable but may be useful when calculating the heating of dust grains in the vicinity of H II regions by Ly$\alpha$ photons.

It is worthwhile to note that there is a spiky dip feature at the line center of the dust-absorption spectra for the gas of $T_{{\rm K}}=10$ K. The feature indicates that Ly$\alpha$ photons are less absorbed at the very line core than at the outside of the core. Verhamme et al. (2006) also found a similar dip feature, as shown in Figure 3 in their paper, but it is much broader and deeper than ours. We found that the dip is noticeable only when the optical depth and temperature are rather low. As the optical depth and temperature increase, the dip feature is found to be smeared out and eventually disappear in the absorbed spectra. The dip feature is discussed in more detail in Appendix 32.

We also estimated the escape fraction of Ly$\alpha$ photons from a dusty, infinite slab. Neufeld (1990) derived an analytic equation for the escape fraction at a limit of large optical depths, as follows:

$$f_{{\rm esc}}=\left.1\right/\cosh\left[\frac{\sqrt{3}}{\pi^{5/2}\zeta}\left\{\left(a\tau_{0}\right)^{1/3}\tau_{{\rm abs}}\right\}^{1/2}\right].$$

(34)

He estimated $\zeta\approx 0.525$ by comparing this equation with numerical results of Hummer & Kunasz (1980), obtained for $\tau_{0}=5\times 10^{7}$ and $a=4.7\times 10^{-2}$ (or equivalently $T_{{\rm K}}=1$ K). Figure 6 shows the escape fraction as a function of $(a\tau_{0})^{1/3}\tau_{{\rm abs}}$ calculated for $\tau_{0}=10^{4}-10^{8}$ and $T_{{\rm K}}=10$, $10^{4}$ K. The dust optical thickness was varied so that $(a\tau_{0})^{1/3}\tau_{{\rm abs}}$ ranges from $10^{-2}$ to $10^{2}$. To compare our results with those of others, we assumed the dust scattering albedo of 0.5 and isotropic scattering ($g=0$). It appears that Equation (34) with $\zeta=0.525$ is indeed a good approximation for the cases of $\tau_{0}>10^{7}$. However, $\zeta=0.525$ overpredicts the escape fraction when the optical depth is lower than $10^{7}$. From the figure, it is also clear that $\zeta\approx 0.525$ is more adequate for lower temperature. This is because the equation was derived under the condition of $(a\tau_{0})^{1/3}>\tau_{{\rm abs}}$ and the width ($a$) of the Voigt profile is proportional to $T_{{\rm K}}^{-1/2}$. However, adopting a different value of $\zeta$ for different temperature $T$ and optical depth $\tau_{0}$, we found that the equation still provides a reasonably good representation of the escape fraction even when the condition required for the equation is not satisfied. The $\zeta$ values suitable for different temperatures $T_{{\rm K}}$ and optical depths $\tau_{0}$ are shown in Figure 6.

We present the scattering rate obtained for an unit generation rate of Ly$\alpha$ photons. Figure 7 compares the scattering rate calculated through a Monte-Carlo RT simulation and that obtained using an analytic approximation for the slab geometry. The results for the sphere geometry are shown in Figure 8. In the figures, the solid black lines represent the scattering rate calculated by directly counting the scattering events in the simulation. The red dots were estimated using the Ly$\alpha$ spectrum at the line center $J_{x}(x=0)$ and Equation (27). The green dashed lines denote the analytical, approximate solutions at the limit of a large optical thickness. The approximate formulae for the scattering rate in slab and sphere geometries are, respectively, obtained by combining Equations (B) and (C) for $J_{x}(0)$ with Equation (27). Note that the analytic solutions overlap with the simulation results for $\tau_{0}=10^{7}$ and are barely distinguishable in the figures.

The functional shape of $P_{\alpha}$ as a function of $\tau/\tau_{0}$ is more or less independent of the gas temperature and total optical thickness $\tau_{0}$. However, its absolute strength is proportional to $T_{{\rm K}}^{-1/2}$, as shown in Equation (27). Therefore, in Figures 7 and 8, the scattering rate $P_{\alpha}$ at $T_{{\rm K}}=10$ K is greater than that at $T_{{\rm K}}=10^{4}$ K by a factor of $\sim 31$.6 for a given optical depth $\tau_{0}$. We also note that the functional form of $P_{\alpha}$ in a sphere is much steeper than that in a slab. This tendency is due to photons in an infinite slab being able to escape the medium only through one direction that is perpendicular to the slab.

As $\tau_{0}$ increases, the scattering rate $P_{\alpha}(\tau/\tau_{0})$ gradually decreases and converges to an asymptotic form, which is well represented by the analytic formulae denoted by the green dashed lines in Figures 7 and 8. The convergence is achieved at a lower optical thickness when the medium has a lower temperature. As denoted in the figures, good agreement between the formula and simulation was obtained at $\tau_{0}\gtrsim 10^{5}$ for $T_{{\rm K}}=10$ K, but at a slightly higher optical depth of $\tau_{0}\gtrsim 10^{6}$ for $T_{{\rm K}}=10^{4}$ K. A similar trend can be also found in Figure 4, where the average number of scatterings begins to approach the asymptotic solution at a lower optical depth when the medium is in a lower temperature ($T_{{\rm K}}=10$ K) than in a higher temperature ($T_{{\rm K}}=10^{4}$ K).

The scattering rate $P_{\alpha}$ is affected by the presence of dust grains in the medium as well as by the systematic, bulk motion of the gas. Figure 9 shows the dependence of $P_{\alpha}$ on the dust-absorption optical depth in a slab. In the figure, the gas temperature and optical depth are $T_{{\rm K}}=10^{4}$ K and $\tau_{0}=10^{7}$, respectively, and the dust-absorption optical depth varies from $\tau_{{\rm abs}}=0$.0 (no dust) to $\tau_{{\rm abs}}=0.24$. These results were adopted from the models calculated for Figure 6. As expected, when $\tau_{{\rm abs}}$ gradually increases, $P_{\alpha}$ begins to drop. We note that, if the gas-to-dust ratio of the Milky Way is assumed, the absorption optical depth of the medium with $T_{{\rm K}}=10^{4}$ K and $\tau_{0}=10^{7}$ is $\tau_{{\rm abs}}=0.14$. This value is similar to the second-highest $\tau_{{\rm abs}}$ in Figure 9, indicating that dust grains effectively absorb Ly$\alpha$ photons and reduce the number of Ly$\alpha$ scatterings at a location of $\tau/\tau_{0}=0.5$ by a factor of $\sim 10$.

Figure 10 shows the effect of outflow on the scattering rate. In the figure, we assumed a Hubble-like, expanding sphere with a temperature of $T_{{\rm K}}=100$ K and an optical depth of $\tau_{0}=10^{3}$. The maximum velocity is increased from $v_{{\rm max}}=0$ to 10 km s^-1. As expected, the scattering rate rapidly decreases as the expansion velocity rises. Figures 9 and 10 illustrate how dramatically the number of scatterings can be reduced by the presence of dust and the bulk motion of the gas.

We first examine the line profile of Ly$\alpha$ formed within a spherical medium when the recoil effect is ignored, i.e., $g_{*}=0$ in Equation (15). Figures 11 and 12 show the line profiles at various radial locations of a spherical medium with a temperature of $T_{{\rm K}}=10$ and $10^{4}$ K, respectively. In each figure, the results for two different optical depths ($\tau_{0}=10^{5}$ and $10^{7}$ in Figure 11, and $\tau_{0}=10^{6}$ and $10^{7}$ in Figure 12) are shown. The black lines in the figures denote the simulation results, while the red lines show the analytic solutions given in Equation (C3). It is clear that the approximate formula reproduces the simulation results reasonably well, in particular, when the optical depth is large enough, i.e., $\tau_{0}\gtrsim 10^{7}$ and the gas temperature is low ($T_{{\rm K}}=10$ K). For $T_{{\rm K}}=10^{4}$ K and $\tau_{0}<10^{6}$, however, the analytic solution significantly underestimates the intensity. We also note that the analytic solution tends to underestimate the Ly$\alpha$ intensity in the central portion and overestimate in outer parts. The same trend was found by Higgins & Meiksin (2012) who solved the moment equations numerically. We also found a similar trend for a slab, but not shown in this paper. It should also be mentioned that the Ly$\alpha$ spectrum evolves progressively to a typical double-peak shape as upon approaching the system boundary ($\tau/\tau_{0}\gtrsim 0.95$), but not clearly shown in the figures.

We next consider the case in which the atomic recoil effect is taken into account in the simulation. Figure 13 shows the Ly$\alpha$ spectra in a static, infinite slab for six different sets of physical parameters. The gas temperature is assumed to be $T_{{\rm K}}=10$ K in Figures 13(a) and (b), while it is $T_{{\rm K}}=10^{2}$ K in Figures 13(c) and (d), and $T_{{\rm K}}=10^{4}$ K in Figures 13(e) and (f). The optical thickness for the system is $\tau_{0}=10^{2}$ for Figures 13(a), (c), and (e), and $\tau_{0}=5\times 10^{2}$ for Figures 13(b), (d), and (f). In each figure, the left panels show the Ly$\alpha$ spectra measured in an optical depth bin of $0.2<\tau/\tau_{0}<0.5$; the right panels show the spectra in $0.5<\tau/\tau_{0}<0.8$.

As noted in Section 2.3, the frequency gap between the fine-structure levels of the $n=2$ state can be comparable to the line width formed in the medium with a relatively low optical depth and low temperature (i.e., $T_{{\rm K}}\lesssim 10^{2}$ K). In Figures 13(a) to (d), we, therefore, show also the Ly$\alpha$ spectra obtained when the fine structure is considered. In the figures, the black and blue solid lines denote the Ly$\alpha$ spectra calculated when the fine structure is ignored and taken into consideration, respectively. In Figure 13(a) for $\tau_{0}=10^{2}$, we see that the spectra, denoted in blue, show two steps due to the fine structure. However, as the optical thickness increases higher than $\tau_{0}\approx 5\times 10^{2}$, the step-like shape disappears, as shown in Figure 13(b) for $\tau_{0}=5\times 10^{2}$. The intensity ratio of the doublet line becomes 1:1 at high optical depths as a result of a considerable number of resonance scatterings, and thus the step signature of the doublet disappears. The step shape disappeared at a lower optical depth as temperature increases. For instance, in the case of $T_{{\rm K}}=10^{2}$ K (Figure 13(c)), the step shape is not found already at $\tau_{0}=10^{2}$.

The most notable thing in Figures 13(a) to (d) is that the Ly$\alpha$ spectra are tilted in the central part, unlike those shown in Figures 11 and 12 in which no recoil effect was considered. To compare with the expectation in the WF effect theory, we overlay an exponential function $\exp\left[-h(\nu-\nu_{\alpha})/k_{{\rm B}}T_{{\rm K}}\right]$ for the gas temperature $T_{{\rm K}}$ with a red dashed line. As shown in the figures, the spectra clearly appears to be consistent with the function expected for the WF effect. However, in Figures 13(e) and (f) for the highest temperature of $T_{{\rm K}}=10^{4}$ K, the spectral shape appears to be nearly constant at the central part. In fact, the exponential function for a high $T_{{\rm K}}$ is almost constant. Hence, at first glance, it seems that it is not possible to confirm whether the Ly$\alpha$ spectra within the medium at a high temperature follow an expected functional form or not.

However, we note that the spectrum can be expressed by

$$\frac{k_{{\rm B}}T_{{\rm K}}}{h\Delta\nu_{{\rm D}}}\ln\left(J_{x}/J_{0}\right)=-x.$$

(35)

Thus, if the spectrum is drawn in logarithmic scale, as in the left side of the above equation, and compared with $-x$, we can readily illustrate that the spectrum follows the expected function. Figure 14 compares the Ly$\alpha$ spectra represented in logarithmic scale with $-x$, which is indicated by the dashed red line, for all cases shown in Figure 13. Consequently, we conclude that the Ly$\alpha$ spectrum is well described by an exponential function with a slope of the gas temperature, even at the highest temperature ($T_{{\rm K}}=10^{4}$ K) presented here. Interestingly, Figures 13(a) and 14(a) show that the double steps (blue lines) have the same exponential slope corresponding to the gas temperature, although they are disconnected.

Here, we note that the exponential shape should remain valid across the frequency interval corresponding to the width of the resonance profile $\phi_{x}$, i.e., $|x|\lesssim x_{*}=0.84\ (=\ln 2)$. In this paper, we take this width as the full-width at half maximum (FWHM) of the Voigt function for $a\ll 1$, or equivalently the FWHM of the Gaussian function with a standard deviation of $1/\sqrt{2}$. If we additionally consider the line splitting due to the fine structure, the exponential shape should be kept over the whole range of the frequency gap of the doublet plus the width of the resonance profile, i.e., $|x|\lesssim x_{*}=0.84+|\nu_{\nicefrac{{3}}{{2}}}-\nu_{\nicefrac{{1}}{{2}}}|/(2\Delta\nu_{{\rm D}})$. Otherwise, the 21-cm spin temperature could not approach the gas temperature. Note that $|\nu_{\nicefrac{{3}}{{2}}}-\nu_{\nicefrac{{1}}{{2}}}|/(2\Delta\nu_{{\rm D}})=5.83\times 10^{-2}\left(T_{{\rm K}}/10^{4}\,\text{K}\right)^{-1/2}$. In order to check if this is the case, the vertical blue dashed lines are overlaid in Figures 13 and 14 (and in other figures as well where necessary) to denote the frequency interval $|x|\leq x_{*}$.

In the case of $T_{{\rm K}}=10$ K and $\tau_{0}=10^{2}$ shown in Figure 14(a), it appears that the spectral shape does not satisfy the requirement because of the discontinuous step feature. On the other hand, we found that the spectra obtained for $\tau_{0}\gtrsim 3\times 10^{2}$ and $T_{{\rm K}}=10$ K follows the required exponential shape over the frequency interval $|x|\leq x_{*}$, for instance, as shown in Figure 14(b). In the case of $T_{{\rm K}}=10^{2}$ K ($\tau_{0}=10^{2}$ and $\tau_{0}=5\times 10^{2}$), the desired profile is more clearly seen in Figures 14(c) and (d). Figures 14(e) and (f) also show that the line profiles follow the expected shape very well, over the frequency interval $|x|\leq x_{*}$, even for the highest temperature $T_{{\rm K}}=10^{4}$ K. We examined many different cases and found that the necessary condition for the Ly$\alpha$ spectral profile is fulfilled at an optical depth as low as $\tau_{0}\simeq(1-5)\times 10^{2}$. In general, as the temperature is low, it appeared that the requirement is satisfied at a lower optical depth. If the optical depth is lower than $\tau_{0}\simeq 10^{2}$, the Ly$\alpha$ spectrum was, in general, not fully thermalized in the temperature range of $10\,\text{K}\leq T_{{\rm K}}\leq 10^{4}\,\text{K}$, examined in this paper.

The above results were obtained in static media. It would, therefore, be worthwhile to examine whether the Ly$\alpha$ spectrum will have a shape of the required exponential function in a medium in motion. For instance, we show the spectra within a Hubble-like, expanding medium in Figures 15 and 16. Both figures were obtained for an expanding medium with $\tau_{0}=10^{3}$, $T_{{\rm K}}=10^{2}$ K. Figure 15 shows the Ly$\alpha$ spectra at several locations, expressed in terms of optical depth, in a spherical medium with a maximum expanding velocity of $v_{{\rm max}}=5$ km s^-1. In Figure 15, the spectra near the line center ($x=0$) seem flat. However, more detailed views in Figure 16 clearly illustrate that the spectra calculated at various optical depth bins follow the exponential function for $T_{{\rm K}}=10^{2}$ K. We, therefore, can conclude that the Ly$\alpha$ line profile follows the desired exponential function with a slope corresponding to the gas temperature even in an expanding medium, if the expanding velocity is not too fast. If a medium rapidly expands or, more precisely, has a velocity gradient, across a region with an optical thickness of $\tau_{0}\approx 100-500$, larger than the thermal velocity dispersion, Ly$\alpha$ photons will escape the region before they fully develop the spectral shape required for the WF effect. The relevant issue in moving media is further addressed in the next section.

The initial frequency of Ly$\alpha$ was fixed at $\nu=\nu_{\alpha}$ ($x=0$) for the above discussion. The same conclusion was obtained even when we assumed a Voigt profile or a flat continuum. The Ly$\alpha$ line is broadened or adjusted very quickly in a medium with an optical thickness of $\tau_{0}\gtrsim 10^{2}$ so that the initial frequency distribution did not significantly alter the resulting spectra; only the spectra in the vicinity of the source showed a weak dependence on the initial frequency distribution.

In the above calculation of the Ly$\alpha$ RT, only the thermal motion and systematic bulk flow of gas were taken into account. Aside from the thermal broadening, the line width of Ly$\alpha$ may also be broadened by the presence of the turbulent velocity field. It is convenient to define the Doppler temperature $T_{{\rm D}}$ as in Liszt (2001), to incorporate the turbulence effect (see also Draine, 2011):

$$2k_{{\rm B}}T_{{\rm D}}/m_{{\rm H}}=2k_{{\rm B}}T_{{\rm K}}/m_{{\rm H}}+v_{{\rm turb}}^{2},$$

(36)

where $v_{{\rm turb}}$ is the velocity dispersion of turbulent motion. The Doppler parameter $b$ defined by

$$b=\sqrt{v_{{\rm th}}^{2}+v_{{\rm turb}}^{2}}$$

(37)

is also useful to describe effects caused by the turbulent motions, as in Verhamme et al. (2006). The Doppler temperature $T_{{\rm D}}$ (or Doppler parameter $b$) has been used in place of the thermal kinetic temperature (or thermal velocity) to describe the Ly$\alpha$ spectral shape. However, it has not been verified whether this prescription of incorporating the turbulent fluid motion is appropriate in the context of Ly$\alpha$ RT, especially with regards to the WF effect. Hydrodynamic simulations for galaxies, where the turbulence effects are included, have been adopted to calculate the emergent Ly$\alpha$ spectra. However, in these studies, the turbulence effect could not be separated from the pure thermal motion effect.