EFFECT OF DUST ON Lyα PHOTON TRANSFER IN AN OPTICALLY THICK HALO

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

$$\begin{equation} \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}, \end{equation} \tag{ 1 }$$

where ϕ(x, a) is the normalized Voigt profile (Hummer 1965). As usual, the photon frequency ν in Equation (1) is described by the dimensionless frequency x ≡ (ν − ν₀)/Δν_D, with ν₀ = 2.46 × 10¹⁵ s⁻¹ being the resonant frequency, Δν_D = ν₀(v_T/c) = 1.06 × 10¹¹(T/10⁴)^1/2 Hz the Doppler broadening, $v_T=\sqrt{2k_BT/m}$ the thermal velocity, and T the gas temperature of the halo. σ₀/π^1/2 is the cross section of scattering at the resonant frequency ν₀. The parameter a in Equation (1) is the ratio of natural to Doppler broadening. For the Lyα line, a = 4.7 × 10⁻⁴(T/10⁴)^−1/2. The optical depth of Lyα photons with respect to H i resonant scattering is $\tau _s(x) = n_{{\rm H\,\mathsc{i}}}R \sigma (x) = \tau _0\phi (x,a)$, where $\tau _0=n_{{\rm H\,\mathsc{i}}}\sigma _0R$.

If the absorption and scattering of dust are described by the effective cross section per hydrogen atom σ_d(x), the total optical depth is given by

$$\begin{equation} \tau (x) = \tau _0\phi (x,a)+ \tau _{d}(x), \end{equation} \tag{ 2 }$$

where the dust optical depth $\tau _{d}(x)=n_{{\rm H\,\mathsc{i}}}\sigma _{d}(x)R$. This is equal to assuming that dust is uniformly distributed in the intergalactic medium. 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α photons in a spherical halo with dust is given by

$$\begin{eqnarray} &&\lefteqn{ {\partial I\over \partial \eta } + \mu \frac{\partial I}{\partial r}+\frac{(1-\mu ^2)}{r}\frac{\partial I}{\partial \mu } - \gamma \frac{\partial I}{\partial x} = - \phi (x;a)I } \nonumber \\ & & \qquad + \int \mathcal {R}(x,x^{\prime };a)I(\eta, r, x^{\prime },\mu ^{\prime })dx^{\prime }d\mu ^{\prime }/2-\kappa (x) I \nonumber\\ &&\qquad + 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,\nonumber\\ \end{eqnarray} \tag{ 3 }$$

where I(t, r_p, x, μ) is the specific intensity, which is a function of time t, radial coordinate r_p, frequency x, and the direction angle, μ = cos θ, with respect to the radial vector r.

In Equation (3), we use the dimensionless time η defined as $\eta =cn_{{\rm H\,\mathsc{i}}}\sigma _0 t$ and the dimensionless radial coordinate r defined as $r=n_{{\rm H\,\mathsc{i}}}\sigma _0 r_p$. That is, η and r are, respectively, in the units of mean-free flight time and mean-free path of photon ν₀ 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 = η + const. With the dimensionless variable, the size of the halo R is equal to τ₀.

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

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

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

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

The redistribution 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α 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 Equations (4) or (5). The effect of recoil actually is under control (Roy et al. 2009c, 2010). We will address this in the next section.

The absorption and scattering of dust are described by the term κ(x)I of Equation (3), where κ(x) = σ_d/σ₀, which is of the order of 10⁻⁸(T/10⁴)^1/2 (Draine & Lee 1984; Draine 2003). The term with A of Equation (3) describes albedo, i.e., A ≡ σ_s/σ_d, where σ_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 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', and is simply given by a phase function as

$$\begin{eqnarray} \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}}\nonumber\\ &=& \sum _{l=0}^{\infty }\frac{(2 l+1)}{2} g^l P_l(\mu)P_l(\mu ^{\prime }), \end{eqnarray} \tag{ 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 Equation (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α photon, we will take g = 0.73 for realistic dust scattering (Li & Draine 2001). The integral of Equation (6) is performed in Appendix A.

In Equation (3), the term with the parameter γ is due to the expansion of the universe. If n_H is equal to the mean of the number density of cosmic hydrogen, we have γ = τ⁻¹_GP, and τ_GP is the Gunn–Peterson optical depth. Since the Gunn–Peterson optical depth is of the order of 10⁶ at high redshift (e.g., Roy et al. 2009c), the parameter γ is of the order of 10⁻⁵ to 10⁻⁶. Therefore, if the optical depth of halos is equal to or less than 10⁶, the term with γ of Equation (3) can be ignored.

In Equation (3), we neglect the effect of collision transition from the H(2p) state to the H(2s) state, which can significantly affect the escape of Lyα photons when the H i column density is higher than 10²¹ cm⁻² and the dust absorption is very small (Neufeld 1990). This is generally 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).

Equation (6) indicates that the transfer equation (3) can be solved with the Legendre expansion I(η, r, x, μ) = ∑_lI_l(η, r, x)P_l(μ). If we take only the first two terms, l = 0 and 1, it is the Eddington approximation as

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

where

$$\begin{eqnarray} J(\eta,r,x) &=& \frac{1}{2}\int _{-1}^{+1}I(\eta,r,x,\mu)d\mu,\nonumber\\ &&[-6pt] \\ [-6pt] F(\eta,r,x) &=& \frac{1}{2}\int _{-1}^{+1}\mu I(\eta,r,x,\mu)d\mu.\nonumber \end{eqnarray} \tag{ 8 }$$

They are, respectively, the angularly averaged specific intensity and flux. Defining j = r²J and f = r²F, Equation (3) yields the equations of j and f as

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

$$\begin{eqnarray} \frac{\partial f}{\partial \eta } + \frac{1}{3} \frac{\partial j}{\partial r} - \frac{2}{3}\frac{j}{r} & = & -(1-Ag)\kappa f +\gamma \frac{\partial f}{\partial x}-\phi (x;a)f.\qquad \end{eqnarray} \tag{ 10 }$$

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

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

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

where S₀ and ϕ_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 a given S₀ = S are equal to Sj₁(x) and Sf₁(x), where j₁(x) and f₁(x) are the solutions of S₀ = 1. On the other hand, Equation (3) is not linear with respect to the function ϕ_s(x). The solution f(x) for a given ϕ_s(x) is not equal to ϕ_s(x)f₁(x), where f₁(x) is the solution of ϕ_s(x) = 1.

In the range outside the halo, r > R, no photons propagate in the direction μ < 0. The boundary condition at r = R given by ∫⁻¹₀μJ(η, R, x, μ)dμ = 0 is then (Unno 1955)

$$\begin{equation} j(\eta, R,x)=2f(\eta, R,x). \end{equation} \tag{ 12 }$$

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

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

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

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 κ(x) is on the order of 10⁻⁸, its effect will be significant only for halos with optical depth τ₀ ⩾ 10⁶, and ignorable for τ₀ ⩽ 10⁵. To illustrate the dust effect, we use halos of R = τ₀ ⩽ 10⁴ and take larger κ to be ≃ 10⁻⁴ to 10⁻². We also assume that κ is frequency-independent. We consider below only the case of gray dust, i.e., κ is independent of the frequency x. This is certainly not realistic dust. However, the frequency given in the solutions below mostly is in the range |x| < 4. Therefore, the approximation of gray dust would be proper if the cross section of dust is not significantly frequency-dependent in the range |x| < 4.

As the first example of numerical solutions, we show the W–F effect, which requires that the distribution of Lyα photons in the frequency space should be thermalized near the resonant frequency ν₀. The W–F effect illustrates the difference between the analytical solutions of the F–P approximation and those of Equations (9) and (10). The former cannot 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 (Equation (3)).

Figure 1 presents a solution of mean intensity j(η, r, x) at time η = 500 and radial coordinate r = 10² for a halo with size R ≫ r = 10². 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₀ = 1. The solutions of Figure 1 are actually independent of R, if R ≫ 10². 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.

Zoom In Zoom Out Reset image size

A remarkable feature shown in Figure 1 is that all j(η, r, x) have a flat plateau in the range |x| ⩽ 2. This gives the frequency range of the W–F local thermalization (Roy et al. 2009b, 2009c). The range of the flat plateau |x| ⩽ 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 κ term of Equation (3) changes the frequency distribution of photons. The redistribution function (Equation (6)) also does not change the frequency distribution of photons. This point can also be seen from Equations (9) and (10), in which the κ 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| ⩽ 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₀ϕ_s(x) with arbitrary ϕ_s(x) that is non-zero within |x| < 2 (Roy et al. 2009b, 2009c). This property can also be used as a test of the simulation code. It requires that the simulation results of the flat plateau should hold, regardless of whether the source is monochromatic or has finite width around ν₀.

To study the effects of dust scattering on the Lyα photon escape, we show in Figure 2 the flux f(η, r, x) of Lyα photons emergent from halos at the boundary r = R = 10² for Model I. The three panels of Figure 2 correspond to κ = 10⁻⁴, 10⁻³, and 10⁻² from left to right, respectively. The source starts to emit photons at η = 0 with a stable luminosity S₀ = 1 and with a Gaussian profile $\phi _s(x)=(1/\sqrt{\pi })e^{-x^2}$.

Zoom In Zoom Out Reset image size

Figure 2 clearly shows that the time evolution of f(η, r, x) is κ-independent. Although the cross section of dust scattering increases by about 100 times from κ = 10⁻⁴ to κ = 10⁻², the curves of the left and right panels in Figure 2 are actually almost identical.

According to the scenario of "single longest excursion," the photon escape is not a process of Brownian random walk in the spatial space, but rather a transfer in the frequency space (Osterbrock 1962; Avery & House 1968; Adams 1972, 1975; Harrington 1973; Bonilha et al. 1979). A 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 equation (Equation (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 ≠ 0 with respect to the direction of the incident particle, the angular distribution will keep isotropic undergoing a g ≠ 0 scattering. Hence, dust scattering also has no effect on the angular distribution.

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 η = 2000 and 3000 are the same. This means the flux f(η, R, x) at the boundary R is already stable, or saturated at the time η ⩾ 2000. The small difference between the curves of η = 1000 and η ⩾ 2000 of the left panel indicates that the flux is still not yet completely saturated at the time η = 1000. However, comparing the middle and right panels of Figure 3, we see that for κ = 10⁻³, the flux has already saturated at η = 1600, while it has saturated at η = 800 for κ = 10⁻². That is, the stronger the dust absorption, the shorter the saturation timescale. The timescales of escape or saturation do not increase by dust absorption, and even decrease with respect to the medium without dust. Stronger absorption leads to a shorter timescale of saturation.

Zoom In Zoom Out Reset image size

Obviously, dust absorption does not help in producing photons for the "single longest excursion." Therefore, dust absorption cannot make the timescale of producing photons for the "single longest excursion" smaller. However, dust absorptions are effective in reducing the number of photons trapped in the state of local thermalized equilibrium |x| < 2 (see also Section 4.2). This leads to the fact that the higher the value of κ, the shorter the timescale of saturation.

Since Lyα photons underwent a large number of resonant scattering before escaping from the halo with optical depth τ₀ ≫ 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α photons in R = τ₀ ≫ 1 halos.

To calculate the effective optical depth, we first give the total flux of Lyα photons emergent from the halo of radius R, which is defined as F(η) = ∫f(η, R, x)dx. Figure 4 plots F(η) as a function of time η for a halo with sizes R = τ₀ = 10² and 10⁴. 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², and the lower panels are of R = 10⁴. In each panel of R = 10², we have three curves corresponding to κ = 10⁻⁴, 10⁻³, and 10⁻², respectively. In cases of R = 10⁴, we take κ = 10⁻⁴ and 10⁻³.

Zoom In Zoom Out Reset image size

The left panel of Figure 4 shows that the three curves of κ = 10⁻⁴, 10⁻³, and 10⁻² are almost the same. This is consistent with Figure 2 for model I in that the time evolution of f is κ-independent for the pure scattering dust. For the pure absorption dust (the right panel of Figure 4), the saturated flux is smaller for larger κ. We can also see from Figure 4 that the timescale of approaching saturation is smaller for larger κ. 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 τ_effect ≡ −(1/κ)ln f_S. The results are shown in Table 1, in which τ_a is the dust absorption depth. It is interesting to see that the effective absorption optical depth is always equal to about a few times of the optical depth of resonant scattering τ₀, regardless of whether τ_a is less than 1. Namely, the effective absorption depth τ_effect of dust is roughly proportional to τ₀.

According to the random walk scenario, if a medium has optical depths of absorption τ_a and scattering τ_s, the effective absorption optical depth should be equal to $\tau _{\rm effect}=\sqrt{\tau _a(\tau _a+\tau _s)}$ (Rybicki & Lightman 1979, p. 393). 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 the optical depth of dust is lower than the optical depth of resonant scattering τ₀, the timescale of photon escaping basically is not affected by the dust, but is proportional to τ₀; therefore, the absorption is also proportional to τ₀.

With the total flux, we can define the escaping coefficient of the Lyα photon as f_esc(η, τ₀) ≡ F(η)/F₀, where F₀ is the flux of the center source. Figure 5 shows f_esc(η, τ₀) at three times η = 5 × 10³, 10⁴, and 3.2 × 10⁴ for Model II and κ = 10⁻³. At η = 5 × 10³, the flux of halos with τ₀ ⩽ 10³ is saturated. At η = 10⁴, halos with τ₀ ⩽ 3 × 10³ are saturated, and all halos of τ₀ ⩽ 10⁴ are saturated at η = 3.2 × 10⁴.

Zoom In Zoom Out Reset image size

A remarkable feature of Lyα photon emergent from optically thick medium is the double-peaked profile. Figures 1–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 τ₀ larger than 10². Figure 6 presents the double peak frequency |x_±| as a function of aτ₀, where the parameter a is taken to be 10⁻² (left) and 5 × 10⁻³ (right). Comparing the curves with and without dust in Figure 6, we can say that the dust effect on |x_±| is very small until aτ₀ = aR = 10².

Zoom In Zoom Out Reset image size

In the range aτ₀ < 20, the |x_±|–τ₀ relation is almost flat with |x_±| ≃ 2. This 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₋, are basically 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 is weakly dependent on τ₀. Thus, we always have x_± ≃ ±2. When the optical depth is larger, aτ₀ ∼ 10², 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_±| will increase with τ₀.

Figure 6 also shows a line x_± = ±(aτ₀)^1/3, which is given by the analytical solution of the F–P approximation, in which the Doppler broadening core in the Voigt profile is ignored (Harrington 1973; Neufeld 1990; Dijkstra et al. 2006). The numerical solutions of Equations (3) or (9) and (10) deviate from the (aτ₀)^1/3-law at all parameter range of Figure 6. The deviation at aτ₀ < 20 is due to the fact that the Doppler broadening core in the Voigt profile is ignored in the F–P approximation; no locally thermal equilibrium can be reached. Therefore, in the range aτ₀ < 20, |x_±| of the WENO solution is larger than the (aτ₀)^1/3-law. In the range of aτ₀ > 20, the F–P approximation yields a faster diffusion of photons in the frequency space. This point can be seen in the comparison between an F–P solution with Field's analytical solution (Figure 1 in Rybicki & Dell'Antonio 1994). In this range, the numerical results of |x_±| are less than the (aτ₀)^1/3-law.

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–3, neither the narrowing nor the widening of the two peaks.

If dust absorption can cause narrowing, the absorption should be weaker at |x| ∼ 0 and stronger at |x| ⩾ 2. Similarly, if dust absorption can cause widening, the absorption should be weaker at |x| ∼ 2 and stronger at |x| ∼ 0. To test these assumptions, Figure 7 plots ln [f(η, r, x, κ = 0)/f(η, r, x, κ)] as a function of x. It measures the x(frequency) dependence of the flux ratio with and without dust absorption. We take large η, 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 in |x| > 2; therefore, the assumption of the narrowing is ruled out. Figure 7 also shows that the curves of ln [f(η, r, x, κ = 0)/f(η, r, x, κ = 10⁻³)] are almost flat in the range |x| < 2. Therefore, the assumption of widening of the two peaks can also be ruled out.

Zoom In Zoom Out Reset image size

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

If the radiation from the sources has a continuum spectrum, the effect of neutral hydrogen halos is to produce an absorption line at ν = ν₀. The profile of the absorption line can also be found by solving Equations (9) and (10), replacing the boundary equation (Equation (11)) by

$$\begin{equation} j(\eta, 0, x)=0, \quad f(\eta, 0, x)=S_0. \end{equation} \tag{ 14 }$$

That is, we assume that the original spectrum is flat in the frequency space. The spectrum of the flux emergent from the halo of R = 10² and 10⁴ with the central source of Equation (14) for dust models I, II, and III is shown in Figure 8. All curves are for large η, i.e., they are saturated.

Zoom In Zoom Out Reset image size

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

The flux profiles of Figure 8 are typically 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 that at the flat wing. This is because more photons are stored in the frequency range |x| < 2. According to the redistribution function equation (Equation (4)), the probability of transferring an x' photon to an |x| < |x'| photon is greater than that from |x'| to |x| > |x'|. 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 at |x| < 2 are thermalized; 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. This puts the shoulder at |x| ∼ 2.

As expected, for model I (left panels of Figure 8), the double profile is completely κ-independent. Dusty scattering does not change the flux and its profile. For models II and III, the higher the κ, 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 κ-independent. This once again shows that dust absorption and scattering cause neither narrowing nor widening of the double-peaked profile. However, for higher κ 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.

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

We first give the WENO reconstruction procedure in approximating ∂j/∂x,

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

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

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

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

$$\begin{displaymath} \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)}, \end{displaymath}$$

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

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

The nonlinear weights ω_m are given by

$$\begin{displaymath} \omega _m=\frac{\breve{\omega }_m}{\sum _{l=1}^3\breve{\omega }_l}, \quad \breve{\omega }_l=\frac{\gamma _l}{(\epsilon +\beta _l)^2}, \end{displaymath}$$

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

$$\begin{displaymath} \gamma _1=\frac{3}{10},\quad \gamma _2=\frac{3}{5},\quad \gamma _3=\frac{1}{10}, \end{displaymath}$$

and the smoothness indicators β_l are given by

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

To approximate the r-derivatives in the system of Equations (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

$$\begin{displaymath} {\bf u}_t+A{\bf u}_r, \end{displaymath}$$

where u = (j, f)^T and

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

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

To implement the boundary condition (Equation (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 v = R⁻¹u. Denote v = (v₁, v₂)^T; now u_t + Au_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}; \quad \frac{\partial v_2}{\partial t}+\lambda _2\frac{\partial v_2}{\partial r},$$

where $\lambda _1={\sqrt{3}}/{3}$ and $\lambda _2=-{\sqrt{3}}/{3}$. The characteristic line starting from the boundary r = r_max for the first equation is pointing outside the computational domain while the one for the second equation it 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 α and β. We can calculate the values of α and β 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 β = 0. We use extrapolation to obtain the value of v₁ and then compute the value v₂. In the end, we transfer v back to the physical space by u = Rv.

To evolve in time, we use the third-order total variation diminishing Runge–Kutta time discretization (Shu & Osher 1988). For the system of ordinary differential equations u_t = L(u), the third-order Runge–Kutta method is

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