LIGHT-CONE EFFECT OF RADIATION FIELDS IN COSMOLOGICAL RADIATIVE TRANSFER SIMULATIONS

Abstract

We present a novel method to implement time-delayed propagation of radiation fields in cosmological radiative transfer simulations. Time-delayed propagation of radiation fields requires construction of retarded-time fields by tracking the location and lifetime of radiation sources along the corresponding light-cones. Cosmological radiative transfer simulations have, until now, ignored this "light-cone effect" or implemented ray-tracing methods that are computationally demanding. We show that radiative transfer calculation of the time-delayed fields can be easily achieved in numerical simulations when periodic boundary conditions are used, by calculating the time-discretized retarded-time Green's function using the Fast Fourier Transform (FFT) method and convolving it with the source distribution. We also present a direct application of this method to the long-range radiation field of Lyman-Werner band photons, which is important in the high-redshift astrophysics with first stars.

1. INTRODUCTION

Cosmological radiative transfer simulations are used to study astrophysical processes which occur on a large, cosmic scale. A notable example of such processes is the process of cosmic reionization, in which individual H II regions take up large volumes, of order (~ 20 comoving Mpc)3 or larger, after ~ 50% of all the baryons are reionized by astrophysical radiation sources (Furlanetto & Oh 2005; Furlanetto et al. 2006; Iliev et al. 2007; Zahn et al. 2007). The typical radiation sources long after the recombination epoch are stars and quasars, which emit predominantly ultra-violet (UV) and X-ray photons, respectively. The mean free path of hydrogenionizing (H-ionizing) UV photons is simply the average size of H II regions,1 while the mean free path of Xray photons is much larger than that of UV photons due to the relatively small optical depth (the scattering cross section of UV photons is much larger than that of X-ray photons: e.g., Mesinger et al. 2013; Xu et al. 2014). Therefore, the mean free path of X-ray photons is truly cosmological, surpassing several hundred Mpc when the rest-frame photon energy is larger than a few keV (Xu et al. 2014). Another example of long mean free paths is the Lyman-Werner (LW) band photons. These photons lie at the energy range of 11 – 13.6 eV, and thus traverse cosmological distances freely until they redshift into hydrogen Lyman resonance lines and are scattered off neutral hydrogen atoms (Haiman et al. 2000; Ahn et al. 2009). Practically all LW band photons are scattered multiple times and reprocessed into lower-energy photons when they traverse the “LymanWerner horizon” rLW ~ 100 [21/(1 + z)]0.5 Mpc, where z is the redshift at which the photons are emitted (Ahn et al. 2009).

The light-cone effect, which denotes the effect related to fields travelling at the speed of light and thus tracing past-time events along the corresponding light cones, is inherent in any physical processes on cosmic scales. Radiation fields and gravitational fields are the obvious examples governed by the light-cone effect: both fields travel at the speed of light (see Hwang et al. 2008 showing that the gravitational field travels at the speed of light in the form of the electric part of the Weyl tensor). The relevant time or length scale over which the light-cone effect matters can be estimated by an effective time-scale tlc ≡ f /(df /dt), where f is any physical quantity relevant to the problem of interest and df /dt is the change rate of f. For astrophysical problems, the time scale may be naively estimated by the lifetime of radiation sources of interest. In contrast, gravitational fields under non-relativistic motion (velocity ≪ c) will be corrected from the Newtonian field, which is constructed in an action-at-a-distance manner, only very slightly (about 10−6 – 10−4 or somewhat larger if cumulative effects are considered) by general relativistic effects including this light-cone effect (Hwang et al. 2008).

Light-cone effect of radiation fields traversing cosmological distances has been usually approximated by a uniform, spatially-averaged quantity in semi-analytical studies (e.g., Haiman et al. 2000 for LW and H-ionizing backgrounds), or calculated by computationally expensive, brute-force methods of ray-tracing with finite speed of light (e.g., Bryan et al. 2014 for H-ionizing background). The first account of the inhomogeneous LW background in numerical simulations was taken by Ahn et al. (2009), whose consistency with semianalytical and semi-numerical approach was confirmed elsewhere (Dijkstra et al. 2008; Holzbauer & Furlanetto 2012). The method of Ahn et al. (2009) was later used in simulating the process of cosmic reionization under the influence of the first stars (Ahn et al. 2012; Fialkov et al. 2013). This method has been improved for faster computation using FFT. Its application to other radiation fields affected by the light-cone effect is straightforward, as long as the field strength from the source is only a function of distance but not of direction. Even though this improved method was already used in the simulation by Ahn et al. (2012), its formalism and generality have never been described explicitly. Thus we lay out the generic formalism, a numerical scheme and the actual application in this paper. A similar account for X-ray background (Xu et al. 2014) used the same machinery of Ahn et al. (2009) with modifications to implement the opacity due to hydrogen and helium atoms.

This paper is organized as follows. In Section 2, we describe the generic formalism and the numerical scheme. In Section 3, we show how the LW background is calculated using this scheme, and also quantify the possible error from ignoring the light-cone effect. In Section 4, we conclude our findings and discuss further prospects and limitations.

 

2. RETARDED RADIATION FIELD: FORMALISM AND NUMERICAL METHOD

2.1. Green’s Function Formalism

A generic radiation field F(x, t) at position x and time t in a homogeneous and isotropic background spacetime satisfies a d’Alembertian-type equation. In the trivial case of a Minkowski background, electromagnetic (EM) waves (see Jackson 1998) follow

where F(x, t) is either the scalar potential Φ or the vector potential A, and f(x, t) is the corresponding source term (the charge density or the current density, respectively). In this case, F(x, t) can be cast into an convoluted integral form

where the retarded Green’s function2 (RGF henceforth) G(x, t; x′, t′) satisfies

and is given in a closed form as

Note that Equation (2) is in the most generic form. The light-cone effect or causality is apparent in Equations (2) and (4), as the field at observing time t is caused by a source along the light cone at past time t′ = t−|x − x′| /c. Even though Equations (1), (3) and (4) are applicable both to Φ and A only in the Lorentz gauge, it is well known that the electromagnetic field, as a physical quantity derived from both Φ and A, still conserves this causality under any gauge choice (Jackson 2002; see also a nice example of Problem 6.20 in Jackson 1998).

Let us consider a generic field F(x, τ) in the spacetime described by the Robertson-Walker metric with null curvature,

where x is the comoving spatial coordinate, a(t) is the scale factor at conformal cosmic time τ with normalization a = 1/(1 + z), and τ is defined by

It is convenient to use τ and x because the null geodesic is simply given by d |x| /dτ=1. One can then expect that the generic RGF will take the following form:

where 𝓕(|x − x′|, τ, τ′) is a generic, geometrical dilution factor. Note that 𝓕(|x − x′|, τ, τ′), in general, reflects cosmological effects such as photon-redshifting, and therefore depends also on τ and τ′. The corresponding integral form will become (cf., Equation 2)

where f(x′, τ′) is the source term properly defined in the (x, τ) coordinate system.

We now decompose Equation (8) into discrete terms which are affected by different time slices, and also apply a periodic boundary condition to assimilate the usual simulation set-up. Equation (8) can be rewritten

where f(x′, τn) approximates the source term during τ′ = [τn, τn+1] as a fixed 3D quantity, W(x) is a tophat window function such that

and the RGF and the partial contribution to F(x, τ) from the sources during τ′ = [τn, τn+1] are defined as 𝓕n(|x − x′| ; τ) and Fn(x, τ), respectively. It is important to choose small enough time gaps such that f(x′, τ′) does not evolve much during each time step. Now, we apply the periodic boundary condition: the cubical domain with size L (volume V ≡ L3) is infinitely repeated side to side and thus satisfies

where u, v, w are integer indices running from −∞ to ∞ along the three axes of the cubical domain, both x and x′ are confined to one domain, the second identity is guaranteed by the periodicity of the source term, and the last identity defines the “effective” RGF 𝒢n(x − x′). Note that for a given observing point (x, τ) and the past time-slice at [τn, τn+1], the (u, v, w) points contributing to 𝒢n(x − x′) are only those confined in the spherical shell bounded by comoving radii τ − τn+1 and τ − τn.

We have just cast the full integral (Equation 8) for the retarded field into the summation (Equation 9) of spatially convoluted integrals (Equation 10). In the next section, we will describe how we calculate the full retarded field under a uniform Eulerian grid.

2.2. Numerical Method

We now develop a numerical method to calculate the retarded fields. A scheme for calculating 3D fields in numerical simulations becomes the most transparent under a uniform grid, and thus we restrict the description also to the case of a uniform grid. This choice is further justified by the fact that a retarded field originates from sources located far from the observing point, whose impact on each observing point is greatly averaged out and thus can be accurately calculated even on a uniform grid.

If we discretize relevant quantities on a uniform grid, Equation (10) turns into

where the domain is divided into cells with size l ≡ L/N, {i, j, k} and {i′, j′, k′} are integer indices for real-space positions, and we absorbed l3 into fn(i′, j′, k′) by letting ∫ d3x′ f(x′, τn) ↦ ∑ i′, j′, k′ fn(i′, j′, k′). If (p, q, r) is the Fourier component of some field A(i, j, k) given by

where the imaginary unit i is shown in bold face to avoid confusion with the index i, then the convolution theorem gives

where {p, q, r} are the indices of the Fourier-space coordinates.

The merit of Equations (11) – (13) is that the retarded field can now be calculated by FFT as long as periodicity is assumed. The algorithm then becomes straightforward, as follows:

For given (observing) time τ, divide the lookback time into discrete time slices with τ′ = [τn, τn+1]. For a given (past) time slice n, locate a unit source at the origin (i′, j′, k′) = (0, 0, 0) of the realspace domain (box), and attach identical boxes – each containing one unit source at its origin – side by side. The total number of boxes to attach may be limited by the “effective horizon”, beyond which the contribution to 𝒢n(i, j, k) (or to the net Fn(i, j, k) due to a rapid decrease of fn(i′, j′, k′) in lookback time) becomes too weak to include. Calculate 𝒢n(i, j, k) from this configuration by a direct summation (see its definition in Equation 10) over the contributions from these sources, when the functional form of 𝓕(|x − x′|, τ, τ′) in Equation (9) is known. Note that 𝒢n(i, j, k) at each observing point {i, j, k} is affected only by the unit sources inside the spherical shell given by τ − τn+1 ≤ ros ≤ τ − τn, where ros is the comoving distance between x = (i, j, k)l and x′ = (u, v, w)L and the index notation follows those of Equations (10) and (11). Get (p, q, r) by applying FFT to 𝒢n(i, j, k) calculated from step 2. Get (p, q, r) by applying FFT to fn(i, j, k). Get (p, q, r) from Equation (13), and then apply FFT to get Fn(i, j, k). Finally, sum over time slices to get F(i, j, k) = ∑n Fn(i, j, k).

Note that once the box size, the grid resolution and the output times are determined, one can precalculate 𝒢n(i, j, k) for each set of observing time (τ) and lookback time (τn–τn+1) and store all (p, q, r)’s regardless of the source distribution, which is a merit of the convolution theorem. The number of data files for 𝒢n(i, j, k) is ~ NobserveNlookback, where Nobserve is the number of observing redshifts and Nlookback is the number of lookback times. Nlookback can depend on the observing redshift, depending on the change of tlc over time. When the box size (L) is much smaller than the size of the effective horizon (rh), the number of unit sources repeated around the computation domain will be Nunit−source ~ (rh/L)3, and the computation of each 𝒢n(i, j, k) will require ~ Nmesh(r)Nunit−source operations.

 

3. APPLICATION: LYMAN-WERNER BACKGROUND

As an application, we apply the method from Section 2.2 to the calculation of LW background. LW bands are composed of many excitation levels of hydrogen molecules (H2), and at high intensities the LW photons can photo-dissociate H2. Because first stars are born with the help of H2-cooling, calculating LW background is crucial in the study of first-star formation.

The functional form of the RGF for LW background was first laid out by Ahn et al. (2009). The band-averaged and angle-averaged intensity JLW(x, z) observed at comoving location x and redshift z due to a source at (x′, zs) in units of erg s−1 cm−2 Hz−1 sr−1 is given by

where 〈Lν(x′, zs)〉 d(hpν)Lν(x′, zs)/(2.1eV) is the rest-frame, frequency(ν)-averaged source lumi- nosity of the emitted LW band photons at the source redshift zs and hp is the Planck constant,

and

with the Hubble constant (in units of 100 km/s/Mpc) h and the present ratio of matter density to the critical density Ωm. Of course, the light-cone effect is implicit here: |x − x′| = τ(z)−τ(zs). The LW horizon is rLW ≡ 97.39α, beyond which no sources contribute to JLW. Comparing Equation (14) to Equations (7) and (8), we find that

if we take the source term to be a sum of individual source luminosities such that f(x′, τ′) = .

When visualized, the real-space RGF 𝒢n(i, j, k) can give some insight on how JLW is constructed. Figure 1 depicts RGFs for JLW, which clearly show how the region of influence is determined by the light-cone effect from a unit source at one corner of the box3. Following the steps described in Section 2.2, we have calculated JLW and used this to study the formation and suppression of the first stars during the epoch of reionization (Ahn et al. 2009, 2012). Ahn et al. (2012) simulated the process of cosmic reionization including the rise and fall of first stars, in a box of size 114/h Mpc, on which Figures 1 – 3 are based. Lifetime of first stars, which are believed to be as massive as M∗ 100 M⊙ on average (e.g., see a recent result by Hirano et al. 2014), is only about a few Myrs or less, and thus the life and death of these stars render the source function to change rapidly in time and space. In Figure 2, the rapid evolution of 〈Lν〉 every 2 Myrs due to first stars is shown, together with the partial contributions Fn(p, q, r) to JLW. rLW resides within the lookback time of 10 Myrs. One can see that contributions to JLW becomes weaker and more homogeneous as lookback time increases. Nevertheless, contribution of those time slices with lookback time tlookback = [2 − 10] Myrs (2nd to 5th panels from left in Figures 1 and 2) to JLW is not completely negligible.

Figure 1.RGFs for LW intensity at observing redshift z = 27.90, in a periodic box of size 114/h (comoving) Mpc and a uniform grid of 2563 cells. Both the real-space (top row) and Fourier-space (bottom row) RGFs (𝒢n and respectively) are shown in log scales, on a slice containing 4 corners of the box. For which has large dynamic range, we took the absolute value of its real part, ⎢Re() ⎢, for better visualization. From left to right, RGFs of progressively longer lookback-time slices are shown, in units arbitrarily chosen but consistent throughout time slices. Lookback-time range is shown in terms of redshift, below each set of 𝒢n and . Each redshift range corresponds to 2 Myrs of cosmic time. Note that the right-most panels correspond to the time slice containing the horizon rLW.

Figure 2.Source distribution and LW background from different time slices, with the same configuration as in Figure 1. The top row shows the 2D maps of 〈Lν〉 (assumed to be fixed at and beyond the redshift specified below each panel until the adjacent future redshift) on a slice of thickness (114/256)/h Mpc, and the bottom row the “partial” JLW’s on the same slice contributed from the look-back time ranges specified below each panel. The net JLW observed at the observing redshift z = 27.90 is given by the summation of all the partial JLW’s, which is shown in Figure 3.

Figure 3.Comparison of JLW’s with the light-cone effect and without. JLW under the correctly implemented light-cone effect (left), JLW under an action at a distance (middle; see Equation 18) and the fractional difference between these two fields (right) are plotted in the same 2D slice as in Figures 1 and 2.

How important is the light-cone effect in the LW background from first stars? To answer this question, we now calculate JLW in an action-at-a-distance way, and then compare this to JLW correctly calculated including the light-cone effect. For this, we now freeze the source distribution at the observing redshift, while attaching identical boxes side by side for periodicity. Basically, we now calculate this “instantaneously affected” JLW as

where 〈Lν〉 is now frozen at the observing redshift z (e.g., the source term corresponding to the top-leftmost panel in Figure 2), but we keep all other terms intact including the implicit condition |x − x′| = τ(z)−τ(zs), in order to quantify only the impact of the temporal evolution of the source term. Figure 3 shows that the rapid evolution of source term indeed imprints its signature in JLW. The fractional difference fdiff ≡ (JLW, instant−JLW)/JLW is about ~ 10%, which is quite significant. Therefore, astrophysical processes governed by radiation sources of short lifetime should be affected by the light-cone effect in general.

 

4. DISCUSSION

We have formulated a generic method to numerically calculate radiation fields which are affected by the lightcone effect. This method requires discretizing both time and space coordinates, and is analogous to the usual mesh scheme used in N-body gravity simulations. The light-cone effect is realized by the retarded-time Greeen’s function, and its calculation becomes efficient when a periodic-box condition is used. Once the RGFs are calculated, the process of obtaining a radiation field of interest is performed by FFT, which greatly expedites computation compared to direct summation over all sources inside and outside the simulation box.

The relevance of this method is clear in astrophysical problems occurring on cosmological scales. Once the size of the simulation box is chosen and the functional form of the RGF is known, one can follow the algorithm in Section 2.2, as in the example of LW background in Section 3. Similar large-scale astrophysical processes which require an accurate account of the light-cone effect are the reionization of baryons by UV photons and the reionization and heating of gas by X-ray photons. It becomes more important at very high redshift when most astrophysical sources are first stars and their byproducts, because the lifetime of first stars are believed to be very short and thus the source function evolves rapidly in time (Section 3).

There are, of course, limitations to this method. The method is based on the fact that the RGF is a function only of distance to the source but not of direction to the source. In reality, radiative transfer is usually affected by the actual path of the ray in the form of the locallyvarying optical depth. For example, formation of an H II bubble from a source does not occur in a spherically symmetric way, even though light travels at the speed of light. This is due to the inhomogeneous distribution of hydrogen density (and also the optical depth), which then regulates the H-ionizing background to be dependent on the direction of the ray. For H-ionizing background, our method can then be used only after the IGM is fully ionized and the optical depth becomes negligible everywhere. In case of LW background, we could safely use the method because the modulation factor fmod in Equation (15), derived from the cosmic redshifting of LW photons and the frequency of hydrogen Lyman resonance lines, is valid even when the IGM is ionized, because simulations of cosmic reionization find that the IGM leaves trace amount of neutral hydrogen atoms due to recombination (e.g., Iliev et al. 2008). In this case, very-high opacity of Lyman resonance lines guarantees that once LW band photons are redshifted into one of these lines, even a trace amount of hydrogen atoms can successfully scatter off those photons, and then fmod becomes the only modulation factor in addition to the usual geometrical dilution factor (Equations 14 and 17). Another limitation is that the method cannot treat the gravitational lensing effect, which depends on the actual mass distribution around each ray’s path. If a problem of interest is found to be strongly affected by the lensing effect, one should resort back to the ray-tracing method. In addition, if source locations are resolved beyond the mesh resolution, one can adaptively choose to use the particle-mesh (PM) scheme to calculate the near-source influence by direct summation over nearby-source contributions while calculating far-source influence by our method.