# Effect of Dust on Lyman-alpha Photon Transfer in Optically Thick Halo

Yang Yang11affiliation: Division of Applied Mathematics, Brown University,
Providence, RI 02912, USA , Ishani Roy22affiliation: Computing Laboratory,
University of Oxford, Oxford, OX1 3QD, United Kingdom , Chi-Wang
Shu11affiliation: Division of Applied Mathematics, Brown University,
Providence, RI 02912, USA and Li-Zhi Fang33affiliation: Department of Physics,
University of Arizona, Tucson, AZ 85721, USA

###### Abstract

We investigate the effects of dust on Ly$\alpha$ photons emergent from an
optically thick medium by solving the integro-differential equation of the
radiative transfer of resonant photons. To solve the differential equations
numerically we use the Weighted Essentially Non-oscillatory method (WENO).
Although the effects of dust on radiative transfer is well known, the resonant
scattering of Ly$\alpha$ photons makes the problem non-trivial. For instance,
if the medium has the optical depth of dust absorption and scattering to be
$\tau_{a}\gg 1$, $\tau\gg 1$, and $\tau\gg\tau_{a}$, the effective absorption
optical depth in a random walk scenario would be equal to
$\sqrt{\tau_{a}(\tau_{a}+\tau)}$. We show, however, that for a resonant
scattering at frequency $\nu_{0}$, the effective absorption optical depth
would be even larger than $\tau(\nu_{0})$. If the cross section of dust
scattering and absorption is frequency-independent, the double-peaked
structure of the frequency profile given by the resonant scattering is
basically dust-independent. That is, dust causes neither narrowing nor
widening of the width of the double peaked profile. One more result is that
the time scales of the Ly$\alpha$ photon transfer in the optically thick halo
are also basically independent of the dust scattering, even when the
scattering is anisotropic. This is because those time scales are mainly
determined by the transfer in the frequency space, while dust scattering,
either isotropic or anisotropic, does not affect the behavior of the transfer
in the frequency space when the cross section of scattering is wavelength-
independent. This result does not support the speculation that dust will lead
to the smoothing of the brightness distribution of Ly$\alpha$ photon source
with optical thick halo.

cosmology: theory - intergalactic medium - radiation transfer - scattering

##  1 Introduction

Ly$\alpha$ photons have been widely applied to study the physics of luminous
objects at various epochs of the universe, such as Ly$\alpha$ emitters,
Ly$\alpha$ blob, damped Ly$\alpha$ system, Ly$\alpha$ forest, fluorescent
Ly$\alpha$ emission, star-forming galaxies, quasars at high redshifts as well
as optical afterglow of gamma ray bursts (Haiman et al. 2000; Fardal et al.
2001; Dijkstra & Loeb. 2009; Latif et al. 2011). The resonant scattering of
Ly$\alpha$ photons with neutral hydrogen atoms has a profound effect on the
time, space and frequency dependencies of Ly$\alpha$ photons transfer in an
optically thick medium. Ly$\alpha$ photons emergent from an optically thick
medium would carry rich information of photon sources and halo surrounding the
source of the Ly$\alpha$ photon. The profiles of the emission and absorption
of the Ly$\alpha$ radiation are powerful tools to constrain the mass density,
velocity, temperature and the fraction of neutral hydrogen of the optically
thick medium. Radiation transfer of Ly$\alpha$ photons in an optically thick
medium is fundamentally important.

The radiative transfer of Ly$\alpha$ photons in a medium consisting of neutral
hydrogen atoms has been extensively studied either analytically or
numerically. Yet, there have been relatively few results which are directly
based on the solutions of the integro-differential equation of the resonant
radiative transfer. Besides the Field solution (Field 1959, Rybicki &
Dell’Antonio 1994), analytical solutions with and without dust mostly are
based on the Fokker-Planck (P-F) approximation (Harrington 1973, Neufeld 1990,
Dijkstra et al. 2006). The P-F equation might miss the detailed balance
relationship of resonant scattering (Rybicki 2006), and therefore, the
analytical solutions cannot describe the formation and evolution of the
Wouthuysen-Field (W-F) local thermalization of the Ly$\alpha$ photon frequency
distribution (Wouthuysen 1952, Field 1958), which is important for the
emission and absorption of the hydrogen 21 cm line (e.g. Fang 2009). The
features of the Ly$\alpha$ photon transfer related to the W-F local
thermalization are also missed. An early effort (Adams et al. 1971) trying to
directly solve the integro-differential equation of the resonant radiative
transfer with numerical method. It still is, however, of a time-independent
approximation.

Recently, a state-of-the-art numerical method has been introduced to solve the
integro-differential equation of the radiative transfer with resonant
scattering (Qiu et al. 2006, 2007, 2008, Roy et al. 2009a). The solver is
based on the weighted essentially non-oscillatory (WENO) scheme (Jiang & Shu
1996). With the WENO solver, many physical features of the transfer of
Ly$\alpha$ photons in an optically thick medium (Roy et al. 2009b, 2009c,
2010), which are missed in the Fokker-Planck equation approximations, have
been revealed. For instance, the WENO solution shows that the time scale of
the formation of the W-F local thermal equilibrium actually is only about a
few hundred times of the resonant scattering. It also shows that the double
peaked frequency profile of the Ly$\alpha$ photon emergent from an optically
thick medium does not follow the time-independent solutions of the P-F
equation. These results directly indicate the needs of re-visiting problems
which have been studied only via the F-P time-independent approximation.

We will investigate, in this paper, the effects of the dust on the Ly$\alpha$
photons transfer in an optically thick medium. Dust can be produced at epochs
of low and moderate redshifts, and even at redshift as high as 6 (Stratta et
al. 2007). Absorption and scattering of dust have been used to explain the
observations on Ly$\alpha$ emission and absorption (Hummer & Kunasz 1980),
such as the escaping fraction of Ly$\alpha$ photons (Hayes et al. 2010, 2011,
Blanc et al. 2010); the redshift-dependence of the ratio between Ly$\alpha$
emitters and Lyman Break galaxies (Verhamme et al. 2008); and the “evolution”
of the double-peaked profile (Laursen et al. 2009).

Nevertheless, it is still unclear whether the time scale of photon escaping
from optically thick halo will be increasing (or decreasing) when the halo is
dusty. It is also unclear whether the effects of dust absorption can be
estimated by the random walk picture (Hansen & Oh 2006). As for the dust
effect on the double-peaked profile, the current results given by different
studies seem to be contradictory: some claims that the dust absorption leads
to the narrowing of the double-peaked profile (Lauresen et al 2009), while
others result that the width between the two peaks apparently should be
increasing due to the dust absorption (Verhamme et al. 2006). We will focus on
these basic problems, and examine them with the solution of the integro-
differential equation of radiative transfer.

This paper is organized in the following way: section 2 presents the theory of
the Ly$\alpha$ photon transfer in an optically thick medium with dust. The
equations of the intensity and flux of resonant photons in a dusty medium are
given. We will study three models of the interaction between dust and photons:
(1) dust causes only scattering with photons; (2) dust causes both scattering
and absorption; and (3) dust causes only absorption of photons. Section 3
gives the solutions of Ly$\alpha$ photons escaping from an optically thick
spherical halos with dust. The dusty effect on the double-peaked profile will
be studied in Section 4. The discussion and conclusion are given in Section 5.
Some mathematical derivations of the equations and numerical implementation
details are given in the Appendix.

##  2 Basic theory

###  2.1 Radiative transfer equation of dusty halo

We study the transfer of Ly$\alpha$ photons in a spherical halo with radius
$R$ around an optical source. The halo is assumed to consist of uniformly
distributed HI gas and dust. The optical depth of HI scattering over a light
path $dl$ is $d\tau=\sigma(\nu)n_{\rm HI}dl$, where $n_{\rm HI}$ is the number
density of HI, and $\sigma(\nu)$ is the cross section of the resonant
scattering of Ly$\alpha$ photons by neutral hydrogen, which is given by

$\sigma(x)=\sigma_{0}\phi(x,a)=\sigma_{0}\frac{a}{\pi^{3/2}}\int^{\infty}_{-\infty}dy\frac{e^{-y^{2}}}{(x-y)^{2}+a^{2}}$
(1)

where $\phi(x,a)$ is the normalized Voigt profile (Hummer 1965). As usual, the
photon frequency $\nu$ in eq.(1) is described by the dimensionless frequency
$x\equiv(\nu-\nu_{0})/\Delta\nu_{D}$, with $\nu_{0}=2.46\times 10^{15}$ s-1
being the resonant frequency, $\Delta\nu_{D}=\nu_{0}(v_{T}/c)=1.06\times
10^{11}(T/10^{4})^{1/2}$ Hz the Doppler broadening, $v_{T}=\sqrt{2k_{B}T/m}$
the thermal velocity, and $T$ the gas temperature of the halo.
$\sigma_{0}/\pi^{1/2}$ is the cross section of scattering at the the resonant
frequency $\nu_{0}$. The parameter $a$ in eq.(1) is the ratio of the natural
to the Doppler broadening. For the Ly$\alpha$ line, $a=4.7\times
10^{-4}(T/10^{4})^{-1/2}$. The optical depth of Ly$\alpha$ photons with
respect to HI resonant scattering is $\tau_{s}(x)=n_{\rm
HI}R\sigma(x)=\tau_{0}\phi(x,a)$, where $\tau_{0}=n_{\rm HI}\sigma_{0}R$.

If the absorption and scattering of dust are described by effective cross-
section per hydrogen atom $\sigma_{d}(x)$, the total optical depth is given by

$\tau(x)=\tau_{0}\phi(x,a)+\tau_{d}(x)$ (2)

where the dust optical depth $\tau_{d}(x)=n_{\rm HI}\sigma_{d}(x)R$. This is
equal to assume that dust is uniformly distributed in IGM. The effects of
inhomogeneous density distributions of dust (Neufeld 1991; Haiman & Spaans
1999) will not be studied in this paper.

The radiative transfer equation of Ly$\alpha$ photons in a spherical halo with
dust is given by

$\displaystyle{\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}=$
$\displaystyle-\phi(x;a)I+\int\mathcal{R}(x,x^{\prime};a)I(\eta,r,x^{\prime},\mu^{\prime})dx^{\prime}d\mu^{\prime}/2$
$\displaystyle-\kappa(x)I+A\kappa(x)\int\mathcal{R}^{d}(x,x^{\prime};\mu,\mu^{\prime};a)I(\eta,r,x^{\prime},\mu^{\prime})dx^{\prime}d\mu^{\prime}+S$

where $I(t,r_{p},x,\mu)$ is the specific intensity, which is a function of
time $t$, radial coordinate $r_{p}$, frequency $x$ and the direction angle,
$\mu=\cos\theta$, with respect to the radial vector ${\bf r}$.

In eq.(3), we use the dimensionless time $\eta$ defined as $\eta=cn_{\rm
HI}\sigma_{0}t$ and the dimensionless radial coordinate $r$ defined as
$r=n_{\rm HI}\sigma_{0}r_{p}$. That is, $\eta$ and $r$ are, respectively, in
the units of mean free flight-time and mean free path of photon $\nu_{0}$ with
respect to the resonant scattering without dust scattering and absorption.
Without resonant scattering, a signal propagates in the radial direction with
the speed of light, the orbit of the signal is then $r=\eta+{\rm const}$. With
dimensionless variable, the size of the halo $R$ is equal to $\tau_{0}$.

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

$\displaystyle\mathcal{R}(x,x^{\prime};a)=$
$\displaystyle\frac{1}{\pi^{3/2}}\int^{\infty}_{|x-x^{\prime}|/2}e^{-u^{2}}\left[\tan^{-1}\left(\frac{x_{\min}+u}{a}\right)-\tan^{-1}\left(\frac{x_{\max}-u}{a}\right)\right]du$

where $x_{\min}=\min(x,x^{\prime})$ and $x_{\max}=\max(x,x^{\prime})$. In the
case of $a=0$, i.e. considering only the Doppler broadening, the re-
distribution function is

$\mathcal{R}(x,x^{\prime})=\frac{1}{2}{\rm erfc}[{\rm
max}(|x|,|x^{\prime}|)].$ (5)

The re-distribution function of equation (5) is normalized as
$\int_{-\infty}^{\infty}\mathcal{R}(x,x^{\prime})dx^{\prime}=\phi(x,0)=\pi^{-1/2}e^{-x^{2}}$.
With this normalization, the total number of photons is conserved in the
evolution described by equation (3). That is, the destruction processes of
Ly$\alpha$ photons, such as the two-photon process (Spitzer & Greenstein 1951;
Osterbrock 1962), are ignored in equation (3). The recoil of atoms is also not
considered in equation (4) or (5). The effect of recoil actually is under
control (Roy et al. 2009c, 2010). We will address it in next section.

The absorption and scattering of dust are described by the term $\kappa(x)I$
of eq.(3), where $\kappa(x)=\sigma_{d}/\sigma_{0}$, which is of the order of
$10^{-8}(T/10^{4})^{1/2}$ (Draine & Lee 1984; Draine 2003). The term with $A$
of eq.(3) describes albedo, i.e. $A\equiv\sigma_{s}/\sigma_{d}$, where
$\sigma_{s}$ is the cross section of dust scattering. Generally, $A$ lies
approximately between 0.3 and 0.4 (Pei 1992; Weingartner & Draine 2001).

Since dust generally is much heavier than a single atoms, the recoil of dust
particles can be neglected when colliding with a photon. Under this “heavy
dust” approximation, photons do not change their frequency during the
collision with dust. The redistribution function of dust $\mathcal{R}^{d}$ is
independent of $x$ and $x^{\prime}$, and is simply given by a phase function
as

$\mathcal{R}^{d}(\mu,\mu^{\prime})=\frac{1}{4\pi}\int^{2\pi}_{0}d\phi^{\prime}\frac{1-g^{2}}{(1+g^{2}-2g\bar{\mu})^{3/2}}=\sum_{l=0}^{\infty}\frac{(2l+1)}{2}g^{l}P_{l}(\mu)P_{l}(\mu^{\prime}),$
(6)

where $\bar{\mu}=\mu\mu^{\prime}+\sqrt{(1-\mu^{2})(1-\mu^{\prime
2})}\rm{cos}\phi^{\prime}$ and $P_{l}$ is the Legendre function. The factor
$g$ in eq.(6) is the asymmetry parameter. For isotropic scattering, $g=0$. The
cases of $g=+1$ and -1 correspond to complete forward and backward scattering,
respectively. Generally, the factor $g$ is a function of the wavelength. For
the Ly$\alpha$ photon, we will take $g=0.73$ for realistic dust scattering (Li
& Draine 2001). The integral of eq.(6) is performed in Appendix A.

In eq. (3), the term with the parameter $\gamma$ is due to the expansion of
the universe. If $n_{\rm H}$ is equal to the mean of the number density of
cosmic hydrogen, we have $\gamma=\tau_{GP}^{-1}$, and $\tau_{GP}$ is the Gunn-
Peterson optical depth. Since the Gunn-Peterson optical depth is of the order
of $10^{6}$ at high redshift (e.g. Roy et al. 2009c), the parameter $\gamma$
is of the order of $10^{-5}-10^{-6}$. Therefore, if the optical depth of halos
is equal to or less than 106, the term with $\gamma$ of eq.(3) can be ignored.

In eq.(3) we neglect the effect of collision transition from $H(2p)$ state to
$H(2s)$ state, which can significantly affect on the escape of Ly$\alpha$
photons when HI column density is higher than $10^{21}$ cm-2 and dust
absorption is very small (Neufeld, 1990). This generally is out of the
parameter range used below. We are also not considering the effects of bulk
motion of the medium of halos (e.g. Spaans & Silk 2006, Xu & Wu, 2010).

###  2.2 Eddington approximation

Eq.(6) indicates that the transfer equation (3) can be solved with the
Legendre expansion $I(\eta,r,x,\mu)=\sum_{l}I_{l}(\eta,r,x)P_{l}(\mu)$. If we
take only the first two terms, $l=0$ and 1, it is the Eddington approximation
as

$I(\eta,r,x,\mu)\simeq J(\eta,r,x)+3\mu F(\eta,r,x)$ (7)

where

$J(\eta,r,x)=\frac{1}{2}\int_{-1}^{+1}I(\eta,r,x,\mu)d\mu,\hskip
14.22636ptF(\eta,r,x)=\frac{1}{2}\int_{-1}^{+1}\mu I(\eta,r,x,\mu)d\mu.$ (8)

They are, respectively, the angularly averaged specific intensity and flux.
Defining $j=r^{2}J$ and $f=r^{2}F$, Eq.(3) yields the equations of $j$ and $f$
as

$\displaystyle{\partial j\over\partial\eta}+\frac{\partial f}{\partial r}$
$\displaystyle=$ $\displaystyle-(1-A)\kappa
j-\phi(x;a)j+\int\mathcal{R}(x,x^{\prime};a)jdx^{\prime}+\gamma\frac{\partial
j}{\partial x}+r^{2}S,$ (9) $\displaystyle\frac{\partial
f}{\partial\eta}+\frac{1}{3}\frac{\partial j}{\partial
r}-\frac{2}{3}\frac{j}{r}$ $\displaystyle=$ $\displaystyle-(1-Ag)\kappa
f+\gamma\frac{\partial f}{\partial x}-\phi(x;a)f.$ (10)

The mean intensity $j(\eta,r,x)$ describes the $x$ photons trapped in the
position $r$ at time $\eta$ by the resonant scattering, while the flux
$f(\eta,r,x)$ describes the photons in transit.

The source term $S$ in the equations (3) and (9) can be described by a
boundary condition of $j$ and $f$ at $r=r_{0}$. We can take $r_{0}=0$. Thus,
the boundary condition is

$j(\eta,0,x)=0,\hskip 28.45274ptf(\eta,0,x)=S_{0}\phi_{s}(x),$ (11)

where $S_{0}$, and $\phi_{s}(x)$ are, respectively, the intensity and
normalized frequency profile of the sources. Since equation (3) is linear, the
solutions of $j(x)$ and $f(x)$ for given $S_{0}=S$ are equal to $Sj_{1}(x)$
and $Sf_{1}(x)$, where $j_{1}(x)$ and $f_{1}(x)$ are the solutions of
$S_{0}=1$. On the other hand, the equation (3) is not linear with respect to
the function $\phi_{s}(x)$. The solution $f(x)$ for a given $\phi_{s}(x)$ is
not equal to $\phi_{s}(x)f_{1}(x)$, where $f_{1}(x)$ is the solution of
$\phi_{s}(x)=1$.

In the range outside the halo, $r>R$, no photons propagate in the direction
$\mu<0$. The boundary condition at $r=R$ given by $\int_{0}^{-1}\mu
J(\eta,R,x,\mu)d\mu=0$ is then (Unno 1955)

$j(\eta,R,x)=2f(\eta,R,x).$ (12)

There is no photon in the field before $t=0$. Therefore, the initial condition
is

$j(0,r,x)=f(0,r,x)=0.$ (13)

We will solve equations (9) and (10) with boundary and initial conditions
eqs.(11) - (13) by using the WENO solver (Roy et al. 2009a, b, c, 2010). Some
details of this method is given in Appendix B.

###  2.3 Dust models

We consider three models of the dust as follows:

I. pure scattering, $A=1,\ g=0.73$: dust causes only anisotropic scattering,
but no absorption;

II. scattering and absorption. $A=0.32,\ g=0.73$: dust causes both absorption
and anisotropic scattering.

III. pure absorption. $A=0$: dust causes only absorption, but no scattering;

Models I and III do not occur in reality. They are, however, helpful to reveal
the effects of pure scattering and absorption on the radiative transfer.

Since $\kappa(x)$ is on the order of $10^{-8}$, its effect will be significant
only for halos with halos with optical depth $\tau_{0}\geq 10^{6}$, and
ignorable for $\tau_{0}\leq 10^{5}$. To illustrate the dust effect, we use
halos of $R=\tau_{0}\leq 10^{4}$, and take larger $\kappa$ to be $\simeq
10^{-4}-10^{-2}$. We also assume that $\kappa$ is frequency-independent. We
consider below only the case of grey dust, i.e. $\kappa$ is independent of
frequency $x$. This certainly is not realistic dust. Yet, the frequency range
given in solution below mostly are in the range $|x|<4$. Therefore, the
approximation of grey dust would be proper if cross section of dust is not
significantly frequency dependent in the range $|x|<4$.

###  2.4 Numerical example: Wouthuysen-Field thermalization

As the first example of numerical solutions, we show the Wouthuysen-Field
(W-F) effect, which requires that the distribution of Ly$\alpha$ photons in
the frequency space should be thermalized near the resonant frequency
$\nu_{0}$. The W-F effect illustrates the difference between the analytical
solutions of the Fokker-Planck approximation and that of eqs. (9) and (10).
The former can not show the local thermalization (Neufeld 1990), while the
latter can (Roy et al. 2009b). All problems related to the W-F local thermal
equilibrium should be studied with the integro-differential equation (3).

Figure 1: The mean intensity $j(\eta,r,x)$ at $\eta=500$ and $r=100$ for dust
models I (left panel), II (middle panel) and III (right panel). The source is
$S_{0}=1$ and $\phi_{s}(x)=(1/\sqrt{\pi})e^{-x^{2}}$. The parameter
$a=10^{-3}$. In each panel, $\kappa$ is taken to be 0, 10-4, 10-3 and 10-2.

Figure 1 presents a solution of mean intensity $j(\eta,r,x)$ at time radial
$\eta=500$ coordinate $r=10^{2}$ for halo with size $R\gg r=10^{2}$. The three
panels correspond to dust models I (left panel), II (middle panel) and III
(right panel). The source is taken to have a Gaussian profile
$\phi_{s}(x)=(1/\sqrt{\pi})e^{-x^{2}}$ and unit intensity $S_{0}=1$. The
solutions of Figure 1 actually are independent of $R$, if $R\gg 10^{2}$. The
intensity of $j$ is decreasing from left to right in Figure 1, because the
absorption is increasing with the models from I to III.

A remarkable feature shown in Figure 1 is that all $j(\eta,r,x)$ have a flat
plateau in the range $|x|\leq 2$. This gives the frequency range of the W-F
local thermalization (Roy et al, 2009b, c). The range of the flat plateau
$|x|\leq 2$ is almost dust-independent, either for model I or for models II
and III. This is expected, as neither the absorption nor scattering given by
the $\kappa$ term of eq.(3) changes the frequency distribution of photons. The
redistribution function (6) also does not change the frequency distribution of
photons. This point can also be seen from eqs.(9) and (10), in which the
$\kappa$ terms are frequency-independent. The evolution of the frequency
distribution of photons is due only to the resonant scattering.

Since thermalization will erase all frequency features within the range
$|x|\leq 2$, the double-peaked structure does not retain information of the
photon frequency distribution within $|x|<2$ at the source. That is, the
results in Figure 1 will hold for any source $S_{0}\phi_{s}(x)$ with arbitrary
$\phi_{s}(x)$ which is non-zero within $|x|<2$ (Roy et al. 2009b, c). This
property can also be used as a test of the simulation code. It requires that
simulation results of flat plateau should be hold, regardless of the source to
be monochromatic or with finite width around $\nu_{0}$.

##  3 Dust effects on photon escape

###  3.1 Model I: scattering of dust

To study the effects of dust scattering on the Ly$\alpha$ photon escape, we
show in Figure 2 the flux $f(\eta,r,x)$ of Ly$\alpha$ photons emergent from
halos at the boundary $r=R=10^{2}$ for Model I. The three panels of Figure 2
correspond to $\kappa=10^{-4}$, 10-3, and 10-2 from left to right,
respectively. The source starts to emit photons at $\eta=0$ with a stable
luminosity $S_{0}=1$, and with a Gaussian profile
$\phi_{s}(x)=(1/\sqrt{\pi})e^{-x^{2}}$.

Figure 2: Flux $f(\eta,r,x)$ of Ly$\alpha$ photons emergent from halos at the
boundary $R=10^{2}$, and for the dust model I $A=1,\ g=0.73$. The parameter
$\kappa$ is taken to be 10-4 (left), 10-3(middle) and 10-2 (right). The source
is $S_{0}=1$ and $\phi_{s}(x)=(1/\sqrt{\pi})e^{-x^{2}}$. The parameter
$a=10^{-3}$.

Figure 2 clearly shows that the time-evolution of $f(\eta,r,x)$ is
$\kappa$-independent. Although the cross section of dust scattering increases
about 100 times from $\kappa=10^{-4}$ to $\kappa=10^{-2}$, the curves of the
left and right panels in Figure 2 actually are almost identical.

According to the scenario of “single longest excursion”, photon escape is not
a process of Brownian random walk in the spatial space, but a transfer in the
frequency space (Osterbrock 1962; Avery & House 1968; Adams, 1972, 1975;
Harrington 1973; Bonilha et al. 1979). Photon will escape, once its frequency
is transferred from $|x|<2$ to $|x|>2$, on which the medium is transparent. On
the other hand, dust scattering given by the redistribution function eq.(6)
does not change photon frequency. Dust scattering has no effect on the
transfer in the frequency space.

Moreover, photons with frequency $|x|<2$ are quickly thermalized after a few
hundred resonant scattering. In the local thermal equilibrium state, the
angular distribution of photons is isotropic. Thus, even if the dust
scattering is anisotropic $g\neq 0$ with respect to the direction of the
incident particle, the angular distribution will keep isotropic undergoing a
$g\neq 0$ scattering. Hence, dust scattering also has no effect on the angular
distribution.

###  3.2 Model III: absorption of dust

Similar to Figure 2, we present in Figure 3 the flux of Model III, i.e. dust
causes only absorption without scattering. All other parameters of Figure 3
are the same as in Figure 2. In the left panel of Figure 3, the curves at the
time $\eta=2000$ and $3000$ are the same. It means the flux $f(\eta,R,x)$ at
the boundary $R$ is already stable, or saturated at the time $\eta\geq 2000$.
The small difference between the curves of $\eta=1000$ and $\eta\geq 2000$ of
the left panel indicates that the flux is still not yet completely saturated
at the time $\eta=1000$. However, comparing the middle and right panels of
Figure 3, we see that for $\kappa=10^{-3}$, the flux has already saturated at
$\eta=1600$, while it has saturated at $\eta=800$ for $\kappa=10^{-2}$. That
is, the stronger the dust absorption, the shorter the saturation time scale.
The time scales of escape or saturation do not increase by dust absorption,
and even decrease with respect to the medium without dust. Stronger absorption
leads to shorter time scale of saturation.

Figure 3: Flux $f(\eta,r,x)$ of Ly$\alpha$ photons emergent from halos at the
boundary $r=R=10^{2}$. The parameters of the dust are $A=0$ and $\kappa=$ 10-4
(left), 10-3 (middle) and 10-2 (right). Other parameters are the same as in
Figure 2.

Obviously, dust absorption does not help in producing photons for the “single
longest excursion”. Therefore, dust absorption can not make the time scale of
producing photons for “single longest excursion” to be smaller. However, dust
absorptions are effective in reducing the number of photons trapped in the
state of local thermalized equilibrium $|x|<2$ (see also §4.2). This leads to
the fact that the higher the value of $\kappa$, shorter the time scale of
saturation.

###  3.3 Effective absorption optical depth

Since Ly$\alpha$ photons underwent a large number of resonant scattering
before escaping from halo with optical depth $\tau_{0}\gg 1$, it is generally
believed that a small absorption of dust will lead to a significant decrease
of the flux. However, it is still unclear what the exact relationship between
the dust absorption and the resonant scattering is. This problem should be
measured by the effective optical depth of dust absorption of Ly$\alpha$
photons in $R=\tau_{0}\gg 1$ halos.

To calculate the effective optical depth, we first give the total flux of
Ly$\alpha$ photons emergent from halo of radius $R$, which is defined as
$F(\eta)=\int f(\eta,R,x)dx$. Figure 4 plots $F(\eta)$ as a function of time
$\eta$ for halo with sizes $R=\tau_{0}=10^{2}$ and $10^{4}$. The curves
typically are the time-evolution of growing and then saturating. The three
panels correspond to the dust models I, II and III from left to right. The
upper panels are of $R=10^{2}$, and lower panels for $R=10^{4}$. In each panel
of $R=10^{2}$, we have three curves corresponding to $\kappa=$ 10-4, 10-3 and
10-2, respectively. In cases of $R=10^{4}$, we take $\kappa=$ 10-4 and 10-3.

Figure 4: The time evolution of the total flux $F(\eta)$ at the boundary of
halos with $R=\tau_{0}=10^{2}$ (upper panels), and $R=\tau_{0}=10^{4}$ (lower
Panels). The source of $S_{0}=1$ and $\phi_{s}(x)=(1/\sqrt{\pi})e^{-x^{2}}$
starts to emit photons at time $\eta=0$. The parameters of dust are
$(A=1,g=0.73)$ (left); $(A=0.32,\ g=0.73)$ (middle) and $A=0$ (right). In each
panel of $R=10^{2}$, $\kappa$ is taken to be 10-4, 10-3 and 10-2. In the cases
of $R=10^{4}$, $\kappa$ is taken to be 10-4, 10-3.

The left panel of Figure 4 shows that the three curves of $\kappa=10^{-4}$,
10-3 and 10-2 are almost the same. This is consistent with Figure 2 that for
Model I, the time-evolution of $f$ are $\kappa$-independent for the pure
scattering dust. For the pure absorption dust (the right panel of Figure 4),
the saturated flux is smaller for larger $\kappa$. We can also see from Figure
4 that the time scale of approaching saturation is smaller for larger
$\kappa$. The result of model II is in between that for models I and III.

With the saturated flux of Figure 4, one can define the effective absorption
optical depth by $\tau_{\rm effect}\equiv-(1/\kappa)\ln f_{S}$. The results
are shown in Table 1, in which $\tau_{a}$ is the dust absorption depth. It is
interested to see that the effective absorption optical depth is always equal
to about a few times of the optical depth of resonant scattering $\tau_{0}$,
regardless whether $\tau_{a}$ is less than 1. Namely, the effective absorption
depth $\tau_{\rm effect}$ of dust is roughly proportional to $\tau_{0}$.

Table 1. Effective absorption optical depth $\tau_{\rm effect}$

| Model II | Model III
---|---|---
$R=\tau_{0}$ | $\kappa$ | $\tau_{a}$ | $f_{S}$ | $\tau_{\rm effect}$ | $\tau_{a}$ | $f_{S}$ | $\tau_{\rm effect}$
$10^{2}$ | $10^{-4}$ | 0.0068 | 0.978 | $2.2\times 10^{2}$ | 0.01 | 0.963 | $3.8\times 10^{2}$
$10^{2}$ | $10^{-3}$ | 0.068 | 0.760 | $2.7\times 10^{2}$ | 0.10 | 0.670 | $4.0\times 10^{2}$
$10^{2}$ | $10^{-2}$ | 0.68 | 0.116 | $2.2\times 10^{2}$ | 1.00 | 0.057 | $2.9\times 10^{2}$
$10^{4}$ | $10^{-4}$ | 0.68 | $6.28\times 10^{-2}$ | $2.8\times 10^{4}$ | 1.00 | $3.02\times 10^{-2}$ | $3.5\times 10^{4}$
$10^{4}$ | $10^{-3}$ | 6.8 | $4.07\times 10^{-7}$ | $1.5\times 10^{4}$ | 10.0 | $2.87\times 10^{-9}$ | $1.97\times 10^{4}$

According to the random walk scenario, if a medium has optical depths of
absorption $\tau_{a}$ and scattering $\tau_{s}$, the effective absorption
optical depth should be equal to $\tau_{\rm
effect}=\sqrt{\tau_{a}(\tau_{a}+\tau_{s})}$ (Rybicki & Lightman 1979).
However, the results of the last line of Table 1 show that the random walk
scenario does not work for the dust effect on resonant photon transfer. This
result is consistent with Figures 2 and 3. When optical depth of dust is lower
than the optical depth of resonant scattering $\tau_{0}$, the time scale of
photon escaping basically is not affected by the dust, but is proportional to
$\tau_{0}$, and therefore, the absorption is also proportional to $\tau_{0}$.

###  3.4 Escape coefficient

With the total flux, we can define the escaping coefficient of Ly$\alpha$
photon as $f_{\rm esc}(\eta,\tau_{0})\equiv F(\eta)/F_{0}$, where $F_{0}$ is
the flux of the center source. Figure 5 shows $f_{\rm esc}(\eta,\tau_{0})$ at
three times $\eta=5\times 10^{3}$, 104 and 3.2$\times 10^{4}$ for Model II and
$\kappa=10^{-3}$. At $\eta=5\times 10^{3}$, the flux of halos with
$\tau_{0}\leq 10^{3}$ is saturated. At $\eta=10^{4}$, halos with $\tau_{0}\leq
3\times 10^{3}$ are saturated, and all halos of $\tau_{0}\leq 10^{4}$ are
saturated at $\eta=3.2\times 10^{4}$.

Figure 5: Escaping coefficient $f_{\rm esc}(\eta)$ as a function of the
optical depth $\tau_{0}$ of halo at time $\eta=5\times 10^{3}$, 104, and
3.2$\times 10^{4}$ from bottom to up. Dust is modeled by II, $A=0.32,\
g=0.73$, and $\kappa=10^{-3}$.

##  4 Dust effects on double-peaked profile

###  4.1 Dust and the frequency of double peaks

A remarkable feature of Ly$\alpha$ photon emergent from optically thick medium
is the double-peaked profile. Figures 1, 2 and 3 have shown that the double
peak frequencies $x_{+}=|x_{-}|$ are almost independent of either the
scattering or the absorption of dust. In this section, we consider halos with
size $R$ or $\tau_{0}$ larger than $10^{2}$. Figure 6 presents the double peak
frequency $|x_{\pm}|$ as a function of $a\tau_{0}$, where the parameter $a$ is
taken to be 10-2 (left) and $5\times 10^{-3}$ (right). Comparing the curves
with dust and without dust in Figure 6 we can say that the dust effect on
$|x_{\pm}|$ is very small till $a\tau_{0}=aR=10^{2}$.

Figure 6: The two-peak frequencies $x_{+}=|x_{-}|$ as a function of
$a\tau_{0}$. The parameter $a$ is taken to be 10-2 (left) and $5\times
10^{-3}$ (right). Dust model III (pure absorption) is used, and $\kappa$ is
taken to be $10^{-3}$. The dashed straight line gives $\log x_{\pm}$-$\log
a\tau$ with slope 1/3, which is to show the $(a\tau)^{1/3}$-law of $x_{\pm}$.

In the range $a\tau_{0}<20$, the $|x_{\pm}|$-$\tau_{0}$ relation is almost
flat with $|x_{\pm}|\simeq 2$. It is because the double-peaked profile is
given by the frequency range of the locally thermal equilibrium. The positions
of the two peaks, $x_{+}$ and $x_{-}$, basically are at the maximum and
minimum frequencies of the local thermalization. The frequency range of the
local thermal equilibrium state is mainly determined by the Doppler
broadening, and weakly dependent on $\tau_{0}$. Thus, we always have
$x_{\pm}\simeq\pm 2$. When the optical depth is larger, $a\tau_{0}\sim
10^{2}$, more and more photons of the flux are attributed to the resonant
scattering by the Lorentzian wing of the Voigt profile. In this phase,
$|x_{\pm}|$ will increase with $\tau_{0}$.

Figure 6 shows also a line $x_{\pm}=\pm(a\tau_{0})^{1/3}$, which is given by
the analytical solution of the Fokker-Planck approximation, in which the
Doppler broadening core in the Voigt profile is ignored (Harrington 1973,
Neufeld 1990, Dijkstra 2006). The numerical solutions of eqs (3) or (9) and
(10) deviate from the $(a\tau_{0})^{1/3}$-law at all parameter range of Figure
6. The deviation at $a\tau_{0}<20$ is due to that the Doppler broadening core
in the Voigt profile is ignored in the Fokker-Planck approximation, and then,
no locally thermal equilibrium can be reached. Therefore, in the range
$a\tau_{0}<20$, $|x_{\pm}|$ of the WENO solution is larger than the
$(a\tau_{0})^{1/3}$-law. In the range of $a\tau_{0}>20$, the Fokker-Planck
approximation yields a faster diffusion of photons in the frequency space.
This point can be seen in the comparison between a Fokker-Planck solution with
Field’s analytical solution (Figure 1 in Rybicki & Dell’Antonio 1994). In this
range, the numerical results of $|x_{\pm}|$ is less than the
$(a\tau_{0})^{1/3}$-law.

###  4.2 No narrowing and no widening

The dust effect has been used to explain the narrowing of the width between
the two peaks (Laursen et al. 2009). Oppositely, it is also used to explain
the widening of the width between the two peaks (Verhamme et al. 2006).
However, Figures 1, 2, 3 and 6 already show that the width between the two
peaks of the profile is very weakly dependent on dust scattering and
absorption. This result supports, at least in the parameter range considered
in Figures 1, 2, 3, neither the narrowing nor the widening of the two peaks.

Figure 7: $\ln[f(\eta,r,x,\kappa=0)/f(\eta,r,x,\kappa)]$ as function of $x$
for model II (up), and III (bottom), and $\kappa=10^{-3}$ (left) and 10-2
(right). Other parameters are the same as in Figure 2.

If dust absorption can cause narrowing, the absorption should be weaker at
$|x|\sim 0$, and stronger at $|x|\geq 2$. Similarly, if dust absorption can
cause widening, the absorption should be weaker at $|x|\sim 2$, and stronger
at $|x|\sim 0$. To test these assumptions, Figure 7 plots
$\ln[f(\eta,r,x,\kappa=0)/f(\eta,r,x,\kappa)]$ as a function of $x$. It
measures the $x$(frequency)-dependence of the flux ratio with and without dust
absorption. We take large $\eta$, and then the fluxes in Figure 7 are
saturated. Figure 7 shows that the absorption in the range $|x|<2$ is much
stronger than that of $|x|>2$, and therefore, the assumption of the narrowing
is ruled out. Figure 7 shows also that the curves of
$\ln[f(\eta,r,x,\kappa=0)/f(\eta,r,x,\kappa=10^{-3})]$ are almost flat in the
range $|x|<2$. Therefore, the assumption of widening of the two peaks can also
be ruled out.

Since the cross sections of dust absorption and scattering are assumed to be
frequency-independent. Eqs. (9) and (10) do not contain any frequency scales
other than that from resonant scattering. However, either narrowing or
widening would require to have frequency scales different from that of
resonant scattering. This is occurence is not possible if the dust is gray.

###  4.3 Profile of absorption spectrum

If the radiation from the sources has a continuum spectrum, the effect of
neutral hydrogen halos is to produce an absorption line at $\nu=\nu_{0}$. The
profile of the absorption line can also be found by solving equations (9) and
(10), but replacing the boundary equation (11) by

$j(\eta,0,x)=0,\hskip 28.45274ptf(\eta,0,x)=S_{0}.$ (14)

That is, we assume that the original spectrum is flat in the frequency space.
The spectrum of the flux emergent from halo of $R=10^{2}$ and $10^{4}$ with
central source of eq.(14) for dust models I, II and III are shown in Figure 8.
All curves are for large $\eta$, i.e. they are saturated.

Figure 8: The spectrum of the flux emergent from halo of $R=10^{2}$ (upper
panels) and 104 (lower panels) with central source of eq.(14) for the dust
model I (left), II (middle) and III (right). Other parameters are the same as
in Figure 2.

The optical depths at the frequency $|x|>4$ are small, and therefore, the
Eddington approximation might no longer be proper. However, those photons do
not strongly involve the resonant scattering, and hence they do not
significantly affect the solution around $x=0$. The solutions of Figure 8 is
still useful to study the profiles of $f$ around $x=0$.

The flux profile of Figure 8 typically are absorption lines with width given
by the double peaks similar to the double peaked structure of the emission
line. The flux at the double peaks is even higher than the flat wing. It is
because more photons are stored in the frequency range $|x|<2$. According to
the redistribution function eq.(4), the probability of transferring a
$x^{\prime}$ photon to a $|x|<|x^{\prime}|$ photon is larger than that from
$|x^{\prime}|$ to $|x|>|x^{\prime}|$. Therefore, if the original spectrum is
flat, the net effect of resonant scattering is to bring photons with frequency
$|x|>2$ to $|x|<2$. Photons stored $|x|<2$ are thermalized, and therefore, in
the range $|x|<2$, the profile will be the same as the emission line, and the
double peaks can be higher than the wing. It makes the shoulder at $|x|\sim
2$.

As expected, for model I (left panels of Figure 8), the double profile is
completely $\kappa$-independent. Dusty scattering does not change the flux and
its profile. For models II and III, the higher the $\kappa$, the lower the
flux of the wing, because the dust absorption is assumed to be frequency-
independent. The positions of the double peaks, $x$, in the absorption
spectrum are also $\kappa$-independent. This once again shows that dust
absorption and scattering causes neither narrowing nor widening of the double-
peaked profile. However, for higher $\kappa$ the flux of the peaks is lower.
When the absorption is very strong, the double-peaked structure might
disappear, but will never be narrowed or widened.

##  5 Discussions and conclusions

The study of dust effects on radiative transfer has had a long history related
to extinction. However, dust effects on radiative transfer of resonant photons
actually have not been carefully investigated. Existing works are mostly based
on the solutions of the Fokker-Planck approximation, or Monte Carlo
simulation. These results are important. We revisited these problems with the
WENO solver of the integro-differential equation of the resonant radiative
transfer, and have found some features which have not been addressed in
previous works. These features are summarized as follows.

First, the random walk picture in the physical space will no longer be
available for estimating the effective optical depth of dust absorption. For a
medium with the optical depth of absorption and resonant scattering to be
$\tau_{a}\gg 1$, $\tau(\nu_{0})\gg 1$ and $\tau_{s}(\nu_{0})\gg\tau_{a}$, the
effective absorption optical depth is found to be almost independent of
$\tau_{a}$, and to be equal to about a few times of $\tau_{s}(\nu_{0})$.

Second, dust absorption will, of course, yield the decrease of the flux of
Ly$\alpha$ photons emergent from optical thick medium. However, if the
absorption cross-section of dust is frequency independent, the double-peaked
structure of the frequency profile is basically dust-independent. The double-
peaked structure does not get narrowed or widened by the absorption and
scattering of dust.

Third, the time scales of Ly$\alpha$ photon transfer basically are independent
of dust scattering and absorption. It is because those time scales are mainly
determined by the kinetics in the frequency space, while dust does not affect
the behavior of the transfer in the frequency space if the cross section of
the dust is wavelength-independent. The local thermal equilibrium makes the
anisotropic scattering to be ineffective on the angular distribution of
photons. Dust absorption and scattering do not lead to the increase or
decrease of the time of storing Ly$\alpha$ photons in the halos.

The differences between the time-independent solutions of the Fokker-Planck
approximation, or Monte Carlo simulation and the WENO solution of eq.(3) is
mainly related to the W-F effect. Therefore, all above-mentioned features can
already be clearly seen with halos of $\tau_{0}\sim 10^{2}$, in which the W-F
local thermal equilibrium has been well established.

In this context, most calculation in this paper is on holes with
$\tau_{0}<10^{5}$. This range of $\tau_{0}$ certainly is unable to describe
halos with column number density of HI larger than 1017 cm-2 (e.g. Roy et al.
2010). Nevertheless, the result of $\tau_{0}<10^{5}$ would already be useful
for studying the 21 cm region around high-redshift sources, of which the
optical depth typically is (Liu et al 2007; Roy, et al. 2009c).

$\tau_{0}=3.9\times 10^{5}f_{\rm HI}\left(\frac{T}{10^{4}{\rm
K}}\right)^{-1/2}\left(\frac{1+z}{10}\right)^{3}\left(\frac{\Omega_{b}h^{2}}{0.022}\right)\left(\frac{R_{\rm
ph}}{10{\rm kpc}}\right),$ (15)

where $f_{\rm HI}$ is the fraction of HI. All other parameters in eq. (15) is
taken from the concordance $\Lambda$CDM mode. For these objects the relation
between dimensionless $\eta$ and physical time $t$ is given by

$t=5.4\times 10^{-2}f^{-1}_{\rm HI}\left(\frac{T}{10^{4}{\rm
K}}\right)^{1/2}\left(\frac{1+z}{10}\right)^{-3}\left(\frac{\Omega_{b}h^{2}}{0.022}\right)^{-1}\eta,\
\ {\rm yr}.$ (16)

The 21 cm emission rely on the W-F effect. On the other hand, the time-scale
of the evolution of the 21 region is short. The effect of dust on the time-
scales of Ly$\alpha$ evolution should be considered.

We have not considered the Ly$\alpha$ photons produced by the recombination in
the ionized halo. If the halo is optical thick, photons from the recombination
will also be thermalized. The information of where the photon comes from will
be forgotten during the thermalization. Therefore, photons from recombination
should not show any difference from those emitted from central sources. Only
the photons formed at the place very close to the boundary of the halo will
not be thermalized, and may yield different behavior.

This research is partially supported by ARO grants W911NF-08-1-0520 and
W911NF-11-1-0091.

##  Appendix A Integral of the phase function [eq.(6)]

Eq.(6) can be rewritten as

$\mathcal{R}^{d}(\mu,\mu^{\prime})=\frac{1}{4\pi}\int^{2\pi}_{0}d\phi^{\prime}\frac{1-g^{2}}{|{\bf
I}-g{\bf I^{\prime}}|^{\frac{3}{2}}}$ (A1)

where ${\bf I}$ and ${\bf I^{\prime}}$ are unit vector on the direction of
polar angle $\theta$ and $\theta^{\prime}$, and azimuth angle $\phi$ and
$\phi^{\prime}$, respectively. That is ${\bf I}\cdot{\bf I}={\bf
I^{\prime}}\cdot{\bf I^{\prime}}=1$ and ${\bf I}\cdot{\bf
I^{\prime}}=\cos\gamma=\cos\theta\cos\theta^{\prime}+\sin\theta\sin\theta^{\prime}\cos(\phi-\phi^{\prime})$,
and $\mu=\cos\theta$, $\mu^{\prime}=\cos\theta$. We have

$\displaystyle\frac{d}{dg}\frac{1}{|{\bf I}-g{\bf
I^{\prime}}|^{1/2}}=\frac{1-g^{2}}{2g|{\bf I}-g{\bf
I^{\prime}}|^{3/2}}-\frac{1}{2g|{\bf I}-g{\bf I^{\prime}}|^{1/2}},$ (A2)

and therefore,

$\frac{1-g^{2}}{|{\bf I}-g{\bf
I^{\prime}}|^{3/2}}=2g\frac{d}{dg}\frac{1}{|{\bf I}-g{\bf
I^{\prime}}|^{1/2}}+\frac{1}{|{\bf I}-g{\bf I^{\prime}}|^{1/2}}.$ (A3)

The expansion with Legendre functions $P_{l}(\cos\gamma)$ gives

$\frac{1}{|{\bf I}-g{\bf
I^{\prime}}|^{1/2}}=\sum_{l=0}^{\infty}g^{l}P_{l}(\cos\gamma),$ (A4)

and then

$\frac{1-g^{2}}{|{\bf I}-g{\bf
I^{\prime}}|^{3/2}}=\sum_{l=1}^{\infty}2lg^{l}P_{l}(\cos\gamma)+\sum_{l=0}^{\infty}g^{l}P_{l}(\cos\gamma).$
(A5)

Since
$\cos\gamma=\cos\theta\cos\theta^{\prime}+\sin\theta\sin\theta^{\prime}\cos(\phi-\phi^{\prime})$,
we have the following identity for the Legendre function $P_{l}(\cos\gamma)$
as

$P_{l}(\cos\gamma)=P_{l}(\cos\theta)P_{l}(\cos\theta^{\prime})+2\sum_{m=1}^{m=l}\frac{(l-m)!}{(l+m)!}P^{m}_{l}(\cos\theta)P^{m}_{l}(\cos\theta^{\prime})\cos[m(\phi-\phi^{\prime})].$
(A6)

The integral of $\phi^{\prime}$ in eq.(A1) kills the second term of eq.(A6),
we have then

$\displaystyle\mathcal{R}^{d}(\mu,\mu^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}2\pi\left[\sum_{l=1}^{\infty}2lg^{l}P_{l}(\cos\theta)P_{l}(\cos\theta^{\prime})+\sum_{l=0}^{\infty}g^{l}P_{l}(\cos\theta)P_{l}(\cos\theta^{\prime})\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\sum_{l=1}^{\infty}2lg^{l}P_{l}(\mu)P_{l}(\mu^{\prime})+\sum_{l=0}^{\infty}g^{l}P_{l}(\mu)P_{l}(\mu^{\prime})\right].$

Using the orthogonal relation
$\int_{-1}^{1}P_{l}(\mu)P_{l^{\prime}}(\mu)d\mu=\frac{2}{2l+1}\delta_{l,l^{\prime}}$,
we have

$R_{0}(g)=\frac{1}{2}\int_{-1}^{1}d\mu\int_{-1}^{1}d\mu^{\prime}R^{d}(\mu,\mu^{\prime})=1,$
(A8)

for which only the term $l=0$ in eq.(A7) has contribution. Similarly,

$R_{1}(g)=\frac{1}{2}\int_{-1}^{1}d\mu\int_{-1}^{1}d\mu^{\prime}\mu
R^{d}(\mu,\mu^{\prime})=\frac{1}{2}\int_{-1}^{1}d\mu\int_{-1}^{1}d\mu^{\prime}\mu^{\prime}R^{d}(\mu,\mu^{\prime})=0,$
(A9)
$R_{2}(g)=\frac{1}{2}\int_{-1}^{1}d\mu\int_{-1}^{1}d\mu^{\prime}\mu\mu^{\prime}R^{d}(\mu,\mu^{\prime})=\frac{g}{3}.$
(A10)

These results are used in deriving eqs.(9) and (10).

##  Appendix B Numerical algorithm

To solve Equations $(9)$ and $(10)$ as a system, our computational domain is
$(r,x)\in[0,r_{\textrm{max}}]\times[x_{\textrm{left}},x_{\textrm{right}}],$
where $r_{\textrm{max}},x_{\textrm{left}}$ and $x_{\textrm{right}}$ are chosen
such that the solution vanishes to zero outside the boundaries. We choose mesh
sizes with grid refinement tests to ensure proper numerical resolution. In the
following, we describe numerical techniques involved in our algorithm,
including approximations to spatial derivatives, numerical boundary condition,
and time evolution.

###  B.1 The WENO Algorithm: Approximations to the Spacial Derivatives

The spacial derivative terms in Equation $(9)$ and $(10)$ are approximated by
a fifth-order finite difference WENO scheme.

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

$\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})$

with fixed $\eta=\eta^{n}$ 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. Denote

$h_{j}=j(\eta^{n},r_{i},x_{j}),\quad\quad j=-2,-1,\cdots,N+3$

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

$\hat{h}_{j+1/2}=\omega_{1}\hat{h}^{(1)}_{j+1/2}+\omega_{2}\hat{h}_{j+1/2}^{(2)}+\omega_{3}\hat{h}_{j+1/2}^{(3)}$

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

$\displaystyle\hat{h}_{j+1/2}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{1}{6}h_{j-1}+\frac{5}{6}h_{j}+\frac{1}{3}h_{j+1},$
$\displaystyle\hat{h}_{j+1/2}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{3}h_{j}+\frac{5}{6}h_{j+1}-\frac{1}{6}h_{j+2},$
$\displaystyle\hat{h}_{j+1/2}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{11}{6}h_{j+1}-\frac{7}{6}h_{j+2}+\frac{1}{3}h_{j+3}.$

and the nonlinear weights $\omega_{m}$ are given by

$\omega_{m}=\frac{\breve{\omega}_{m}}{\sum_{l=1}^{3}\breve{\omega}_{l}},\quad\breve{\omega}_{l}=\frac{\gamma_{l}}{(\epsilon+\beta_{l})^{2}}$

where $\epsilon$ is a parameter to avoid the denominator to become zero and is
taken as $\epsilon=10^{-8}$. The linear weights $\gamma_{l}$ are given by

$\gamma_{1}=\frac{3}{10},\quad\gamma_{2}=\frac{3}{5},\quad\gamma_{3}=\frac{1}{10},$

and the smoothness indicators $\beta_{l}$ are given by

$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\frac{13}{12}(h_{j-1}-2h_{j}+h_{j+1})^{2}+\frac{1}{4}(h_{j-1}-4h_{j}+3h_{j+1})^{2},$
$\displaystyle\beta_{2}$ $\displaystyle=$
$\displaystyle\frac{13}{12}(h_{j}-2h_{j+1}+h_{j+2})^{2}+\frac{1}{4}(h_{j}-h_{j+2})^{2},$
$\displaystyle\beta_{3}$ $\displaystyle=$
$\displaystyle\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}.$

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

${\bf u}_{t}+A{\bf u}_{r},$

where ${\bf u}=(j,f)^{T}$ and

$A=\left(\begin{array}[]{cc}0&1\\\ \frac{1}{3}&0\end{array}\right)$

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 $\Lambda$, $R$ and $R^{-1}$. We then project u to the local
characteristic fields v with ${\bf v}=R^{-1}{\bf u}.$ Now ${\bf u}_{t}+A{\bf
u}_{r}$ of the original system is decoupled as two independent equations as
${\bf v}_{t}+\Lambda{\bf v}_{r}.$ We approximate the derivative ${\bf 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 ${\bf v}_{r}$ back
to the physical space by ${\bf u}_{r}=R{\bf v}_{r}$. We refer the readers to
Cockburn et al. 1998 for more implementation details.

###  B.2 Numerical Boundary Condition

To implement the boundary condition $(12)$, we also need to perform a
characteristic decomposition as discussed above. Using the same notation as
before, we project u to the local characteristic fields v with ${\bf
v}=R^{-1}{\bf u}$. Denote ${\bf v}=(v_{1},v_{2})^{T}$, now ${\bf u}_{t}+A{\bf
u}_{r}$ of the original system is decoupled to two independent scalar
operators given by

$\frac{\partial v_{1}}{\partial t}+\lambda_{1}\frac{\partial v_{1}}{\partial
r};\qquad\frac{\partial v_{2}}{\partial t}+\lambda_{2}\frac{\partial
v_{2}}{\partial r}$

where $\lambda_{1}=\frac{\sqrt{3}}{3}$ and $\lambda_{2}=-\frac{\sqrt{3}}{3}$.
The characteristic line starting from the boundary $r=r_{\textrm{max}}$ for
the first equation is pointing outside the computational domain while the one
for the second equation is pointing inside. For well-posedness of our system,
we need to impose the boundary condition there as

$v_{2}=\alpha v_{1}+\beta$

with constants $\alpha$ and $\beta$. We can calculate the values of $\alpha$
and $\beta$ based on equation $(12)$ and the left and right eigenvectors of
$A$. For example, if we take

$R=\left(\begin{array}[]{cc}\frac{\sqrt{3}}{2}&\frac{\sqrt{3}}{2}\\\
\frac{1}{2}&-\frac{1}{2}\end{array}\right),$

we can compute that $\alpha=7-4\sqrt{3}$ and $\beta=0$. We use extrapolation
to obtain the value of $v_{1}$ and then compute the value $v_{2}$. In the end,
we transfer ${\bf v}$ back to the physical space by ${\bf u}=R{\bf v}$.

###  B.3 Time Evolution

To evolve in time, we use the third-order TVD Runge-Kutta time discretization
(Shu & Osher 1988). For system of ODEs $u_{t}=L(u)$, the third order Runge-
Kutta method is

$\displaystyle u^{(1)}$ $\displaystyle=$ $\displaystyle u^{n}+\Delta\tau
L(u^{n},\tau^{n}),$ $\displaystyle u^{(2)}$ $\displaystyle=$
$\displaystyle\frac{3}{4}u^{n}+\frac{1}{4}(u^{(1)}+\Delta\tau
L(u^{(1)},\tau^{n}+\Delta\tau)),$ $\displaystyle u^{n+1}$ $\displaystyle=$
$\displaystyle\frac{1}{3}u^{n}+\frac{2}{3}(u^{(2)}+\Delta\tau
L(u^{(2)},\tau^{n}+\frac{1}{2}\Delta\tau)).$

## References

  * (1)
  * (2) Adams, T.F. 1972, ApJ, 174, 439
  * (3)
  * (4) Adams, T.F. 1975, ApJ, 201, 350
  * (5)
  * (6) Adams, T.F., Hummer, D.G. & Rybicki, G.B. 1971, J. Quant. Spectrosc. Radiat. TransferApJ, 11, 1365
  * (7)
  * (8) Avery, L.W., & House, L.L. 1968, ApJ, 152, 493
  * (9)
  * (10) Blanc, G. A, et al. 2010, arXiv:1011.0430
  * (11)
  * (12) Bonilha, J. R. M., Ferch, R., Salpeter, E. E., Slater, G., Noerdlinger, P. D. 1979, ApJ233, 649
  * (13)
  * (14) Cockburn, B., Johnson, C., Shu, C-W. & Tadmor, E. 1998, Lecture Notes in Mathematics, 1697, 450
  * (15)
  * (16) Dijkstra, M., Haiman, Z. & Spaans, M. 2006, ApJ, 649, 14
  * (17)
  * (18) Dijkstra, M. & Loeb, A. 2009, MNRAS400, 1109
  * (19)
  * (20) Draine, B.T. 2003, ApJ, 598, 1017
  * (21)
  * (22) Draine, B. T.; Lee, H. M. 1984, ApJ, 285, 89
  * (23)
  * (24) Fang, L.Z., 2009, Inter. J. Mod. Phys. D18, 1943
  * (25)
  * (26) Fardal et al. 2001, ApJ, 526, 505
  * (27)
  * (28) Field, G.B., 1958, Proc. IRE, 46, 240
  * (29)
  * (30) Field, G.B. 1959, ApJ, 129, 551
  * (31)
  * (32) Haiman Z. et al. 2000, ApJ, 537, L5
  * (33)
  * (34) Haiman Z. & Spaans, M. 1999, ApJ, 518, 138
  * (35)
  * (36) Hansen, M., & Oh, S. P. 2006, MNRAS, 367, 979
  * (37)
  * (38) Harrington, J.P. 1973, MNRAS, 162, 43
  * (39)
  * (40) Hayes, M., et al. 2010, Nature, 464, 562
  * (41)
  * (42) Hayes, M., Schaerer, D., Östlin, G., Mas-Hesse, J. M., Atek, H. & Kunth, D. 2011, ApJ, 730, 8
  * (43)
  * (44) Henyey, L. G. & Greestein, J. L. 1941, AJ, 93, 70
  * (45)
  * (46) Hummer, D.G. 1962, MNRAS, 125, 21
  * (47)
  * (48) Hummer, D.G. 1965, Mem. R. astr. Soc., 70, 1
  * (49)
  * (50) Hummer, D.G. 1969, MNRAS, 145, 95
  * (51)
  * (52) Hummer, D.G. & Kunasz, P.B. 1980, ApJ, 236, 609
  * (53)
  * (54) Jiang, G. & Shu, C.-W. 1996, J. Comput. Phys., 126, 202
  * (55)
  * (56) Latif, M. et al. 2011, MNRAS413, L33
  * (57)
  * (58) Laursen, P., Sommer-Larsen, J., & Andersen, A., 2009, ApJ, 704, 1640
  * (59)
  * (60) Li, A., & Draine, B. T. 2001, ApJ, 550, 213
  * (61)
  * (62) Liu, et al. 2007, ApJ, 663, 1
  * (63)
  * (64) Neufeld, D. A. 1990, ApJ, 350, 216
  * (65)
  * (66) Neufeld, D. A. 1991, ApJ, 370, L85
  * (67)
  * (68) Osterbrock, D.E. 1962, ApJ, 135, 195
  * (69)
  * (70) Pei, Y. C. 1992, ApJ, 395, 130
  * (71)
  * (72) Qiu, J.-M., Feng, L.-L., Shu, C.-W. & Fang, L.-Z. 2006, New Astronomy, 12, 1
  * (73)
  * (74) Qiu, J.-M., Feng, L.-L., Shu, C.-W. & Fang, L.-Z. 2007, New Astronomy, 12, 398
  * (75)
  * (76) Qiu, J.-M., Shu, C.-W., Liu, J.-R. & Fang, L.-Z. 2008, New Astronomy, 13, 1
  * (77)
  * (78) Rybicki, G.B. 2006, ApJ, 674, 709
  * (79)
  * (80) Rybicki, G.B. & Dell’Antonio, I.P. 1994, ApJ, 427, 603
  * (81)
  * (82) Rybicki G.B. & Lightman, 1979 Radiative Processes in Astrophysics, (J. Wiley New York.
  * (83)
  * (84) Roy, I., Qiu J.-M., Shu C.-W. & Fang L.-Z., 2009a New Astronomy 14, 513
  * (85)
  * (86) Roy, I. Xu, W., Qiu J.-M., Shu C.-W. & Fang L.-Z., 2009b ApJ, 694, 1121
  * (87)
  * (88) Roy, I. Xu, W., Qiu J.-M., Shu C.-W. & Fang L.-Z., 2009c ApJ, 703, 1992
  * (89)
  * (90) Roy, I, Shu, C.-W. & Fang, L. Z. 2010, ApJ, 716, 604
  * (91)
  * (92) Shu, C.-W. & Osher, S., 1988, J. Comp. Phys., 77, 439
  * (93)
  * (94) Spaans, M. & Silk, J., 2006, ApJ, 652, 902
  * (95)
  * (96) Spitzer, L. & Greenstein, J.L. 1951, ApJ, 114, 407
  * (97)
  * (98) Stratta, G., Maiolino, R., Fiore, F. & D’Elia, V., 2007, 661, 9
  * (99)
  * (100) Unno, W. 1955, Publ. Astron. Soc. Japan. 7, 81
  * (101)
  * (102) Verhamme, A., Schaerer, D. & Maselli, A. 2006, AA, 460, 397
  * (103)
  * (104) Verhamme, A.; Schaerer, D., Atek, H. & Tapken, C. 2008, AA, 491, 89
  * (105)
  * (106) Weingartner, J. C.; Draine, B. T. 2001, ApJ, 548, 296
  * (107)
  * (108) Wouthuysen, S. A. 1952, AJ, 57, 31
  * (109)
  * (110) Xu, W. & Wu, X.-P. 2010, ApJ, 710, 1432
  * (111)