Late Reheating of the IGM by Quasars: A Radiation Hydrodynamical Simulation of Helium II Reionization
Pascal Paschos, Michael L. Norman, James O. Bordner, Robert Harkness
aa r X i v : . [ a s t r o - ph ] N ov Late Reheating of the IGM by Quasars: A Radiation HydrodynamicalSimulation of Helium II Reionization
Pascal Paschos , Michael L. Norman , James O. Bordner & Robert Harkness ABSTRACT
We study the ionization and thermal evolution of the intergalactic medium duringthe epoch of He II reionization by means of radiation hydrodynamical cosmological sim-ulations. We post-process baryonic density fields from a standard optically-thin IGMsimulation with a homogeneous galaxy-dominated UV background (UVB) which reion-izes H I and He I at z=6.5 but does not have any contribution to the ionization of He II .Therefore, we suppress the He II photoheating contribution to the gas temperature dueto the homogeneous UVB. Quasars are introduced as point sources throughout the 100Mpc simulation volume located at cold dark matter (CDM) density peaks consistentwith the Pei luminosity function. We assume an intrinsic quasar spectrum J ( ν ) ∝ ν − . and a luminosity proportional to the halo mass. We evolve the spatial distribution ofthe He II ionizing radiation field at h ν = 4, 8, and 16 Ryd using a time-implicit variabletensor Eddington factor radiative transfer scheme. Simultaneously, we also solve for thelocal ionization of He II to He II and the associated photoheating of the gas includingopacity effects. We find that the percolation of the He III regions is essentially completeby z=2.5. When comparing to a self-consistent optically thin simulation at the sameredshift, in which He II is also ionized by the uniform UVB, we find that inclusion ofopacity effects results in higher IGM temperature by a factor of approximately 1.7 atthe mean gas density level. We construct synthetic absorption line spectra from whichwe derive statistical parameters of the He II Ly α forest . We use 300 long (∆ z = 0 . z = 2 . ± . II Ly α line transmissionof ¯ F = 0 . ± . ≃
11% the meanvalue. The opacity effect on the gas temperature is shown by comparing the broadeningwidth of the H I and He II Ly α lines to the results from the self-consistent optically thinsimulation. We find a shift by approximately 1.25 km/s to higher b-parameter valuesfor both H I and He II . Finally, we estimate the relative broadening width between the Laboratory for Computational Astrophysics,Center for Astrophysics and Space Sciences,University of California at San Diego, La Jolla, CA 92093, U.S.A.Email: ppaschos, mlnorman, [email protected] Physics Department, UCSD San Diego SuperComputer Center II median b-parameter is about 0.8 times the medianH I broadening width. This implies that the He II absorbers are physically extendedconsistent with conclusions from observed lines of sight. Subject headings: cosmology: hydrodynamical simulations: radiation transfer
1. Introduction
The thermal evolution of the IGM after reionization is governed chiefly by photoheating (Efs-tathiou 1992; Miralda-Escud´e & Rees 1994). Models for the thermal evolution of the IGM (Miralda-Escud´e & Ostriker 1990; Giroux & Shapiro 1996; Abel & Haehnelt 1999) generally assume thatH I and He I are nearly fully ionized by z ∼ II ionizes somewhat later at z ∼ II ’s higher ionization po-tential and recombination rate (Sokasian, Abel & Hernquist 2002; hereafter SAH). Observationalsupport for late He II reionization is summarized in SAH. Abel & Haehnelt (1999) emphasized theimportance of opacity effects during reionization in establishing the post-reionization temperatureof the gas. They showed that models using the optically thin expression for the photoheating rateduring He II reionization underestimate the IGM temperature at mean density by a factor of ∼ α forest (Cen et al. 1994; Zhang, Anninos& Norman 1995; Hernquist et al. 1996; Miralda-Escud´e et al. 1996; Zhang et al. 1997, 1998;Theuns et al. 1998) generally adopt the optically thin expression for photoheating for simplicityand computational economy. This is a reasonable assumption for the H I and He I photoheatingat z ∼ II photoheating as He II is in theprocess of being ionized by the percolation of He III spheres (SAH). These simulations thereforeunderestimate the temperature of the IGM during the epoch of helium reionization. The stan-dard approach taken in these simulations is to assume a homogeneous photoionizing backgroundwhich evolves with redshift consistent with observed quasar and galaxy counts, such as that byHaardt & Madau (2001). The ionization and thermal state of the baryonic gas is then computedself-consistently with its dynamics by solving the equations of hydrodynamic cosmology includingradiative heating and cooling (Cen 1992; Katz, Weinberg & Hernquist 1996; Anninos et al. 1997).It is found that the temperature of the low density IGM is determined almost entirely by photo-heating balancing adiabatic cooling due to cosmic expansion. This results in a tight relationshipbetween gas temperature and density, the so-called equation of state of the IGM (Hui & Gnedin1997): T = T ∆ β (1)where ∆ = ρ/ ¯ ρ is the gas overdensity and β is redshift dependent but in the range 0 ≤ β ( z ) ≤ . II reionization raises T and reduces β relative to an optically thin calculation. A technique often employed to include opacity effectswithin hydrodynamic simulations is to simply multiply the He II photoheating rate by a constant 3 –factor X He II ≥
1, with a value of 2-4 being sufficient to match observed temperatures (e.g., Bryan& Machacek 2000; Jena et al. 2005).The combination of high spectral resolution quasar absorption line observations and hydro-dynamical cosmological simulations provide a means for measuring the thermal evolution of theIGM. The thermal state of the gas is generally deduced from H I Ly α linewidths (b-values) in highresolution spectra (Rauch et al. 1997; Schaye et al. 1999, 2000; Bryan & Machacek 2000; Theunset al. 2000; Bolton et al. 2005; Jena et al. 2005) although the flux power spectrum has also beenemployed (Zaldarriaga, Hui & Tegmark 2001). If the temperature of the IGM alone determinedthe b-values, then it would be straightforward to measure the temperature from high resolutionspectra. However, Hubble broadening is always of the same order as the thermal broadening inLy α forest absorption lines (Schaye 1999), hence the need for comparison with simulations. The-uns, Schaye & Haehnelt (2000) used the b-parameter distribution to measure the temperature ofthe IGM at z=3.25, finding T ≥ ,
000 K. Schaye et al. (2000) used the lower cutoff in thelinewidth–column density scatter diagram to measure the temperature evolution of the IGM overthe redshift range 2-4.5. They found evidence of late reheating at z ∼
3, which they ascribed to lateHe II reionization by quasars. Bryan & Machacek (2000) independently explored the same diagnos-tics and found that temperature estimates were sensitive to the assumed cosmology, in particularthe amplitude of mass fluctuations on a few Mpc scales. Zaldarriaga, Hui & Tegmark (2001) usedthe falloff of the Ly α forest flux power spectrum at small scales to measure the IGM temperature,and found T ∼ × K, dependent on their assumed β .A proper calculation of He II reionization would treat quasars as point sources within thesimulation volume and solve the equation of radiative transfer for the spatial distribution of theionizing background coupled self-consistently to the dynamical, thermal, and ionization evolutionof the gas. This was done in an approximate way by SAH, who post-processed a series of densityfields taken from a SPH hydrodynamical simulation of the IGM using the GADGET code. Thehydrodynamical simulation used the optically thin prescription for ionizing and heating the gas,assuming all ionization states of hydrogen and helium were in ionization equilibrium with the UVBof Haardt & Madau (1996). In the post-processing step, the helium reionization calculation wasrecomputed for an evolving quasar source population treated as point sources within the volume.Quasars were assumed to have a constant lifetime of 10 yr. Every 10 yr, peaks within thedensity field of a suitably chosen data dump from the hydrodynamic simulation were populatedwith quasars consistent with an empirical luminosity function. Around each point source the staticequation of radiative transfer for He II ionizing photons was solved using a photon-conserving ray-casting scheme. The ionization state of He II was then updated ignoring thermal feedback to thegas. SAH found that for reasonable parameter choices, He II reionization occurred in the range3 ≤ z ≤ III spheres percolate and eventually merge. Our hydrodynamic simulation wasperformed on an Eulerian grid of 512 cells versus SAH’s 224 in the same volume, giving us ∼ ∼ . I and He I due tostellar sources only using the homogeneous UVB of Haardt & Madau (2001) arising from galaxies”GAL”. In the post-processing step, quasars are treated as point sources, which ionize He II toHe III . Although our radiative transfer scheme is based on completely different spatial and angulardiscretizations as SAH, the important difference for the purpose of this paper is that we solve theRT equation at three frequencies hν =4, 8 and 16 Ryd in order to evaluate, albeit crudely, the localHe II photoheating rate taking the processed QSO spectrum into account.The outline of this paper is as follows. In § I and He I reionization described in §
3. In § II reionization and also describe our treatmentof quasars in the simulation volume. In § II reionization. In § § II , and in § I and He I , followed by the inhomogeneous reionization of He II , our calculation is not fullyself-consistent. Although we keep track of the late reheating, we do not modify the underlyingdensity fields nor do we alter temperature-dependent recombination rates, which will affect thedetailed ionization state of the gas. In § II reionization, as well as the importance ofneglecting these coupling effects. We show that these effects, while present, are small, and do notseriously undermine our estimate of the reheating. In § II reionization based on synthetic H I and He II Ly α absorption line spectra derived from oursimulation. In § II photoheating on the ionization state of the gas (Appendix B), anddocument the reaction rate coefficients we use (Appendix C).
2. Simulations2.1. Cosmic Realization
In this work, we present the results from post-processing the redshift evolution of a cosmicrealization computed with the Eulerian cosmological hydrodynamic code Enzo (Bryan & Norman1997; Norman & Bryan 1999; O’Shea et al. 2004). The box of size 67 h − Mpc comoving wasevolved in a flat (Ω = 1) ΛCDM cosmology on a unigrid mesh of 512 grid cells and 512 dark 5 –matter particles from z=99 to z ≃
2. We used the initial power spectrum of matter fluctuationsby Eisenstein & Hu (1999) to initialize the calculation which was then computed forward underthe Zel’dovich approximation (Bertschinger & Gelb 1991). Our choice of cosmological parametersis σ = 0 .
8, Ω Λ = 0 .
7, Ω b = 0 . , n s = 1 with a present day Hubble constant of h = 0 .
67 inunits of 100 km/s/Mpc. With these parameters our cosmic realization has a mesh resolution of130 h − ≃
200 kpc within the 100 Mpc cube, and ≃ . h − × M ⊙ dark matter particle mass. InFigure (1) we show a volumetric rendering of the cosmic gas distribution as mapped by the baryonoverdensity in our simulation at z=2.6.One of the fundamental simplifications used here, which places limits on the validity of ourresults is that the ionization of H I and He I are treated entirely differently than He II . The moti-vation is that of numerical simplicity and calculation speed. From the IGM cosmological evolutionstandpoint, hydrogen and neutral helium are believed to be globally ionized at an earlier cosmicepoch than He II . In the simulation discussed here, the photo–ionization/heating rates of H I andHe I are computed self-consistently during the hydrodynamical calculation in the optically thinregime. The premise of our argument is that, if one adopts a picture where neutral hydrogen andneutral helium reionization is completed by z ≃
6, such species will have large ionization fractionsby the time He II reionization occurs, which theoretical models and observations place at z ≃ − I /He I species is achieved by the evolving uniform metagalactic fluxdue to stellar sources as computed in Haardt & Madau (2001). The uniform ultra-violet backgroundphoto-ionizes and photo-heats the IGM, however it is prevented by hand to radiatively alter theHe II abundance. We do so by suppressing the ionization/heating rates in Enzo that control theHe II ↔ He III chemistry. The latter is computed separately by our inhomogeneous point-sourcedistribution of QSOs which is discussed in § II reionization proceeds, under the limitations discussed above, via the mergersof individual He III
I-fronts during the cosmic epoch that spans the redshift interval z = 6 − II reionization In this section we introduce some definitions and basic results related to He II reionization thatwe will refer to later on. The softness parameter of the cosmic radiation is defined as S = Γ HI Γ HeII . Inthe optical thin limit under a single power law profile for the radiation spectrum, J ν = J ( νν ) − α q ,where J denotes the volume averaged mean intensity at the hydrogen ionization threshold, wecan estimate the photoionization rates to be Γ HI ≡ Γ = πJ σ oHI h α q and Γ HeII ≡ Γ = πJ σ oHeII h α q − α q . Dividing the last two relations yields an estimate for the softness parameterin the optical thin limit equal to S = 4 α q . Therefore, the ability of a background radiation 6 –field to ionize the H I and He II species depends on the spectral slope. A soft radiation field wouldhave large spectral slope and would primarily favor hydrogen ionization. He II ionization thereforerequires a hard ionizing spectrum.The available sources of radiation at redshifts z ≤
10 are of two types. Stellar sources, as-sociated with radiation from galaxies forming in the cosmic medium at such redshifts, have largesoftness parameter values of S stellar ∼ × . QSOs on the other hand, have much smaller softnessparameter values of S QSO ∼
50 which makes them ideal for He II ionization. However, the numberdensity of QSOs, as measured by observations sharply rises in the interval 6 > z > z ≈ z em = 6 .
56 (Hu,Cowie & McMahon 2002). Direct observations of QSO in Ly α spectra are very difficult becauseof the hydrogen Gunn-Peterson effect, the optical depth manifestation of hydrogen neutrality athigh redshifts (Fan, Carilli & Keating 2006). This suggests that either the quasars have a wellestablished population by z ≈ II ionizer, He II reionization is expected to occurlater; at times when a sufficient number hard photons becomes available.The Gunn-Peterson optical depth of He II for a uniform IGM can be expressed in the sameway as the corresponding H I and that leads to the ratio τ HeII ( z ) τ HI ( z ) = n HeII n HI σ HeII σ HI = n HeII n HI . Theratio of optical depths is known as the R-factor and is used in observations to measure the relativeproperties between the He II and H I Ly α forest spectra. Observations find that in transmissionspectra longward of quasars seen at emission redshifts z ∼ R ≈ −
100 within δz ∼ . z & .
5. In theoptical thin limit the R-factor and the softness parameter relation can be derived.In ionization equilibrium the ratio of n HeII /n HI can be computed if we assume that the higherionization states of both species (He III , H II ) dominate. In that case, for hydrogen we get χ Γ = n e α (1 − χ ) where χ = n HI /n H and α is H II the recombination coefficient. Similarly for He II wecan also write ψ Γ = n e α (1 − ψ ), where ψ = n HeII /n He and α is the He II recombinationcoefficient. The above two relations yield χ /ψ ≃ Γ Γ α α , for χ ≪ ψ ≪
1. The lastequation allows for the determination of the ratio n HI /n HeII = χ n H ψ n He ≃ × Γ α Γ α . Therefore, for α α ∼ T ≈ K) the ratio of number densities becomes n HeII n HI ≃ S . From the last relation itfollows that R ≃ . S . For S ≃ × (stellar radiation) the relation predicts R ≃ S ≃ R ≃
5. In addition, we can compute that since the quantity n e α ( T ) / Γ = τ ion /τ rec , the ratio of the ionization to recombination time scales, τ ionHeII τ recHeII = 5 S τ ionHI τ recHI Typical softness parameter values S ≥
100 then yield τ ionHeII τ recHeII ≥ × τ ionHI τ recHI . The last relation 7 –shows that, due to the larger recombination coefficient and the fewer number of photons available forionization, it is more difficult to ionize He II compared to H I . The shorter (longer) recombination(ionization) time scale effectively restricts He II reionization to take place at a slower rate.Large values in the R-factor is suggestive of large optical depths in He II and/or small opticaldepths in neutral hydrogen. Comparable optical depths at late redshifts, which correlate to smallcosmic neutral hydrogen fraction and therefore small cosmic fraction in He II , yield small valuesfor the R-factor. The observations can therefore infer a range in the observed softness parameterfrom measuring the ratio of optical depths in the He II and H I line forest. The observed range ofR=10-100 consequently suggests that S=100-1000 between z=2-3 although the upper limit is anoverestimate because we assumed ψ ≪
1. Because such observed values are sampled in trans-mission spectra that probe the ionization phase transition of He II to He III , we can infer that atthat epoch the galaxy dominated ultraviolet background is gradually being replaced by a quasardominated type.In §
6, we show that the R-factor evolves with redshift from large to small values as it is com-puted in transmission through the computational volume. The quantity that is actually computedis the η parameter defined as the column density ratio between He II and H I , η = N HeII N HI . TheR-factor evolution and value follows directly as R = η . The gradual evolution of the R-factor and η parameter from large values is indicative of the cosmic evolution of the softness parameter froma stellar to a quasar dominated type. For the quasar placement and evolution in the simulation we use the quasar luminosity functionby Pei et al. (1995), shown in Equation (2). This luminosity function was also used by Haardt& Madau (1996,2001) to derive the quasar emissivity and the volume average photoionization andheating rates used in the simulation for every species other than He II . φ ( L, z ) = φ ∗ /L ∗ ( z )[ L/L ∗ ( z )] β + [ L/L ∗ ( z )] β (2) L ∗ ( z ) = L ∗ (0)(1 + z ) α q − e ζz (1 + e ξz ∗ ) e ξz + e ξz ∗ Our first requirement is the placement of a single QSO in the computational volume at z=6.5.The redshift is for all practical purposes a matter of choice, since there is no accurate prediction ofwhen the first QSO appears in the universe. We adopt a scenario where quasars become visible inobservations after the epoch of hydrogen reionization is completed by z ≃ . ≃ × M ⊙ . Therefore, quasar point sources have to be placed in the centers of halos above thiscutoff. We do so by identifying the evolution of a list of halo centers through the computationalvolume from z = 6 . z = 2. The number of quasars per redshift interval is determined by theluminosity function in Equation (2) and the functional form of the mass-halo to quasar luminosityrelation.We assume an intrinsic quasar spectrum J ν = J ( νν ) − α q for the LyC part of the spec-trum with a spectral index α q = 1 .
8. Therefore, the number flux of emitted LyC photons, in photons/s/cm , is then ˙ n ph = πJ h R ∞ ǫ − αq ǫ dǫ . In this relation, ǫ = hνhν . The integrationyields ˙ n ph = πJ hα q . Similarly the LyC energy flux is l LyC = πν J α q − , which yields l LyC =( hν ) ˙ n ph α q − α q ⇒ L = ( hν ) ˙ N ph α q − α q . In the last equation, L and ˙ N ph represent the totalLyC luminosity (ergs/s) and emitted photon rate (photons/s) respectively per point source (QSO)above the H I ionization threshold of 1 Ryd (13.6 eV). This allows for a parametrization of theemitted ionizing flux based on the number of LyC photons rather than energy. The motivationis entirely for consistency with the photon-conserving schemes in simulating reionization (Abel &Haehnelt 1999; Sokasian, Abel & Haehnelt 2001, 2002; Ciardi et al. 2003; Whalen, Abel & Norman2004). The total luminosity for the He II ionizing radiation above 4 Ryd (54.4 eV) is then obtainedthrough L = 4 − ( α q +1) × L . The luminosity of each source is determined by the dark matter massof the halo that initially creates it. For a dark matter halo of mass M halo we put in a UV sourceemitting ˙ N ph = 10 × M halo H I ionizing photons/s. This is equivalent to placing a ˙ N ph = 10 ph/s mini-quasar inside a 10 M ⊙ dark matter halo which is the prescription used in Abel &Haenhelt (1999). Sokasian et al. (2001) investigated an array of QSO placement methods in thecomputational volume and found that the results are largely insensitive to the choice. We adoptthe linear relation for convenience. In this work, the most massive halo computed at at z=2.5 hasmass M halo ≈ × M ⊙ . If a QSO source is placed there then it will emit LyC photons at a rateof ˙ N ph = 6 × s − . For an input spectrum with slope α q = 1 . L = 2 . × ergs/s and L = 6 . × ergs/s for H I and He II respectively.The placement of the point sources in the volume is dynamical in nature. The list of darkmatter halos is assigned a quasar source with a luminosity value determined by our phenomeno-logical prescription. The location of the quasar from that point on is locked to the position ofthe dark matter particle closest to the center of the halo. The distribution of luminosities is thenintegrated from the higher value to the smallest up to the point where the average luminosity perunit volume reproduces the distribution fit given by Equation (2) at each redshift. For simplicity,we do not evolve the luminosity in each dark matter halo, which remains the same at initialization.As the luminosity function increases with decreasing redshift additional sources are spawned in thesimulation, leading to an overall increase of their numbers. At z . § f ij ( x ) at t n is an interpolatedvalue between the local values of the two data dumps whose redshifts bound the instantaneoustime ( z ≤ t n ≤ z ). As we will discuss in the next sections, the cosmic evolution of the heliumreionization is decoupled from the corresponding hydrogen one. Consequently, the placement of aHe II ionizing source on a density peak with a precomputed hydrogen and therefore electron densitycan yield significant discrepancy between our calculation and a self-consistent one, particularly inclose proximity to the source.
3. Homogeneous Hydrogen Ionization
As mentioned above, in the numerical and physical setup of this present work, we treat theionization of H I and He I separately from He II . The metagalactic flux we use from Haardt &Madau (2001) tabulates the contributions due to galaxies (GAL) and quasars (QSO) separately.In Enzo we have the option of running a simulation with GAL only, QSO only, or GAL+QSO. In astandard optically thin simulation the GAL+QSO UVB would be used. We have carried out sucha simulation, hereafter called Simulation A, for comparison with the inhomogeneous reionizationsimulation. For the latter, we first run a hydro simulation using the GAL UVB to ionize H I andHe I . Then He II reionization is accomplished by treating quasars as point sources as detailedin §
4. Since our quasar population follows the same luminosity function and intrinsic spectrumas assumed by Haardt & Madau, we are able to compare the homogeneous and inhomogeneoussimulations directly. We consider the effect of the ǫ ≥ . II ⇔ He III chemistry while the H I andHe I ionization states remain unaffected. The treatment is by all measures an approximation thatallows the problem to be solved as a single species ionization problem only and therefore it remainsconceptually simple. In reality, the ionization of He II has two inter-dependent effects; it releasesadditionally one electron per He II ionization which in turn, when thermalized, raises the meantemperature of the IGM. These two effects combined would in principle shift the ionization balanceof the H I /H II and He I /He II species. However, we show in this section that because by the timethis takes place hydrogen and helium have already large ionization fractions the aforementionedeffects are small.Starting with the rate equation for hydrogen and the cosmic mean density, we can safelyignore the collisional contributions and write the chemical balance in the proper frame of referenceas follows: ˙ n HI = − H ( z ) n HI − n HI Γ + n e n HII α ( T ) . (3)In Equation (3), n e = n HII + n HeI + 2 n HeII , Γ is the integrated H I photoionization rate,and α ( T ) is the radiative recombination coefficient in the H II + e → H I + γ reaction and is afunction of temperature. In a cosmic medium, n He = f n H , where f ≃ for a mass of fraction 10 –of ρ H = 0 . ρ and ρ He = 0 . ρ respectively. The He II number density can then be rewritten as n HeII = n He − n HeI − n HeIII which in turn allows the electron density to be expressed as follows: n e = n H [ χ HII + (1 − χ HeI + χ HeIII )12 ] (4)Changing slightly the notation, we can rewrite Equation (3) according to Appendix (A) in termsof the ionization fractions χ = n HII /n H , ψ = n HeI /n He , ψ = n HeII /n He and ψ = n HeIII /n He as: ˙ χ = Γ (1 − χ ) − χ n H α ( T )[ χ + 112 (1 − ψ + ψ )] ⇒ ˙ χ ≈ Γ (1 − χ ) − χ n H α ( T )[ χ + 1 + ψ
12 ] (5)In the last equation, we assumed that almost all of the helium is highly ionized to the He II state ψ →
0. In ionization equilibrium ( ˙ χ = 0), after setting A ≡ (1 + ψ ) /
12 we can rewrite Equa-tion (5) as (1 + χ A ) χ − χ = Γ / ( An H α ). On one hand, the ionization of the He II would increasequantity A due to the increase in the He III abundance. On the other hand, the increased tem-perature would decrease the recombination coefficient which can be approximated as α ∝ T − β ,for temperatures in the range of T ≃ − K and β = 0 .
51. The approximation is based onexpression fits by Bugress (1964).Therefore, we will investigate two extreme cases. Case (I) represents zero He
III abundance, ψ → T I . Case (II) represents a limitof almost full He II ionization, ψ →
1, that would correspond to a new temperature T II . In bothcases, the hydrogen ionization rate is the same because our treatment of the individual QSO sourcesplacement has a distribution that is statistically identical to the global average used for H I andHe I ionization. Therefore, we can form the ratio:(1 + χ II /A II ) χ II (1 − χ I )(1 + χ I /A I ) χ I (1 − χ II ) = ( T II T I ) β A I A II (6)Equation (6) is solved in Appendix (B), for χ II in the range of 1 − χ I = 10 − − − and for A II = 1 / A I = 1 / T II /T I = 1 . −
2. The ratio of temperatures is based on estimates of thetemperature increase due to the photoelectrons injected in the IGM from the ionized He II atoms(Haehnelt & Steinmetz 1998; Abel & Haehnelt 1999) and will be discussed in more detail in § II ionization results in a decrease of the neutral hydrogen fractionprimarily due to the decrease of the ionized hydrogen’s recombination coefficient.The ≈ −
25% reduction depends on the temperature increase, where the largest increase(a factor of 2) produces the biggest shift. The fractional decrease in the neutral hydrogen frac-tion is smallest at low fractions of ionized hydrogen. For a typical neutral hydrogen fraction of 11 – χ ≃ − at post-reionization redshifts, we can then estimate that our treatment (of not updatingthe hydrogen abundances) would overestimate the neutral hydrogen fraction by about 25% if thegas temperature increase is between 1.5-2 times the non-perturbed value. Consequently, we an-ticipate an overestimate by the same amount in the optical depth of hydrogen Ly α radiation andan underestimate in the mean transmitted flux. The problem is similar to the one described inZhang et al. (1997), Bryan et al. (1999) and Jena et al. (2005). An unmodified UVB by Haardt& Madau (1996), applied in those simulations, would not yield an agreement with the observedb-parameter of the Ly α forest. An adjustment by a factor of 1.5-2.0 in the He II photoheating ratewas then necessary to match the observed results. Therefore, we caution against a strict interpre-tation of the resulting hydrogen Ly α forest in the original calculation in which we suppressed theHe II photoionization and photoheating processes altogether.
4. Inhomogeneous Helium II Reionization
We compute the inhomogeneous He II reionization due to an evolving distribution of local QSO-type sources by solving for the time and spatial evolution of the ionizing radiation energy density E ν ( r , t ) as shown in following section. In § § § We consider simulation volumes of box length L much smaller than the horizon scale L ≪ L H = c/H ( z ). Also, prior to bubble overlap, the time between emission and absorption of a randomHe II ionizing photon will be much shorter than a Hubble time. In this limit, the cosmologicalradiative transfer equation reduces to the familiar one (Norman, Paschos & Abel 1998):1 c ∂I ν ∂t + ˆ n · ∇ I ν a = η ν − χ ν I ν (7)where I ν is the monochromatic specific intensity, η ν , χ ν are emission and extinction coefficients,and ν is the instantaneous, comoving frequency. In equation (7) the gradient is comoving, andhence we divide by the cosmological scale factor a to convert to proper distances. The zeroth andfirst angular moments of equation (7) yield ∂E ν ∂t + 1 a ∇ · F ν = ǫ ν − ck ν E ν (8)and 1 c ∂ F ν ∂t + 1 a ∇ · P ν = − c k ν F ν , (9)where E ν , F ν and P ν are the radiation energy, flux vector, and pressure tensor, respectively. Allquantities are measured in the fluid (proper) frame, where ǫ ν = 4 πη ν is the emissivity, and k ν
12 –is the absorption coefficient. The UV photons are scattered by the free electrons and for therange of electron number densities in the IGM (in the redshift interval we are interested in) wecan estimate that λ s (the scattering mean free path) = n e σ Th > L H . Therefore, we ignore thescattering coefficient.The radiation pressure tensor is coupled to the energy density P ν → E ν in the momentequations through the tensor Eddington factor, f ν = P ν E ν . The latter guarantees the correct directionof the flux vector (Mihalas & Mihalas, 1984). It can be shown that the time derivative term inequation (9) is small compared to the rest if we integrate using a timestep long compared to thelight crossing time of a computational cell, as we do. Therefore, dropping the time derivative inequation (9) and combining it with equation (8) we get ∂E ν ∂t = 1 a ∇ · [ ck ν ∇ · ( f ν E ν )] + ǫ ν − ck ν E ν (10)where ∇ · ( f ν E ν ) ≡ ∂∂x i ( f ijν E ν ) . In Equation (10), ǫ ν is the spatially discrete monochromatic emissivity at the locations of theemitting sources and k ν = n HeII σ HeIIν is the local opacity, where the functional form of σ HeIIν isgiven by Osterbrock (1989) and in Appendix C. Equation 10 can be solved for E ν ( r , t ) for a givensource distribution provided the spatially-dependent Eddington tensor is known. Formally, f ij isobtained from angular quadratures of the specific intensity (e.g., Hayes & Norman 2003). However,we wish to avoid solving the full angle– and frequency–dependent equation of radiative transfer.Instead we employ a geometric closure introduced by Gnedin & Abel (2001) in which we calculatethe radiation pressure tensor assuming the medium is optically thin. In this limit P ij ν ( r ) = 14 πc N X k L k ν | r − s k | (ˆ n i ( r − s k ))(ˆ n j ( r − s k )) | r − s k | (11)where s k are the positions of the ionizing sources, and ˆ n i,j are the direction vectors (basis) at thepoint r . Equation 11 describes pure radial streaming radiation from a collection of point sourcesin a transparent medium. Until He III bubbles begin to overlap, this is an excellent approximationinside the He
III regions but a poor one outside. However, since there is very little ionizing radiationin the He II regions, it makes little difference what one chooses for f . As discussed by Gnedin &Abel, the greatest error is when two bubbles begin to overlap and the two ionizing sources begin to“see one another.” If one source is much more luminous that the other, this can lead to a ∼ R V g E ν dV = V g E ν , R V g ǫ ν dV = V g ǫ ν , R V g k ν dV = 13 – V g k ν , where V g is the volume of a grid cell and E ν , ǫ ν , k ν are now understood to be cell averages.Equation (10) then becomes: ∂E ν ∂t + 1 aV g I cell F ν · d S cell surf = ǫ ν − ck ν E ν (12) F iν = − cak ν ∂∂x j ( f ijν E ν )Our time-implicit discretization scheme will be discussed in a follow up paper. Briefly, equation (12)is discretized on a uniform cartesian mesh and integrated using backward Euler time differencing.Spatial discretization of the RHS of equation (12) yields a 19-point stencil. The resulting sparse-banded system of linear equations is solved using the stabilized biconjugate gradient (BiCGstab)algorithm implemented in the MGMPI package developed by the Laboratory for ComputationalAstrophysics (Bordner 2002).For this problem, we compute the radiative energy density at three frequency values abovethe ionization threshold value, E , E , E , corresponding to photon energies ǫ = 4, ǫ = 2 × ǫ and ǫ = 4 × ǫ in Rydberg units ( ǫ = hνhν ) respectively. The three points plus a fourth oneat ǫ = 32 Ryd with E = E ( ǫ ǫ ) − α q , are used to infer the interpolated profile of a 4th degreepolynomial E ( ǫ ) /E = P k =4 k =0 c k ( ǫ ) − k between ǫ = 4 −
16 Ryd. At ǫ ≥
16 Ryd, we assume that E ( ǫ ) = E ( ǫǫ ) − α q , which follows from the reasonable assumption that at four times the ionizationthreshold energy, ≈ . keV already in the soft X-ray energy band, there is little effect in theattenuation of the radiative energy due to opacity. The interpolation is necessary in order to be ableto compute the ionization and heating rates which involve an integration in frequency space (photonenergy). The use of three frequencies and locally computing c k bypasses the requirement more manyfrequency bins or the rewrite of the moment equations in frequency groups. Our motivation is thatour calculation does not aim in computing the reprocessing of the radiation field spectrum, butrather in a reasonable and spatially inhomogeneous estimate of the photoionization/heating ratesthat will give rise to the 3D He II reionization process.Upon obtaining the solution E .. , E hereafter, we proceed to compute the local ionizationand heating rates as described in § c k are no longer necessaryand are discarded. The obvious advantage of this method is computational speed in the derivationof the photo-ionization/heating rates through an analytical formula. The disadvantage is that theaccuracy is as only good as the cubic interpolation scheme. However, we note that our interpolationscheme does a reasonably good job of describing the processed quasar spectrum obtained with thefull multifrequency calculation of Abel & Haehnelt (1999). In fact, our choice of frequency pointswas strongly guided by the inset spectra in their Figure 1.The distribution of local sources determine the point source emissivity (source function) and the3D distribution of the Eddington factor. Because in our numerical setup all sources emit radiation lca.ucsd.edu/portal/software/mgmpi
14 –with identical spectral slope α q , the calculation of the Eddington tensor via equation (11) yields aquantity that does not depend on frequency ( f ijν ≡ f ij ). The pre-processing involves keeping trackof the sources position and luminosity during the calculation and updating the Eddington factorfunctional form and emissivity source function at the beginning of each data dump calculation.The Eddington factor is initialized at f ijsource = 1 / δ ij within each grid cell containing a source.Between time steps, assumed to be the redshift interval between the data dumps ( δz = 0 . E ( r , t n ) which we use to update the He II photo-ionization rate Γ rad described in § Our simple single species chemistry determines the time step of the evolution. We update theionization fraction ψ ( x , t ) = n HeIII n He through the rate equation:˙ n HeIII = − H ( z ) n HeIII + n HeII Γ − n e n HeIII α ( T ) (13)In Equation (13), Γ = Γ rad + Γ col is the sum of the radiative and collisional ionization rates. Thecollisional ionization is due to collisions of the He II ion with electrons, He II + e → He III + 2 e ,and therefore is proportional to the electron density n e . The collisional ionization rate per electronΓ col /n e is a function of the gas temperature. An analytic fit to the temperature dependence isprovided in Appendix C where we also provide the functional form for the recombination coefficient(both quantities are measured in cm s − . The photoionization rate per baryon is then givenby the equation (Osterbrock 1989) Γ rad ≡ ∞ R ν c E ν hν σ HeIIν dν . Because we find a parametric fit ofthe energy density in frequency space the photoionization can be directly computed as follows.Γ rad = cσ HeII h [ E α q +3 + E P k =0 c k k +3 (1 − − ( k +3) )], where we approximate the ionization crosssection with the power law σ ǫ = σ HeII ( ǫ ) − .To determine the time-step δt n for the advancement of the radiation solution between t n and t n +1 , we follow the procedure below. We collect the photoionization time-scales, τ chem = (Γ rad ) − ,from grid cells that lie in the vicinity of the ionization front and then calculate δt through δt = max ( τ chem , τ lc ), where τ lc is the cell light crossing time-scale, τ lc = δxc − . Cells close the I-frontinterface are “captured” by the criterion ∆ n HeIII n HeIII ( t n ) & .
1. Comparison between the light-crossingand chemical time scales is necessary in the initial moments of the evolution, because in proximityto the source the I-front propagates close to the speed of the light. Evolving the energy density withthe light-crossing time-scale constrains the I-front expansion to subluminal speeds. In addition, theuse of a chemical time-scale close to the source would yield very small time-steps, due to the largenumber of photons, that would in turn slow down the overall calculation.Before we proceed, we need to make the following very important clarification. In this 15 –work, we do not consider contributions to the radiative energy density by He II recombinationsdespite including them in the update of the He II abundance. The significance of the contributingHe II recombinations to the energy density are important when a consistent ionization calculationof all species is performed. In that case, radiative recombinations to energies ǫ < . III ionization frontyield photons with energies ǫ < . I and H I and assist in the advancement of the He II and H II fronts respectively. In oursetup, there are no He II and H II Str¨omgen regions, only He
III ones. Therefore, the photon fluxfrom recombinations on the He
III
I-fronts are ignored and we only consider the chemical effectsresulting from the attenuation and percolation of the individual He II ionizing flux.However, recombinations are included in the abundance update primarily because it can bean important effect for the diffuse helium in proximity to a local source that shut off. The lack ofdirect photons could lead to a rapid recombination of the He III bubble, if no additional radiativeflux reaches that region from another QSO source. Equation (13) can then be rewritten as follows.˙ ψ = (1 − ψ ) τ ion − ψ τ rec (14)In Equation (14), τ ion = (Γ ) − is total ionization time scale and τ rec = ( α ( T ) n e ) − is the localHe III recombination time scale and Γ is the sum of all ionization processes that lead to the forwardreaction He II → He III and we have assumed ψ ≈
0. These processes in our scheme involve thecombination of the direct photoionizations from the QSO sources Γ rad and collisional ionizationΓ col = n e k col ( T ). We integrate equation 14 using backward Euler time differencing where all sourceterms are computed at the advance time t n +1 . ψ n +13 = ψ n + δt n / ( τ ion ) n +1 δt n [ τ ion ) n +1 + τ rec ) n +1 ] (15)In Equation (15), only the photoionization rate Γ rad is available at the advanced time. There-fore, we are forced to initialize the abundance update by computing ( τ ion ) n +1 ≃ rad ) n +1 + n ne ( k col ( T )) n and substitute ( τ rec ) n +1 → ( τ rec ) n . However, upon obtaining the He III number density at t n +1 , n HeIII = ψ n +13 n He we can update the electron density at t n +1 and insert it back in Equation (15).We therefore can improve upon the original estimate by iterating Equation (15) until ∆ n e n e ≤ . §
5, is crude and isnot used in the iterative scheme. The local helium number density at any time is computed fromthe gas density as n He = ρ He m H ≃ ρ m H . The electron density is given by the charge conservationequation n e = n HII + n HeII + 2 n HeIII . Our calculation assumes no change in the ionized hy-drogen density between the value in the original simulation ( o ), and the post-processed value ( ). 16 –Therefore, n e − n oe = ( n HeII − n oHeII ) + 2( n HeIII − n oHeIII ) = ( n HeIII − n oHeIII ). which allows thecalculation of the local electron density at any time from n e = n oe + ( n HeIII − n oHeIII ). In Figure (2), we show volumetric renderings of n HeIII ( r , z ) at two redshift instances. The 3Dvisualization shows the expanding ionized bubbles filling up the cosmic volume due to the combinedeffect of the radiative energy transport and sources being turned on at different parts of the volumeat later redshifts. Individual bubbles of He III may stagnate as they reach their Str¨omgren radiidue to recombinations. However, overall the volume filling factor (VFF) of He
III increases asmore quasars are placed in the computational volume and the percolation between the I-frontsincreases the mean free path of the ionizing photon flux. In left panel of Figure (3), we showa slice through the cosmic volume at z = 2 . II , He III density distributions. Ionizedregions have percolated through the cosmic medium to “open up” the IGM to & . II reionization by such redshifts. In the right panel of Figure (3), we showthe redshift evolution of the VFF, as measured by the fraction of the grid cells with ionized heliumat He III abundance of ψ ≥ − . As the ionized regions begin to merge, assisted by the increasein the QSO number density, the VFF(z) rapidly increases, leading to a value of >
68% at z ≤ . V F F ≈ .
90 is achieved by z ≈ . II reionization. This redshift value compares well with the observeddetermination by Kriss et al. (2001). The solid line in the VFF evolution figure is a spline fitthrough the computed data ( δz = 0 . δz = 0 . δx ) /V , where δx is the grid resolution and V is the volumeof the computational box. The evolution of the VFF shows a rapid increase in He III at redshifts z .
4, following an earlier epoch of apparent stagnation.An alternative way to illustrate He II reionization is to plot the redshift evolution of the volumeaveraged abundance fraction. In the left panel of Figure (4), we show that the mean mass fractionin He II ρ He II /ρ , drops significantly at z .
4. For reference, we also show the mean fraction inH I ρ H I /ρ , undergoes a steep drop at z ≃ . II reionization epoch is much more extended than that of H I . Arapid drop in the H I fraction occurs between z = 7 − . I mass fraction drops by 4 dex. By visualinspection of the mean fraction in He II from Figure (4) we can determine that it drops by ∼ z = 3 . − .
5. When the difference in the redshift interval is converted to cosmic timeit yields that a similar drop in the mass fraction takes approximately eight times longer to occurin the case of He II when compared to H I . This is consistent with the conclusion reached in § II reionization is attributed to higher recombination time scales andless available ionizing photons per baryon for a stellar dominating UV background.In the right panel of Figure (4), we plot the He II mass fraction versus the local overdensityby logarithmically binning the latter quantity and computing the mean and median He II from thecells with gas density within the bin. Such a graph is intended to show the trend between the twoquantities. For reference, we plot the trend between the two quantities from a standard cosmologicalsimulation, where all ionizations are computed self-consistently due a uniform UV background inthe optically thin approximation.The simulation is terminated at z ≃ . V F F ≃ .
9. The reason for suspending thecalculation is that when the cosmic volume is effectively transparent to the ionizing radiation,the assumptions underlying equation (7) no longer apply. Specifically, free streaming photonsmay cross the volume unimpeded by absorption and their mean free path can become comparableto the size of the horizon. The latter effect requires keeping the cosmological dilution term inEquation (12), which was ignored. From the numerical perspective, solving a parabolic equationwhen the conditions call for a hyperbolic one, can create local superluminal speeds of the ionizingfront if the chemical time scale becomes smaller than the cell light crossing time. Concluding thecalculation at the end of the opaque phase of He II still addresses the main question investigated inthis work; what is the primary mechanism that leads to the He II reionization. We conclude, that arising population of QSO sources, assisted by the gas dilution due to cosmic expansion, clumpinessdue to structure formation and the pre-ionization of neutral hydrogen which allows for the almostexclusive usage of the He II ionizing radiation are a set of physical conditions that reproduce theepoch of He II reionization by redshifts z ≤ .
5. Late Heating due to He II Reionization5.1. Physical Considerations
Photoionization of He II at an epoch later than that of H I releases an additional one photo-electron per ionization. We can estimate the mean energy of such electrons due to absorption ofphotons with energies ≥ . ≡ II and compare it to the value obtained in the case of ≥ . ≡ I atoms. For simplicity, we will ignore the geometric attenuationand optical depth effects in the radiative flux due to local sources and gas opacity and assume thatthe local mean intensity of the radiation is the same as the emitted, with a spectrum J ǫ = J ǫ − α q .As before, in this notation ǫ = hνhν . The photo-heating rate, in ergs/s, due to electrons ejectedfrom He II atoms, can then be computed from ¯ G HeII = ν R ∞ (4 πJ ǫ /ǫ )( ǫ − σ ǫ dǫ . Substitutingfor the power law of the mean intensity, for σ ǫ = σ oHeII ǫ − we get ¯ G HeII = πν J σ oHeII − αq ( α q +2)( α q +3) . 18 –In a similar manner, we can derive ¯ G HI = πν J σ oHI ( α q +2)( α q +3) . If the sources responsible for the He II andH I ionization are the same then the J amplitude and the α q spectral slope are identical in bothequations which allows for the derivation of the following ratio:¯ G HeII ¯ G HI = σ oHeII σ oHI − α q (16)Substituting for σ oHeII = σ oHI we get ¯ G HeII ¯ G HI = 4 − α q . The total photoheating rate per unitvolume is proportional to the number density of the absorbers which yields the ratio betweenHe II and H I photoheating rates to be G HeII /G HI = n HeII n HI − α q . The ratio of number densities isestimated in § n HeII /n HI = S , where S is the softness parameter. The ratio of thephotoheating rates then becomes G HeII /G HI = S − α q ≃ α q +1 − α q = 5 /
3. The latter factoris indicative of the degree of temperature increase due to the thermalization of the photoelectronsejected by He II ionizations when compared to the temperature inferred by HI ionization alone.In deriving the above value, we made the assumption that hydrogen and singly ionized heliumare photoionized simultaneously by the same ultraviolet field. The effects of distinct ionizationepochs can however be modeled in the above relation if we adopt the scenario of hydrogen ionizationby a soft radiation background ( α (2) ∼
5) and of He II ionization by a hard one ( α (1) ∼ . z ≃ G HeII G HI =
53 ( α (2) q +2)( α (1) q +2) ≃ for α (1) q = 1 . α (2) q = 5. In conclusion, the thermalization of the photoelectron in the He II + γ → He III + e − reaction can be an important determinant of the intergalactic mediumtemperature.Photoionization models based on hydrogen ionization alone predict a temperature of the inter-galactic medium of T IGM ≈ . × K. However, Ly α forest observations yield lines with medianbroadening widths of b = 26 −
36 km/s (Carswell et al. 1987, 1989; Zhang et al. 1997; Dav´e et al.1997) in the intermediate redshift range of z=2-4. Because the thermal and differential Hubble flowcomponents in the total HI line broadening are of the same order of magnitude (Zhang et al., 1998;Aguirre, 2002), one can infer a thermal width range in the HI Ly α forest between b th = 13 − b th would require the temperature in the intergalactic medium tobe T IGM = 20 , − ,
000 K, a factor of ∼ . − . II ionization photoheating due to a hard ultraviolet spectrum. 19 – The post-processing of the simulation data dumps involves an update of the local gas temper-ature by solving for the perturbed temperature in the thermal energy equation. In the simplestapproximation, we will assume that the change in the net rate of the cell thermal energy densityis due to the balance between the heating by the ejected photoelectrons from He II ions and theHe III recombination cooling, γ − kρµm p ∆ δTδt = G − Λ, where γ = is the gas adiabatic index, ρ thelocal gas density and µ the mean molecular weight. All quantities have local proper values. Inaddition, in the notation followed here G = n HeII G r is the total photoheating rate, measured in ergs/s/cm , which is the radiative rate G r ≡ ¯ G mentioned above. The cooling rate, Λ, is a functionof the local gas temperature and density and is equal to Λ ≡ n e n He III L ( T ). L denotes the coolingfunction due to He III recombinations and is given by L ( T ) = 3 . × − T / T − . (1 + T . ) − inunits ergs cm /s (Cen 1992), where T n = T / n .Recombination He III cooling is only one out of several processes that cool the cosmic gas,the list of which is described in detail in Anninos et al. (1997). The cooling rate coefficients,in parametric fits that depend on the local gas temperature, due to excitation, ionization andrecombination of the primary chemical species along with bremsstrahlung and Compton coolingare fully incorporated into the ENZO code. We expect that He
III cooling would dominate thecooling processes inside the He
III bubbles when compared to the corresponding He II recombinationcooling, because the He II abundance is significantly reduced there to fractions ψ ≃ − − − .In addition, the He III recombination cooling rate, along with H II recombination, dominate at thelow temperatures found in underdense cosmic regions T . K over all other types. In a schemewhere hydrogen is already almost completely ionized at δ . χ ≈ ψ ≈
1) the thermal balance would be controlled by the heating and cooling ofthe latter species’ ionization.However, our ability to accurately post-process the temperature field in our scheme is limitedby the fact that excitation and collisional ionization cooling from processes involving H I , He I andHe II species dominate the cooling curves at temperatures T & K, while exhibiting very steepprofiles in the temperature range T ≃ − K. An increase in the gas temperature by thephotoelectrons ejected in the He II radiative ionization would cause an increase of the cooling co-efficients. Combined with the increase in the electron number density this shift in the thermalbalance could significantly increase the cooling rate even though the fractional abundances in H I ,He I and He II are small. The end result is that, by only including He III recombination cooling inrecomputing the gas temperature, we may overestimate the value of the adjusted temperature. Thisupper limit in the temperature estimate implies that the recombination time-scale in Equation (15)is also an upper estimate and that the overall propagation of the cumulative He
III ionization frontis faster than in the case of a self-consistent calculation. However, such calculation would requirecomputing the ionization and thermal balance of all species self-consistently and is beyond thescope of this paper. A numerical scheme is under development that will allow us to do this in thenear future (Reynolds et al., in prep ). Therefore, we conclude that the epoch of He II reionization 20 –may be placed at a later redshift than the one we calculated here ( z reion ∼ .
5) even though wenote that such adjustments could be reversed or may not be necessary if the diffuse recombina-tion radiation is added in the calculation. Such radiation would further ionize HI and HeI speciessuppressing their contributions to the cooling curve. That may explain why, even under all of theassumptions and approximations that we allowed and followed in this work, the end result of ourpredicted redshift evolution of the He II opacity correlates well with observations of the He II Ly α forest as we shall show in § T ( o ) and T (1) denote thecell temperatures before and after the presence of He II ejected photoelectrons, then we approximatethe change in the thermal energy equation as follows: δT (1) δt ≃ δT ( o ) δt + ( γ − m p µ (1) kρ ( G (1) − Λ (1) ) (17)In an explicit, time-discretized form, where the original temperature at time t n is obtained byinterpolation between the logarithm of the local temperature in the two data dumps with redshiftsthat bound the time evolution ( z ≤ t n ≤ z ), Equation (17) becomes: T (1) n +1 − T (1) n = T ( o ) n +1 − T ( o ) n + ( δt ) n ( γ − k ( m p µρ ) n + ( n n + HeII G n + r − n n + e n n + He III L n ) (18)In Equation (18), time-centered quantities are computed as X n + = 0 . X n +1 + X n ) andonly for variables that we know their value at the forward time t n +1 . Therefore, L is computedat time t n because it is a function of temperature which is unknown at t n +1 . The update intemperature occurs after all principal quantities of E n +1 , n HeII , n He III and n e are computed at t n +1 by advancing the local solutions at the implicit time-step of the radiation field discussed in § m p µρ = P n i is equal to total number density of the cosmic species plus electronsand is computed as follows: m p µρ = [ ρm p + n e ] − ≃ [0 . · − Ω b h (1 + z ) + n e ] − , where theexpression for n e was provided in § III bubble at t n then T (1) n = T ( o ) n and n nHe III →
0. If no direct ionizing radiation reaches that grid cell by t n +1 then n n +1 He III → T (1) n +1 = T ( o ) n +1 . In a cell where thermal equilibrium was reached duringthe original calculation between t n +1 and t n then during post-processing T (1) n +1 = T (1) n is achievedonly when the heating and cooling terms balance out.A point of concern is that there was no physical reason to suppress the collisional ionization ofHe II in the original simulation. Even if there are no ionizing sources capable of radiatively ionizingHe II , a local temperature of T ≥ . · K, found in overdense regions, may be enough to 21 –collisionally eject the 54.4 eV bound electron in the He II atom. Consequently, the temperatureevolution T ( o ) ( z ) includes such an effect and can skew the post-processed evolution T (1) ( z ). This isevident if we assume that due to collisional ionization in the pre-processed data n He III = 0 and at t n , X (1) n = X ( o ) n Equation (18) would then predict T (1) n +1 = T ( o ) n +1 + ( δt ) n W ( n ( o ) HeII , n ( o ) He III , n ( o ) e , T ( o ) )where W ≡ ( γ − k ( m p µ o ρ ) n + × ( − n n + e,o n n + He III ,o L n ( T ( o ) )The last equation would unnecessarily recompute the temperature in regions which lack ra-diative input but have significant collisional rates in He II ↔ He III and therefore, n He III ,o = 0.To account for such discrepancy in the temperature evolution, each updated local temperatureis adjusted as T (1) n +1 = T (1) n +1 + ( δt ) n W ( n ( o ) HeII , n ( o ) He III , n ( o ) e , T ( o ) ). Although in doing so we improveupon the temperature evolution, the collisional effect can never fully be readjusted because ofinterdependency between all of the physical quantities. In Figure 5, we show the results of our temperature calculation. In the left panel, we plot theevolution of the mean temperature in overdensities log ( δ ) = 0 − z . . II photoheating due to the rising population of QSO’sin the cosmic volume. The overdensity interval was chosen in order to show that He II reionizationis a cosmic event that primarily affects the diffuse and mildly overdense IGM where the the tem-perature can get increased by about a factor of 2 at the end of the calculation. One should considertwo effects that support this conclusion. In regions of significant overdensity and therefore darkmatter potential, gas is heated by the gravitationally controlled free fall compression to temper-atures T & K. Therefore, the effects of photoheating due to the photoelectrons ejected byradiative ionizations of He II only result in an insignificant fractional change. In addition, the largeelectron density and gas density can be a source of large optical depths and radiation trappingdue to increased recombination time scales for a modest increase in the temperature. As a resultthe radiative energy density can significantly drop in such regions which furthermore reduces theefficiency of He II ionization.The optical depth effects are demonstrated in the right panel of Figure (5), where we show(z=2.5) the median temperature vs. density (solid line) in the scatter plot between the two quan-tities on the simulation grid. We obtain the curve by binning the gas overdensity in logarithmicintervals and computing the median gas temperature within each bin. For reference, the relationin the original calculation is also shown (dashed lines). On the left panel of Figure (6), we plot theratio between the optically thick optically thin calculations. The dashed lines shows the effect ofour postprocessing on the gas temperature. He II photoionization contributes primarily to the gastemperature at the mean and low gas densities. The effect is diminished at higher overdensitieswhere collisional ionization dominates. In addition to comparing to the original calculation we 22 –also compare to a simulation where the chemical species abundunces are self-consistently computedduring the simulation due to the same homogeneous UV background and in the optically thin limit.The dot-dashed curves on the right panel of Figure (5) and the left panel of Figure (6) trace the tem-perature density relation in that case and the ratio to our optically thick calculation respectively.We find that our postprocessed temperature is higher by a factor of ≃ . II ionizations in the redshift interval that correspondsto the rise in the number density of sources emitting hard radiation. The physical implication isthat the late increase in temperature may well be the reason why the galaxy luminosity functiondecreases at z .
4. A raise in the mean temperature due to the cosmic evolution of QSO sourceswould increase the Jeans mass threshold by a factor of 2.2-5 if we adapt an average increase inthe temperature between 1.7-3 ( M jeans ∝ T ). The latter would in turn suppress the furtherformation of dwarf galaxies in the cosmic volume. However, we find that, according to the rightpanel of Figure (5), the fractional increase in local gas temperature is on average not that large athigher overdensities. If we exclude cosmic neighborhoods that are in close proximity to the localUV sources, where the temperature increase is significantly greater due to the high radiative energydensity, then in collapsed structures at δ ≥
100 which are increasingly self-shielding, photoheatingis not as effective as shock heating due to gravitational collapse. Nonetheless, a firm conclusionon that effect is not possible in this work, due to the very coarse grid resolution that does notadequately resolve the aforementioned structures.
6. Signatures of He II Reionization6.1. Synthetic Flux Spectra
The updated He II density and gas temperature can be used to study the effects of this in-homogeneous reionization scheme on the transmission of the intergalactic medium transmission inthe He II Ly α restframe wavelength. The objective is to compare with the He II transmissivity ob-tained from analyzing the currently available observed lines of sight. We synthesize spectra of theHe II Ly α absorption at the rest wavelength of 304 ˚A along 300 random lines of sight (LOS). Thenumber of LOS was chosen in order to yield a less than 1% fluctuation to the mean transmittedflux by the end of the calculation.The synthesis method is described in detail in Zhang et al. (1997). In addition to He II , wealso compute the transmitted flux of the corresponding H I forest at the rest wavelength of Ly α II and H I and therefore maps onto 23 –the same grid location for the two redshifted wavelengths. The redshift interval of the spectrumoutput is set at ∆ z = 0 .
2. In addition to this interval being the redshift output interval of thehydrodynamic simulation it is a long redshift path that minimizes the sightline to sightline varianceby forcing the transmission path to cross and wrap through the volume boundaries along the samedirectional vector about ∆ z ∆ z cube times. ∆ z cube = L H ( z ) /c is the linear size of the cube in redshiftunits where H ( z ) = 100 h ((1 + z ) Ω M + Ω Λ ) is the Hubble constant at redshift z and L=100 Mpccomoving. At z = 2 . H ( z = 2 .
5) = 249 . z cube = 0 .
083 and ∆ z ∆ z cube = 2 .
4, the number of wraps through the computational volume.The output from the hydrodynamic simulation is saved every δz dump = 0 . z and z . The output spectrum istherefore centered at 0 . × ( z + z ), where z = z − δz dump . The input fields are species density(H I and He II ), gas temperature and the three peculiar velocity components, all assumed to befrozen in the comoving frame of reference. However, we allow proper evolution along the redshiftpath of the sightline for the densities ( ∝ (1 + z ) ) and velocities ( ∝ (1 + z ) − ). The spectra arecomputed along lines of sight that sample the computational volume continuously for z ≤ .
1. Wedo this by restarting the calculation in the new redshift interval from the mesh location wherethe previous calculation stopped and continue along the same directional vector. The initial pointof origin is randomly selected. After the completion of each ∆ z = 0 . I (top left panel) and He II (bottom left panel)Ly α transmission at z=2.5 and compare with the results from an Enzo simulation that computesabundances self-consistently in the optical thin limit (top and bottom right panels) using a UVBfrom a mix of evolving quasar and galaxy populations (Haardt & Madau 2001). Both transmis-sion lines are casted from the same initial point and along the same directional vector. We usethe updated temperature and the unperturbed neutral hydrogen abundance for calculating theH I transmission spectra. Although there are few visible differences at first glance, the H I spectrumin the postprocessed temperature case has a continuum flux level below that of the optically thin,self-consistent case. Since the mean transmitted flux is sensitive to the number of pixels close tothe continuum at low redshifts the line of sight average transmitted flux in the post processedcalculation is below the value obtained from the self consistent calculation by about ≈ ≈ . α He II transmission versus observed wavelength. Below, we plot the instantaneouscomoving distance of each velocity pixel to the closest UV source versus the comoving path length ofthe sightline. The last curve is a joint ensample of parabolic curves each having a minimum, markedwith crosses, that corresponds to the position of nearest distance to the ionizing source. Each timeanother UV source becomes closest to the trajectory, the curve is marked by the beginning of ananother parabola. We note that at the location of closest proximity for that sightline, at ≈ ≈ ≈
550 km/s). Due to negative skewness, the mean is shiftedto the right at a value of ≈
14 Mpc ( ≈
750 km/s). The negative skewness is most likely due tothe fact that the distribution of the UV point sources on the grid is not isotropic but are clusteredaccording to the clustering properties of the host dark matter halos. From the distribution in theright panel of Figure (8) we can compute a typical range of impact parameters in the sightlines of ≈ −
22 Mpc comoving contained within the values of the distribution at 1/e the peak value. Weare interested to determine whether the transmission properties of our sightline sample is in anyway biased by the proximity to the ionizing sources or if we are mostly sensitive to the values ofthe ionizing flux at the mean level.Detailed analysis of the proximity effect around each source is beyond the scope of this paper.However, we can impose an upper bound limit based on the highest luminosity quasar in thecomputational volume at z=2.5, L max = 2 × ergs/s where L designates continuum ionizingluminosity above 4 Ryd. An estimate of a proximity size distance can be obtained from the 25 –volume averaged energy density above the He II ionization threshold, ¯ E . In our calculation wehave computed ¯ E = 7 . × − ergs/cm − which corresponds to a volume averaged He II ionizingflux equal to ¯ F = c ¯ E = 2 . × − ergs/s/cm − . The condition that the geometrical attenuationof the central emission equals the total flux through the surface of a sphere centered at the sourceyields an upper bound for this proximity distance of d properHeII ≤ ( L max / π ¯ F ) ≃ .
86 Mpc properor d HeII . ≈ II Ly α transmission value.On the left panel of Figure (9), we show a quantitative estimate of that effect. There, wescatter plot the average impact parameter of the sightline versus the mean transmission in eachsightline (cross symbols). The degree of their linear correlation, measured by the Pearson correlationcoefficient, is rather small, r = − .
24. The figure does show that there is non-negligible negativecorrelation between mean transmission and the average proximity to the distribution of ionizingsources for each sightline. A power law fit to the data in the form ¯ F los ∝ ¯ d − sip where ¯ F los is thesightline specific mean flux and ¯ d ip is the average impact parameter yields s = − . ± .
04 and isplotted over the individual points. The curve lies within the two standard deviations levels of themean transmission. Therefore, we conclude that the proximity effects have no statistical weight onthe sightline transmissions. That is consistent with our estimate of the upper limit on the proximitydistance ( ≃ II Ly α absorbed flux at that distance fromthe point source (points). The data show a tail of low absorption at small impact parameters atdistances . . & . σ level. To get a clear view of the trendwe bin the horizontal axis into constant logarithmic bins and compute the median (red histogram)and mean (blue histogram) per bin. A Fermi-Dirac function in the form f ( x ) = exp ( ( do − x ) µ is thenfitted to the median histogram with do = 3 . µ = 0 .
63 Mpc (skinwidth) for impact parameters up to 10 Mpc. The point of half maximum is consistent with theprevious estimate of the proximity effect distance based on the luminosity and mean intensity of theUV field. The figure clearly shows that the effect of source proximity on the local absorbed flux arerelevant at comoving distances of . II Ly α absorption in the IGM and becomes significant at impact parameters . . σ level of the mean absorption. 26 – The straightforward average of the flux, ¯ F = Npx · nlos Σ nlosj Σ Npxi F ij , from all pixels and lines ofsight per redshift interval, is a measure of the opacity of the cosmic volume due to the absorption bythe particular chemical species (He II and H I in this case). In the notation followed, N px = 30 , R z = δz =0 . / , ≃ . × − , and nlos = 300 is the number of random lines of sight. The redshiftresolution is fixed in our calculation, which results in a redshift dependent spectral resolution of R λ = λ ∆ λ = zδz = 1 . × (1 + z ). At z = 3 this yields a spectral resolution of 50 times higherthan the designed value of the spectrograph aboard FUSE.The evolution of the mean transmitted flux is typically represented by the effective opticaldepth, defined as τ eff = − ln ( ¯ F ). As discussed in Paschos & Norman (2005) (PN), the effectiveoptical depth is biased by high transmission gaps in the pixel flux distribution and therefore willsystematically yield lower values when compared to the mean optical depth. The latter is definedas the raw average of the pixel optical depth per redshift interval, τ mean = Npx · nlos Σ nlosj Σ Npxi τ ij ,where τ ij is the pixel optical depth at redshift z i in the line of sight index LOS = j . In Figure (10),we show the redshift evolution of the the optical depth of the 304 ˚A line using both representations,mean (left panel) and effective (right panel), which suggests a rather smooth evolution leading toHe II reionization at z ≃ .
5. Error bars show the 1 σ standard deviation of the mean flux andoptical depth values between different lines of sight. In the effective optical depth case, standarddeviation at more than 100% the mean flux value at z ≥ F min = 10 − was therefore imposed in order to be able to compute thenatural logarithm.The error bars also indicate a large degree of variance in the data at redshifts z & . − I transmission during hydrogen reionization.There, the large degree of variance during reionization was attributed to the presence of hightransmission gaps associated with underdense regions in the IGM which ionize first. In this pictureof inhomogeneous He II reionization, the high transmission segments at high redshifts are associatedwith He III bubbles that the lines of sight intersect as they are cast through the simulation box.However, this is not inconsistent with the conclusions reached in PN because the He
III bubblesprimarily extend in underdense to mean density cosmic regions.Estimates of the observed He II effective optical depth, shown in Figure (10), are comprisedof the six known and analyzed sightlines todate: HS1700+6416 (Fechner at al. 2006; Davidsenet al. 1996), HE2347-4342 (Zheng et al. 2004b; Kriss et al. 2001), HS1157+3143 (Reimers etal. 2005), PKS1935-692 (Anderson et al. 1999), Q0302-003 (Heap et al. 2000) and SDSSJ2346-0016 (Zheng et al. 2004a). Although very few in number they suggest a lack of He II Ly α foresttransmission at redshifts above z ≃ II reionization. Upper bound estimates suggest a steeply rising optical depthabove z ≃ . z ≃ . z ≃ .
7. In thiswork, we have calculated that between z = 2 . − .
4, at ¯ z = 2 .
5, a mean transmitted flux value at¯ F = 0 . ± .
002 which corresponds to τ eff = 1 . ± . σ/ √ N los . Observed effective optical depthestimates are τ effobs = 0 . ± .
01 for HE2347-4342 averaged over z ≃ . − . τ effobs = 0 . ± .
34 and τ effobs = 1 . ± .
18 at ¯ z ≃ .
45 for HS1700+6416 (Fechner et al. 2006).If reionization is completed by z ≃ .
5, as suggested by the observed lines of sight, that wouldcorrespond to an optical depth of τ eff ≃
1. Our calculation is about 1.5 σ above that value atthat redshift. Furthermore, the computed redshift evolution of the optical depth is smoother thansuggested by observations. We attribute such differences to the limitations of our treatment andmost notably ignoring diffuse emission due to recombinations to the ground state of He II . Suchrecombinations increase the ionization of He II and therefore lower the opacity. The slope of theredshift evolution is affected by such omission because of the diffuse emission’s dependency on thegas temperature, ∝ T − . At earlier redshifts, the temperature is lower and therefore recombinationsto the ground state may contribute a significant amount of He II ionizing radiation. Nonetheless,due to the lack of a self-consistent calculation we do not know how large that effect would be onthe slope of the optical depth redshift evolution. We do anticipate that the effect diminishes withredshift due to photoheating which raises the IGM temperature.From the right panel of Figure (10) we can infer a general agreement between the observed andcomputed values at redshifts that evidently sample the tail of the He II reionization epoch. Thatleads to the conclusion that the final opacity distribution in He II is largely insensitive to the detailsof our numerical setup. Direct photo-ionizations due to a rising number density of QSO sourcesappears to be adequate to yield a value of the He II abundance close to the observed one at lateredshifts. The largely ionized He II by z=2.5 in our computational volume and the update in the gastemperature due to photo-heating, allows for the derivation of two standard statistical propertiesof a Ly α forest, the column density and b-parameter distribution. Our grid resolution is too coarseto draw any absolute conclusions about these distributions. In order to resolve the H I Ly α forest,a grid resolution of ∼
40 kpc is required (Bryan et al. 1999). A larger grid resolution may not berequired for He II lines if they primarily form in underdense regions. However, comparison with aresolved H I forest is desirable. Our intent in this work is to look at the relative differences betweenthe two forests and draw conclusions that may stand the test of an improved grid resolution in afuture simulation.Between z=2.6-2.4, at ¯ z = 2 .
5, we identify the He II and H I Ly α lines and compute theircolumn density and b-parameter width. The line identification was performed using the method 28 –described in Zhang et al. (1998). We count lines per constant logarithmic column density bins(∆ N X = 0 .
25) and normalize per column density and redshift interval, dz = 0 .
2, for the twospecies, d NdzdN X (X = H I , He II ). The result is shown in the left panel of Figure (11). TheH I distribution flattens at LogN HI ∼ . LogN
HeII ∼ .
5. In addition, the He II column density distribution shows adecrease in the slope and large error bars as the column density increases. The error bars refer tothe 1 σ standard deviation of the identified line counts per logarithmic bin. By LogN
HeII & . II that are not fully ionized by z=2.5. Such regions occupy about 10% of the computationalvolume and are a source of large optical depths and therefore column densities in lines of sight thatpass through0 them. The two profiles suggest that for column densities in the range 13.5 - 15.5we identify 10-100 times more He II than H I lines. This result is consistent with a high resolutioncosmological simulation that resolves the two forests (Zhang et al. 1997).Power law fits to the column density distributions in Figure (11), f X ( N ) dN X = βN − α X X dN X ,yield α HI = 1 . ± .
14 and α HeII = 1 . ± .
12 in the logarithmic column density range of 13.5 -15.5. The error estimates are derived from the least squares fit and do not take into considerationany propagated error in the individual column density bins which is sensitive to the bin size. TheLy α forest is well studied in the literature, however results that are pertinent to the slope of thecolumn density distribution typically refer to either the range of logN HI ( cm − ) = 12 . − . N HI = 10 . − . cm − of α HI = 1 . − .
85 forredshifts between z=2-3. The slope for the H I column density distribution obtained in this work isslightly above. That is primarily due to the poor grid resolution that underestimates the numberof H I absorbers with high column densities. To a lesser extent, we also expect a steepening inthe distribution because we overestimate the H I opacity for lower column density absorbers in ournumerical setup. However, we note that within the error estimate of the power law fit, the valueof the slope we computed here overlaps with the range of published values.The estimated slope of the He II Ly α column density distribution, in the optical thin limit,using uniform photoionization models, ranges between α HeII ≈ . − . logN HeII = 13 . − . z = 2 . α HeII = 1 .
41 which is below the previous estimates. As we shallshow below the He II absorbers correspond to physically extended IGM structures and thereforethe He II column density distribution is less sensitive to our coarse grid resolution than hydrogen.The smaller slope is due to the opacity effects of the inhomogeneous ionizing radiation. Overdenseregions will have higher opacity to He II ionizing radiation compared to a uniform photoionizationmodel. That in turn results in lines migrating to higher column density bins which softens theslope of the distribution compared to an optically thin calculation. 29 –We now proceed to investigate the effects of non-uniform photoheating in the distribution ofthe line broadening widths which is the observable that closely traces the thermal state of the IGM.The b-parameter distribution, computed here as the fraction of lines per km/s is shown for the twospecies on the right panel of Figure (11) (solid lines). The distributions were computed from lineswith column densities in the range LogN X = 13 . − . b HImed = 34 . − and b HeIImed = 28 .
16 kms − for the H I and He II lines respectively. For comparison, the medianvalues in the uniform UVB case are b HImed ( u ) = 32 .
98 kms − and b HeIImed ( u ) = 26 .
89 kms − . When wecompare the non-uniform and uniform results we note that the peaks of the distributions are offsetby ≃ .
25 km/s in both species. This is expected according to the right panel of Figure (5) wherethe median temperature per logarithmic overdensity bin is plotted. Differences in the overdensitydependent median temperature distribution arise at overdensities ∆ .
3. If we adopt a relationbetween H I column density and overdensity in the form of ∆ ≥ N HI ) / ( z ) − (Schaye et al.2003) then we can compute that at z=2.5 a H I column density of 10 . cm − is due to overdensitiesof ∆ & .
5. At such overdensities we predict the biggest difference in the temperature between thenon-uniform versus the uniform cases. The increase in temperature by a factor ≃ . II lines have about ≈
82% the broadening width of the H I lines. For pure thermal broadening,the heavier by a factor of 4 helium atoms would yield a width only half the corresponding size ofthe hydrogen line, b HeIIth = b HIth . Our calculation suggests that the broadening due to the Hubbledifferential flow and peculiar velocities dominate the line formation in the He II forest. That in turnwould indicate that the He II lines primarily form in underdense and cosmic mean density regions.Our result is consistent with the conclusions in Zheng et al. (2004) where they calculated thatalong the HE2347-4342 transmission line b HeII ≃ . b HI which also supports a Hubble dominatedabsorption line broadening for He II . η -Parameter Evolution We conclude this section, by computing the redshift evolution of the η -parameter throughthe R-factor, which are defined in § η -parameter in the range 3 . ≤ z ≤ . δz = 0 .
2, averaged over all lines of 30 –sight. Direct computation of the quantity not only requires knowledge of the column density foreach individual line but also knowledge of the local association between the H I and He II absorptionfeatures. However, because the transmission of the He II Ly α line exhibits a trough at redshifts z &
3, obtaining the column density value for individual He II lines is not feasible there due tothe inability to deblend the absorption features in the spectrum. In addition, whatever correlationbetween the H I and He II lines exists at z . I lines canbe associated with He II absorption features and vice versa (Kriss et al. 2001). Our biggest problemthough in this calculation lies in the separate treatment of hydrogen and helium ionization. Thelatter coupled to a low resolution simulation can be a source of significant bias against the truecorrelation between the He II and H I absorbers.Therefore, we will approximate here the line of sight and redshift interval average of the η -parameter as < log ( η ) > ≃ < log (4 × τ HeII τ HI ) > = < log (4 × R ) > where R stands for the R-factor.Through this approximation, the identification of individual lines and the derivation of columndensity through their gaussian profiles are not required. Instead, in each redshift interval wecompute the average of the ratio between the local He II and H I pixel optical depths along thetransmission line. The values for each line of sight are then averaged to obtain a mean value for the η -parameter. We plot the redshift evolution of such calculation in Figure (12). Because the errorsdue to pixel averaging are large, they are ignored. Instead, the error bars in the figure represent the1 σ standard deviation due to line of sight averaging. For reference, we overplot the computed valueof the η -parameter in the spectrum of HE2347-4342 from Kriss et al. (2001) and also data fromFechner et al. (2006) (blue circle) towards HS1700+6416. We obtained the latter by computingthe mean logη from the published data in the range z=2.3-2.75. At z=2.5 our computed value logη = 2 . ± . II Ly α optical depth at that redshift. As with the optical depth calculation, the limitations ofour treatment result in an overestimate. However, we came close enough to the observed valuesat redshifts later than the He II transmission trough, to reinforce our conclusion that our limitedcalculation is able to reproduce the much of the observed signatures of He II reionization.
7. Summary & Conclusions
We have simulated the late inhomogeneous reionization of He II by quasars and the attendantphotoheating of the IGM including opacity effects. We post-process baryonic density fields froma standard optically thin IGM simulation with a homogeneous galaxy-dominated UV backgroundwhich reionizes H I and He I at z=6.5. Quasars are introduced as point sources throughout the100 Mpc simulation volume located at density peaks consistent with the Pei luminosity function.We assume an intrinsic quasar spectrum J ( ν ) ∝ ν − . and a luminosity proportional to the halomass. We evolve the spatial distribution of the He II ionizing radiation field at h ν = 4, 8, and16 Ryd using a time-implicit variable tensor Eddington factor radiative transfer scheme which wedescribe. Simultaneously, we also solve for the local ionization of He II to He III and the associated 31 –photoheating of the gas. This distinguishes our calculation from the work of Sokasian, Abel andHernquist (2002) who simulated inhomogeneous He II reionization but did not study the thermalevolution of the gas.We find that the cosmic evolution of the QSO population causes the individual He III regions tooverlap and subsequently ionize 90 % of the volume at He
III ionization fractions ψ = n HeIII n He & − by z ≃ .
5. In addition, to the He II and He III number density calculation, we also update the localgas temperature due to the thermalization of the photoelectrons ejected by the HeII ionizationprocess. As expected, the temperature is higher compared to the unprocessed simulation withthe difference increasing rapidly at redshifts z .
3. Relative to a self-consistent optically thinsimulation where He II is also photoionized by a homogeneous UV background, we find that opticaldepth effects result in an increase in the temperature of the intergalactic medium at the cosmicmean density level by a factor of ≈ . z = 2 .
5. The results of our temperature calculation areconsistent with analytic and numerical predictions of the HeII heating effect (Abel & Haehnelt 1999;Schaye et al. 2000; McDonald et al. 2000). Finally, we trace the redshift evolution of the He II Ly α transmission using randomly casted synthetic spectra through the simulated volume. The analysisof the mean transmission allows for the derivation of the effective optical depth of the 304 ˚A line.We have calculated that at ¯ z = 2 . ± .
1, the average of pixel flux among all lines of sight yieldsa value for the mean transmission of ¯ F = 0 . ± .
002 which corresponds to τ eff = 1 . ± . z = 0 .
2) is longer than the size of the simulation in order tominimize the transmission variance across the redshift path. We compute the error only due to thesightline to sightline variance and find that at z=2.5 the 1 σ standard deviation is 11% the meanLy α flux. When we compare to estimates from He II forest spectra observed with the FUSE (FarUltraviolet Space Explorer) satellite and HST at the same redshift interval we find that the ourvalue of the effective depth is comparable but slightly above the observed values. We attribute thedisagreement to our approximate treatment of the inhomogeneous ionizing radiation that ignoresthe diffuse component and only focuses on the point sources’ input.In addition to the mean transmission estimate, we also compute the column density and b-parameter distribution of the identified He II Ly α lines and compare them to the correspondingstatistical properties of the HI Ly α lines. We find that an optically thick calculation results inextra heating that shifts the b-parameter distributions to higher broadening widths. In additionand within the limitations of our coarse numerical resolution, which does not resolve the H I forest,our calculation shows that the median He II b-parameter at ¯ z = 2 . II lines primarily form in underdense to cosmic meandensity regions.ACKNOWLEDGEMENTSThis work was supported by NSF grants AST-9803137, AST-0307690 and AST-0507717. Thesimulations were carried out on the TITAN IA64 Linux cluster at the National Center for Super- 32 –computing Applications and on the IBM Blue Horizon system at San Diego Supercomputing Centerunder LRAC allocation MCA098020. A. Species Concentration Equation in an Expanding Universe
For a single species medium and in proper coordinates the rate equation that governs theabundance of ionization state i is given by the following equation:˙ n i = − aa n i + Σ j β ij n j (A1)In Equation (A1), n j denotes the number density of ionization states of the species in ionizationor recombination coupling to state i. The ratio ˙ aa is a function of redshift and is equal to the redshiftvalue of the Hubble constant.Dividing by the total number density of the species, n T , we get:˙ n i /n T = − aa n i /n T + Σ j β ij n j /n T ⇒ ˙ n i /n T = − aa y i + Σ j β ij y j (A2)The LHS of Equation (A2), can be rewritten as follows:1 n T ∂∂t n i = ∂∂t ( n i n T ) − n i · ∂∂t ( 1 n T ) = ∂∂t y i + 1 n T ˙ n T n i = ˙ y i − aa y i (A3)Plugging the result back into Equation (A2), we get ˙ y i = Σ j β ij y j , which shows that whensolving the rate equation using ionization fractions, instead of number densities, the equationis independent of the Hubble expansion term. However, coefficients β ij that are pertinent torecombinations to state i depend on the electron density which is computed on the proper frameof reference. B. Effect of Late Reheating on Hydrogen Ionization Balance
In this calculation, we estimate the neutral hydrogen fraction if an additional source of ion-ization and heating is introduced due to the full photoionization of He II by a local UV radiation 33 –followed and the ejection of an additional electron per helium atom, assumed to be already singlyionized. The unmodified neutral fraction is labeled as Case (I). Case (II) represents the updatedhydrogen abundance.In Figure (13), Equation (6) was solved for two values of temperature increase, T II /T I =1 . −
2. On the left panel we show the neutral hydrogen fractions in Case (II) for a range ofCase (I) inputs. The open squares (circles) represent the smallest (largest) temperature increase.The solid line stands for no difference between the two cases (slope of 1). The calculation showsthat the largest temperature increase results in a greater departure from the straight line, althoughthe degree of such departure becomes smaller at larger neutral fractions.This is shown on the right panel of Figure (13) where we plot the percentage change in thehydrogen neutral fraction. We see a fixed degree of change between χ I = 10 − − − . and a rapiddecrease at χ I & − . . For the range of temperature increase used in the calculation we can thendetermine that a total He II photoionization would lead to an adjustment of the hydrogen neutralfraction by 12-24%. This corresponds to an overestimate of χ , in the absence of He II photoheating,by ≈ − T ≃ K where the slope becomes a function of temperature itself. Theresults shown in Figure (13) are derived with a constant slope of β = 0 . T ≃ K. To explore the effects of the slope on the adjustment of the hydrogen neutral fraction,due to He II photoheating, we show on the left panel of Figure (14), the results of the calculationfor an input neutral fraction of χ I = 10 − , β = (0 . , . , .
6) and for a range of temperature ratios(0.5-2.0).The upper dashed line at T II /T I ≥ β value (0.6) and indicates thata steepening of the temperature exponent leads to an increase in the neutral hydrogen fractionadjustment. This trend is reversed at T II /T I < II photoheating effects but it was a numerical investigation of thesolution. However, one can make the argument, depending on whether the UVB drops after acertain redshift, that additional cooling terms, such as collisional ionization cooling ( He II + e − → He III + 2 e − ) and recombination radiation cooling ( He III + e − → He II + γ ), might lower thetemperature below the unperturbed value. However, such possibility is small at least up to z ≃ − χ II χ I ≥ β = 0 . − χ II χ I ≤
0. This a result of the increased electron number density due to the additionalphotoelectron ejected from the He II atom. In this narrow range of temperature increase the greaterelectron number density dominates the recombination rate, which consequently increases to shiftthe ionization balance towards more neutral hydrogen. 34 –This effect depends on the initial neutral fraction, as shown on the right panel of Figure (14).There, we plot the family of solutions of Equation (6) ( β = 0 .
5) for the range of neutral fractions(10 − ≤ χ I ≤ − ) and temperature adjustment factors (0 . ≤ T II /T I ≤ . χ I = 10 − from the left panel. The plot shows that as theneutral fraction increases the adjustment decreases for T II /T I ≥
1. As discussed in the previousparagraph, the range of T II /T I where the electron density dominates the recombination rate shiftsto the right towards lower ionization fractions. C. Radiative and Collisional Rates
In this section we document the radiative and collisional rates we used in our simulation. Forthe full range of chemical reaction rates incorporated into the cosmological hydro code Enzo pleasesee Anninos et al. (1997) and Abel et al. (1997).
C.1. Collisional Ionization and Radiative Recombination of Singly Ionized Helium
Fits are accurate within 1% in the temperature range 1 − K. He + + e − → He ++ + 2 e − Abel et al. (1997)In the following fit temperature T is in units of eV: T = T K Γ ion /n e ≡ k = exp [ − . . ln ( T ) − . ln ( T ) + 4 . ln ( T ) − . ln ( T ) +8 . × − ln ( T ) − . × − ln ( T ) +1 . × − ln ( T ) − . × − ln ( T ) ] cm s − He ++ + e − → He + + γ Cen (1992) α R ≡ k = 3 . × − T − / ( T K ) − . (1 + ( T K ) . ) − cm s − Temperature is expressed in K.
C.2. Photoionization cross section of Hydrogen and Helium II H + γ → H + + e − He + + γ → He ++ + e −
35 –The cross section for both rates are expressed through the same functional form because singlyionized helium is a hydrogen like atom (from Osterbrock 1974). σ = AZ ( νν Z ) exp [4 − arctanǫ ) /ǫ ]1 − exp ( − π/ǫ ) A = 6 . × − cm , ǫ = ( νν Z − / and ν Z = 13 . × Z ( Z = 1 , ν ≡ ν ( Z =1) , ν ≡ ν HeI ( Z =2) and ν ≡ ν HeII ( Z =2) .For numerical convenience we use an approximation to the full He II cross section in the form σ HeII = A (= σ oHeII )( hνhν ) − , where hν ≡ . ≈
3% (Haehnelt & Abel 1999).
REFERENCES
Abel, T., & Haehnelt, M. 1999, ApJ, 520L, 13Aguirre, A., Schaye, J. & Theus, T. 2002, ApJ, 576, 1Anderson, S. F., Hogan, C. J., Williams, B. F. & Carswell, R. F. 1999, AJ, 117, 56Anninos, P., Zhang, Y., Abel, T. & Norman, M. L. 1997, NewA, 2, 209Bertschinger, E. & Gelb, J. M. 1991, ComPh, 5, 164Bolton, J., Meiksin, A. & White, M. 2004, MNRAS, 348L, 43Bryan, G. L. & Norman, M. L. 1997, in IMA workshop in Mathematics and Applications, Vol.117, p. 165 (Kluwer Academic)Bryan, G. L., Machacek, M., Anninos, P. & Norman, M. L. 1999, ApJ, 517, 13Bryan, G. L. & Machacek, M. E. 2000, ApJ, 534, 57Carswell, R. F., Webb, J. K., Baldwin, J. A. & Atwood, B. 1987, ApJ, 319, 709Carswell, R. F. 1989, Nature, 341, 692Cen, R. 1992, ApJS, 78, 341Cen, R., Miralda-Escud´e, J., Ostriker, J. P., & Rauch, M. 1994, ApJ, 437, L9Ciardi, B., Ferrara, A. & White, S. D. M. 2003, MNRAS, 344, 7Dav´e, R., Hernquist, L., Weinberg, D. H. & Katz, N. 1997, ApJ, 477, 21 36 –Davidsen, A. F., Kriss, G. A. & Zheng, W. 1996 , Nature, 380, 47Efstathiou, G. 1992, MNRAS, 256P, 43Eisenstein, D. & Hut 1998, ApJ, 498, 137Eisenstein, D. & Hu, W. 1999, ApJ, 511, 5Fechner, C., Reimers, D., Kriss, G. A. et al. 2006, A&A, 455, 91Giroux, M. L. & Shapiro, P. R. 1996, ApJS, 102, 191Gnedin, N. Y., & Abel, T. 2001, NewA, 6, 437Haardt, F. & Madau, P. 1996, ApJ, 461, 20Haardt, F. & Madau, P. 2001, in XXI Moriond Astrophysics Meeting, Galaxy Clusters and theHigh Redshift UniverseHaehnelt, M. G., & Steinmetz, M. 1998, MNRAS, 298, 21Heap, S. R. et al. 2000, ApJ, 534, 69Hernquist, L., Katz, N., Weinberg, D. H. & Miralda-Escud´e, J. 1996, ApJ, 457L, 51Hu, E. M., Cowie, L. L., McMahon, R. G., et al. 2002, ApJ, 568, L75Hui, L. & Gnedin, N. Y. 1997, MNRAS, 292, 27Jacobsen, P., et al. 1994, Nature, 370, 35Jacobsen, P., Jansen, R. A., Wagner, S. & Reimers, D. 2003, A&A, 397, 891Jena, T., Norman, M., L., Tytler, D., Kirkman D., et al. 2005, MNRAS, 361, 70Katz, N., Weinberg, D. H. & Hernquist, L. 1996, ApJS, 105, 19Kirkman, D. & Tytler, D. 1997, ApJ, 484, 672Kriss, G. A., et al. 2001, Science, 293, 1112McDonald, P., Miralda-Escud´e, J., Rauch, M., Sargent, W. L. W. et al. 2000, ApJ, 543, 1Mihalas, D. & Mihalas, B. W., 1984, ”Foundation of Radiation Hydro dynamics”, Oxford UniversityPress, New YorkMiralda-Escud´e & Ostriker 1990Miralda-Escud´e, J. & Rees, M. J. 1994, 266, 343 37 –Miralda-Escud´e, J., Cen, R., Ostriker, J. P, & Rauch, M. 1996, ApJ, 471, 582Norman, M. L., Paschos, P., & Abel, T. 1998, H in the Early Universe, ed. F.Palla, E. Corbelli,& D. Galli (Mem. Soc. Astron. Italiana), 69, 271Osterbrock, D. E., 1989, “Astrophysics of Gaseous Nebulae and Active Galactic Nuclei”, UniversityScience Books, Mill Valley, CAPaschos, P., Norman, M. L., 2005, ApJ, 631, 59Pei, Y. C. 1995, ApJ, 438, 623Petitjean, P., Webb, J. K., Rauch, M., Carswell, R. F. & Lanzetta, K. 1993, MNRAS, 262, 499Rauch, M., Miralda-Escud´e, J., Sargent, W. L. W., Barlow, T. A., Weinberg, D. H., Hernquist, L.,Katz, N., Cen, R. & Ostriker, J., P., 1997, ApJ, 489, 7Razoumov, A., Norman, M. L., Abel, T. & Scott, D. 2002, ApJ, 572, 695Reimers, D., Koller, S., Wisotski, L., Groote, D., Rodriguez-Pascual, P. & Wamsteker, W. 1997,A & A, 327, 890Reimers, D., Fechner, C., Hagen, H.-J. et al. 2005, A&A, 442, 63Schaye, J., Theuns, T., Leonard, A. & Efstathiou, G. 1999, MNRAS, 310, 57Schaye, J., Theuns, T., Rauch, M., Efstathiou, G. & Sargent, W. L. W. 2000, MNRAS, 318, 817Sokasian, A., Abel, T. & Hernquist, L. 2002, MNRAS, 332, 601Theuns, T., Leonard, A. & Efstathiou, G. 1998, MNRAS, 297L, 49TTheuns, T., Schaye, J. & Haehnelt, M., G. 2000, MNRAS, 315, 600TTytler, D., Kirkman, D., O’Meara, J. M., Suzuki, N., Orin, A., Lubin, D., Paschos, P., Jena, T.,Lin, W., Norman, M. L. & Meiksin, A., 2004, ApJ, 617, 1Wadsley, J. W. & Bond, J. R. 1997, ASPC, 123, 332WWhalen, D., Abel, T. & Norman, M. L. 2004, ApJ, 610, 14WZaldarriaga, M., Seljak, U. & Hui, L. 2001, ApJ, 551, 48ZZaldarriaga, M., Hui, L. & Tegmark, M. 2001, ApJ, 557, 519ZZhang, Y., Anninos, P. & Norman, M. L. 1995, ApJ, 453L, 57ZZhang, Y., Anninos, P., Norman, M. L. & Meiksin, A. 1997, ApJ, 485, 496Zhang, Y., Meiksin, A., Anninos, P. & Norman, M. L. 1998, ApJ, 495, 63Z 38 –Zheng, W., Chiu, K., Anderson, S. F., et al. 2004a, AJ, 127, 656Zheng, W., Kriss, G. et al. 2004, ApJ, 605, 631Z This preprint was prepared with the AAS L A TEX macros v5.2.
39 –Fig. 1.— Left Panel: Volumetric rendering of the logarithm of baryon overdensity at z=2.6. RightPanel: Evolution in the number of QSO sources versus redshift in this simulation 40 –Fig. 2.— Volumetric renderings of the He
III distribution mapped by the Log n HeIII at z=3 andz=2.6. The ionization cutoff in this visualization is set to ψ ≡ n HeIII n He = 10 − (dark regions). 41 –Fig. 3.— Left Panel: Slice of the He III (darkened regions) and He II
3D mass distribution at z=2.6.At that redshift, the volume filling fraction of the He
III is V F F ≃ .
90. Right Panel: RedshiftEvolution of the VFF for χ HeIII ≥ − . The error bars result from the uncertainty in the locationof the cumulative I-front due to finite grid cell size. 42 –Fig. 4.— Left Panel: Redshift evolution of the volume averaged mass fraction of singly helium.Overplotted is the corresponding neutral hydrogen mass fraction. Right Panel: The scatter plotbetween He II fraction and local gas overdensity at z=2.5 yields the mean and median value of Logf
HeII per logarithmic bin of gas overdensity. Solid (dashed) curves refer to the inhomogeneous(homogeneous) case. The median profiles correspond to the curves that peak at
Logδ ≃ δ = 1 −
10. Right Panel: The median temperature distribution in each logarithmicoverdensity interval at z=2.5. It shows that the processed (solid) and unprocessed (dashed) cases.The curves show that He II reionization primarily increases the temperature of the IGM in low tomildly overdense range. 44 –Fig. 6.— Left Panel: From temperature-density plot in Figure (5) we compute the ratio of thepostprocessed result over the input simulation (dashed) and the optically thin simulation (dashed-dot). Right Panel: We plot the distribution of the logarithmic slope dlnTdlnρ = γ − II and H I Ly α spectra between z = 2 . − .
4, at ¯ z = 2 .
5, along a randomly castedsightline through the computational volume. The upper panels refer to the HI Ly α transmission.The lower panels refer to the HeII Ly α transmission. The horizontal axis is converted to observedwavelength λ = λ o (1 + z ), where λ o is the restframe wavelength of the resonant Ly α scatter.Left Panels: Uniform UVB ionizes and heats the IGM in the optically thin approximation. HeIIionization and heating is included in this calculation. Right Panels: The uniform UVB only ionizesH I to H II and He I to He II . Photoheating therefore due to the uniform UVB only refers to theabove processes. The non-uniform point source radiative input from QSO type sources is used tophotoionize He II to He III and to calculate the subsequent photoheating. Although the H I spectraseem similar between the two sets of calculations, ignoring the feedback on H I abundance by theHe II ionization results in underestimating the mean H I Ly α transmission, compared to the fullsimulation, by ≈ z = 2 .
5. Ateach point along the path of the sightline, we compute the distance to the closest point source. Asthe trajectory crosses the computational volume, the closest source changes depending on position.On the lower left part of the figure, we plot this proximity distance for a sightline that comes theclosest to a point source, about one grid resolution element. On the top left panel, we show theHe II transmitted flux for that sightline. Crosses mark the locations along the path that are closestto the point source. The alternating proximity sources are indicated by the alternating parabolicprofiles. On the right panel, we show the probability distribution of the minimum distance to thesources, which we define as the impact parameter, from all random sightlines. The asymmetry inthe distribution is small, which is due to the near isotropic distribution of high density peaks atthe scale of the simulated volume. 47 –Fig. 9.— On the left panel, we show the scatter plot between the average impact parameter foreach sightline in respect to the distribution of the ionizing point sources versus the sightline-specificmean flux. The red line is power law fit to the data with a slope of s = − . ± .
04. The blue solid(dashed) line shows the average (2 σ ) He II Ly α transmission from all sightlines at ¯ z = 2 .
5. Thesmall negative correlation r = − .
24 suggests that the mean flux of a sightline is largely insensitiveto the average proximity in respect to the distribution of ionizing sources along it’s path. On theright panel, we show the locally absorbed flux (y-axis) at the point of closest proximity to theionizing source shown on the x-axis. Individual points correspond to the absorbed flux at the pixelof closest distance. The red (blue) histogram is the median (mean) absorbed flux per constantlogarithmic bins. Solid and dashed blue horizontal lines are the level of the mean absorption inour sightline sample and the 2 σ standard deviation respectively. The solid line is a Fermi-Diracfunction that fits the median (red) histogram for impact parameters less that 10 Mpc comoving.The half point of the Fermi-Dirac function is at d o = 3 . . . σ level of the meanabsorption, show a steep sensitivity in the absorbed flux due to the quasar proximity. 48 –Fig. 10.— We compute the mean optical depth and effective optical depth at selected meanredshift points, from ¯ z = 6 − . δz = 0 .
25 about themean redshift. Left Panel: Redshift evolution of the mean pixel Ly α He II optical depth averagedover all lines of sight. The error bars are 1 σ standard deviation to the sightline average opticaldepth. Right Panel: Redshift evolution of the Ly α He II effective optical depth. The error barsin the simulation results are sightline-to-sightline errors estimated from the standard deviation tothe average transmission from sightline-specific mean flux. Overplotted are sightline measurementsalong HE2347-4342 (red triangles), HS1700+6416 (blue circles), HS1157+3143 (blue diamond),PKS1935-692 (blue cross), Q0302-003 (green square) and SDSSJ2346-0016 (orange star). Errorbars to these points refer to error estimates of the pixel flux along the sightline specific to thequasar. The HE2347-4342 data from Zheng et al. (2004b) have been resampled to a larger binsize to make it comparable to our data. The redshift profiles of the optical depth point statisticssuggest a smooth evolution of the He II reionization epoch. 49 –Fig. 11.— Left Panel: Column Density Distributions, f ( N ) = d NdzdN X of H I and He II absorbersbetween z = 2 . z = 2 .
4. Circle and square points correspond to the H I and He II columndensity distribution respectively. Solid (dashed) straight lines (in log-log) show the power law fitto the H I (He II ) distribution between column densities N X = 10 . − . cm − . The errorbars show the 1 σ standard deviation to the number of lines per logarithmic column density binequal to 0.25. Right Panel: Broadening width distributions at ¯ z = 2 . II and H I . Thedistributions are shown as the fraction of absorption lines per km/s in bins of 2 km/s. Black andred colors correspond to H I and He II distributions respectively. The dashed lines show a referenceself-consistent calculation in the optical thin limit using an ionizing uniform UVB. The peaks ofnon-uniform calculations are offset by ≃ .
25 km/s to the right of the uniform results. 50 –Fig. 12.— Redshift evolution of the η parameter in this simulation. Results are shown as < logη > per redshift interval and are obtained as the line column density ratio of He II over H I at low z( z <
3) and estimated from the ratio of optical depths at higher redshifts. The filled triangle showsthe average estimate between z = 2 . − . z = 2 .
5. 51 –Fig. 13.— Left Panel:
Logχ I (input) vs. Logχ II (output) for T II /T I = 1 . T II /T I =2 . Logχ I . 52 –Fig. 14.— Left Panel: Percentage change vs. T II /T I for three values of the exponent in α (1) R ∝ T − β .Solid line: β = 0 .
5. Upper dashed line at T II /T I ≥ β = 0 .
6. Lower dashed line at T II /T I ≥ β = 0 .
4. Right Panel: Percentage change vs. T II /T I for 10 − ≤ χ I ≤ −1