Table Of ContentA WENO algorithm for radiative transfer
with resonant scattering and the
Wouthuysen-Field Coupling
9
0
0 Ishani Roya, Jing-Mei Qiub, Chi-Wang Shua, Li-Zhi Fangc
2
n aDivision of Applied Mathematics, Brown University, Providence, RI 02912
a
J bMathematical and Computer Science, Colorado School of Mines, Golden, CO
5 80401
1
cDepartment of Physics, University of Arizona, Tucson, AZ 85721
]
O
C
. Abstract
h
p
-
We develop a numerical solver for the integral-differential equations, which de-
o
r scribe the radiative transfer of photon distribution in the frequency space with
t
s resonant scattering of Lyα photons by hydrogen gas in the early universe. The
a
[ time-dependent solutions of this equation is crucial to the estimation of the ef-
fect of the Wouthuysen-Field (WF) coupling in relation to the 21 cm emission and
1
v absorption at the epoch of reionization. However, the time-dependent solutions of
1 this equation have not yet been well performed. The resonant scattering leads to
7
the photon distribution in the frequency space to be piecewise smooth containing
1
2 sharp changes. The weighted essentially nonoscillatory (WENO) scheme is suitable
.
1 to handle this problem, as this algorithm has been found to be highly stable and
0 robust for solving Boltzmann equation. We test this numerical solver by 1.) the
9
analytic solutions of the evolution of the photon distribution in rest background;
0
: 2.) the analytic solution in expanding background, but without resonant scattering;
v
3.) the formation of local Boltzmann distribution around the resonant frequency
i
X
with the temperature to be the same as that of atom for recoil. We find that the
r
a evolution of the photon distribution due to resonant scattering with and without
recoil generally undergoes three phases. First, the profile of the photon distribution
is similar to the initial one. Second, an extremely flat plateau (without recoil) or
local Boltzmann distribution (with recoil) form around the resonant frequency, and
the width and height of the flat plateau or local Boltzmann distribution increase
withtime.Finally,thedistributionaroundtheresonantfrequencyissaturatedwhen
thephotonsfromthesourceisbalanced bytheredshiftoftheexpansion.Thisresult
indicates that the onset of the W-F coupling should not be determined by the third
phase, but by the time scale of the second phase. We found that the time scale of
the W-F coupling is equal to about a few hundreds of the mean free flight time
of photons with resonant frequency, and it basically is independent of the Sobolev
parameter if this parameter is much less than 1.
Preprint submitted to Elsevier 22 January 2009
Key words: cosmology: theory, radiation, hydrodynamics, methods: numerical
PACS: 95.30.Jx, 07.05.Tp, 98.80.-k
1 Introduction
It is generally believed that detecting redshifted 21 centimeter signals from
early universe is one of the next frontiers in observational cosmology, because
it would be able to provide the information of the first generation of light
sources in the cosmic dark ages. Many studies have been done on the 21 cm
emission and absorption from the halo of individual first stars (Chuzhoy et
al. 2006, Cen 2006, Liu et al. 2007). A common conclusion of these works
is that the configurations of the 21 cm emission and absorption regions is
strongly time-dependent. The reason is simple. The necessary conditions of
21 cm emission and absorption are that 1. the fraction of neutral hydrogen
(HI) is still high; 2, Lyα photons are available for the Wouthuysen-Field (W-
F) coupling. The region of 21 cm emission around first stars is then a thin
shell just outside the I-front (ionization-front). Hence the 21 cm emission shell
should move with a speed higher than the speed of the I-front v , which is
f
rather high, even comparable to the speed of light. Therefore, the time scale
of the formation and evolution of the regions of 21 cm signal can roughly be
estimated by d/v , d being the thickness of the 21 cm emission and absorption
f
shell. This time scale is found to be of the order of 1 Myr and less (Cen 2006,
Liu et al. 2007).
In all of the above-mentioned works, the spin temperature is calculated with
the assumption that the Wouthuysen-Field (W-F) coupling (Wouthuysen,
1952; Field, 1958, 1959) is effective. That is, the color temperature of photons
around Lyα frequency is assumed to be the same as the kinetic temperature
of hydrogen gas. The W-F coupling locks the color temperature of Lyα pho-
ton to the kinetic temperature of hydrogen gas, and then links the internal
(spin-) degree of freedom with the kinetic temperature of gas. Obviously, the
W-F coupling would be available if the time scale of the onset of the W-F
coupling is less than the time scale of the formation and evolution of the
21 cm emission and absorption shells. However, most calculations on the W-
F coupling are based on the time-independent Fokker-Planck approximation
(Chen & Miralda-Escude, 2004; Hirata, 2006; Furlanetto & Pritchard 2006;
Chuzhoy & Shapiro 2006) of the radiative transfer with resonant scattering
of Lyα photons by hydrogen gas. The time-independent solution is acceptable
only if the time-independent state is approached in a scale shorter than that
of the evolution of the 21 cm signal regions. Unfortunately, this assumption is
not obvious. Time-dependent solutions are necessary.
2
The W-F coupling is due to the resonant scattering of Lyα photons by hydro-
gen atoms. It is described by a Boltzmann-like integro-differential equation of
radiative distribution in the phase space. The W-F coupling is sensitive to the
evolution of the photon distribution in the frequency space. The analytical
solutions of this equation with resr background has revealed that the fre-
quency distribution of photons under the resonant scattering of Lyα photons
by gaseous hydrogen atoms without recoil has two features: 1. an extremely
flat top of the photon profiles in the frequency space; 2. a sharp boundary
of the flat top range(Field 1958). Currently, the numerical results are still far
from precise to match these features (e.g. Meiksin, 2006).
In this paper, we develop a numerical solver with the weighted essentially
nonoscillatory (WENO) scheme. The WENO method has high order of ac-
curacy and good convergence in capturing discontinuities as well as to be
significantly superior over piecewise smooth solutions containing discontinu-
ities (Shu 2003). WENO schemes have been widely used in applications. It
is also effective in solving Boltzmann equations (Carrillo et al. 2003, 2006)
and radiative transfer (Qiu et al. 2006, 2007, 2008). Therefore, one can expect
that the integral-differential equations of resonant scattering can be properly
handled numerically by the WENO scheme.
The paper is organized as follows. Section 2 presents the basic equations of
the resonant scattering of radiation. Section 3 gives the numerical solver of
the WENO scheme. Section 4 presents the tests of the numerical solver. The
time scale of the W-F coupling is briefly addressed in Section 5. A discussion
and conclusion are given in Section 6.
2 Basic equations
2.1 Radiative transfer equations with resonant scattering
Consideringaspatiallyhomogeneousandisotropicallyexpandinginfinitemedium
consisting of neutral hydrogen with temperature T, the kinetics of photons in
the frequency space is described by the radiative transfer equation with reso-
nant scattering (Hummer & Rybicki, 1992; Rybicki & Dell’antonio 1994)
∂J(x,t) cH ∂J(x,t)
+2HJ(x,t) =
∂t − v ∂x
T
kcφ(x)J(x,t)+kc (x,x′)J(x′,t)dx′ +C(t)φ(x) (1)
− R
Z
3
where J is the specific intensity in terms of the photon number, H(t) =
R˙(t)/R(t)istheHubbleparameter,R(t)isthecosmicfactor,v = (2k T/m)1/2
T B
is the thermal velocity of hydrogen atom, the dimensionless frequency x is re-
lated to the frequency ν and the resonant frequency ν by x = (ν ν )/∆ν ,
0 0 D
−
and ∆ν = ν v /c is the Doppler broadening. The parameter k = χ/∆ν ,
D 0 T D
and the intensity of the resonant absorption χ is given by χ = πe2n f /m c,
1 12 e
where n the number density of neutral hydrogen HI at ground state, and
1
f = 0.416 is the oscillation strength. The cross section at the line center is
12
πe2
σ = f (∆ν )−1. (2)
0 12 D
m c
e
In eq.(1), C(t) is the source of photons with the frequency distribution φ(x),
which is the Voigt function of the frequency profile with the center at ν , i.e.
0
a ∞ e−y2
φ(x) = dy (3)
π3/2 (x y)2+a2
−Z∞ −
where a is the ratio of the natural to the Doppler broadening. We have a =
A /(8π∆ν ), and A = 6.25 108 Hz is the Einstein spontaneous emission
21 D 21
×
coefficient. φ(x) is normalized with φ(x′)dx′ = 1. When a 0, we have
→
pure Doppler broadening as Z
1
φ (x) = e−x2. (4)
D
√π
Theredistributionfunction (x,x′)givestheprobabilityofaphotonabsorbed
R
atthefrequencyx′,andre-emittedatthefrequencyx.Itdependsonthedetails
of the scattering (Henyey 1941; Hummer 1962; Hummer, 1969). If we consider
coherent scattering without recoil, it is
(x,x′) = (5)
R
∞
1 x +u x u
e−u2 tan−1 min tan−1 max − du
π3/2 a − a
|x−Zx′|/2 (cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21)
∞
where x = min(x,x′) and x = max(x,x′). Obviously, (x,x′)dx′ =
min max
R
Z
−∞
φ(x). In the case of a = 0, i.e. considering only the Doppler broadening, eq.(5)
becomes
1
(x,x′) = erfc[max( x , x′ )]. (6)
R 2 | | | |
4
Considering the recoil of atoms, the redistribution function for the Doppler
broadening is (Field, 1959, Basko, 1981)
1
(x,x′) = e2bx′+b2erfc[max( x+b , x′ +b )], (7)
R 2 | | | |
where parameter b = hν /mv c = 2.5 10−4(104/T)1/2.
0 T
×
2.2 Re-scaling the equations
We use the new time variable τ defined as τ = cn σ t, which is in units
1 0
of the mean free flight time of photons at resonant frequency. The number
density of neutral hydrogen atoms n = f n , with f being the fraction of
1 HI H HI
neutral hydrogen. For the concordance ΛCDM model, we have n = 1.88
H
×
10−7(Ω h2/0.02)(1+z)3 cm−3. The factor 0.75 is from hydrogen abundance.
b
Therefore,
T 1/2 10 3 0.022
t = 0.054τf−1 yrs (8)
HI 104 1+z Ω h2
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) b (cid:19)
We re-scale the eq.(1) by the following new variables
J′(x,t) = [R(t)/R ]2J, C′(t) = [R(t)/R ]2C. (9)
0 0
Thus, eq.(1) becomes
∂J′(x,τ)
= φ(x)J′(x,τ)
∂τ −
∂J′
+ (x,x′)J′(x′,τ)dx′ +γ +C′(τ)φ(x). (10)
R ∂x
Z
We will use J for J′ and C for C′ in the equations below. The parameter
γ is the so-called Sobolev parameter, γ = (H/v k) = (8πH/3A λ3n ) =
T 12 1
(Hm ν /πe2n f ), where λ is the wavelength for Lyα transition. γ is simply
e 0 1 12
related to the Gunn-Peterson optical depth τ by
GP
0.25 1/2 Ω h2 1+z 3/2
γ−1 = τ = 4.9 105h−1f b . (11)
GP HI
× (cid:18)ΩM (cid:19) 0.022!(cid:18) 10 (cid:19)
The redshift evolution of f is dependent on the reionization models. Before
HI
reionization f 1; after reionization f 10−5 in average. Therefore, the
HI HI
≃ ≃
5
parameter γ has to be in the range from 1 (z 7) to 10−7 (z 10). For static
≤ ≥
background, γ = 0.
3 Numerical solver: the WENO scheme
3.1 Computational domain and computational mesh
The computational domain in the case of static background is x [ 6,6].
∈ −
The initial condition is J(x,0) = 0. The boundary condition is
J(x,τ) = 0, at x = 6. (12)
| |
In the case of expanding background, i.e.γ = 0, the computational domain is
6
bigger than x [ 6,6] depending on the value of the Sobolev parameter γ.
∈ −
The domain (x ,x ) is chosen such that for the particular value of γ we
left right
have
J(x ,τ) 0, J(x ,τ) 0. (13)
left right
≈ ≈
For example, the domain is taken to be x [ 100,6] for the case of γ = 10−3.
∈ −
Also for different values of γ the solutions reach saturation at different time.
For example, we see in our numerical results that the solution J(ν,t) of eq.
(9) reaches saturation at time τ = 100 for γ = 10−1, and reaches saturation
at time τ = 104 for γ = 10−3. The computational domain is discretized into a
uniform mesh in the x direction,
x = x +i∆x, i = 0,1,2, ,N ,
i left x
···
where∆x = (x x )/N ,isthemeshsize.WealsodenoteJn = J(x ,τn),
right− left x i i
the approximate solution values at (x ,τn).
i
3.2 Algorithm of the spatial derivative
To calculate ∂J, we use the fifth order WENO method (Jiang & Shu, 1996).
∂x
That is,
∂J(x ,τn) 1
i ˆ ˆ
(h h ) (14)
j+1/2 j−1/2
∂x ≈ ∆x −
6
ˆ
where the numerical flux h is obtained by the procedure given below. We
j+1/2
use the upwind flux in the fifth order WENO approximation because the wind
direction is fixed (negative). First we denote
h = J(x ,τn), i = 2, 1, ,N +3 (15)
i i x
− − ···
where n is fixed. The numerical flux from the WENO procedure is obtained
by
ˆ ˆ(1) ˆ(2) ˆ(3)
h = ω h +ω h +ω h , (16)
i+1/2 1 i+1/2 2 i+1/2 3 i+1/2
ˆ(m)
where h are the three third order fluxes on three different stencils given
i+1/2
by
1 5 1
ˆ(1)
h = h + h + h ,
i+1/2 −6 i−1 6 i 3 i+1
1 5 1
ˆ(2)
h = h + h h ,
i+1/2 3 i 6 i+1 − 6 i+2
11 7 1
ˆ(3)
h = h h + h ,
i+1/2 6 i+1 − 6 i+2 3 i+3
and the nonlinear weights ω are given by,
m
ωˇ γ
m l
ω = , ωˇ = , (17)
m 3 l (ǫ+β )2
l
ωˇ
l
l=1
X
where ǫ is a parameter to avoid the denominator to become zero and is taken
as ǫ = 10−8. The linear weights γ are given by
l
3 3 1
γ = , γ = , γ = , (18)
1 2 3
10 5 10
and the smoothness indicators β are given by
l
13 1
β = (h 2h +h )2 + (h 4h +3h )2,
1 i−1 i i+1 i−1 i i+1
12 − 4 −
13 1
β = (h 2h +h )2 + (h h )2,
2 i i+1 i+2 i i+2
12 − 4 −
13 1
β = (h 2h +h )2 + (3h 4h +h )2.
3 i+1 i+2 i+3 i+1 i+2 i+3
12 − 4 −
7
3.3 High order numerical integration
The integration of the resonance scattering term is calculated by a fifth order
quadrature formula (Shen et al. 2007)
xright Nx
f(x)dx = ∆x w f(x )+O(∆x5), (19)
j j
xlZeft jX=0
where the weights are defined as,
251 299 211 739
w = , w = , w = , w = ,
0 1 2 3
720 240 240 720
739 211 299 251
w = , w = , w = ,w = ,
Nx−3 720 Nx−2 240 Nx−1 240 Nx 720
and w = 1 otherwise. Notice that this numerical integration is very costly
j
and uses most part of the CPU time. For the equation with recoil and the
redistribution function with a = 0, we have used a grouping of the numerical
integration operations at different x locations so that the computational cost
can be reduced to order O(N) rather than O(N2), where N is the number of
grid points in x, without changing mathematically the algorithm and its accu-
racy. The O(N) algorithm for numerical integration is given in the appendix
which further highlights the speed and accuracy of the numerical algorithm
proposed in this paper. Unfortunately, this grouping technique does not work
for the case a = 0, hence the CPU cost for the case with a = 0 is much larger.
6 6
3.4 Time evolution
To evolve in time, we use the third-order TVD Runge-Kutta time discretiza-
tion (Shu & Osher, 1988). For systems of ODEs u = L(u), the third order
t
Runge-Kutta method is
u(1)=un +∆τL(un,τn),
3 1
u(2)= un + (u(1) +∆τL(u(1),τn +∆τ)),
4 4
1 2 1
u(3)= un + (u(2) +∆τL(u(2),τn + ∆τ)).
3 3 2
8
4 Tests of the WENO solver
4.1 Static background
We first test the WENO solver with two analytical solutions of eq. (10) with
static background, H = γ = 0 and Doppler broadening a = 0 (Field, 1958).
The first analytical solution is for the initial radiative field J(x,τ) = 0 and
the constant source C(τ) = 1. It is
J(x,τ) = π−1/2[1 exp( τe−x2)] (20)
− −
∞
+ ew2[1 (1+τe−w2)exp( τe−w2)]erf(w)dw.
− −
Zx
The second analytical solution of eq.(10) is also for H = γ = 0, a = 0, but the
source C = 0, while the initial radiative field is
J(ν,0) = π−1/2e−ν2. (21)
The solution is
∞
J(x,τ) = π−1/2e−x2exp( τe−x2)+τ e−w2exp( τe−w2)erf(w)dw. (22)
− −
Zx
The analytical solutions (20) and (22) are shown in Figures 1 and 2 respec-
tively. We also plot the numerical solutions given by our algorithm in same
figures. The numerical results show very small deviation from the analytical
solutions.
Acommon featureofFigures 1 and2is that theoriginallyDoppler peak atthe
center of the frequency profile gradually becomes a flat plateau. The width
of the plateau increases with time. This is because the resonant scattering
makes a non-uniform distribution in the frequency space to a uniform one. It
issimilar tothatinthephysical space, diffusiongenerally leads toanevolution
from non-uniform distribution to a uniform one. The height of the plateau of
Figure 1 is increasing with time, while in Figure 2 it is decreasing, because for
the case C = 1, the number of photons increases, while for the case of C = 0,
the total number of photons is conserved.
9
τ=10000
103
τ=1000
102
τ=100
101
τ=10
x)
J(100 τ=1
10-1 τ=0.1
10-2 τ=0.01
10-3
-4 -2 0 2 4 6
x
Fig. 1. Static solutions (γ = 0) for pureDoppler redistribution (a = 0) of eq.(10), in
which C = 1 and J(x,0) = 0. The analytical solutions are shown by dashed lines,
while the numerical results are shown by solid lines.
τ=0.01,0.1
τ=1
τ=10
τ=100
τ=1000
τ=10000
10-1
x)
J(
10-2
10-3
-4 -2 0 2 4
x
Fig. 2. Static solutions (γ = 0) for pure Doppler redistribution (a = 0) of eq.(10),
in which C = 0 and J(x,0) = π−1/2exp( x2). The analytical solutions are shown
−
by dashed line, while the numerical results are shown by solid lines.
4.2 Expanding background
The second test is given by considering expanding background. Without ab-
sorption and scattering, eq.(10) becomes
∂J ∂J
= γ +C(τ)φ(x). (23)
∂τ ∂x
10
