RESONANT SCATTERING AND Lyα RADIATION EMERGENT FROM NEUTRAL HYDROGEN HALOS

The property of the halo around individual luminous object depends on the luminosity, the spectrum of UV photon emission, and the time evolution of the center object. For luminous objects at high redshift, the halos generally consist of three spheres (Cen 2006; Liu et al. 2007). The innermost region is the highly ionized Strömgren sphere, or the H ii region, which is optically thin of Lyα photons. The second region, which is just outside the H ii region, is optically thick of Lyα photons. The temperature of the baryon gas is about 10⁴ K, which is due to the heating of the UV photons. The third region is outside the heated region. It is unheated and therefore the temperature of the baryon gas can be as low as 10² K.

In this context, we will study the transfer of Lyα photons in a radius R spherical halo of neutral hydrogen with temperature T in the range of 10²–10⁴ K. Assuming that the uniformly distributed H i gas has number density $n_{\rm H\,{\mathsc i}}$, the optical depth over a light path dl is $d\tau (\nu)=\sigma (\nu)n_{\rm H\,{\mathsc i}}{\it dl}$, where σ(ν) is the cross section of the resonant scattering of Lyα photons by hydrogen, given as

$$\begin{equation} \sigma (x)=\sigma _0 \phi (x,a), \end{equation} \tag{ 1 }$$

where x is the dimensionless frequency defined by x ≡ (ν − ν₀)/Δν_D, ν₀ being the resonant frequency. Δν_D = ν₀(v_T/c) is the Doppler broadening and $v_T=\sqrt{2k_BT/m}$. Therefore, x measures the frequency deviation Δν = |ν − ν₀| in units of Δν_D = 1.06 × 10¹¹(T/10⁴)^1/2 Hz. σ₀ = πe²f/m_ecΔν_D = 1.10 × 10⁻² cm² is the cross section of the resonant scattering at the frequency ν₀ = 2.46 × 10¹⁵ s⁻¹. The function ϕ(x, a) in Equation (1) is the normalized profile given by the Voigt function as (Hummer 1965)

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

The parameter a in Equation (2) is the ratio of the natural to the Doppler broadening. For the Lyα line, a = 4.7 × 10⁻⁴(T/10⁴)^−1/2. The optical depth of the halo with column number density of neutral hydrogen $N_{\rm H\,{\sc i}}=n_{\rm H\,{\sc i}}R$ is

$$\begin{equation} \tau (x) = N_{\rm H\,{\sc i}} \sigma (x) = \tau _0\phi (x,a). \end{equation} \tag{ 3 }$$

Since $\phi (0,a)=1/\sqrt{\pi }$ when a ≪ 1, the line-center optical depth is then $\tau (0)=\tau _0/\sqrt{\pi }$ and

$$\begin{equation} \tau _0 = 1.04\times 10^{7}\left(\frac{T}{10^4}\right)^{-1/2} \left(\frac{N_{\rm H\,{\sc i}}}{10^{20}\,{\rm cm^{2}}}\right). \end{equation} \tag{ 4 }$$

The radiative transfer of Lyα photons in spherical halo is described by the equation of the specific intensity I(η, r, x, μ) as

$$\begin{eqnarray} {\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 \\ &+& \int \mathcal {R}(x,x^{\prime };a)I(\eta, r, x^{\prime },\mu ^{\prime })dx^{\prime }d\mu ^{\prime } + S, \end{eqnarray} \tag{ 5 }$$

where we use dimensionless time η defined as ${\eta =cn_{\rm H\, \mathsc i}\sigma _0 t}$ and dimensionless coordinate r defined as ${r=n_{\rm H\,\mathsc i}\sigma _0 r_p}$, with r_p being the physical radial coordinate. That is, η and r are, respectively, in the units of mean free flight time and mean free path of photon ν₀. For a signal propagated in the radial direction with the speed of light, we have r = η + const. In Equation (5), μ = cos θ is the direction relative to the radial vector $\hbox{\bf r}$.

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 1941; Hummer 1962, 1969). If we consider coherent scattering without recoil, the redistribution function with the Voigt profile (Equation (2)) is

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

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{ 7 }$$

The redistribution function of Equation (7) 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 (5). That is, the destruction processes of Lyα photons, such as the two-photon process (Spitzer & Greenstein 1951; Osterbrock 1962), are ignored in Equation (5). In Equations (6) and (7), we also do not consider the recoil of atoms. It is equal to assuming that the mass of atom is very large. The effect of recoil actually is under control (Roy et al. 2009c). We will address it in next section.

In Equation (5), 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 $\gamma =\tau _{\rm {GP}}^{-1}$, and $\tau _{\rm {GP}}$ is the Gunn–Peterson optical depth. Since Gunn–Peterson optical depth is of the order of 10⁶ at high redshift (e.g., Roy et al. 2009c), we will take the parameter γ = 10⁻⁵–10⁻⁶.

When the optical depth is large, we can take 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{ 8 }$$

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

$$\begin{eqnarray} {\partial j\over \partial \eta } + \frac{\partial f}{\partial r} & = & - \phi (x;a)j + \int \mathcal {R}(x,x^{\prime };a)j dx^{\prime }+ \gamma \frac{\partial j}{\partial x}+ r^2S,\nonumber \\ \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} & = & - \phi (x;a)f. \end{eqnarray} \tag{ 10 }$$

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

The source term S in the Equations (9) and (10) can be described by a boundary condition of j and f at r = r₀. We can take r₀ = 0, as r₀ is not important if the optical depth of the halo is large. Thus, we have

$$\begin{equation} j(\eta, 0, x)=0, \hspace{1cm} 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 Equations (9)–(11) are linear, the intensity S₀ can be taken as any constant. That is, the solution f(x) of S₀ = S is equal to Sf₁(x), where f₁(x) is the solution of S = 1. On the other hand, Equations (9) and (10) are not linear with respect to the function ϕ_s(x), i.e., the solution f(x) of ϕ_s(x) is not equal to ϕ_s(x)f₁(x), where f₁ is the solution of ϕ_s(x) = 1.

In the range outside the halo, r>R, no photons propagate in the direction μ < 0. Therefore, 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 solve Equations (9) and (10) with the numerical method developed recently (Roy et al. 2009a, 2009b, 2009c). Some details of this method are given in the Appendix. We first solve the problems when the sources S is steady.

First we consider steady sources. That is, the parameter S₀ in Equation (11) is time independent after the switch-on of the sources at η = 0 (Equation (13)). A typical solution of the flux f(η, r, x) given by Equations (9) and (10) is shown in Figure 1, for which the source is taken to be S₀ = 1 and $\phi _s(x)=(1/\sqrt{\pi })e^{-x^2}$, i.e., the emission line width is equal to the Doppler broadening. The parameters a and γ are taken to be 10⁻³ and 10⁻⁵, respectively. The effect of γ = 10⁻⁵ actually is ignorable in these solutions. The left panel is the solutions f(η, r, x) at radius r = 10² and time η = 500, 1000, 2000, and 3000 with the boundary condition (Equations (11) and (12)). The solutions approach a stable state at time η ⩾ 2000.

Zoom In Zoom Out Reset image size

The right panel of Figure 1 is the total flux F_t(η) ≡ ∫f(η, r, x)dx. It again shows that total flux approaches a stable state with F_t = 1 at time η ⩾ 2000. If the Lyα photon transfer is due only to spatial diffusion, the time scale of a spatial transfer with an optical depth τ₀ = 10² should be as large as η ∼ τ²₀ = 10⁴. However, Figure 1 shows that photons have already escaped from the τ₀ = 10² hole within the time η ∼ 2000, which is much less than the time scale of a purely Brownian diffusion. Therefore, the transfer of photons should not be a process of purely Brownian diffusion. This point is very well known in early studies on the escape of resonant photons from opaque clouds (Osterbrock 1962; Harrington 1973; Avery & House 1968). The time scale of escape shown in Figure 1 is also consistent with the estimation given by Monte Carlo simulations (Adams, 1975; Bonilha et al. 1979). However, these works are mainly based on the time scale of the escape of photons with frequency, at which the photon can take a "single longest excursion" (Adams 1972). Figure 1 shows that the escape time scale η = 2000 is available not only for photons that can take a "single longest excursion" but also for all photons with frequency around ν₀ or x = 0.

The right panel of Figure 1 also plays the role of testing our algorithm. Because ∫ϕ(x, a)jdx = ∫R(x, x'; a)jdxdx', Equation (9) yields

$$\begin{equation} F_t(\eta, r)=F_t(\eta, 0)=S_0 \end{equation} \tag{ 14 }$$

when the solution approached the stable state. Equation (14) is the conservation of photon number. Since numerical errors which have accumulated over a long time evolution could be huge, Equation (14) is useful to check the reliability of the code. The right panel of Figure 1 shows perfectly F_t(η, r = 10²) = F_t(η, 0) = S₀ = 1 at stable state.

Figure 2 shows a solution of the mean intensity j(η, r, x) of Equations (9) and (10) at r = 10² and time η = 200, 300, and 500. The source is the same as Figure 1, namely S₀ = 1 and $\phi _s(x)=(1/\sqrt{\pi })e^{-x^2}$. Other parameters are also the same as Figure 1.

Zoom In Zoom Out Reset image size

A remarkable feature of the solutions is to show a flat plateau around x = 0. As has been shown by Roy et al. (2009c), the flat plateau actually is the Boltzmann statistical equilibrium distribution around x = 0 when the atomic mass is infinite. If the mass is finite, i.e., considering the recoil in the redistribution functions (6) or (7), the flat plateau will become e^−2bx, where b = hν₀/mv_Tc. This is the local Boltzmann distribution required by the Wouthuysen–Field effect (Roy et al. 2009b).

Without resonant scattering, the pure absorption will lead to the flux of ν₀ photons at r = 10² to be attenuated by a factor e^τ(0) ≃ 10⁻²⁵. Therefore, Figure 2 shows that a major effect of resonant scattering is to restore ν₀ photons in optically thick halos. According to the redistribution function $\mathcal {R}(x,x^{\prime })$ (Equation (5)), the probability of transferring a x' photon to a |x| < |x'| photon is larger than that from x' to |x|>|x'|. Therefore, the net effect of resonant scattering is to bring photons with frequency x' ≠ 0 to the central range x ∼ 0 of frequency space. Photons of x ∼ 0 are then effectively restored. Moreover, the restored photons are thermalized. We see from Figure 2 that in the time range from η = 200 to 500, the mean intensity j at x = 0 quickly increases by a factor larger than 10. At the same time, the flat plateau or the local thermalization around x = 0 is perfectly held. Therefore, the time scales t_ther of ν₀ photon restoring and thermalization should be less than η = 200, which is much smaller than the time scale t_escape of escaping from a halo of r = 10².

The result of t_ther ≪ t_escape in optically thick halos is very important. One can conclude that the photons of the flux f around x = 0 emergent from optically thick halos actually come from the thermalized photons, regardless of the initial distribution of the photons. The initial conditions of the photon source have already been forgotten during the thermalization. Therefore, one can expect that the profile of flux f has to be independent of the sources. We demonstrate this point with Figure 3.

Zoom In Zoom Out Reset image size

Figure 3 presents the flux f given by Equations (9) and (10) with the same parameters as the solution of Figure 1, and the source profiles are $\phi _s(x)=(1/\sqrt{\pi })e^{-x^2}$ (left panel) and $\phi _s(x)=(1/\sqrt{\pi /2})e^{-2x^2}$ (right panel). That is, the line widths are, respectively, 1 and $1/\sqrt{2}$. We see that the left and right panels of Figure 3 are almost identical within |x| ⩽ 3. Therefore, the initial distribution of photon frequency is already forgotten, and the photons of the left and right cases at η>200 actually are from the same thermalized sources. The profiles of the flux are also held if the line width is broader than 1. We will show this point when the source has a continuant spectrum (Section 3.4).

The profiles of either j or f shown in Figures 1–3 have two-peak structure. The flux f is dominated by photons with frequency x_± ≃ ±(2–3). The two-peak structure has been recognized in all the time-independent solutions of the Fokker–Planck approximation (Harrington 1973; Neufeld 1990; Dijkstra et al. 2006) and Monte Carlo simulations (Lee 1974; Zheng & Miralda-Escude 2002; Ahn et al. 2002; Cantalupo et al. 2005; Verhamme et al. 2006). A point we would like to emphasize is that this structure is independent of the profile of the source S. It is because the initial properties of the central sources have been forgotten during the local thermalization.

Since the shape of the locally thermalized distribution is time independent, the frequency of the two peaks, |x_±|, at a given r is also time independent. When the photons of the flux f mainly come from the locally thermalized photons, the frequency of the two peaks, |x_±|, should not be described by the relation |x_±| = (aτ)^1/3, because this relation is from a solution of the time-independent Fokker–Planck equation, which does not describe the thermalization (Neufeld 1990).

Figure 4 presents the peak frequency |x_±| as a function of ar for solutions given by Equations (9) and (10) with a = 10⁻². We consider only r ⩾ 10², as in the case r ⩽ 10² photons do not undergo a large enough number of scattering and therefore are not completely thermalized yet. With the dimensionless variables, ar is equal to aτ. A line |x_±| = (ar)^1/3 is also shown in Figure 4. It shows that our numerical result of |x_±| is significantly different from the (aτ₀)^1/3 law if ar < 30. That is, in the range ar < 30, photons of the flux f are dominated by the locally thermalized photons. The frequency x_± actually is dependent on the width of the flat plateau or the locally thermalized distribution of j. Therefore, x_± is always larger than two. This feature has also been addressed in Bonilha et al (1979). The curve of Figure 4 is approximately available for a = 10⁻³ and 10⁻⁴. Thus, |x_±| is larger than (ar)^1/3 if r < 3 × 10⁵ and 3 × 10⁶ for a = 10⁻³ and 10⁻⁴, respectively.

Zoom In Zoom Out Reset image size

Figure 4 shows a slow increase of |x_±| with r in the range ar ⩽ 30, and then it approaches (ar)^1/3 for larger ar. When r is large, more photons of the flux are attributed to the resonant scattering by the Lorentzian wing of the Voigt function. |x_±| is then approaching (ar)^1/3. It should be emphasized once again that all these results are independent of the intrinsic width of Lyα emission from the central source.

If the radiation from the sources has a continuant 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), but replacing the boundary Equation (11) by

$$\begin{equation} j(\eta, 0, x)=0, \hspace{1cm} f(\eta, 0, x)=S_0. \end{equation} \tag{ 15 }$$

That is, we assume that the original spectrum is flat in the frequency space.

A solution of the time evolution of j and f at r = 10² with the source Equation (15) is shown in Figure 5. The optical depths at the frequency |x|>4 are small, and therefore the Eddington approximation would 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. Therefore, the solution is still useful to study the profiles of j and f around x = 0. The profile of f typically is an absorption line with width given by the position of the two-peak structure. As expected, the profile in the range |x| < 4 is the same as the left panel of Figure 1. It shows again that the two-peak structure is independent of the frequency spectrum of the central source.

Zoom In Zoom Out Reset image size

The flux f of Figure 5 has a flat wing in both sides of x>3 and x < −3, because the effect of resonant scattering is negligible for photons with frequency |x|>4, and they can freely transfer from r = 0 to 10². On the other hand, the mean intensity j, which is the isotropic component of the intensity J (Equation (8)), does not have wings at |x|>4. That is, resonant scattering cannot store photons with |x|>4 in the halo. Nevertheless, the mean intensity j still has a flat plateau around x = 0. It means that in the period from η = 200 to η = 500, the halo trapped and stored more and more photons of |x| < 4. These photons are in the locally thermalized state.

As an application of the absorption spectrum in Figure 5, we study the profile of the red damping wing of the optically thick IGM at high redshifts. Since the variable x is independent of redshift, the profile of a red damping wing is directly given by the flux f(η, r, x) at the wing of low frequency x ⩽ 0.

If the hydrogen clouds are located far from the Lyα sources, the red damping wing can be modeled as the absorption of optically thick halos. The red damping wing of damped Lyα system (DLA) model is then given by f(x) = e^−τ(x) and x < 0, where τ(x) is from the Voigt function Equation (2) as

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

The column number density of neutral hydrogen atoms, $N_{\rm H\,{\mathsc i}}$, generally is estimated with fitting profile Equation (16) with observation, and then, $N_{\rm H\,{\sc i}}$ can be found from τ₀ by Equation (4). If the light source is located in a hydrogen cloud, the column number density given by the fitting of Equation (16) should be underestimated, because resonant scattering helps the transfer of resonant photons.

An example Figure 6 gives (a) the red damping wings of Equation (16) with τ₀ = 2 × 10² and (b) a solution f of Section 3.4 at r = 10⁴, or τ₀ = 10⁴ and η = 9 × 10⁴. The profiles of panels (a) and (b) are very different from each other. The former is quickly dropping when |x| is less than about 3, while the later at |x| < 3 is softly dependent on x. The curve of panel (b) approaches a much higher bottom at x = 0.

Zoom In Zoom Out Reset image size

More importantly, Figure 6 shows that the curve (a) with τ₀ = 2 × 10² is more or less close to the curve (b). That is, the DLA model at τ₀ = 2 × 10² may give a similar data fitting as the solution with resonant scattering. Therefore, for optical depth τ₀ = 10⁴, the DLA model of Equation (16) will underestimate the column number density of neutral hydrogen by about two orders.

If the light source is significantly time dependent, like the optical afterglow of GRBs, we can treat the source as a flash described by boundary conditions as

$$\begin{equation} j(\eta, 0, x)=0, \quad f(\eta, 0, x) = \left\lbrace \begin{array}{@{}ll} S_0 \phi _s(x), & \mbox{$0<\eta <\eta _0$} \\ \ms 0, & \mbox{$\eta >\eta _0$}. \end{array} \right. \end{equation} \tag{ 17 }$$

It describes a flash within a time range 0 < η < η₀ or the time duration is Δη = η₀. We consider the case of r ⩾ η₀. That is, the size of the halo is larger than the spatial range lasted by the flash, as photon spatial transfer in an optically thick medium cannot be faster than the speed of light. The initial condition is still the same (Equation (13)).

Figure 7 presents two time-dependent solutions of the mean intensity j and flux f: one is for a flash of Equation (17) with η₀ = 100 at r = 10²; the other is for a flash with η₀ = 500 at r = 10³. Other parameters are S₀ = 1, a = 10⁻², and $\phi _s(x)=(1/\sqrt{\pi }) \, e^{-x^2}$. We still see the typical flat plateau of j at all times, even when the original time duration of the flash is as short as Δη = 100.

Zoom In Zoom Out Reset image size

The time dependence of j has a rising phase and a decaying phase. The thermalization of j is held in both rising and decaying phases. We see from Figure 7 that the rising and decaying phases are frequency dependent. Photons at the two peaks reach their maximum at about η = 200 (top right panel) and η = 4000 (bottom right panel), while the flat plateau reaches their maximum at about η = 500 (top right panel) and η = 6000 (bottom right panel). That is, the halo holds a locally thermalized photons for a much longer time than the original time durations Δη = 100 (top panel) and Δη = 500 (bottom panel).

The time evolution of f also consists of a rising phase and a decaying phase. Therefore, the flux emergent from the halo is also a flash. However, the time duration is very different from the original one. For the top panel, we see that the profile of f is almost time independent from η = 300 to 500. That is, the time duration 500 − 300 = 200 is much larger than the original one Δη = 100. For the bottom panel, the time independent flux is from η = 3000 to 5000, and therefore the time duration of the flash is about 2000, which is also much larger than the original time duration η = 500.

To further study the feature of time dependence of flash, we plot Figure 8, which gives the light curve of the flux f at x = 0 (top panels) and the total flux (bottom panels). The top panels show the light curve of rising and decaying phases of the ν₀ photon flux f. With these light curves, one can find the maximum of the flux f at time η_max; the rising phase is then η < η_max and decaying phase is η>η_max. For each curve, one can also define a time duration Δη = η₂ − η₁, where η₁ < η_max and η₂>η_max, and both are given by the condition f(η_1,2, r, 0) = 0.9f(η_max, r, 0). The time duration is then Δη = η₂ − η₁.

Zoom In Zoom Out Reset image size

The two top panels of Figure 8 are for a flash source with original time duration Δη = 50. Figure 8 shows that the time duration of the flash, Δη, is r-dependent. At r = 0, i.e., at the source, Δη is 50. At r = 10² (top left of Figure 8), Δη is about 200, while at r = 10³, we have Δη ≃ 2000. Therefore, the time duration Δη increases with r. This result shows that the time duration seems to be dependent mainly on the size r (or optical depth τ) of the halo, regardless the initial time duration Δη = η₀.

The bottom two panels of Figure 8 are the light curves of the total flux F_t(η, r) = ∫f(η, r, x)dx. The peak times η_max of total flux generally are less than that of ν₀ photon flux, as the ν₀ photon needs longer restoration time. The time durations given by the total flux are also a little less than that of ν₀ photon flux, but it still shows that the time duration is mainly dependent on the size r of the halo.

Figure 9 presents the light curves of flash sources with time duration Δη = η₀ = 1, 5, and 20. The halo size is r = 10². For the case of η₀ = 1, we have η₀ ≪ r. Therefore, the source can be considered as a pulse. The top panels of Figure 9 are for flux of ν₀ photons, while the bottom panels are the corresponding total flux. Although the three flash sources have very different time durations at r = 0, their light curves of f at x = 0 are very similar. The maximum values of f for η₀ = 1, 5, and 20 can even be described by relations as f₂₀ ≃ 4f₅ ≃ 20f₁. The coefficients 4 and 20 are from the ratio of the total numbers of photons of the three flashes. The three light curves of F_t at the bottom panels of Figure 9 are also similar to each other, although they are not as good as the top three curves f. This is because the light curves are frequency dependent. Nevertheless, we still see that the three maximum values of F_t also roughly satisfy the relation F₂₀ ≃ 4F₅ ≃ 20F₁.

Zoom In Zoom Out Reset image size

Either Figure 8 or Figure 9 reveals that the time scales of the propagation of a flash in halos are mainly dependent on the size r, regardless the original time duration. This feature indicates that the spatial transfer of a flash essentially is a diffusion process. As mentioned in Section 3, the spatial transfer of resonant photons cannot be described as a purely Brownian diffusion process, by which the time duration Δη should increase with r² or τ². On the other hand, if the spatial transfer of a photon can be, on average, described by a constant speed, the time duration Δη of a flash should be, on average, equal to the original time duration or, at least, dependent on the original time duration. However, Figures 8 and 9 show that the time duration Δη is roughly proportional to r and independent of the initial time duration. Therefore, the spatial transfer of resonant photons essentially should still be a diffusion process.

The effect of trapping and storing photons in an optically thick halo can be more clearly revealed with a flash source. Let us consider a flash source with a continuant spectrum, which is described by Equation (17), but we take S₀ϕ_s(x) = 1. The time evolutions of j and f for η₀ = 50 at r = 10² are presented in Figure 10, in which we take the time to be η = 200, 300, and 500.

Zoom In Zoom Out Reset image size

The top panels show the evolution of the mean intensity j in the halo. In an early time η = 200, there are photons in the central range |x| < 4 as well as in wings |x|>4. At later time η = 300, wing photons disappear, because all wing photons from the flash source have already escaped from the r halo. At time η = 500, the mean intensity j of |x| < 4 is still about the same as j at η = 300. The flat plateau of j is shown at all times. Therefore, the locally thermalized photons are stored in the halo at least from time η = 200 to 500, which is much longer than the η₀ = 50.

The evolution of the flux f given by the bottom panels of Figure 10 is more interesting. At the time η = 200, the flux shows a typical absorption line at ν₀. However, at η = 300, the flux f becomes a typical emission line with two-peak profile. At time η = 500, f is still a two-peak emission line. The flux of the emission at η = 500 is as strong as that at η = 300. Its light curve is similar to that of Figures 8 and 9. Note that Figures 8 and 9 are for a source of emission line, while Figure 10 is for a continuant spectrum. The similarity of the light curves of Figure 10 with Figures 8 and 9 is again due to the local thermalization and the diffusion in the physical space, both of which lead to the initial frequency spectrum and the time dependence of the photon sources being forgotten.

The emission at η>300 is a delayed emission, as the flash of source has already ceased. The photons of the delayed emission are provided by the |x| < 4 photons stored in the halo r < 10². The time duration of the delayed emission is about the same as the time scale of the decaying phase of Figures 8 and 9. Therefore, it is also proportional to r, regardless the original time duration. Thus, at large r, a flash with a continuant spectrum and very short original time duration η₀ can produce a two-peak emission with time duration proportional to r.

The spatial 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 $\frac{\partial j}{\partial x}$,

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

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

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

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

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

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

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

and the nonlinear weights ω_m are given by

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

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

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

and the smoothness indicators β_l are given by

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

To approximate the r-derivatives in the system of Equations (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{equation} {\hbox{\bf u}}_t + A \hbox{\bf u}_r, \end{equation} \tag{ A6 }$$

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

A fifth-order conservative finite difference WENO scheme can only be applied to a uniform grid or a smoothly varying grid. A smooth transformation,

$$\begin{equation} \xi = \xi (r), \end{equation} \tag{ A7 }$$

gives us a uniform grid in a new variable ξ. In this case, ξ is sufficiently smooth, i.e., it has as many derivatives as the order of accuracy of the scheme. Therefore, the left-hand side of Equations (9) and (10) can be written as

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

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

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

is transformed to

$$\begin{equation} \hbox{\bf u}_t + A \xi _{r}\hbox{\bf u}_{\xi } \end{equation} \tag{ A9 }$$

and the WENO derivative approximation is now applied to $\hbox{\bf u}_{\xi }$.