 |
| Fig.1: Hydrogen reionization must have occurred somewhere in the redshift interval z = 40 - 5. |
Understanding details of reionization would help us to study the role of feedback mechanisms on subsequent structure formation, as well as learn something about the first generation of objects by studying their imprint on the IGM. By the time the first halos collapsed, the gaseous component of the Universe must had been very clumpy. Unfortunately, propagation of ionization fronts into an inhomogeneous medium is an unsolved problem. Naively, one would expect that ionizing radiation would first stream into low-density voids where absorption by neutral hydrogen is the lowest. However, sources which cause reionization most likely reside inside virialized dark matter halos within highly overdense knots at the intersection of large-scale cosmological filaments. Assuming that most ionizing photons are not very energetic, that is reionization is not caused by soft X-rays from quasars or star forming regions, reionization is bound to demonstrate a very patchy character. A question then is what is the topology of ionized regions given the spectra, luminosities and the spatial distribution of ionizing sources. The answer is not trivial since it depends on the details of radiative transfer (RT) in a highly inhomogeneous IGM.
Most numerical studies in the field either rely on some simplified description of RT (Gnedin 1999), or simply neglect density inhomogeneities around each source (Ciardi et al. 2000a) assuming that photodissociation and photoionization regions around each source are spherically symmetric. A fast way to compute a radius of the ionization front propagating from a point source into a non-homogeneous medium is patch together solutions at different azimuthal angles assuming that ionizing photons propagate just along radial vectors originating at the source. However, ionizing photons can come from recombinations in the medium, for one thing, and the recombination rate is especially high in dense regions. One obvious example where pure radial RT would fail is ionization of shadows formed behind dense regions (Fig. 2).
 |
| Fig.2: Cross-sections of the dense clump (shaded) and the neutral shadow region behind it, for different fractions q = a1 / atot of recombinations to the ground state. The clump is being illuminated with a parallel flux of direct ionizing photons entering through one side of the box (marked with arrows), and the shadow region behind it is being photoionized by recombination Lyman continuum photons from the surrounding HII gas. Each model has been evolved to an equilibrium state, corresponding to many passages of the wavefront across the computational region. The nine contour levels correspond to ionization fractions of xe = 0.1, 0.2, ..., 0.9. |
Numerical scheme
Recently, there has been a surge of interest in constructing an efficient 3D solver for cosmological RT (Abel et al. 1998, Nakamoto & Umemura 1998, Gnedin 1999, Ciardi et al. 2000b). At UBC, over the last few years, we have developed a numerical scheme for 3D, time-dependent RT in a clumpy medium where every cell can be a source of radiation. This is a very interesting problem because in 3D each point can affect every other point via the radiation field, like in gravity except that radiation is also attenuated along the line connecting these two points, due to absorption and scattering. In addition, the optical properties of the medium can change very quickly (chemical reaction rates are sometimes much faster than the light crossing time through the computational cell). Finally, the intensity of radiation depends on frequency. Therefore, from the computational standpoint, at the same spatial resolution RT is a much more challenging task than solving an N-body problem (Fig. 3).
Fig.3: The radiative transfer equation. |
The code we have developed is based on a fast explicit advection of radiative intensities in the multi-dimensional phase space (3 spatial coordinates, 2 angles in the direction of photon propagation, and frequency resulting in 6D transport) and an implicit solution of non-equilibrium rate equations independently at each point in space. Main requirements for this method were:
- fast advection in 5D (without frequency) or 6D (with frequency) phase space,
- high spatial resolution of 3D fronts (little numerical dispersion),
- correct front speed in 3D, and
- reasonable angular resolution for a given spatial mesh.
There is a large variety of methods for explicit multidimensional advection (for a brief listing, see Fig. 4). Our scheme implements one of the most straightforward techniques, time-dependent transport of radiation intensities along a large number of rays in 3D (Fig. 5). Unlike in other areas of astrophysics, explicit transport of radiation at the speed of light does not necessarily yield a prohibitively large number of time steps. For instance, the evolution of a cosmological volume from z = 20 to z = 3 with radiation at the resolution 20 kpc (comoving) would require one to take ~ 105 time steps. Frequency dependence is taken care of with the multigroup approach, that is, the RT equation is solved numerically to find integrated radiative quantities inside individual frequency groups (Fig. 6).
 |
| Fig.4: Methods for explicit multidimensional advection. |
 |
| Fig.5: The volume is covered uniformly and isotropically with rays. Radiative intensities are being transported explicitly along each ray at the speed of light. At every time step these intensities are projected to the 3D Cartesian grid to reconstruct the energy density at each point which is then used to evolve chemistry rate equations. |
 |
| Fig.6: Schematic representation of the three frequency groups used in our calculations, along with some other radiation thresholds. |
Uncertainties
Apart from obvious difficulties in implementing numerical RT, there are a number of major uncertainties which we face in modelling reionization.
- We do not know the star formation efficiency, that is, what fraction of baryons in a collapsed halo actually goes into stars.
- We can only guess the initial mass function of the primordial stellar population. By far, the most efficient feedback would come from most massive stars. Among other physical processes, the stellar mass scale is determined by the cooling properties of contracting protostellar clouds. There are a number of mechanisms which can lead to fragmentation on scales smaller than 100 solar masses (e.g., 3-body formation of molecular hydrogen H2, or rotation). On the other hand, Norman et al. (2000) recently performed 3D simulations of the first structure formation using adaptive meshes to cover the dynamical range 3 × 107, taking a comoving box of 128 kpc on a side and zooming in on dense regions to evolve them to nearly stellar densities. Despite very high final densities (1012 cm-3) achieved in their simulation, they conclude that there is no evidence of fragmentation on scales smaller than ~ 100 solar masses.
- The effect of the first generation of stars on their surroundings depends on the escape fraction of ionizing photons. If most of the Lyman continuum radiation is absorbed within the galaxy, then simply many more stars would be needed to reionize the IGM.
Quasar reionization
Quasar reionization is distinctively different from stellar reionization. With quasars, one usually has just a few bright sources inside a cosmological volume. Their luminosities are so high that - at least at the resolution level of the simulations - they ionize everything in their immediate neighbourhood producing very large HII regions (Fig. 7). Their spectra are much flatter than those from stars, hence, ionization fronts are wider since there are many more high-energy photons travelling ahead of the front (Fig. 8).
| Fig.7: This figure shows a cross-section through a 3D hydrogen ionization front propagating into the neutral IGM from a virialized halo hosting a mini-quasar around z = 12, at four different output times (0.048, 0.53, 1.07 and 2.15 Myrs after the quasar had been turned on, respectively). The quasar has a luminosity equivalent of 1054 hydrogen Lyman limit photons s-1 distributed over the entire spectrum above 13.6 eV as a power law with the spectral index n = 1.5. The simulation employs the full 3D radiative transfer in three frequency groups (for both direct ionizing and diffuse photons) combined with the solution of the rate equations for nine chemical species, at the spatial resolution 643. The size of the volume is 2.4 Mpc (comoving). To compare the ionization front speed to the speed of light, the light travel time across the volume is 0.6 Myrs. |
|
|
| Fig.8a: Same as Fig. 7, but for singly ionized helium. Note that the HeII region is a relatively thin shell bounded by neutral helium in the ambient medium and by doubly ionized helium closer to the quasar. | Fig.8b: Same model as in Fig. 8a, except that now the spectral index is changed to n = 5, to mimic much softer radiation. The singly ionized helium region occupies a considerably larger volume sice there are fewer photons capable of doubly ionizing helium. The bolometric luminosity of the point source is the same as in Fig. 8a. |
Stellar reionization
On the other hand, if reionization is caused by stellar radiation, there are many more not so luminous sources hidden inside virialized dark matter halos. It is likely that stars did not ionize their host environment at once, and ionization fronts might have first propagated into the low-density IGM via narrow ionized channels (Gnedin 1999). Fig. 9 shows an example of such calculation. There are about 800 stellar sources (stellar groups) in the volume. Since in our calculations the angular-dependent radiation field is computed at each point in 3D, the models are very-resolution limited at the moment (643 spatial × 82 angular). At the moment we are preparing a run with a much higher angular resolution. Of course, the goal is to eventually connect every two points in the volume.
 |
| Fig.9: This diagram shows the position of the hydrogen ionization front for stellar reionization at four consecutive redshifts (z = 11.6, 11.2, 10.8 and 10.3), for a starburst at z = 12. The resolution is 643 (spatial) × 82 (angular). The star formation rate is normalized to that in the highest resolution run in Gnedin (1999), stellar population synthesis is taken from Ciardi (2000a). |
Summary
Most models of galaxy formation only considered the force of gravity plus some gas physics, mostly because until now dealing with realistic radiative transfer has been well beyond our computational capabilities. Methods for high-resolution numerical RT are being currently developed, to solve transfer problems on both parallel computers and stand-alone workstations. One application of cosmological radiative transfer is in understanding details of how the Universe went from a neutral state to being ionized by the light from the first stars and galaxies as we observe it today.
References
Abel, T., Norman, M. L., Madau, P., 1998, astro-ph/9812151
Ciardi, B., Ferrara A., Governato, F., Jenkins, A., 2000a, MNRAS, 314, 611
Ciardi, B., Ferrara, A., Marri, S., Raimondo, G., 2000b, astro-ph/0005181
Gnedin, N. Y., 1999, astro-ph/9909383
Nakamoto T., Umemura M., 1998, presentation at the Ann. Meeting of the Astron. Soc. of Japan
Norman, M. L., Abel, T., Bryan, G., 2000, astro-ph/0005246
