3D numerical cosmological radiative transfer in an inhomogeneous medium
Abstract
A numerical scheme is proposed for the solution of the threedimensional radiative transfer equation with variable optical depth. We show that timedependent ray tracing is an attractive choice for simulations of astrophysical ionization fronts, particularly when one is interested in covering a wide range of optical depths within a 3D clumpy medium. Our approach combines the explicit advection of radiation variables with the implicit solution of local rate equations given the radiation field at each point. Our scheme is well suited to the solution of problems for which line transfer is not important, and could, in principle, be extended to those situations also. This scheme allows us to calculate the propagation of supersonic ionization fronts into an inhomogeneous medium. The approach can be easily implemented on a single workstation and also should be fully parallelizable.
keywords:
radiative transfer – methods: numerical – cosmology: theory – intergalactic medium1 Introduction
Understanding the effect of the radiation field on the thermal state of interstellar and intergalactic gas is important for many areas of astrophysics, and in particular for star and galaxy formation. Of special interest for cosmological structure formation are the epoch of reionization of the Universe (Gnedin & Ostriker 1997, MiraldaEscudé et al. 1998), the effects of selfshielding on the formation of disk and dwarf galaxies (Navarro & Steinmetz 1997, Kepner et al. 1997) and the absorption properties of Ly clouds (Katz et al. 1996, Meiksin 1994). While our knowledge of the physics of largescale structure and galaxy formation has benefited significantly from numerical Nbody and gasdynamical models (see, e.g., Zhang et al. 1998, and references therein), there is very little that has been done to include radiative transfer (RT) into these simulations. Challenges seem to abound, not least of all, the fact that the intensity of radiation in general is a function of seven independent variables (three spatial coordinates, two angles, frequency and time). While for many applications it has been possible to reduce the dimensionality (e.g. to build realistic stellar atmosphere models), the clumpy state of the interstellar or intergalactic medium does not provide any spatial symmetries. Moreover, coupled equations of radiation hydrodynamics (RHD) have very complicated structure, and are often of mixed advectiondiffusion type which makes it very difficult to solve them numerically. Besides that, the radiation field in optically thin regions usually evolves at the speed of light, yielding an enormous gap of many orders of magnitude between the characteristic timescales for a system.
One way to avoid the latter problem is to solve all equations on the fluidflow timescale. While there are arguments which seem to preserve causality in such an approach [Mihalas & Mihalas 1984], even then the numerical solution is an incredibly difficult challenge [Stone et al. 1992]. On the other hand, several astrophysical problems allow one to follow the system of interest on a radiation propagation timescale, without imposing a prohibitively large number of time steps. In the context of cosmological RT, we would like to resolve the characteristic distance between the sources of reionization. Cold dark matter (CDM) cosmologies predict the collapse of the first baryonic objects as early as , with the typical Jeans mass of order (see Haiman & Loeb 1997, and references therein). The corresponding comoving scale of fragmenting clouds is . If stellar sources reside in primordial globular clusters of this mass, then for the average separation between these objects is (comoving). For explicit schemes the Courant condition imposes a timestep
(1) 
Thus evolution from to takes (assuming density parameter , and defining the Hubble constant to be ). Substituting and into eq. (1) yields the total number of timesteps in the range over the entire course of evolution. In other words, to resolve reionization by stellar clusters, we need to compute few tens of thousands of timesteps on average, and for clusters the number of steps required will be ten times smaller. A grid with the required resolution will result in computational boxes of several on a side. A full cosmological radiative transfer simulation with boxes at least this big and for all the required timesteps has not been feasible in the past.
This is far from the only challenge we face. Since any two points can affect each other via the radiation field, even for a monochromatic problem, we must describe the propagation of the radiation field anisotropies in the full fivedimensional space. Standard steadystate RT solvers, which have been widely used in stellar atmosphere models, are not efficient in this case. Nonlocal thermodynamic equilibrium (NLTE) steadystate radiative transport relies on obtaining the numerical solution via an iterative process for the whole computational region at once, and is usually effective only for very simplified geometries. Any refinement of the discretization grid and/or increase in the number of atomic rate equations to compute NLTE effects will necessarily result in an exponential increase in the number of iterations required to achieve the same accuracy. On the other hand, the 3D solution of the steadystate transfer equation in the absence of any spatial symmetries can often be obtained with Monte Carlo methods [Park & Hong 1998]. However, these methods demonstrate very slow convergence at higher resolutions and are hardly applicable if one is interested in following a timedependent system.
The change in the degree of ionization in a lowdensity environment occurs on a radiation propagation timescale . To track ionization fronts (Ifronts) in this regime, it is best to apply a highresolution shockcapturing scheme similar to those originally developed in fluid dynamics. One possible approach is the direct numerical solution of the monochromatic photon Boltzmann equation in the 5dimensional phase space (Razoumov, in preparation). To allow for a tradeoff between calculational speed (plus memory usage) and accuracy, a more conventional approach is to truncate the system of angleaveraged radiation moment equations at a fixed moment and to use some closure scheme to reconstruct the angle dependence of the intensity at each point in 3D space. The method of variable Eddington factors first introduced by Auer & Mihalas (1970) has been shown to produce very accurate closure for timedependent problems in both 2D [Stone et al. 1992] and 3D [Umemura et al. 1998]. However, to the best of our knowledge, all schemes employed so far for calculating the timedependent variable Eddington factor were based on a steadystate reconstruction of the radiation field through all of the computational region at once, given the thermal state of material and level populations at each point. Since advection (or spatial transport) of moments is still followed on the timescale of typical changes in the ionizational balance of the system, this approximation certainly provides physically valid results, assuming that the reconstruction is being performed often enough. However, the steadysteady closure relies on the iterative solution of a large system of nonlinear equations, which becomes an exceedingly difficult problem, from the computational point of view, as one moves to higher spatial and angular resolution and to the inclusion of more complicated microphysics.
The goal of the present paper is to demonstrate that in the cosmological context it is possible and practical to solve the whole RT problem on a radiation propagation timescale – as opposed to the fluidflow timescale – and we present a simple technique which gives an accurate solution for the angledependent intensity in three spatial dimensions. The scheme can track discontinuities accurately in 3D and is stable up to the Courant number of unity. Since all advection of radiation variables is being done at , the scheme is well tailored to the numerical study of the propagation of Ifronts into a nonhomogeneous medium with any optical depth, and gives very accurate results for scattering. Rather than solve the proper chemistry equations applicable to cosmological structures, we adopt a simple toy model described below. Similarly, in the current paper we do not make any attempt to model the thermal properties of the gas, concentrating just on efficient multidimensional advection techniques.
This paper is organized as follows. In Section 2 we briefly review the state of numerical RT in the study of reionization. We then concentrate on methods for 5D numerical advection. In Section 3 we describe our numerical algorithm and we present the results of numerical tests in Section 4. Finally, in Section 5 we discuss the next steps towards a realistic 3D RT simulation.
2 Formulation of the problem
It is believed that light from the first baryonic objects at led to a phaselike transition in the ionizational state of the Universe. This process of reionization significantly affected the subsequent evolution of structure formation (Couchman & Rees 1986). In detail reionization did not happen at a single epoch, with details of ‘preheating’, percolation, helium ionization and other physical processes having been studied in great detail over the last decade (some recent contributions include Madau et al. 1997, Haiman & Loeb 1997, Gnedin & Ostriker 1997, Shapiro et al. 1998, Tajiri & Umemura 1998).
It now seems clear that the full solution of the problem requires a detailed treatment of the effects of RT. To complicate matters, by the time of the first star formation, the smallscale density inhomogeneities had entered the nonlinear regime [Gnedin & Ostriker 1997], and the medium was filled with clumpy structures. The success of cosmological Nbody and hydrodynamical models in quantifying the growth of these objects (e.g., Zhang et al. 1998) suggests that the next step will be to include the effects of global energy exchange by radiation. Indeed, there is a need for timedependent 3D RT models as numerical tools for understanding the effect of inhomogeneities in the dynamical evolution of the interstellar/intergalactic medium. For instance, the ability of gas to cool down and form structures depends crucially on the ionizational state of a whole array of different chemical elements, which in turn directly depends on the local energy density of the radiation field.
The hydrogen component of the Universe is most likely ionized by photons just above the Lyman limit, because (1) the crosssection of photoionization drops as at higher frequencies, and (2) the medium will be dominated by softer photons, even in the case of quasar reionization (when ionizing photons come mostly from diffuse H ii regions). Therefore, we argue that either monochromatic or frequencyaveraged transfer will be a fairly good approximation in our models.
Recently, the problem of simulating 3D inhomogeneous reionization with realistic radiative transfer has attracted considerable interest in the scientific community. Umemura et al. (1998) calculated reionization from to , solving the 3D steadystate RT equation along with the timedependent ionization rate equations for hydrogen and helium. The radiation field was integrated along spatial dimensions using the method of short characteristics [Stone et al. 1992]. The steadystate solution implies the assumption that the radiation field adjusts instantaneously to any changes in the ionization profile. One drawback of this approach, however, is in lowdensity voids where there are probably enough Lyman photons to ionize every hydrogen atom, so that the velocity of Ifronts is simply equal to the speed of light. Then the rate equations still have to be solved on the radiation propagation timescale. Besides, implicit techniques in the presence of inhomogeneities will become exponentially complicated, if we want to solve timedependent rate equations for multiple chemical species.
Norman et al. (1998) and Abel et al. (1998) present a scheme for solving the cosmological radiative transfer problem by decomposing the total radiation field into two parts: highly anisotropic direct ionizing radiation from point sources such as quasars and stellar clusters, and the diffuse component from recombinations in the photoionized gas. In their method the direct ionizing radiation is being attenuated along a small number of rays, each of which is forced to pass through one of the few point sources within the simulation volume. The diffuse part of the radiation field is found with a separate technique which can benefit from the nearly isotropic form of this component, for instance, through the use of the diffusion approximation. Both solutions are obtained neglecting the time dependent term in the radiative transfer equation, with the default timestep dictated by the speed of the atomic processes.
If reionization by quasars alone is ruled out (Madau 1998, however see Haiman & Loeb 1998), then Ifronts will be caused by Lyman photons from lowluminosity stellar sources at high redshifts. In this case the pressure gradient across the ionization zone is more likely to become important before the front is slowed down by the finite recombination time. In the present paper we ignore hydrodynamical effects, concentrating on an efficient method to track supersonic Ifronts. Our approach is to solve the timedependent RT coupled with an implicit local solver for the rate equations. This method gives the correct speed of front propagation and it also quickly converges to a steadystate solution for equilibrium systems. However, we should note that until a detailed comparison is made between explicit advection (at the speed of light) and the implicit reconstruction (through an elliptic solver), it is difficult to judge which approach works best in simulating inhomogeneous reionization in detail.
Although radiation propagates with the speed of light and the intensity of radiation depends on five spatial variables, plus frequency and time, the RT equation is inherently simpler than the equations of compressible hydrodynamics, since its advection part is strictly linear. Nonlinearities are usually introduced when we are trying to reduce the dimensionality of the problem. Much of the difficulty thus comes from inability to get decent numerical resolution in the 5D (or 6D with frequency) space with presentday computers.
In the current work we have attempted to develop an efficient method to describe the anisotropies in the monochromatic radiation field propagating through an inhomogeneous medium, which we now describe.
3 The numerical technique
The classical RT equation (without cosmological terms) reads
(2) 
where is the intensity of radiation in direction and and are the local emissivity and opacity. This equation is valid for those problems in which the lightcrossing time across the computational volume is much smaller than the Hubble time, and we can neglect the redshift effects within the simulation volume. With proper boundary conditions one can easily account for the Doppler shift in the background radiation.
The basic idea of our technique is to solve eq. (2) directly for the angledependent intensity at each point. The total radiation field at each point is divided into the direct (from ionizing sources) and diffuse (due to recombinations in the gas) components:
(3) 
where are the three indices in a rectangular grid. The energy density due to direct photons coming from point sources of ionization can be easily calculated on the 3D grid via summation over all sources (assuming that there are not too many of them within the volume):
is the optical depth, is the physical distance between the current point and the source and is the age of the source. Note that the rate of emission of photons can be modified to allow for variability of sources on short timescales. For the diffuse component we use an upwind monotonic scheme to propagate 1D wavefronts along a large number of rays in 3D at the speed of light. Following Stone & Mihalas (1992), we apply an operator split explicitimplicit scheme, in which advection of radiation variables is treated explicitly and the atomic and molecular rate equations are solved implicitly and separately at each point. Unlike Abel et al. (1998), we use rays to track the diffuse component of the radiation field and the direct ionizing background radiation (streaming into the computational volume). Since we need to draw rays essentially through every grid point in the 3D volume, at first glance this approach appears to have very large memory requirements. However, efficient placing of the rays can significantly reduce the computational effort.
3.1 A uniform and isotropic grid of rays
At each timestep we are interested in getting a solution for the mean radiation energy density and material properties on a 3D rectangular grid. Instead of shooting rays though each grid node in 3D, we choose to cover the whole computational volume with a separate grid of rays which is uniform both in space and in angular directions. Assuming that the computational volume corresponds to the range , we first construct a 2D rectangular base grid containing nodes with coordinates with the origin at the centre of the cube
(4) 
where the factor ensures that the entire volume is covered with rays, and we shoot rays normal to the base grid through all of its grid points. Note that the separation between 2D nodes is allowed to be larger than . To cover the whole volume with rays, we then rotate the base grid by an angle around the axis and by around the axis, where the rotation angles are discretized to mimic an isotropic distribution of rays (with fewer azimuthal angles close to the poles)
(5) 
Only those rays which pass through the volume (and not all of them do!) are stored in memory. Let the total number of all possible orientations of the base grid be
(6) 
The resulting number of rays is significantly smaller than . In fact, with and (), we only require rays.
Now, that we have the grid of rays and the 3D rectangular mesh, we have to specify the rules of interpolation between them. Before we start our simulations, for each 3D grid point we also store an array of the closest four rays going through its neighbourhood in the direction . Since the rays do not pass exactly through 3D grid points, we use the values of the intensity on the four closest rays in that direction to compute the angulardependent intensity. Assume that for the point the distances to the four closest rays going in the direction are , , and , respectively. We project the point onto these rays and read the values of the intensities, which we write as , , and . We then calculate the intensity at in the direction according to
(7) 
where the weights for are set to zero if rays were found in the immediate neighbourhood of . This might be the case close to the edges of the computational volume; we shall comment more on this while discussing the boundary conditions. The form of eq. (7) was chosen specifically because: (1) if a ray labeled happens to pass exactly through the point , then ; (2) if all four rays encompass the point, ; and (3) if happens to be far from all four rays, then the resulting intensity will just be an average of the four ’s.
At each point on our 3D rectangular mesh we assume a piecewise linear dependence of the intensity on two angles, and ,
(8)  
within a spherical rectangular element (or a spherical triangle adjacent to either of the poles) bounded by the angles , in and and in , with the rectangular grid defined as
(9) 
and
We then integrate the intensity over with appropriate weights to get the scalar radiation energy density at each point of the rectangular spherical grid
(10) 
where the quadrature terms for the integration
are modified to allow for the nonorthogonal angular grid (eq. 5).
Since the advection part on the lefthand side of eq. (2) is strictly linear, the simplest way to propagate intensities is just to shift wavefronts by one grid zone at each timestep, accounting for sources and sinks of radiation. Assuming that all discretization points along each ray are strictly equidistant, the intensity at a point is updated simply as
(11) 
Alternatively, one could take special care of the length of each ray segment contained within a particular 3D grid cell, and use a scheme similar to the thirdorderaccurate piecewise parabolic advection method (PPA) of Stone & Mihalas (1992). In either case we can track sharp discontinuities in 1D with very little numerical diffusion, and, therefore, our approach is well suited to the calculation of Ifronts.
3.2 Local chemistry equations
Since in our calculation all advection of radiation variables is performed explicitly, we can solve NLTE rate equations separately at each point. This makes it relatively easy to implement an implicit solver for all atomic and molecular processes.
To demonstrate the capabilities of explicit advection, instead of solving the proper chemistry equations for multiple species with primordial chemical composition, we have here adopted a simple toy model with just photoionization and radiative recombination in a pure hydrogen medium. The implicit solution of possibly stiff rate equations described below can be implemented in a similar manner for more realistic chemistry models.
The time evolution of the degree of ionization is given simply by
(12) 
with opacity
(13) 
Here is the photoionization coefficient, the energy density of ionizing radiation and the recombination coefficient. The correct expression for the total emissivity of the gas can be obtained by considering conservation of the thermal energy density for matter:
(14) 
where all symbols represent the change of variables during one time step. The full recombination coefficient
(15) 
is the sum of recombination coefficients to the ground state () and to all levels above the ground state (, the ‘case B’ recombination coefficient), is the frequency just above the Lyman limit, and we assume that recombinations in Lyman lines occur on a short timescale compared to . Similarly, the full emissivity (or gas energy loss through recombinations) is
(16) 
where . This simple notation ensures radiation energy conservation in eq. (2) for pure scattering of Lyman continuum photons (i.e., when ). Eq. (12) does not account properly for the number of photons entering the volume, so that a large photoionization coefficient might lead to overproduction of ions. To compensate for this, the number of photoionizations inside a 3D grid cell per unit time is not allowed to be larger than the number of photons actually absorbed inside this cell, i.e.
(17) 
Eq. (12,17) are solved separately at each point, given the local radiation energy density . Discretization of eq. (12) in time yields
(18) 
where is just the righthand side of eq. (12) and for stability. This equation can almost always be solved via Newton’s method for small enough . Linearizing eq. (18), we get the th approximation to the value of at time :
(19)  
which can then be iterated.
3.3 The algorithm
We start calculations by specifying the initial conditions (temperature, degree of ionization, and in the simplest cases no radiation field inside the volume) and boundary conditions (the intensity of radiation entering the volume – all outward flux at the edges can freely escape the computational box). The inward flux at the boundaries is isotropic within and is simply
where is the average background radiation energy density. Since each ray within the volume starts and ends at the sides of the volume, we automatically have boundary conditions for each of the 1D advection problems. The density field is kept static for all tests in this study, but since the radiation field is being evolved explicitly at the speed of light, one could easily evolve the underlying density distribution on a much bigger fluid flow timescale if desired.
The course of the algorithm at each time step can be divided into the following steps:

project , from 3D grid to rays,

evolve wavefronts along individual rays by ,

project intensities from rays to 3D to reconstruct as a function of five variables ,

calculate the radiation energy density ,

solve all chemical equations via an implicit iterative technique, separately at each point, and

compute new , .
At the beginning of each timestep we advect 1D intensities according to eq. (11) along each ray. The numerical resolution along each ray is simply set to the resolution of the 3D rectangular mesh, so that along the ray the point has a coordinate
where is the length of the ray segment inside the cube. Note that one can have much coarser 1D grids along rays, speeding up the advection but sacrificing both spatial and angular resolution. The advected intensities are then projected onto a 3D grid to reconstruct the mean energy density at each point using eq. (7 – 10). This operation is one of the most demanding from the computational point of view, since at each of our points we have to deal with the angular dependence of the radiation field. When this update is done, we solve the matterradiation interaction equations implicitly to compute the local level populations (eq. (19)). This gives us the 3D distributions of emissivity and opacity (eq. (13 – 16)) which are then mapped back to the rays and used in the advection scheme at the next timestep.
This simple scheme which we will refer to as timedependent ray tracing can be used as a standalone solver, or as a closure scheme for the system of moment equations through the use of variable Eddington factors (as in Stone et al. 1992). In the absence of any sinks and sources of radiation, the intensity is conserved exactly along each ray. Since the number of rays does not vary with time, our method guarantees exact conservation of the radiation energy in 3D.
Note that if Ifronts do not propagate fast, for instance, in the cosmological context where the evolution on the lightcrossing timescale will normally require many thousand time steps, it is possible to update the radiation energy density once every few tens or few hundred time steps, while still evolving the intensities and solving all rate equations properly at the speed of light. We found that in practice this shortcut leads to an increase in speed by a factor of 10 or even higher without any loss of accuracy. However, it is important to compute the energy density properly along the edges of Ifronts to guarantee the correct rate of growth of ionized regions.
To further accelerate this computation, we use adaptive time stepping to put higher time resolution on the 3D cells in the vicinity of Ifronts. Since wavefronts cannot propagate further than one grid zone during one time step, only those cells which just experienced large change in and their immediate adjacent neighbours need a proper update of . Depending on the width of Ifronts, typically only a few percent of all 3D cells require the new, exact value of .
One advantage of the use of angleaveraged moments of radiation is that the advection mechanism is essentially reduced to 3D, and it is relatively straightforward to implement the multidimensional conservation scheme for the linear advection part of the moment equations. Then one could use a much denser spatial grid for the solution of the moment equations, and a relatively course grid for the angular reconstruction of the intensity of radiation via ray tracing. In practice, however, we have found that the mismatch between the spatial resolutions of the moment solver and of the ray tracing usually leads to numerical instabilities. In what follows, we consider ray tracing only as a standalone solver.
Another – perhaps, a better – way of coupling angular and spatial variations of the intensity may be an extension of the Spherical Harmonics Discrete Ordinate Method [Evans 1998]. For steadystate transfer problems, instead of storing the radiation field, this method keeps track of the source function as a spherical harmonic series at each point. Although the direct implementation of this technique for timedependent problems is probably not realistic, due to the lookback time (i.e. the finite speed of light propagation), the spherical harmonic representation of the radiation field might require less storage and might result in smoother angular dependence as compared with a pure ray tracing approach.
4 Tests
For all of our test runs, except the study of the shadow behind a neutral clump in Sec. 4.4, we set up a numerical grid with dimensions . The angular resolution was chosen to match the equivalent resolution of data points for 5D advection. There are rays passing in the immediate neighbourhood (th of the total computational volume) of each 3D grid cell, each ray containing grid nodes. Thus, in 5D we obtain the equivalent resolution of data points.
Also, for all simulations we take a (comoving) cosmological volume, scaled down to . All densities (from lowdensity optically thin ambient gas to dense clumps) fall in the range , assuming a Hubble constant of . This low density contrast is chosen to cover the typical range encountered in cosmological hydrodynamics, and is ideally suited to demonstrate transient features during patchy ionization. The absorption coefficients are taken to mimic complete selfshielding (assumed to be ) at neutral hydrogen column densities higher than , except where we probe different regimes, specifically in Sec. 4.3 where corresponds to , and Sec. 4.4 where corresponds to .
4.1 An isolated spherical expanding Ifront
A simple problem that would test the ability of the method to track Ifronts properly is that of a single, isolated Strömgren sphere expanding around a source of ionizing radiation (see, e.g., Abel et al. 1998). One difficulty of the current approach is that rays are drawn in a way to cover the whole computational volume uniformly and isotropically. For a single point source of radiation this would mean that only a small number of rays pass through its neighbourhood. Although intensities along individual rays are strictly conserved, there is no guarantee that the energy density has the right value far from a source of a specified luminosity.
At time we turn on a point source, which starts to blow an expanding H ii bubble around it. Since our algorithm conserves 1D intensities exactly, independently of spatial or angular resolution, the Ifront must propagate at the right speed, which can be obtained analytically by equating the number of direct ionizing photons to the flux of neutral atoms crossing the front. In the absence of radiative recombinations (), the H ii bubble grows indefinitely with the radius
where . Parameters used for this calculation are the diffuse neutral gas density (at ) and the central source luminosity (all emitted in photons above the hydrogen Lyman limit). The comparison between the numerical speed of the Ifront and the exact solution from eq. (4.1) is plotted in Fig. 1. The difference between the two always stays within one grid zone.
4.2 An isolated Strömgren sphere in the presence of a density gradient
For this test we put a point source of radiation into a density gradient along one of the principal axes of the cube. In the absence of diffuse radiation from H ii regions () the only ionizing photons come directly from the source in the centre, in which case the shape of the ionized bubble would be a simple superposition of Strömgren spheres with radii varying with the azimuthal angle and given by the classical solution [Spitzer 1968]
(20) 
where is the photon production rate of the central source. For an exponential density gradient along the yaxis
(21) 
(, being the hydrogen densities on the opposite faces of the cube), the equilibrium Strömgren radius is given by a simple equation
(22) 
where
We take the physical densities on the opposite faces of the volume to be and , and the luminosity is . Since we do not solve any realistic atomic rate equations in this paper, we take the value of the total hydrogen recombination rate from Hummer (1994) for some fiducial temperature (). In Fig. 2 we plot a time sequence of models, for ionization by a central source, with no scattering of Lyman photons (i.e. ). The numerical solution at appears to be very close to the exact one for an equilibrium Strömgren sphere. The sharp transition layer between the ionized and the neutral regions in the high optical depth regime indicates that, indeed, the scheme introduces very little numerical diffusion even when extended to 3D.
4.3 Ionization in the presence of a UV background
The uniform coverage of the whole volume with rays implies that extended sources of radiation will be represented statistically much better than point sources. A simple test mimicking the evolution of dense clouds in the presence of ionizing radiation is to enclose the computational region in an isotropic bath of photons. The simplest way to accomplish this is just to set up a uniform, isotropically glowing boundary at the edges of the cube at . An effective demonstration of timedependent ray tracing would be its ability to deal with any distribution of state variables within the simulation volume. For this test, we set up a density condensation shaped as the acronym for ‘radiation hydrodynamics’ (RHD), with a density times that of the ambient homogeneous medium. The ambient medium has a constant density of , and the energy density of the background radiation is . Fig. 4 shows the result of this run. Most of the lowdensity environment is ionized on the radiation propagation timescale. It takes somewhat longer for ionizing photons to penetrate into the dense regions. Whether these regions can be ionized on a timescale of interest, depends on the ratio of the recombination timescale to the flux of background radiation. One can easily see ionization ‘eating in’ to the neutral zone, e.g. in the disappearance of the serifs on the letters at late times. Note that the width of the ionization fronts does not usually exceed one grid zone (Fig. 5).
4.4 Diffuse radiation from H ii regions: shadows behind neutral clouds
Part of the ionizing radiation at high redshifts comes in the form of hydrogen Lyman continuum photons from recombinations in diffuse ionized regions. The following test, simulating the formation of shadow regions behind dense clouds at the resolution , was adapted from Canto, Steffen & Shapiro (1998). A neutral clump of radius (comoving) is being illuminated by a parallel flux ) from one side; here is Planck’s constant. A shadow behind the clump is being photoionized by secondary recombination photons from the surrounding H ii region (Fig. 6) of physical density . Neglecting hydrodynamical effects, the width of the shadow region can be estimated using a simple twostream approximation [Canto et al. 1998]: of stellar ionizing photons (just above
(23) 
and the dimensionless parameter is defined as
(24) 
For , recombination Lyman continuum photons from the illuminated region will eventually photoionize the shadow completely. For , radiative losses through lowenergy cascade recombination photons will stop the Ifront, forming a neutral cylinder behind the dense clump. Strictly speaking, equations (23)–(24) are valid only for a shadow completely photoionized by secondary photons, and should be viewed as an approximation to Ifronts driven by scattering.
For this run we modify boundary conditions to include recombination photons originating outside the box. Each of the rays – starting on any face except the upper side of the volume (which goes through the neutral shadow) – carries the additional intensity of , where the angular brackets denote the average throughout the currently photoionized gas inside the box.
In Fig. 7 we plot the radius of the shadow neutral region as a function of in our 3D numerical models. To play with the size of the neutral core, we fix the total recombination coefficient but we vary the portion of recombinations into the ground state (which would produce more Lyman continuum photons capable of ionizing the medium). The width of the Ifront driven by secondary photons depends on the assumed opacity (the optical depth of corresponding to the column density of ). This low opacity was chosen to reach equilibrium quicker – equilibrium itself does not depend on the opacity – but it cannot be too low otherwise there will be an unnecessary spread of the transition zone in the Ifront over many grid cells.
We find a remarkably good agreement between the results of our models and the analytic solution, taking into account that ionization of the shadow is due to scattering in the medium.
4.5 Diffuse radiation from H ii regions: ionization of a central void
To demonstrate the ability of our scheme to handle scattering in more complicated situations, we also set up a model with ionization of a central lowdensity void by secondary, recombination photons. The void region is surrounded by two nested cubes with opposite faces open. The walls of the cubes are set to be much denser than the rest of the medium, to screen completely the central void from direct ionizing photons, and the ionizing UV flux is introduced at all faces of the computational volume. The total hydrogen recombination coefficient is again taken for the temperature . Similar to Sec. 4.4, we vary the amount of scattering in the medium by changing the fraction of atoms recombining into the ground state. In reality at the value [Hummer 1994] gives a solution in between our extreme values of .
Similar to the test problem of Section 4.3, if , then the medium will be ionized completely, since there is a constant flux of primordial ionizing photons. The speed of ionization depends on the values of , , and . Note, however, that if is too high, the Ifront will be very slow, since a large portion of the original ionizing photons are scattered back. On the other hand, if is too low, the Ifront will propagate much faster in those regions where ionization is driven by primordial photons, but in shadowed regions there will be too few recombination photons. Thus, it seems that the speed of ionization of the central void will be the highest at some intermediate .
In Fig. 8 we demonstrate ionization of the void region for models with complete scattering and with no scattering at all. The parameters used for this model are . and the diffuse neutral gas density
As expected, for the noscattering model () the central region remains neutral, since there is no direct path for the ionizing photons. However, for the model which includes scattering (), at least part of the central region becomes ionized. This demonstrates that our scheme is perfectly capable of dealing with rescattering of the ionizing photons. Beyond this simple test case, there are many astrophysical situations where progress can be made via numerical radiative transfer. For example, analytic solutions are often used, which are steadystate, and which assume a sharp boundary between the neutral and ionized zones. Using our numerical techniques it should be possible to follow general systems with complex density inhomogeneities as well as regions of partial ionization.
5 Conclusions
In this paper we have demonstrated that, with existing desktop hardware, it is possible to model cosmological inhomogeneous reionization on a lightcrossing time in three spatial dimensions. Since the photoionization timescale in the low optical depth regime () is of order of the lightcrossing time , explicit advection might be a faster method in covering at least these regions. Compared to the elliptictype solvers on the fluidflow timescale or the timescale of atomic processes, explicit radiative advection produces very accurate results without the need to solve a large system of coupled nonlinear elliptic equations. The computing requirements with explicit advection grow linearly with the inclusion of new atomic and molecular rate equations, which is certainly not the case for quasistatic solvers (although it is feasible that the development of multigrid techniques for elliptic equations might actually approach similar scaling).
Using eq. (1) we can see that the entire history of reionization can be modeled with – timesteps (depending on the required resolution), which makes explicit advection an attractive choice for these calculations. However, the efficiency of the explicit radiative solver has still to be explored. Future work should include a detailed comparison between explicit advection and implicit reconstruction (through an elliptic solver), to demonstrate which method works best for calculating inhomogeneous reionization.
As we have demonstrated here, for certain problems, including the propagation of supersonic Ifronts, the Courant condition does not seem to impose prohibitively small timesteps. In this case the biggest challenge is to accurately describe anisotropies in the radiation field, i.e. to solve for inhomogeneous advection in the 5D phase space, in the presence of nonuniform sources and sinks of radiation. Strictly speaking, the storage of one variable at, say, data points requires about GB of memory, which stretches the capabilities of topend desktop workstations. One attractive possibility for future exploration is to directly solve the monochromatic photon Boltzmann equation in 5D. To demonstrate the feasibility of the numerical solution, however, among different methods, we here chose to concentrate on simple ray tracing at the speed of light. The numerical approach we have used is completely conservative and produces very little numerical dissipation.
We want to conclude that despite the high dimensionality of the problem, with a reasonable expenditure of computational resources, of the type available today, it is possible to numerically model many different aspects of the full 3D radiative transfer problem. Furthermore we feel that the methods described here represent a significant and realizable step towards the goal of full cosmological RHD.
Acknowledgments
We wish to thank Taishi Nakamoto for providing us with the results of reionization models from the University of Tsukuba ahead of publication. A.R. would like to thank Jason R. Auman for numerous enlightening discussions, and Gregory G. Fahlman for constant encouragement on this project, as well as Randall J. LeVeque for help with numerical methods for multidimensional conservation laws. This work was supported by the Natural Sciences and Engineering Research Council of Canada.
References
 [Abel et al. 1997] Abel T., Anninos P., Zhang Y., Norman M.L., 1997, New Astronomy, 2, 181
 [Abel et al. 1998] Abel T., Norman M.L., Madau P., 1998, astroph/9812151
 [Anninos et al. 1997] Anninos P., Zhang Y., Abel T., Norman M.L., 1997, New Astronomy, 2, 209
 [Auer & Mihalas 1970] Auer L.H., Mihalas D., 1970, MNRAS, 149, 65
 [Canto et al. 1998] Canto J., Steffen W., Shapiro P.R., 1998, ApJ, 502, 695
 [Couchman & Rees 1986] Couchman H.M.P., Rees M.J., 1986, MNRAS, 221, 53
 [Evans 1998] Evans F., 1998, J. Atmospheric Sciences, 55, 429
 [Gnedin & Ostriker 1997] Gnedin N.Y., Ostriker J.P., 1997, ApJ, 486, 581
 [Haiman & Loeb 1997] Haiman Z., Loeb A., 1997, ApJ, 483, 21
 [Haiman & Loeb 1998] Haiman Z., Loeb A., preprint astroph/9710208
 [Hummer 1994] Hummer D.G., 1994, MNRAS, 268, 109
 [Katz et al. 1996] Katz N., Weinberg D.H., Hernquist L., MiraldaEscudé J., 1996, ApJ, 457, L57
 [Kepner et al. 1997] Kepner J.V., Babul A., Spergel D.N., 1997, ApJ, 487, 61
 [Madau et al. 1997] Madau P., Meiksin A., Rees M.J., 1997, ApJ, 475, 429
 [Madau 1998] Madau P., 1998, preprint astroph/9807200
 [Meiksin 1994] Meiksin A., 1994, ApJ, 431, 109
 [Mihalas & Mihalas 1984] Mihalas D., Mihalas B., 1984, Foundations of Radiation Hydrodynamics. Oxford University Press, New York
 [MiraldaEscudé 1998] MiraldaEscudé J., Haehnelt M., Rees M.J., 1998, astroph/9812306
 [Navarro & Steinmetz 1997] Navarro J.F., Steinmetz M., 1997, ApJ, 478, 13
 [Norman et al. 1998] Norman M.L., Paschos P., Abel T., 1998, astroph/9807282
 [Park & Hong 1998] Park Y.S., Hong S.S., 1998, ApJ, 494, 605
 [Shapiro 1986] Shapiro P.R., 1986, PASP, 98, 1014
 [Shapiro et al. 1998] Shapiro P.R., Raga A.C., Mellema G., 1998, preprint astroph/9804117
 [Spitzer 1968] Spitzer L., Jr., 1968, Diffuse matter in space. Interscience publishers, New York
 [Stone & Mihalas 1992] Stone J.M., Mihalas D., 1992, JComputPhys, 100, 402
 [Stone et al. 1992] Stone J.M., Mihalas D., Norman M. L., 1992, ApJS 80, 819
 [Tajiri & Umemura 1998] Tajiri Y., Umemura M., 1998, preprint astroph/9806046
 [Umemura et al. 1998] Umemura M., Nakamoto T., Susa H., 1998, in Miyama S. M., Shibata K., eds, Numerical Astrophysics 1998. Kluwer, in press
 [Zhang et al. 1998] Zhang, Y., Meiksin, A., Anninos, P., Norman, M.L., 1998, ApJ, 495, 63