Metal Cooling in Simulations of Cosmic Structure Formation
aa r X i v : . [ a s t r o - ph ] J a n Mon. Not. R. Astron. Soc. , 1–19 (2007) Printed 2 November 2018 (MN L A TEX style file v2.2)
Metal Cooling in Simulations of Cosmic StructureFormation
Britton Smith ⋆ , Steinn Sigurdsson ⋆ , and Tom Abel ⋆ Astronomy & Astrophysics, The Pennsylvania State University, 525 Davey Laboratory, University Park, PA 16802, U. S. A. Kavli Institute for Particle Astrophysics and Cosmology, Stanford Linear Accelerator Center, 2575 Sand Hill Road, MS 29,Menlo Park, CA 94025, U. S. A.
In original form May 23, 2007
ABSTRACT
The addition of metals to any gas can significantly alter its evolution by increasingthe rate of radiative cooling. In star-forming environments, enhanced cooling can po-tentially lead to fragmentation and the formation of low-mass stars, where metal-freegas-clouds have been shown not to fragment. Adding metal cooling to numerical sim-ulations has traditionally required a choice between speed and accuracy. We introducea method that uses the sophisticated chemical network of the photoionization soft-ware, Cloudy, to include radiative cooling from a complete set of metals up to atomicnumber 30 (Zn) that can be used with large-scale three-dimensional hydrodynamicsimulations. Our method is valid over an extremely large temperature range (10 K T K), up to hydrogen number densities of 10 cm − . At this density, a sphere of1 M ⊙ has a radius of roughly 40 AU. We implement our method in the adaptive meshrefinement (AMR) hydrodynamic/N-body code, Enzo. Using cooling rates generatedwith this method, we study the physical conditions that led to the transition fromPopulation III to Population II star formation. While C, O, Fe, and Si have been pre-viously shown to make the strongest contribution to the cooling in low-metallicity gas,we find that up to 40% of the metal cooling comes from fine-structure emission by S,when solar abundance patterns are present. At metallicities, Z > − Z ⊙ , regions ofdensity and temperature exist where gas is both thermally unstable and has a coolingtime less than its dynamical time. We identify these doubly unstable regions as themost inducive to fragmentation. At high redshifts, the cosmic microwave backgroundinhibits efficient cooling at low temperatures and, thus, reduces the size of the doublyunstable regions, making fragmentation more difficult. Key words: stars: formation methods: numerical
The first luminous objects in the universe formed from pri-mordial gas, comprised solely of H and He, with only traceamounts of D an Li. The relatively simple chemistry ofmetal-free gas, combined with tighly constrained cosmolog-ical parameters (Spergel et al. 2007), has allowed the for-mation of the first stars to be simulated with extremelyhigh precision, from the hierarchical growth of their hostdark matter halos through to the point where the denseproto-stellar cores becomes optically thick (Abel et al. 2002;Bromm et al. 2002; Bromm & Loeb 2004; Yoshida et al.2006; O’Shea & Norman 2007; Gao et al. 2007). With thedeaths of these stars came the creation of the first heavy ⋆ E-mail:[email protected] (BDS); [email protected](SS); [email protected] (TA) elements. Core-collapse and pair-instability supernovae cre-ated metals in copious amounts (Heger & Woosley 2002)and ejected them into the IGM (Madau et al. 2001).The presence of metals alters the dynamics of collaps-ing gas-clouds by increasing the number of available atomicand molecular transitions, allowing the gas to lose its inter-nal energy more quickly than in case of no metals (Omukai2000; Bromm et al. 2001; Bromm & Loeb 2003b). The in-troduction of metals adds a new level of complexity to theproblem of simulating the formation and evolution of cos-mic structure. Abel et al. (1997) identified a minimal setof 21 chemical reactions necessary for accurately followingthe non-equilibrium evolution of a gas consisting solely ofspecies of H and He, including H . Galli & Palla (1998)showed that 33 total reactions were required when includingD and Li species to the gas. Omukai (2000) performed oneof the first numerical studies of collapsing gas-clouds to con- c (cid:13) Britton Smith, Steinn Sigurdsson, and Tom Abel sider the contribution of metals. Their chemical network ofH, He, C, and O included 50 atomic and molecular speciesand 478 reactions. While theirs was not a minimal model,the above examples illustrate the great expense associatedwith the expansion of chemical networks to include addi-tional elements. Other works have studied the effect of met-als on star-forming gas using similar methodologies to thatof Omukai (2000), e.g., Schneider et al. (2002, 2003, 2006);Omukai et al. (2005). The complexity of the chemical net-works used in these studies limited their treatment of gasevolution to one-zone, semi-analytical models. In the earli-est work to incorporate metal cooling into three-dimensionalhydrodynamic simulations to study metal-enriched star for-mation, Bromm et al. (2001) used a small set of the mostdominant atomic transitions of C, N, O, Fe, Si, and S,as decribed by Ricotti et al. (1997). Their method also ig-nored the cooling from H , which was justified within theirstudy by the assumption of a very large photo-dissociatingUV background, but is otherwise an extremely importantcoolant in low-metallicity environments. For high temper-ature gases, Sutherland & Dopita (1993) computed metalcooling functions that included 14 heavy elements over arange of metallicities, with solar abundance patterns. Thesecooling functions are useful for simulating the IGM andother hot, ionized environments, but a minimum tempera-ture of 10 K makes them inapplicable to studies of the cold,neutral gas associated with star-formation. These coolingfunctions assume collisional equilibrium of the species andas such cannot capture the important role of UV and X-rayradiation.We introduce a new method for including the cool-ing from heavy elements in large-scale hydrodynamic sim-ulations that is valid over a wide range of physical con-ditions, covers a great number of elemental species, andis fast enough to be used in large-scale numerical sim-ulations. We have utilized the established photoioniza-tion software, Cloudy (Ferland et al. 1998) to constructlarge grids of metal cooling data. We have developed amethod to include both the cooling from heavy elementsand the non-equilibrium cooling from H in hydrodynamicsimulations. This method has been used successfully inthe numerical simulations of star formation performed bySmith & Sigurdsson (2007). In §
2, we describe our methodfor creating the metal cooling data, including a new code toexpedite the process. We, then, present two implementationsof the cooling method in the AMR, hydrodynamic/N-bodycode, Enzo (Bryan & Norman 1997; O’Shea et al. 2004). In §
3, we focus on the application of metals to low-temperatureenvironments, identifying the dominant cooling mechanisms,and studying the possibility of fragmentation and thermalinstability in metal-enriched gas. Finally, we end with a dis-cussion in § At the current time, it is still too computationally expen-sive and memory intensive to follow the non-equilibriumchemistry for a large set of heavy elements in a three-dimensional hydrodynamic simulation. The exact mass of the first massive stars is not known (Abel et al. 2002;Tan & McKee 2004; Yoshida et al. 2006). Also unknown arethe exact yields of early supernovae (Heger & Woosley 2002;Maeder et al. 2005; Nomoto et al. 2006; Rockefeller et al.2006). Similarly, in many astrophysical systems one mightwant to model computationally the exact metal distribu-tions. Consequently, it is not clear a priori what level ofsophistication of cooling model is needed to adequately cap-ture the hydro and thermodynamic evolution of the gas un-der consideration. Note that uncertain grain physics alsoincreases the potentially important parameter space. In ourapproach, we assume ionization equilibrium, which allows usto calculate, in advance, the cooling rate for a parcel of gaswith a given density and temperature, with incident radia-tion of known spectral shape and intensity. For this problem,we find the photoionization code, Cloudy (Ferland et al.1998), especially apt. Cloudy is conventionally used to modelthe transmitted spectrum from a cloud of gas with a givenchemical composition, being irradiated by a specified source.The code must calculate an equilibrium solution by balanc-ing the incident heating with the radiative cooling from afull complement of atomic and molecular transitions, as wellas continuum emission from dust. The chemical network ofCloudy covers all atomic species from H to Zn, as well as amultitude of molecular species. Each elemental abundancecan be specified individually, giving us the ability to modelthe cooling from a gas with any composition. Since Cloudypermits the use of virtually any input spectrum, we are ableto create cooling data that is suitable for any radiation en-vironment. Instead of allowing the code to cycle throughtemperatures until converging on a thermodynamic equi-librium solution, we use the constant temperature com-mand to fix the temperature externally, allowing us to utilizeCloudy’s sophisticated machinery to calculate cooling ratesout of thermal equilibrium. In this manner, we create a gridof heating and cooling values as a function of temperature,gas density, chemical composition, and incident spectrum.The cooling rates presented in this work were created usingversion 07.02.01 of the Cloudy software.To automate the process of data production and orga-nization, we have created a code, called ROCO (RecursivelyOrganized Cloudy Output.) ROCO uses a recursive algo-rithm to process user-specified loop parameters, making itpossible to create data-grids of any dimension. Commandsthat are to be issued to Cloudy are given to the ROCO codein either one of two formats - loop commands with a setof parameters through which the code will iterate, and con-stant commands that are to be issued with the same valueduring each iteration over the loop commands. Since mostuses of Cloudy involve the creation of large grids of modelsconstructed by looping over a set of input parameters, thecapabilities of ROCO give it the potential to be useful toa broader community of Cloudy users than just those whowould use it to create the cooling tables discussed here. Tothis end, ROCO is structured in such a way that the post-Cloudy data analysis routines can be easily interchanged tosuit the needs of different users. The code features an extrarunning mode that simply runs Cloudy over the specifiedparameter-space with no further processing of the data, aswell as a template designed to help users create new run-ning modes suited to their specific needs. ROCO also hasthe ability to run multiple instances of Cloudy simultane- c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation ously, greatly reducing runtime. The parallel feature workswell on individual machines with multiple processors, as wellas Beowulf clusters using the MPI (Message Passing Inter-face) framework. A copy of the ROCO code will be madeavailable upon request to the authors.In Figure 1, we display the resulting cooling function forgas with n H = 1 cm − at metallicities, from Z = 0 (metal-free) to 10 Z ⊙ , over the temperature range, 50 T K. For these cooling rates, we use the coronal equilibrium command in Cloudy to simulate an environment free of ra-diation, where all ionization is collisional. We also neglectthe cooling from H , so as to better illustrate the coolingcontribution from metals at temperatures less than 10 K.We accomplish this by issuing the Cloudy command, no H2molecule . We implement our metal cooling method in the Eule-rian adaptive mesh refinement hydrodynamic/N-body code,Enzo(Bryan & Norman 1997; O’Shea et al. 2004). When asimulation is initialized, Enzo reads in the Cloudy/ROCOdata-grid, storing the heating and cooling values as functionsof temperature, H number density, and any other parame-ters, such as spectral intensity, depending on the nature ofthe simulation. The heating, Γ, and cooling, Λ, are storedwith code units corresponding to [ergs s − cm ]. During thesimulation, Enzo stores the mass density and internal energyfor each grid cell in the box. At each hydrodynamic time-step, the radiative cooling solver cools the gas by loweringthe internal energy via a simple Euler update, u n +1 i,j,k = u ni,j,k + ˙ u ni,j,k × δt, (1)where u n +1 i,j,k denotes the internal energy of the grid cell with(x,y,z) coordinates, (i,j,k), at the (n+1)’th time-step, ˙ u is thecooling rate in code units corresponding to [ergs s − ], and δt is the adopted time-step. For every hydrodynamic time-step, the code subcycles through Equation 1, selecting fromthree possible time-steps, until one full hydrodynamic time-step has been completed. The time-step, δt , in Equation 1,adopts the minimum of the following three values: (1) halfof the hydrodynamic time-step, (2) 10% of the cooling time,(u/ ˙ u ), or (3) the time remaining to have integrated over onefull hydrodynamic time-step. The Euler update is the stan-dard method used for updating the internal energy in all ofthe established cooling routines in the Enzo code. Since thecooling rate is such a nonmonotonic function of temperature,this approach of using substeps with an explicit integrationmethod has been found to yield the best combination ofaccuracy and speed (Anninos & Norman 1994). The inter-nal energy and mass density for each grid cell are convertedto temperature and number density and the heating andcooling values are calculated by linearly interpolating overthe Cloudy/ROCO data-grid. The change in internal energyfrom the Cloudy/ROCO cooling rates is expressed as˙ u C/R = (Γ − Λ) n H , (2)where n H is the H number density.We implement two distinct versions of the method de-scribed above. In the first and simplest version, the cooling iscalculated solely from the Cloudy/ROCO data, as in Equa-tion 2. The total change in internal energy is ˙ u tot = ˙ u C/R . (3)When converting the internal energy to temperature, it isnecessary to know the value of the mean molecular weight, µ . In this implementation, we assume µ to be a constantwith the value 1.22. For high temperatures, T ∼ > K, thismethod is sufficient for providing accurate gas cooling withinhydrodynamic simulations. This implementation is not suit-able, however, when T ∼ < K and the formation of H be-comes important. Disregarding formation on grain surfacesand three-body formation, H primarily forms through thefollowing channels:the H − channel,H + e − → H − + γ ,H − + H → H + e − , (4)and the H +2 channel,H + H + → H +2 + γ ,H +2 + H → H + H + . (5)When a significant electron fraction exists, these reac-tions proceed to form H very quickly, with the H − chan-nel typically dominating, except in the very high redshiftuniverse (z > − is readily destroyed bythe CMB (Abel et al. 1997; Bromm et al. 2002). Recently,Hirata & Padmanabhan (2006) have suggested that forma-tion of H +2 viaHe + H + → HeH + + γ ,HeH + + H → H +2 + He, (6)is responsible for more H than the H +2 channel in Equa-tion 5. If, however, the ionization fraction is low, H formsvery slowly, with equilibrium timescales that can exceed thecurrent age of the universe. The consequence is that H for-mation is so sensitive to the thermal history of the gas thatH fractions cannot be known without explicitly followingthe non-equilibrium chemistry during the simulation. Wefind this to be the case when using Cloudy to compute thecooling rate from H . In searching for ionization equilibrium,Cloudy integrates over timescales that are unphysically long,leading to an overcalculation of the H fraction, producingcooling rates that are too high.We solve this problem in our second implementation byfirst removing the H molecule from Cloudy’s chemical net-work with the no H2 molecule command. We, then, run twoCloudy/ROCO data-grids: one with the full set of elementsand the other with H and He only. We obtain a metals-onlydata-grid by subtracting the metal-free data-grid from thecomplete data-grid. Using the established H/He network inEnzo (Anninos et al. 1997; Abel et al. 1997), we follow thenon-quilibrium fractions of H, H + , H − , H , H +2 , He, He + ,He ++ , and e − , and directly calculate the associated atomic(Black 1981; Cen 1992) and molecular (Galli & Palla 1998)cooling rates. We provide the option to use the H cool-ing rates of Lepp & Shull (1983), which are obsolete, butallow a means of comparison to older simulations. We alsoinclude the cooling, or heating, from electrons scattering offthe CMB asΛ Comp = 5 . × − (1 + z ) n e ( T − T CMB ) , (7)where T CMB = 2.7 (1 + z) K (Bromm et al. 2002). We pre-vent the metals from cooling the gas below the CMB temper- c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel ature by subtracting the metal cooling rate at T = T
CMB , asin Bromm et al. (2001). Including the CMB explicitly in cos-mological simulations would require adding an extra dimen-sion in the Cloudy/ROCO data-grid to account for the evo-lution of the CMB with redshift. While this is certainly pos-sible, interpolating over an extra dimension to calculate thecooling during the simulation would be significantly slowerthan the approximation described above. To test the va-lidity of this approximation of the CMB temperature-floor,we create a Cloudy/ROCO metals-only data-grid, explicitlyincluding the CMB. At low densities (n ∼ − ), the val-ues of (Λ - Γ) from the data with the CMB included differfrom the values of (Λ - Λ( T CMB )) from the data without theCMB by a factor of roughly 2 near T
CMB . At higher tem-peratures and densities, the two values are nearly identical.The total rate of energy loss applied to the simulation gasin the second implementation is˙ u tot = ˙ u H,He + ˙ u Comp + ˙ u C/R , (8)where ˙ u H,He is the total atomic and molecular cooling fromthe H/He network, and ˙ u C/R is the metals-only cooling ratetaken from the Cloudy/ROCO data in the manner describedabove. The value of µ is calculated directly from the H/Hespecies fractions, neglecting any addition from the metals.In low-metallicity gases, this approach is reasonable, as theincrease in µ from the metals only reaches ∼ − for Z =10 − Z ⊙ . Since the Cloudy/ROCO data-grids also store thevalues of µ for each point, this can be added to the valuecalculated without the heavy elements when the metallicityis very high.In Figure 2, we display low-temperature cooling func-tions for gases with varying density and metallicity, con-structed with the third implementation of the metal coolingmethod. To produce the data for Figure 2, we set up an un-physical, two-dimensional grid in Enzo that varies smoothlyover density and temperature. We iterate the reaction net-work for a time equivalent to that between z = 99 and20 ( ∼
160 Myr), with hydrodynamics disabled, then com-pute the cooling with the third implementation of our metalcooling method, using the H cooling rates of Galli & Palla(1998). Since the first stars are predicted to form at the cen-ters of ∼ M ⊙ dark matter halos at z ∼
20 (Tegmark et al.1997; Yoshida et al. 2003), integrating the rate equationsover this time interval places each of the species in the rel-ative abundances in which they would be found during theepoch of first-star formation. Hence, Figure 2 provides a di-rect comparison of the cooling rate of the gas that formedthe first and successive generations of stars. The H frac-tions used to create the cooling rates for Figure 2 resultfrom integrating the H/He chemical network for the periodof time between z = 99 and 20. As such, Figure 2 shouldnot be used as a general cooling function, except in the con-text mentioned above. However, Figures 3–8 do not includethe cooling from H and may, therefore, be used as generalpurpose cooling functions when the cooling from H is notneeded. Much attention has been given recently to the role of thefirst heavy elements in transitioning from the singular,high-mass mode of star formation of the first stars to themode producing stars with a Salpeter initial mass function(IMF). Analytical studies by Bromm & Loeb (2003b) andSantoro & Shull (2006) have focused on the contributions ofindividual elements toward triggering fragmentation in star-forming clouds. Bromm & Loeb (2003b) suggest C and O tobe the dominant coolants in low-metallicity gas, in the pres-ence of an H dissociating UV background created by thefirst stars (Bromm & Loeb 2003a). By calculating the cool-ing rate necessary to equate the cooling time to the free-falltime at n = 10 cm − and T = 200 K, the point where H cooling becomes inefficient (Abel et al. 2002; Bromm et al.2002), Bromm & Loeb (2003b) predict individual criticalabundances of C and O to be [C/H] crit ≃ -3.5 and [O/H] crit ≃ -3.05, where [A/H] = log (N A /N H ) - log (N A /N H ) ⊙ .Santoro & Shull (2006) consider the cooling from Fe and Si,in addition to C and O, and take into account the densitydependence of metal cooling. In doing so, they find thatthe critical abundance of each element varies with density,reaching a minimum at a critical density that is different ineach case. They also note that different elements dominatedifferent density and temperature regimes.In Figs. 3–6, we plot the individual cooling contribu-tions for number density, n H = 10 cm − and metallicities,Z = 10 − Z ⊙ and 1 Z ⊙ . In each case, we create a set ofcooling data with the full complement of elemental species,from H through Zn, neglecting H . We plot only the coolantswhose contributions reach, at least, 10 − of the total coolingwithin the temperature range, 10 K T § I , instead of C II , as in Bromm & Loeb(2003b) and Santoro & Shull (2006). For Z = 10 − Z ⊙ , themost important coolants are the fine structure transtions at369.7 µ m and 609.2 µ m from C I and at 63.2 µ m from O I . Athigher temperatures (T >
100 K), the Si II transition at 34.8 µ m becomes important as well. At higher metallicities, COreplaces atomic C (Fig. 5) and emission from Si I at 129.6 µ mbecomes dominant (Fig. 6). Fe cooling is relatively unimpor-tant up to n H = 10 cm − , but is completely dominant forT >
200 K by n H = 10 cm − , with O I strongest slightlybelow that temperature, and CO the most important be-low 50 K (Fig. 7). In addition to the elements covered bySantoro & Shull (2006), we note that cooling from neutralS reaches roughly the 10% level at n H = 10 cm − . The [S I ]25.2 µ m transition peaks at 40% of the total cooling at n H = 10 cm − and T ∼ − of the total cooling isthe 60.6 µ m transition of [P II ] at n H = 10 cm − . The num-ber of coolants that reach 10 − of the total grows quicklywith density. We observe 23 distinct coolants contributingat that level at n H = 10 cm − , and 32 by 10 cm − . If welower the threshhold to 10 − , there are a total of 84 coolantsat n H = 10 cm − , illustrating the strength of Cloudy andour cooling method. c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation Recently, studies have suggested that dust cooling at highdensities can trigger fragmentation for metallicities as lowas 10 − Z ⊙ (Omukai et al. 2005; Schneider et al. 2006;Tsuribe & Omukai 2006; Clark et al. 2007). Schneider et al.(2004) have claimed that between 15 and 30% of the massof the progenitor of a pair-instability supernova is con-verted in dust. However, observations of the Crab nebula(Green et al. 2004) and the Cassiopeia A supernova rem-nant (Krause et al. 2004) have returned little or no signs ofdust, suggesting that type II supernova may not actuallybe large dust producers. Given the controversy surroundingthe existence of dust grains in the formation environments ofsecond-generation stars, we do not include the cooling fromdust in the analysis, but leave it for a separate work. Toprovide an example of the effect of dust on the cooling rate,we run a simple model including dust in Cloudy. We use amodel designed to simulate the dust within the ISM, usingthe Cloudy command, grains ISM . The dust physics used inCloudy is described in detail by van Hoof et al. (2004). TheISM dust model in Cloudy consists of both graphite and sil-icates with sizes ranging from 5 × − µ m to 0.25 µ m anda power-law size-distribution with a power-law index of -3.5.For solar metallicity, the total grain abundances per H are10 − . for graphite and 10 − . for silicates. In Fig. 8, weplot the cooling rate from metals at n H = 10 cm − , withand without dust grains, for metallicities, Z = 10 − Z ⊙ , 10 − Z ⊙ , and 10 − Z ⊙ . Since it is inappropriate to think of thedust and gas-phase metal abundances as independent, wedirectly scale the dust abundances with the gas-phase metalabundancs. For number densities lower than 10 cm − , theadditional cooling from dust is almost negligible. At higherdensities, dust becomes more important. If dust is, in fact,produced in the supernovae of the first stars, it is likely tobe as important as previous studies have claimed. In order to study the ability of a collapsing gas-cloud tofragment, we first identify the regions of density and tem-perature where the classical fragmentation criterion, t cool < t dyn (Field 1965), is met, with the dynamical time expressedas t dyn = r π Gρ . (9)We limit this analysis to solar abundance patterns. This rep-resents the first step of an incremental approach to studyingthe general criteria that lead to fragmentation in collapsingclouds. The use of solar abundance patterns will allow usto begin to quantify the chemical abundance required forfragmentation. In a following paper, we will study nonsolarabundance patterns motivated by predicted yields of pri-mordial supernovae. In this future work, we will vary theabundances of individual elements, which will require theexploration of a much larger parameter-space. We createa Cloudy/ROCO data-grid with the following parameters:50 K T δ log(T) ≃ − n H cm − with δ log(n H ) = 0.1, and 10 − Z ⊙ Z − Z ⊙ with δ log(Z) = 1. We, then, followthe same procedure used to produce Fig. 2, as described in § (t dyn /t cool ) >
0. As expected, there is no density atwhich metal-free gas can fragment for T <
200 K. As themetallicity increases, the value of log (t dyn /t cool ) slowly in-creases, first in the low-temperature regime, where H cool-ing is inefficient, so even a small amount of metals has aneffect. For gas with n H = 10 cm − at T = 200 K, the frag-mentation criterion is nearly met by Z = 10 − Z ⊙ and wellsatisfied one order of magnitude higher in metallicity. Oncethe metallicity reaches 10 − Z ⊙ , the entire parameter-spaceis fragmentable. At high densities, however, fragmentationwill be curtailed as the cloud becomes optically thick toits own radiation (Low & Lynden-Bell 1976; Rees 1976). Aswas also reported by Santoro & Shull (2006), the efficiencyof the metal cooling peaks at n H ∼ cm − , significantlylowering the critical metallicity required for fragmentation.The addition of metals to a gas also has the potentialto trigger thermal instabilities during cloud-collapse. As inAbel et al. (2002) (Field 1965), we define a parcel of gaslosing energy at a rate, L, to be thermally unstable if ρ (cid:16) ∂L∂ρ (cid:17)(cid:12)(cid:12)(cid:12) T − T (cid:16) ∂L∂T (cid:17)(cid:12)(cid:12)(cid:12) ρ + L ( ρ, T ) > , (10)where L is expressed in terms of the cooling rate, Λ, as L ( ρ, T ) = ρ Λ( ρ, T ) . (11)The cooling rate, Λ, can be locally approximated by a power-law in both temperature and density asΛ( ρ, T ) ∝ (cid:16) TT (cid:17) α (cid:16) ρρ (cid:17) β . (12)The partial derivatives of Eqn. 10 become (cid:16) ∂L∂ρ (cid:17)(cid:12)(cid:12)(cid:12) T = ( β + 1)Λ( ρ, T ) (13)and (cid:16) ∂L∂T (cid:17)(cid:12)(cid:12)(cid:12) ρ = ρα Λ( ρ, T ) T . (14)The thermal instability criterion simplifies to α − β < . (15)In Fig. 10 we plot the value of the instability parameter,( α - β ), for the same cooling data used for Fig. 9. Formetal-free gas (Fig. 10, top-left), the instability parameter isgreater than 3 over nearly the entire parameter space, andalways greater than 4 at high densities. Abel et al. (2002)and Yoshida et al. (2006) both arrive at the same conclu-sion using this analysis for metal-free gas, with the addedassumption that the cooling function is independent of den-sity. As the metallicity increases, a region of thermal in-stability forms at low density and temperature. When themetallicity reaches 10 − Z ⊙ (Fig. 10, middle-right), a sec-ond unstable region exists for 10 cm − ∼ < n H ∼ < cm − ,at a temperature of a few hundred K. The second unstableregion coincides with the increase in cooling efficiency to itsmaximum value, illustrated in Fig. 9.Fragmentation is more likely to occur when the gas isboth thermally unstable, and can cool faster than the dy- c (cid:13)000
200 K. As themetallicity increases, the value of log (t dyn /t cool ) slowly in-creases, first in the low-temperature regime, where H cool-ing is inefficient, so even a small amount of metals has aneffect. For gas with n H = 10 cm − at T = 200 K, the frag-mentation criterion is nearly met by Z = 10 − Z ⊙ and wellsatisfied one order of magnitude higher in metallicity. Oncethe metallicity reaches 10 − Z ⊙ , the entire parameter-spaceis fragmentable. At high densities, however, fragmentationwill be curtailed as the cloud becomes optically thick toits own radiation (Low & Lynden-Bell 1976; Rees 1976). Aswas also reported by Santoro & Shull (2006), the efficiencyof the metal cooling peaks at n H ∼ cm − , significantlylowering the critical metallicity required for fragmentation.The addition of metals to a gas also has the potentialto trigger thermal instabilities during cloud-collapse. As inAbel et al. (2002) (Field 1965), we define a parcel of gaslosing energy at a rate, L, to be thermally unstable if ρ (cid:16) ∂L∂ρ (cid:17)(cid:12)(cid:12)(cid:12) T − T (cid:16) ∂L∂T (cid:17)(cid:12)(cid:12)(cid:12) ρ + L ( ρ, T ) > , (10)where L is expressed in terms of the cooling rate, Λ, as L ( ρ, T ) = ρ Λ( ρ, T ) . (11)The cooling rate, Λ, can be locally approximated by a power-law in both temperature and density asΛ( ρ, T ) ∝ (cid:16) TT (cid:17) α (cid:16) ρρ (cid:17) β . (12)The partial derivatives of Eqn. 10 become (cid:16) ∂L∂ρ (cid:17)(cid:12)(cid:12)(cid:12) T = ( β + 1)Λ( ρ, T ) (13)and (cid:16) ∂L∂T (cid:17)(cid:12)(cid:12)(cid:12) ρ = ρα Λ( ρ, T ) T . (14)The thermal instability criterion simplifies to α − β < . (15)In Fig. 10 we plot the value of the instability parameter,( α - β ), for the same cooling data used for Fig. 9. Formetal-free gas (Fig. 10, top-left), the instability parameter isgreater than 3 over nearly the entire parameter space, andalways greater than 4 at high densities. Abel et al. (2002)and Yoshida et al. (2006) both arrive at the same conclu-sion using this analysis for metal-free gas, with the addedassumption that the cooling function is independent of den-sity. As the metallicity increases, a region of thermal in-stability forms at low density and temperature. When themetallicity reaches 10 − Z ⊙ (Fig. 10, middle-right), a sec-ond unstable region exists for 10 cm − ∼ < n H ∼ < cm − ,at a temperature of a few hundred K. The second unstableregion coincides with the increase in cooling efficiency to itsmaximum value, illustrated in Fig. 9.Fragmentation is more likely to occur when the gas isboth thermally unstable, and can cool faster than the dy- c (cid:13)000 , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel namical time. We indicate the regions where both the frag-mentation and thermal instability criteria are met in whitein Fig. 11. No doubly unstable realm exists for metallicities,Z − Z ⊙ . The CMB creates a temperature floor, below which gas can-not cool radiatively. We study the influence of the CMB onthe evolution of star-forming gas by applying a CMB floor atz = 20 in the manner described in § § α , from Equation 15, at low temperatures, making thegas thermally stable. In Fig. 12, we illustrate the influenceof the CMB on the doubly unstable regions, shown previ-ously. The unstable region that existed in the low-density,low-temperature regime is completely eliminated. There re-mains, however, a small area of instability for metallicitiesas low as 10 − Z ⊙ . We have introduced a new method for including the ra-diative cooling from metals in large, three-dimensional hy-drodynamic simulations. In addition to its implementationin the AMR code, Enzo, this method has also been usedby Bogdanovi´c et al. (2006) in numerical simulations withthe smoothed particle hydrodynamics (SPH) code, Gadget(Springel et al. 2001; Springel 2005). Our technique takesadvantage of the extremely complex chemical reaction net-work of the preexisting radiative transfer code, Cloudy,which includes a full elemental coverage from H to Zn, alongwith a variety of molecular species and dust. With the sin-gular assumption of ionization equilibrium for the heavy el-ements, we are able to precalculate cooling rates for gaseswith any chemical abundance in all manners of radiation en-vironments over a temperature range of 10 to 10 K. Withcooling rates computed in advance, we eliminate the bar-rier that has classically prevented large chemical modelsfrom being incorporated into three-dimensional numericalsimulations. Because our cooling scheme is valid over sucha large range of density and temperature, and features somany coolants, it can be applied to a huge variety of as-trophysical problems, such as the evolution of the ISM andIGM, normal star formation, planetary nebulae, accretiondisks, and protoplanetary disks.One advantage of the large chemical network of Cloudyis that we are able to determine the dominant coolants froma complete sample of atomic species up to an atomic numberof 30. Fine-structure transitions of C and O are the greatestcontributors to the cooling up to number densities of about10 cm − , where Fe cooling becomes significant and C ismarginalized. The importance of these three elements, alongwith Si, in triggering the formation of the first low-massstars has been studied in great detail by Santoro & Shull(2006). The cooling models we have constructed using so-lar abundance patterns reveal S cooling to be important for 10 cm − ∼ < n H ∼ < cm − . S is produced only slightlyless than Si in type II (Woosley & Weaver 1995) and pair-instability supernovae (Heger & Woosley 2002), and shouldbe taken into account when considering the metals respon-sible for the transition from Population III to Population IIstar formation. The ability to specify individual abundancesin our cooling method makes it straightforward to simulatethe evolution of gas with non-solar abundance patterns.In addition to elevating the cooling rate of a gas tosatisfy the classical fragmentation criterion, metals also in-crease the potential for fragmentation by creating thermalinstabilities. We have identified regions in temperature anddensity in which both the classical fragmentation and ther-mal instability criteria are met to be the physical conditionsmost likely to see fragmentation occur. We observe thesedoubly unstable regions to exist for metallicities as low as10 − Z ⊙ . If fragmentation cannot occur outside these re-gions, then the fate of a star-forming gas-cloud will be deter-mined by the path taken through density-temperature spaceas it collapses. If we consider the doubly unstable regionsin Fig. 11, appropriate for star-formation in current epoch,there will almost certainly be a period of double instabilitywhen Z > − Z ⊙ . Interestingly, the density-temperaturetracks shown in Fig. 1 of Omukai et al. (2005) indicate thatgas with Z = 10 − Z ⊙ will pass right through the stableregion that separates two instabilities. Omukai et al. (2005)also find that star formation at that metallicity only pro-duces high-mass fragments. At high redshift, the CMB sig-nificantly reduces the size of the doubly unstable regions.Thermal instabilities, though, are extremely sensitive to theslope of the cooling rate as a function of density and tem-perature. Every element has distinct cooling properties, andwill, therefore, produce different thermal instabilities. Assuch, the key to uncovering the nature of the first PopulationII stars will be in the determination of the mass function oftheir precursors. In papers to follow, we will extend our study to gases withnon-solar abundance patterns. We will explore thermal anddouble instabilities created by individual elements, as well asabundance patterns produced in Population III supernovae.Future studies will also investigate the effects of backgroundradiation on the evolution of star-forming gas. The finalword, however, will only come from three-dimensional, nu-merical simulations. The simulations by Smith & Sigurdsson(2007), employing the methods described here with solarabundance patterns, have confirmed that fragmentation oc-curs for metallicities, Z > − Z ⊙ . We intend to pair allfuture predictions made from analysis of thermal instabili-ties with full numerical simulations.Although dust physics has been implemented in Cloudy,we currently treat only metals in the gas-phase in our anal-ysis. In the future, we will study the effects of dust coolingin more detail. We will also couple the dust chemistry to theH/He chemical network in Enzo, so as to properly model theformation of H on grain surfaces. A great strength of thismethod is the use of the ever-expanding chemical network ofCloudy. As more physical processes are incorporated into theCloudy software, the utility of this method will increase as c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation well. One major constraint of this work, however, is that itsvalidity is confined to the optically thin limit. The approx-imations made here break down at opacities of order unity.For higher opacities, our method provides a core modulefor flux-limited diffusion schemes. Complex geometries in-troduce problems of self-heating and shadowing, which willrequire full, three-dimensional radiative transfer. ACKNOWLEDGMENTS
Based on observations made with the NASA/ESA Hub-ble Space Telescope, obtained from the data archive at theSpace Telescope Institute. STScI is operated by the associ-ation of Universities for Research in Astronomy, Inc. underthe NASA contract NAS 5-26555. This research was sup-ported by Hubble Space Telescope Theory Grant HST-AR-10978.01, a grant from NASA’s ATP NNG04GU99G, andNSF CAREER award AST-0239709 from the National Sci-ence Foundation. We are grateful to an anonymous refereefor sharing many insightful comments and suggestions. Wealso thank Brian O’Shea and Mike Norman for useful dis-cussions and support. SS also thanks KIPAC and StanfordUniversity for their hospitality.
REFERENCES
Abel T., Anninos P., Zhang Y., Norman M. L., 1997, NewAstronomy, 2, 181Abel T., Bryan G. L., Norman M. L., 2002, Science, 295,93Anninos P., Zhang Y., Abel T., Norman M. L., 1997, NewAstronomy, 2, 209Anninos W. Y., Norman M. J., 1994, ApJ, 429, 434Black J. H., 1981, MNRAS, 197, 553Bogdanovi´c T., Smith B. D., Eracleous M., Sigurdsson S.,2006, in Laser Interferometer Space Antenna: 6th Inter-national LISA Symposium Vol. 873 of American Instituteof Physics Conference Series, Electromagnetic Signaturesof Massive Black Hole Binaries. pp 257–263Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23Bromm V., Ferrara A., Coppi P. S., Larson R. B., 2001,MNRAS, 328, 969Bromm V., Loeb A., 2003a, ApJ, 596, 34Bromm V., Loeb A., 2003b, Nature, 425, 812Bromm V., Loeb A., 2004, New Astronomy, 9, 353Bryan G., Norman M. L., 1997, in Chrisochoides N., ed.,Workshop on Structured Adaptive Mech Refinement GridMethods IMA Volumes in Mathematics No. 117, A Hy-brid AMR Application for Cosmology and Astrophysics.Springer-VerlagCen R., 1992, ApJS, 78, 341Clark P. C., Glover S. C. O., Klessen R. S., 2007, ArXivAstrophysics e-prints 0706.0613Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W.,Kingdon J. B., Verner E. M., 1998, PASP, 110, 761Field G. B., 1965, ApJ, 142, 531Galli D., Palla F., 1998, A&A, 335, 403Gao L., Yoshida N., Abel T., Frenk C. S., Jenkins A.,Springel V., 2007, MNRAS, 378, 449 Green D. A., Tuffs R. J., Popescu C. C., 2004, MNRAS,355, 1315Heger A., Woosley S. E., 2002, ApJ, 567, 532Hirata C. M., Padmanabhan N., 2006, MNRAS, 372, 1175Krause O., Birkmann S. M., Rieke G. H., Lemke D., KlaasU., Hines D. C., Gordon K. D., 2004, Nature, 432, 596Lepp S., Shull J. M., 1983, ApJ, 270, 578Low C., Lynden-Bell D., 1976, MNRAS, 176, 367Madau P., Ferrara A., Rees M. J., 2001, ApJ, 555, 92Maeder A., Meynet G., Hirschi R., 2005, in Barnes IIIT. G., Bash F. N., eds, Cosmic Abundances as Records ofStellar Evolution and Nucleosynthesis Vol. 336 of Astro-nomical Society of the Pacific Conference Series, ChemicalAbundances and Yields from Massive Stars. pp 79–+Nomoto K., Tominaga N., Umeda H., Kobayashi C., MaedaK., 2006, Nuclear Physics A, 777, 424Omukai K., 2000, ApJ, 534, 809Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005,ApJ, 626, 627O’Shea B. W., G. B., Bordner J., Norman M. L., Abel T.,Harknes R., Kritsuk A., 2004, in Plewa T., Linde T., WeirsV. G., eds, Adaptive Mesh Refinement - Theory and Ap-plications Vol. 41 of Lecture Notes in Computational Sci-ence and Engineering, Proceedings of the Chicago Work-shop on Adaptive Mesh RefinementO’Shea B. W., Norman M. L., 2007, ApJ, 654, 66Rees M. J., 1976, MNRAS, 176, 483Ricotti M., Ferrara A., Miniati F., 1997, ApJ, 485, 254Rockefeller G., Fryer C. L., Li H., 2006, ArXiv Astrophysicse-prints 0608028Santoro F., Shull J. M., 2006, ApJ, 643, 26Schneider R., Ferrara A., Natarajan P., Omukai K., 2002,ApJ, 571, 30Schneider R., Ferrara A., Salvaterra R., 2004, MNRAS,351, 1379Schneider R., Ferrara A., Salvaterra R., Omukai K.,Bromm V., 2003, Nature, 422, 869Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006,MNRAS, 369, 1437Smith B. D., Sigurdsson S., 2007, ApJ, 661, L5Spergel D. N. et al., 2007, ApJS, 170, 377Springel V., 2005, MNRAS, 364, 1105Springel V., Yoshida N., White S. D. M., 2001, New As-tronomy, 6, 79Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253Tan J. C., McKee C. F., 2004, ApJ, 603, 383Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T.,Palla F., 1997, ApJ, 474, 1Tsuribe T., Omukai K., 2006, ApJ, 642, L61van Hoof P. A. M., Weingartner J. C., Martin P. G., VolkK., Ferland G. J., 2004, MNRAS, 350, 1330Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003,ApJ, 592, 645Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ,652, 6This paper has been typeset from a TEX/ L A TEX file preparedby the author. c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel
Figure 1.
Cooling functions, excluding cooling from H , for gaseswith n H = 1 cm − and metallicities, Z = 0 (solid), 10 − Z ⊙ (dot),10 − Z ⊙ (dash), 10 − Z ⊙ (dot-dash), 1 Z ⊙ (dot-dot-dash), and10 Z ⊙ (dot-dash-dash). c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation Figure 2.
Cooling functions, including H cooling, for gases withn H = 1 cm − (top-left), n H = 10 cm − (top-right), n H = 10 cm − (bottom-left), and n H = 10 cm − (bottom-right). Metal-licities are Z = 0 (solid), 10 − Z ⊙ (dot), 10 − Z ⊙ (dash), 10 − Z ⊙ (long dash), 10 − Z ⊙ (dot-dash), and 10 − Z ⊙ (dot-long dash).c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel
10 100 1000T [K]10 -32 -30 -28 Λ [ e r g s - c m ] Total C I 369.7 µ mO I 145.5 µ m C II 157.6 µ mO I 63.2 µ m C I 609.2 µ mC I 985 nm Figure 3.
Cooling contributions from C and O species that reachat least 10 − of the total cooling for gas with with n H = 10 cm − and Z = 10 − Z ⊙ . The total cooling (solid, black) includesall species contained within the Cloudy chemical network. Com-ponents shown are [O I ] 145.5 µ m (dot-dot-dash), [O I ] 63.2 µ m(dash-dash-dot), [C I ] 369.7 µ m (dash), [C I ] 609.2 µ m (long dash),C I
985 nm (dash-dot), and [C II ] 157.6 µ m (long dash-dot). Thesolid, grey line indicates the sum of all the components plotted. c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation
10 100 1000T [K]10 -34 -32 -30 -28 Λ [ e r g s - c m ] Total Si I 129.6 µ mFe IIS I 25.2 µ mS I 56.3 µ m Si I 68.4 µ mSi II 34.8 µ m Figure 4.
All other coolants not plotted in Fig. 3 that reachat least 10 − of the total cooling for gas with n H = 10 cm − and Z = 10 − Z ⊙ . The total cooling (solid, black) includes allspecies contained within the Cloudy chemical network. Compo-nents shown are [Fe II ] (dot), [S I ] 25.2 µ m (dash), [S I ] 56.3 µ m(long dash), [Si I ] 129.6 µ m (dash-dot), [Si I ] 68.4 µ m (long dash-dot), and [Si II ] 34.8 µ m (dash-dot-dot). The solid, grey line indi-cates the sum of all the components plotted.c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel
10 100 1000T [K]10 -32 -30 -28 -26 Λ [ e r g s - c m ] Total C I 369.7 µ mCOO I 145.5 µ mO I 63.2 µ m C I 609.2 µ mC I 985 nm Figure 5.
Cooling contributions from C and O species that reachat least 10 − of the total cooling for gas with with n H = 10 cm − and Z = 1 Z ⊙ . The total cooling (solid, black) includesall species contained within the Cloudy chemical network. Com-ponents shown are CO (dot), [O I ] 145.5 µ m (dot-dot-dash), [O I ]63.2 µ m (dash-dash-dot), [C I ] 369.7 µ m (dash), [C I ] 609.2 µ m(long dash), and C I
985 nm (dash-dot). The higher dotted linerepresents CO emission from C O , while the lower dotted lineshows emission from C O . The solid, grey line indicates thesum of all the components plotted. c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation
10 100 1000T [K]10 -32 -30 -28 -26 Λ [ e r g s - c m ] Total S I 25.2 µ mFe IIFe I 24.0 µ mFe I 34.7 µ m S I 56.3 µ mSi I 129.6 µ mSi I 68.4 µ m Figure 6.
All other coolants not plotted in Fig. 5 that reach atleast 10 − of the total cooling for gas with n H = 10 cm − and Z= 1 Z ⊙ . The total cooling (solid, black) includes all species con-tained within the Cloudy chemical network. Components shownare [Fe II ] (dot), [Fe I ] 24.0 µ m (dash), [Fe I ] 34.7 µ m (long dash),[S I ] 25.2 µ m (dash-dot), [S I ] 56.3 µ m (long dash-dot), [Si I ] 129.6 µ m (dash-dot-dot), and [Si I ] 68.4 µ m (dash-dash-dot). The solid,grey line indicates the sum of all the components plotted.c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel
10 100 1000T [K]10 -34 -32 -30 -28 Λ [ e r g s - c m ] Total O I 63.2 µ mCOFe I 24.0 µ mFe II O I 630 nm Figure 7.
Subset of the most dominant coolants at n H = 10 cm − and Z = 1 Z ⊙ . The total cooling (solid) includes allspecies contained within the Cloudy chemical network. Compo-nents shown are CO (dot), [Fe I ] 24.0 µ m (dash), [Fe II ] (longdash), [O I ] 63.2 µ m (dash-dot), and O I
630 nm. c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation
10 100 1000T [K]10 -40 -38 -36 -34 -32 -30 Λ [ e r g s - c m ] Figure 8.
Total cooling rate from metals for gas at n H = 10 cm − , with metallicities Z = 10 − Z ⊙ (solid), 10 − Z ⊙ (dashed),and 10 − Z ⊙ (dash-dot). The black lines indicate the total cool-ing from gas-phase metals only. The grey lines show the coolingwith gas-phase metals and dust grains, created with the ISM dustgrain model in Cloudy, using the command, grains ISM . The ISMdust model in Cloudy consists of both graphite and silicates withsizes ranging from 5 × − µ m to 0.25 µ m and a power-law size-distribution with a power-law index of -3.5. For solar metallicity,the total grain abundances per H are 10 − . for graphite and10 − . for silicates. In each case shown, the dust grain abun-dances have been scaled to the gas-phase abundances.c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel log (n/cm −3 ) T [ K ] log (n/cm −3 ) T [ K ] log (n/cm −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] Figure 9.
Contours of log (t dyn /t cool ) over number density andtemperature for gases with metallicities, Z = 0 (top-left), 10 − Z ⊙ (top-right), 10 − Z ⊙ (middle-left), 10 − Z ⊙ (middle-right),10 − Z ⊙ (bottom-left), and 10 − Z ⊙ (bottom-right). H coolingis included. c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation log (n/cm −3 ) T [ K ] log (n/cm −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] Figure 10.
Contours of the instability parameter, ( α - β ) overnumber density and temperature. The medium is unstable forvalues less than 2. The metallicities are Z = 0 (top-left), 10 − Z ⊙ (top-right), 10 − Z ⊙ (middle-left), 10 − Z ⊙ (middle-right), 10 − Z ⊙ (bottom-left), and 10 − Z ⊙ (bottom-right). At Z = 10 − Z ⊙ ,two separate thermally unstable regions exist. These two regionsmerge by 10 − Z ⊙ .c (cid:13) , 1–19 Britton Smith, Steinn Sigurdsson, and Tom Abel log (n/cm −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] Figure 11.
The white patches indicate regimes of density andtemperature where log (t dyn /t cool ) > α - β ) < − Z ⊙ (top), 10 − Z ⊙ (middle), and 10 − Z ⊙ (bottom). As in Fig. 10, there are two individual doubly unstableregions at Z = 10 − Z ⊙ that have merged by 10 − Z ⊙ . c (cid:13) , 1–19 etal Cooling in Simulations of Cosmic Structure Formation log (n/cm −3 ) T [ K ] −3 ) T [ K ] −3 ) T [ K ] Figure 12.
Doubly unstable regions for the same metallicities asin Fig. 11, but with a CMB temperature floor at z = 20 included.c (cid:13)000