Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Calculating the inhomogeneous reionization of the universe Razoumov, Alexei O. 1999

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2000-486982.pdf [ 6.39MB ]
Metadata
JSON: 831-1.0085451.json
JSON-LD: 831-1.0085451-ld.json
RDF/XML (Pretty): 831-1.0085451-rdf.xml
RDF/JSON: 831-1.0085451-rdf.json
Turtle: 831-1.0085451-turtle.txt
N-Triples: 831-1.0085451-rdf-ntriples.txt
Original Record: 831-1.0085451-source.json
Full Text
831-1.0085451-fulltext.txt
Citation
831-1.0085451.ris

Full Text

C A L C U L A T I N G T H E I N H O M O G E N E O U S R E I O N I Z A T I O N O F T H E U N I V E R S E By Alexei O. Razoumov B. Sc., M. Sc. (Astronomy) Moscow State University A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E STUDIES PHYS ICS & A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UN IVERS ITY O F BR IT ISH C O L U M B I A October 1999 © Alexei O. Razoumov, 1999 1 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Physics & Astronomy The University of British Columbia 129-2219 M a i n M a l l Vancouver, Canada V6T 1Z4 Date: bJoviJUc {S 4333 * — Abstract A numerical scheme for the solution of the three-dimensional, frequency- and time-dependent ra-diative transfer equation with variable optical depth is developed for modelling the reionization of the Universe. Until now, the main difficulty in simulating the inhomogeneous reionization has been the treatment of cosmological radiative transfer. The proposed approach is drastically different from previous studies, which either resorted to a very simplified, parametric treatment of radiative transfer, or relied on one-dimensional models. The algorithm presented here is based on explicit multidimensional advection of wavefronts at the speed of light, combined with the implicit solution of the local chemical rate equations separately at each point. I implement the ray-tracing version of this algorithm on a desktop workstation and check its performance on a wide variety of test problems, showing that explicit advection at the speed of light is an attractive choice for simulation of astrophysical ionization fronts, particularly when one is interested in covering a wide range of optical depths within a 3D clumpy medium. This scheme is then applied to the calculation of time-dependent, multi-frequency radiative transfer during the epoch of first object formation in the Universe. In a series of models, the 2.5 Mpc (comoving) simulation volume is evolved between the redshifts of z = 15 and z = 10 for different scenarios of star formation and quasar activity. The highest numerical resolution employed is 643 (spatial) x 102 (angular) x 3 (frequency), and at each point in space I calculate various stages of hydrogen and helium ionization accounting for nine chemical species altogether. It is shown that at higher numerical resolution these models of inhomogeneous reionization can be used to predict the observational signatures of the earliest astrophysical objects in the Universe. At present, the calculations are accurate enough to resolve primordial objects to the scale typical of globular clusters, 1O6M0. ii Table of Contents Abstract " List of Tables vii List of Figures viii Acknowledgements x 1 Introduction to Cosmological Radiative Transfer 1 2 Numerical Radiative Transfer in an Inhomogeneous Medium 12 2.1 Radiative Transfer in Moving Media 12 2.2 The Equation of Radiative Transfer in Stationary Material 14 2.3 The Eddington Approximation 15 2.4 Time-Dependent vs. Steady-State Techniques 18 2.5 Formal Solution of the Transfer Equation 22 2.5.1 Long Characteristics 23 2.5.2 Short Characteristics 25 2.6 Time-dependent 3D Moment Solver 26 2.7 Truncating the System of Moment Equations at a Higher Order 28 2.8 Solution of the Multi-dimensional Photon Boltzmann Equation 30 2.8.1 Normal Flux Update 32 2.8.2 First-order Transverse Fluxes in 3D 33 2.8.3 First-order Corner Fluxes in 3D 36 iii 2.8.4 Second-order Updates 37 2.8.5 Summary 38 3 Time-Dependent Ray Tracing 39 3.1 Descr ipt ion of the M e t h o d 39 3.1.1 A Uni fo rm and Isotropic G r i d of Rays 40 3.1.2 Angula r Reconstruction 44 3.1.3 T h e A l g o r i t h m 46 3.1.4 Adap t ive T i m e Stepping 49 3.2 Tests 49 3.2.1 L o c a l Chemis t ry Equations 49 3.2.2 Test Parameters and Resolution 52 3.2.3 A n Isolated Spherical Expand ing I-front 53 3.2.4 A n Isolated Stromgren Sphere in the Presence of a Densi ty Gradient . . . 54 3.2.5 Ionization in the Presence of a U V Background 56 3.2.6 Diffuse Radia t ion from H n Regions: I. Shadows behind Neut ra l Clouds . 56 3.2.7 Diffuse Radia t ion from H I I Regions: II. Ionization of a Cent ra l V o i d . . . 58 4 Including the Physics of Matter-Radiation Interaction 68 4.1 Frequency Dependence: Mul t i g roup Transfer 68 4.2 Rate Equat ions 70 4.2.1 Number Density of Absorpt ion-Induced Processes 74 4.3 Absorp t ion Coefficients 75 4.4 Photoioniza t ion Rates i n the L o w Opt i ca l Dep th Regime 76 4.5 Photoioniza t ion Rates: the General Case 76 4.6 Emiss iv i ty Inside Each Frequency Group 77 4.6.1 Hydrogenic Ions 77 iv 4.6.2 He l ium-Like Ions 78 4.6.3 Mat te r Energy Conservation 79 4.7 Testing the Chemis t ry M o d e l vs. Cosmic Recombinat ion at z ~ 1000 80 5 Modelling Inhomogeneous Reionization 84 5.1 The Density Dis t r ibu t ion 84 5.2 Temperature of the I G M before Reheating 85 5.3 Preheat ing 86 5.4 T i m e Stepping wi th Variable Redshift • 88 5.5 M o d e l l i n g the F i rs t Objects 88 5.5.1 M o d e l A : Reionizat ion by Shor t -Lived Mini -Quasars 90 5.5.2 M o d e l B : Reionizat ion by Br ight Quasars 93 '5.5.3 M o d e l C : Reionizat ion by an Ex te rna l U V F i e l d 94 5.5.4 Cr i t e r i a for Reionizat ion by Popula t ion I I I Stars 95 5.6 Observational Properties of Reionizat ion 100 5.7 A t o m i c Hydrogen 21cm Emiss ion 101 5.7.1 21-cm Emiss ion and Absorp t ion 102 5.7.2 High-z Absorp t ion Features 104 5.7.3 High-z Emiss ion Features 104 5.8 Discussion and Conclusions 108 6 Concluding Remarks 128 Appendices 131 A High-Resolution Schemes for One-Dimensional Hyperbolic Systems 131 B An Unsplit Method for 3D Moment Transfer 134 v C Observing Programs to Detect the Epoch of Reionization 137 D Glossary 1 3 9 References 141 v i List of Tables 5.1 Number of vi r ia l ized halos vs. redshift 89 5.2 Reioniza t ion models 90 C . l Observational probes of the epoch of reionization 137 C.2 Observat ional probes of the epoch of reionization (continued) 138 v i i List of Figures 1.1 E p o c h of reionization 3 2.1 Methods to solve the R T E numerically 19 2.2 T h e effect of poor angular resolution 24 2.3 Smoothing wi th poor angular resolution 25 2.4 L o n g Characterist ics 26 2.5 Short Characteristics 27 2.6 Transverse propagation in 2D ' 35 2.7 Transverse and corner propagation i n 3D 36 3.1 Rays shooting through the base grid i n the direction 41 3.2 Piece-wise linear angular variat ion of the intensity 43 3.3 Numer ica l speed of a spherically symmetric I-front 54 3.4 I-front expanding around a point source into a non-homogeneous med ium . . . . 60 3.5 Nested 3D rectangular and angular grids 61 3.6 Ionizat ion of the ' R H D ' acronym: 3D isosurfaces 62 3.7 Ionizat ion of the ' R H D ' acronym: contour plot 63 3.8 Shadows behind neutral clouds 64 3.9 The w i d t h of the neutral shadow behind a dense c lump 65 3.10 Geometry of the central void region 66 3.11 Ionizat ion of the void region 67 4.1 Frequency groups 69 v i i i 4.2 T h e test of the photochemistry model against cosmic recombination 82 5.1 The adopted density field at eight different redshifts I l l 5.2 Hydrogen I-front i n M o d e l A 112 5.3 H e l i u m I-front in M o d e l A 113 5.4 Temperature of the I G M i n M o d e l A 114 5.5 Hydrogen I-front i n M o d e l B 115 5.6 H e l i u m I-front in M o d e l B 116 5.7 Color map of the hel ium I-front i n M o d e l B 117 5.8 T h e energy density of the radiat ion field i n three frequency groups vs. fractional ionizat ion i n M o d e l B 118 5.9 Hydrogen I-front in M o d e l C 119 5.10 Hydrogen I-front for a different cross-section through M o d e l C 120 5.11 Hydrogen I-front i n M o d e l C at a sl ightly lower redshift 121 5.12 3D dis t r ibut ion of neutral hydrogen i n M o d e l C 122 5.13 3D hydrogen I-front in M o d e l C at different spatial resolutions 123 5.14 Sp in temperature at high redshift 124 5.15 21-cm absorption or emission at high redshifts 125 5.16 21-cm radiat ion skymap from M o d e l B 126 5.17 21-cm radiat ion skymap from M o d e l B (continued) 127 ix Acknowledgements I would like to thank Jason A u m a n , Gregory Fahlman and Douglas Scott for constant encour-agement on all my projects. I am part icular ly grateful to my research supervisor Douglas Scott for a great deal of valuable insight and for his seemingly limitless patience. I would also like to thank Rodr igo Ibata for numerous enjoyable discussions and for encouraging me to try out new ideas, and B e l l Chen for a very friendly atmosphere i n our office. Special thanks go to my parents and my sister for their constant love and emotional support. x Chapter 1 Introduction to Cosmological Radiative Transfer One of the most fundamental problems in astronomy is the nature of the very first generation of astrophysical objects at high redshifts, which formed when the Universe was only 5 — 10% of its present age. T w o facets — the dynamics of dark matter and the thermal properties of the high-redshift gas — shaped the physical evolution at those early epochs. Over the last 15 years the theory of large-scale structure formation has benefited significantly from numerical N - b o d y and gas-dynamical simulations, which successfully describe the emergence of hierarchical structures out of p r imord ia l density fluctuations, as well as the collapse of i nd iv idua l objects. These techniques work quite well when properties of the gas are not coupled strongly to the spat ial gradients i n the radiat ion field; e.g., i n the smal l opt ical depth regime one can approximate the energy balance between matter and radiat ion v i a a combinat ion of photoheating by a uniform background flux and a parametric cooling function associated w i t h this radia t ion (Efstathiou 1992, Mira lda-Escude & Rees 1994, K a t z et al. 1996b, Q u i n n et al. 1996, Navarro & Steinmetz 1997, Haehnelt & Steinmetz 1998, B r y a n et al. 1999), or even derive an effective equation of state of the intergalactic medium ( I G M ) i n the presence of a uniform ioniz ing background ( H u i & Gned in 1997, Schaye et al. 1999, R i c o t t i et al. 1999). T h i s approximat ion is sufficient i n the low-density, completely ionized I G M which is exposed to the ba th of ultraviolet ( U V ) photons at lower redshifts, but it breaks down i n the cores of collapsed halos or i n the neutral , high-redshift I G M which is s t i l l opaque to ionizing photons. A t present, the majori ty of cosmological gas-dynamical simulations lacks a proper description of radiative transfer (RT) mak ing these models inapplicable to more realistic states of opaque and inhomogeneous I G M . 1 o Chapter 1. Introduction to Cosmological Radiative Transfer 2 A n approach t radi t ional ly popular in cosmology is to solve the radiative transfer equation ( R T E ) i n a space-averaged way (Miralda-Escude & Ostriker 1990, 1992, G i r o u x & Shapiro 1996) assuming that the matter dis t r ibut ion is homogeneous, the radiat ion field is uniform and isotropic, and neglecting R T effects w i t h i n ind iv idua l clouds. The mean effect, of clumps is usually smoothed over a uniformly distr ibuted ambient gas. These techniques can describe the overall R T i n the I G M reasonably well i f the effective opt ical depth of the Universe is smaller than unity; however, even then, they fail woefully inside dense structures. The simplest approximat ion which takes into account the effects of R T wi th spat ial ly vari-able opt ical depth is probably that of an infinite slab (Haiman et al. 1996, Haardt & M a d a u 1996) or a spherically symmetric cloud (Kepner et al. 1997, Ta j i r i & U m e m u r a 1998, A b e l & Haehnelt 1999, Loeb & R y b i c k i 1999, H a i m a n et al. 1999) irradiated by a uniform background field. T h i s model can be a good approximat ion to a vi r ia l ized halo immersed into a U V back-ground after the epoch of reionization (Kepner et al. 1997, Ta j i r i & Umemura 1998), or to a smal l ( 1 O 4 _ 5 M 0 ) cloud i l luminated by low-energy (below 13.6 e V ) photons before reionization (Haiman et al. 1999), but it does not account accurately for the dumpiness of the I G M and for the spatial d is t r ibut ion of sources. A s we shal l see, there are reasons to believe that the radiat ion field has an important effect on the thermal state of interstellar and intergalactic medium, which can be accounted for only w i t h a detailed model for cosmological R T . W i t h i n currently popular cosmological scenarios, the Universe recombined around redshift z ~ 1000, as i t was expanding and cool ing down after the p r imord ia l hot phase. A t some stage, after decoupling from radiat ion at a redshift of few hundred, the gas began to collapse slowly under its own gravity forming the very first bound objects (Peebles 1993, p. 176). It is believed that the light from these objects later reionized most of the gaseous component of the Universe (Fig . 1.1). Observations of the cosmic microwave background ( C M B ) fluctuations on the sky, on the one hand, and of absorption features i n the spectra of high-redshift quasars, on the other hand, place the epoch of reionizat ion somewhere Chapter 1. Introduction to Cosmological Radiative Transfer i n the redshift interval 6 & z <^  40 (Haiman & K n o x 1999). 3 u neutral IGM first objects fc fc c, ionized hydrogen ^ fc fc growth of gravitational instabilities hydrogen reionization reheating quasars galaxies -1000 (recombination) -300 (decoupling) 35-40 -6 0 redshift Figure 1.1: A schematic diagram showing the epoch of reionization. The best upper l imi t on the reionization redshift comes from the analysis of C M B anisotropies. The scattering of C M B photons off free electrons appearing i n abundance dur ing reionization has several important effects on the C M B (see, e.g., Peebles 1993, p. 581, also Tegmark et al. 1994). Mos t notably, it changes the paths of the photons, suppressing pr imary C M B fluctu-ations on degree and sub-degree scales. Besides that, free electrons Doppler shift C M B photons leading to secondary anisotropies. F r o m a compila t ion of a l l existing data of C M B anisotropies, Griffiths et al. (1998) recently derived a model-dependent upper value of the redshift of reion-izat ion, w i t h the best fit around 2reion ~ 35 — 40. O n the other hand, the opt ical depth of the neutral I G M to L y a absorption (Haiman & K n o x 1999) is very high, destroying radiat ion from a l l sources beyond the redshift z at observed wavelengths shorter than A L y Q ( l + z). The absence of strong absorption blueward of the L y a line i n the spectra of high-redshift objects suggests that hydrogen is almost completely ionized by z = 5 (Gunn & Peterson 1965). Ind iv idua l L y a forest clouds reveal themselves i n absorption at distinct redshifts. However, the opt ical depth of the I G M is T L y a ^ 1, y ie ld ing the average neutral hydrogen fraction <> 10~ 5 . Recent discoveries of L y a emit t ing galaxies at even higher Chapter 1. Introduction to Cosmological Radiative Transfer 4 redshifts (z = 6.68, Chen et al. 1999; z = 5.74, H u et al. 1999; z = 5.60, Weymarm et al. 1998) move the lower l imi t for the epoch of reionization even further back i n time. Patchy reheating and reionization are interesting for a number of reasons. It is l ikely that observational data on the earliest luminous objects i n the Universe wi l l first come from the s tudy of the thermal state and the dis t r ibut ion of the inhomogeneous I G M around these objects (e.g., see Tables C . l - C.2 in A p p e n d i x C for a list of observational programs). A number of fundamental questions could be addressed i n a comparison between the predicted and observable features of this high-redshift gas: • W h a t is the nature of the sources of reionization? Was the Universe ionized by stars or quasars, or both? • How were the first luminous objects distr ibuted i n space? W h e n d id they form? W h a t were their spectra and luminosities? • How long d id reionization last? • W h a t is the connection between reheating and reionization? • Was reheating efficient i n providing the negative feedback on subsequent quasar ac t iv i ty or star formation? • How does c lumping of the I G M affect the timescale of reionization? W h a t is the result ing cooling function? • How large were inhomogeneities at the t ime of the first structure formation? It is important to understand how many photons are actual ly needed to reionize the U n i -verse. T h e most probable source of ionizing radiat ion is the popula t ion of stars and quasars formed dur ing the gravitat ional collapse of the earliest structures. It has been claimed (Shapiro & G i r o u x 1987, Mira lda-Escude & Ostriker 1990, Shapiro et al. 1994, M a d a u et al. 1999) that Chapter 1. Introduction to Cosmological Radiative Transfer 5 since the observed populat ion of quasars i n both opt ical and radio wavelengths declines steadily beyond z ~ 3 (Warren et al. 1994, Pe i 1995, Schmidt et al. 1995), the U V radia t ion from quasars alone is not sufficient for reionization. However, i n the past few years, there has been a grow-ing interest i n models which predict many more low-luminosity quasars at ul t ra-high redshifts (z ~ 10 — 15), based on a suggested correlation between the masses of host dark matter halos and central black holes (Haiman &; Loeb 1998, Haehnelt et al. 1998). O n the observational front, recently, some bright quasars were detected to as far as z ~ 5 (Fan et al. 1999, and references therein), w i th a dozen or so objects known at z > 4.5. T h e recently begun Sloan D i g i t a l Sky Survey (SDSS) has the potential to increase this number to a few hundred over the next five years, possibly, wi th the detection of bright objects in the redshift interval 5.5 <> z <^  7. In addi t ion, there might be some evidence for a large popula t ion of low-luminosi ty quasars at h igh redshifts from recent X - r a y data (Miya j i et al. 1998) which do not show any decline i n the quasar number density at z > 2.7 (Haiman & Loeb 1999). O n the other hand, i f reionizat ion is caused by a generation of stars not long before z = 5, the rate of star formation must be similar to or bigger than that inferred from galaxy observations at z = 3 (Madau et al. 1999). A d d i t i o n a l constraints on the efficiency of star formation at higher redshifts come from stud-ies of metal enrichment of the I G M . T h e detection of heavy elements i n L y a absorption sys-tems shows a clear correlation between metallicities and the degree of collapse of these objects. D a m p e d L y a systems (those wi th neutral hydrogen column densities J V ( H i ) ^ 1 0 2 0 c m - 2 ) demonstrate metal abundances i n the range 3 x 10~ 3 < Z/ZQ < 1 0 _ 1 ( L u et al. 1996) indicat-ing that these are young galaxies w i th strong local star formation. L o w column density clouds (N(Hi ) ~ 3 x 1 0 1 4 c m ~ 2 ) which can be associated w i t h the diffuse I G M and which exhibi t metallicit ies of order Z/ZQ ~ 10~ 2 (Cowie et al. 1995), might have been enriched by Popula-tion III stars formed i n vir ia l ized objects w i t h masses 10 5 — 10 7 M Q at much higher redshifts (30 z <^  15, H a i m a n et al. 1996a). However, recent observations of very low column-density clouds (those w i t h 3 x 1 0 1 3 i V ( H l ) <> 1 0 1 4 c m ~ 2 , L u et al. 1998) seem to have ru led out the Chapter 1. Introduction to Cosmological Radiative Transfer 6 possibi l i ty that Popula t ion III stars could have pol luted the entire Universe to a uniform metal-the amount of heavy elements which one would expect to detect i f these elements were produced by the same early generation of stars that reionized the Universe and were later ejected into the low-density I G M is of order Z/Z@ ~ 3 x 10~'3 (or lower, depending on the dumpiness of the gas and star formation efficiency, Gned in & Ostr iker 1997), hence, the observed metall icit ies are more or less consistent w i t h our theoretical insight on Popula t ion III objects ( L u et al. 1998, C i a r d i et al. 1999). Energetically, it would be sufficient to convert only a smal l ( <> 10~ 5 —10~ 3 , H a i m a n & Loeb 1997) fraction of baryons into quasars and stars to reionize the I G M . However, photoionizat ion is balanced by recombinations, the efficiency of which is a sensitive function of the local matter density, thus depending on the redshift and on the dumpiness of the gas. In a homogeneous universe the ratio of recombination t ime to the Hubble t ime is (Haiman & K n o x 1999) that is, recombinations are much more important at higher redshifts (z <; 10). Simi lar ly , the efficiency of recombinations rises sharply i n high-density clouds. There is no reason to think that reionization was an instantaneous, phase-like t ransi t ion. To the contrary, the gaseous component must have been c lumpy enough to allow for the collapse of ind iv idua l objects." Therefore, one should expect that the ionizing radiat ion w i l l first stream into the low-density voids, where absorption by neutral hydrogen is the lowest. One complicat ion arises though. Sources which caused reionization most l ikely reside inside v i r ia l ized dark matter halos w i t h i n highly overdense spherical knots at the intersection of large-scale filaments ( A b e l et al. 1998b). A question is then whether low-density voids or high-density regions get ionized first (Ha iman & K n o x 1999). The answer is not t r iv i a l since it depends on the details of R T i n a highly inhomogeneous I G M . Gned in & Ostriker (1997) calculated reionizat ion numerical ly l ic i ty level of Z/Ze ~ 10~ 2 — instead, the inferred value is closer to Z/ZQ < 10 3 . Note that 3/2 Chapter 1. Introduction to Cosmological Radiative Transfer 7 concluding that c lumping of the gas increases the global recombination rate, mak ing it harder to reionize the Universe. In this scenario, at the beginning, a large fraction of the ioniz ing photons does not travel far from sources, since these are located inside dense structures. It is i n these clumps that the ionization proceeds first, and at some stage it is followed by an almost instantaneous decrease i n the neutral hydrogen fraction i n low-density regions. In this case, the Universe becomes transparent to the ioniz ing cont inuum on a very short timescale, since most of the neutral hydrogen inside dense filaments must have been ionized by then. However, in a more detailed model (albeit w i t h an approximate treatment of the R T ) G n e d i n (1998) concluded that reionization proceeds i n the reverse direction: first, ionizat ion fronts (I-fronts) propagate into voids and quickly ionize the low-density I G M . In the second stage the ioniz ing radia t ion eats into high-density filamentary structures un t i l they get ionized and the universe as a whole becomes transparent. Similar ly. Mira lda-Escude et al. (1998) argue that — despite the dumpiness of the gas — recombinations are not very important globally since the high-density filamentary gas is not ionized un t i l a later time. In this scenario, ioniz ing photons escape sources w i th in dense structures through a smal l solid angle, without causing complete loca l ionizat ion. In fact, the ionizat ion occurs "outside-in". s tar t ing i n voids and gradual ly eating into overdense regions. For this reason, the complete overlap of H I I regions occurs at much later times after reionization. U n t i l then the Universe stays globally opaque to ioniz ing photons. T h e reionizat ion of the gas has three dramatic consequences for the subsequent evolution of the Universe. F i r s t of all, the photoionized gas is much warmer, and is therefore is far less l ikely to collapse into bound structures. In other words, reionization raises the overall Jeans mass. Second, the cooling of a neutral hydrogen/hel ium medium w i t h p r imord ia l abundance i n the absence of any ionizing flux is dominated by line cooling. If the same med ium is submerged i n the b a t h ' o f photoionizing radiat ion, line cooling is far less efficient (the gas is ionized), and cooling is dominated by recombinations (Efstathiou 1992). O n the other hand, objects Chapter 1. Introduction to Cosmological Radiative Transfer 8 which reionize the Universe w i l l also enrich it w i th heavy elements, enhancing cooling in the temperature range 10 4 —10 7 K . F ina l ly , the ionized gas is v i r tua l ly transparent to high-frequency photons streaming from nearby quasars or star forming regions — unlike the neutral med ium which blocks a l l radia t ion wi th frequencies above the L y a line. Unfortunately, the evolution of I-fronts i n a general 3D medium is an unsolved problem. It now seems clear that the full solution requires a detailed treatment of the effects of R T . To com-plicate matters, by the t ime of the first star formation, the small-scale density inhomogeneities have entered the non-linear regime, and the medium was filled w i t h c lumpy structures (Gned in & Ostr iker 1997). In recent years, hydrodynamical simulations of structure formation in a universe dominated by cold dark matter ( C D M ) and often inc luding an ioniz ing background, have been very successful in quantifying the density d is t r ibut ion and the ionizat ion state of the I G M (e.g., C e n et al. 1993, Hernquist et al. 1996, K a t z et al. 1996, Zhang et al. 1998). The degree of sophistication of these simulations suggests that the next step w i l l be to include the effects of global energy exchange by radiat ion. Indeed, there is a need for time-dependent 3D R T models as numerical tools for understanding the effect of inhomogeneities on the dynamica l evolution of the interstellar/intergalactic medium. For instance, the abi l i ty of gas to cool down and form structures depends crucial ly on the ionizat ional state of a whole array of different chemical elements, which i n tu rn directly depends on the local energy density and spectrum of the radiat ion field. Mos t of the previous studies which included some elements of radiative transfer into cosmo-logical simulations had to rely on a number of very simplifying assumptions, such as a unique dependence of the radiat ion field at each point on the local matter density and the loca l density gradient (the so called ' local opt ical depth approximation ' , Gned in & Ostriker 1997), or adopt an opt ical ly t h in Universe w i t h cooling rates depending on the strength and spectrum of a uniform U V background field w i t h complete self-shielding of gas from the external radia t ion field beyond some cr i t ica l density (Efstathiou 1992, C h i b a & N a t h 1994, Navarro & Steinmetz Chapter 1. Introduction to Cosmological Radiative Transfer 9 1997, Weinberg et al. 1997). It was only dur ing the last couple of years that it has become possible to numerically model some aspects of the full 3D R T problem. Challenges seem to abound, not least of a l l , the fact that the intensity of radia t ion i n general is a function of seven independent variables (three spatial coordinates, two angles, frequency and t ime). W h i l e for many applications it has been possible to reduce the dimensionali ty (e.g., i n classical stellar atmosphere models), the c lumpy state of the interstellar or intergalactic med ium does not provide any spatial symmetries. Moreover, coupled equations of radia t ion hydrodynamics ( R H D ) have very complicated structure, and are often of mixed advection-diffusion type which makes it very difficult to solve them numerically. Besides that, the radia t ion field i n opt ical ly th in regions usually evolves at the speed of light, y ie lding an enormous gap of many orders of magnitude between the characteristic time-scales for a system. One way to avoid the latter problem is to solve a l l equations on the fluid-flow time-scale. W h i l e there are arguments which seem to preserve causality i n such an approach (Mihalas & Miha las 1984, p. 342), even then the numerical solution is an incredibly difficult challenge (Stone et al. 1992). Since any two points can affect each other v i a the radiat ion field, even for a monochromatic problem, we must describe the propagation of the radiat ion field anisotropies i n the full five-dimensional space. Standard steady-state R T solvers, which have been widely used i n stellar atmosphere models, are not efficient i n this case. Non- local thermodynamic equi l ib r ium ( N L T E ) steady-state radiative transport relies on obtaining the numerical solution v i a an iterative pro-cess for the whole computat ional region at once, and is usually effective only for very simplif ied geometries. A n y refinement of the discretization gr id and /or increase i n the number of a tomic rate equations to compute N L T E effects w i l l necessarily result i n an exponential increase i n the number of iterations required to achieve the same accuracy. O n the other hand, the 3D solution of the steady-state transfer equation i n the absence of any spatial symmetries can often be ob-tained w i t h Monte Car lo methods (Park & Hong 1998). However, these methods demonstrate very slow convergence at higher resolutions and are hardly applicable i f one is interested i n Chapter 1. Introduction to Cosmological Radiative Transfer 10 following a time-dependent system. Recently, the problem of s imulat ing 3D inhomogeneous reionization wi th realistic radiative transfer has attracted considerable interest i n the scientific community. Umemura et al. (1998) calculated reionizat ion from z = 9 to z = 4, solving the 3D steady-state R T equation along wi th the time-dependent ionizat ion rate equations for hydrogen and hel ium. T h e radia t ion field was integrated along spatial dimensions using the method of short characteristics (discussed i n Sec. 2.5.2). The steady-state solution implies the assumption that the radiat ion field adjusts instantaneously to any changes i n the ionizat ion profile. One draw-back of this approach, however, is in low-density voids where there are probably enough L y m a n cont inuum photons to ionize every hydrogen atom, so that the velocity of I-fronts is s imply equal to the speed of light. T h e n the rate equations s t i l l have to be solved on the radiat ion propagation timescale. Besides, impl ic i t techniques in the presence of inhomogeneities w i l l become exponential ly complicated, i f we want to solve time-dependent rate equations for mult iple chemical species. N o r m a n et al. (1998) and A b e l et al. (1998a) present a scheme for solving the cosmological radiative transfer problem by decomposing the total radiat ion field into two parts: h ighly anisotropic direct ioniz ing radiat ion from point sources such as quasars and stellar clusters; and the diffuse component from recombinations i n the photoionized gas. In their method the direct ionizing radiat ion is being attenuated along a smal l number of rays, each of which is set up to pass through one of the few point sources w i th in the s imulat ion volume. T h e diffuse part of the radia t ion field is found w i t h a separate technique which might benefit from a nearly isotropic form of this component, for instance, through the use of the diffusion approximat ion. B o t h solutions are obtained neglecting the t ime dependent term i n the radiative transfer equation, w i t h the default time-step dictated by the speed of the atomic processes. If reionizat ion by quasars alone is ruled out (Madau 1998, however see H a i m a n & Loeb 1998), then I-fronts are most l ikely to be driven by ionizing photons from low-luminosi ty stellar sources at h igh redshifts. In this case, the pressure gradient across the ionizat ion zone is more Chapter 1. Introduction to Cosmological Radiative Transfer 11 l ikely to become important before the front is slowed down by the finite recombinat ion time. In this thesis I ignore hydrodynamical effects, concentrating on an efficient method to track supersonic I-fronts. M y approach is to solve the time-dependent R T coupled w i t h an impl ic i t local solver for the rate equations. Th i s method gives the correct speed of front propagation, and it also quickly converges to a steady-state solution for equi l ibr ium systems. However, I should note that un t i l a detailed comparison is made between explici t advection (at the speed of light) and the impl ic i t reconstruction (through an el l ipt ic solver), it is difficult to judge which e approach works best i n s imulat ing inhomogeneous reionization in detai l . As ide from their relevance to numerical cosmology, high-resolution mul t i -d imensional con-servation techniques and numerical R T i n part icular represent a formidable computa t ional prob-lem, the solut ion of which w i l l be useful i n many areas of computat ional astrophysics. T h e m a i n question I want to address i n this thesis is: can we numerical ly solve the problem of global en-ergy transfer i n a general 3D highly inhomogeneous medium v i a the radia t ion field? T h e present work describes one of the first attempts to attack this problem. T h i s thesis is organized as follows. In Chapter 2, I review techniques for solut ion of the time-dependent R T E . In Chapter 3, a method for accurate advection of radiat ion intensity i n 3D for monochromatic problems is developed. In Chapter 4, this scheme is extended to include frequency dependency of the radiat ion field, and I also introduce the detailed p r imord ia l chemistry equations into the model . In Chapter 5, I compute a number of models for different scenarios of cosmic reionization, and quantify some of the observational predictions of the epoch of reionization. F ina l ly , I present concluding remarks in Chapter 6. Chapter 2 Numerical Radiative Transfer in an Inhomogeneous Medium 2.1 Radiative Transfer in Moving Media For a moving medium the radiative transfer equation (RTE) in the inertial frame is (Mihalas & Mihalas 1984) iai^v)_ + n _ = _ c ot where n is the direction of photon propagation, I the intensity of radiation (per unit frequency per unit time per unit square area per unit solid angle), e the emissivity and K = K t n + KSC the total (absorption and scattering) extinction coefficient. For simplicity, we have omitted the variable denoting the dependence on spatial coordinates. Since the background material through which the light propagates is moving, e and K are anisotropic in general. For numerical simulations it is convenient to rewrite eq. (2.1) in the form where all matter properties are isotropic. Historically, the comoving RTE was first derived by Castor (1972), using a full general relativity formalism in application to spherically symmetric flows. Mihalas (1980) was the first to derive the comoving RTE exactly to all orders in v/c within special relativity, again in spherical symmetry. Through conservation of the number of photons one can show (Castor 1972) that the specific intensity, emissivity and opacity can be written in the frame comoving with the fluid using the following invariants: 12 Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 13 I(n,i>) = [^~J I 0 ( n 0 , ^ o ) , e(n,^) = (j^J e 0 ( f 0 ) , where the emissivity en and the absorption KQ are isotropic in the material 's frame, and the Doppler shift and the angle of aberration are obtained wi th the Lorentz transformations n 0 = n • v 7 n • v (2.3) (7 + l ) c j After some algebraic work eq. (2.1) can be wri t ten to order 0(v/c) i n the comoving frame (Buchler 1979): ld_ c dto n 0 ; v \ I, 1 + ^— — ^ 0 JL - n x • n 0 + IM cj v§ J ( n 0 • V 0 ) v IM c VQ . + d CM ( n 0 • V 0 ) ( v • n 0 ) IM + KM-IjM "0 ' (2.4) Here n x is a unit vector along the projection of the fluid velocity onto the plane perpendicular to no (in the comoving frame). In wr i t ing this equation we have already neglected the acceleration term, assuming that the fluid frame is inert ial . E q . (2.4) is not considerably easier to solve numerical ly than eq. (2.1): a l l matter properties are isotropic — but now we have to solve the equation i n the comoving frame, which means that the discretization grid has to move w i t h the fluid (the Lagrangian grid). E q . (2.1) or eq. (2.4) coupled w i t h the equations of hydrodynamics Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 14 mass conservation (continuity) : Dp/Dt + pV • v = 0, momentum conservation : pDv/Dt = — V p — p V $ + ^ X F F , ^ ^ matter energy conservation : pD(e/p)/Dt = —pV • v — AnksS + ck^E, equation of state : 1 p = p[p,Te), (writ ten i n the Lagrangian frame) constitute a general problem of radiat ion hydrodynamics ( R H D ) . Here p, p, e and v are the f luid density, pressure, energy density and velocity, respec-tively; <J> the external (gravitational) potential; E the energy density of the radiat ion field; F the radiat ive flux; Te the electron temperature; S the radiat ion source function; c the speed of light; XF, ks, various opacities (see eq. (2.9)-(2.11) for exact definitions); and D/Dt denotes the Lagrangian derivative. The numerical solution of a coupled R H D problem is, perhaps, the most fundamental challenge of computat ional astrophysics. T h e topic of this thesis is radiative transfer i n static media. In what follows, we w i l l neglect Doppler shifts w i th in the volume (unless indicated otherwise), the effects of aberration, and radiat ion pressure. Each of these elements can be included expl ic i t ly into the method we are developing here, since we w i l l solve the problem on the smallest of a l l advection times — the light crossing timescale. T h e goal of this s tudy is develop a fast and accurate scheme for energy transfer v i a radiat ion in general 3D inhomogeneous media. 2.2 T h e E q u a t i o n o f R a d i a t i v e T r a n s f e r i n S t a t i o n a r y M a t e r i a l T h e R T E (without cosmological terms) i n the static medium is (Chandrasekhar 1950) + n • V J „ = ev - KUIU, (2.6) c dt where I „ is the intensity of radiat ion i n direction n and e^'and KV are the local emissivity and opacity. E q . (2.6) is a 7-dimensional par t ia l differential equation, since the intensity lv Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 15 is a function of three spat ia l coordinates r, two angles i n the direct ion of photon propagat ion n, frequency v and time. To get a numerical solution to eq. (2.6), we have to discretize the intensity on a gr id i n the six-dimensional (without time) phase space which immediately results i n large amounts of data and the corresponding need for the computat ional power to deal w i t h these data. A gr id of say 64 points along each of the six independent variables would yie ld 512 G B of data to store just one time snapshot of the intensity i n double precision. Needless to say, these memory requirements are well beyond the capabilities of modern computers. A number of different techniques have been suggested to reduce the dimensionali ty of this problem. Historical ly, the radiative transfer problem was of crucial importance for model l ing stellar atmospheres. In the plane parallel geometry or i n spherical symmetry a l l differential equations were usually wri t ten as functions of two variables: the radius (or the mass w i t h i n that radius) and the frequency. The angle dependence was taken care of using various approx-imations, such as, historically, the Schwarzschild-Schuster approximat ion (Schuster 1905), or the Edd ing ton approximation, commonly used even today. 2.3 The Eddington Approximation The approach leading to the Edding ton approximat ion consists of rewri t ing eq. (2.6) as a system of angle-averaged moment equations (Krook 1955, Uhno & Spiegel 1966). In a static med ium integration over a l l directions and frequencies inside some range v\ < v < 1/2, w i t h corresponding angular weights, yields an infinite number of moment equations, the first two of which are dE —— = — V • F + ATVKSS — CKEE (radiation energy density conservation), (2.7) dt 1 <9F 1 and -TT^T- = — V • P Y F F (flux conservation). (2.8) <r ot c Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 16 Here we have introduced the frequency-averaged opacities X * = ' « * F W A > ; (2-9) KS = \ P K„S{v)dv; (2.10) KE = ^ P nvE(v)dv. (2.11) T h e energy density E, the flux F and the radiat ion pressure tensor PQ /g are the first three moments (zero, first and second rank tensors, respectively) of the specific intensity defined as angle-averaged integrals of the specific intensity (Unno & Spiegel 1966) E = - f IdQ, (2.12) c JA-K F = f Indfi, • (2.13) JA-K PaB = - f Inan0dn. (2.14) C JA-K K r o o k (1955) studied the hierarchy of moment equations (2.7) - (2.8) and the val id i ty of different approximations for one-dimensional problems (stellar atmospheres). Unno & Spiegel (1966) were the first to investigate the accuracy of moment equations for a general 3D inhomogeneous medium. Clearly, one can form an infinite hierarchy of moment equations. Unfortunately, like w i t h any system of differential moment equations, there is always one more unknown compared to the number of equations, thus requir ing an addi t ional equation of state to close this system. Over the years a number of different closure schemes have been suggested. T h e most common approach is to stop at the second moment equation and add an approximate relat ion between E Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 17 and P a /3; often it is wri t ten i n the form of the so-called dimensionless tensor Eddington factor (Auer & Miha las 1970) fa/3 = (2-15) Original ly , Edd ing ton derived a relation between E and Pap assuming a part icular form for 1(6, <j)). Since a number of different dependencies 1(9, 4>) can give the same value of f a / g , i t is the final form of the Edd ing ton tensor that is usually referred to as the Eddington approxima-tion. T h e factors fap essentially determine the degree of anisotropy of the radia t ion f i e ldbu t , unfortunately, can be calculated exactly only for the two extreme cases: (1) for free streaming (optical depth T « 1, i.e. the mean free path of a photon \ p is much larger than the typica l length L of the system) for a beam propagating i n the direct ion n the flux F = cEn and the Edding ton factor is highly anisotropic fQ/3 = nanp; (2.16) (2) in the diffusion l imi t ( r 2> 1) the flux |F | = cE • 0(\p/l) « O a n d the Edd ing ton factor is isotropic fa/3 = ^ a / 3 , ( 2 ' 1 7 ) where 5ap is the Kronecker delta symbol . A s for its value i n between, the degree of anisotropy has to be computed numerically. Clearly, it has to be done accurately; e.g., app ly ing the diffu-sion l imi t value of fap i n the near free-streaming regime would result i n overestimating the flux of radiat ion, effectively giving the wrong speed of energy propagation. Very often i n numerical calculations the Edd ing ton factors are evaluated using the flux-limited diffusion approximation, a modif icat ion of the diffusion theory which gives the right flux of radia t ion i n the opt ical ly th in Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 18 regime as well (see, e.g., A l m e & W i l s o n 1974 or Levermore & Pomran ing 1981). W h i l e this method gives the correct values i n bo th l imits (optically th in and optical ly thick media), it ef-fectively results i n dropping the time-derivative i n the radiat ion momentum equation (eq. (2.8), see Miha la s h Miha las 1984), which can lead to serious errors for intermediate opt ical depths. T h e best existing solution to recover the Eddington tensor is to resort to formal closure, reconstructing the angle-dependent intensity v i a a separate technique (Mihalas & Miha las 1984, Stone et al. 1992). Th i s is usually done by integrating intensities over a l l sources w i t h i n the computa t ional volume, taking into account boundary conditions and the opacity inside the volume. 2.4 Time-Dependent vs. Steady-State Techniques Another complicat ion i n obtaining a meaningful numerical solution of eq. (2.6) arises from the fact that radia t ion at least i n opt ical ly th in regions propagates at the speed of light, requir ing one to take very smal l time steps to avoid propagation of unphysical oscillations. A t the same time a l l fluid flows evolve on much slower, hydrodynamical timescales tf. A n extreme example of such a si tat ion would be formation of winds at the base of stellar atmospheres. Consider a low-mass ma in sequence (MS) star i n a low-mass X - r a y binary which is being irradiated by a hard flux of X-rays from a companion accreting neutron star. To compute the reaction of the M S star to this i rradiat ion, one would have to bu i ld a self-consistent dynamica l model of the system evolving at the smallest of the characteristic timescales: the light crossing timescale £R, through one pressure scale-height is of order 10~ 4 — l O " 3 s; the dynamica l timescale for w i n d acceleration is typical ly a few seconds; the timescale of energy transport v i a horizontal surface motions is of order a few days; while the global convection (mixing) timescale approaches one month. Thus the range of timescales, each of which is of crucia l interest for the behaviour of the irradiated star, covers more then"ten orders of magnitude, making it hopeless Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 19 Methods to solve the R T E numerically \ quasi-static solvers (multigrid techniques might be the best; all equations are solved implicitly) solution of the first two moment equations (energy & flux conservation) + some closure scheme (Eddington approximation) explicit advection (combined with an implicit chemistry solver) solution of as many moment equations as possible, with some "symmetric" closure scheme, e.g-diffusion aproximation flux-limited diffusion formal closure time-dependent ray tracing (discrete ordinates) A Monte Carlo methods (very low resolution) solve the photon Boltzmann equation directly in 5D (or 6D with frequency) phase space using modem, high-order multidimensional conservation schemes long characteristics isotropic radiation pressure tensor short characteristics diffusion isotropic higher-order moments (n>30) good approximation to free streaming Figure 2.1: Methods to solve the RTE numerically. Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 20 to talk about any numerical model at first glance. A common approach widely used in the literature (Mihalas & Miha las 1984, Stone et al. 1992, A b e l et al. 1998a) is to assume that the radiat ion field adjusts instantenously to any changes i n the matter density dis t r ibut ion, and the R T E can be solved impl i c i t l y on much longer timescales, e.g., on the f luid flow timescale tf. T h e arguments support ing this approach were or iginal ly put forward by D . Miha las (Mihalas & Miha la s 1984, p. 342) and are based on an order-of-magnitude analysis of different terms i n the equations of R H D (see also Stone et al. 1992). For example, in opt ical ly th in regions, the ratio of the light crossing timescale £R to ti is 0(v/c), which is a small number in most astrophysical situations, typical ly i n the range 10~ 4 — 10~ 2 . Therefore, in most cases one can safely ignore a l l terms i n the R H D equa-tions which are 0(v/c) and smaller, unless one is dealing wi th relat ivist ic flows. For processes occurr ing w i t h timescales larger than £R, one can s imply neglect the explici t variat ion of the radia t ion field w i t h t ime, assuming that it is always frozen to the current density d is t r ibut ion . T h i s radiat ion field is often referred to as quasi-static. The terms describing the interaction of radiat ion w i t h matter (emissivity and opacity) depend on the current microphysical state which is i n tu rn governed by atomic and molecular rate equations. These can have very steep t ime dependence evolving on timescales sometimes even faster than t^. Accordingly, a l l microphys-ica l popula t ion rate equations must be solved impl ic i t ly as well . Hence, a l l computations can be followed on the timescale of typical changes i n material propeties inside the computa t ional volume. W h i l e this approach is very powerful, the major uncertainty here is the accuracy at which the matter-radiat ion interaction is treated. Also , since the great majori ty of equations is solved impl ic i t ly , the numerical implementat ion of quasi-static techniques usually involves the solution of a large system of coupled non-linear el l ipt ic equations. The computer requirements for such calculations are known to grow steeply w i t h higher numerical resolution and the in t roduct ion of more complicated physics. Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 21 Mathemat ical ly , it is possible that the usage of fast mul t igr id techniques might actually give algorithms which scale linearly w i t h the amount of data processed at each t ime step. However, at the t ime of wr i t ing of this thesis mul t ig r id techniques have been successfully applied only to those systems of conservation laws which demonstrate a 'fair' degree of el l ipt ici ty, e.g., to the steady state Navier-Stokes equations (Sidilkover 1989). On-the other hand, several astrophysical problems allow one to follow the system of interest on a radia t ion propagation time-scale, without imposing a prohibi t ively large number of t ime steps. In the context of cosmological R T , we would like to resolve the characteristic distance between the sources of reionization. C o l d dark matter ( C D M ) cosmologies predict the collapse of the first baryonic objects as early as z = 30 — 50, w i t h the typica l Jeans mass of order 10 5 M 0 (see H a i m a n & Loeb 1997, and references therein). The corresponding comoving scale of fragmenting clouds is ~ 7kpc . If stellar sources reside i n p r imord ia l globular clusters of this mass, then for fih = 0.05 the average separation between these objects is A a ; ~ 20kpc (comoving). For explicit schemes the Courant condi t ion imposes a time-step A T 10 5 At < tR = ~ ^ h - 1 yrs. (2.18) C(l + Z) 1 + 2 Thus evolution from z = 20 to z = 3 takes ~ 7.5 x 10 8 h~l yrs (assuming the density parameter S7o = 1) and defining the Hubble constant to be HQ = 100 h k m s~l M p c - 1 ) . Subs t i tu t ing z = 3 and z = 20 into eq. (2.18) yields the total number of time-steps i n the range 20,000 — 150,000 over the entire course of evolution. In other words, to resolve reionization by 10 5 M 0 stellar clusters, we need to compute ~ few tens of thousands of time-steps on average, and for 10 8 M 0 clusters the number of steps required w i l l be even ten times smaller. A 100 3 gr id w i t h the required resolution w i l l result i n computat ional boxes of several M p c on a side. A ful l cosmological radiative transfer s imulat ion w i t h boxes at least this big and for a l l the required timesteps has not been feasible i n the past. Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 22 The rest of this chapter is devoted to the methods for solution of the time-dependent equa-t ion (2.6) on the light-crossing timescale £R. These are the ma in requirements for such a method: • fast advection i n 5D or 6D; • high spat ial resolution of 3D fronts (little dispersion); • correct front speed in 3D; • reasonable angular resolution for a given spatial mesh. There now exist a large variety of different techniques for the explici t solution of the t ime-dependent R T E i n 3D on the light propagation timescale. The major ones are l isted in F i g . 2.1. Because of the inherent complexity of mult i-dimensional advection, there is no single method which is best suited for cosmological R T . In sections 2.6 - 2.8 we w i l l go i n more detai l through some of the techniques we have tr ied to implement numerically. It seems likely that each of the methods covered — (1) an explici t 3D moment solver, (2) a combinat ion of a large number (n < 30) of moment equations w i th a 'symmetric ' closure scheme, or (3) a high-resolution con-servation scheme for the mult i -dimensional photon B o l t z m a n n equation — could be developed to be equally efficient i n studies of inhomogeneous I-fronts. Each of these techniques relies on ex-pl ic i t advection of wavefronts across the computat ional volume at the speed of light, hence, the timesteps are always l imi ted by the Courant condit ion. We have opt imized a further method, (4) time-dependent ray tracing to the degree that it can be used for model l ing reionizat ion at spat ial resolution of 64 3 , and describe it separately i n Chapter 3. 2.5 Formal Solution of the Transfer Equation T h e first step i n comput ing a numerical solution to the R T E is to b u i l d a gr id on which it can be discretized. Since the intensity of radiat ion depends bo th on the spat ial coordinates and on Chapter 2, Numerical Radiative Transfer in an Inhomogeneous Medium 23 the direct ion of photon propagation, one has to discretize simultaneously i n 3D space and i n two angles. Ideally, we would like to put enough angular data points NQN^ to resolve a l l 3D gr id points along the boundary of the computat ional volume. Therefore, i f N is the number of spat ia l g r id points i n each direction, the number of angles would scale roughly as NgN^ ~ N2 ensuring that even smal l clumps wi th in the volume are taken into account. Clearly, i f the average absorption inside the volume is large, the information about the precise locat ion of the sources w i l l be lost on a diffusion timescale, making it unnecessary to draw a lot of rays. In this case even the diffusion approximation w i l l do the job. However, generally one has to worry about the precise geometry of the field inside or close to opt ical ly th in regions. If the opacity is known beforehand (which is not usually the case i n the study of steep I-fronts), a simple rule is that the resolved number of angles should be a monotonical ly increasing function of Xp/L, the ratio of the mean free path of a photon to the size of the volume. F i g . 2.2 demonstrates a typical si tuation arising wi th poor angular resolution. T h e wave-fronts are transported along a finite number of rays going through the source i n the centre. A t large distances from the source the separation between the rays in the adjacent directions is several times larger than the 3D mesh size. In general, the business of interpolat ion between rays to get non-zero intensities is a t r icky one, since it can easily lead to the wrong I 7front speed i n 3 D . Smooth ing of the source improves the s i tuat ion slightly but clearly loses i n resolution (F ig . 2.3). Once the discretization of angular variables is done, one has to devise an efficient way of t ransport ing the angle-dependent intensity. For both quasi-static and explici t advection schemes two methods have been t radi t ional ly employed (Stone et al. 1992, Nakamoto & U m e m u r a 1998). 2.5.1 Long Characteristics Since i n a flat spacetime the paths of photons are straight lines, one can difference the R T E along mul t ip le lines i n 3D. The idea of long characteristics is to draw, at each point on a 3D mesh, a set Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 24 Figure 2.2: The effect of poor angular resolution around a point source of radiation. This figure shows a particular isosurface of equal ionization for a front expanding from a point source of radiation into a homogeneous, pure hydrogen medium. The true solution is an expanding HII spherical bubble. of rays covering all directions (Fig. 2.4). Reconstruction of the energy and radiation moments at each point then involves a simple integration at that point over all directions. Although it is possible to match the discretization in angle and space in such a way that invidual rays will pass exactly through more than one point, in general one has to draw the same set of rays individually through each point. Hence, the memory requirements for this method scale as 0(N4NgNfj)), since one has of order NgN^ rays for each cell and N points along each ray. Clearly, this method is the easiest to implement numerically, but it is also by far the most inefficient in CPU time and memory usage. A modification of this method relies on evaluation of the formal solution of eq. (2.6) only along the main axes and principal diagonals of the 3D mesh (NEWS, see Norman et al. 1998 for details). However, its angular resolution is only TT/4, placing its accuracy next to the diffusion Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 25 Diargy Enargy Figure 2.3: A model s imilar to the one in F i g . 2.2 except that here we show the projected (integrated column) energy density w i t h no smoothing and w i t h smoothing of the source, re-spectively. T h e colour gradient represents the energy density i n relative units. approximation. 2.5.2 Short Characteristics The goal of any formal closure scheme is to minimize the number of data points preserving the numerical accuracy in a mult i -dimensional phase space. The idea of short characteristics (see, e.g., Stone et al. 1992) is to calculate the intensity along short ray segments bounded by the edges of a 3D cell (F ig . 2.5). T h e values of the intensity at the ends of each ray segment are interpolated along pr inc ipa l spatial axes from the closest 3D grid points. T h e n one has to store only the angle-dependent intensity on the 3D gr id , w i th the overall memory requirements of order 0(N^N$N^). However, since the very nature of this method is based on iterations i n the plane not coinciding w i t h the direction of photon propagation, consecutively at each 3D mesh layer along the beam, it leads to large angular dispersion (see the numerical domain of Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 26 i . t 1 BJ ui ui, k Ik Figure 2.4: In the method of long characteristics rays are drawn to cover the entire sphere independently at each point. The intensity at point A in a given direction depends on the intensity advected from point B (and earlier from C) in that direction, corrected for emissivity and absorption along the ray segment BA (based on the diagram from Nakamoto & Umemura 1998). dependence in Fig. 2.5, based on Nakamoto & Umemura 1998; this effect was also demonstrated in the searchlight beam test in Stone et al. 1992). 2.6 Time-dependent 3D Moment Solver Assume that we can estimate Eddington factors via a stand-alone fast closure scheme, be it an accurate formal solution or a variant of the flux-limited diffusion approximation. The Eddington factors set the direction and speed of photon propagation. However, they do not give an explicit update to radiation variables — we still have to solve moment equations — •= —V • F dt Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 27 A / B2 C3 / /, Bl CI C2 . > Figure 2.5: The numerical domain of dependence i n the method of short characteristics. The intensity at the point B l is interpolated from the intensities at B 2 and B 3 , which in tu rn are based on interpolations at previous t ime steps (based on the diagram from Nakamoto & Umemura 1998). dF — = -c2V-(f^E). (2.19) In an attempt to perform fast explici t advection of radiat ion variables i n 3D on the light-crossing timescale, we have developed an unsplit mul t idimensional solver of the truncated moment equa-tions (2.19). The details of this method are described i n A p p e n d i x B . We have carried out a number of tests concluding that this solver gives a very accurate description for an isotropic radiat ion field (e.g., i n the diffusion l imi t ) . O n the other hand, for free streaming it tends to quickly lead to negative radiative energy densities E, since it does not have a bu i l t - in posi t ivi ty requirement. In practice, negative values of E appear i n those regions where there are large spatial gradients i n the Edd ing ton factors fap. How can we construct a scheme the advection part of which always gives positive energy densities? It seems likely that such a requirement has to be t ied to the one-dimensional R iemann Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 28 solver (eq. (B.7)). However, the full development of this technique is beyond the scope of this thesis. 2.7 Truncating the System of Moment Equations at a Higher Order The main l imi ta t ion of the Edd ing ton approximat ion is that it requires a separate, often very computat ional ly intensive closure scheme. Another — perhaps more mathematical ly elegant — idea which has been developed pr imar i ly for relativist ic R T is to extend the system of moment equations to higher orders. Anderson & Spiegel (1972) showed that, i n a med ium w i t h relat ivist ic differential motions, the Eddington approximat ion does not necessarily work even for opt ical ly thick cases. They develop a next-order approximat ion w i t h an arbi t rary number of moments of the type — where 1 is a 3D unit vector, and is its projection on the /j,a-axis — and suggest the use of a ' symmetr ic ' closure scheme which assumes that the highest order non-zero moment is isotropic. To i l lustrate the idea, let us expand the angle-dependent intensity i n terms of moments defined i n eq. (2.20). The expansion polynomials are generalized spherical functions (Anderson & Spiegel 1972) which are orthogonal, (2.20) for n ^  m, (2.22) and traceless, w i t h the first ones being Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium T h e intensity can then be wri t ten as n=0 where the coefficients A^ltl2..^n are expressed through the moments: 29 Y = 1, (2.23) Ya = n,a, Yo.fi = nanp - ^ 6 a j 3 , Yap1 = n-anpn-f - ^(Saprij + 5ianp + 8pjna), W) = E A = M , (2.25) A a = - 3 M a , Aap = y (M a / 3 - ^M5ap^j , T h e proposed closure scheme simply states ^ i W . . . / i „ = 0 for n > " m a x . (2-26) Choos ing n m a x = 1 gives the classical diffusion approximation w i t h an isotropic radia t ion field. Tak ing n m a x = 2 gives the next-order approximation. Th i s formalism was further developed by Thorne (1981) and Struchtrup (1997), bo th i n the covariant form i n the context of relativist ic R T . In this theory, i nd iv idua l moment equations are Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 30 coupled by matrices of mean absorption and scattering coefficients, which reduce to the well-known Rosseland mean values for local thermodynamic equi l ib r ium ( L T E ) . It is interesting to note that Struchtrup (1997) derived a closure scheme from maximiza t ion of entropy i n the rest frame of the f luid. More recently, Struchtrup (1998) applied this formalism to the case of homogeneous matter at rest, showing that the number of moments n m a x needed to represent a par t icular geometry of the radiat ion field varies w i th the opt ical depth. In fact, fimax = 30 is a good approximat ion for free streaming. For a system of moment equations, the speed of propagation i n a part icular direction (the characteristic speed) is given by the eigenvalue of the transport mat r ix i n the corresponding direction (see.eq. (B.3) i n A p p e n d i x B ) . It is easy to check that the highest characteristic speed (given by the largest eigenvalue) approaches the speed of light asymptot ical ly as n m a x — ¥ oo. E .g . , i n the diffusion approximat ion ( n m a x = 1), the advection speed is always c / \ / 3 , whereas at n m a x = 30 it is w i th in 1% of the speed of l ight. Here, we suggest that i t might be possible to use this property of moment equations to develop an automatic closure scheme which w i l l adjust itself to the current opt ical depth. T h e moment of t runcat ion n m a x can be given a simple functional dependence on the local absorption coefficient, whereas the numerical solution of the moment equations up to n m a x w i l l set the preferential direction of photon propagation at the corresponding characteristic speed. In i t i a l attempts proved to be computat ional ly time-consuming. However, ful l development of this approach is beyond the scope of this thesis, but future work could incorporate this idea into a numerical method and compare its performance to other closure schemes. 2.8 Solution of the Multidimensional Photon Boltzmann Equation Instead of solving the truncated system of moment equations, an alternative possibi l i ty is to employ a mult i -dimensional , high-resolution conservation scheme to solve the photon B o l t z m a n n equation di rect ly i n the five- (without frequency) or six-dimensional (wi th frequency) phase Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 31 space. T h e advantage of this technique is that one uses a single gr id for a l l calculations, w i t h high-order corrections reducing numerical dissipation, and it can, i n principle, be implemented i n curved spacetimes and/or non-inertial frames. For simplici ty, to il lustrate the basic idea of the method, let us consider monochromatic transfer i n flat spacetime. The photon Bo l t zmann equation is the conservation law for / R , the photon dis t r ibut ion function, defined as the number of photons per unit volume per unit momentum p = nhu/c, and related to the intensity of radia t ion v i a J ( r , n , M ) = ~ / R ( r , n , M ) - (2-28) c z Above we also introduced the photon four-momentum Ma — ——(l,n). (2.29) c For flat spacetimes — i n which photon paths are straight lines (Ma = 0) — eq. (2.27) can be easily reduced to the standard R T E (eq. (2.6)). In 3D Cartesian coordinates the monochromatic version of eq. (2.6) can be wri t ten as a 5D scalar advection problem w i t h coefficients constant i n the first three variables: It + A(e,<b)Ix + B{0,(t>)Iy + C(e,<l>)Iz = c(e-Kl), (2.30) where the intensity is a function of six independent variables I — I(r, n , t), the speed of ad-vection along the pr inc ipa l axes is (A, B,C) = c- (cos # cos <?!>, cos # s i n s i n ( 9 ) , and indices are short-hand notat ion for t ime and coordinate derivatives. Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 32 A high-resolution, mult i-dimensional wave-propagation algori thm for general time-dependent hyperbol ic systems of conservation laws was developed by R . LeVeque and collaborators ( L e V -eque 1997, 1998, Langseth & LeVeque 1997). A l though this technique is not based on the solu-t ion of the mult i -dimensional R iemann problem (evolution of the discontinuity at the interface between adjacent cells) but rather combines the first-order solut ion of I D R i e m a n n problems (in the directions normal to the interfaces between adjacent cells on a rectangular spat ial grid) w i t h special second-order mult idimensional corrections, it achieves a b ig step forward i n bo th s tabi l i ty and accuracy of the calculation. 2.8.1 Normal Flux Update We start by wr i t ing the first-order Godunov's scheme (LeVeque 1997) for fluxes normal to the cell interfaces. For example, an update along the x-axis involves the solution of the scalar linear wave equation It + A{0,<l>)Ix = O. (2.31) The corresponding Godunov's (1959) flux update is 7 n+i = i?-!L [ A + A J ^ / a + A- A J i + l / 2 ] , (2.32) where i labels cells, n labels timesteps, and A+ and A" are the positive and negative values of the wave speed A and A i j _ i y 2 = h ~~ li-i- Here h is the mesh size and k is the t ime step. T h e term A+AIi_1/2 = A+It — A + 7 , ; _ i is the difference of two waves travell ing to the right from points Xi and X i - \ , respectively. Alternat ively, A + A / j _ 1 / 2 + A~AIi+i/2 = [A+Il + A~Ii+\) - (A+Ii~i + A~l{) can be interpreted as the difference of net waves travelling to the right and to the left from Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 33 Xi. In comparison, expressions for A± for systems of equations or non-linear equations are considerably more complicated (see A p p e n d i x A ) . The full normal first-order update in 3D is a combination of updates along a l l three pr inc ipa l axes: fijtlm = lijklm ~ ft [^im^i-l^Jklm + Am^i+l/2,jklm\ ~f [ 5 f a A / » r l / 2 , W m + 5 i m A / i j + l/2,Wm] ( 2 " 3 3 ) ~TL \ChrAIijk-l/2,lm + Clm&Iijk+1/2,lm\ , where the subscripts denote discretization of the intensity on the numerical gr id: I(rijki n l m i t n ) = lijklm-T h e stabil i ty condi t ion of this scheme is given by the sum of the I D Courant conditions, that is the timestep k i n eq. (2.33) cannot exceed the t ime it takes a wave travelling at the sum of the m a x i m u m advection speeds along ind iv idua l axes [max(A) + m a x ( S ) + max(C)] to cross one computat ional gr id zone. 2.8.2 First-order Transverse Fluxes in 3D In flat spacetimes w i t h no scattering the trajectories of photons are straight lines, therefore, i n the 5D phase space ( r y ^ n j m ) information can only propagate to cells adjacent i n the first three indices (i,j,k). The normal flux update (eq. (2.33)) takes care of waves moving normal to cell interfaces. If the speed of advection i n the y-direction (given by the operator B i n eq. (2.30)) is not identical ly zero, some part of the wave moving from the cell (i — k, I, m) to (i, j, k, /, m) w i l l be advected to the cell (i, j + l,k,l, m) dur ing the same t ime step (F ig . 2.6). To take this transverse flux into account, we have to include spatial cross-derivative terms like Ixy (LeVeque 1997). Later on these terms w i l l be important for extending the method to second-order accuracy. Also , the transverse terms w i l l improve the stabil i ty of the mult i -dimensional Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 34 numerical scheme to the full Courant number of 1, i n units of the m a x i m u m advection speed i n any given direct ion (Langseth & LeVeque 1997). Consider the advection terms i n eq. (2.30): It = -(AIx + BIy + CIz). (2.34) A simple Taylor series expansion of the solution at t ime t + At around the already known solution at t ime t yields: /(r, n, t + At)= I(r, n, t) + J t(r, n, t)At + ~Itt(r, n, t)At2 + jUM(r, n, t)At3 + 0(At4), (2.35) where the th i rd and fourth R H S terms can be obtained by differentiating eq. (2.34): Itt = A{Alxx + BIyX + CIZX)+B{AIXy+BIyy+CIZy) + C(AIxz + BIyz + CIzz), (2.36) — Ittt = A2(AIXXX + BIyXX + CIZXX) + AB(AIXyX + BIyyX + CIzyx) + AC'(AIXzx "I- BIyzx + CIZZX) + BA(AIXXy + BIyxy + CIZXy) + B2(AIXyy + BIyyy + CIZyy) + BC(AIXZy + BIyzy + CIZZy) + 'CA(AIXXZ + BIyxz + CIZXZ) + CB(AIxyz + BIyyz + CIzyz) + C2(AIXZZ + BIyzz + CIZZZ). The first-order normal flux update outl ined i n the previous section takes into account only the first two terms on the R H S of eq. (2.35). Par t of the normal wave transported i n the rr-direction w i l l almost always (since B ^ 0 i n general, see F i g . 2.6) be advected along the y-and z T axes dur ing the same time step into neighbouring cells. The 0(At2) t e rm i n eq. (2.35) corresponds to the transverse propagation yielding an update to eq. (2.33) which can be wr i t ten schematically i n terms of contributions to transverse fluxes: Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 35 ij+l/2,fc(m' ,6 ij —1/2,Mm' ij+l/2,klmi ,b ij—l/2,klm' (2.37) ij+l i+lj+1 £ + A"Aq i + | / 2 j BA A q i + ] / 2 j i.J- 1 i+lj-l Figure 2.6: Propagat ion of transverse fluxes in 2D (based on LeVeque 1997). In this notation, the superscript b on the R H S labels the transverse wave update to the normal flux along the y-axis. Since eq. (2.30) is a single P D E , a l l wave spectra are scalars, and wave strengths are given by the positive and negative parts of advection matrices A, B and C. Other 0(At2) spatial cross-derivative terms (A±B±, A±C±, B±C±, C±A±, and C ± B ± , contr ibuting to $ a , $ f t and $ c ) are formed s imi lar ly to eq. (2.37). The resulting transverse flux update is then 1 k ~7iTBlmAtm^Ii-l/'2,jklm Ik 2h' 1 fcj A T ~2h l m l m + 1 l 2 J K L M Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium jn+l ijklm jn+l _ * (<f.a 1ijklm L \^ i+l /2, , — <ba \ Jklm i-l/2,jklmj (2.38) klm ij — l/2,klrn ®ijk-l/2,lm) 2.8.3 First-order Corner Fluxes in 3D Accordingly, i n three spatial dimensions, part of the transverse flux w i l l transported into the corner cell , a l l dur ing the same time step (Fig . 2.7). The corresponding update accounts for cross-derivative 0(Ats) terms i n eq. (2.35): -1 (ABCIzyx + ACBIyzx + BACIzxy + BCAIxzy + CABIyxz + CBAIxyz) At3. (2.39) comer flux transverse flux Figure 2.7: Propagat ion of transverse and corner fluxes i n 3D (based on LeVeque 1997). Since eq. (2.30) is scalar, the first-order corner update consists s imply of first-order upwind terms along each of the coordinate axis. The corner update corresponding to the first term i n Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 37 eq. (2.39) is 1 2 6 1 fk^ 2 6 U /2,jklmi k\* foj Atm BhnCL AIijk-l/2,lm > \ 2 _ g Atm BhnCL. AIijk-l/2,lm ~ * $ i + l / 2 , . 6 U 1 / / c N 2 6 6 U 1 2 im^lrn^lm^ijk-l^lm ' *i-l/2,jklm> 1 / ^ V Br.C,t A L j L i / 2 , m —> AtmBlTrfilm^ijk+\./2,lm > * ? + l / 2 J A : / m > AlmBlm^lmAUjk+l/2,lm > ®i-\/2,jklnv 6 U (2.40) Other terms contribute s imilar updates to the transverse fluxes $ a , $ f t and $ c in eq. (2.38). 2.8.4 Second-order Updates Second-order corrections can be performed separately for each of the normal , transverse and corner fluxes. Note that the transverse advection of normal waves i n eq. (2.38) already accounts for 0(At2) cross-derivative terms i n the Taylor expansion (eq. (2.35)). T h e ful l second-order accurate update formally has to include the terms A2IXX, B2Iyy and C2IZZ which is carried out by propagating the second-order correction wave for the normal update (eq. (2.33)) using techniques described i n A p p e n d i x A . However, as pointed out by LeVeque (1997) and Langseth & LeVeque (1997), to keep good stabil i ty properties, it is also useful to perform a transverse propagation of the correction wave. The implementat ion details can be found i n Langseth & LeVeque (1997). Chapter 2. Numerical Radiative Transfer in an Inhomogeneous Medium 38 2.8.5 Summary There are several clear advantages of the numerical solution of the photon Bolzmann ' s equation over the formal closure along a large number of rays. F i r s t of a l l , the direct B o l z m a n n solver uses a single 5D gr id for a l l calculations, e l iminat ing the need for interpolat ion between the 3D mesh and the rays. It can also natural ly deal w i t h mult iple point sources w i t h i n the volume. Since the wavefronts streaming from these sources are tracked explici t ly, at the speed of light, there is no need to calculate the opt ical depth along the line connecting every single gr id cell and the point source, as we do for direct ioniz ing photons (eq. (3.1)). However, a finite angular resolution i n the direct Bo l t zman solver w i l l s t i l l result i n a numerical diffusion s imi lar to the one shown i n F i g . 2.2, imp ly ing that the number of angles resolved should roughly match the 3D spatial resolution. Note that normal, transverse and corner propagation of bo th the first-order and the second-order correction waves have to be performed at each t ime step. We have numerically implemented the first-order monochromatic Bo lzmann solver in three spatial dimensions. Since the Bolzmann 's eq. (2.30) is a scalar P D E w i t h coefficients constant i n spat ial variables, and since the intensity is discretized on a single gr id i n the 5 D phase space, without the need to draw rays, its numerical implementat ion does not present much of a challenge on the programming side. P re l iminary tests have demonstrated large dissipation across the wavefront, indicat ing that second-order updates are essential. O n the other hand, experience has shown that the explici t ray tracing method described i n the next chapter can get s imi lar results for the diffuse part of the radiat ion field — although w i t h some tradeoffs — at a fraction of the computat ional cost. In the rest of the thesis we shal l concentrate on the ray-tracing algori thm. Chapter 3 Time-Dependent Ray Tracing 3.1 Description of the Method T h e basic idea of our technique is to solve eq. (2.6) directly for the angle-dependent intensity I(r, n, v, t) at each point. The total radiat ion field at each point is d iv ided into the direct (from ioniz ing sources) and diffuse (due to recombinations i n the gas) components: where k are the three indices i n a rectangular grid. T h e energy density ESJ^- k due to direct photons coming from point sources of ionizat ion can be easily calculated on the 3D gr id v i a summat ion over a l l sources (assuming that there are not too many of them w i t h i n the volume): huNph _T... I 1 ^ n,j, k,s < cta s * " t - ' i , j , f c , s I 0 if r - i j ) f c , s > cts, fri,j,k,s where ruuk,s = / OITV Jo is the opt ica l depth, r - j j ^ g is the physical distance between the current point and the source and ts is the age of the source. Note that the rate of emission of photons A ^ p n can be modified to allow for var iabi l i ty of sources on short timescales. Clearly, this calculat ion is very computer intensive since for each cell where we want to reconstruct the intensity we have to integrate the opacity along the line connecting this cell and the source, and the summat ion over a l l sources w i l l result i n C ( A r 4 A ? s r c ) operations. Fortunately, this update has to be done either on the 39 Chapter 3. Time-Dependent Ray Tracing 40 timescale of typica l changes i n the opacity i n the volume, or the var iabi l i ty timescale of sources. T h i s is achievable only for a smal l number of sources. For the diffuse component we use an upwind monotonic scheme-to propagate I D wavefronts J(r, n, v, t) along a large number of rays i n 3D at the speed of light. Fol lowing Stone & Miha la s (1992), we apply an operator split expl ic i t - impl ic i t scheme, i n which advection of radia t ion variables is treated expl ic i t ly and the atomic and molecular rate equations are solved impl i c i t l y and separately at each point. Unl ike A b e l et al. (1998a), we use rays to track the diffuse component of the radiat ion field and the direct ionizing background radiat ion (streaming into the computat ional volume). Since we need to draw rays essentially through every gr id point i n the 3 D volume, at first glance this approach appears to have very large memory requirements. However, efficient placing of the rays can significantly reduce the computat ional effort. 3.1.1 A Uniform and Isotropic Grid of Rays A t each time-step we are interested i n getting a solution for the mean radia t ion energy density and mater ial properties on a 3D N3 rectangular gr id . Instead of shooting rays though each gr id node i n 3D, we choose to cover the whole computat ional volume wi th a separate gr id of rays which is uniform bo th in space and in angular directions. Assuming that the computat ional volume corresponds to the range 0 < x, y, z < 1, we first construct a 2D rectangular base gr id containing N2 nodes (where N can be the same as the number of data points on each side of the computat ional volume) w i t h coordinates w i t h the origin at the centre of the cube / - / * — 1/2 1\ 1 y i j = V 3 { - H - - 2 ) + 2\ zv = \> with i,j = 1,...,N, (3.3) where the \ /3 factor ensures that the entire volume is covered w i t h rays, and we shoot rays normal to the base gr id through a l l of its gr id points. Note that the separation between 2D Chapter 3. Time-Dependent Ray Tracing 41 nodes is allowed to be larger than 1/N. To cover the whole volume with rays, we then rotate the base grid by an angle #/ — ir/2 around the y-axis and by 0;m around the z-axis, where the rotation angles are discretized to mimic an isotropic distribution of rays (with fewer azimuthal angles close to the poles) Figure 3.1: Rays shooting through the base grid in the direction (6,(j)). 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 Chapter 3. Time-Dependent Ray Tracing 42 ^ a n g l e s - -Ne(NB - 1). (3.5) 7T T h e to ta l number of rays scales as O(N2Aangies), which for large N is significantly smaller than N3 x A a n g ] e s (long characteristics, Sec. 2.5.1). In fact, w i t h N = 64 and i V a n g i e s = 110 (Ng = 10), we only require 218,242 rays. Since there are ~ 2/3 x JV points along each ray on average, the to ta l memory requirements of the method are of order 0(N3NgN(p), the same scaling as for the method of short characteristics (Sec. 2.5.2). However, time-dependent ray t racing does not have the angular dispersion which is typica l of short characteristics. Now, that we have the grid of rays and the 3D rectangular mesh, we have to specify the rules of interpolat ion between them. Th i s task is complicated, since we have to ensure that 3D wavefronts w i l l be reasonably sharp (the ideal s i tuat ion is to preserve 2nd- or 3rd-order accuracy from I D rays after interpolation to three spatial dimensions), as well the correct speed of I-fronts i n 3 D . A m o n g a large variety of possible schemes we tr ied, here we present the two most ly used in the tests. Since our pr imary goal is to minimize the number of rays (and thus greatly reduce, the computer memory requirements), using more than one ray to reconstruct the intensity at the point (i,j, k) i n the direct ion (#/, </>;m) implies that I D data points on the rays might fall outside the cell (i,j,k). Therefore, we need to use some smoothing which, however, w i l l preserve the correct front speed. A n alternative is to use a one-to-one correspondence between the I D data points along the rays and 3D mesh points. Th i s later method guarantees the right shock velocity (given that it is correct i n ID) but requires several times more rays. A : Smooth Interpolation Before we start our simulations, for each 3D grid point (i,j,k) we also store an array of the closest four rays going through its neighbourhood in the direction (0i,4>im)- Since the rays do not pass exactly through 3D gr id points, we use the values of the intensity on the four closest Chapter 3. Time-Dependent Ray Tracing 43 Figure 3.2: We calculate the radiation energy density assuming piece-wise linear angular variation of the intensity at each point in 3D. See text for details. rays in that direction to compute the angular-dependent intensity. Assume that for the point (i, j, k) the distances to the four closest rays going in the direction (9i, 4>im) axe d\, d2, c?3 and d^, respectively. We project the point k) onto these rays and read the values of the intensities, which we write as 1%, I2, h and I4. We then calculate the intensity at (i,j,k) in the direction (9i,cj)lm) according to d^dzdi /J^ did2d3d4 hMm = 2^ — J «V«/ 2- — 1  WQ' (3-6) q=l 1 ' q=l 1 Chapter 3. Time-Dependent Ray Tracing 44 where the weights wq for q > qo are set to zero i f qo < 4 rays were found i n the immediate neighbourhood of (i,j,k). Th i s might be the case close to the edges of the computat ional volume; we shall comment more on this while discussing the boundary conditions. T h e form of eq. (3.6) was chosen specifically because: (1) i f a ray labeled q happens to pass exactly through the point (i,j,k), then Iij}k,i,m = Iq', (2) i f a l l four rays encompass the point, m in ( J 9 ) < h,j,k,l,m < niax(Jg); and (3) i f (i,j,k) happens to be far from a l l four rays, then the result ing intensity w i l l just be an average of the four T^'s. B: One-to-One Correspondence W i t h the interpolation scheme i n eq. (3.6) there is always the possibi l i ty that the emissivi ty and opacity calculated i n a 3D cell w i l l be using I D intensities from rays i n adjacent cells. Th i s - unless the interpolat ion scheme is specially debugged for a fixed set of grids (spatial and angular) at a fixed resolution - w i l l cause the wrong speed of advection i n 3D, even though we might be using the correct high-resolution R iemann solver along I D rays. T h e cure is simple: the interpolat ion between 3D cells and I D data points on the rays has to be on a one-to-one correspondence. The energy density reconstructed from intensities from a given set of I D ray points has to contribute to the emissivity and opacity of that and only that same set of points. Therefore, i n a given direction (#;, </f>;m), the separation between rays should be smal l enough to guarantee at least one ray data point w i th in that 3D cell. The consequence is a few times (3 — 4 i n practice) larger number of rays required to cover the entire computat ional volume. 3.1.2 Angular Reconstruction A t each point on our 3D rectangular mesh we assume a piece-wise linear dependence of the intensity I on two angles, 9 and (j), I(M) = ( l-fr)( l-£*)/l-l,n-l + J Chapter 3. Time-Dependent Ray Tracing 45 (1 - &)€<l>Il-l,n + t,e£<pll,n + & ( 1 - i<t>)h,n-l, (3.7) w i t h i n a spherical rectangular element (or a spherical triangle adjacent to either of the poles) bounded by the angles 9U 9L-i i n 9 and 4>n and 0n_i i n (p, w i t h the rectangular gr id defined as - 1 r 9, = 7T 7 V 0 - 1 2 1 < I < Ne, 2im ~Ne~' 1 < n < Ng, (3.8) and e-et. - <Pn-l 9i - 9l-i <Pn - 9n-l We then integrate the intensity over 4ir w i t h appropriate weights to get the scalar radiat ion energy density at each point (l,n) of the rectangular spherical gr id C JAn C ; n n - l / 2 , / 2 (3.9) where the quadrature terms = L i<e<ei iij,k{0,<!>)dsi after integration inside each spherical rectangle or triangle (F ig . 3.2) can be wr i t ten as n - 1 / 2 A#, &4>n-l/2 X < 2-1/2 (hn + Iin-i){&9i-\/2 sin6>; + A c o s 6 > n _ 1 / 2 ) -- [h-i,n-i + ^ - i , n ) ( A 6 » ; _ 1 / 2 s i n ^ - x + A c o s 0 j _ 1 / 2 ) for rectangles, i i n ( A 0 j _ 1 / 2 s i n 0 j + A c o s 0 n _ 1 / 2 ) / 2 -- {h-\,n-i + ^ - i , n ) ( A 6 ' ; _ 1 / 2 s i n 0 ; _ i + A cos 9i_l/2)/2 for triangles. (3.10) The R H S of this expression is then mul t ip led by appropriate coefficients to allow for the non-orthogonal angular gr id (eq. 3.4), to yie ld the final weights. Chapter 3. Time-Dependent Ray Tracing 46 Since the advection part on the left-hand side of eq. (2.6) is s t r ic t ly linear, the simplest way to propagate intensities is just to shift wavefronts by one gr id zone at each time-step, accounting for sources and sinks of radiat ion. Assuming that a l l discretization points along each ray are s t r ic t ly equidistant, the intensity at a point j is updated s imply as - = / ; _ i e - ^ ' c A t + £jcAt. (3.11) Al ternat ively , one could take special care of the length of each ray segment contained w i t h i n a part icular 3D gr id cell , and use a scheme similar to the third-order-accurate piecewise parabol ic advection method ( P P A ) of Stone & Mihalas (1992). In either case we can track sharp discon-tinuities i n I D w i t h very l i t t le numerical diffusion, and, therefore, our approach is well suited to the calculat ion of I-fronts. Note that our placing of rays is functionally s imi lar to the method of of accelerated ray tracing employed by Nakamoto & Umemura (1998) for the quasi-static reconstruction (formal solution) of the angle-dependent intensity. However, here we are solving the advection problem explici t ly, propagating wavefronts i n the mult idimensional phase space at the speed of light, i n a sense making radiative transfer a local problem. 3.1.3 The Algorithm We start calculations by specifying the in i t i a l conditions (temperature, degree of ionizat ion, and i n the simplest cases no radiat ion field inside the volume) and boundary conditions (the intensity of radia t ion entering the volume - a l l outward flux at the edges can freely escape the computa t ional box). T h e inward flux at the boundaries is isotropic w i t h i n 2n and is s imply I + = £*• where Et, is the average background radiat ion energy density. Since each ray w i t h i n the volume starts and ends at the sides of the volume, we automatical ly have boundary conditions for each Chapter 3. Time-Dependent Ray Tracing 47 of the I D advection problems. The density field is kept static for a l l tests in this study, but since the radiat ion field is being evolved expl ic i t ly at the speed of light, one could easily evolve the under lying density dis t r ibut ion on a much slower fluid flow timescale i f desired. T h e course of the a lgori thm at each t ime step can be d iv ided into the following steps: • project the emissivity e„ and opacity nu from 3D grid to rays, • evolve wavefronts along ind iv idua l rays by At, • project intensities from rays to 3D to reconstruct Iw as a function of five variables • calculate the radiat ion energy density E(r), • solve a l l chemical equations v i a an impl ic i t iterative technique to get new state variables, separately at each point, and • compute new eu and KV. A t the beginning of each time-step we advect I D intensities according to eq. (3.11) along each ray. T h e numerical resolution along each ray is s imply set to the resolution 1/JV of the 3D rectangular mesh, so that along the ray i the point j has a coordinate Ai = )^A j = l,...,Nit Ni = liN, where Zj is the length of the ray segment inside the cube. Note that one can have much coarser I D grids along rays, speeding up the advection but sacrificing bo th spat ial and angular resolution. The advected intensities are then projected onto a 3D grid to reconstruct the mean energy density at each point using eq. (3.6 - 3.9). Th i s operation is one of the most demanding from the computat ional point of view, since at each of our N3 points we have to deal w i t h the angular dependence of the radiat ion field. W h e n this update is done, we solve Chapter 3. Time-Dependent Ray Tracing 48 the matter-radiat ion interaction equations impl i c i t ly to compute the local level populat ions (eq. (3.20)). T h i s gives us the 3D distr ibutions of emissivity and opacity (eq. (3.14 - 3.17)) which are then mapped back to the rays and used i n the advection scheme at the next time-step. T h i s simple scheme which we w i l l refer to as time-dependent ray tracing can be used as a stand-alone solver, or as a closure scheme for the system of moment equations through the use of variable Edd ing ton factors (as i n Stone et al. 1992). In the absence of any sinks and sources of radiat ion, the intensity is conserved exactly along each ray. Since the number of rays does not vary w i t h t ime, the advection part of our algori thm - i n other words, the proper solut ion of the I D R T equation (2.1) w i t h the zero right-hand side - can guarantee exact conservation of the to ta l number of photons i n the volume, once these photons have been mapped to the rays at t ime t = 0. O n the other hand, since emissivities and opacities entering the transfer equation as sources and sinks on the right-hand side are interpolated quantities, the total energy stored i n the radiat ion field at any given moment is subject to the interpolat ion error, which is a function of the spat ial and angular resolution. One advantage of the use of angle-averaged moments of radiat ion is that the advection mechanism is essentially reduced to 3 D , and it is relatively straightforward to implement the mult i -dimensional conservation scheme for the linear advection part of the moment equations. T h e n one could use a much denser spatial gr id for the solution of the moment equations, and a relatively course gr id for the angular reconstruction of the intensity of radia t ion v i a ray tracing. In practice, however, we have found that the mismatch between the spat ial resolutions of the moment solver and of the ray tracing usually leads to numerical instabili t ies. In what follows, we consider ray tracing only as a stand-alone solver. Another - perhaps, a better - way of coupling angular and spatial variations of the inten-sity may be an extension of the Spherical Harmonics Discrete Ordinate M e t h o d (Evans 1998). For steady-state transfer problems, instead of storing the radiat ion field, this method keeps track of the source function as a spherical harmonic series at each point. A l t h o u g h the direct Chapter 3. Time-Dependent Ray Tracing 49 implementat ion of this technique for time-dependent problems is probably not realistic, due to the lookback time (i.e. the finite speed of light propagation), the spherical harmonic represen-ta t ion of the radiat ion field might require less storage and might result i n smoother angular dependence as compared wi th a pure ray tracing approach. 3.1.4 Adaptive Time Stepping Note that i f I-fronts do not propagate fast compared to the speed of light, for instance, i n the cosmological context where the evolution on the light-crossing timescale w i l l normal ly require many thousand t ime steps, it is possible to update the radiat ion energy density E(r) once every few tens or few hundred t ime steps, while s t i l l evolving the intensities and solving a l l rate equations properly at the speed of light. We found that i n practice this shortcut leads to an increase i n 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 I-fronts to guarantee the correct rate of growth of ionized regions. To further accelerate this computat ion, we use adaptive t ime stepping to put lower t ime resolution on the 3D cells far from I-fronts. Since wavefronts cannot propagate further than one gr id zone dur ing one t ime step, only those cells which just experienced large change in E(r) and their immediate adjacent neighbours need a proper update of E(r). Depending on the w i d t h of I-fronts, typica l ly only a few percent of a l l 3D cells require the new, exact value of E(r). 3.2 Tests 3.2.1 Local Chemistry Equations Since i n our calculat ion a l l advection of radiat ion variables is performed expl ici t ly , we can solve N L T E rate equations separately at each point. Th i s makes it relatively easy to implement an impl ic i t solver for a l l atomic and molecular processes. Chapter 3. Time-Dependent Ray Tracing 50 To demonstrate the capabilities of explici t advection, instead of solving the proper chemistry equations for mult iple species w i t h p r imord ia l chemical composit ion, we have here adopted a simple toy model w i th just photoionization and radiative recombination i n a pure hydrogen medium. T h e impl ic i t solution of possibly stiff rate equations described below can be imple-mented i n a s imilar manner for more realistic chemistry models. One of these models is later developed i n Chapter 4. Let us assume that the change i n the degree of ionizat ion dxe/dt due to photoionizations is s imply propor t ional to the energy density of the monochromatic radiat ion field E of frequency v\, w i t h the proport ional i ty factor (1 — xe) * gm, where the photoionizat ion coefficient gui has the dimensionali ty c m 3 e r g - 1 s _ 1 . The total change i n fractional ionizat ion is ^ = gmE(l - xe) - a , ' e 2 n H a , (3.12) where a is the recombination coefficient. For an arbitrary opt ical depth, the finite difference approximat ion to eq. (3.12) can be wri t ten i n the conservation form, xe(t + At) - xe(t) E 1 - e~™At 2 A+ = I A + xe ma, 3.13 At hpvi nuAt w i t h opacity K = (1 - xe)nnhPvigm/c. (3-14) Instead of comput ing the true emissivity of the hydrogen medium (which depends on recombina-t ion rates to different atomic energy levels and a Maxwel l i an velocity d i s t r ibu t ion of recombining electrons), we obta in the total emissivity e of the gas by considering conservation of the thermal energy density Eih for matter, A £ t h = E ( l - e - C K A t ) - 47reAi = hPuinilAxe, (3.15) Chapter 3. Time-Dependent Ray Tracing 51 where a l l A symbols represent the change of variables dur ing one t ime step. T h e full recombi-nat ion coefficient a = ct\ + a s (3.16) is the sum of recombination coefficients to the ground state (ct\) and to a l l levels above the ground state (CUB, the 'case B ' recombination coefficient), v\ is the frequency just above the L y m a n l imi t , and we assume that recombinations i n L y m a n lines occur on a short timescale compared to ( i e r i H Q : B ) - 1 - Similar ly , the full emissivity (or gas energy loss through recombina-tions) is e = e i + e B , (3.17) where ei/e = a\/a. Th i s simple notat ion ensures radiat ion energy conservation i n eq. (2.6) for pure scattering of L y m a n continuum photons (i.e., when « B = 0)- E q . (3.13) does not account properly for the number of photons entering the volume, so that a large photoionizat ion coefficient (/HI 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. dr F, 1 — p - c r e A i ^ 1 < h 1 e , (3.18) dt hpvi rift At E q . (3.13,3.18) are solved separately at each point, given the local radia t ion energy density E. Discret izat ion of eq. (3.13) in t ime yields = (1 - e)h(xen) + 0fi(xen+1), (3.19) A t where / i (x e ) is just the right-hand side of eq. (3.13) and 1/2 < 6 < 1 for stability. Such equations can almost always be solved v i a Newton's method for smal l enough A t . L inear iz ing Chapter 3. Time-Dependent Ray Tracing 52 eq. (3.19), we get the (i + l ) - t h approximat ion to the value of xe at t ime t' d / i ( z e " + M dx^1'1 - l 1 + x 1-9 h(xen) + h(xe n+l , i (3.20) 9 9At which can then be iterated. 3.2.2 Test Parameters and Resolution For a l l of our test runs, except the study of the shadow behind a neutral c lump i n Sec. 3.2.6, we set up a numerical gr id w i t h dimensions 64 3 x 10 2 . We have experimented w i t h different angular resolutions finding that i V 2 n g l e s ~ 10 2 produces reasonably accurate results for the amount of gas c lumping typical ly found i n cosmological simulations. M u c h higher resolutions (Angles 3> 10 2) do not provide any significant gain, since, even at 10 2 , rays which are adjacent i n angles resolve a l l 3D gr id cells wi th in the mean free pa th of a photon. It is interesting to point out that - / V 2 n g l e s ~ 10 2 also matches the equivalent resolution of 64 5 data points for 5D advection. There are N2ng]es rays passing i n the immediate neighbourhood ( l / 6 4 3 t h of the to ta l computat ional volume) of each 3D gr id cell , each ray containing Np ~ 2/3 x 64 gr id nodes. Thus , in 5D we obtain the equivalent resolution of iVp x 64 3 x i V 2 n g l e s ~ 1.1 x 10 9 ~ 64 5 da ta Lower angular resolutions ( i V 2 n g l e s < 10 2 ) would probably suffice for those problems where the radia t ion field is nearly isotropic (that is, the mean free path of a photon is comparable to or smaller than the grid resolution), or where the angular variat ion is dominated by the lower frequency harmonics. For a l l test runs we take a 2.5 M p c (comoving) cosmological volume, w i t h the physical scale corresponding to z = 10. A l l densities (from low-density opt ical ly t h i n ambient gas to dense clumps) fall i n the range Clg&s = 0.001 — 0.05 (in units of the cr i t ica l density of the points. Chapter 3. Time-Dependent Ray Tracing 53 Universe), assuming a Hubble constant of 50 k m s - 1 M p c - 1 . Th i s low density contrast is chosen to cover the typica l range encountered i n cosmological hydrodynamics, and is ideally suited to demonstrate transient features during patchy ionizat ion. The absorption coefficients are taken to mimic complete self-shielding (assumed to be r J> 10) at neutral hydrogen co lumn densities higher than 1 0 1 8 c m - 2 (for which the opt ical depth at the L y m a n l imi t exceeds unity, see, e.g., Zhang et al. 1997), except where we probe different regimes, specifically in Sec. 3.2.5 where r = 10 corresponds to 1 0 1 7 c m - 2 , and Sec. 3.2.6 where r = 10 corresponds to 1 0 2 1 c m - 2 . 3.2.3 A n Isolated Spherical Expanding I-front A simple problem that would test the abi l i ty of the method to track I-fronts properly is that of a single, isolated Stromgren sphere expanding around a source of ionizing radiat ion (see, e.g., A b e l et al. 1998a). One difficulty of the current approach is that rays are drawn i n a way to cover the whole computat ional volume uniformly and isotropically. For a single point source of radia t ion this would mean that only a smal l number of rays pass through its neighbourhood. Hence, al though intensities along ind iv idua l rays are str ict ly conserved, there is no guarantee that the energy density has the right value far from a source of a specified luminosity. Instead, i n the first two tests (Sec. 3.2.3 and 3.2.4), the propagation of ionizing photons is computed on the 3D gr id , according to eq. (3.2). A t t ime t = 0 we turn on a point source, which starts to blow an expanding H n bubble around it . T h e speed of the I-front can be obtained analyt ical ly by equating the number of direct ioniz ing photons to the flux of neutral atoms crossing the front. In the absence of radiative recombinations (a = 0), the H l l bubble grows indefinitely (Spitzer 1968): where tc = ( A p h / 4 7 r c 3 n H ) • Parameters used for this calculat ion are the diffuse neutral gas density ftb = 0.01 (at z = 10) and the central source luminosi ty J V p h = 3 x 1 0 5 3 s - 1 (all emit ted Chapter 3. Time-Dependent Ray Tracing 54 i n photons above the hydrogen L y m a n l imi t ) . The comparison between the numerical speed of the I-front and the exact solution from eq. (3.2.3) is plotted i n F i g . 3.3. T h e difference between the two always stays w i t h i n one grid zone. t ime , s Figure 3.3: The analyt ical radius of the spherical H n bubble (solid line) compared to the numerical results (filled data points). The error bars indicate the resolution of the grid (±Ax), while the dotted lines give the radius at which the speed of the I-front drops below the speed of light. The horizontal dashed line at the top shows the radius at which the I-front reaches the boundary of the computat ional volume. 3.2.4 A n Isolated Stromgren Sphere in the Presence of a Density Gradient ' Consider a point source of radiat ion put into a density gradient along one of the p r inc ipa l axes of the cube. In the absence of diffuse radiat ion from H n regions (a% = 0) the only ioniz ing photons come directly from the source i n the centre, i n which case the shape of the ionized bubble would be a simple superposition of Stromgren spheres w i t h rad i i Rs{<fr) varying w i t h the az imuthal angle <f> and given by the classical solution (Spitzer 1968) Chapter 3. Time-Dependent Ray Tracing 55 S o = 4 ™ B I V n H ( r , < £ ) d r , (3.21) Jo where S'o is the photon product ion rate of the central source. For an exponential density gradient along the y-axis lognn = lognH.i + r cos i yo l Q g I n H ,2 (3.22) Ay ~ \ n H , i ( n H i , n H , 2 being the hydrogen densities on the opposite faces of the cube), the equ i l ib r ium Stromgren radius Rs (</>) is given by a simple equation So = [ea+bRsW(b2Rl(<t>) - 2bRs((f>) + 2) - 2ea (3.23) where a = 2 l n n H , i + -—-In Ay \ n H , i ^ = 2 c o s ^ ^ / n H ,2 and Ay \ n H , i y We take the physical densities on the opposite faces of the volume to be O g a S ; i = 0.001 and ^gas,2 = 0.05, and the luminosi ty to be A^ ph = 1 0 5 1 s _ 1 . To simplify the calculat ion by avoiding any realistic atomic rate equations, we take the value of the total hydrogen recombinat ion rate from H u m m e r (1994) for some fiducial temperature (T = 1 0 4 K ) . In F i g . 3.4 we plot a t ime sequence of models, for ionizat ion by a central source, w i t h no scattering of L y m a n photons (i.e. a± = 0 ) . The numerical solution at t = 4.1 Gyrs appears to be very close to the exact one for an equ i l ib r ium Stromgren sphere. T h e sharp t ransi t ion layer between the ionized and the neutral regions i n the high optical depth regime indicates that, indeed, the scheme introduces very l i t t le numerical diffusion even when extended to 3D (wi th in the accuracy of interpolat ion between the 3D mesh and the rays, see F i g . 3.5 for details). Chapter 3. Time-Dependent Ray Tracing 56 3.2.5 Ionization in the Presence of a U V Background The uniform coverage of the whole volume w i t h rays implies that extended sources of radia t ion w i l l be represented statist ically much better than point sources. A simple test mimick ing the evolution of dense clouds i n the presence of ionizing radiat ion is to enclose the computat ional region i n 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 t = 0. A n effective demonstrat ion of time-dependent ray tracing would be its abi l i ty to deal w i t h any d i s t r ibu t ion of state variables w i th in the s imulat ion volume. For this test, we set up a density condensation shaped like the acronym for ' radiat ion hydrodynamics ' ( R H D ) , w i t h a density 50 times that of the ambient homogeneous medium. The ambient medium has a constant density of f i g a s = 10~ 3 , and the energy density of the background radiat ion is E = 5 x 1 0 ~ 2 3 ^ i erg c m - 2 H z - 1 s r - 1 . F i g . 3.6 shows the result of this simulation. Most of the low-density environment is ionized on the radiat ion propagation timescale. It takes somewhat longer for ioniz ing photons to penetrate into the dense regions. Whether these regions can be ionized on a timescale of interest, depends on the rat io of the recombination timescale to the flux of background radiat ion. One can easily see ionizat ion 'eating i n ' to the neutral zone, e.g. i n the disappearance of the serifs on the letters at late times. Note that the w id th of the ionizat ion fronts does not usually exceed one gr id zone (F ig . 3.7). 3.2.6 Diffuse Radiation from H n Regions: I. Shadows behind Neutral Clouds Par t of the ioniz ing radiat ion at high redshifts comes in the form of hydrogen L y m a n cont inuum photons from recombinations i n diffuse ionized regions. The following test, s imula t ing the formation of shadow regions behind dense clouds at the resolution 32 3 x 10 2 , was adapted from Canto et al. (1998). A neutral c lump of radius Rc = 1.25 M p c (comoving) is being i l lumina ted by a paral lel flux F* = 6.34 x 1 0 ~ 1 6 hp1 erg c m - 2 s _ 1 H z - 1 of stellar ion iz ing photons (just Chapter 3. Time-Dependent Ray Tracing 57 above 13.6 e V ) from one side; here hp is P lanck ' s constant. A shadow behind the c lump is being photoionized by secondary recombination photons from the surrounding HII region (Fig . 3.8) of physical density np = 3.7 x 10~ 2 c m - 3 . Neglecting hydrodynamical effects, the w i d t h Rj of the shadow region can be estimated using a simple two-stream approximat ion (Canto et al. 1998) ? = TlTWri iT' w h e r e r ' = TT' <3-24) 1 + r/(21nr-j - 1) Rc and the dimensionless parameter £ is defined as (3-25) For £ 2 < 2, recombination L y m a n cont inuum photons from the i l lumina ted region w i l l eventually photoionize the shadow completely. For £ 2 > 2, radiative losses through low-energy cascade recombinat ion photons w i l l stop the I-front, forming a neutral cyl inder behind the dense c lump. St r ic t ly speaking, equations (3.24)-(3.25) are val id only for a shadow completely photoionized by secondary photons, and should be viewed as an approximat ion to I-fronts driven by scattering. For this r un we modify boundary conditions to include the effect of recombinat ion photons originat ing 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 addi t ional intensity of (ei) / ( « ) , where the angular brackets denote the average throughout the currently photoionized gas inside the box. Overa l l , we find a remarkably good agreement between the results of our models and the analyt ic solut ion (Canto et al. 1998), taking into account that . ionizat ion of the shadow is due to scattering i n the medium. In F i g . 3.9 we plot the radius R\ of the shadow neutral region as a function of £ i n our 3D numerical models. To vary the size of the neutral core, we fix the to ta l recombinat ion coefficient but we change the por t ion of recombinations into the ground Chapter 3. Time-Dependent Ray Tracing 58 state (which would produce more L y m a n continuum photons capable of ioniz ing the medium). The w i d t h of the I-front driven by secondary photons depends on the assumed opacity (the opt ica l depth of r = 10 corresponding to the column density of 1 0 2 1 c m - 2 ) . T h i s low opacity was chosen to reach equi l ib r ium quicker - equi l ibr ium itself does not depend on the opacity -but it cannot be too low, otherwise there w i l l be an unnecessary spread of the t ransi t ion zone i n the I-front over many grid cells. 3.2.7 Diffuse Radiation from H l l Regions: II. Ionization of a Central Void To demonstrate the abi l i ty of our scheme to handle scattering i n more complicated situations, we also set up a model w i t h ionizat ion of a central low-density void by secondary, recombinat ion photons. T h e void region is surrounded by two nested cubes w i t h opposite faces open (F ig . 3.10). T h e walls of the cubes are set to be much denser than the rest of the medium, to screen completely the central void from direct ioniz ing photons, and the ioniz ing U V flux is introduced at a l l faces of the computat ional volume. The total hydrogen recombination coefficient <x\ + an is again taken for the temperature T — 1 0 4 K . Similar to Sec. 4.4, we vary the amount of scattering i n the medium by changing the fraction q of atoms recombining into the ground state. In reality at T — 1 0 4 K the value q ss 0.38 (Hummer 1994) gives a solut ion i n between our extreme values of q. Simi lar to the test problem of Section 3.2.5, i f « B = 0, then the medium w i l l be ionized completely, since there is a constant flux of p r imord ia l ionizing photons. The speed of ionizat ion depends on the values of a\, OB , gm and nn- Note, however, that i f a.\ is too high, the I-front w i l l be very slow, since a large por t ion of the or iginal ioniz ing photons are scattered back. O n the other hand, i f a\ is too low, the I-front w i l l propagate much faster i n those regions where ionizat ion is driven by pr imord ia l photons, but i n shadowed regions there w i l l be too few recombinat ion photons. Thus, it seems that the speed of ionizat ion of the central void w i l l be the highest at some intermediate a\. Chapter 3. Time-Dependent Ray Tracing 59 In F i g . 3.11 we demonstrate ionization of the void region for models w i t h complete scat-tering and w i t h no scattering at a l l . The parameters used for this model are E = 5 x 1 0 _ 2 0 i v i e r g c m - 2 H z - 1 s r - 1 and the diffuse neutral gas density Q g a s = 0.01. A s expected, for the no-scattering model (q = 0) the central region remains neutral , since there is no direct path for the ionizing photons. However, for the model which includes scattering (q = 1), at least part of the central region becomes ionized. T h i s demonstrates that, our scheme is perfectly capable of dealing w i t h re-scattering of the ionizing photons. Beyond this simple test case, there are many astrophysical situations where progress can be made v i a numerical radiative transfer. For example, analytic solutions are often used, which are steady-state, and which assume a sharp boundary between the neutral and ionized zones. Us ing our numerical techniques it should be possible to follow general systems w i t h complex density inhomogeneities as well as regions of par t ia l ionizat ion. Chapter 3. Time-Dependent Ray Tracing 60 11.6 Myrs 290 Myrs 0 0.2 04 0 8 1 Figure 3.4: The cross-section of the numerical 3D ionizat ion front going through the source of ionizat ion for a model w i t h no scattering, shown for six output times. The density of the gas follows an exponential gradient (eq. 3.22) increasing by a factor of 50 from the upper to the lower side of the cube. The dashed line gives the analytic solution (eq. 3.23) for an equi l ib r ium Stromgren sphere without scattering. The wid th of the transi t ion layer between the neutral and fully ionized regions is consistent everywhere wi th the mean free path of ionizing photons. The ionizing source of luminosi ty TVph = 1 0 5 1 s _ 1 (all emitted i n photons just above the hydrogen L y m a n l imi t ) is marked by the asterisk. The nine contour levels correspond to ionizat ion fractions of xe = 0 . 1 , 0 . 2 , 0 . 9 . Chapter 3. Time-Dependent Ray Tracing 61 y , - y - \ r y y i y v -Figure 3.5: In order to cover the whole computat ional volume w i t h rays uniformly and isotropically, we set the number of rays inside a la t i tud ina l angle interval [6,6 + d6] propor t ional to cos 9. A n odd number of rays w i l l necessarily introduce slight left-right asymmetry for interpolat ion between the 3D rectangular gr id and the radial gr id for any non-central cross-section through the volume. In this figure we show a 2D schematic representation of the two grids employed i n our method. To simplify the plot, the source of radiat ion is assumed to be s i t t ing in the centre, and only rad ia l rays are d rawn here. Due to the finite angular resolution, interpolat ion between the 3D mesh and the rays is not symmetric w i th respect to the centre. Th i s shows up i n the asymmetry between the upper and the lower sides of the ionizat ion contours i n F i g . 3.11. Chapter 3. Time-Dependent Ray Tracing 62 ntime=400 ntime=500 I ntime=600 Figure 3.6: The isosurface of ionizat ion level xe = 0.5 is plotted at six t ime intervals for the model w i t h an ionizing background coming from outside of the cube. The density contrast between the ambient medium and the high-density acronym ' R H D ' (radiation hydrodynamics) is 50. Th i s s imulat ion demonstrates 5D advection on the radiat ion propagation time-scale. The numerical resolution is 64 5 . The effects of shielding are clearly visible dur ing par t ia l ionizat ion. Chapter 3. Time-Dependent Ray Tracing 63 Figure 3.7: Contour plot of ionizat ion i n a cross-section, through the data presented i n F i g . 3.6. There are 11 contour lines spaced linearly from = 0 to x» = 1. Note the sharpness of the ionizat ion boundaries, the shielding, and the fact that the most deeply embedded regions take the longest to ionize. Chapter 3. Time-Dependent Ray Tracing Figure 3.8: Cross-sections of the dense c lump (shaded) and the neutral shadow region behind it, for different fractions q = a\/a of recombination into ground states. The c lump is being i l luminated w i t h a parallel flux F* of direct stellar ionizing photons entering through one side of the box (marked wi th arrows), and the shadow region behind it is being photoionized by recombination L y m a n cont inuum photons from the surrounding HII gas. Each model has been evolved to an equi l ib r ium state, corresponding to many passages of the wavefront across the computat ional region. The nine contour levels correspond to ionizat ion fractions of xe = 0 . 1 , 0 . 2 , 0 . 9 . Th i s test is adapted from Canto et al. (1998). Chapter 3. Time-Dependent Ray Tracing 65 1 1 1 1_| 0 2 4 6 8 1 0 '£ Figure 3.9: The wid th R\ of the shadow behind a neutral c lump of radius Rc as a function of the dimensionless parameter £ of eq. (3.25). T h e solid l ine represents the two-stream analyt ic solution (Canto et al. 1998), assuming an infinitely sharp boundary. The points show the results of our 3D numerical model , plotted expl ic i t ly at xe = 0.5, w i t h the errorbars g iv ing the w i d t h of the I-front formed by secondary photons. Since the w i d t h of the neutral core i n our numerical model depends on mul t ip le scatterings i n 3D, we would not expect a detailed match; however, the agreement between the approximate analytic and the numerical solutions is remarkably good. Chapter 3. Time-Dependent Ray Tracing 66 Figure 3.10: The central void region in the test model in Sec. 3.2.7 is completely screened from direct ionizing photons by two nested cubes with opposite faces open. The only way to ionize the void is to allow for scattering of secondary photons in the side tunnels. Chapter 3. Time-Dependent Ray Tracing 67 Figure 3.11: Diagram demonstrating ionizat ion of the void region, which consists of two nested cubes w i t h opposite sides open. The central void and the side ' tunnels ' have the same density as the surrounding medium, and the density of the walls is set to be much higher to completely block a l l external ioniz ing photons. In this figure we plot one of the central cross-sections showing contours of ionizat ion at five different output times for the model w i th complete scattering (q — 1), and a quasi-equi l ibr ium configuration for the model w i th no scattering (q — 0) after n = 500 timesteps. The shaded arears represent neutral regions. Chapter 4 Including the Physics of Matter-Radiation Interaction 4.1 Frequency Dependence: Multigroup Transfer The numerical approach developed i n the previous section describes an efficient method for advection of monochromatic intensity i n the 5D phase space of three spat ial coordinates and two angles. E x a c t l y the same technique w i l l work for broad frequency intervals i f a l l terms w i t h subscript v i n eq. (2.6) are replaced w i t h frequency-integrated quantities. However, since cross-sections of photoionizat ion for most chemical species drop approximately as v~3 above the photoionizat ion threshold (e.g.j Tucker 1975, p. 236), higher frequencies have very l i t t le power i n frequency-integrated values of E, e and n. Therefore, a scheme which is devised to be sensitive to photoionizat ion of hydrogen w i l l give only a very crude description for ionizat ion rates of species w i t h higher ionizat ion potentials, as well as emissivities of those elements at higher frequencies. Therefore, we would like to include as many frequency intervals as possible. For the s tudy of cosmic reionizat ion it is sufficient to consider three frequency groups, above the hydrogen L y m a n l imi t (13.6 eV) and the single (24.6 eV) and double (54.4 eV) ionizations of hel ium, respectively, w i t h a l l frequency-dependent terms being averaged inside those groups. Besides hydrogen and single and double he l ium photoionization and recombination, there are several other chemical processes which could be easily included. However, at this stage we feel that in t roducing more than three frequency groups w i l l result i n an unnecessary complexity of the numerical model . We neglect entirely the effects of preheating of the neutral I G M by photons below 13.6 e V — which is a very complex process i n itself, acting through absorption i n a l l a tomic L y m a n lines w i t h 2 < n <^  150, and v i a absorption i n L y m a n and Werner lines of molecular H-2 (Ha iman 68 Chapter 4. Including the Physics of Matter-Radiation Interaction 69 et al. 1999) — and we consider radiative processes only above 13.6eV. The three frequency groups together w i th some other radiat ion thresholds are displayed schematically in F i g . 4.1. H Lyman limit H I He I He II oo Figure 4.1: Schematic representation of the three frequency groups, along w i t h some other radiat ion thresholds. We define group-averaged values as JT%+» /•"g+i Iv Ig= Iudv, eg = evdv, K i g = / Kv—dv, (4.1) Jisg Jug Jus i g where g = 1,2, 3 and 1/4 = 00. Integrating eq. (2.6) over frequency inside the interval [ i / g , i / g + i ] we get 1 dl - ~of + n • V / g = e g - K I g/ g. (4.2) Let us assume that we know the cross-section of some photochemical reaction (e.g., photoion-ization) o-*(v). Its contr ibut ion to the g-group average absorption is s imply « I g < =- / o*(v)Ivdv, (4.3) Chapter 4. Including the Physics of Matter-Radiation Interaction 70 and the total photoionizat ion rate is k>= I dQ / o*{v)~dv. (4.4) JA-K JV\ hv To calculate « i g and /e*, we are free to chose any spectral shape Iv{v) to satisfy the definit ion of Ig. For a simple power law W o g ( £ ) (4.5) we have h = h = h = hi a i - 1 1 OL2 — 1 ^ 0 3 a 3 - 1 V2 a 3 > 1. OL2 a i ^ 1, "2 7^  1, (4.6) (4.7) (4.8) In what follows, we w i l l be solving three scalar eq. (4.2) for g — 1, 2, 3 using the techniques described i n Chapter 3. 4.2 Rate Equations We use the atomic data from A b e l et al. (1997), inc luding 28 chemical rate equations for nine species: H , H + , H ~ , H2, H ^ , He, H e + , H e + + and e~. It is not entirely clear i f elements heavier than he l ium are present i n the intercluster environment during the epoch of reionization. W h i l e the absorption spectra of ul tra-high redshift (up to z ~ 5) quasars seem to indicate at least some presence of heavy elements in intergalactic clouds (e.g., see Cowie et al. 1995), there is a consensus that these elements were ejected into the intergalactic medium after the epoch of the first star formation. If most of cosmic reionization occurs v i a supersonic, R- type shock waves propagating outward from the point sources which made the first light i n the Universe, it seems Chapter 4. Including the Physics of Matter-Radiation Interaction 71 safe to assume that it takes some time for heavy elements to f i l l in the Universe more or less uniformly. Metals , i f present at that epoch, would supply extra free electrons, which would change the cooling function of the medium (e.g., facil i tat ing cooling i n the temperature range 10 4 — 10 7 K ) and affect some of the observable properties of patchy ionizat ion. However, at this stage we do not include metals into our models assuming uniform pr imord ia l composi t ion. Here we give a brief list of the chemical reactions which we take into account. For de-tails of actual microphysics and numerical values of reaction cross-sections .we refer the reader to A b e l et al. (1997) and Anninos et al. (1997), as well as to the on-line database of p r i -mord ia l gas chemistry made available by T . A b e l (the P r i m o r d i a l Gas Chemis t ry page, at http://zeus. ncsa. uiuc. edu:8080/^ abel/P Gas). C o l l i s i o n a l : (4.9) 1) H + e" -> H+ + 2e" 2 ) H + + e " - + H + 7 3) He + e" He+ + 2e~ 4) He+ + e" - » He + 7 5) He+ + e" -> H e + + + 2e~ 6) H e + + + e" -> H e + + 7 7) H + e~ ->• H ~ + 7 8) H ~ + H -4 H 2 + e~ 9) H+ + H -> H ^ + 7 10) H 2 + H —>• H 2 + H + 11) H2 + H + - ^ H 2 f + H 12) H 2 + e " - > 2 H + e" 13) H 2 + H -»• 3 H 14) H - + e" -> H + 2e" 15) H ~ + H ->• 2 H + e~ 16) H + + H ~ —> 2 H 17) H+ + H " ->• H ^ + e~ 18) H j + e" ->• 2 H 19) H ^ + H " -> H + H 2 R a d i a t i v e : (4.10) 20) H + 7 H + + e - 21) He + 7 H e + + e" 22) He+ + 7 -»• H e + + + e" 23) H " + 7 -> H + e~ 24) H 2 + 7 -> H+ + e" 25) H+ + 7 H + H+ 26) + 7 -»• H+ + e~ + H+ 27) H 2 + 7 -> H 2 * -> 2 H 28) H 2 + 7 -> 2 H We use the following definitions for chemical abundances: Chapter 4. Including the Physics of Matter-Radiation Interaction '<P = ™He total/™H total, • xi = n H o / n H total, x2 X4 = n H 2 / « H total, X5 x7 = « H e + / ' ( / ' n H total, ^8 where tp = 0.075 is hel ium abundance, electrons reads X\ + X2 + £3 + 2^4 + 2^5 = 1, x6 + x7 + xs = l, (4-12) X2 - x-j + X5 + tp(xr + 2xg) = Xg. We write down the chemical rate equations i n their most general form ^ = ^ E E <*ijkil>jil>kXjXk + E Mixi + i = 1 , . . . , 9, (4.13) where the coefficients a ^ , fiij and 7* are functions of reaction rates from A b e l et al. (1997) and Anninos et al. (1997), w i t h the first term on the R H S of eq. (4.13) accounting for col l is ional reactions, the second term representing radiative processes, the last term (which might give an impl ic i t conservation of species) being almost always zero, and ^ = (i,i,i,i,i,V,</>,V>, i f -Equat ions 4.13 describe evolution on drastically different timescales. Photo ioniza t ion can proceed on timescales much shorter than the light-crossing (of one gr id zone) t ime £R. O n the other hand, many species w i l l quickly acquire equ i l ib r ium abundances, and any further change i n their number density w i l l progress at the Hubble rate. The timescales entering eq. (4.13) span over ten orders of magnitude. For this reason, the numerical solut ion of chemical rate equations should be based on impl ic i t techniques. A l t h o u g h the routines for solution of = ™H+ /™H total, x3 = nH~ / n H total, (4-11) = n H + / n H total, XQ = nne/lpnu total, 2 = nHe++/V' nH total, ^9 = ne/nR t o t a i , and conservation of hydrogen and he l ium nuclei and Chapter 4. Including the Physics of Matter-Radiation Interaction 73 the systems of stiff ordinary differential equations are generally available (e.g., the popular Livermore Solver, L S O D A R , at http://gams.nist.gov; or the routine S T I F B S i n Numer ica l Recipes, Press et al. 1992), inside each solver it is necessary to update dynamica l ly a large variety of parameters, such as the chemical reaction rates which depend on the temperature, or the strength and spectrum of the radiat ion field. For this reason, a numerical solver of rate equations opt imized for automatic handl ing of terms wi th variable stiffness on the R H S of eq. (4.13) was developed specifically for this task. To simplify calculations, we assume that a^-fc, Pij and 7i do not change much wi th temperature during one integration t ime step 8t({Sxj}) (which can be much smaller than At; as discussed below), al though the impl ic i t solver is s t i l l temperature-dependent on the At timescale. We use l inearizat ion to construct an impl ic i t , first-order accurate chemistry solver. If we rewrite eq. (4.13) in the form dt and replace it w i t h a finite difference equation x?+l - x? * = ( i - ejMx?}) + oMixp1}) - 1 A t \ (4.15) where 0 > 1/2 for stability, then l inearizing eq. (4.15) we can obtain a system of equations dx?+1 | | ^ | | = | | ^ | | , (4.16) where the Jacobian on the L H S is composed of ^ = 2*nH [E^r1+^?+1] +e&-jtf (4"17) dxi \m J Here 8u is the Kronecker delta symbol and Chapter 4. Including the Physics of Matter-Radiation Interaction 74 Note that At is the time step for the advection problem. As far as the solution of the chemical reactions is concerned, At is the total period of time during which these reactions have to be evolved. The time step 5t({5xj}) for integration of the rate reactions can be much smaller. We obtain the implicit solution to eq. (4.13) using the iterations x?+1'k+l = x?+1>k + 6xi, with x?+1'° = x?. (4.18) Since the full change of variables, x™+1 — xf, during one advection time step At is of order O(At), we can always guarantee convergence simply by reducing At. Mathematically, the implicit solution (eq. (4.16)—(4.18)). is always stable; the stability is not degraded by the fact that the function <f>i{{xj}) is not very smooth, since it depends on the temperature and is obtained by interpolation from the tabular values. To guarantee high accuracy of the solution, even in the presence of very stiff photoionization terms, we use automatic integration stepping with the time step St({Sxj}) determined by the largest Jacobian element in eq. (4.17) for At -> oo. 4.2.1 Number Density of Absorption-Induced Processes Imagine a case where the photoionization rate in eq. (4.4) is very high. Since our chemistry solver is not completely implicit in the strictest sense — we assume that the kinetic temperature of the gas does not change considerably during an integration time step St (derived from the largest element in the Jacobian matrix in eq. (4.17)), hence, the absorption coefficient is constant within St — we might encounter a situation when the total number density of absorption-induced processes in a 3D cell exceeds the number of photons absorbed within that cell. To prevent this, we have to limit the total number of radiative reactions, similarly to eq. (3.18). But contrary to the toy chemistry model in Sec. (3.2.1), we now include multiple processes, both collisional and radiative, so that the number of photons absorbed due to some photochemical Chapter 4. Including the Physics of Matter-Radiation Interaction 7 5 reaction does not give a correct upper limit for the abundance of a particular species. However, for very steep photochemical reactions, we can subtract already thermalized (absorbed) photons (number density of absorbed photons) « ^ ( / ^ ) (1 ~ e CK*At)i (4-19) from the radiation energy density Eg at that point, within the chemistry solver. In the equation above, the quantity E \ _ Eg ^ - ^g+i [ j ^ j q e - l ( 4 2 0 ) is the local number density of photons in the frequency group g. Note that for large local absorption, the energy density Eg —> 0 on a timescale much shorter than At, and new photons have to be supplied to that cell in the advection part of the algorithm. 4.3 Absorption Coefficients To estimate the corresponding absorption coefficient, we can rewrite eq. (4.3) for a contribution from a given photochemical reaction, concluding that it is a constant within each frequency group: «ig <— m1-^ °*{v) ( - ) 8 du = nn1^ (a*) , (4.21) where both hg/Ig and (cr*) depend on the spectrum of the radiation field, but not on its intensity. The total absorption coefficient accounting for all photochemical reactions from eq. (4.10) in a given cell in 3D is then «ig 0g nu ( < 7 2 0 ) 5 + n H e (<72l) f l + « H e + (^22)9 + " H - (^23> f l + « H 2 ( ^ 2 4 > g + (4-22) + nn+ (a 2 5 ) g + + %+ (cr26)5 + m2 {°27)g + m2 (cr28)c Chapter 4. Including the Physics of Matter-Radiation Interaction 76 4.4 Photoionization Rates in the Low Optical Depth Regime A g a i n , consider some photochemical reaction wi th known cross-section <T*(I/). Tak ing into account the spectrum of the radiat ion field i n eq. (4.4), we can rewrite eq. (4.5) as where f dn £ i0g(e, ft f 8 + 1 ^ ( l ) = G*gcE0g, (4.23) E0G = -f i0g(e,<i>)<m, G*g = [Vg+l (-) 8'dv- (4-24) T h e coefficients G*g do not depend on the strength of the radia t ion field and can be com-puted once and for a l l at the start of the calculation. Once we know EG from the solut ion of eq. (4.2), we can find the radiat ion field strength EQS E0S = E g - f K ~ V ^ ~ " S + 1 ' " l * U 9 = 1 0 1 % • (4.25) [ ( "g - l)/"g> « g > 1,' 9 = 3-Note that this approximation w i l l work only for smal l opt ical depths, that is when the opt ical depth across a 3D cell is smaller than unity. 4.5 Photoionization Rates: the General Case If a computat ional cell on the 3D mesh is opaque to some photochemical reaction (which we denote w i t h "*") , one has to use the general expression for the frequency-dependent photoion-iza t ion rate TP 1 _ p-CKtu^-t where K*V is one of the ind iv idua l absorption coefficients from eq. (4.22). After integration over a specific frequency group we get Chapter 4. Including the Physics^ of Matter-Radiation Interaction 77 which can be wr i t ten approximately as E„ 1 cto , a e —1 g 4.6 E m i s s i v i t y I n s i d e E a c h F r e q u e n c y G r o u p T h e rates of radiative recombinations w i l l give us the total energy loss by gas through radiat ion. Not a l l of this energy w i l l be emitted i n L y m a n continuum photons (hv > 13.6 e V ) . Some might be scattered away close to the L y a line, and some por t ion of this energy w i l l be lost i n low-energy cascade recombination photons. In other words, knowing the rates for the radiative recombinations is not the same as knowing the emissivity inside each of our frequency groups. To calculate emissivities as defined in eq. (4.1), we have to use models for hydrogenic and helium-like ions to compute the energy loss as a function of frequency of the emit ted photons. 4.6.1 H y d r o g e n i c I o n s The energy loss coefficient due to recombination to the n- th level i n a hydrogenic ion is (Hummer 1994) UTe, Z, v) = ^ W 2 n - 2 e ( l + n \ f e ^ • (4.29) where v Z2hui e=zh7i ~ kBTe Chapter 4. Including the Physics of Matter-Radiation Interaction 78 a is the fine-structure constant, v\ is the hydrogen L y m a n l imi t , Z is the atomic number, k-Q is the Bo l t zmann constant and on(Z,e) is the cross-section of photoionizat ion for level n . Similar ly , the free-free energy loss coefficient MTe,Z,v) = - ^ ( f ) ' a2X2cZXl/2e-ug{f(u,X)^. (4.30) Here we have introduced the free-free Gaunt factor gfi(u,X), the C o m p t o n wavelength A c , and u = hf/k^Tf, (Hummer 1994). T h e n the frequency-dependent total energy loss due to recombinations and free-free emission for hydrogenic ions is kBTe Y<Pn(Te,Z,v)+MTe,Z,v) 47T where n0 is either n H + for hydrogen or nHe++ for doubly ionized hel ium. We take the numer-ical values for the hydrogenic photoionization cross-sections an(Z,e) from Storey & H u m m e r (1991) and the free-free Gaunt factor gs(u,X) from Hummer (1988), and integrate eq. (4.31) numerically to get e„dv for each frequency group. 4.6.2 Helium-Like Ions In he l ium atoms, at least at lower pr inc ipa l quantum numbers n\ and n,2, the nls configuration is dominated by both the electron-nucleus and electron-electron interactions. A s a result, the model for emissivity due to recombination of singly ionized he l ium H e + + e~ — • He(nZS) + 7 (4.32) is somewhat more complicated. We calculate the energy loss coefficient due to recombinations to the nlS state of neutral hel ium using (Hummer & Storey 1998) Chapter 4. Including the Physics of Matter-Radiation Interaction 79 , l s ( T e , * ) = ^ ( 2 Z + l ) f + A ) e " A e ^ - (4-33) Here, again, e = hu/Wyd = v/v] is the free electron energy i n R y d , a is the fine-structure constant, \i is the effective quantum number of state nlS, and _ R y d _ 1 5 7 ' 8 9 0 K T h e n the total recombination and free-free emissivity due to singly ionized he l ium is kBTe -n7ng oo n—1 E E E Vnis{%, V) + MTe, Z = 1,U) n = l ;=o S=i S = 3 (4.34) 477 which, again, is averaged wi th in each frequency group at the beginning of simulations. In our calculations we neglect the emissivity from di-electronic recombinations of hel ium. 4.6.3 Matter Energy Conservation Simi la r ly to the energy conservation in eq. (3.15), the change i n the matter energy density is determined by a l l time-dependent heating and cooling processes occurr ing inside a gas cell dur ing one timestep: Tn+i_Tn = [r(T) - A(T) ] A i • n H _ 1 [fM^i + 3:2 + 3:3) + ^kB{xA + x5) + + \kB(x& + x7 + x8)ip + \kBx$ where the heating (4.35) At g is given by photoionizat ion heating, and cooling I _ - C K g A t A ( T ) = 4vr E € g + A c i + A C o m + A c e + A m o i (4.37) Chapter 4. Including the Physics of Matter-Radiation Interaction 80 includes cooling due to radiative recombination and bremsstrahlung (combined into the first term, the sum of e g), coll isional ionizat ion cooling A c ; , Compton cooling A c 0 m (or heating — depending on whether the kinet ic temperature of the gas is higher or lower than the C M B temperature), coll isional excitat ion cooling A c e , and cooling A m 0 ] due to molecular hydrogen. The expressions for the latter four terms are taken from Anninos et al. (1997). A l so , since we are solving the R T E for a given, precomputed density field without a hydrodynarnical module, we also include adiabatic cooling due to expansion expl ic i t ly into our model . However, at the redshifts of interest to us (namely, at z > 5) Compton cooling is the dominant cooling process (Mira lda-Escude & Rees 1994), whereas adiabatic cooling is more important at lower redshifts. 4.7 Testing the Chemistry Model vs. Cosmic Recombination at z ~ 1000 A simple test for the chemistry solver described above is to follow the recombinat ion of the Universe at z ~ 1000. We assume that the matter d is t r ibut ion is uniform and the radiat ion field is bo th uniform and isotropic. For this test we take the "standard" model of G a l l i & P a l l a (1998; G P hereafter) w i th H0 = 6 7 k m s _ 1 , fi0 = 1 and fib = 0.0367. A l t h o u g h we could calculate the energy density of the cosmic radiat ion field i n a self-consistent way v i a for this par t icular test this would be a very bad approximation, since most of thermal radiat ion w i t h a temperature of hundreds or thousands of K e l v i n would be inside or below the first frequency group 5 = 1. Instead, we s imply use the blackbody Planck spectrum w i t h and the gas temperature is calculated from eq. (4.35). The physics of recombinat ion at z ~ 1000 E n + 1 = Ene-CKlsAt + 4 7 r e A * , (4.38) T(z) = 2.726(1 +z)K, (4.39) Chapter 4. Including the Physics of Matter-Radiation Interaction 81 is somewhat more complicated than radiative processes at lower redshifts (50 <; z <: 10), pr imar-i ly due to a much higher temperature and abundance of C M B photons dur ing recombinat ion (Seager et al. 1999, and references therein). Since our chemical model does not include any excited states of neutral hydrogen atoms — just the ground state and the cont inuum — we would not expect close quantitative agreement wi th other models (Seager et al. 1999, G P ) . Par t icular ly , we neglect L y a resonance photons, which keep most of the neutral hydrogen i n the first excited state after the epoch of recombination, either un t i l the expansion of the U n i -verse makes the I G M optical ly th in to Lya photons (they get redshifted out of the resonance at 10.2eV before they get absorbed), or un t i l there are sources of reheat ing/reionizat ion which would change the ionizat ional balance altogether. Compared to detailed studies of cosmic re-combinat ion, our photochemical model features a greatly reduced set of reactions for hydrogen and he l ium species (compare to G P ) , and also we do not include any heavy elements other than he l ium (e.g., deuter ium and l i th ium) . However, even for this crude model , we expect to get approximately the right redshift for the epoch of reionization, as well as roughly the correct abundance of p r imord ia l molecular hydrogen. We start evolving our photochemical model at z = 2000, assuming in i t i a l complete ionizat ion of bo th hydrogen and hel ium: X l = ( 0 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 1 + 2 ^ ) T . The result of this calculat ion is shown i n F i g . 4.2. Overa l l , we get good agreement wi th the adopted value of the redshift of recombination. Our residual ionizat ion fraction n e - /nn ~ 1.9 x 1 0 - 3 is overabundant by a factor of six compared to G P , or by a factor of three compared to some earlier studies (Stancil et al. 1996). Consequently, we get a s imi lar overabundance of molecular hydrogen, n n 2 / n H ~ 2.5 x 10~ 5 , which at high redshifts (10 3 > z i> 400) forms p r imar i ly through H ^ and at lower redshift (z < 110 - 200) through negative hydrogen ions H ~ (which are not destroyed by the C M B at these lower redshifts), and thus strongly depends Chapter 4. Including the Physics of Matter-Radiation Interaction 82 iooo Figure 4.2: Abundances of four chemical species vs. redshift. The vert ical axis is l o g x ; for % = 2 ( H + , solid line), i = 4 ( H 2 , dashed line), i = 7 ( H e + , dash-dot line) and i = 8 ( H e + + , dotted line), respectively. Chapter 4. Including the Physics of Matter-Radiation Interaction 83 on the amount of free electrons ( G P ) . Another shortcoming in this model is that we have not accounted for the popula t ion of excited states of w i th an excitat ion temperature equal to the temperature of the C M B just after recombination, art if icially suppressing the photodissociat ion of at 10 3 > z ^ 400 and, therefore, getting an overproduction of H2. In addi t ion, we do not include any reactions which form molecular combinations w i t h H e + (most important ly, H e H + , and the subsequent formation of H e ^ ) , hence, its relatively high residual abundance n H e +/nn «s 3.3 x 10~ 8 . Chapter 5 Modelling Inhomogeneous Reionization 5.1 The Density Distribution To be able to confront R T dur ing the epoch of reionization, we have to construct a realistic density field and evolve it through redshift. Since the R T E is solved on the light-crossing timescale, it is relatively straightforward to add a much slower evolving density field to the model . One way to do that is to use an N - b o d y code working paral lel to the R T solver. However, the subject of this thesis is the numerical solution of the R T E . For the sake of s implici ty, here we use a series of precomputed N-body density fields k ind ly supplied by M . W h i t e (see, e.g., M e i k s i n , W h i t e & Peacock 1999). The size of the volume is 2.5 M p c (in units comoving w i t h the Hubble flow), which at z = 15 corresponds to ~ 1/500 of the Hubble radius. A t the present epoch, 2.5 M p c is s imilar to the size of a typica l r ich cluster of galaxies. One interesting property of reionization is that it probes the formation of objects i n the mass range 10 5 — 10 9 M 0 (Ha iman et al. 1999), w i t h the boundaries set by the Jeans mass and the free-fall timescale, respectively. C M B - n o r m a l i z e d cosmological models, which also give approximately the right number of galaxies at low redshifts, can produce very different power spectra at z ^ 10. One way to discriminate between the competing models is to match the predicted patterns of reionization to observed quantities. However, at the moment we do not have enough computer power (both C P U time and memory) to run R T models for different cosmologies. Instead, we focus on a part icular cosmological model exploring different scenarios of quasar act ivi ty and star formation. 84 Chapter 5. Modelling Inhomogeneous Reionization 85 We interpolate the density field from a series of eight output times calculated for a cos-mological constant dominated C o l d Dark Mat te r cosmology (Oo = 0.3, OA = 0.7, 0 D = 0.04), and covering the redshift interval 30 > z > 8. Since for a "no evolution" model the comoving density 1 + Sijk(z) = riijk(z)/(l + z ) 3 = constant, to obta in a comoving density for an evolving field we use linear interpolation i n redshift of 1 + b~ijk{z) between adjacent precomputed fields. 5.2 Temperature of the I G M before Reheating Due to the expansion and adiabatic cooling after decoupling of matter and radiat ion at Zdec ~ few x 100, and before reheating at z J> 30 — 50, the kinetic temperature of the matter drops wel l below the temperature TCMB of the C M B : 2 726K TC M B = 2.726(1 + z)K and r u = — (1 + zf. (5.1) 1 "r 2dec I f there had been no coupl ing between gas and radiat ion after recombination at z r e c , then zdec = zrec, and prior to the epoch of reionization at z ~ 10 — 30 the matter temperature T^(z) would have dropped down to a fraction of a K e l v i n . However, even i n the absence of sources of preheating there is some residual ionization, w i t h free electrons exchanging energy wi th C M B photons (before decoupling at redshift of few hundred). In the expanding Universe the scattering of the C M B off free electrons (Compton heating) w i l l raise the electron temperature Te (Peebles 1993, p. 177) by dTe = xe 8 a T a R r C M B 4 _ ^ + 2Te(z) dt 1 + xe 3 m e c e 1 + z where xe is the fractional ionizat ion, and the second term on the R H S represents adiabatic cooling due to expansion. Coll is ions among electrons and neutral particles drive -> T e . The fractional ionizat ion can be calculated analytically, taking into account recombinat ion to excited states of hydrogen atoms (Peebles 1993, p. 167-173): Chapter 5. Modelling Inhomogeneous Reionization 86 d ne dt n fan-Is c-hva/kT n 1 + K(A + / 3 e ) n i s 1 + KAnls (5.3) where a e is the recombination coefficient to a l l excited states, (3e is the coefficient of ionizat ion and A = 8.23 s - 1 is the two-photon decay rate from the metastable 2s level to the ground state. Slight modifications to these formulae can be found i n Seager et al. (1999), but are unimportant for our purposes. Solut ion of eq. (5.2) - (5.3) shows that heating by C M B through residual free electrons keeps the kinetic temperature around ~ 6.5 K at z = 13 (as opposed to ~ 0.5 K w i t h no fractional ionization at al l) . However, the temperature of the clumps wi th in the I G M just pr ior to reionization is l ikely to be somewhat higher, since i n a non-homogeneous Universe denser structures w i l l not expand at the average rate of the Hubble flow, therefore, main ta in ing warmer temperatures. In addit ion, shocks from gradually collapsing objects w i l l contribute to the overall heating of the I G M (Sunyaev & Zeldovich 1972), not to mention radiative preheating by the very first stars and quasars. The calculation of the propagation of I-fronts which are going to heat the medium to many thousands of K e l v i n is not l ikely to be very sensitive to whether the in i t i a l gas temperature was a few K or a few tens of K . For the purpose of this study we adopt an in i t i a l ly uniform temperature T g a s = 30 K . We also assume that prior to the U V irradiat ion a l l the gas is neutral , and that the H 2 fraction is 3 x 10~ 6 ( A b e l et al. 1998b). 5.3 Preheating Ligh t coming from the first collapsed structures gradually builds up a cosmic radia t ion back-ground. Low-energy (below 13.6 eV) photons do not contribute to photoionizat ion traveling much further i n a neutral gas. The photons in the L y m a n and Werner bands of molecular hydrogen (in the range 11.2 — 13.6 eV) can photo-dissociate H 2 bo th i n the I G M and i n other collapsing objects (Haiman et al. 1999, C i a r d i et al. 1999). H a i m a n et al. (1997) showed that Chapter 5. Modelling Inhomogeneous Reionization 87 the flux which is required for complete H 2 photo-dissociation of the Universe is several orders of magnitude smaller than the flux needed for reionization. Since molecular hydrogen is the m a i n coolant for a gas of p r imord ia l composit ion, the U V background below 13.6 e V can strongly affect subsequent small-scale structure formation. O n l y halos w i t h masses high enough for efficient hydrogen L y a line cooling can continue further collapse. In addi t ion, at lower redshifts (z < 100), the density of the C M B is not sufficient to keep the major i ty of neutral hydrogen atoms above the ground state. T h e opt ical depth of the I G M inside atomic hydrogen L y m a n lines is 'extremely high (Haiman et al. 1999), and a photon redshifted to the frequency of any hydrogen L y m a n line w i t h 2 < n <^  150 gets immediately absorbed. L y a photons (n = 1) can be neglected, since they get re-emitted at the same frequency, while a l l h igh (n J> 150) L y m a n lines are not important because they have smal l opt ical depths. P r i o r to reionization, absorption i n L y m a n lines w i l l completely block a l l sources beyond some redshift Z m a x given by 1 "f" ^obs ^obs (Ha iman et al. 1999), where ui is the frequency of the L y m a n line closest from above to the observation frequency f 0 bs at the redshift z 0bs-In this thesis, we do not compute radiative transfer below 13.6 e V , neglecting a l l low-energy photons. In other words, one can argue that we do not conserve the radia t ion energy density since we allow low-frequency photons to escape the volume. The details of preheating, which are extremely important for subsequent structure formation, depend on the bui ld-up of the isotropic background below 13.6 e V as a function of redshift. Here, we only give a tool for calculat ing 3D R T i n an inhomogeneous medium, and compute propagation of ioniz ing photons. R T below 13.6 e V could be performed i n a similar way, subdiv id ing the interval 11.2 — 13.6 e V into i nd iv idua l frequency groups. Chapter 5. Modelling Inhomogeneous Reionization 88 5.4 Time Stepping with Variable Redshift We start our simulations at the redshift z = 15 when the first ioniz ing photons appear i n the volume. For explici t advection of wavefronts at the speed of light, the t ime steps are bound by the Courant condi t ion + 1 / i A t t ( t " ) l + * ( * n + 1 ) * ~~~c 1 + *(*») ' ( 5 - 5 j where fj, is the dimensionless Courant number, Ax(tn) is the 3D mesh size, and /J, < 1 for stabili ty. Hav ing adopted a part icular cosmological model , i.e. knowing the function z{t), we can easily solve eq. (5.5) numerically to determine tn+1. 5.5 Modelling the First Objects We have performed two runs of patchy reionization by sources residing inside the computat ional volume, as well as a run wi th reionization by an external background field. C o l d Dark Mat t e r cosmologies predict the collapse of objects w i t h masses ~ 10 5 — 10 6 M Q around the redshift z ~ 30 (Ha iman & Loeb 1998, H L 9 8 hereafter; Tegmark et al. 1997) w i th more massive structures forming at lower redshifts. It is only this relatively dense gas that is able to cool further and collapse into halos, leading to star formation or to quasar activity. To estimate the fraction of baryons that reside i n collapsed halos, the dark matter Jeans mass is often wr i t ten as ( A b e l et al. 1998b, H L 9 8 , O h 1999) M m i n = 10 8 [(1 + z)/10}-3/2 M Q , (5.6) which is s imply the cr i t ica l mass needed to reach the v i r i a l temperature of 1 0 4 K to be dense enough for efficient atomic hydrogen cooling. A l l halos w i th masses greater than M m i n are assumed to be vi r ia l ized. Note that this v i r ia l iza t ion threshold assumes that a l l molecular Chapter 5. Modelling Inhomogeneous Reionization 89 N U M B E R OF VIRIALIZED HALOS VS. REDSHIFT redshift number of mini-quasars 15 1 13 • 1 12 2 10 2 8 8 Table 5.1: Number of collapsed halos i n a 2.5 M p c (comoving) volume vs. redshift. The v i r ia l iza t ion mass Mm-m is given by eq. (5.6). hydrogen has been previously destroyed dur ing preheating — otherwise Mmm can be several orders of magnitude lower (HL98) . We assume that each halo wi th Mhaio > Mm\n hosts a mini-quasar, the luminosi ty of which is uniquely determined by the mass of the halo. In Table 5.1 we give the number of such mini-quasars inside our computat ional volume as a function of redshift. To estimate the pho-toionizat ion rate of the I G M by early quasars, it is essential to know their luminosi ty function ( L F ) . T h e number of observed quasars i n both opt ical and radio wavelengths declines steadily beyond z ~ 2.8 (Warren et al. 1994, Pe i 1995), al though some bright quasars were recently detected to as far as z ~ 5 (Fan et al. 1999, and references therein). A l so , it is possible that state-of-the-art X- r ay observations w i l l provide better constraints on the high-redshift quasar L F than currently inferred from optical data (Haiman & Loeb 1999). Several analyt ical models extrapolat ing the observed low-redshift L F of quasars to higher redshifts have been proposed i n the past (HL98, Pe i 1995, and references therein). F i t t i n g the observed quasar L F to semi-analytical models of structure formation introduces a number of uncertainties. The quasar L F is more or less established only at lower redshifts (z < 2.2 — 2.9), where most phenomenological models predict a double power law L F (for h igh and low luminosit ies, respectively), w i t h a redshift dependence of the form oc (1 + z)1 (Boyle, Shanks & Peterson 1988, Pe i 1995). A t higher redshifts, assuming that each black hole ( B H ) accretes matter at a rate close to the Edd ing ton l imi t (the mass infall is l imi ted by radia t ion pressure), Chapter 5. Modelling Inhomogeneous Reionization 90 REIONIZATION MODELS model description -Zstart ^finish ag ( R T module) A short-lived mini-quasars 15.0 14.4 1.5 1.5 1.5 B a single bright quasar 10.0 9.9 1.8 1.8 1.8 C an external U V field 13.0 12.6 1.5 1.5 1.5 Table 5.2: Reionizat ion models. Sma l l &L B landford (1992) have used a universal quasar L F to infer the bir thrate of quasar B H s . Haehnelt & Rees (1993), applying the Press-Schechter theory of halo formation and relat ing the B H - h a l o mass dependence to both the B H formation efficiency and the accretion luminosity, fitted the observed quasar L F . H L 9 8 used a s imi lar approach to derive the ultra-high-redshift (z > 10) quasar L F i n the context of reionization, finding that the early popula t ion of min i -quasars reionizes the Universe at z = 12. We take their quasar luminosi ty model to explore the effects of inhomogeneous reionization. Note that, i n general, we do not assume the same spectral energy d is t r ibu t ion for sources of radia t ion — these can be point-l ike quasars, star forming regions w i th in the volume or ioniz ing photons coming from the outside — and for the specific intensity i n eq. (4.5). Instead, we adopt a simple power-law spectrum for the calculation of R T of both direct and diffuse radiat ion, while sources are allowed to exhibit any spectral dependence. E .g . , for the U V background (Sec. 5.5.3) we assume a double power-law spectrum accounting separately for stellar and quasar photons. In Table 5.2 we list the reionization models we use for our computations. It is worth stressing again that direct ionizing photons from point sources w i th in the volume i n models A and B are computed on the 3D gr id using eq. (3.1)-(3.2), whereas photons reprocessed i n the I G M are followed w i t h ray tracing. Let us now examine these models i n more detail . 5.5.1 Mode l A : Reionization by Short-Lived Mini-Quasars For this model we adopt the high-redshift quasar luminosi ty function derived by H L 9 8 . The i r model is normalized to the observed quasar luminosi ty function at low redshifts (4.5 > z > 2.6) Chapter 5. Modelling Inhomogeneous Reionization 91 and assumes a constant ratio of the central black hole mass to the host halo mass M B H + . e = — — = constant. M h a l 0 For halos w i t h masses M h a ] 0 > Mmjn the luminosi ty of a quasar w i t h a B H mass M B H at a t ime t after its formation is L ( i ) = e M h ^ x l 4 x l 0 3 8 e - t / t o ^ g . M 0 s T h e fit to the observed luminosi ty function produces a model i n which quasars have very short lifetimes, of order to = 1 0 5 8 y r s , and the mass ratio is e = 1 0 ~ 3 ' 2 . To obtain the number of photons above 13.6eV, H L 9 8 use the median spectrum of E lv i s et al. (1994), to get Nph = Nph(> 13.6eV) = 6.6 x l O ^ s " 1 x ^ j ^ g - * / * f o r M h a J o > M m i n . Tak ing into account the variabil i ty of quasars and neglecting the lookback effects for the cal-culat ion of the opt ical depth, for a cell {i,j, k} inside our cosmological volume, the summat ion over a l l sources gives the energy density of direct ionizing photons (eq. (3.2)) =E (5.7, and i V p h ( i ) = 6 . 6 x l 0 4 7 s - 1 x ^ x < ; e 'o i f t > ts; 0 i f * < t s . (5.8) where ts is the t ime when the quasar turns on. Integration of eq. (5.7) w i t h i n each frequency group gives an approximate expression for « ^ E —2 / - rij^Mdu. (5.9) Chapter 5. Modelling Inhomogeneous Reionization 92 For model l ing quasar spectra, we adopt the same power law as for the intensity w i t h i n i nd iv idua l frequency groups (eq. (4.5)). Since the energy density ^ i V „ ) P n oc ^ Q g for v% < v < f g + i , we can write A ^ p h = A r0g,Ph — and NgtPh = / NVtPhdv. F r o m continuity of i V j , ) P r i over a l l frequencies the photon energy product ion w i t h i n each group is s imply where ( g and £ g are constants defining the adopted spectrum of quasar radiat ion: (5.10) , » 2 + l c. = s^Th-^fer]. e. = i I- - "2 ( S ) ° ' " & = (g) w + I h - •» (a)"] •« » = A te)"+1 h - (s)' + l ^ w ^ C*2 + l ^ I \ a i + l I V^ \ Ct2 +1 For this model we adopt the power law spectrum w i t h ag = 1.5 for a l l frequencies shortward of the hydrogen L y m a n l imi t (Haardt & M a d a u 1996, Haehnelt & Steinmetz 1998). T h e result of this run is presented i n F i g . 5.2 - 5.3. In F i g . 5.2 we show time evolution of the 3D hydrogen I-front propagating into the inhomogeneous I G M i n a cross-section through the quasar (the bright spot in the upper right corner), at four different redshifts. N ine contour levels of hydrogen ionizat ion are plotted on top of the color map representing the logar i thm of the to ta l hydrogen number density. For this s imulat ion we employ the full 64 3 (spatial) x 10 2 (angular) x 3 (frequency) resolution. The quasar turns on i n the densest cell in the volume at z = 15, and its luminosi ty is decaying on the timescale of to = 1 0 5 ' 8 yrs according to eq. (5.8), un t i l i t becomes negligible around z0g ~ 14.92. T h e in i t i a l expansion of the hydrogen I-front' is dr iven entirely by direct ionizing photons from the quasar. Close to the quasar, the diffuse radiat ion from the H11 region builds up slowly, dominat ing the radiat ion field near the 1-front Chapter 5. Modelling Inhomogeneous Reionization 93 after z0g. A r o u n d this t ime the I-front starts to show patchiness, due to percolat ion of the radia t ion field through the density inhomogeneities. In a separate run we have verified that this patchiness has nothing to do w i t h the finite angular resolution. T h e intensity of secondary photons w i t h i n the HII volume is not large enough to keep up w i t h the in i t i a l fast expansion of the front. A s it slows down, its thickness clearly reflects the mean free pa th of photons, being wider in the low-density regions and dropping below the resolution length (2.5/64 M p c comoving) inside the dense filament below the quasar (F ig . 5.2). Singly ionized he l ium shows a th in , shell-like structure bounded by neutral he l ium on the outer side and by doubly ionized hel ium closer to the quasar (F ig . 5.3). T h e H e l l front has the shape s imilar to the hydrogen I-front, bo th emphasizing the under lying density profile. The size of the H e l l void is smaller due to the higher ionizat ion potent ial of hel ium. We do not compute directly heating of the neutral I G M by photons wi th wavelengths greater than the hydrogen L y m a n edge, accounting only for hv > 13.6 e V radiat ion. Thus , the result ing temperature ahead of the hydrogen I-front might be somewhat different from more detailed y studies considering R T of L y a photons. In F i g . 5.4 we show the temperature i n the same cross-section at z = 14.785. Note that there is a prominent region where the neutral m e d i u m is heated to 10 3 —10 4 K ahead of the I-front, whereas hydrogen ionizat ion occurs at (3 —4) x 10 4 K . 5.5.2 Model B : Reionization by Bright Quasars T h e major uncertainty i n the previous model is the number, d is t r ibut ion and the luminosi ty function of quasars at z > 5. It is l ikely that preheating of the I G M by the earliest structures at z ~ 15 - either by light (soft X-rays or L y a ) or through hydrodynamica l shocks - w i l l raise the Jeans mass throughout the Universe, preventing further collapse i n many regions. In this case one might expect the medium to stay neutral un t i l bigger halos form closer to z ~ 10. In this section we compute a model for the reionization of a periodic 2.5 M p c (comoving) volume by a single bright (Nph = 1 0 5 4 s _ 1 ) quasar switching on at z — 10. A l so , for this calculat ion Chapter 5. Modelling Inhomogeneous Reionization 94 we adopt a steeper spectral index, ag = 1.8, throughout a l l frequency groups, mimick ing soft reprocessed radiat ion originating i n the circumnuclear region of the quasar. F igure 5.5 illustrates the resulting profile of the time-dependent hydrogen I-front at the beginning of the run. The dense filaments near the host halo are ionized almost on the light-crossing time, due to the high luminosi ty of the source. Not affected by the quasar finite .lifetime, ionizat ion proceeds quickly, and by z ~ 9.6 most of the low-density gas i n the volume is ionized. In F i g . 5.8 we show the energy density of the radiat ion field Eg i n three frequency groups (g = 1, 2, 3), i n a cross-section through M o d e l B at z = 9.473. 5.5.3 Model C: Reionization by an External U V Field Let us now consider a model w i t h direct ionizing photons coming exclusively from outside the volume. T h e common parameterization of the uniform U V background is a simple, redshift-dependent power-law 1 r ( \ J„(z) = — / Iv(n,z)dSl = J-2i{z) x 1 (T 2 1 . I — I erg c m " 3 &~ l H z " 1 s r _ 1 , (5.11) 47T H-K \U%) V„ < V < V, g+1-For this scenario, we have taken an ionizing background evolving w i t h redshift as proposed by A b e l & Haehnelt (1999), except that we shift it towards higher redshift. O u r form is *w - { r + i ^ e " , ( ' " 9 > / , , * ( ^ ) " 5 + 8 (5-12) + j e - K * - 9 ) / 2 - 5 ] 3 fiiV1"8} x 1 0 - 2 1 e r g c m ~ 3 s" 1 H z " 1 s r " 1 , 1 + [7/(* - 8)] 4 U / J w i t h the first term accounting for stellar reionization (wi th a soft component, a g = 5) around z ~ 12 — 14, and the second term representing the harder quasar radiat ion reaching its m a x i m u m Chapter 5. Modelling Inhomogeneous Reionization 95 closer to z ~ 11 — 12. For the processed radiat ion (eq. (4.5)) we assume the spectral power law a g = 1.5 (g = 1 , . . . ,3) . Figures 5.9 - 5.12 show the results of this run. To speed up the evolution, we start calculations when the background field has been already established, at z = 13. In addi t ion, the spat ial resolution is reduced to 32 3 . A l l low density hydrogen gas is ionized wi th in about ten passages of the wavefront across the volume, by z ~ 12.8. Denser filaments (p/p ~ 10 — 50) survive much longer, at least through z ~ 12.6 when we stopped the calculations. Note that this s imulat ion does not have the spatial resolution required to compute the evolution of neutral regions to much lower redshifts. In F i g . 5.13 we present a cross-section through the hydrogen I-front computed at three different spat ia l resolutions (8 3 x 8 2 , 16 3 x 8 2 and 32 3 x 8 2 , respectively), at the same output time. It can be seen that the posi t ion of the I-front converges as we go to higher resolution. 5.5.4 Criteria for Reionization by Population III Stars In the case of numerous and more uniformly distr ibuted faint stellar sources, w i t h softer spectra than for quasars, one would expect to witness a much more homogeneous, and, perhaps, slower reionization. In this section, we discuss how to set up conditions for numerical reionizat ion by stellar sources inside the computational volume. Because of the inherent complexi ty of the associated physics, there is s t i l l no standard theory of star formation. Inside v i r ia l ized dark matter halos, the gas is supported by the pressure gradient. Its further collapse and fragmentation into stars is possible only i f there is addi t ional cooling which w i l l loosen pressure support. For a gas of pr imordia l composit ion, the most efficient coolant is molecular hydrogen (H2), through coll is ional excitat ion of v ibra t ional and rotat ional levels and subsequent radiative transitions (Lepp & ShuII 1983, 1984; H a i m a n et al. 1996a; Tegmark et al. 1997). T h e post-recombination fraction of molecular hydrogen is ~ 3 x 10~ 6 . Inside dense halos, i n the absence of the U V background, this fraction can rise to ~ 10~ 4 (Haiman et al. 1996b, Tegmark et al. 1997) Chapter 5. Modelling Inhomogeneous Reionization 96 faci l i tat ing fast cooling. In fact, H a i m a n et al. (1996b) concluded that it is possible, contrary to naive expectations, that the increased abundance of free electrons due to a U V i r radia t ion of cold gas clouds can catalyze the formation of molecular hydrogen, accelerating line cooling and triggering the eventual collapse and fragmentaion of these clouds. However, it is more l ikely that i n lower density neutral regions, just prior to the epoch of reionization, U V photons below the ionizat ion threshold corning from the very first collapsed objects w i l l propagate ahead of I-fronts easily dissociating H 2 molecules (Haiman et al. 1997). D u r i n g this preheating the soft U V flux required to completely photo-dissociate the Universe is several orders of magnitude smaller than the flux needed for reionization (Haiman et al. 1997). Now, it seems f i rmly established that this preheating of the I G M — often referred to as a "negative feedback" — w i l l suppress molecular hydrogen cooling, effectively ha l t ing any star formation i n objects w i t h masses below 1 O 5 ~ 6 M 0 . The former w i l l remain neutral , pressure-supported gas clouds which later either get slowly evaporated by the U V background or merge into larger objects (Haiman et al. 1999). Therefore, one could assume that if reionizat ion is caused by the early popula t ion of stars, these must reside only i n v i r ia l ized halos w i t h masses of at least 1 0 7 - 8 M 0 (or v i r i a l temperatures in excess of 1 0 4 K , H a i m a n & Loeb 1997) inside which molecular line cooling is efficient before preheating. O n the other hand, recently, H a i m a n et al. (1999) performed detailed calculation of the bui ld-up of the U V background assuming that the first stellar ac t iv i ty occurs i n vir ia l ized halos appearing near z ~ 30. T h e corresponding Jeans mass is of order 1 0 4 - 5 M 0 g iving a v i r i a l temperature of a few hundred K e l v i n . Star formation inside smal l clouds suppresses the H 2 abundance i n these halos and shuts itself off before reionizat ion. H a i m a n et al. (1999) conclude that their calculations confirm the existence of the negative feedback, but find that it depends sensitively on the spectrum of the first sources. C i a r d i et al. (1998) found that the soft U V background in the L y m a n and Werner bands builds up to ~ 1 0 - 3 0 — 1 0 _ 2 7 e r g c m ~ 2 s _ 1 H z - 1 which is not sufficient to provide the negative feedback on structure formation i n the redshift interval 30 > z > 20. Chapter 5. Modelling Inhomogeneous Reionization 97 In view of these uncertainties, it is quite questionable to identify star forming regions w i t h a density cri ter ion alone, without relating them to a threshold molecular hydrogen fraction, above which gas clouds can efficiently cool, fragment and form stars. Tegmark et al. (1997) calculated the evolution of the H2 abundance, solving a system of kinet ic rate equations w i th in I D spherically symmetr ic models of halo profiles for different halo masses, and concluded that this fraction is close to 5 x 10~ 4 — 1 0 - 3 . A b e l et al. (1998b) performed fully 3D numerical simulations of the first bound objects, confirming that the gas can cool efficiently inside v i r ia l ized dark matter halos (which form inside spherical knots at the intersection of filaments) once the H2 fraction reaches 5 x 10~ 4 . Unfortunately, the spatial resolution in our R T models does not al low to resolve correctly regions of efficient H2 formation. These regions are barely resolved even i n the state-of-the-art, high-resolution hydro dynamical simulations ( A b e l et al. 1998b). A variety of techniques has been suggested to implement starbursts i n numerical simulations. To bypass rigours of not very well constrained physics of star formation, H a i m a n & Loeb (1997, H L 9 7 hereafter) s imply assume that a constant fraction / s t a r of the v i r ia l ized gas is converted into stars, cal ibrat ing this value against the inferred carbon abundance i n the I G M from observations of the L y a forest. The central assumption here is that heavy elements detected i n high-redshift L y a absorbers are produced by the same generation of stars that reionized the Universe. However, the estimate based on the rate of carbon ejection into the I G M presents a number of uncertainties (HL97) , since it is not entirely clear how much carbon gets actual ly ejected into the surrounding medium, or exactly what fraction of carbon burns into nitrogen at the base of convective envelopes of stars w i t h masses i n the range 3 — 8 M Q ("hot bo t tom burning") . In addi t ion, carbon is produced most ly by 3 — 6 M Q stars (HL97) , while most ioniz ing photons come from 10 — 12 MQ stars — hence, due to different lifetimes of stars w i t h different masses, there must be a time delay between stellar reionizat ion and metal enrichment of the I G M . H L 9 7 take the standard model as / s t a r = 4%, w i t h i n a broader possible range 1% < / s t a r < 99%. More recently, H a i m a n & Loeb (1998) (and later O h 1999) adopt a Chapter 5. Modelling Inhomogeneous Reionization 98 starburst model which is normalized to the observed metal l ic i ty range 10~ 3 < Z/ZMQ < 10~ 2 of the I G M at z ~ 3, inferred from carbon and sil icon abundances, y ie ld ing the star formation efficiency 1.7% < / s t a r < 17%. Note that the observed value of / s t a r i n star forming regions (defined as the fraction of H2 gas converted into stars) falls i n the range from ~ 1% to 30% (C ia rd i et al. 1999). C i a r d i et al. (1999) introduce an addi t ional parameter, / b , which represents the fraction of already vir ia l ized baryons which are able to cool and form stars, taking the average value fh = 8% (given by the 3D numerical simulations of star forming regions i n A b e l et al. 1998b). In their models, / s t a r = 15%, which we assume in our simulat ion. F ina l ly , not al l U V photons produced i n a star forming region can escape it easily. T h e escape probabi l i ty / e s c of ionizing photons (HL97) depends on a number of factors: the size and density of the ionized HII region associated w i t h the site of star formation, the dumpiness of the gas w i t h i n that region, the recombination timescale, and the rate of heating (by bo th U V and winds from young stars or SNe) . H L 9 7 calculated I D ionizat ion profiles of halos of different masses and redshifts, finding that / e s c is almost uniquely determined by the redshift, and giv ing an approximate fitt ing formula log / e sc = < 1.92 1510 — 1 for z > 10, (5.13) 0 for z < 10, which we adopt. Aga in , the value of / e s c inferred from observations of nearby star forming regions is not very well established, but rather falls into a range from few % to ~ 20% (see H L 9 7 and C i a r d i et al. 1999 for references). The to ta l stellar component inside a vi r ia l ized halo of mass M b can then be wr i t ten as (C ia rd i et al. 1999) M s t a r = M b / b / s t a r - (5.14) Chapter 5. Modelling Inhomogeneous Reionization 99 We have attempted to perform a R T run w i t h star formation inside the volume. We have assumed that once conditions for star formation are met (the vir ia l ized cell reaches the threshold H2 fraction), there is an instantaneous starburst w i t h a Salpeter i n i t i a l mass function. We have approximated the composite stellar spectrum (derived i n C i a r d i et al. (1999), for a metal l ic i ty Z/ZQ = 10~ 2 ) w i t h a tr iple power law, ag = (0.95,5,12). Immediately after the outburst at z = 30, the spectrum is nearly flat between the H I and He II ionization edges, becoming much softer on a timescale of ~ 10 7 yrs. A t t ime t = 0 the luminosi ty at the hydrogen L y m a n edge can be wri t ten expl ic i t ly as (Cia rd i et al. 1999) M * ) « 4.8 x 1 0 » ( ^ - ) (4) (&g) e r g H . - . (5.15) O f these photons, only L „ ( z ^ ) / e s c escape into the I G M . We adopt the values /b = 8%, / s t a r = 15% and / e s c = 2%. Unfortunately, the numerical resolution of our model d id not allow us to form any significant fraction of molecular hydrogen above the background value of 3 x 10~ 6 . A d d i t i o n a l complicat ion arises from the fact that the first stars can form wi th in ~ 10 5 M© clouds as early as z ~ 30 — 50. Since stellar spectra are much softer than those for quasars, star formation might be able to easily shut itself off due to photodissociation of H2, without leading to any significant hydrogen ionizat ion. Thus reionization might proceed v i a sporadic starbursts occurr ing i n different parts of the volume and separated by many redshifts units. We started our model at z = 30 observing only margina l star formation. We conclude that the present spatial resolution (64 3) is not sufficient for model l ing stellar reionization w i t h sources inside the computat ional volume v ia full R T . Chapter 5. Modelling Inhomogeneous Reionization 100 5.6 Observational Properties of Reionization A k i n to the density and temperature fluctuations i n the I G M at h igh redshifts (1 < z < 5) responsible for the L y a absorption forest i n the spectra of background quasars, the density, temperature and ionizat ion structure fluctuations at ul t ra-high redshifts (z > 5) should be visible i n absorption (either against more distant objects, e.g., quasars, galaxies or 7-ray burst afterglows, or against the C M B ) or emission above the C M B . Below we have l isted the major possibly detectable manifestations of the epoch of reionization: • atomic hydrogen 21cm emission (Hogan & Rees 1979, Oort 1984, Scott & Rees 1990, M a d a u et al. 1997, Gned in & Ostriker 1997, Tozzi et al. 1999, M e i k s i n 1999, Shaver et al: 1999); • L y a emission dur ing patchy reionization (Gou ld & Weinberg 1996, Mi ra lda -Escude & Rees 1998, Loeb & R y b i c k i 1999); • H a emission prior to and dur ing reionization (Oh 1999); • complete blanketing of a l l flux shortward of A L y ] j m j t ( l + z) i n spectra of a l l quasars beyond some cr i t ica l redshift, as well as observation of the damping wing of L y a absorption i n the neutral I G M around the same redshift; • ' radio-lines' of atomic hydrogen; • damping of p r imary C M B fluctuations and generation of secondary anisotropies; • free-free emission from ionized halos i n L y m a n l imi t systems (Af (H1 ) ~ 1 0 1 7 c m - 2 ) du r ing and after reionization (Loeb 1996, H a i m a n & Loeb 1997, O h 1999); and • the absorption properties of He II clouds Chapter 5. Modelling Inhomogeneous Reionization 101 Huge datasets which are going to be available soon (see Tables C . l - C.2 i n A p p e n d i x C for the list of currently implemented or planned observational programs) w i l l require comparison to high-resolution quantitative theoretical models. Absorp t ion or emission features from the high-redshift gas represent a convolution of 3D distr ibutions of different quantities, such as the temperature, the fractional ionization, the absorption coefficient at a specific frequency, and so on. In case of the inhomogeneous I G M , cosmological hydrodynamical simulations coupled w i t h 3D R T models present an accurate way of predict ing observable skymaps. In the next section, for the example of the 21-cm line radiat ion from atomic hydrogen, we demonstrate how these skymaps can be generated. 5.7 Atomic Hydrogen 21cm Emission Very-high-redshift sources (z > 5) can be detected through emission or absorption of 21-cm radia t ion by the surrounding neutral I G M . The physics of the 21-cm emission and absorption in a diffuse med ium was well established i n the 1950s (see, e.g., Wouthuysen 1952, F i e l d 1959). A number of researchers at the t ime were looking for 21-cm features from the intergalactic gas. Sunyaev & Zeldovich (1972) were the first to suggest that a large fraction of the gravi ta t ional energy released dur ing collapse of protoclusters of galaxies i n the "pancake" scenario w i l l result i n shock heating of the I G M , w i t h the spin temperature Ts of the gas being dr iven above the background radiat ion temperature T C M B J potential ly leading to detectable 21-cm radia t ion from neutral hydrogen. Reheating by shocks and light from the first collapsed objects was further investigated by Sunyaev & Zeldovich 1975, Hogan & Rees 1979, Oort 1984, Scott & Rees 1990, and recently by M a d a u et al. 1997, Gned in k Ostriker 1997, Tozz i et al. 1999, M e i k s i n 1999, Shaver et al. 1999. It has been shown that it is possible to detect sources appearing before reionization, through the reheating of the neutral I G M and subsequent emission or absorption against the C M B spectrum at the rest wavelength of the atomic hydrogen 21 c m line. Due Chapter 5. Modelling Inhomogeneous Reionization 102 to inhomogeneities i n the gas at high redshift these spectral signatures w i l l result i n angular fluctuations across the sky, as well as showing some structure i n redshift space (Hogan & Rees 1979, Scott & Rees 1990). T h i s potential ly opens up the possibil i ty of reconstructing the 3D dis t r ibut ion of gas dur ing patchy reionization. 5.7.1 21-cm Emission and Absorption T h e 21-cm line i n atomic hydrogen corresponds to the transi t ion between the singlet and triplet n = 1 hyperfine states which exist due the electron-proton spin coupling. A tr iplet a tom radiates spontaneously w i t h probabi l i ty A i o = 2.85 x 1 0 ~ 1 5 s - 1 (Fie ld 1959). The intensity of the 21-cm radia t ion depends on the spin temperature defined by the Bo l t zmann equation s imply as the exci ta t ion temperature of the hyperfine levels (e.g., M e i k s i n 1999) ! i i = 3 exp ( , T* = « 0.068K, (5.16) where no and ri\ are the populations of the upper and lower states, and the numeric coefficient is just the ratio of their statist ical weights. The coupling between T s and the temperature of the background radiat ion TC M B is pretty tight, so that i n the absence of any other pumping mechanism T s - » TC M B on a timescale T* / (TC M BA I O ) ~ 2 X 1 0 4 y r s (Meiks in 1999), and there w i l l be no 21-cm absorption or emission features on top of the C M B . Fortunately, there are at least two mechanisms which can decouple the spin temperature from the C M B (F ie ld 1959; also, see F i g . 5.14). F i r s t of a l l , i n dense regions collisions between hydrogen atoms can change the col l is ional / radiat ive balance between the two levels, essentially d r iv ing the spin temperature T s towards the gas temperature T^. Before there are any heating sources at z <> 30 — 50, the coll is ion-induced decoupling of T S and TC M B w i l l be the only mechanism leading to 21-cm spectral features on top of the C M B (Scott & Rees 1990). However, its efficiency is very smal l Chapter 5. Modelling Inhomogeneous Reionization 103 and might be important only for the densest regions at Sp/p J> 30[(1 + z)/10] 2 (Madau et al. 1997). A more l ikely mechanism is the exci tat ion to the triplet a tom by absorption of a L y a photon to the n = 2 level, a so-called Wouthuysen-Field effect (Wouthuysen 1952, F i e l d 1959). In this process a singlet hydrogen atom i n the lrjSi/2 sta'te (where the notat ion is npLj, w i th L being the electron orbi ta l angular momentum, J = S + L the total electron angular momentum, and F = (proton spin) + J the total angular momentum) absorbs a L y a photon and makes a t ransi t ion to one of the n — 2 states, 2 i P i y 2 or 2i_P 3 / / 2 , decaying i n some cases to the tr iplet 1] S\/2- The opposite reaction, the de-excitation of the hyperfine triplet state through scattering of a L y a photon is also possible. Since i n the neutral medium the optical depth i n the L y a line is very large, mult iple scatterings of a L y a photon drive the temperature i n the line close to the kinetic temperature of the gas Tk, and, accordingly, the spin temperature Ts might be par t ia l ly coupled to T^. W i t h these two mechanisms, the spin temperature can be wri t ten as (Fie ld 1959) rp _ TCMB + {yc + ya)Tk i + yc + ya where yc is the coll is ional coupling parameter (see Scott & Rees 1990 for a f i t t ing formula) C T, y. Aw 2 1 ' w i t h C being the coll isional rate, and ya is the L y a pumping coefficient (Fie ld 1959) (5.18) Aw T k ' w i t h Pw being the de-excitation rate of the triplet state v i a absorption of a L y a photon to the n = 2 level. Chapter 5. Modelling Inhomogeneous Reionization 104 5.7.2 High-z Absorption Features A t very high redshifts, before any discrete sources of heating, the deviat ion of the spin temper-ature from the background radiat ion temperature is entirely due to collisions between neutral atoms (Scott & Rees 1990). However, this effect is very smal l (see F i g . 5.15), wi th Ts « TC M B , and one would not expect to detect any absorption features on top of the C M B . After the first sources of reheating appear, L y a photons w i l l de-excite the upper hyperfine (triplet) state, lowering the spin temperature closer to the gas temperature. D u r i n g this short epoch Tk < T s <C TC M B , y ielding stronger absorption features. Accord ing to an estimate i n Tozz i et al. (1999), the I G M w i l l be visible in the 21-cm absorption over a timescale of ~ 10 7 yrs, wi th brightness temperature variations of order A T b ~ 40 m K over a frequency range of ~ 5 M H z . 5.7.3 High-z Emission Features Very soon (~ 1 — 10 Myrs ) after the I G M becomes visible i n absorption, the heat input from the first cosmological objects w i l l raise the gas temperature above that of the C M B . L y a photons of increasing number density w i l l couple the spin temperature Ts to the kinet ic temperature Tk, d r iv ing the former considerably above the temperature of the background radia t ion and rendering neutral hydrogen observable i n 21-cm emission. Let us see how the L y a pumping by photons originating close to H I I regions (which i n tu rn are dr iven by U V and soft X-rays from the first cosmological sources) results i n angular and redshift-space structures. We determine the temperature of the gas by comput ing the R T of photons above the hydrogen L y m a n edge. In reality, the thermal state of the I G M ahead of the hydrogen I-front is strongly affected by photons below 13.6 e V . Here, we neglect transfer i n L y a and other lines, effectively, lowering the temperature of the neutral medium. However, we s t i l l have significant preheating i n our models. The higher energy (above 13.6 eV) U V and soft Chapter 5. Modelling Inhomogeneous Reionization 105 X - r a y photons propagate further i n neutral hydrogen since the opacity drops w i t h frequency as oc and they are able to heat the gas ahead of the I-front very efficiently. In radioastronomy, the specific intensity Iu is often replaced by the frequency-dependent brightness temperature T b ( ^ ) , which in the Rayleigh-Jeans approximat ion (hu <S U-QT) is pro-por t ional to the intensity: T h e n the R T E • *k = -lv + Bv{T) (5-20) through a med ium w i t h some physical temperature T of emission can be wr i t ten as ^ = - T b + T . (5.21) drv If T = constant along the line of sight, the solution for the brightness temperature is T b = Thoe-T" + T ( l - e - T - ) . (5.22) For a high-redshift neutral cloud, T b 0 = TOMB, T = TS and TV = r is the total 21-cm line opt ical depth across the region. The resulting brightness temperature can be lower or higher than that of the C M B , hence, the cloud can be visible either i n absorption or emission against the background radiat ion. The differential temperature between the 21-cm signal and the C M B is ' A T = T b - T C M B = (T s - T C M B ) ( 1 - e~T). (5.23) We w i l l assume that there are enough Lya : photons to couple the spin temperature T s to the kinetic temperature of the gas Tk throughout the entire volume, and that T s 3> T* w 0.068 K Chapter 5. Modelling Inhomogeneous Reionization 106 everywhere. In this case, the ratio of the atomic hydrogen densities i n the upper and lower hyperfine states no go is equal to the ratio of their statistical weights, gi/go = 3, independently of the spin temper-ature. T h e n the fluctuations i n the 21-cm radiat ion w i l l be caused only by the gradient of the local neutral hydrogen density, shaped by the edges of H n regions. The 21-cm absorption coefficient is c2Aw 1. 9 l ( g0\ 3c2Aw 1 T* « = g 2 2 A re0 - n x — ~ Q 2 2 A ~ n ° ^ r " ( 5 - 2 4 ) STT^U^ Au go \ giJ STTZU^ AU 1s (see, e.g., Scott & Rees 1990). In this equation, Au is the effective w i d t h of the 21-cm line profile (approximated by a rectangle), and A%o is the spontaneous rate for the 21-cm transi t ion. Since n i / n o —> 3, the lower level popula t ion is no = (1 — xe)nn/4, where n n is the number density of atomic hydrogen, and xe is the fractional ionizat ion. The 21-cm line opt ica l depth across the cloud, r = K X (path length), is approximately r « 2.59 x 1 0 ~ 1 9 c m 2 H z K ^ H 1 ) , (5.25) AuTs where 7Y(H I ) is the neutral hydrogen column density along the line of sight. For smal l opt ical depths ( r 1), the differential temperature becomes A « T s ~ r ° M B ffiLl x 2.59 x 1 0 - 1 9 c m 2 Hz K . (5.26) % Au y ' W i t h i n a large cloud of geometrical depth I, the effective line w i d t h is most strongly affected by the Hubb le flow velocity gradient Av across the region, and can be wri t ten schematically as Au RJ u2icm — , where Av = H0(l + z)3/2l, c Chapter 5. Modelling Inhomogeneous Reionization 107 whereas the co lumn density is n H ( l - xe)l3 iV(Hl I2 T h e final expression for the differential brightness temperature (its sign indica t ing absorption or emission against the C M B ) takes the form A T ^ T ^ - T C M B n H ( l - x e ) C _ x 2 ^ x 1 Q _ 1 9 ^ H z T s ^21cm#o(l + Z)b'1 where the addi t ional 1/(1 + z) factor accounts for the change i n the observed brightness tem-perature due the expansion of the Universe. O u r numerical simulations give the 3D distributions of the fraction of hydrogen ioniza t ion and temperature. In order to. obtain the observed profile of 21-cm radiat ion on the sky, we integrate the brightness temperature along mult iple lines of sight, solving the transfer equation dTh = (T s - Th)dr, w i t h T b ; i n = T C M B , and obtaining the signal ( T b j O U t - T C M B ) / ( 1 + z). Figures 5.16 - 5.17 show a skymap of 21-cm emission originat ing from the wa rm neutral I G M around the bright quasar i n M o d e l B'. We assume that an efficient pumping mechanism (e.g., L y a photons) drives the spin temperature of the neutral gas very close to the kinet ic temperature everywhere i n the volume. T h i s temperature, as well as the ionizat iohal structure of the I G M are calculated v i a detailed R T above 13.6 e V . T h e adopted power law for quasar radiat ion, w i t h the spectral index a g = 1.8 i n a l l frequency groups, results i n a relatively flat spectrum, quickly raising the temperature of the neutral I G M . T h e 21-cm signal on the sky is a convolution of the gas temperature and the density of atomic hydrogen, showing an abrupt cut-off near the edge of the expanding H I I region. Note, that other theoretical maps of 21-cm emission from the inhomogeneous medium at a h igh redshift have been recently published (Meiks in 1999, Tozz i et al. 1999). However, to the best of our knowledge, none of them include detailed numerical R T as yet, but are based on an approximate estimate of the radius of a I D Stromgren sphere which is superimposed on the Chapter 5. Modelling Inhomogeneous Reionization 108 exist ing 3D density dis t r ibut ion. In this thesis, we have shown that such calculations can be performed w i t h detailed time- and frequency-dependent numerical R T . 5.8 Discussion and Conclusions T h e speed of I-fronts depends, among other variables, on the ambient neutral hydrogen density. Thus , the timescale on which the ionizing background 'eats i n ' to dense clumps and filaments is determined by the typical density contrast i n the volume. T h e latter quanti ty depends crucial ly on the spat ial resolution. In the simulations presented i n Sec. 5.5.1 - 5.5.3, the predominant density contrast is of order ~ 100 at 64 3 , and closer to ~ 50 at 32 3 . Note that the method developed i n this thesis does not treat R T on scales smaller than one cell. Therefore, one has to pay special attention to modell ing the photochemical evolution of, say, a single overdense cell submerged i n a bath of ionizing photons. In the present study, we s imply take out the number of photons absorbed wi th in that cell from the radiat ion field w i t h i n the chemistry solver, given that this field w i l l recomputed again at the next t ime step, after one light-crossing time of the cell . However, imagine a cell harbouring a dense halo. The overall transmission properties of this cell (frequency-dependent opacity and emissivity) w i l l differ from those of the same volume element i f the halo is not resolved. Here, we can only speculate that, depending on the astrophysical context, it might be possible to introduce a self-shielding fudge factor accounting for c lumping of the gas on sub-resolution scales. In this Chapter , we have performed several runs of cosmic reionizat ion by bo th point sources and the background radiat ion field. We solve self-consistently the coupled equations of radiative transfer and non-equi l ibr ium chemistry. It is shown that i n general the 3D d i s t r ibu t ion of ionized regions follows the gas dis t r ibut ion. In a l l models low density voids are ionized first, and dense filaments are destroyed much more slowly, on the timescale determined by the ratio of the photoionizat ion to recombination rates at those high densities. Chapter 5. Modelling Inhomogeneous Reionization 109 In models w i t h quasars inside the volume, the spectrum of the direct radia t ion does indeed become harder w i t h t ime due to the wavelength dependence of the opt ical depth. We have not included any Doppler shifts into the model . However, having per iodic boundary conditions allows direct ioniz ing photons from quasars to travel over many computat ional volumes and experience the Doppler shift due to the Hubble flow, given that these photons are not absorbed before they get redshifted. T h e Hubble speed on the scale of 2.5 M p c is 2 5 0 / z k m s - 1 . E v e n w i t h periodic boundary conditions, the error due to neglect of the Doppler shift is not large compared to other sources of uncertainties. The most advanced s imulat ion is model A , wh ich has been evolved from z = 15.0 to z = 12.3, w i t h sl ightly under 4000 timesteps. I f we assume that there are some direct high-frequency photons which manage to travel over the entire length of the model without being absorbed, they cross the volume ~ 60 times, which corresponds to the Hubb le speed of ~ 15 x 1 0 3 / i k m s - 1 . T h e resulting error is most certainly smaller than 5% or 10% at most and does not exceed the error coming from the replacement of frequency-dependent quantities w i t h integrals inside each frequency group. In addi t ion, par t ia l ly due to the l imi ted frequency resolution above 54.4 e V , the vast majori ty of direct photons gets absorbed i n the neutral med ium (which s t i l l occupies most of the volume i n model A at the end of the s imulat ion) . Hence, the error resulting from the neglect of the Doppler shift due to the Hubble flow is much smaller than other uncertainties i n the model . In future work, we w i l l increase the frequency resolution to calculate the propagation of high energy ( U V and X-ray) photons w i t h much greater care. Since these photons can travel much further i n a neutral medium, they might be largely responsible for the preheating of the I G M , as well as for "blurr ing" I-fronts. We have made no attempt to put observational constraints on different scenarios of reion-izat ion, instead, concentrating on developing tools to accomplish this task i n the future. O n the example of 21-cm line emission i n Sec. 5.7, we have demonstrated that fu l l , RT-enabled simulations of this k i n d w i l l produce maps which could quanti tat ively show the observational signatures of reionization. Chapter 5. Modelling Inhomogeneous Reionization 110 In future work, we w i l l increase the resolution of these simulations to 128 3 or 256 3 , w i t h the corresponding update i n angular resolution. Time-dependent ray t racing (as well as any explici t R T solver) is wel l suited to clusters of dis t r ibuted memory machines, since remote volume elements cannot communicate on timescales shorter than the light-crossing time, and therefore can be followed independently. Chapter 5. Modelling Inhomogeneous Reionization 111 Figure 5.1: The isosurfaces of the density d is t r ibut ion i n a 2.5 M p c (comoving) volume used for the R T model i n the present work are plotted at the overdensity S = 2 at eight different redshifts. The density field was k ind ly supplied by M . Whi t e . Chapter 5. Modelling Inhomogeneous Reionization 112 Figure 5.2: Th i s figure shows a cross-section through a 3D hydrogen I-front propagating into the neutral I G M from a vi r ia l ized halo hosting a mini-quasar around z = 15, at four different output times. T h e quasar luminosi ty is 4.2 x 1 0 5 3 s _ 1 M h a i o / 1 0 9 M 0 , and its lifetime is 1 0 5 8 y r s (model A , see text i n Sec. 5.5.1 for details). The s imulat ion employs the full 3D radiative transfer i n three frequency groups (for bo th direct ionizing and diffuse photons) combined wi th the solution of the rate equations for nine chemical species, at the spatial resolution 64 3 . The size of the volume is 2.5 M p c (comoving). The quasar turns off around z0g ~ 14.92, and after that the I-front is driven entirely by secondary (recombination) photons coming from the H I I bubble. To compare the I-front speed to the speed of light, it 's wor th point ing out that by the t ime of the lower right corner plot, the wavefront would have crossed the volume four times. Chapter 5. Modelling Inhomogeneous Reionization 113 Figure 5.3: Same as F i g . 5.2, but for singly ionized hel ium. Note that the H e l l region is a relatively th in shell bounded by neutral hel ium i n the ambient medium and by doubly ionized hel ium closer to the quasar. Chapter 5. Modelling Inhomogeneous Reionization 114 Figure 5.4: Temperature for a cross-section through the I G M i n model A at z = 14.785. The quasar is at the upper right corner. T h e long spikes going outward from the quasar are separated by shadows forming behind denser regions (these spikes have nothing to do w i t h the finite angular resolution displayed i n F i g . 2.2). Hydrogen ionizat ion corresponds to the temperature range (3 — 4) x 10 4 K . Chapter 5. Modelling Inhomogeneous Reionization 115 Figure 5.5: A cross-section through a 3D hydrogen I-front at the beginning of the s imulat ion for M o d e l B , at four different output times. The quasar (in the upper left part of the plot) turns on at z = 10 and emits ionizing photons a constant luminosi ty 1 0 5 4 s - 1 throughout the entire run. Chapter 5. Modelling Inhomogeneous Reionization 116 Figure 5.6: Same as F i g . 5.5, but for singly ionized hel ium. The H e l l region is bounded by the He i n bubble on the inside, and by neutral hel ium further away from the quasar. To better illustrate the structure of the front, the fraction of ionizat ion from the lower right plot is displayed i n color i n F i g . 5.7. Chapter 5. Modelling Inhomogeneous Reionization 117 Chapter 5. Modelling Inhomogeneous Reionization 118 Figure 5.8: Th i s figure displays the energy density of the radiat ion field E G i n three frequency groups (g = 1,2,3 top to bottom) vs. fractional ionizat ion of hydrogen and helium, in a cross-section through M o d e l B at z = 9.473. Chapter 5. Modelling Inhomogeneous Reionization 119 Figure 5.9: Th i s figure illustrates hydrogen reionization by a U V background field (model C , see Sec. 5.5.2 for details) start ing at z = 13. The spatial resolution is 32 3 , the box size is 2.5 M p c . The speed of complete ionizat ion is a function of the spatial resolution which affects the max imum density contrast w i t h i n the volume. A t 32 3 , the densest neutral filaments survive at least un t i l z ~ 12.6 when we stopped the s imulat ion. Figure 5.11 Cont inued from F i g . 5.10. Chapter 5. Modelling Inhomogeneous Reionization 122 Figure 5.12: 3D dis t r ibut ion of neutral (atomic and molecular) hydrogen from the simulat ion on F i g . 5.9 (ionization by the U V background, model C ) . Note that this is not a plot of the integrated co lumn density but rather a 3D rendering of the actual neutral hydrogen number density. The colorbar on the left of each plot is approximate, since it is affected by shielding i n the visual izat ion routine, and is given here only for reference. Chapter 5. Modelling Inhomogeneous Reionization 123 Figure 5.13: T h i s figure shows a cross-section through the hydrogen I-front in M o d e l C com-puted at three different spatial resolutions (8 3 x 8 2 , 16 3 x 8 2 and 32 3 x 8 2 , top to bottom), at the same output t ime. To ease the comparison, the underlying density field is visualized at 3 2 3 for a l l three models. Chapter 5. Modelling Inhomogeneous Reionization 124 adiabatic expansion ~ (1+z) -CMB Zt collisions (important in dense clouds at lower z) Is X L y alpha photons adiabatic expansion: no direct effect K adiabatic expansion — ( l+z) A 2 Figure 5.14: Energy contributors to the spin temperature TS i n an adiabatically expanding universe. Arrows indicate the direction of energy transfer. In the absence of coll isional and L y a pumping the spin temperature comes to equi l ib r ium w i t h the C M B temperature on a timescale T*/(TCMBAW) ~ 5 X 1 0 4 y r s (Meiks in 1999). Chapter 5. Modelling Inhomogeneous Reionization 125 T «T <T A K S CMB very close, no 21-cm features TK < Ts <K XCMB ... bigger gap, 21 -cm absorption ' larger coupling T <T<T andT ~T CMB S K -""CMB K very close, no 21-cm features < < anH « 21-cm emission CMB S K CMB K disappearance of neutral H time Figure 5.15: Emergence of 21-cm emission or absorption features at h igh redshifts. Chapter 5. Modelling Inhomogeneous Reionization 126 Figure 5.16: Temperature of the I G M i n a cross-section through M o d e l B (upper panel), at two output times, and the resulting map of 21-cm emission against the C M B (lower panel), i n terms of the brightness temperature. The spin temperature of the 21-cm transi t ion is assumed to be the same as the kinet ic temperature of the gas throughout the entire volume, which is close to the temperature of the C M B far from the quasar. We do not take into account the effects of R T below 13.6 e V . In addi t ion, the adopted spectrum of the quasar is rather flat, resulting in an artificially high temperature i n the neutral I G M ahead of the hydrogen I-front. Chapter 5. Modelling Inhomogeneous Reionization 127 Figure 5.17: Cont inued from F i g . 5.16. The snapshot at z = 9.967 (the column on the right) is close to complete overlapping of HII regions. Chapter 6 Concluding Remarks In this thesis, I have described a method for the numerical solution of the frequency- and time-dependent radiative transfer equation i n a three-dimensional inhomogeneous medium. It is demonstrated that, w i th existing desktop hardware, it is possible to model cosmological inhomogeneous reionization on a light-crossing time £R i n three spatial dimensions. Since the photoionization time-scale in the small optical depth regime (T < 1) is of order of the light-crossing time explicit advection is an efficient method in covering at least these regions. Compared to el l ipt ic-type solvers on the fluid-flow time-scale or the time-scale of typ ica l changes i n the density and ionization structure dis t r ibut ion, explici t techniques produce very accurate results without the need to solve a large system of coupled non-linear el l ipt ic equations. The comput ing requirements wi th explicit advection grow linearly wi th the inclusion of new atomic and molecular rate equations, which is certainly not the case for quasi-static solvers (although it is feasible that the development of mul t igr id techniques for el l ipt ic equations might actually approach s imilar scaling). Us ing eq. (2.18) one can see that the entire history of reionization can be modeled w i t h ~ 1 0 4 - 1 0 5 time-steps (depending on the required resolution), which makes explici t advection an attractive choice for these calculations. However, the efficiency of this method has s t i l l to be explored. Future work should include a detailed comparison between explici t advection and impl ic i t reconstruction (through an ell ipt ic solver), to demonstrate which method works best for calculat ing inhomogeneous reionization. A s I have demonstrated here, for certain problems, inc luding the propagation of supersonic 128 Chapter 6. Concluding Remarks 129 I-fronts, the Courant condi t ion does not seem to impose prohibi t ively smal l time-steps. In this case, the biggest challenge is to accurately describe anisotropies i n the radiat ion field, that is, to solve for inhomogeneous advection i n the five-dimensional (without frequency) phase space, i n the presence of non-uniform sources and sinks of radiat ion. St r ic t ly speaking, the storage of one variable at, say, 64 5 data points requires about 9 G B of memory, which stretches the capabilit ies of top-end desktop workstations. One attractive possibil i ty for future explorat ion is to direct ly solve the monochromatic photon Bo l t zmann equation i n 5D. To demonstrate the feasibility of the numerical solution, however, among different methods, I have chosen to concentrate on simple ray tracing at the speed of light. The numerical approach being used is completely conservative and produces very l i t t le numerical dissipation. I have computed a series of reionization models for different strengths and spectra of the background radiat ion, wi th in a fixed cosmological framework. To be able to dis t inguish between different scenarios of the first object formation, I w i l l attempt to go to a higher spat ial resolution (128 3 or 256 3 ) in the future work, as well as include explici t radiative transfer below 13.6 e V dur ing preheating. O u r simulations show that, depending on the luminosi ty of the first objects, expanding H l l regions could overlap very soon after these sources tu rn on (model B ) . In models w i t h lower photon product ion rates (models A and C ) , the complete overlapping happens later, a l though the densest clouds in the I G M stay neutral by the time we stop calculations. It is wor th point ing out that al though HII regions completely overlap at z > 5 (the exact epoch depends on the scenario of reionization), there is s t i l l a significant fraction of residual neutral hydrogen i n dense L y a clouds and Lyman- l imi t systems. M a d a u et al. 1999 define the so-called breakthrough redshift at z^ ~ 1.6 when the Universe becomes opt ical ly t h in to Lyman-con t inuum photons. In future work we shal l evolve our models to much lower redshifts to study how z D r is affected by the details of reionization at earlier stages. Time-dependent ray tracing (or another explicit advection technique) could be used to solve Chapter 6. Concluding Remarks 130 the problem of R T in an inhomogeneous medium in many other areas of astrophysics, especially, when the Courant condi t ion does not result i n prohibi t ively smal l time steps. E x p l i c i t techniques w i l l work well bo th for monochromatic (line) radiat ion and for frequency-integrated quantities. Examples of problems w i t h R T i n ind iv idua l lines include cooling i n star forming regions, preheating at h igh redshifts, and L y a emissivity of par t ia l ly neutral clouds submerged in a bath of ioniz ing photons. The global exchange of energy v i a the radiat ion field at a fixed wavelength is probably one of the easiest problems to solve i n numerical R T . In the cosmological context, it seems l ikely that the progress i n computer technology, a few years from now, w i l l allow routine solut ion of this problem in three spatial dimensions w i th the sort of resolution obtained i n modern, state-of-the-art hydro dynamica l simulations (10 7 -10 9 three-dimensional data points). O n the other hand, the full solut ion of the R H D equations, retaining a l l 0(v/c) terms, is a much more com-plicated problem. In this case one deals wi th frequency-dependent R T , and issues such as line transfer, broadening effects and spatial motions w i th in the s imulat ion box become important . Nevertheless, despite the high dimensionality of the problem, w i t h a reasonable expenditure of computat ional resources (of the type available today), it is possible to numerical ly model many different aspects of the full 3D radiative transfer problem. The method described i n this thesis represents a significant and realizable step towards the goal of full cosmological R H D . Numer ica l model l ing of cosmological R T w i l l need to continue to develop i n the near future, i n concert w i t h the rapid observational progress in understanding the end of the "Dark Ages" . In conclusion, i n Tables C . l - C.2 I have listed some of the major observational programs currently being planned or already implemented which could potential ly probe the reionizat ion epoch and the formation of the first objects. Appendix A High-Resolution Schemes for One-Dimensional Hyperbolic Systems The goal of this appendix is to describe a high-resolution numerical scheme for the solut ion of I D non-linear advection problems. The linear version of this a lgor i thm is used to calculate first-order flux updates and higher-order corrections for normal , transverse and corner waves i n the numerical Bo l t zmann solver i n Sec. 2.8. The full non-linear version is appl ied to the derivation of an unsplit 3D moment solver i n Append ix B . F i rs t , consider a linear system qt + Aqx = 0, q=(Qi,Q2,...,QN)T ( A . l ) (LeVeque 1997) of source-free hyperbolic equations in I D . B y "linear" we mean that the operator A does not depend on the coordinate x. To find the solution of eq. ( A . l ) , we have to separate variables Qi,. •., QN, i-e. to diagonalize the mat r ix A. We begin by in t roducing a new variable q = Pu, so that eq. ( A . l ) becomes ut + P~lAPux = 0. (A.2) We require that P~lAP = A is a diagonal matr ix . In other words, A is the diagonal mat r ix composed of eigenvalues Aj of A , and P consists of right eigenvectors (pp) of A . E q . (A.2) yields qt + PAP~lqx = 0. (A.3) 131 Appendix A. High-Resolution Schemes for One-Dimensional Hyperbolic Systems 132 We continue by d iv id ing A into positive and negative parts A = A + + A " , A ± = d i a g ( A ± , . . . , A ± ) , A+ = max(A,0) and A ~ = min (A,0 ) . (A.4) Since A is linear, we can also write A± = PA±P ± D - 1 The quanti ty A+Aqi_1/2 = A+(qi - qi-x) = PA+P^Aq^^ = PA+AUi_1/2 = / XfAUl ^ = P A+AUAT 0 (A.5) (LeVeque 1997) where Aup is a linear combination of A Q 1 ) j _ 1 / 2 , • • • ^ Q j v , j - i /2 i c a n be thought of as the difference between all right-going waves to the left (i — l —> i) and to the right (i —> i + 1) from point xi (waves are numbered (5 = 1 , . . . , N). The eigenvectors A^" are s imply the wave speeds, Aup the wave amplitudes and pp the wave spectra. The ful l Godunov ' s update to the numerical solution of eq. ( A . l ) is then q"+1 =q?~l [ A + A ? i - l / 2 + ^ " A f t + l / 2 ] • (A.6) Now, instead of eq. ( A . l ) consider a non-linear system qt + A(x)qx = 0, q={Qi,Q2,...,QN)T- (A-7) The wave speeds and spectra now vary w i t h x, and eq. (A.5) is not val id any more. Instead, we can write Appendix A. High-Resolution Schemes for One-Dimensional Hyperbolic Systems 133 A+Aqi-i/2 = E A / W * - i / 2 (P0i) > 0 A'Aqi_1/2 = E A ^ - l a / 3 * - l / 2 (P/3t-l) • ' ( A - 8 ) To get the correct wave amplitudes a ^ i - i ^ we have to use the proper j u m p conditions for the R i e m a n n problem at the interface between Xi-i and xf. E l Ppi-l P/3i I , . <*0i-l/2 t 0 r t | ' (A- 9) i.e., the total j u m p A g ^ ^ is just the sum of jumps across ind iv idua l waves. E q . (A.6 , A . 8 -A.9) ensure the full conservation for original hyperbolic systems; however, they provide only a first-order solution (LeVeque 1997). Appendix B A n Unsplit Method for 3D Moment Transfer In this appendix we develop an explicit (on the light-crossing timescale) solver of the first two radiat ion moment equations for Sec. 2.6 given a separate closure technique. Let us consider the advection part of the truncated system of 3D moment equations (2.7) - (2.8): ^ = -c2V-(fa0E). ( B . l ) Assume, that the Edd ing ton factor coefficients fafj{r) are already known from some closure scheme (the formal solution for the angle-dependent intensity). Cont ra ry to the 5D Bo l t zmarm equation w i t h constant coefficients (eq. (2.30)), the system ( B . l ) can be very non-linear, for instance, at the transi t ion layer between optical ly th in and opt ical ly thick regions. T h e solut ion of eq. ( B . l ) on a 3D rectangular mesh begins wi th the construction of the first-order accurate scheme. It is useful to write eq. ( B . l ) i n the matr ix form Qt + [A(x, y, z)q]x + [B(x, y, z)q] + [C(x, y, z)q]z = 0, where / E \ Fi F2 \F3J ( \ A = 0 1 0 0 c 2 / n 0 0 0 c 2 / 2 i 0 0 0 V c 2 / 3 i 0 0 0 J 134 Appendix B. An Unsplit Method for 3D Moment Transfer 135 B = ( 0 0 1 0 x c 2 / i 2 0 0 0 c 2 / 2 2 0 0 0 V c 2 / 3 2 0 0 0 / and C = 0 0 0 1 \ c 2 / i 3 0 0 0 c 2 / 2 3 0 0 0 V c 2 / 3 3 0 0 0 ) (B.2) The first-order normal update i n the ^-direct ion is given by the left- and right-going flux differences from eq. (A.8). The speeds of waves along the rr-axis are given by eigenvalues of the mat r ix A Xl,ijk = 0, A2,ijfe = 0, \-A,ijk — Cyjfn/tjk, ^4,ijk — —c\fhi,ijk The corresponding wave spectra (right eigenvectors of A) are (B-3) / 0 \ Pl,ijk V 1 J P3,ijk = ( \f fll,ijk cf 21,ijk fll,ijk f21,ijk 1 . f31,ijk \ f21,ijk J P2,ijk 0 1 V 1 I Pi,ijk — / \/hl,ijk \ cf 21,ijk fll,ijk f21,ijk 1 V fzi,ijk , f21,ijk ' The j u m p condi t ion (eq. (A.9)) then takes the form (B.4) &Qi-l/2,jk = al,i-l/2,jkPhiJk + a2,i-l/^3kP2,ijk + a3,i-l/2,jkPS,ijk + + ali-l/2,jkP4,i~lJ,k (B-5) yielding the strengths of the waves. Since waves #1 and #2 have zero velocities, the only amplitudes we are interested i n are Appendix B. An Unsplit Method for 3D Moment Transfer 136 f21,ijk ( • v / / l l , i - l , i f e^92 , ? : - l/2 , j f e + c / i i , i - i 1 j f c A Q 1 ) j _ 1 / 2 j f c ) 3,i l/2,jk y/ fn,ijkfll,i-l,jk + fll,ijk\/ fll,i-ljk f21,i-l,jk ( \ / / l l , i j f c ^ 9 2 , i - l / 2 j f c ~ cfn,ijk&<ll,i-l/2,jk) «4 , i - l / 2 j f e = / f f / f . • T h e first-order normal update for advection i n the x-direct ion is (B.6). = " I XhjkaP,i-l/2jkPl3,ijk + E Xf>jjk<$,i+l/2jkPP,iJkj • (B.7) N o r m a l updates along the y- and z-axes are performed i n a similar way. Appendix C Observing Programs to Detect the Epoch of Reionization Probing the epoch of reionization program description first l ight Sloan Digital Sky =4> five-color imaging and spectroscopic survey of one 1998, dura t ion Survey (SDSS) quarter of the sky ~ 5 years possibility of detecting up to a few hundred bright quasars and galaxies at z > 5, with the brightest objects in the redshift interval 5.5 <> z 7 http://www. sdss. org Chandra X-ray Ob- => 17' x 17' field of view; sensitivity threshold of ~ 2 x servatory (CXO) 1 0 - l 6 e r g s - 1 c m - 2 might detect of order a hundred quasars per field of view in the redshift interval 5 <i z <> 8; IR spectroscopy of X-ray selected quasars with N G S T could probe the reionization history up to z ° http .-//snail, msfc .nasa.gov/axaf/ax af.html 1999 Giant Me-trewave Radio Tele-scope ( G M R T ) 1996 (?; an array of 30 fully steerable 45m parabolic dishes spread over the area wi th a baseline of up to 25 k m (wi th a. compact central area enclosed i n a larger "Y"-shaped ! configuration) could detect 21-cm emission or absorption against the C M B from high-redshift neutral hydrogen, if reionization occurred at 6 z <> 10 http://158.144-11.2/page/grnrt/grnrt. html Square Kilometer Array (SKA) a project to bu i ld a giant radiotelescope w i t h the to- I ~ 2010 tal collecting area of at least a square kilometer; several different options currently under study could detect 21-cm emission or absorption against the C M B from high-redshift neutral hydrogen, if reionization occurred at z <J 15 http://www. nfra. nl/skai http://www. atnf. csiro. au/lkT Table C . l : Current and future observing programs potential ly able to probe the epoch of reionization (continued on the next page). 137 Appendix C. Observing Programs to Detect the Epoch of Reionization 138 Probing the epoch of reionization program description first light Next Genera-tion Space Telescope (NGST) • possibly, will allow the dete medium-resolution spectros Peterson optical depth TQP http://ngst.gsfc.nasa.gov http:/'/www. ngst. stsci. edv =>• expected imaging sensitivity in the infrared better than 1 nJy at wavelengths 1 — 3 pm ction of individual galaxies and mini-quasars directly to z ~ 10, copy to z ~ 8; in the case of gradual reionization might measur > 1; might detect Lya emission from the reionization epoch ~ 2008 and perform e the Gunn-M A P Planck all-sky differential map of C M B fluctuations i n five fre-quency bands from 22 to 90 G H z w i t h an angular resolu-t ion of ;> 0.21° - 0.93°; w i l l be launched on a stationary orbit around the Sun-Ear th Lagrange ( L 2 ) point =4> high-angular resolution mapping of the C M B over, the whole sky, simultaneously over a wide frequency range; designed to provide 10 times the sensitivity and more than 50 times the angular resolution of the C O B E satellite ~ F a l l 2000 ~ 2007 measuring small scale C M B anisotropies with sensitivity to the electron scattering optical depth as small as a few %; should be able to determine the exact (if any) redshift of reionization and probe secondary effects during patchy reionization http://map.gsfc.nasa.gov http://astro. estec.esa.nl/SA-general/Projects/Planck Table C.2: Cont inued from table C . l . Appendix D Glossary B H black hole C D M cold dark matter CFL, or the Courant condition Courant-Friedrichs-Levy condi t ion T h e stabil i ty condit ion in explicit advection techniques requir ing the time step to be less than the t ime taken for a sound wave (or an advection wave for general hyperbol ic equations of non-hydrodynamical nature) to cross one computa t ional gr id zone. C M B cosmic microwave background D M dark matter FD finite difference I-front ionizat ion front IGM intergalactic medium L T E local thermodynamic equi l ibr ium A n assumption that a smal l (local) volume element is i n thermodynamic equi-l i b r i u m - al though the whole system need not be i n equ i l ib r ium and might experience large fluxes of energy transport - and consequently the local a tomic 139 Appendix D. Glossary 140 popula t ion numbers depend only on the local electron temperature and density v i a the Saha-Bol tzmann equations. L T E is enforced by coll isional processes. L F luminosi ty function M S ma in sequence N L T E n o n - L T E (absence of L T E ) In N L T E the local level populations depend not only on the local electron tem-perature and density but on other parameters as well , e.g., the energy density of the radiat ion field. Models are usually computed by replacing the equ i l ib r ium Saha-Bol tzmann equations w i t h time-dependent rate equations. N L T E effects are larger at higher temperatures and lower densities when radiative processes might be dominant over collisions. P D E par t ia l differential equation P P A piecewise parabolic advection R H D radia t ion hydrodynamics R T radiative transfer R T E radiative transfer equation SDSS Sloan D i g i t a l Sky Survey U V ultraviolet R e f e r e n c e s [1] A b e l T . , Ann inos P. . Zhang Y . , Norman M . L . , 1997, New Astronomy, 2, 181 [2] A b e l T . , Haehnelt M . G . , 1999, astro-ph/9903102 [3] A b e l T . , N o r m a n M . L . , M a d a u P., 1998, astro-ph/9812151 [4] A b e l T . , Ann inos P., Norman M . L . , Zhang Y . , 1998, A p J , 508, 518 (astro-ph/9705131) [5] A l m e M . L . , W i l s o n J .R . , 1974, A p J , 194, 147 [6] Anderson J . L . , Spiegel E . A . , 1972, A p J , 171, 127 [7] Anninos P . , Zhang Y . , A b e l T . , Norman M . L . , 1997, New Astronomy, 2, 209 [8] Arons J . , Winger t D W . , 1972, A p J , 177, 1 [9] A u e r L . H . , Miha las D . , 1970, M N R A S , 149, 65 [10] B o n d J .R . , Kofmar i L . , Pogosyan D . , Wadsley J . , 1998, astro-ph/9810093 [11] B o n d J .R . , Myers S.T., 1996, A p J S , 103, 1 [12] Boyle B . J . , Shanks T . , Peterson B . A . , 1988, M N R A S , 235, 935 [13] B r y a n G . L . , Machacek M . , Anninos P. , Norman M . L . , 1999, A p J , 517, 13 [14] Buchler J .R . , 1979, J Q S R T , 22, 293 [15] Canto J . , Steffen W . , Shapiro P .R . , 1998, A p J , 502, 695 [16] Castor J.I., 1972, A p J , 178, 779 [17] Cen , R . , G n e d i n N . Y . , Ostriker J .P . , 1993, A p J , 417, 387 141 References 142 [18] Chandrasekhar S., 1950, Radiat ive Transfer. Oxford Univers i ty Press, Oxford [19] C h e n H . - W . , Lanzet ta K . M . , Pascarelle S., 1999, astro-ph/9904161 [20] C h i b a M . , N a t h B . B . , 1994, A p J , 436, 618 [21] C i a r d i B . , Ferrara A . , A b e l T . , 1998, astro-ph/9811137 [22] C i a r d i B . , Ferrara A . , Governato F . , Jenkins A . , 1999, astro-ph/9907189 [23] Couchman H . M . P . , Rees M . J . , 1986, M N R A S , 221, 53 [24] Cowie L . L . , Songaila A . , K i m , T . -S . , H u E . M . , 1995, A J , 109, 1522 [25] Efs ta thiou G., ' 1992, M N R A S , 256, 43P [26] E l v i s M . , Wi lkes B . J . , M c D o w e l l J . C . , Green R . F . , Bechtold J . , W i l l n e r S.P., Oey M . S . , Po lomski E . , C u t r i R . , 1994, A p J S , 95, 1 [27] Evans F . , 1998, J . Atmospheric Sciences, 55, 429 [28] Fan X . et al. (SDSS Col laborat ion) , 1999, astro-ph/9903237 [29] F i e l d G . B . , 1959, A p J , 129, 536 [30] G a l l i D . , P a l l a F . , 1998, A & A , 335, 403 (astro-ph/9803315) [31] G i r o u x M . L . , Shapiro P .R . , 1996, A p J S , 102, 191 [32] G n e d i n N . Y . , Ostriker J . R , 1997, A p J , 486, 581 [33] Gned in N . Y . , 1998, in Proceedings of the 19th Texas Sympos ium on Rela t iv is t ic As t ro-physics and Cosmology, held i n Paris , France, Dec. 14-18, 1998. J . P a u l , T . Montmerle , a n d . E . A u b o u r g , eds. [34] Godunov S . K . , 1959, Matematichesky Sbornik, 47, 271 [35] G o u l d A . , Weinberg D . H . , 1996, A p J , 468, 462 References [36] Griffiths L . M . , Barbosa D . , L idd le A . R . , 1998, astro-ph/9812125 [37] G u n n J . E . , Peterson B . A . , 1965, A p J , 142, 1633 [38] Haardt F . , M a d a u P. , 1996, A p J , 461, 20 [39] Haehnelt M . G . , Rees M . J . , 1993, M N R A S , 263, 168 [40] Haehnelt M . G . , Steinmetz M . , 1998, M N R A S , 298, L21 [41] Haehnelt M . G . , Natara jan P., Rees M . J . , 1998, M N R A S , 300, 817 [42] H a i m a n Z. , A b e l T . , Rees M . J . , 1999, astro-ph/9903336 [43] H a i m a n Z. , Loeb A . , 1999, astro-ph/9904340 [44] H a i m a n Z. , K n o x L . , 1999, astro-ph/9902311 [45] H a i m a n Z. , Loeb A . , 1997, A p J , 483, 21 [46] H a i m a n Z. , Rees M . J . , Loeb A . , 1997, A p J , 476, 458 [47] H a i m a n Z. , Loeb A . , 1998, A p J , 503, 505 (astro-ph/9710208) [48] H a i m a n Z . , T h o u l A . A . , Loeb A . , 1996, A p J , 464, 523 [49] H a i m a n Z . , Rees M . J . , Loeb A . , 1996, A p J , 467, 522 [50] Hernquist L . , K a t z N . , Weinberg D . H . , Mira lda-Escude J . , 1996, A p J , 457, L51 [51] Hogan C . J . , Rees M . J . , 1979, M N R A S , 188, 791 [52] H u E . M . , M c M a h o n R . G . , Cowie L . L . , 1999, astro-ph/9907079 [53] H u i L . , G n e d i n N . Y . , 1997, M N R A S , 292, 27 [54] H u m m e r D . G . , 1988, A p J , 327, 477 [55] Hummer D . G . , 1994, M N R A S , 268, 109 References 144 [56] Hummer D . G . , Storey P . J . , 1998, M N R A S , 297, 1073 [57] K a t z N . , Weinberg D . H . , Hernquist L . , Mira lda-Escude J . , 1996, A p J , 457, L57 [58] K a t z N . , Weinberg D . H . , Hernquist L . , 1996, A p J S , 105, 19 [59] Kepner J . V . , B a b u l A . , Spergel D . N . , 1997, A p J , 487, 61 [60] Kofman L . , Pogosyan D . , Shandarin S., 1990, M N R A S , 242, 200 [61] K o f m a n L . , Pogosyan D . , Shandarin S.F. , Melot t A . L . , 1992, A p J , 393, 437 [62] K r o o k M . , 1955, A p J , 122, 488 [63] Kunasz P., Aue r L . H . , 1988, J Q S R T , 39, 67 [64] Langseth J .O . , LeVeque R . J . , 1997, A wave propagation method for three-dimensional hyperbol ic conservation laws, preprint [65] Lepp S., Shul l J . M . , 1983, A p J , 270, 578 [66] Lepp S., Shul l J . M . , 1984, A p J , 280, 465 [67] LeVeque R . J . , 1997, J . Comput . Phys. , 131, 327 [68] LeVeque R . J . , 1998, Nonlinear Conservation Laws and F in i t e Volume Methods for A s -trophysical F l u i d F low, i n Computa t iona l Methods for Ast rophys ica l F l u i d F low, 27th Saas-Fee Advanced Course Lecture Notes, O . Steiner and A . Gautschy, eds., Springer-Verlag [69] Levermore C D . , Pomran ing G . C . , 1981, A p J , 248, 321 [70] Loeb A . , 1996, A p J , 459, 5 [71] Loeb A . , R y b i c k i G . B . , 1999, A p J , i n press, astro-ph/9902180 [72] L u L . , Sargent W . L . W . , Bar low T . A . , Church i l l C . W . , Vogt S.S., 1996, A p J S , 107, 475 References 145 [73] L u L . , Sargent W . L . W . , Bar low T . A . , Rauch M . , 1998, astro-ph/9802189 [74] M a d a u R , M e i k s i n A . , Rees M . J . , 1997, A p J , 475, 429 [75] M a d a u R , 1998, astro-ph/9807200 [76] M a d a u R , Haardt F . , Rees M . J . , 1999, A p J , 514, 648 (astro-ph/9809058) [77] M e i k s i n A . , 1994, A p J , 431, 109 [78] M e i k s i n A . , 1999, astro-ph/9902384 [79] M e i k s i n A . , W h i t e M. , Peacock J . A . , 1999, M N R A S , 304, 851 (astro-ph/9812214) [80] Miha las D . , 1980, A p J , 237, 574 [81] Miha las D . , Miha las B . , 1984, Foundations of Rad ia t ion Hydrodynamics . Oxford Univers i ty Press, New York [82] Mira lda-Escude J . , Ostriker J .P. , 1990, A p J , 350, 1 [83] Mira lda-Escude J . , Ostriker J .P. , 1992, A p J , 392, 15 [84] Mira lda-Escude J . , Rees M . J . , 1994, M N R A S , 266, 343 [85] Mira lda-Escude J . , Haehnelt M . , Rees M . J . , 1998, astro-ph/9812306 [86] Mira lda-Escude J . , Rees M . J . , 1998, A p J , 497, 21 [87] M i y a j i T . , Hasinger G . , Schmidt M . , 1998, astro-ph/9809398 [88] Navarro J .F . , Steinmetz M . , 1997, A p J , 478, 13 [89] Nakamoto T . , U m e m u r a M . , 1998, A New Numer ica l M e t h o d to Solve the M u l t i -Dimensional Radia t ive Transfer Equa t ion , poster at the A n n . Meet ing of the As t ron . Soc. of Japan [90] N o r m a n M . L . , Paschos P. , A b e l T . , 1998, astro-ph/9807282 References 146 [91] O h S . R , 1999, astro-ph/9904255 [92] Oort J . H . , 1984, A & A , 139, 211 [93] Pa rk Y . - S . , Hong S.S., 1998, A p J , 494, 605 [94] Peebles P . J . E . , 1993, Pr inciples of physical cosmology Pr ince ton Univers i ty Press, Pr ince-ton, New Jersey [95] Pe i Y . G . , 1995, A p J , 438, 623 [96] Press W . H . , Teukolsky S .A. , Vet ter l ing W . T . , F lannery B . P . , 1992, Numer ica l Recipes i n For t ran. 2nd ed., Cambridge Universi ty Press, Cambridge [97] Q u i n n T . , K a t z N . , Efstathiou G . , 1996, M N R A S , 278, L49 [98] Razoumov A . , Scott D . , 1999, M N R A S , accepted, astro-ph/9810425 [99] R ico t t i M . , Gned in N . Y . , Shul l J . M . , 1999, astro-ph/9906413 [100] Sathyaprakash B . S . , Sahni V . , M u n s h i D . , Pogosyan D . , Melo t t A . L . , 1995, M N R A S , 275, 463 [101] Seager S., Sasselov D . D . , Scott D . , 1999, A p J , i n press [102] Schaye J . , Theuns T . , Leonard A . , Efstathiou G . , 1999, astro-ph/9906271 [103] Schmidt M . , Schneider D . P . , G u n n J . E . , 1995, A J , 110, 68 [104] Schuster A . , 1905, A p J , 21, 1 [105] Scott D . , Rees M . J . , 1990, M N R A S , 247, 510 [106] Shapiro P .R . , 1986, P A S P , 98, 1014 [107] Shapiro P .R . , G i r o u x M . L . , 1987, A p J , 321, L107 [108] Shapiro P .R . , G i r o u x M . L . , B a b u l A . , 1994, A p J , 427, 25 References '. [109] Shapiro P .R . , Raga A . C . , Mellerna G . , 1998, astro-ph/9804117 [110] Shaver P . A . , Windhors t R . A . , M a d a u P., de B r u y n A . G . , 1999, astro-ph/9901320 [111] Sidilkover D . , 1989, P h . D . thesis, the Weizmann Institute of Science, Israel [112] Smal l T . A . , Blandford R . D . , 1992, M N R A S , 259, 725 [113] Spitzer L . , Jr . , 1968, Diffuse matter i n space. Interscience Publishers, New Y o r k [114] S tanci l P . C . , Lepp S., Dalgarno A . , 1996, A p J , 458, 401 [115] Stone J . M . , Miha las D . , 1992, J . Comput . Phys. , 100, 402 [116] Stone J . M . , Miha las D . , Norman M . L . , 1992, A p J S 80, 819 [117] Storey P . J . , Hummer D . G . , 1991, Comput . Phys. Commun. , 66, 129 [118] Struchtrup H . , 1997, Anna l s of Physics , 257, 111 [119] Struchtrup H . , 1998, Annals of Physics, 266, 1 [120] Sunyaev R . A . , Zeldovich Y a . B . , 1972, A & A , 189 [121] Sunyaev R . A . , Zeldovich Y a . B . , 1975, M N R A S , 171, 375 [122] Ta j i r i Y . , U m e m u r a M . , 1998, A p J , 502, 59 (astro-ph/9806046) [123] Tegmark M . , Si lk J . , Blanchard A . , 1994, A p J , 420, 484 [124] Tegmark M . , S i lk J . , Rees M . J . , B lanchard A . , A b e l T . , P a l l a F . , 1997, A p J , 474, 1 [125] Tozz i P. , M a d a u P., M e i k s i n A . , Rees M . J . , 1999, astro-ph/9905199 [126] Thorne K . S . , 1981, M N R A S , 194, 439 [127] Tucker W . H . , 1975, Rad ia t ion Processes i n Astrophysics. T h e M I T Press, Cambridge References 148 [128] U m e m u r a M . , Nakamoto T . , Susa H . , 1998, i n M i y a m a S. M . , Shibata K . , eds, Numer ica l Astrophysics 1998. Kluwer , in press [129] Unno W . , Spiegel E . A . , 1966, P A S J , 18, 85 [130] Valageas P. , Si lk J . , 1999, astro-ph/9903411 [131] Warren S.J . , Hewett P . C , Osmer P.S. , 1994, A p J , 421, 412 [132] Weinberg D . H . , Hernquist L . , K a t z N . , 1997, A p J , 477, 8 [133] Weymann R . J . , Stern D . , Bunker A . , Spinrad H . , Chaffee F . H . , Thompson R . I . , Storrie-L o m b a r d i L . J . , 1998, A p J , 505, L95 [134] Wouthuysen S .A. , 1952, A J , 57, 31 [135] Zhang Y . , Anninos P., Norman M . L . , M e i k s i n A . , 1997, A p J , 485, 496 [136] Zhang Y . , M e i k s i n A . , Anninos P. , N o r m a n M . L . , 1998, A p J , 495, 63 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085451/manifest

Comment

Related Items