Mon.Not.R.Astron.Soc.001,1–19(2010) Printed4January2012 (MNLATEXstylefilev2.2) Radiative transfer of energetic photons: X-rays and helium 2 ionization in C -RAY 2 Martina M. Friedrich1 ⋆, Garrelt Mellema1, Ilian T. Iliev2 and Paul R. Shapiro3 1 1DepartmentofAstronomy&OskarKleinCentre,AlbaNova,StockholmUniversity,SE-10691Stockholm,Sweden 0 2AstronomyCentre,DepartmentofPhysics&Astronomy,PevenseyIIBuilding,UniversityofSussex,Falmer,BrightonBN19QH 2 3DepartmentofAstronomyandtheTexasCosmologyCenter,TheUniversityofTexasatAustin,TX78712,USA n a J 3 4January2012 ] O ABSTRACT C We presentan extensionto the short-characteristicray-tracingand non-equilibriumphoton- h. ionizationcodeC2-RAY.Thenewversionincludestheeffectsofheliumandimprovedmulti- p frequency heating. The motivation for this work is to be able to deal with harder ionizing - spectra,suchasforexamplefromquasar-likesourcesduringcosmicreionization.Wereview o the basic algorithmic ingredients of C2-RAY before describing the changes implemented, r t which include a treatment of the full on the spot (OTS) approximation, secondary ioniza- s tion,andmulti-frequencyphoto-ionizationandheating.Weperformedaseriesoftestsagainst a [ equilibriumsolutionsfromCLOUDYaswellascomparisonstothehydrogenonlysolutions by C2-RAY in the extensive code comparison in Ilievetal. (2006). We show that the full, 1 coupled OTS approximation is more accurate than the simplified, uncoupled one. We find v thatalso with heliumand a multi-frequencyset up, long timesteps(up to ∼ 10%of the re- 2 combinationtime) still give accurate results for the ionization fractions.On the other hand, 0 accurateresultsforthetemperaturesetstrongconstrainsonthetimestep.Thedetailsofthese 6 constraintsdependhoweverontheopticaldepthofthe cells. We use thenew versionofthe 0 . codetoconfirmthattheassumptionmadeinmanyreionizationsimulations,namelythathe- 1 lium is singly ionized everywherewere hydrogenis, is indeed valid when the sourceshave 0 stellar-likespectra. 2 1 Keywords: methods:numerical–radiativetransfer–galaxies:intergalacticmedium–HII v: regions i X r a 1 INTRODUCTION warmandionized.Thisprocessisknownasreionizationandwas theUniverse’slastglobalphasetransition(seeforexamplethere- Photo-ionizationisoneofthemajorradiativefeedback processes viewinBarkana2006). inastrophysics. Theextremeultraviolet(EUV)photonsproduced by massive stars in star formation regions are capable of heating Traditionallyphoto-ionizationcodesconcentratedonsolving the gas to temperatures around 104 K and the heated ions pro- equilibrium situations, as for example CLOUDY (Ferlandetal. ducecopiousamountsofcollisionallyexcitedlineradiation,lead- 1998), MAPPINGS (Sutherland&Dopita 1993) and the three- ing tothewell-known and sometimes spectacular images of H II dimensional code MOCASSIN (Ercolanoetal. 2008) do. The regions(suchasforexampletheOrionnebula,O’Dell2001).Ac- mainaimofthesecodesistoaccuratelycalculatelinestrengthsfor creting black holes and neutron stars also produce ionizing radi- comparison tospectroscopic observations. However, sincethein- ation and, depending on their mass, can ionize smaller or larger creaseinpressurecandrivepowerfulflowsinthegas,thereisalso regionsaround themselves. Ongalacticscales, theemissionfrom aneedtocouple photo-ionization calculationstohydrodynamics. supermassive central black holes are observed to produce ioniza- Thisnecessitatesasimpler versionof theradiativetransfer, since tionconesstretchingintothegalaxy’simmediateenvironment(see ithastobecalculatedinstepwiththehydrodynamicsandfordy- e.g.Pogge1989).Evenatthelargestscalesphoto-ionizationisim- namiccalculationstheindividuallinestrengthsarelessinteresting. portant. Some time before redshift 6, ionizing radiation escaped Thehistoryofthesetypesofcalculationsgoesbackquitefar,see from the first generations of galaxies and percolated through the forexampleYorke(1986)foranoverviewoftheworkdonebefore intergalacticmedium(IGM),changingitfromcoldandneutralto 1986. Forsomeoftheapplicationsmentionedabove,alowerdimen- sional approach (one or two dimensional) is sufficient. However, ⋆ e–mail:[email protected] other applications, such ascosmic reionization, require thetrans- 2 Friedrichet al. port tobeperformed inthefullthreedimensions. Because ofthe anoverviewofthebasicalgorithmicideasbehindC2-RAYtothen higherdimensionalityhereaswellasimplerversionoftheradiation proceedinSection3withanoverviewofhowweextendeditscapa- and photo-ionization physics is typically implemented, although bilitiestohandleharderphotons.Section4containsthedescription thelevelofsophisticationvariesbetweenmethods(seee.g.thefirst ofaseriesofoneandthree-dimensionaltestsforthenewmethod, generation methods of Razoumov&Scott 1999; Sokasianetal. alsoevaluatingtheeffectsofthevariousnewelementssuchassec- 2001;Nakamotoetal.2001;Ciardietal.2001).Insomecasesthe ondary ionizationsand thecoupled on-the-spot approximation. A couplingtothedynamicswasalsodoneforcosmologicalapplica- seriesof appendices describe several important elements in more tions(e.g.,Ricottietal.2002;Tracetal.2008;Wise&Abel2008). detail. One of the simplifications which is often employed is to consider hydrogen as the only element being photo-ionized. Since close to 10% of the gas is helium, this approxima- tion is crude. However, as shown for example in figure 2.4 of Osterbrock&Ferland(2006),fortypicalO-starspectratheioniza- tionofheliumfollowslargelythatofhydrogen.So,ifoneisonly 2 REMINDEROFBASICSTEPSOFTHEORIGINAL interestedinthestructureoftheionizedregions,assumingthathe- C2-RAYALGORITHM liumfollowshydrogenisnotabadapproximation. However,asonemovestoharderspectra,thisassumptionbe- TheC2-RAYmethodwasdevelopedtobeatime-dependentphoto- comeslessandlessvalid.Notonlydoesoneneedtotakeintoac- ionization algorithm that could be efficiently combined with a count that helium can become doubly ionized, also the fact that hydrodynamics calculation, and not impose impractically short thecrosssectionofthetwotypesofheliumstartcontributingsub- timesteps and small cell sizes on the latter. This is achieved by stantially to the opacity of the gas becomes an issue. This prob- assuming that the ionization evolution of individual cells follows lem becomes especially important when the ionizing spectrum is anexponential decay tothe equilibriumsolution and that a time- powerlaw-like, as one expects from hot accretion disks around averagedvalueoftheopticaldepthcanbeusedtodescribetheef- blackholes.Specificallyforthecaseofcosmicreionization,where fect of a cell on the radiative transfer during the entire timestep. there is a possible contribution from powerlaw-like sources such Thisapproachisabletocorrectlytracktheprogressofionization asquasarsandmini-quasars,anyproperstudyoftheircontribution frontsover many cellsduring onetimestep. Inaddition, optically shouldconsiderbothhydrogenandhelium. thickcellsaredealtwithbydefiningthephoto-ionizationratesuch This has motivated us to extend the capabilities of the code thatitisconsistent withthenumberofphotonsabsorbedinsidea C2-RAY to include the effects of helium. C2-RAY is a photon- cell.C2-RAYwasdescribedandtestedindetailinM+06.Herewe conserving radiative transfer code that uses short characteristic summarize some of the ideas in order to define our notation and ray tracing and is described in detail in Mellemaetal. (2006b) provideanintroductiontotheextensionsdescribedinSect.3. (hereafter M+06) . It has been used extensively for reionization The evolution of the ionized hydrogen fraction derives from simulations(e.g.Mellemaetal.2006c;Ilievetal.2007,2008, to thesetofchemicalevolutionequations: name a few). It was also combined with a hydrodynamics code (CAPREOLE-C2, Mellemaetal.2006a),andanidealmagnetohy- d xHI = −ΓHI+CHIne αBHIIne xHI , drodynamicscode(PHAB-C2, deColle&Raga2006;Arthuretal. dt(cid:18) xHII (cid:19) (cid:18) ΓHI+CHIne −αBHIIne (cid:19) (cid:18) xHII (cid:19) 2011), to investigate galactic H II regions. Furthermore, it was (1) tested against other non-equilibrium radiative transfer codes in Ilievetal.(2006)(hereafterI+06)andinconjunctionwithagrid- wherexHIistheneutralhydrogenfraction,xHIIistheionizedhy- basedhydrodynamiccode(fordetails,seeMellemaetal.2006a)in drogen fraction, ne is the electron density, ΓHI is the hydrogen a second comparison project, including gas dynamics (Ilievetal. photo-ionizationrate(seebelow),CHIisthecollisionalionization 2009). rateandαBHIIistherecombinationrate.SincexHI+xHII=1,we Adding helium implies introducing another source of obtain frequency-dependent opacity, thus making a multi-frequency ap- d proach inevitable. In addition the on-the-spot approximation be- dtxHII=−(ΓHI+CHIne+αBHIIne) xHII+(ΓHI+CHIne) . comes morecomplicated asone hastotake intoaccount how re- AH gH combination photons fromheliumaffect thehydrogen ionization. | {z } | {z }(2) Asonemovestohigherphotonenergies,oneshouldalsotakeinto accountthesecondaryionizationscausedbythesuperthermalelec- Here we introduce the notation, AH and gH, which allows us to trons produced when high energy photons ionize the atoms and writetheequationinthegeneralvectorform,usefullateron, ions. Including EUV and soft X-ray (SX) photons therefore is a d non-trivialextensionofthephoto-ionizationcalculationswhichwe x=Ax+g. (3) dt describeandtestinthispaper. Asiswellknown,thegeneralsolutionx(t)toasetofequationsof In terms of the physical processes included the method we thistypeisthesumofthesolutiontothehomogeneouscase,x (t), present here is similar to a number of others published in re- h whereg=,andaparticularsolutionx : centyears,butdiffersinthealgorithmsused.Thecodes CRASH p (Masellietal.2009)andLICORICE(Baeketal.2010)useMonte n Carlo techniques. The codes SPHRAY (Altayetal. 2008) and x(t)=xh(t)+xp with xh= cixietλi. (4) TRAPHIC (Pawlik&Schaye 2011) implement ray tracing on Xi=1 particledatawhereasRADAMESH(Cantalupo&Porciani2011) Here,nistherankof A,i.e.thenumber of coupled equations (1 usesanadaptivemeshapproach. in the case of hydrogen), λi are the eigenvalues of A, xi are the The lay-out of the paper is asfollows. In Section2 we give corresponding eigenvectors and ci are coefficients which can be Radiativetransferofenergeticphotons:X-raysandheliumionizationin C2-RAY 3 calculatedfromtheboundaryconditionx(t=0)=x : 0 n x0 = cixi+xp (5) i=1 X Insubsequent timesteps,x isthestateattheendoftheprevious 0 timestep.Inthecaseofaconstantg,theparticularsolutioncanbe theequilibriumsolutiongivenby Axp+g=. Forthesimplest,hydrogenonlycase,thecoefficientsarethus λH = −(ΓHI+CHIne+αBHIIne) xH = 1 xH = ΓHI+CHIne p ΓHI+CHIne+αBHIIne cH = x −xH (6) 0 p Here,weaddedthesuperscriptHforhydrogenandweskippedthe subscript1sinceforhydrogen,n=1. The photon-conserving photo-ionization rate Γ in each cell usedinAiscalculatedusing Γ = ∞ Lνe−hτνi1−e−h∆τνidν, (7) HI Zνth hν hnHIiVshell iFnigMu+re061,.fiItgeurraetio2nschemeofC2-RAYforasinglecell,conceptionallyas whereh∆τνiisthetimeaveraged opticaldepthover thecelland V isthevolumeoftheshellthecellbelongsto.Thisquantity shell can be calculated from the time evolution of the neutral fraction cies higher than the ionization threshold of HeII, the ionization (Eq.4).Bysolvingthesetwoequations(4and7)inalternatingor- cross-sectionsofHeIandHeIIareroughlyanorderofmagnitude deroneiteratestoconvergenceasillustratedinFig.1.Thisiteration larger than the HI ionization cross-section. Therefore, neglecting alsoinvolvestheelectrondensityne whichiscalculatedfromthe heliuminthecaseofsourceswithhardspectrawillunderestimate timeaveragedionizedfraction.Thetimeaveragedopticaldepthto the optical depth substantially. We therefore have to add helium the cell hτνi iscalculated by short characteristic ray-tracing over chemistrytoC2-RAY. thesolutionsfoundforcellslyingnearertothesource.Thismakes Inordertoincludehelium,boththechemicalevolutionequa- thealgorithmcausal. tion,Eq.(2),andthecalculationoftheionizationrate,Eq.(7),have Inthecaseofmultiplesources,theiterationasshowninFig.1 tobechanged.Additionally,asshownbelow,theiterationscheme is split up in two parts: The first part, including step three (find- fromFig.1hastobemodified.Wedescribeeachofthesechanges ingtheionizationand heatingratesineach cell)isdonefor each here. source, looping through theentirecomputational gridusing short characteristicraytracing.Foreachcell,theratesfromallsources areadded.Thesetotalratesareusedintheremainingtwostepsin 3.1 Chemicalevolutionequation theiteration.SeealsoM+06foraflowchart andadescriptionof Theprocedurefor addingheliumphoto-ionization toour calcula- theimplementationhowtoloopthroughthesourcelist. tionsisbyitselfrelativelystraightforwardasthebasicalgorithmic NotehoweverthattheflowchartinM+06forthesinglesource idea described in Sect. 2 provides the framework for this. How- loop(figure4)incorrectlyincludesthelasttwostepsofthesingle ever,acomplicatingfactoristhepresenceofionizingrecombina- cell loop. The electron density and the time averaged ionization tionphotonssincetheycoupletherateequationsofhydrogenand fractionsareinfactnotupdatedinthesourceloopbutareupdated helium. Here we present two approaches for dealing with these, firstafterthephoto-ionizationratesfromallsourcesaresummedto wherethefirstoneislessaccurate,butsimplerandmoresimilarto aglobalphoto-ionizationrate. whatotherauthorshaveused.Wecomparetheresultsofthesetwo approachesinSect.4. 3 EXTENDINGC2-RAY 3.1.1 Simplerecombination:nocouplingofspecies The original C2-RAY methodology works well for soft, stellar spectra. In this case, the IGM can be considered to be hydrogen Whenionsrecombine,photonsareemitted.Incaseofrecombina- only since there are not many photons capable of ionizing He II tionofhydrogenions,onlyrecombinationstothegroundstatere- andheliumcanbeassumedtobesinglyionizedeverywherewhere sultinphotonsenergeticenoughtoionizehydrogen.Ifoneassumes hydrogen isionized. InSect. 4.3.1weshow that themorphology thesephotonstoionizeimmediatelyanotherhydrogenatomclose oftheionizationfractionfieldinacosmologicalreionizationsim- by, thisiscalled theon the spot (OTS)approximation for hydro- ulationwithstellarsources(only),isindeedhardlyaffectedbythe gen (e.g. Osterbrock&Ferland 2006). In this approximation, the inclusionofhelium.However,whenthespectrumhasasignificant recombinationcoefficienttoallstatesofhydrogen,αAisreplaced amountofSXphotons,heliumcontributessignificantlytotheop- bytherecombinationcoefficienttoallstatesbutthegroundstate, ticaldepthandamulti-frequencyapproachisrequired:atfrequen- αB.Foramixofhydrogenandhelium,theOTSapproximationis 4 Friedrichet al. morecomplicatedasheliumrecombinationphotonscanionizeboth HeII,HeIandHI,dependingontherelativeopticaldepthoverthe hydrogenandhelium.However,asafirststep,weassumethatpho- cellinquestionatνHeII.Ofthoserecombinationphotons,afraction th tons fromrecombinations to the ground statecan only ionize the yagoesintoHeIIionization,afractionybgoesintoHeIionization 2 2 samespeciesfromwhichtheyoriginateandthatinrecombinations andafraction1−ya−ybgoesintoHIionization.Here,thefrac- 2 2 tootherstatesthanthegroundstatenoionizingphotonsareemit- tionsy,ya,andybaredependentontherelativeopticaldepthofthe 2 2 ted.Thatmeans,weusetheαB recombination-coefficientsforall speciesatthethresholdfrequenciesofHeIandHeII,asdescribed species.Inthefollowingwerefertothisasthe“uncoupledon-the- below. spotapproximation”(U-OTS).Inthisapproximationhydrogenand RecombinationsofHeIIItootherstatesthanthegroundstate helium can be treated separately. For helium, theset of chemical contributetoionizationofHIandHeI.Thoserecombinationslead evolutionequations(inanalogytoEq.1forhydrogen)is: eithertotwophotonemission(inafractionv ofthecases,where vistemperaturedependent,Hummer&Seaton1964),ofwhichon x d HeI average a fraction l is energetic enough to ionize HI and a frac- x = dt HeII  tion m is energetic enough to ionize HeI; therefore, a fraction x HeIII vw = v(l −m+ my) goes into HI ionization and a fraction   −UHeI neαHeII 0 xHeI v(m(1−y)) goes into HeI ionization. The remaining fraction, UHeI −neαHeII−UHeII neαHeIII xHeII 1−v,leadstoemissionofaHeLymanαphoton.Thosephotons     0 UHeII −neαHeIII xHeIII areabsorbedbyanyspeciestoafractionf.Bylettingescapesome    (8) of the helium Ly α photons (f 6= 1) we lose them since we do not include those at larger distances from the source. Of the ab- wherex istheneutralheliumfraction,x isthesinglyion- HeI HeII sorbedHeLymanαphotons, afractionz goesintoHIionization izedheliumfractionandx isthedoublyionizedheliumfrac- HeIII andtheremainingfraction1−zintoHeI.Additionally,theBalmer tion.Wegroupedtheionizing‘uprates’intooneterm continuumemissionphotons(α2 )canionizehydrogen.Table1 HeIII Un ≡Γn+neCn. (9) summarizesthephotonemittingrecombinationprocessesincluded in the on the spot treatment. The numerical parameters used are ThesubscriptsonC,αandΓindicateonwhichspeciestheyact,so listedinTable2. HeIIα⇌HeIIIHeIII.UsingthefactthatxHeI+xHeII+xHeIII =1, We can now introduce these terms in the general equation, ΓHeII theequivalentequationtoEq.(2)canbewrittenas Eq.(3): d xHeII =AHe· xHeII +gHe, (10) xHII UHI dt(cid:18) xHeIII (cid:19) (cid:18) xHeIII (cid:19) x= xxHeII  g= UH0eI  (12) HeIII wherethevectorgHeandthematrixAHehavethefollowingforms     and U gHe= HeI 0 (cid:18) (cid:19) nHene× AHe =(cid:18) −neαHeII−UHeUIIH−eIIUHeI ne−αnHeeIαIIH−eIIUI HeI (cid:19)(11) +−α(UBHHnIe+) (yα1HneHInIeH+npeα×BHeII) α2H((efIIzI(+1(−1v−)ny+H2av−wy)2bα)BHαe1HIIeII+II))    A= −UHeI+  Thesolution of this set of linear differential equations for the U-  −UHeII−UHeI− ne(α1HeIIIy2b+  OTScasecanbefoundinAltayetal.(2008)orwiththehereintro-  0 (αAHeII−(1−y)α1HeII)ne (αAHeIII−α1HeIIIy2a)+  ducednotationinAppendixA.  αBHeIII(f(1−z)(1−v)+(l−w)v))    0 UHeII (−ne(αAHeIII−y2aα1HeIII))    3.1.2 Onthespotapproximation:couplingofspecies (13) Inreality,recombinationphotonsfromheliumcanionizeeitherhy- ThedensityfractionsnHe/nH inA12 andA13 areduetothe drogenorhelium,introducingtheneedtocoupletherateequations factthattheequationsarewrittenintermsoffractions,notinterms for the two elements. Thisis the proper OTS approximation. Ta- ofdensities.Thecompletesolutionforthissetofequationsispre- ble 1 gives an overview of the recombination processes affecting sentedinAppendixB. hydrogenandheliumfractionsintheOTSapproximation. To implement the OTS approximation we follow 3.2 Extendingtheiterationmechanism Osterbrock&Ferland (2006) for dealing with the recombina- tions of HeII to HeI and Flower&Perinotto (1980) for those of In Sect. 3.1.2, we treat the set of chemical evolution equations HeIII to HeII (except for recombinations to the ground-state). dx/dt = f(x) as a set of linear differential equations, Eq. (3), Mostlywealsousetheirnotation. anduseiterationstoobtainthecorrectelectrondensityandioniza- Forhydrogen,wetaketheαB-recombinationcoefficientαB . tionrates.However,furtherdependencesonthedifferentionization HII ForHeII,thephotonsfromrecombinationstothegroundstateare fractionsarehiddenintheparametersy,z,w,yaandyb.Wefound 2 2 distributed between helium and hydrogen depending on the frac- thatthishiddenextranon-linearitycancomplicatetheconvergence tionofopticaldepthattheheliumionizationthresholdfrequency, oftheiteration.Wethereforeextendedouriterationschemetomake νHeI:afractiony goesintohydrogenionization,afraction1−y itmorerobust.Theextensionconsistsofupdating theparameters th goesintoheliumionization.Thephotonsfromstatesotherthanthe y,z,w,ya andyb afterthechemicalevolutionequationhasbeen 2 2 ground state contribute with a fraction p to hydrogen ionization. solvedandthensolvingthechemicalevolutionequationforasec- Similarly,for HeIIItherecombinations totheground stateionize ondtimewiththisnew setofparameters. Wethentakethemean Radiativetransferofenergeticphotons:X-raysandheliumionizationin C2-RAY 5 Table1.SummaryofrecombinationprocessesincludedintheOTStreatment HII α−1H→IIHI groundstaterecomb ⇒ HI→HIIionization HeII α−BH→eIIHeI deexciationsfrom ⇒p HI→HIIionization recombinationston≥2 HeII α−1H→eIIHeI groundstaterecomb ⇒y HI→HIIionization 1−y ⇒ HeI→HeIIionization HeIII α−1H→eIIIHeII groundstaterecomb ⇒y2a HeII→HeIIIionization ⇒y2b HeI→HeIIionization 1−y⇒2a−y2b HI→HIIionization HeIII α2H→eIIIHeII recombton=2 ⇒ HeIIBalmercontinuum ⇒ HI→HIIionization HeIII αBH→eIIIHeII deexcitationsfrom ⇒v 2photondecay ⇒w HI→HIIionization m(1−y) recombinationston≥2 ⇒ HeI→HeIIionization 1−v f(z−1) ⇒ HeIILyαphoton ⇒ HeI→HeIIionization fz ⇒ HI→HIIionization Table2.OverviewofthenumericalparametersusedintheOTSapproximation p = 0.96or0.66dependingonne(Osterbrock&Ferland2006),inourcasesofinterest(cosmologicalsimulations)always0.96 y = τH/(τH+τHeI)atνtHheI(ionizationthresholdofHeI) y2a = τHeII/(τH+τHeI+τHeII)atνtHheII(ionizationthresholdofHeII) y2b = τHeI/(τH+τHeI+τHeII)atνtHheII(ionizationthresholdofHeII) z = τH/(τH+τHeI)athν =40.8eV(HeILyα) f = 1to0.1(”escape”fractionofLyαphotons)dependingontheneutralfraction v = temperaturedependentcoefficient(Hummer&Seaton1964) w = (l−m)+my l = 1.425,fractionofphotonsfrom2-photondecay,energeticenoughtoionizehydrogen(Flower&Perinotto1980) m = 0.737,fractionofphotonsfrom2-photondecay,energeticenoughtoionizeneutralhelium(Flower&Perinotto1980) ofthesefirstandsecondsolutionsforthetime-averagedionization inthebins,itisusefultofurthersubdividebin2andbin3,hence- fractionsandusethismeantocalculatethephoto-ionizationrates. forth referred to as (frequency) sub-bins. To refer to a particular These extrasteps mean that the iterationscheme now consists of sub-bin,weusethefollowingnotation:n.mreferstosub-binmof seveninsteadoffoursteps,seeFig.2. binn(n=2,3).InAppendixCweshowourpower-lawfitstothe cross-sectiondatafromVerneretal.(1996). The general treatment in every such sub-bin follows M+06 3.3 Calculatingthephotonrates andisschematicallyillustratedinFig.3.Toavoidexpensiveinte- For a pure hydrogen medium Eq. (7) gives the photo-ionization grationsforeachphoto-ionizationcalculation,wetabulatethetotal rate,i.e.Γ=Γ andτ =τ .Inthecaseofamediumconsisting ionizationrateΓ asafunctionoftotalopticaldepthτ atthe HI HI tot tot ofhydrogenandhelium,thissimpletreatmentcanonlybeusedfor minimumfrequencyofthesub-bininquestion,assumingforallrel- photons with frequencies below the ionization threshold for neu- evantspeciesthesamefrequencydependenceforthecross-section, tralhelium,νHeI.Abovetheionizationthresholdfrequencyofion- followingtheconceptoutlinedinTenorio-Tagleetal.(1985).The th izedhelium,νHeIIallthreespecies,HI,HeIandHeIIIcontribute assumption of having the same frequency dependence allows the th totheopticaldepth.Inthefrequencybinbetweenthosethreshold summationoftheopticaldepthofeachspeciesandthetabulation frequencies, HI atoms and HeI atoms contribute to the total op- ofthephoto-ionizationrateasafunctionofthissingletotaloptical ticaldepth. Therefore,theminimumnumber offrequencybinsto depth.Sincethefrequencydependencesofthecross-sectionsofthe considerseparatelywhencalculatingthephotoionizationratesfor differentspeciesareinfactnotidentical(seeforexampleFig.C1), eachspeciesisthree,henceforthreferredtoas(frequency)bins.As thisisanapproximation. Themoresub-binsweusethemoreac- we explain below, we use a single frequency dependence for all curatelywefollowtheactualshapeofthe(frequencydependence) speciesineachbin.Sincethethedifferentfrequencydependences curves of the different cross-sections. The frequency dependence oftheionizationcrosssectionsσ ofthespeciesareverydifferent imposedonallspeciesisthatofσ forallsub-binsofbin2and HeI 6 Friedrichet al. eachsub-bin.Doingsooverestimatesthecontributionofhydrogen totheopticaldepthinfrequencybin2sincetheHIionizationcross section drops faster withfrequency than theionization cross sec- tionofHeI.ThisleadsthereforetoanoverestimateofΓ /Γ . HI HeI Infrequencybin3,theslopesofthecross-sectioncurvesareingen- eralmoresimilar(seeAppendixC)butstillsteeperforHIwhich resultsinanoverestimateofΓ /Γ .Themoresub-binsused, HI HeII thesmaller theoverestimateswill be. Wetest theconvergence of thisinAppendixE.Wefindthatincreasingthenumberofsub-bins inbin2improvestheslopeoftheionizationfrontwhileincreasing thenumber ofsub-binsinbin3improvestheionizationfractions insideandoutsidethefront. 3.4 Heatingandsecondaryionizations OneofthemotivationsforextendingC2-RAY,istoinvestigatethe effectsthatheliumandhardspectrahaveonthetemperatureevo- lutionoftheIGMduringtheEoR.Inthissectionwefirstdescribe theimplementationofthetemperaturecalculation.Thisisfollowed byremarksontheinclusionofsecondaryionizationsandconsider- ationsabouttheadditionaltimesteprestrictionscausedbythetem- peraturecalculation(inadditiontotheionizationcalculation). WeusetheidealgaslawPV = NkBT,wherekB isBoltz- mann constant, to calculate the pressure P from the temperature T and the gasparticlenumber density N/V = nH +nHe +ne in each cell of volume V. From the pressure we then calculate the internal energy per volume V (in each cell), using the di- mensionless heat capacity (per particle) at constant volume, cV: uint = Uint/V = cVP Foramonatomicgas,cV = 3/2.Sothe equationtoconverttemperatureintointernalenergydensityreads: 3 uint = 2kBT(nH+nHe+ne). (15) Thisenergydensityisaffectedbytheheating(H)andcooling(C) Figure 2.Extended iteration scheme for C2-RAY including helium with ofthegasper(cell-)volumeandperunittime.Inordertofollow couplingofthespecies. thetemperatureevolutionofthegas,wethereforehavetosolvethe equation ∂uint thatofσHeIIforallsub-binsofbin3.Thetableentriesarethein- ∂t =H−C. (16) tegralsofionizationrateoverthefrequencyrangeofthesub-binin As the only contribution to the heating rate H, we consider question,thusthetotalionizationrateduetophotonswithfrequen- photo-ionizationheating.Ineveryphoto-ionizationofspeciesi,the ciesinthatfrequencysub-bin. excessenergy,(ν−ν (i))istransferredtotheelectronreleased.If Inordertousethispre-calculatedtableanddeterminethetotal th oneassumesthatthecellsareopticallythick,(almost)allphotons photo-ionizationrateinagivencell,weneedtodeterminetheion- areabsorbedandtheaverageenergyperphotoionizationissimply izationcross-sectionsforeveryspeciesattheminimumfrequency theaverageexcessenergyover thewholespectrum. Thiswasthe of every sub-bin. We use the fitting formula from Verneretal. approachusedbytheoriginalC2-RAYinthetestswithtemperature (1996) for the cross-sections and compute the total optical depth evolutioninI+06.1 τtotasthesumoftheopticaldepthsofallspeciesattheminimum To properly take into account the effects of different optical frequencyofeachsub-bin.Thisopticaldepthisusedtodetermine depthontheheatingwecancalculateHinanalogytotheioniza- thetotalphotoionizationrateinthetable. tionrateΓ(Eq.7).Itshouldberememberedthat,althoughwewrite Next, this total photo-ionization rate has to be split up into thisasone equation, for thephoto-ionization rate itisin fact the ionizationratesforHI(Γ ),HeI(Γ )andincaseofbin3,HeII HI HeI differencebetweentheingoingphoto-ionizationrate(calculatedon (Γ ).Thisisdoneaccordingtotherelativefractionsofoptical HeII thebasisofhτi)andtheoutgoingphoto-ionizationrate(calculated depthineachfrequencyinterval. onthebasisofthehτi+h∆τi).Fortheheating,thesequantities Γi =Γtotττtoit, τtot = τi, i=HI,HeIand HeII (14) athreefroatrhmeroafbpshtroatcotnssinwcheentheeyntesyrimngbothliezecethlleaenxdctehsesaemneorugnytsotifllexin- i X cessenergystillintheformofphotonsleavingthecell.Neverthe- InAppendix D weevaluate other approaches for thisdistribution thathavebeenproposedintheliterature. Wecalculatethesefractionsofopticaldepthattheminimum 1 The H-only C2RAY used in Mellemaetal. (2006a) and Arthuretal. frequencyofeachsub-bin.Thisimpliesthatwealsohereassumea (2011)actuallyusedaonespecies,onefrequencybinversionoftheheating singlepowerlawfitforthefrequencydependenceofallspeciesin methoddescribedbelow. Radiativetransferofenergeticphotons:X-raysandheliumionizationin C2-RAY 7 Figure3.Sketchshowinghowtocalculatetheionizationratecontributionforbin3.Instep1,thetotalopticaldepthineachsub-biniscalculated.Instep 2,thepre-calculatedphotoionizationratetableisusedtocalculatethetotalphotoionizationrateineachsub-bin.Onthebasisoftheopticaldepthofeach speciesattheminimumfrequencyνminineachsub-bin,thefractionsfigoingintohydrogen-,helium-andionizedhelium-ionizationarecalculatedinstep 3.Finally,theproductsofthefractionsandtotalphotoionizationratesaresummedoverallsub-binstoresultintheionizationrateforeachspecies,ΓH,ΓHeI andΓHeII. less, what is tabulated as function of optical depth is this excess ferred energy is greater than the binding energy of this electron, energyasafunctionofopticaldepth.Thecorrespondingequation this electron is released. This process is called secondary ioniza- toEq.(7)isthen tion. The probability of such an event depends on the energy of the primary electron and on the ionization state of the gas. The H(i)= h(ν−νth(i))Lνe−hτνi1−e−h∆τνidν fullprocessrequirescarefulmodelling,butShull&vanSteenberg Zν(sub−bin) hν Vshell (1985),Ricottietal.(2002)andValde´s&Ferrara(2008),toname (17) a few, published separable functional relations dependent on pri- maryelectronenergyandhydrogenionizationfractiontodescribe forspeciesi.Sincetheexcessenergyisdifferentforeachspecies, what fraction of the energy deposited in primary electrons goes instead of having one table of heating rates in each sub-bin, we into secondary ionizations of either HI or HeI, fHI/HeI, and what ion need three tables, where the excess energy is with respect to the intoheat,f .Intheserelationsitisassumedthatx = x heat HII HeII thresholdfrequencyofνth(HI),νth(HeI)andνth(HeII).Theto- andsecondaryionizationsofHeIIareneglected.Weimplementthe talheatingineachfrequencysub-binisnaturallyanalogoustothe separable functional relationship as in Ricottietal. (2002) which photo-ionizationgivenby converges for high electron energies to the functional form of H= τi H(i). (18) Shull&vanSteenberg(1985). τ i tot Furlanetto&Stoever(2010)pointedout thatthereisingen- X eralnosimpleseparablefunctionalformfordecidingthefractions Someoftheelectronsthatreceivedexcessenergy(ν−ν (i)) thatgointosecondary ionizationsandheating. Ratherthedepen- th after being released from the bound state of atom i, will collide dence of the relativefractionsvaries withionization fractionina withboundelectronsofotheratoms(orions)ratherthanotherfree complex way, see their figure 5. However, the complex relation electrons.Thistransfersenergytotheboundelectron.Ifthetrans- givenbyFurlanetto&Stoever(2010)stillassumesx = x HII HeII 8 Friedrichet al. andno secondary ionizations of HeII.Giventheselimitationswe Table3.Parametersfortestingthecode. decidedthatatthispointthereisnoclearbenefitinimplementing thesemorecomputationallyexpensiverelations. nH/cm−3 constanthydrogennumberdensity Forthecoolingrate,weincludefree-freeandrecombination nHe/cm−3 constantheliumnumberdensity coolingforHII,HeIIandHeIIIandcollisionalexcitationcooling N˙γ/s rateofphotonsintheenergyinterval for H I, He I and He II. A full overview of the radiative cooling [13.6,5441.6]eV ratesusedisgiveninAppendixG.Incaseofcosmologicalsimu- Teff/K effectivetemperatureoftheblackbodysourcein lations,cosmologicalcoolingduetotheexpansionoftheuniverse caseofBBsource and Compton cooling against the cosmic microwave background β powerlawindexincaseofPLsource photonsareincludedaswell. Tini,x initialtemperatureandinitialionizationstates,HII, HeIIandHeIII In order to numerically solve Eq. (16) we use forward Eu- ∆r thecellsize ler integration. However, the cooling rate depends sensitively on ∆t thetimestep the gas temperature which is changing because f the heating. To (n2/n3) numberofsub-binsinfrequencybins2and3 accuratelyfollowthisbehaviourweusetheforwardEulermethod OTS/U-OTS/αA assumingthe(U-)OTS-approximation,orusing withsub-timestepsdeterminedbyalimitonthetemperaturechange αArecombinationrates (typically10%).WekeeptheheatingtermH(whichrepresentsthe averageheatingrateoverthetimestep)constantbutusethetemper- aturefromtheprevioussub-timesteptocalculatethecoolingrate C. time,anaccuratevalueforthetemperaturerequirestimestepsofthe orderoftheionizationtime.Thisconstraintbecomesmorestrictif Above,wedescribedhowtocalculatetheheatingrateinanal- thecellsareveryopticallythick. ogy to the photo-ionization rate. Specifically, we use the time- averagedopticaldepth(tothecellandinthecell)bothinEq.(7) Given these conclusions it is difficult to provide a simple andEq.(17).Thereishoweveranimportantdifferencebetweenthe recipeforchoosingthetimestep.Wethereforerecommendtesting dependence ofthephoto-ionization ratesandtheheatingrateson fornumericalconvergenceifaccuratetemperaturesarerequired. opticaldepth.AlthoughEq.(7)givesthecorrectrateofabsorbed photons,fortheheatingthefrequencyofthephotonsmatters.Let usfirstconsiderthesourcecell,i.e.,theopticaldepthtothecell,τ, 4 TESTINGTHECODE isconstantwithtimeandequal0.Whileforthephotoionization, Totestthevalidityoftheapproximationswemadeforthesakeof ∆tΓ(τ(t′))dt′ Γ(hτi )= 0 , (19) codeefficiency,weperformedaseriesoftests.Firstofallwevali- ∆t R ∆t datedthenewversionagainsttheoldversionofC2-RAYbysetting fortheheatingingeneral the helium abundance to a very low value. This test showed that thenewversionhasthesamephoton-conserving propertiesasthe ∆tH(τ(t′))dt′ H(hτi )≤ 0 , (20) methodpresentedinM+06.Fortestingthetimedependentsolution ∆t ∆t withhelium,includingtheOTSapproximation,secondaryioniza- R becauseofthe(ν−ν (i))termintheintegralforH.Qualitatively tion and temperature evolution, we would need a fully validated th onecouldsaythatinthebeginningofthetimesteptheheatingper timedependentphoto-ionizationcode,whichasfarasweknowis photo-ionizationistheopticallythickvalue(averagephotonenergy notpubliclyavailable.Wewerethereforeforced,justasforexam- overthewholespectrum),forlatertimesitshiftstolowervalues, pleBaeketal.(2010)andPawlik&Schaye(2011)tocompareour thelimitofwhichisgivenbytheopticallythinvalue(averagepho- results against the one-dimensional photo-ionization equilibrium ton energy over the spectrum weighted withthe photo-ionization code CLOUDY (version 08.00, last described in Ferlandetal. cross-section) 1998). Below we present the results of two sets of such one- Forthisfirstcellwheretheonlydependenceisontheoptical dimensional tests, one set in which the source has a black body depthinsidethecell,itwouldbepossibletotabulateacorrection (BB)spectrumwithaneffectivetemperatureofTeff,theotherfor factorbasedonthechangeoftheopticaldepthduringatimestep. a source with power-law (PL) like spectrum, where the energy- However, inthemoregeneral casewheretheopticaldepthtothe distributioncanbedescribedasL(ν) ∝ ν−β.FortheBBsource cell depends on the time varying optical depth of the other cells wetestvariousaspectsofthecalculationoftheionizationfractions along the ray, a local fix (such as using an energy-average based whilekeeping the temperature constant. For the PL case we also on thechange of optical depth inthecell or sub-timestepping on consider thetemperatureevolution. Wechosethisapproach since a cell-by-cell basis) to the heating rate is no longer possible. We thelattercasehavingawideenergyrangeofhardphotons,consti- exploreddifferentapproachesforcalculatingaccurateheatingrates tutes a more difficult test for the photo-heating. Table 3 gives an usinglargetimesteps,butwewereunabletofindageneralsolution. overviewoftheparametersusedforalltheone-dimensionaltests. InAppendixFweinvestigatehowtheheatingdependsonthe To demonstrate the multi-dimensional and multi-source capabili- choiceofthetimestep.Fromthesetestsweconcludethefollowing: tiesoftheextendedC2-RAYwealsopresenttheresultsoftwotest problemsusingthree-dimensionalcosmologicaldensityfields. • Ifweareinterestedintimescaleslargerthantherecombina- tiontime,wecanstilluselargetimesteps,astheinitialheatingis 4.1 Test-suite1:Blackbodysource nolongerdominating. • Ifthecoolingtimet <∆twecanalsouselargetimesteps, Thefirstsetoftestsconsiders theexpansion ofanionizedregion cool as the temperature is set by the equilibrium heating and cooling produced byasinglesourceintoaconstant densitymediumwith rates.Thisisthecasefortypicalinterstellarmediumconditions. constanttemperature.Undertheassumptionofsphericalsymmetry • If we are interested in time scales below the recombination thiscanbecalculatedinonedimensionfortheradialdistancerto Radiativetransferofenergeticphotons:X-raysandheliumionizationin C2-RAY 9 thesource.Theparametersthatarenotvariedinthissetoftestsare: also plot the ionization profiles after 1010 yr for the latter ap- nH/cm−3 =1.0×10−3,nHe/cm−3 =8.70×10−5,N˙γ/s=5.0× proach (thin black lines). It can be seen that the differences be- 1048,T /K=105.T =10000K,x=(10−40,10−40,10−40) tweentheOTSandU-OTSresultsarelargerthanbetweentheOTS eff ini (i.e.completelyneutral),∆r/pc=150,∆t/yr=107and(n ,n )= approximationandthefullradiativetransferresults,whencompar- 2 3 (10,14). ingthe(almost)equilibriumsolutions.Fromthisweconcludethat Wechoosetheseparameterssothattheheliumfractionf = althoughmorecomplicatedtoimplement,thecoupledOTSapprox- He n /(n +n ) = 0.08andtheremainingphysicalparameters imationshouldbethepreferredapproach. He H He arethesameasinTest 2ofI+06,exceptthatwedonotfollow SinceC2-RAYcannot(now)dealwiththediffusephotonsin thetemperatureevolutionhere.Thetimestep∆tis∼ 10%ofthe another way than using the OTS approximation, we are not able recombinationtimeinafullyionizedmedium,n αA,B,and∼104 to investigate the effects of using this approximation during the H timestheionizationtimeforthefirstcell. growth of the H II region. Cantalupo&Porciani (2011) compare resultsfromwhat wecallthe U-OTSapproximation withfullra- diativetransferofdiffusephotonsusingtheircodeRADAMESH. 4.1.1 Test1A Theyfindratherlargeeffectsattimeswhentheequilibriumsolu- tionhasnotyetbeenreached.Thiswarrantsfurtherinvestigation, Inthefirsttestofthissuiteweassumethatallionizingphotonsfrom probablybestpursuedinmannersimilartoI+06,involvingresults recombinationscanescape,whichmeansthatweareusingtheαA frommultiplecodes. rates.Thistestismeanttoshowhowwellourbasicapproachcan matchtheequilibriumsolutionfromCLOUDY.InFig.4weshow the resultstogether withthat equilibrium solution. It can be seen 4.2 Test-suite2:Apowerlawsource thataftert = 1010 yrthegeneralagreement withthe CLOUDY Inasecondset oftests,weuseUV-photon emittingsourceswith equilibriumsolutionisexcellentfordistancesbelow12kpc,butnot apower-lawspectrum,L(ν) ∝ ν−β,implyingthatthenumberof forlargerdistances.However,sincetherecombinationtimescaleis given by (n α)−1, it follows that for ionization fractions below photonsgoesasf(ν)∝ν−(β+1).Fortheseweconsideredthetwo e ∼ 1%, the recombination timescale is larger than ∼ 1010 yr. In casesβ =1andβ=2.Wepresenttheresultsforsimulationswith temperature evolution using T = 102 K. We only use the full otherwords,theseouterregionshavenotyetreachedequilibrium. ini OTS approximation, apply the secondary ionizations and choose AnothernoticeablefeatureinFig.4isthattheHeIIIfractionshows thenumberofsub-binsinbin2and3tobe(n ,n ) = (26,20). atransientbumparoundr = 4kpc.Thisbumpdisappearsasthe 2 3 TheremainingparametersareasinTEST 1. equilibrium solution is approached, but both, our results and the Wefoundatimestepof∆t=105yrtobesufficienttoobtain equilibriumsolutionstillshow aslopechangeintheHeIIIcurve convergenceinthetemperatureevolution.Choosingatimestepof around thatsamedistance. Fromallthisweconclude that there- the order of the ionization time of the first cell (∆t = 103 yr) sultsofthistestshowthatourbasicmulti-frequencyradiativetrans- gavetemperature results whichfor cellsclose tothe source were ferisconsistentwithCLOUDY’s. only different by at most 7%. For theparameters of this test, the opticaldepthattheionizationthresholdforhydrogenforonecell isτ ∼ 3.2,inbetweenthetwocasestestedinAppendixFandso 4.1.2 Test1B thisresultforthetimestepisconsistentwiththose. Inthesecond test,TEST1 B,weusetheonthespotapproxima- WeshowtheprofilesforHI/HIIandHeI/HeII/HeIIIandT tionasdescribed inSection3.1.2and otherwiseset the samepa- forthreetimesandthetwodifferent power lawspectrainFig.6. rameters as in the previous test. We show the ionization profiles WealsoplottheequilibriumsolutionsfromCLOUDYwithOTS inFig.5togetherwithtworesultsfromCLOUDY,oneusingthe forthosetwospectra.Wefindthatatt = 109 yrthe C2-RAYre- OTS approximation (green dashed line) and the other using full sults are close to the CLOUDY ones for those distances where radiativetransferofrecombinationphotons(bluedashedline).As the equilibrium solution has been reached, although not as close isapparent inFig.5,theagreement betweenthe C2-RAYandthe asfor theresultsof TEST 1.Thesesomewhat larger differences OTSsolutionof CLOUDYisgoodforthosedistanceswherethe withtheCLOUDYresultscannotbeexplainedbyeffectsfromthe equilibriumsolutionhasbeenreached(r<12kpcatt=1010 yr, temperaturecalculationsincewefoundcomparablediscrepancies justasinTEST1 A).Thegeneralpatternofevolutionoftheion- whenimposingaconstanttemperature.However,overallthematch izationfractionsisquitesimilartotheoneinTEST1 A,although isstillreasonable. theOTSapproximationdoesofcoursechangethedetailedvalues The time evolution shows that the two values of β produce ofthefractions. similarresults,withinitiallyquitesteepprofilesfortheionization TheCLOUDYsolutionusingfullradiativetransferofrecom- fractionsforalltimesandarelativelysteepfrontinthetemperature bination photons can be used to evaluate the validity of the OTS evolutionevenforthequitehardβ = 1spectrumattherelatively approximationforthisparticulartestproblem.InspectionofFig.5 latetimeoft = 107 yr.Comparingtheβ = 1andβ = 2results shows that the error introduced by using on the spot approxima- itcanbenoticedthattheformergivesahigherdegreeofhydrogen tion is barely noticeable, but larger for hydrogen and largest for andheliumionizationoutsidethefrontthanthelatter.Furtheritcan theHIfractioninsidetheHIIregion.Closerinspectionshowsthat beseenthattheresidualneutralhydrogenfractioninsidethefront theneutralhydrogenprofileinsidetheHIIregion,uptoaneutral ishigherfortheresultswithβ = 1.JustasinTEST1,thebumps fractionofabout1%,followsthecurvefromtheαArecombination inthe HeIII fractionare transient phenomena, although even the coefficientfromTEST1 A.Thisisnotsurprisingsincemostofthe equilibrium solution shows a slope change around a distance of recombinationphotonsfromthishighlyionizedregionwillnotbe ∼6kpc. absorbedonthespot,butescapetolargerdistances. Tittley&Meiksin(2007)reportedtheHIIfronttobetrailing TogaugetheimportanceofapplyingthecoupledOTSrather behindtheHeIIIfrontforspectrawithβ <1.8.However,evenin thanthepopular U-OTS(i.e.usingαB rates,seeSect.3.1.1),we ourβ =1resultswefindtheHIIfractiontobealwayslargerthan 10 Friedrichet al. 1 t/yr=107 H I H II t/yr=108 t/yr=109 t/yr=1010 0.1 CLOUDY ni o cti a fr 0.01 1 He I He II He III 0.1 ni o cti a fr 0.01 1 5 10 15 1 5 10 15 1 5 10 15 Distance /kpc Distance /kpc Distance /kpc Figure4.ResultsfromTEST1 A.FractionsofHI(upperleftpanel),HII(uppermiddlepanel),HeI(lowerleftpanel),HeII(lowermiddlepanel)andHeIII (lowerrightpanel)attimest/yr=[1×107,1×108,1×109,1×1010]asindicatedbylinethicknessinthelegend.Wealsoshowtheequilibriumsolutionof CLOUDY(greendashed).Notethatduetothelowelectrondensityintheouterregions,theequilibriumsolutionhasnotyetbeenreachedbeyondadistance ofabout12kpcatt=1010yr. theHeIIIfraction.Wewereabletoobtain(double)frontcrossings 4.3.1 TEST3:Effectofheliumonthemorphologyofthe when using the αA recombination rates (no OTS) and disabling hydrogenionizationfractionfieldduringEoRwithout thesecondaryionizations.However,thelocationswheretheHeIII temperatureevolutionforstellartypesources fractionwasfoundtobelargerthantheHIIfractionwereinalim- ited area in the ionization fraction-time space, roughly in the in- Many cosmological reionization simulations only include hydro- terval [0.1, 108 yr] to [0.01, 109 yr]. We therefore conclude that gen and implicitly assume that helium is singly ionized every- trailingHIIfrontsareatbestamarginaleffect. wherewherehydrogenisionized.Theusedhydrogennumberden- In order to evaluate the effects of secondary ionizations we sity is therefore equal to the total number density. In this sec- alsoranTEST 2withoutthem.Wefoundthatthisledtochanges tion, we test if the morphology of H II regions in a reioniza- larger than the differences we found between our results and the tionsimulationwithstellarsourceschangesifheliumisincluded. equilibriumsolutionofCLOUDY.Wethereforeconcludethatitis For this comparison, we use simulation 53Mpc g8.7 130S from worthtoincludesecondaryionizationsinthecalculation. Friedrichetal. (2011) and Ilievetal. (2011). The electron scat- teringopticaldepth produced bythissimulation, τ = 0.083, is es consistent withthe 1–σ range allowed by theseven year WMAP results,τ = 0.088±0.015(Komatsuetal.2011).Inthissimu- es lation,thenumberofionizingphotonsproducedbyadarkmatter haloofmassM isdefinedthrough MΩ 4.3 Cosmologicaltests N˙γ =gγ10Ωmmb p , (21) Inthissection,wetestthe3D-versionoftheextendedC2-RAYon where N˙γ is the number of ionizing photons emitted per Myr, twocosmologicaldensityfields.Thefirsttestisalargerscalecos- Ωb = 0.044, Ωm = 0.27andmp istheprotonmass(Thisequa- mologicaltestwithouttemperatureevolutiontoevaluatetheeffect tionincludedincorrectlyaµinFriedrichetal.2011, equation1). ofheliumatafixedtemperatureinasimulationwithsourceswith Massivehalosareassignedanefficiencygγ =8.7whilelowmass asoftspectrum.Thesecondtestisarathersmallcosmologicalvol- sourceshaveanefficiencyofgγ = 130andaresuppressedinre- umebutincludestemperatureevolution:weredothecosmological gions where the ionization fraction (of hydrogen) is higher than test-problem with multiple sources from I+06 (their test4) to test 10%.Toevaluate theeffectofheliumonthemorphology of HII theeffectofheliumontheheatingandonthehydrogenionization regionsweusethedimensionlesspowerspectrumoftheHIIfrac- fractionfield. tion∆2 InFig.7weshowthepowerspectraforthissimulation xx