Thermal Phases of the Neutral Atomic Interstellar Medium -- from Solar Metallicity to Primordial Gas
DD RAFT VERSION J UNE
21, 2019
Preprint typeset using L A TEX style emulateapj v. 12/16/11
THERMAL PHASES OF THE NEUTRAL ATOMIC INTERSTELLAR MEDIUMFROM SOLAR METALLICITY TO PRIMORDIAL GAS S HMUEL B IALY (cid:63) & A MIEL S TERNBERG
Draft version June 21, 2019
ABSTRACTWe study the thermal structure of the neutral atomic (H I ) interstellar medium across a wide range of metallic-ities, from supersolar down to vanishing metallicity, and for varying UV intensities and cosmic-ray ionizationrates. We calculate self-consistently the gas temperature and species abundances (with a special focus on theresidual H ), assuming thermal and chemical steady-state. For solar metallicity, Z (cid:48) ≡
1, we recover the knownresult that there exists a pressure range over which the gas is multiphased, with the warm ( ∼ K, WNM)and cold ( ∼
100 K, CNM) phases coexisting at the same pressure. At a metallicity Z (cid:48) ≈ .
1, the CNM iscolder (compared to Z (cid:48) =
1) due to the reduced efficiency of photoelectric heating. For Z (cid:48) (cid:46) .
1, cosmic-rayionization becomes the dominant heating mechanism and the WNM-to-CNM transition shifts to ever increas-ing pressure/density as the metallicity is reduced. For metallicities Z (cid:48) (cid:46) .
01, H cooling becomes important,lowering the temperature of the WNM (down to ≈
600 K), and smoothing out the multiphase phenomenon.At vanishing metallicities, H heating becomes effective and the multiphase phenomenon disappears entirely.We derive analytic expressions for the critical densities for the warm-to-cold phase transition in the differentregimes, and the critical metallicities for H cooling and heating. We discuss potential implications on thestar-formation rates of galaxies and self-regulation theories. Subject headings: galaxies: ISM – galaxies: structure – galaxies: star formation – ISM: molecules – galaxies:high-redshift – early universe INTRODUCTIONEmission and absorption line studies of the 21 cm hyperfinespin-flip transition of atomic hydrogen (H I ) show that cold( T <
200 K) dense clumps are commonly embedded withinextended warm ( T ∼ I gasis also clearly observed in numerous hydrodynamical simula-tions of the interstellar medium (ISM) (e.g., Hu et al. 2016;Valdivia et al. 2015; Richings & Schaye 2016; Hill et al.2018).In the standard theory for the thermal structures of interstel- (cid:63) [email protected] Harvard-Smithsonian Center for Astrophysics, 60 Garden street, Cam-bridge, MA, USA School of Physics & Astronomy, Tel Aviv University, Ramat Aviv69978, Israel Center for Computational Astrophysics, Flatiron Institute, 162 5thAve., New York, NY, 10010, USA Max-Planck-Institut f¨ur extraterrestrische Physik (MPE), Giessen-bachstr., 85748 Garching, FRG lar H I , the WNM is cooled by electron impact Ly α emission(Gould & Thakur 1970; Dalgarno & McCray 1972), and thecold medium by metal line fine structure cooling, mainly neu-tral impact [C II ] 158 µ m emission (Dalgarno & Rudge 1964).For a given heating rate and at sufficiently low densities metalline cooling becomes very inefficient and the gas heats up totemperatures ( ∼ α emission limitsfurther heating. At high densities the metal line cooling be-comes effective and the temperature drops to values charac-teristic of the fine-structure energy level spacings ( ∼
100 K).Because of the particular shape of the combined Ly α plusfine-structure cooling function (we discuss this in detail), anintermediate range of densities and associated pressures existsat which the gas can be either warm or cold.Hydrodynamical simulations show that CNM forms at siteswhere flows of WNM collide and compress the gas (e.g.,Hennebelle & P´erault 1999; Kritsuk & Norman 2002; Audit& Hennebelle 2005; V´azquez-Semadeni et al. 2006; Gazol &Kim 2013; Saury et al. 2014; Gazol & Villagran 2016). Theformation of the dense CNM phase and the multiphase struc-ture of the ISM may play a key role in controlling the star-formation properties of galaxies. For example, Elmegreen& Parravano (1994); Schaye (2004) have argued that the lowpressure prevailing at the outskirts of galaxy discs do not al-low the formation of the cold phase resulting in a shut-offof star-formation. The multiphase phenomena where CNMand WNM phases coexist at the same pressure may providea feedback loop for star-formation in galaxies (Corbelli &Salpeter 1988; Parravano 1988, 1989; Ostriker et al. 2010;Kim et al. 2011; Ostriker & Shetty 2011). If the star-formationrate drastically increases, the high heating rate induced by thestellar UV prevents the formation of the CNM phase, whichsuppresses star-formation.Previous computations of the H I WNM/CNM thermalphase structures have generally assumed that metals domi- a r X i v : . [ a s t r o - ph . GA ] J un Bialy & Sternbergnate the CNM cooling. Then, as would be expected, as themetal abundance is reduced higher particle densities are re-quired to compensate for the reduced cooling, and the mul-tiphase behavior is predicted to occur at higher gas densitiesand thermal pressures (Wolfire et al. 1995; Liszt 2002; Walchet al. 2011). The central question we address in this paperis how WNM/CNM multiphase behavior of atomic hydrogen(H I ) gas is affected by H cooling and heating processes asthe metallicity becomes vanishingly small.H cooling in the context of primordial and/or low-metallicity gas and first star-formation has been considered bymany authors over many years (e.g., Saslaw & Zipoy 1967;Lepp & Shull 1984; Shapiro & Kang 1987; Haiman et al.1996; Tegmark et al. 1997; Omukai 2000; Abel et al. 2000;Yoshida et al. 2006; Glover & Jappsen 2007; Schneider et al.2012; Wolcott-Green et al. 2017). It has long been known thatin the absence of metals, just small amounts of H can enablecooling to a few 100 K (see for example the recent study byGlover & Clark 2014) thereby controlling the fragmentationand dynamical properties of collapsing star-forming clouds inthe low-metallicity regime. Less attention has been given tothe specifics of WNM to CNM transitions and possible mul-tiphase behavior at very low-metallicities. For example, intheir study of protogalactic disks Norman & Spaans (1997)presented pressure versus density thermal phase diagrams foratomic gas over a range of redshifts and associated metallic-ities. In their models multiphase behavior disappears at highredshifts (see their Fig. 5). However it remains unclear whatare the individual roles of the metallicity and radiation fieldwhen H processes begin playing a role, and how these pro-cesses affect the phase structure. Similarly, Inoue & Omukai(2015) carried out numerical hydrodynamical simulations ofisobarically contracting gas clouds and studied the effects ofH cooling on thermal instabilities. They found that at suffi-ciently low metallicities H acts to smooth out thermal insta-bilities (see their Fig. 2). But why?The goal of our paper is to study the competing microphys-ical processes that control the phase structures of atomic gas,into the regime of vanishing metallicity. We address the fol-lowing questions. How are the WNM and CNM phases af-fected by the onset of H cooling (and heating)? Is multi-phase gas possible when H replaces metal fine-structure linecooling? If not, why not? In this study we present numeri-cal computations of the (equilibrium) thermal states, with UVradiation intensity, cosmic-ray (or X-ray) ionization rate, andmetallicity, as distinct parameters, and with the addition ofthe H processes in a step by step fashion. This enables ananalytic description of the behavior which we also present.The structure of our paper is as follows. In § §
3. We discuss low metallicity be-havior in § dominatedgas. We characterize the critical densities or pressures abovewhich the gas may cool below ∼ K and condense into in-termediate ( ∼
500 K) or cold ( ∼
50 K) phase (depending onmetallicity), and whether two-phase media can exist, as func-tion of metallicity, UV intensity and CR ionization rate. Wesummarize in § MODEL INGREDIENTSWe compute the temperatures and gas phase chemical com-positions of the neutral atomic hydrogen (H I ) components ofthe interstellar medium (ISM) of galaxies assuming thermal equilibrium and chemical steady state. We use the results toconstruct pressure versus density phase diagrams for the H I gas. We consider standard heating and cooling mechanisms,with the addition of molecular hydrogen (H ) processes thatbecome important and eventually dominant in the limit of lowmetallicity. We investigate the role of the assumed far-UV ra-diation intensities and cosmic ray-ionization rates on the cou-pled thermal and chemical properties of the H I gas. We studythe behavior of the thermal phase structures for standard ISMconditions into the limit of vanishing metallicity.2.1. Basic Equations
In thermal equilibrium, L ( n , T ; x i ) = G ( n , T ; x i ) , (1)where G and L are the total heating and cooling rates per unitvolume (erg cm − s − ). The rates depend on the total hydro-gen volume density n (cm − ), the gas temperature T (K), thedust-to-gas abundance ratio, and the gas phase abundances ofthe chemical species, x i ≡ n i / n , especially for C + , O, elec-trons, and H . For heating we include far-UV photoelectric(PE) ejection of electrons from dust grains and polycyclic aro-matic hydrocarbons (PAHs), and cosmic-ray (CR) and/or X-ray ionization of the hydrogen atoms. For cooling, we includeH I Ly α line-emission excited by electron impacts, and fine-structure [C II ] (157.7 µ m), [O I ] (63.2 µ m, 145.5 µ m), and [C I ](609.1 µ m, 370.4 µ m) line emissions excited primarily by col-lisions with the neutral hydrogen atoms. We also include en-ergy losses due to electron recombination with dust grains andPAHs, and due to gas-dust collisions. As discussed by Wolfireet al. (2003, hereafter, WMHT) these are the dominant heat-ing and cooling mechanisms for thermally stable and unsta-ble warm/cold ( ∼ to ∼
100 K) interstellar H I for standardISM conditions .Importantly, we also incorporate heating and cooling pro-cesses involving H . This includes heating via the release ofbinding energy during H formation, kinetic energy impartedto the H-atoms produced by H photodissociation, and inter-nal H energy released to the gas via collisional deexcitationof UV excited states (we refer to this latter process as ”UVpumping heating”). Cooling is mainly by H ro-vibrationalline-emissions excited by collisions with electrons and neu-trals. We also include cooling by isotopic HD line emissions,and due to dissociation of the H in collisions with thermalelectrons and H-atoms, although these are minor. In the Ap-pendix we present detailed discussions and expressions for thevarious heating and cooling processes we have included in ourcomputations. The roles and effects of H heating and cool-ing processes have not been considered in previous analyticstudies of the multiphase structure of the neutral ISM (Fieldet al. 1969; Wolfire et al. 1995; Liszt 2002; WMHT). We willshow how the H heating and cooling processes influence thethermal structures as the heavy element abundances becomesmall.Eq. (1) is coupled to the equilibrium rate equations for thechemical abundances, ∑ j > l x j x l nk jl , i + ∑ j x j q j , i = . (2)This is a set of N equations, equal to the number of speciesin the chemical network. The first sum is over the formation Metastable line emission (e.g. [OI] 6302 ˚A) also cools the WNM, but itscontribution is subdominant (WMHT, see also Fig 30.1 in Draine 2011). hermal Phases – from Solar Metallicity to Primordial Gas 3and destruction rates (s − ) of species i via two-body reactionsinvolving species j and l , where k jl , i are the temperature-dependent reaction rate coefficients (cm s − ). The secondsum is over the formation and destruction rates of species i due to ionization or dissociation of species j by either cosmic-rays or UV photons, and the q j , i are the rates (s − ) for theseprocesses. For j = i , the coefficients k jl , i and q j , i are neg-ative and represent the destruction of species i in reactionswith species l , and by UV and CR dissociations and ioniza-tions. For j (cid:54) = i the coefficients are positive and represent theformation channels for species i . Eqs. (2) are augmented bythe (dependent) conservation equations for the total elementalabundances and electric charge. These are X m = ∑ i α i , m x i , (3)where X m is the total elemental abundance of element m , and α i , m is the number of atoms of element m in species i . Eq. (3)also describes the electric charge conservation equation, with α i , m representing the number of positive charges in species i , and with X m =
0. In total, there are M + M is the number of elements in the chemicalnetwork. We discuss our chemical network in § Parameters
The gas temperature T , and species abundances x i , dependon the gas density, n , the far-UV radiation intensity, I UV , theCR ionization rate, ζ cr , and the metallicity Z (cid:48) . We define theseparameters below. Additional physical parameters affectingthe thermal and chemical solutions are listed in Table 1.2.2.1. Metallicity and dust-to-gas ratio
The metallicity, i.e. the gas phase abundances of the heavyelements, and the dust abundance, are crucial parameters forthe thermal balance. The gas phase carbon and oxygen abun-dances determine the efficiency of metal-line cooling, and thedust-to-gas mass ratio determines the efficiency of PE heat-ing and dust recombination cooling. We assume that the total gas-phase abundances of the carbon and oxygen species aregiven by X m = A m Z (cid:48) (cid:0) − δ m × ( Z (cid:48) d / Z (cid:48) ) (cid:1) , (4)where m = C, or O. In these expressions Z (cid:48) is the normalizedmetallicity such that Z (cid:48) = A O = . × − and A C = . × − (Asplund et al. 2009), δ C and δ O are dust depletion factors,and Z (cid:48) d is the normalized dust-to-gas ratio. We adopt Z (cid:48) d = (cid:26) Z (cid:48) for Z (cid:48) ≥ Z (cid:48) Z (cid:48) ( Z (cid:48) / Z (cid:48) ) α for Z (cid:48) < Z (cid:48) , (5)with Z (cid:48) = .
2, and α =
3. This broken power-law relation issuggested by recent observations of low metallicity galaxies(R´emy-Ruyer et al. 2013). Here, Z (cid:48) d = . × − (Zubko et al.2003, Table 6, the BARE-GR-S model).For the dust depletion factors, we adopt δ C = .
53 and δ O = .
41. With these values, we recover the characteristicinterstellar Galactic carbon and oxygen gas-phase abundances X C = . × − , X O = . × − at Z (cid:48) = Z (cid:48) d = Z (cid:48) d / Z (cid:48) (cid:28) .
2, the dustdepletion terms vanish, and all the available carbon and oxy-gen remain in the gas-phase, with X C = A C Z (cid:48) , X O = A O Z (cid:48) . 2.2.2. Cosmic-ray and X-ray ionization
Cosmic rays (or X-rays) ionize atomic and molecular hy-drogen, H + cr → H + + e (6)H + cr → H + + e (7)H + cr → H + + H + e (8)producing energetic non-thermal electrons. The electronsgenerated by these primary ionizations lose energy throughcoulomb scatterings that heat the gas, and through furtheratomic and molecular excitations and secondary ionizations(Dalgarno et al. 1999). The total (primary+secondary) ion-ization rate per hydrogen atom is ζ cr = ζ p × ( + φ s ) , where ζ p is the primary ionization rate and φ s ∼ . ζ cr ≈ − s − (Indriolo & Mc-Call 2012; Tielens 2013), and we define the normalized CRionization rate, ζ − ≡ ζ cr / ( − s − ) . Our fiducial value is ζ − =
1. We study the effects of varying ζ − in § ζ cr and φ s .2.2.3. Interstellar far-UV radiation field
The far-UV radiation intensity determines the efficiency ofphotoelectric heating, the H photodissociation rate, and theproduction rate of the C + coolant. Draine (1978, 2011) esti-mated the local UV interstellar radiation field (ISRF), and ob-tained u UV , = .
056 eV cm − , for the energy density over the6-13.6 eV band (=1.69 times the Habing 1968 estimate). Weadopt the Draine spectral shape and allow for variations in theUV intensity by multiplying the overall spectrum with a nor-malized intensity parameter, I UV , such that u UV = u UV , I UV .The UV intensity determines the PE heating rate throughEq. (A3). For a unit Draine ISRF, the free-space (unatten-uated) H photodissociation rate is D = . × − s − (Sternberg et al. 2014). At low metallicity, a harder spectrummay be more typical, since the first generations of stars (pop-III stars) are theorized to be very massive, with stellar effec-tive temperatures as high as 10 K (e.g., Cojazzi et al. 2000;Barkana & Loeb 2001; Abel et al. 2002; Yoshida et al. 2012).Longer wavelength photons, down to 0.75 eV, play a role inremoving the H − anion via photodetachment (Miyake et al.2010). H − is a critical intermediary for H formation at lowmetallicity (see § − photode-tachment rate is D − = . × − s − . In contrast, for a 10 K blackbody (BB) spectrum (normalized such that the photondensity in the 6-13.6 eV band equals that of the Draine field), D − = . × − s − (Bialy & Sternberg 2015; see their Table2). For our fiducial models we seek a compromise betweenthese two extreme cases and extrapolate the Draine (1978)spectrum (see their Eq. 11) down to 0.75 eV. The resultingH − photodetachment rate is D − = . × − s − . This value Bialy & Sternberg TABLE 1P
ARAMETERS notation fiducial value / range meaning comments Z (cid:48) from 10 − to 3 metallicity relative to Galactic see § § I UV § § ζ − = ζ cr − s − § § A C , A O × − the carbon and oxygen solar elemental abundances see § δ C , δ O § Z (cid:48) d the DGR relative to Galactic see § ζ p ζ cr / ( + φ s ) the primary ionization rate per H atom see § φ s ≈ . § D . × − s − the H photodissociation rate at I UV = § § D − . × − s − the H − photodetachment rate at I UV = § § R (cm s − ) the H formation rate coefficient see § x i = n i / n the species abundance, relative to the hydrogen nucleon density see 2.1 and 2.2.4 L , G (erg cm − s − ) the cooling and heating rates per unit volume see § L i ( T ) (erg cm s − ) the cooling efficiency of species i . For n / n crit (cid:28) L i = x i n L i ( T ) see § A.2 Λ i ( T ) (erg s − ) the cooling rate per particle. For n / n crit (cid:29) L i = x i n Λ i ( T ) see § A.2 is close to the 10 K BB value which represents the pristinemetallicity gas, but is ∼ D − may range between these two extremes, and will depend onthe form of the initial mass function (IMF). We explore theeffect of variations in I UV and D − in § D may be stronglyreduced in the interiors of a cloud due to dust absorption andH self-shielding, leading to partial or complete H I -to-H conversion (van Dishoeck & Black 1986; Draine & Bertoldi1996; McKee & Krumholz 2010; Sternberg et al. 2014; Bialy& Sternberg 2016). In this paper we focus on the thermalproperties of the global atomic ISM, rather than the structureof the H I -to-H transition or the interior shielded moleculargas. Thus, we always use the free-space UV field and do notinclude radiation attenuation in our models. For the WNMshielding is in any case not significant (except for extremelylarge shielding columns). For the CNM, our model applies tothe cloud outskirts that are exposed to the free-field UV.2.2.4. H Formation
At sufficiently high metallicity, the H is formed predomi-nantly on dust grains. We adopt R d = × − (cid:18) T
100 K (cid:19) / Z (cid:48) d cm s − . (9)Sophisticated models of the H formation process on dust-grains (e.g., Cazaux & Tielens 2002; Le Bourlot et al. 2012;Bron et al. 2014) derived more complicated dependences ongas temperature, as well as on other physical parameters (e.g.,through the dust temperature which in turn depends on the UVradiation intensity). However for our purposes such a simpledescription is sufficient, because as we show below, H af-fects the thermal structure only at low metallicities where itsformation is dominated through gas-phase processes.At low metallicities ( Z (cid:48) (cid:46) . is formed mainlythrough gas-phase reactions, with the most important onebeing the H − formation route (McDowell 1961; Peebles &Dicke 1968; Hirasawa 1969; Galli & Palla 1998). Although,the H formation rate is obtained through the numerical so-lution to the chemical network (as discussed in § fractional abundance.This formation channel is initiated by radiative attachmentH + e → H − + ν , (10)followed by associative detachmentH − + H → H + e . (11)Assuming that the electron abundance is set by a balancebetween CR ionization of the atomic hydrogen and proton-electron recombination, we have x e = (cid:18) ζα B n (cid:19) / = . × − T / ζ / − (cid:16) n
10 cm − (cid:17) − / , (12)where T ≡ T / ( K ) and α B is the case-B recombinationrate coefficient. The effective H formation rate coefficientthrough this channel is R − ≡ x e k η ≈ . × − T . η ζ / − (cid:16) n
10 cm − (cid:17) − / cm s − (13)where k is the radiative attachment rate coefficient (Eq. 10).The H formation sequence is moderated mainly by photode-tachment, H − + ν → H + e . (14)and also by mutual neutralization with protons,H − + H + → H + H . (15)In Eq. (13) we define the branching ratio η = (cid:18) + D − I UV k n + k x H + k (cid:19) − (16)for H formation via associative detachment versus H − pho-todetachment and mutual-neutralizaion. For D − = . × − s − photodetachmnet becomes important when n / I UV (cid:46) − , while mutual-neutralization is important when The rate coefficients are based on UMIST 2012 (McElroy et al. 2013).See Table 2 for a summary of the most important reactions and rates. Seeonline material for the full network. hermal Phases – from Solar Metallicity to Primordial Gas 5 n / ζ − (cid:46) .
01 cm − . As we show below, H plays an im-portant role in cooling and heating at densities (cid:29) − .Thus, generally, η ≈ .Assuming that H destruction is dominated by photodisso-ciation, we obtain an analytic approximation for the H frac-tional abundance at low metallicities, x H = x e k nD I UV η = . × − T . ηζ / − I − (cid:16) n
10 cm − (cid:17) / . (17)We find Eq. (17) to be in good agreement with the numericalsolution for the chemical network, over the gas densities andtemperatures in our parameter space.2.3. Numerical Method
Given a set of parameters, I UV , ζ cr , and Z (cid:48) , we solve for theequilibrium temperature and species abundances as functionsof the hydrogen density, n . We form a logarithmic grid, span-ning log T = . n = − ×
90 sam-pling points. For each volume density in the grid, we calculatethe species abundances coupled with the net cooling, L − G , asfunctions of the temperature, and find the temperature(s) forwhich the thermal equilibrium condition (Eq. 1) is satisfied.This gives us the T ( n ) curves shown in Figs. 1 and 4 below.We compute the equilibrium thermal pressure (cm − K) bysimply multiplying by the density, i.e., P ( n ) = nT ( n ) . Forany density, heating exceeds cooling for pressures less than P ( n ) , and cooling exceeds heating for pressures greater than P ( n ) .For our chemical network (Eq. 2), we consider a set of 294two-body gas-phase reactions among 34 species, composedof hydrogen, helium, carbon, and oxygen bearing atomic andmolecular species. We adopt the rate coefficients suggestedby the UMIST 2012 database (McElroy et al. 2013), exceptfor the H collisional dissociation for which we adopt theMartin et al. (1996) rate (including both collisional dissoci-ation and dissociative tunneling). We also include the forma-tion of H on dust grains, CR reactions and photo-reactions.The full reaction list is provided in the supplement materials(online version). The dominant photorates are also listed inBialy & Sternberg 2015, see their table 2.In total we obtain a set of 34 coupled non-linear algebraicequations for the 34 species. The equations depend on T and n , as well as on the additional physical parameters, I UV , ζ cr , and Z (cid:48) . We replace five of these equations with the ele-ment conservation (H, He, C, and O), and charge conservationequations. For any density n there is always a single uniquesolution for the temperature T and abundances x i . Rather thansolving the equilibrium algebraic equations (e.g., via Newton-Raphson iteration), we instead integrate the time-dependentrate equations until a time, t , where dx i / dt = THERMAL STRUCTURE AT SOLAR METALLICITY For high D − values, η may become smaller than unity at the coolingpoint, reducing the H abundance. We discuss this in § In this section we first present results for our solar metallic-ity reference models ( Z (cid:48) = Z (cid:48) d = n , assuming I UV = ζ − =
1. The dashedcurves in the left panels, located above and below the solidcurve, are for I UV = ζ − =
10 and 0.1, respectively.First we focus on the behavior at low to intermediate densi-ties, n (cid:46) cm − , where H cooling and heating processesare subdominant. In this regime, our results essentially re-produce those presented by Wolfire et al. (1995) and WMHT.The behavior is controlled by the temperature dependence ofthe Ly α and fine-structure metal-line cooling functions. Wedisplay these cooling functions in the upper panel of Fig. 3.The black curve shows the exponentially steep electron im-pact excitation Ly α cooling rate (erg s − ) per hydrogen nu-cleus. The blue curves show the much flatter metal coolingrates for Z (cid:48) =
1, 10 − , and 10 − , assuming all of the gas-phasecarbon is C + , and the oxygen is in neutral atomic form. Themetal cooling is controlled by electron and neutral-hydrogenimpact excitation of the fine-structure levels. The solid anddashed curves in Fig. 3 are for x e = − and 10 − respec-tively. For the high x e , and when T (cid:46)
500 K, metal line cool-ing is dominated by electron impact excitation. The shapes ofthese curves are of the utmost importance.Cooling becomes inefficient compared to heating as thedensity becomes small. For n (cid:46) . − , the gas becomeswarm, and Ly α dominates the cooling. In this regime thetemperature depends weakly on density and remains close to10 K. This is because for T (cid:46) K Ly α cooling is expo-nentially sensitive to the temperature (see Fig. 3), and anyincrease in the cooling rate due to an increase in gas den-sity requires just a small decrease in temperature to balancethe heating rate. At n ≈ − , the temperature decreasesto ∼ × K, at which point energy losses by the [C II ]and [O I ] emissions become comparable to Ly α , as seen inFig. 1 and Fig. 3. Above this density the metals dominatethe cooling and the temperature falls abruptly. The sharpdrop occurs because at these high temperatures (and as longas T (cid:38)
300 K, i.e. above the transition energies of the fine-structure lines) metal cooling is a weak function of temper-ature (Fig. 3). Any increase in density then requires a largedecrease in temperature to balance the heating rate. We referto the density at which metal cooling takes over as the “cool-ing point density”, and we provide analytic estimates in § T ( n ) curve flattens again.In Fig. 1, the low-density flat part of the T ( n ) curve fol-lowed by the sharp temperature drop gives rise to a local max-imum, P max , in the P ( n ) curve. In general, a local maximumsuch as this occurs if and where the temperature starts de-creasing more rapidly than 1 / n . For our fiducial model with I UV = ζ − =
1, we find P max = . × cm − K occurringat a density n max = . − . For n < n max the gas is warmwith T WNM = K. This is the classical warm neutralmedium (WNM), and P max is the maximum pressure possiblefor the WNM. Starting from cold gas at high densities, thesteep rise of the T ( n ) curve as the density is reduced givesrise to a local minimum, P min , in the pressure versus densitycurve. In general, a local minimum occurs if and where the Bialy & Sternberg F IG . 1.— The temperature, pressure, and cooling and heating rates per hydrogen nucleus as functions of the volume density, assuming solar metallicity and I UV = ζ − =
1. The thin diagonal (horizontal) lines in the left panels indicate the minimum and maximum pressures for a multiphase medium. The dashed redand blue curves in the left panels show models with I UV = ζ − = . temperature starts rising more rapidly than 1 / n as the densityis reduced. For our fiducial model P min = . × cm − Kat n min = . − . For n > n min , the gas remains cold with T CNM (cid:46)
300 K. This is the cold neutral medium (CNM) and P min is the minimum pressure possible for the CNM. We show P min and P max in lower (upper) panels of Fig. 1 as the pair ofthin horizontal (diagonal) lines.Our Fig. 1 reproduces the well-known result (Field et al.1969; Wolfire et al. 1995) that in the transition from Ly α tometal-line cooling the temperature falls more rapidly than 1 / n over a significant density interval from n max to n min . Thisenables multiple solutions for the gas temperature between P min and P max for densities ranging from n min T CNM / T WNM to n max T WNM / T CNM . For densities between n max and n min thegas is thermally unstable and this is the UNM. For isobaricdensity perturbations δ n in the UNM, the gas either cools tothe stable CNM branch (for positive perturbations), or heatsto the stable WNM branch (for negative perturbations). Forthermal pressures between P max and P min the gas may becomemultiphased with CNM condensations coexisting within anenveloping WNM at the same thermal pressure.Fig. 1 shows that for n (cid:38) × cm − , the temperature in-creases with density, as H heating processes start contribut-ing. In this regime the H heating rate increases more steeplythan n (see Eq. A15 in the limit n / n crit <
1) leading to thetemperature increase with n . We stress that the H heatingprocesses play this role so long as the gas is primarily atomic.If any self-shielding induces an HI-to-H transition, the H heating will be strongly suppressed. For solar metallicity atotal (atomic plus molecular) column as small as N ∼ cm − is sufficient to induce conversion to molecular form for n (cid:38) cm − (and I UV = N ∼ cm − , are required (Bialy &Sternberg 2016)At still higher densities (beyond the considered parameterspace) the hydrogen undergoes a conversion into H and theatomic carbon and oxygen convert into CO and H O. Thesemolecular formation processes become efficient as the tem-perature exceeds few hundred Kelvins, which triggers theneutral-neutral chemical network of H O (Sternberg & Dal-garno 1995; van Dishoeck et al. 2013; Bialy et al. 2015). Asour focus is on the (predominantly) atomic ISM, we restrictourselves to densities n (cid:46) cm − .The dashed curves in the left panel of Fig. 1 are theequilibria for I UV = . ζ − / I UV =
1. For all values of I UV in Fig. 1 photoelectricheating dominates and is proportional to I UV . Thus, when I UV is increased (decreased) by a factor of 10, the T ( n ) curvesshift to the right (left), by a factor of 10. As a result, the pres-sure range within which the gas may be multiphased shiftsto higher or lower values in proportion to I UV . For exam-ple, for I UV = ( P max , P min ) = ( . , . ) × cm − K. For I UV = . ( P max , P min ) = ( . , . ) × cm − K.At solar metallicity, H cooling and heating become impor-tant only in the very dense CNM, but as we discuss below,at lower metallicities, H cooling starts to play a role alsohermal Phases – from Solar Metallicity to Primordial Gas 7 n (cm − ) ζ − CRheating PEheating -1 Z ′ -2 -1 ζ − / I UV F IG . 2.— Contours at which CR heating equals PE heating in the I UV / ζ − − Z (cid:48) parameter space. The solid, dashed and dotted curves arefor T = , and 10 K respectively, and the different colors corre-spond to different values of n / ζ − , from n / ζ − =
10 (left most curves)to n / ζ − = . at low-to-intermediate densities, in the WNM, affecting theWNM-to-CNM transition. THERMAL STRUCTURE AT LOW METALLICITYIn this section we investigate the thermal structures as themetallicity is lowered. We first consider the transition fromphotoelectric to cosmic-ray dominated heating as Z (cid:48) is re-duced, for a range of ratios ζ − / I UV . We then solve the equa-tion of thermal equilibrium, in two steps. First, we lower themetallicity and solve for T with the exclusion of the H heat-ing and cooling processes. Then we solve again with the H included. This procedure enables us to distinguish the effectsof diminished metal line cooling from the appearance of theH heating and cooling processes that eventually dominate thebehavior at very low metallicity.4.1. PE versus CR
Near solar metallicity ( Z (cid:48) ∼
1) photoelectric (PE) heatingusually dominates, and the PE heating rate decreases in pro-portion to the dust-to-gas ratio Z (cid:48) d . At sufficiently low Z (cid:48) heat-ing by cosmic-ray ionization becomes competitive and domi-nates. In Fig. 2 we show the boundary lines at which the PEand CR heating rates are just equal, in the I UV / ζ − versus Z (cid:48) plane, for n / ζ − ratios equal to 0.1, 1, and 10 cm − . For ζ − / I UV = Z (cid:48) = . ζ − / I UV , the transition metallicity islowered, but not by much. This is because below Z (cid:48) = . ζ − / I UV ( (cid:38) Step-1: Reducing the Metallicity
We now consider the thermal structures at reduced metallic-ity, with the exclusion of the H heating and cooling processesin Step-1 . Results for T ( n ) and P ( n ) are shown as the dashedcurves in Fig. 4 (top and middle panels) for Z (cid:48) from 3 to 10 − .In these computations we have assumed I UV = ζ − = Z (cid:48) is lowered to 0 .
1, ( P max , P min ) dropto ( . , . ) × cm − K. These drops are due to the more rapid reduction in Z (cid:48) d and associated PE heating rate, com-pared to the diminished metal line cooling rate that decreaseswith Z (cid:48) . The reduced heating relative to cooling at Z (cid:48) = . Z (cid:48) = n min = . − . The maximal densityof the WNM is also lower, with n max = .
71 cm − . As themetallicty is reduced further CR heating takes over, and thegas then remains warm to ever higher densities. This is ex-pected given the reduction in the cooling efficiencies for afixed CR heating rate that is independent of metallicity. Im-portantly, multiphase behavior is maintained, with P min and P max shifted to higher pressures.The behavior in the CR heating regime may be analyzedstraightforwardly as follows. As we discussed in §
3, as thedensity is increased in the warm phase, metal line-emissionenergy losses start contributing at a temperature of ∼ × K. The cooling point temperature is nearly independent of Z (cid:48) , and we denote the corresponding density as n cool , Z . For n > n cool , Z , the temperature drops sharply as Ly α becomesinefficient and metal cooling dominates. To obtain an ana-lytic formula for n cool , Z we equate the heating and metal linecooling rates at T ≡ × K. We assume heating by CRionization. This gives n cool , Z = ζ p E cr A O L Z ( T ) Z (cid:48) ≈ . (cid:18) − Z (cid:48) (cid:19) ζ − cm − , (18)where we have defined the Solar metallicity fine-structureline cooling efficiency (erg s − cm ) per oxygen atom L Z ≡ L [ OI ] + ( x C + / x O ) L [ CII ] . In Eq. (18) ζ p is the primary ion-ization rate, and E cr is the mean deposition energy per ion-ization (see Eq. A8). In the numerical evaluation we used L Z ( T ) = . × − erg s − cm , x C + / x O = A C / A O = . A O = . × − (assuming no depletion), and E cr = . φ s = .
57. These values are appropriate for x e = . × − ,the equilibrium electron fraction at the transition point at Z (cid:48) = .
01. For Z (cid:48) = − Eq. (18) gives, n cool , Z = . − ,decreasing linearly with Z (cid:48) , in good agreement with the tem-perature drop points shown for T ( n ) in Fig. 4 (dashed curves).For Z (cid:48) (cid:38) . processes are neglected (or negligible) thesharp drop in temperature induced by the transition fromLy α to metal line cooling gives rise to the multiphased P ( n ) curves, and then n cool , Z is approximately equal to n max themaximum density for which a WNM is possible. Thus, P max (cid:39) . n cool , Z T = . × (cid:18) − Z (cid:48) (cid:19) ζ − cm − K . (19)We assume a cosmological helium abundance, A He = . is included this simplepicture is greatly altered.4.3. Step-2: Inclusion of H Processes Here n cool , Z is the transition point for Ly α -to-metal cooling. Later wealso introduce n cool , H for the transition point from Ly α -to-H cooling. As shown in Figs. 1 and 5, [CII] always dominates over [CI] cooling.
Bialy & Sternberg F IG . 3.— Top: The various components of the cooling function per hy-drogen nucleus, L / n (erg s − ), Ly α , metal ([CII], [OI]), and H cooling, asfunctions of temperature. We assume n = − , x H + = x e = − (solid) or10 − (dashed), and x H = − . Typically, for T (cid:28) K the density is highand x e is small. For Ly α cooling we always assume x e = − . The threesets of metal cooling curves are for Z (cid:48) = , − , and 10 − . Bottom: The H cooling efficiency L H (erg cm s − ) by atomic hydrogen excitation (basedon data from Lique 2015). Recall, for low densities, L H = x H n L H . Thecontributions from decays of individual ro-vibrationnal H levels, ( v , J ) , andtheir energies above ground (in Kelvin) are indicated. The red curves in the upper panel of Fig. 3 show the H cooling function. In the lower panel we show the contributionof individual H rotational and vibrational levels to H cool-ing, for excitation by atomic hydrogen. As we discuss below,the fact that (unlike metal cooling) the H cooling functionnever saturates with temperature all the way until Ly α be-comes operational, smooths out of the multiphase phenomenaand eventually eliminates it entirely at low metallicity.The solid curves in Fig. 4 show our results for the equi-librium T ( n ) and P ( n ) with the H processes included in our Step-2 . Again, these are heating via H formation, photodis-sociation, and UV-pumping, and cooling via ro-vibrationalline emissions, and collisional dissociation. The various heat-ing and cooling rates for Z (cid:48) = − , − ..., − , are shown inFig. 5.The efficiencies of the molecular processes depend on theH abundance, x H , that we compute self-consistently with thetemperature, for a given gas density and metallicity. In Fig. 4we show x H versus n , for Z (cid:48) from 3 to 10 − . At all metallici-ties, H destruction by photodissociation is very efficient andthe gas remains predominantly atomic. The H abundance ishigher at high metallicities as the H is formed efficiently ondust-grains. For Z (cid:48) (cid:46) .
1, the H is formed much less effi-ciently in the gas phase. The resulting H abundances then F IG . 4.— Temperature, thermal pressure, H abundance, and the C, C + andCO abundances relative to carbon nuclei abundance, as functions of the gashydrogen density, for various metallicity values, from Z (cid:48) = − .The dashed curves in the top panels are for models that exclude H coolingand heating processes. range from 10 − to 10 − , consistent with Eq. (17). Interest-ingly, it is in the low metallicity limit where x H is very low,hermal Phases – from Solar Metallicity to Primordial Gas 9that H plays the most significant role in heating and coolingthe gas. At high Z (cid:48) molecular hydrogen is more abundant butmetal line cooling is very efficient and induces a WNM-to-CNM transition at densities well below the densities at whichH cooling starts playing a role. However, at sufficiently low Z (cid:48) , H cooling must eventually dominate over metals.4.3.1. The onset of H cooling at low metallicity When the H processes are included the T ( n ) and P ( n ) curves in Fig. 4 are altered in several ways. First, at suf-ficiently low metallicity it is H rather than the metals thatcompete with Ly α at the “cooling point” temperature of ∼ × K below which Ly α becomes inefficient. Thisis also illustrated in the left panels of Fig. 5. We define Z (cid:48) cool , H as the critical metallicity at which this switch frommetals to H occurs at the cooling point. Numerically we find Z (cid:48) cool , H ≈ − . The cooling point density is then shifted to alow value of ≈ − , independent of metallicity (see Figs. 4and 5). The H abundance at this point is x H ≈ − .For an analytic estimate we define the gas density for theH cooling point where n L H x H = ζ p E cr . (20)We again assume heating by cosmic-ray ionization and L H is the H cooling efficiency (erg s − cm ). With Eq. (17) for x H with η =
1, we obtain n cool , H = (cid:18)(cid:114) α B ζ D I UV k E cr ζ p L H (cid:19) / ≈ . ζ / − I / cm − . (21)Plugging back into Eq. (17) we get the H abundance at thecooling point, x H , cool = . × − (cid:18) ζ − I UV (cid:19) / . (22)In the numerical evaluations in Eqs. (21 and 22) we use L H = . × − erg cm s − at T , and E cr and φ , as inEq. (18). The corresponding pressure is P cool , H ≈ . n cool , H T ≈ . × ζ / − I / cm − K . (23)Unlike the Ly α -to-metal transition which is accompanied bya maximum in the P − n curve (Eq. 19), the Ly α -to-H tran-sition does not induce such a maximum but only a kink in the P − n curve (see Fig. 4).The critical metallicity below which H cooling becomesimportant is Z (cid:48) cool , H ≡ L H x H L z A O ≈ . × − (cid:18) ζ − I UV (cid:19) / . (24)As we discuss in §
5, this formula is in good agreementwith our numerical results. Importantly, the temperature de-creases more gradually for transitions from Ly α to H cool-ing. Because of the contributions from many lines in the ro-vibrational ladder (as illustrated in the bottom panel of Fig. 3),the H cooling function remains sensitive to temperature inthe entire range 10 < T < K. This is illustrated in theupper panel of Fig. 3 where the red curves show the muchsteeper H cooling function compared to the flat metals func-tion (blue curves). For this reason, the temperature drop at n cool , H is moderated and does not produce a local maximum in the P ( n ) curve. Most generally, the transition from Ly α toH cooling does not result in a multiphased behavior. While H cooling does not induce a multiphase, Fig. 4shows that multiphased behavior continues to be produced bymetal-line cooling, although over a much narrower pressurerange, down to a metallicity as low as 10 − . This multiphaseis induced by the transition from H to metal cooling (ratherthan Ly α to metals) and occurs for a lower cooling point tem-perature of T ≈
600 K. So long as the heating is by CR ion-ization the cooling point density at which the sharp tempera-ture drop occurs is still given by Eq. (18) but evaluated at thelowered T . Thus, in this regime n cool , Z continues to vary as1 / Z (cid:48) , but the density prefactor in Eq. (18) is reduced. For ex-ample, for Z (cid:48) = − , n cool , Z =
40 cm − , leading to a pressuremaximum of P max ≈ . n cool , Z T = . × cm − K for theWNM (see Fig. 4).Because H cooling reduces the temperature of the WNM,the multiphase strip is now narrower than at solar metallic-ity, and continues to shrink with decreasing metallicity. Forexample, while at solar metallicity the multiphase range iswide, with δ n ≡ ( n max − n min ) / n min = . δ P ≡ ( P max − P min ) / P min = .
9, for Z (cid:48) = − we have δ n = . δ P = .
88. For Z (cid:48) = − , the multiphase phenomena is alreadyextremely suppressed, with δ n = . δ P = . The onset of H heating at extremely low metallicity For still lower metallicities ( Z (cid:48) ≤ − ) sharp gas coolingand associated multiphased structures disappear entirely. Thedensity at which metal cooling can balance CR heating in-creases as 1 / Z (cid:48) , and at sufficiently low metallicities this den-sity becomes so high that H heating becomes operationaland dominant. In this limit metal cooling becomes ineffec-tive. This is because H heating includes a dependence onthe molecular formation rate (see Eq. A15) so that the heatingrate increases more rapidly than n . However, metal cool-ing varies no faster than n , and so cannot compete. In thislimit H is the dominant coolant, multiphase behavior doesnot occur, and the gas is driven to intermediate, thermally sta-ble, temperatures of ∼
600 K. This behavior is illustrated inFigs. 4 and 5.We obtain an analytic estimate for the critical metallicity, Z (cid:48) heat , H , below which metal line cooling is no longer effectiveas follows. We first estimate the density n heat , H at which H and CR heating are equal. For n > n heat , H we are in the H heating regime. The critical Z (cid:48) heat , H is then the metallicity atwhich the metal cooling density, n cool , Z , as given by Eq. (18),is at least as large as n heat , H . For n cool , Z (cid:38) n heat , H metal linecooling is ineffective because it cannot compete with H heat-ing.Equating the H and CR heating rates gives, n heat , H = (cid:32) ζ p E cr n crit E H , k (cid:112) ζ / α B (cid:33) / (25) = . × ζ / − T − . cm − , (26)where we used our H heating expression, Eq. (A15), assum-ing UV pumping heating dominates (but still in the n < n crit limit), and used Eq. (13) for the effective gas-phase formationrate coefficient R − . Next, equating the metal cooling density0 Bialy & Sternberg F IG . 5.— Cooling (left) and heating (right) rates per hydrogen nucleus (erg s − ), as functions of the gas hydrogen density, at various metallicities. hermal Phases – from Solar Metallicity to Primordial Gas 11 n cool , Z (Eq. 18) with n heat , H above gives, Z (cid:48) heat , H = . × − ζ / − . (27)In the numerical evaluation we adopted a metal cooling ef-ficiency L z = . × − erg cm s − at n = cm − , T = K. For Z (cid:48) (cid:46) Z heat , H the gas never undergoes a secondmetal-line cooling phase, and there is no multiphase structure.This formula is in agreement with the numerical results shownin Fig. 4.Eq. (27) is a strict requirement on the metallicity, express-ing the metallicity below which a multiphase cannot possiblyoccur. In practice, as discussed above, much before the metal-licity becomes that low, the multiphase region may alreadybecome very narrow, so that any thermal instability effectsmay become negligible. PARAMETER STUDYIn this section, in Figs. 6 and 7, we present results of ourthermal phase computations in the form of equilibrium tem-perature contour plots in the metallicity versus density ( Z (cid:48) − n )plane, as well as temperature and density contour plots in themetallicity versus pressure ( Z (cid:48) − P ) plane. We investigate theeffects of variations in ζ and I UV . We also discuss the timescales relevant for our assumptions of thermal and chemicalequilibrium. 5.1. Contour plots
The top panels of Fig. 6 show the equilibrium gas tem-perature in the metallicity versus density ( Z (cid:48) − n ) plane, for ζ − = I UV =
1. The colors indicate the temperature, rangingfrom cold ∼ K (blue), to intermediate ∼ K (yellow),to warm ∼ K (red). In the upper-right panel we overplotcontours of the cooling times. The cooling times range from ∼ −
10 Myrs in the WNM, depending on the metallicity,down to 0.1-0.4 Myrs in the CNM at low Z (cid:48) , and down to lessthen 0.01 Myrs at high Z (cid:48) . At higher metallicities the increasein the heating rate due to the onset of PE heating, is accompa-nied by appropriately higher cooling rates. For the CNM, thecooling rate naturally increases with increasing metallicity, as[C II ] and [O I ] dominates the cooling. In the WNM, whereLy α dominates, the cooling rate increases via a slight increasein the equilibrium temperature. Since at low Z (cid:48) a large portionof the parameter space is WNM, or the intermediate ≈ Kmedium (at very low Z (cid:48) ), with cooling times 1-10 Myrs, thegas may potentially be out of thermal equilibrium, dependingon the magnitude of the local dynamical time (e.g., the soundcrossing time, turbulent time, free-fall time, etc.). We furtherdiscuss the cooling and chemical time in the Appendix.In the upper-left panel the various curves delineate the dom-inating heating and cooling processes that operate in differentregions of parameter space. For heating, these are (a) photo-electric (PE) heating, above the dashed curve from the WNMinto the CNM, (b) H heating, to the right of the dashed curveat high density in the CNM at moderately reduced metallic-ity, and again at high densities and very low metallicities inthe intermediate temperature zone, and (c) cosmic-ray (CR)heating everywhere else. For cooling these are (a) Ly α to theleft of the solid curve in the low density WNM, (b) metals tothe right of the solid curve into the CNM, and (c) H coolingbelow the solid curve in the lower right portion of parameterspace.For any Z (cid:48) the gas temperature decreases monotonicallywith increasing n , with sharp drops in T , from the WNM to CNM (red to blue), occurring near the Ly α /metal or H /metalcooling point boundaries. The WNM-to-CNM transitionpoints depend on metallicity. The intermediate temperaturezone appears at very low Z (cid:48) and high density.For Z (cid:48) (cid:38) . Z (cid:48) (see Eq. 5), in this regime the Ly α to metalcooling point runs vertically (independent of Z (cid:48) ) at an ap-proximately constant density of ≈ . − . Below Z (cid:48) ≈ . n cool , Z ∝ / Z (cid:48) , consistent with Eq. (19) for acooling point temperature T = × K appropriate for thetransition from Ly α to metal cooling. Below a metallicityof 8 × − , H cooling (rather than metals) takes over fromLy α in the WNM. This is consistent with our analytic esti-mate (Eq. 24) for Z (cid:48) cool , H . The Ly α to H boundary line isvertical (independent of Z (cid:48) ) and occurs well inside the WNMat a density ≈ − consistent with our Eq. 21 for n cool , H .Below Z (cid:48) cool , H ≈ × − the WNM-to-CNM boundary re-mains diagonal with n cool , Z ∝ / Z (cid:48) as given by Eq. (19), butwith a lower T (i.e., T = Z (cid:48) = − . − . to metalcooling.For Z (cid:48) < Z (cid:48) heat , H ≈ − and n > n heat , H ≈ cm − , H heating dominates over CR ionization, and therefore metalcooling becomes ineffective. The critical values Z (cid:48) heat , H and n heat , H occur at the intersection of the CR-to-H heating andH -to-metal cooling boundary curves in Fig. 6, and are con-sistent with our analytic estimates given by Eqs. (25) and (27).In this portion of parameter space an intermediate equilibriumtemperature of ∼ K is set by the balance of H heating andcooling.The lower panels of Fig. 6 show the temperature (left)and density (right) in the pressure versus metallicity ( P − Z (cid:48) )plane (again for ζ − = I UV = zone at very low metallicity andhigh density. The multiphased region where the CNM andWNM overlap occurs between the solid black curves, from P min (to the left) to P max (to the right).The position of the multi-phased zone in Fig. 6 is deter-mined primarily by (a) the switch from PE to CR heating,(b) the dependence of the metal cooling point on Z (cid:48) , and (c)the onset of H cooling and heating at low metallicity. For Z (cid:48) down to ∼ . P min and P max .The width is maximal at Z (cid:48) ≈ . Z (cid:48) (cid:46) . P min and P max increase as the transition to metal cooling requires larger den-sities.As the metallicity is lowered below Z (cid:48) cool , H ≈ × − the multiphased zone become narrower due to H coolingin the WNM. At Z (cid:48) ≈ − there still exists a (small) tem-perature jump at the WNM to CNM boundary, but the pres-sure range of the multiphase zone is already negligible. At Z (cid:48) < Z (cid:48) heat , H ≈ − metal cooling becomes ineffective, thetemperature then decreases smoothly with increasing pres-sure, from the ∼ K Ly α cooled WNM, to the ∼ KH cooled intermediate temperature zone.2 Bialy & Sternberg F IG . 6.— Top: The temperature as a functions of metallicity and gas density. In the left panel, the solid contours delineate the regions where Ly α -, metal-,and H -cooling become dominant. The dashed contours delineate PE-, CR-, and H -heating dominated regions. The contours in the right panel are the coolingtimes in Myrs. Bottom: The temperature (left) and gas density (right) as a functions of metallicity and pressure. The WNM, CNM and multiphase regions areindicated. A 3D interactive surface plot is presented in the online version (see also https://plot.ly/ sb2580/47). Varying ζ , I UV and D − In Fig. 7 we explore the thermal structures for varying UVintensities and cosmic-ray ionization rates. The left columnshows the cases ζ − = . I UV =
1. Themiddle column is for I UV = . ζ − =
1. Theright column shows models in which I UV and ζ − scale to-gether by ± ζ − = I UV = Z (cid:48) − n plots, the PE/CR/H and Ly α /metals/H heating and cooling boundaries are shownas the solid and dashed curves, as in the upper left panel ofFig. 6. For the Z (cid:48) − P plots, the various phases includingthe multiphased zones are indicated as in the lower panels ofFig. 6.In the lower panels we also show, as dotted white curves, acomparison to the WMHT formula for P min as a function ofhermal Phases – from Solar Metallicity to Primordial Gas 13 F IG . 7.— As Fig. 6, but assuming variations in ζ and I UV . The contours, and shaded regions have the same meaning as in Fig. 6. The dotted curves in thelower panels are the WMHT formula for P min (Eq. 28). I UV , Z (cid:48) , Z (cid:48) d , and ζ − , P min = I UV Z (cid:48) d / ( Z (cid:48) f ) + . ( I UV Z (cid:48) d / ζ − ) . cm − K , (28)(their Eq. 33) over its validity range (their Eqs. 35-38). Weinclude the factor f ≡ ( − δ C Z (cid:48) d / Z (cid:48) ) / ( − δ C ) to correct forthe metallicity dependent depletion onto dust (i.e., our Eq. 4).For Z (cid:48) d we use our Eq. (5). The agreement is not perfect due tothe different treatment of cosmic-ray and X-ray heatings, andPAH chemistry. However, overall, the general trends with Z (cid:48) , ζ and I UV are recovered.The left-panels of Fig. 7 show that as ζ increases, theWNM, CNM and the multiphase regions, are all shifted tohigher densities and pressures. Within the CR-heating re-gion (below the almost horizontal dotted line) higher densi-ties are required for cooling to offset an increasing CR heatingrate. Within the PE-heating region, as ζ increases, the elec-tron abundance increases, the dust-grains are less positivelycharged, and the release of photoelectrons from the dust-grains into the gas is more efficient (i.e., see Fig. 1 in Wolfireet al. 1995). In the middle panels, where I UV varies (with ζ kept constant), the shift is observed only at high metallicities( Z (cid:48) (cid:38) .
1) where PE heating dominates, as expected.The critical metallicity at which H cooling becomes im-portant, Z (cid:48) cool , H , increases with ζ and decreases with I UV .When ζ and I UV scale together, Z (cid:48) cool , H remains unchanged.This is because the H abundance (a) increases with ζ , as thegas phase formation is more efficient, and (b) decreases with I UV , as photodissociation is more efficient. The observed de-pendence of Z (cid:48) cool , H on ζ and I UV is in good agreement withEq. (24). The location of the Ly α -to-H transition (the ver-tical white contour) is in good agreement with Eq. (21). Forexample, looking at the left panels of Fig. 7 we see that as ζ − increases from 0.1 to 10 (with I UV = α -to-H transition shifts from n ≈ − . On theother hand, when I UV varies from 0.1 to 10 (and ζ − = ≈ . ≈
20 cm − . Finally, when both ζ and I UV vary to-gether from 0.1 to 10, the transition point shifts by two ordersof magnitude, from ≈ . − .The critical metallicity, Z heat , H , at which H heating takesover, increases with ζ , and is independent of I UV , as predictedby Eq. (27). As discussed in § heating equal CR heating, which inturn must equal metal cooling, all of which are independentof I UV . An increase in the CR intensity increases both theH formation rate and the CR heating rate, but the latter isaffected more strongly, and thus Z heat , H increases with ζ .As discussed in § − photodetachment rate atlow metallicity is uncertain. We consider variations be-tween D − = . × − s − for the extrapolated Draine field(our fiducial models), up to D − = . × − s − for aDraine+Mathis field. To investigate the effect of a varying D − , we have also computed a set of models with the higherDraine+Mathis D − rate. The resulting temperature map isidentical to that in Fig. 6, with the only difference that the H cooling line is shifted to a higher density, n cool , H =
15 cm − ,i.e., an increase by a factor of 3.5. Since the n cool , Z curve is While WMHT use elemental abundances similar to ours for Z (cid:48) =
1, theydo not have a metallicity dependent depletion in their model, and f correctstheir formula for this effect. unaffected by variations in D − (as it is set by metal cooling),the critical metallicity for H cooling is reduced, and Z cool , H is reduced to 10 − . The shift to a higher cooling density is ex-pected as the higher detachment rate reduces the H formationrate, requiring a higher gas density for the onset of H cool-ing. Unlike n cool , H and Z cool , H , the heating points n heat , H and Z heat , H are unaffected by the value of D − , as they occurat very high densities where H − photodetachment is subdom-inant.Quantitatively, the factor 3.5 shift in n cool , H may be un-derstood as follows. For the higher D − rate, photodetach-ment dominates H − removal for densities n < I UV T . cm − . For such densities, η ≈ . T − . I − ( n / − ) (seeEq. 16), and the H abundance is x H = x e k k n D D − I = . × − T . ζ / − I − (cid:16) n
10 cm − (cid:17) / . (29)At T = n =
10 cm − (i.e., near the Ly α -to-H transi-tion) this x H is a factor of ≈
20 smaller than the correspond-ing value in Eq. (17), and has a steeper dependence on den-sity. Plugging this back into Eq. (29), we see that n will nowenter with a power 5/2, thus the factor 20 lower x H , needsto be compensated by an increase in density by a factor of ≈ / ≈ .
3, in agreement with our numerical results.Importantly, while ζ , I UV and D − affect the values of thecritical densities and metallicities, the qualitative behaviour,the shift of the WNM-to-CNM transition to higher pressureswith decreasing metallicity, the quenching of the multiphasezone due to H cooling at low Z (cid:48) , and the final disappearanceof the multiphase phenomenon at very low Z (cid:48) , is robust. SUMMARY AND DISCUSSIONIn this paper we have studied the thermal phase propertiesof atomic gas from high to vanishing metallicity, with spe-cial emphasis on understanding how multiphase behavior isaltered when H heating and cooling processes become domi-nant. As discussed in §
1, the critical role of H cooling at lowmetallicity has long been recognized in many studies. How-ever, multiphased CNM/WNM behavior into the low metal-licity and primordial regimes has received less attention.For solar metallicity, we reproduce the known result that aCNM phase cooled by [C II ] and [O I ] line emission, ( T CNM ∼
100 K) and a WNM phase cooled by Ly α emission ( T ∼ I UV , and/or increasing ionization rate, ζ , asphotoelectric (PE) heating becomes more efficient.As the metallicity decreases, the metal cooling rate de-creases. However, as long as Z (cid:48) (cid:38) .
1, PE heating continuesto dominate so that the heating rate also decreases. For suffi-ciently low metallicities, Z (cid:48) d is expected to scale superlinearlywith Z (cid:48) (in our model this occurs for Z (cid:48) < . Z (cid:48) decreases, and theCNM is then colder and is less dense compared to solar metal-licity models. The pressure range that allows a multiphase isthen larger compared to solar metallicity models.For metallicities Z (cid:48) (cid:46) .
1, the ionization of neutral gas bycosmic-rays (CR) becomes the dominant heating mechanismof the gas (the metallicity value depends on the CR ioniza-tion rate). In this limit as the metallicity decreases, the metalhermal Phases – from Solar Metallicity to Primordial Gas 15cooling rate decreases, while heating remains constant, shift-ing the WNM-to-CNM transition to higher densities and pres-sures. In this limit, and when H heating and cooling pro-cesses are ignored, we find that the Ly α -to-metal transitiondensity (and associated WNM-to-CNM transition) scales as ζ / Z (cid:48) .However, while previous analytic models of the H I phasestructure often ignored H cooling, we find that H cool-ing dominates and modifies the phase structure at low metal-licity (as we show H cooling and heating also play a roleat solar metallicity gas but only at very high densities, ∼ cm − , deep in the CNM). The critical metallicity be-low which H cooling becomes important is Z (cid:48) cool , H = . × − ( ζ − / I UV ) / . Below this metallicity, H cooling re-duces the temperature of the WNM, to temperatures as lowas 600 K (the exact value depending on metallicity). Thislowers P max and results in a narrower multiphase zone, bothin pressure and density.Finally, at extremely low metallicities Z (cid:48) < Z (cid:48) heat , H = . × − ζ / − , the density at the WNM-to-CNM transition is al-ready so high that H heating starts contributing. In this limitthe thermal structure is determined only by H I (Ly α ) and H cooling. There is a single-phase solution at any pressure, asthe temperature decreases smoothly with gas pressure, fromthe ≈ − T ≈
600 K, as determinedby the balance of H cooling and heating.Gravitational collapse is favored in the cold-dense CNM,where the free-fall time is relatively short. Our result thatthe transition from warm to cold gas occurs at ever increas-ing pressures as the metallicity is reduced (over the range10 − (cid:46) Z (cid:48) (cid:46) .
1, for a constant I UV and ζ ), suggests thatgalaxies with metal-poor ISM (i.e., galaxies in early evolu-tionary stages or low mass dwarf galaxies), should have main-tained very high interstellar pressures in order to allow the for-mation of a CNM phase. For example, for our fiducial modelof I UV = ζ − =
1, and for a metallicity Z (cid:48) = − a minimumpressure, P / k ≈ cm − K is required. This requires abnor-mally large ISM column densities (for hydrostatic pressureequilibrium), deep gravitational potential wells of the dark-matter halos, or a highly pressurized hot ambient gas that mayprovide the needed pressure support.Alternatively, a WNM-to-CNM transition and a multiphaseISM may still occur at normal ISM pressures at low metallic-ity if the heating rate is reduced. This requires a reduced ζ and I UV values, suggesting lower star-formation (supernovae)rates in low metallicity galaxies. For example, for standardGalactic conditions ( Z (cid:48) = I UV = ζ − = P ≈ − K (Fig. 6). For Z (cid:48) = − the multi-phase may still occur at a similar pressure if I UV and ζ arereduced by a factor of 10 (Fig. 7, rightmost column, 3 rd row).When connected to a theory of a self-regulated star formingISM (Ostriker et al. 2010), this would imply that the star-formation rate per unit gas column (i.e., the normalization ofthe Kennicutt Schmidt relation) should be lower at very lowmetallicities compared to solar (e.g., a factor of 10 lower for Z (cid:48) = − compared to Z (cid:48) = ≈ K. In such a case, a differ-ent mode of star-formation would be expected, with only verymassive clouds susceptible to gravitational collapse.Previous studies have suggested that the atomic ISM of star-forming galaxies is naturally driven towards the multiphasepressure range, which provides a regulation mechanism forstar-formation (Parravano 1988, 1989; Ostriker et al. 2010).The idea is as follows. If the star-formation rate of a galaxyincreases this leads to an increase in the UV intensity whichin turn increase the PE heating rate. The P − n curve is thenpushed up and to the right (see Fig. 1, red-dashed curve), thegas is driven to the warm (WNM) phase, resulting in a de-crease in the efficiency of star-formation, and vice versa (seeFig. 1 in Ostriker et al. 2010). In the limit where CR domi-nates over PE heating, the regulation is also expected to hold,since CRs are accelerated in supernovae remnants, and thusthe CR ionization rate is also likely proportional to the star-formation rate.If multiphase structure is indeed fundamental for star-formation regulation in galaxies, then our model providesa prediction for the thermal pressure of the atomic ISM ingalaxies, as a function of metallicity, CR ionization rate andUV intensity (i.e., see Figs. 6, 7). This pressure may be thenrelated to the star-formation rate and galaxy gas and dark mat-ter halo, through the requirement of hydrostatic equilibrium.At extremely low metallicities, the multiphase structuredisappears, and the multiphase-regulation mechanism cannothold. This may imply a different (potentially a bursty) modeof star-formation for galaxies in this early evolutionary stage.In this limit the dense gas remains warm, at ≈
600 K, and iscooled by H cooling. This temperature is a factor of ∼ ≈ / ≈
30 higher.Stellar structures forming in this low metallicity regime areexpected to be more massive compared to their high metallic-ity counterparts.Our finding that the cooling of the atomic ISM at low metal-licity (below ≈ .
01) is often dominated by H line emissionfrom the moderately warm phase, may open up a new obser-vational window for the study of the low metallicity ISM, inlocal dwarfs and in galaxies of early evolutionary stages athigh redshift. This is particularly timely with the upcominglaunch of the James Webb Space Telescope (JWST) whichwill be able to observe IR H lines at extreme sensitivity.We thank John Black, Emeric Bron, Franck Le Pettit, Eve-lyne Roueff, Chris McKee, and Eve Ostriker, for fruitful dis-cussion and helpful suggestions. We thank the referee forconstructive comments that improved our manuscript signifi-cantly. SB thanks the Center for Computational Astrophysicsat the Flatiron Institute for hospitality and funding wheresome of this research was carried out with AS. This workwas also supported by the German Science Foundation viaDFG/DIP grant STE 1869/2-1 GE 625/17-1 at Tel Aviv Uni-versity.6 Bialy & SternbergAPPENDIXA. HEATING AND COOLING PROCESSES Line-Emission Cooling
In general, gas cooling through line-emission occurs whenever a particle (ion/atom/molecule) is collisionaly excited to a higherenergy state, and then relaxes back to a lower energy state by (spontaneously) emitting a photon. In this process, part of the kineticenergy of the collider particle is transferred to the emitted photon, leading to gas cooling.While we compute the cooling rates numerically (as outlined below), it is instructive to consider the simple case of a two-levelsystem (e.g., the [C II ] 158 µ m transition), for which an analytic expression for the cooling rate may be obtained. Let ∆ E be thetransition energy, and x i and x j the abundances of the cooling species and the collisional partners, respectively. The cooling rateper unit volume (erg cm − s − ) is L i = x i n L ( T ) (cid:122) (cid:125)(cid:124) (cid:123) ∑ j x j q i j , ↑ ( T ) ∆ E + n / n crit , (A1)where n crit ≡ A / q i , ↓ is the critical density at which the collisional deexcitation and radiative decay rates are equal, A is the Einsteincoefficient for spontaneous decay, and q i j , ↑ and q i j , ↓ the collisional rate coefficients for the upward and downward transitionsrespectively. The latter are related through ( q i j , ↑ / q i , ↓ ) = ( g u / g d ) exp [ − ∆ E / ( k B T )] , where g u , g d , are the quantum degeneracies.The combination ∑ j x j q i j , ↑ ( T ) ∆ E ≡ L i ( T ) is the cooling coefficient . For densities n (cid:28) n crit , the cooling rate L ∝ n , as thedensity of the coolant and the rate of collisions (which populate the upper level) are both ∝ n . At n (cid:29) n crit the energy levels arethermalized and L ∝ n . In this limit, L i (cid:39) x i n n crit ( T ) L ( T ) ≡ x i n Λ LTE ( T ) , (A2)where Λ LTE is cooling rate per atom/molecule (erg s − ) at local thermal equilibrium (LTE). Eq. (A1) is exact for two level systems(assuming optically thin gas), but it may also provide a useful approximation for multilevel systems (such as [O I ] and H ), withthe appropriate choice of an effective ∆ E and n crit .We calculate the line emission cooling for Ly α , [C II ] 158 µ m, [O I ] 63 and 146 µ m, and the ro-vibrational transitions of H and HD. Ly α cooling: We compute the cooling via radiative decay from the 2p and 2s states of neutral hydrogen collisionaly excitedby electrons. For the collisional de-excitation rate coefficient, we adopt q ↓ = . × − Ω / ( g u (cid:112) T / K ) cm s − , g u = ( l u + ) where l u = Ω is the dimensionless collision strength(Gould & Thakur 1970; Spitzer 1978). We adopt a constant Ω = .
26 for the 2s level, and 0 . q ↓ = . × − / √ T cm s − , in excellent agreement with (Gould & Thakur 1970;Spitzer 1978) as well as with the more recent calculations of Callaway et al. (1987), over the temperature range of interest,6000 < T < K. Metal fine-structure cooling:
We calculate the cooling arising from the from the fine-structure transition of [C II ] P o3 / − P o1 / at λ = . µ m, and the fine structure transitions of [O I ] P i − P j where [ i , j ] = ([ , ] , [ , ] , [ , ]) , with the most importantbeing the [ i , j ] = [ , ] at λ = . µ m. We include collisional excitations with electron, hydrogen atoms and para-H and ortho-H . We use the analytic fits presented in Draine (2011, see his Table 6), which are based on Pequignot (1996); Tayal (2008);Barinovs et al. (2005); Flower & Launay (1977); Abrahamsson et al. (2007); Jaquet et al. (1992). We also include cooling from[C I ] (although we find it to be a subdominant coolant), using the analytic fits from Draine (2011, see his Table 6), based on datafrom Roueff & Le Bourlot (1990); Abrahamsson et al. (2007); Schroder et al. (1991); Staemmler & Flower (1991). H ro-vibrational cooling: For the H cooling function we use the polynomial fits provided by Glover & Abel (2008) withupdates from Glover (2015, see their Appendix A1.2), which account for excitations of the H ro-vibrational ladder throughcollisions with H, He, para-H , ortho-H , H + , and e, assuming a constant ortho-to-para ratio o/p=3. We find that our results arenot sensitive to the adopted o/p value. The cooling functions are based on collisional data from Ehrhardt et al. (1968); Cromptonet al. (1969); Linder & Schmidt (1971); Draine et al. (1983); Gerlich (1990); Krsti´c & S. (2002); Wrathmall et al. (2007); Flower& Roueff (1998); Flower et al. (1998); Balakrishnan et al. (2002); Flower & Roueff (1999); Honvault et al. (2011). We havealso calculated the H cooling efficiency assuming the updated H-H collision rates of (Lique 2015). We find that adopting thisalternative collisional dataset results in at most a negligible difference in the total H cooling function (which includes all thecollisional partners). HD ro-vibrational cooling:
For HD cooling we use the cooling function of Lipovka et al. (2005) , based on collisional datafrom Roueff & Flower (1999); Roueff & Zeippen (1999). Heating and Cooling by Dust
Dust grains and PAHs, contribute both to gas heating, through the injection of energetic photoelectrons, and to cooling, thoughdust/PAH-assisted recombination. For the photoelectric (PE) heating we follow Bakes & Tielens (1994) and WMHT and adopt G pe = γ pe ε Z (cid:48) d I UV n , (A3) We spotted typos in the formulae in Lipovka et al. (2005) (confirmedby the authors): (a) T and n in their Eq. (4) and (5), should be replaced with log T and log n . (b) In their Eq. (5), W HD should be replaced with W HD / n . hermal Phases – from Solar Metallicity to Primordial Gas 17where the PE efficiency is ε = . × − + . × − y . + . × − T . + . × − y (A4)and y ≡ I UV ( T / K ) / ( n e / cm − ) φ PAH . (A5)In Eq. A3 and A4, we assume that the PAH abundance scales with the dust-to-gas ratio, Z (cid:48) d , and adopt γ = . × − erg s − and φ PAH = . ). For recombination cooling we adopt L rec = γ rec ( T / K ) . ( . y ) β φ PAH Z (cid:48) d x e n , (A6)with γ rec = . × − erg cm s − .In our model, we do not solve for the PAH chemistry. At high metallicity, the PAHs may affect the electron fraction as theyintroduce an additional recombination channel for the H + (see WMHT). At low metallicity, this effect becomes unimportant asthe abundance of the PAHs falls rapidly with the vanishing DGR. Even at high metallicity, where the PAHs are importnat forrecombination, we find that the net effect on the phase diagrams is mild, as shown in Fig. 7 where we compare to WMHT (seealso Fig. 8 in WMHT for the effect of PAH abundance on the phase diagrams).We include cooling due to collisions of gas particles with dust-grains. We adopt L g − d = γ g − d ( T / K ) / ( T − T d ) (cid:16) − . − / T (cid:17) Z (cid:48) d n , (A7)where T d is the dust temperature, and γ g − d = . × − erg cm s − , (Hollenbach & McKee 1989; Glover & Jappsen 2007).The dust is typically much colder than the gas, and is approximated by T d = . I / K (Draine 2011).
Cosmic-Ray Heating
Cosmic-ray or X-ray ionization liberate energetic electrons into the gas. These energetic electrons lose energy through coulombinteractions, and ionization, dissociation and excitations. The fraction of energy that is lost into coulomb interactions heat thegas. We adopt G cr = ζ p E cr n , (A8)for the CR heating rate, where ζ p = ζ / ( + φ s ) is the primary CR ionization rate, φ s = (cid:16) − x e . (cid:17) . + ( x e / . ) , (A9)is the number of secondary ionizations per each primary electron, and E cr = . (cid:32) + . (cid:20) x e x e + . (cid:21) / (cid:33) eV , (A10)is the characteristic deposition energy per ionization event (Dalgarno & McCray 1972; Draine 2011). The free electrons in thegas can also interact directly with the CR protons, giving an additional heating rate G cr = ζ p E cr nx e . (A11)For CR proton energies of ∼
50 MeV, E cr =
287 eV (Goldsmith et al. 1969). This direct heating channel dominates when theelectron fraction is larger than x e ∼ . H Heating
When LW photons are absorbed by H molecules, they lead to population of the excited electronic (LW) states. ∼
10 % of theradiative decays lead to H photodissociation, whereas in the rest 90 % of the cases, the molecule decays through ro-vibrationaltransitions back to the ground state. If the gas density is sufficiently high, such that collisions rather than spontaneous emissiondeexcite the H , the energy of the excited level is transferred to the colliding particle (H, e, H , etc.) leading to gas heating.Following Burton et al. (1990) and R¨ollig et al. (2006), the pumping heating rate may be approximated with G H , pump = D I UV E pump x H n + n crit / n = Rn x H E pump + n crit / n , (A12)where 9 D is the H pumping rate (9 excitations occur per dissociation), and E pump and n crit are the effective energy and criticaldensity of the pseudo two-level transition (see discussion above). The second equality holds for chemical equilibrium, for which D I UV x H = Rnx H where R is the total H formation rate coefficient. In this limit, the UV pumping heating rate is proportional The numerical values in Eq. A3 and A4 differ from those presented inWMHT because we used I UV as the basic UV intensity unit, where WMHT use the Habing (1968) G units. The conversion is, 1 I UV = . G . formation rate, Rn . At low densities, below n crit , another factor of n / n crit enters, reflecting the fact that when n < n crit spontaneous emission dominates deexcitation, not collisions. Following R¨ollig et al. (2006) we adopt E pump = .
12 eV, and n crit = . × / √ T cm − for the effective energy released per UV pumping event and for the critical density.We also include heating through H photodissociation, and adopt G H , pd = D I UV E pd x H n = R n x H E pd (A13)where E pd = . n / n crit < E pd / ( E pump ) , which occurs at n (cid:46) . × / √ T cm − , for our adopted parameters.When a new H molecule is formed, part of the binding energy (4.5 eV) is converted into heat, either by providing translationalenergy of the H molecules, or by exciting the H ro-vibrational levels which when followed by collisional de-excitation lead togas heating. The latter is effective when the gas volume density is large, of order of or higher than the critical density of H . Theheating resulting from H may be written as G H , form = x H n (cid:40) ∑ i R i (cid:18) E i , form , + E i , form , + n crit / n (cid:19)(cid:41) (A14)where the summation is over the formation channels of H , (i) formation on dust grains, and (ii) the gas-phase H − formationsequence, and R i are given by Equations (9) and (13), for the dust and gas formation routes, respectively. For each of the formationprocesses, the first term in brackets, E form , , is the energy that goes into translational energy of the products (for dust-formation,this is the H . For the gas phase H − formation route, most of the energy goes to the electron), and E form , is the energy thatgoes into internal molecule excitation, which can then heat the gas through collisional deexcitation. We adopt E form , = . E form , = .
48 eV for the dust catalysis (Hollenbach & McKee 1979, see their Appendix VI(c)). For the H − formation wefollow ˇC´ıˇzek et al. (1998) who finds that on average ≈ . E form , = . E form , = .
13 eV.We may obtain a simple form for the total heating rate by all of the H heating processes discussed above. Assuming thatalways a single formation route dominates (either the dust catalysis or the gas phase H − route) and that the H system hasreached chemical equilibrium, we get G H = Rn x H (cid:18) E H , + E H , + n crit / n (cid:19) . (A15)Here we defined E H , = E form , + E pd and E H , = E form , + E pump and assumed that the effective critical densities are the samefor H formation and UV pumping. In practice, since these critical densities are effective two-level system approximations, theirvalues may differ for different processes. With our assumed parameters, E H , = ( . , . ) eV and E H , = ( . , . ) eV forH formation on dust and via H − , respectively. The second term in Eq. (A15) is dominated by H pumping heating, and becomesgreater than the first term when n (cid:38) ( . , . ) × / √ T cm − , for H formation on dust and via H − , respectively.B. REDUCED CHEMICAL NETWORKTo improve computational speed, we consider a reduced chemical network consisting of the species: H, H − , H + , H , He,C + , O, and e. The chemical reactions considered in the reduced network are listed in Table 2, along with the rate coefficients(based on the UMIST 2012 database, McElroy et al. 2013). For a thorough discussion of these chemical reactions see § is formed via dust catalysis and by the H − gas-phase route, and is destroyed by photodissociation. The H − is TABLE 2R
EDUCED N ETWORK
Reaction Rate coefficient (subscript refers to the reaction no. in the text)H + H : dust → H R d = × − T / Z (cid:48) d cm s − H + e → H − + ν k = . × − T . e − . / T cm s − H − + H → H + e k = . × − T − . e − . / T cm s − H − + H + → H + H k = . × − T − . cm s − e + H + → H + ν α B = . × − T − . cm s − H + H → H + H + H Formula from Martin et al. (1996)Reaction Photo/CR-rateH + CR → H + + e ζ = − ζ − s − H + ν → H + H D I UV = . × − I UV s − H − + ν → H + e D − I UV = . × − I UV s − C + , O, He assume all gas-phase carbon/oxygen/helium are in C + , O, He formed via radiative association, and is destroyed by associative detachment, photodetachment, and mutual neutralization. Forthe electrons we include production by CR-ionization of hydrogen, and photoionization of carbon assuming all the carbon is inhermal Phases – from Solar Metallicity to Primordial Gas 19the form of C + . The electrons are removed by radiative recombinations the protons. We assume the oxygen and helium are inneutral atomic form. The proton abundance is set by the requirement of charge neutrality.This reduced network enables us to obtain analytic expressions for the species abundances and thus to improve significantlythe computation speed. The reduced network was used to produce Figs. 6 and 7. We have confirmed the accuracy of the reducednetwork by comparing the equilibrium T ( n ) , and x H ( n ) curves for the reduced and full networks at the various metallicitiesconsidered in Fig. 4. C. TIME-SCALESThermal and chemical equilibrium are good assumptions whenever the dynamical time of interest, t dyn , e.g., the turbulentcrossing time, the free-fall time, the cloud lifetime, etc., is long compared to the cooling time, t cool , and the chemical time, t chem . The cooling/heating time
The cooling time from an initial temperature T i to an equilibrium temperature T is t cool = (cid:90) TT i ( / ) nk B TG ( n , T ) − L ( n , T ) d T , (C1)where we assumed isochoric cooling of predominantly atomic gas (an adiabatic index γ = / t cool ≈ ( / ) nk B TL ( n , T ) ≈ ( / ) nk B TG ( n , T ) . (C2)This is a good assumption when the cooling function increases superlinearly with T , as occurs when H or Ly α emission dominatethe cooling, for T (cid:38) × K, and if Z (cid:48) < Z (cid:48) cool , H down to ∼ §
4, and Fig. 3 below), or at lowtemperatures, T (cid:46)
100 K where metal cooling dominates but is a steep function of the temperature. On the other hand, Whenmetals dominate the cooling at T (cid:38)
100 K, the cooling function depends weakly on T (with a power-index < T replaced by T i .In the second equality in Eq. (C2) we used the fact that the cooling time is dominated by the time spent near the equilibriumtemperature, where G ≈ L . This latter form of t cool is particularly useful when G scales approximately linearly with the densitywhile being weakly dependent on temperature. For example, for CR heating, G = ζ p E cr n and we get t cool = ( / ) k B T ζ p E cr ≈ . T ζ − − Myr . (C3)In the numerical evaluation we used E cr = . φ s = .
57 (the same values as in Eq. 18). We see that for typical ζ − = t cool ∼
20 Myrs, for T = t cool typically 0.1-0.4 Myrs Since the heating rate is mostly independent of temperature (e.g., for CR heating),Eqs. (C2-C3) also express the heating time from T i to T for any arbitrary T i and T . Chemical Times
Generally, the time for the H (and H) abundance to reach steady-state is t H = DI UV + Rn . (C4)where D is the local photodissociation rate. In highly shielded cloud interiors, D (cid:28) Rn , and the H time equal the H formationtime, 1 / ( Rn ) , which may become long. However, in our models we consider the state of the ambient atomic medium exposed tothe free-space UV field, with D = D . In this regime the H time is the H (free-space) photodissociation time, 1 / D = I − yr,and is always very short.Other relevant chemical timescales are the e, C + , and O, timescales. At low metallicity, the steady-state electron abundanceis set by the ionization-recombination equilibrium, as given by Eq. (12). Under the assumption of x e <
1, the time-scale for thisprocess to achieve steady-state is t e = α B n = . T . n − Myr , (C5)where n ≡ n / ( cm − ) . The electron time is already short ( (cid:28) t e ≈ T . / n years. For the C + time, consider a gas that recombines from a higher ionization state, C + + e → C + . Therecombination time is t C + = / ( nx e α C + ) ≈ ( n ζ − ) − / Myr, where we used α C + = . × − cm s − (Nussbaumer & Storey1969, see Table 14.7 in Draine 2011), and used the equilibrium electron abundance (Eq. 12; we can assume equilibrium for x e since t e (cid:28) t C + ). The timescale for the chemical equilibrium of O may be longer than that of C + , as the recombination rate for O + + . Nevertheless, as long as C + has reached equilibrium, the effect on the temperature will not be significant,since both coolants operate at the same tempearture regime. REFERENCESAbel, T., Bryan, G. L., & Norman, M. L. 2000, ApJ, 540, 39—. 2002, Modes of Star Formation and the Origin of Field Populations,Astronomical Society of the Pacific Conference Series, 285Abrahamsson, E., Krems, R. V., & Dalgarno, A. 2007, ApJ, 654, 1171Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARAA, 47, 481Audit, E., & Hennebelle, P. 2005, AA, 433, 1Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822Balakrishnan, N., Vieira, M., Babb, J. F., et al. 2002, AA, 524, 1122Barinovs, G., van Hemert, M. C., Krems, R., & Dalgarno, A. 2005, ApJ,620, 537Barkana, R., & Loeb, A. 2001, Physics Reports, 349, 125Bialy, S., & Sternberg, A. 2015, MNRAS, 450, 4424—. 2016, ApJ, 822, 83Bialy, S., Sternberg, A., & Loeb, A. 2015, ApJ, 804, L29Black, J. H., & Dalgarno, A. 1977, ApJS, 34, 405Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARAA, 51, 207Braun, R., & Walterbos, R. A. M. 1992, ApJ, 386, 120Bron, E., Le Bourlot, J., & Le Petit, F. 2014, AA, 569, A100Burton, M. G., Hollenbach, D. J., & Tielens, A. G. G. M. 1990, ApJ, 365,620Callaway, J., Unnikrishnan, K., & Oza, D. H. 1987, Physical Review A, 36,2576Cardelli, J. A., Meyer, D. M., Jura, M., & Savage, B. D. 1996, ApJ, 467, 334Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29ˇC´ıˇzek, M., Hor´aˇcek, J., & Domcke, W. 1998, Journal of Physics B: Atomic,Molecular and Optical Physics, 31, 2571Cojazzi, P., Bressan, A., Lucchin, F., Pantano, O., & Chavez, M. 2000,MNRAS, 315, L51Colgan, S. W. J., Salpeter, E. E., & Terzian, Y. 1988, AA, 328, 275Corbelli, E., & Salpeter, E. E. 1988, ApJ, 326, 551Crompton, R., Gibson, D., & McIntosh, A. 1969, Australian Journal ofPhysics, 22, 715Dalgarno, A., & McCray, R. A. 1972, ARAA, 10, 375Dalgarno, A., & Rudge, M. R. H. 1964, ApJ, 140, 800Dalgarno, A., Yan, M. I. N., & Liu, W. 1999, ApJ, 125, 237Dickey, J. M., & Brinks, E. 1993, ApJ, 405, 153Dickey, J. M., McClure-Griffiths, N. M., Gaensler, B. M., & Green, A. J.2003, ApJ, 585, 801Dickey, J. M., Mebold, U., Stanimirovic, S., & StaveleySmith, L. 2000, ApJ,536, 756Dickey, J. M., Terzian, Y., & Salpeter, E. E. 1978, ApJ Supplement Series,36, 77Draine, B. T. 1978, ApJS, 36, 595—. 2011, Physics of the Interstellar and Intergalactic MediumDraine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485Ehrhardt, H., Langhans, L., Linder, F., & Taylor, H. S. 1968, PhysicalReview, 173, 222Elmegreen, B. G., & Parravano, A. 1994, ApJ, 435, L121Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, AA, 155, L149Flower, D. R., & Launay, J. M. 1977, Journal of Physics B: Atomic andMolecular Physics, 10, 3673Flower, D. R., & Roueff, E. 1998, Journal of Physics B: Atomic, Molecularand Optical Physics, 31, 2935—. 1999, Journal of Physics B: Atomic, Molecular and Optical Physics, 32,3399Flower, D. R., Roueff, E., & Zeippen, C. J. 1998, Journal of Physics B:Atomic, Molecular and Optical Physics, 31, 1105Galli, D., & Palla, F. 1998, AA, 335, 403Garwood, R. W., & Dickey, J. M. 1989, ApJ, 338, 841Gazol, A., & Kim, J. 2013, ApJ, 765, 49Gazol, A., & Villagran, M. A. 2016, MNRAS, 462, 2033Gerlich, D. 1990, The Journal of Chemical Physics, 92, 2377Glover, S. C. 2015, MNRAS, 451, 2082Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627Glover, S. C. O., & Clark, P. C. 2014, MNRAS, 437, 9Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1Glover, S. C. O., & Mac Low, M. 2007, ApJS, 169, 239Goldsmith, D. W., Habing, H. J., & Field, G. B. 1969, AA, 158, 173Gould, R. J., & Thakur, R. K. 1970, Annals of Physics, 61, 351Habing, H. 1968, BAN, 19, 421Haiman, Z., Thoul, A. A., & Loeb, A. 1996, ApJ, 464, 523Heiles, C., & Troland, T. H. 2003, ApJ, 586, 1067Hennebelle, P., & P´erault, M. 1999, AA, 351, 309Herrera-Camus, R., Bolatto, A., Wolfire, M., et al. 2017, AA, 835, 201Hill, A. S., Mac Low, M.-M., Gatto, A., & Ib´a˜nez-Mej´ıa, J. C. 2018, ApJ,862, 55Hirasawa, T. 1969, Progress of Theoretical Physics, 42, 523Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555—. 1989, ApJ, 342, 306 Honvault, P., Jorfi, M., Gonz´alez-Lezana, T., Faure, A., & Pagani, L. 2011,Physical Review Letters, 107, 023201Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., & Clark, P. C. 2016,MNRAS, 458, 3528Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91Inoue, T., & Omukai, K. 2015, ApJ, 805, 73Jaquet, R., Staemmler, V., Smith, M. D., & Flower, D. R. 1992, Journal ofPhysics B: Atomic, Molecular and Optical Physics, 25, 285Jeans, J. H. 1902, Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences, 199, 1Kim, C.-G., Kim, W.-T., & Ostriker, E. C. 2011, ApJ, 743, 25Kolpak, M. A., Jackson, J. M., Bania, T. M., & Dickey, J. M. 2002, ApJ,578, 868Kritsuk, A. G., & Norman, M. L. 2002, ApJ, 580, L51Krsti´c, P. S., & S., P. 2002, Physical Review A, 66, 042717Kulkarni, S. R., & Heiles, C. 2011, in Interstellar processes; Proceedings ofthe Symposium, Grand Teton National Park, WY, July 1-7, 1986(A88-14501 03-90). Dordrecht, D. Reidel Publishing Co., 1987, p.87-122. NSF-supported research., 87–122Le Bourlot, J., Le Petit, F., Pinto, C., Roueff, E., & Roy, F. 2012, AA, 541,A76Lepp, S., & Shull, J. M. 1984, ApJ, 280, 465Linder, F., & Schmidt, H. 1971, Zeitschrift f¨ur Naturforschung A, 26, 1603Lipovka, A., N´u˜nez-L´opez, R., & Avila-Reese, V. 2005, MNRAS, 361, 850Lique, F. 2015, MNRAS, 453, 810Liszt, H. 2002, AA, 389, 393Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461, 265Marx-Zimmer, M., Herbstmeier, U., Dickey, J. M., et al. 2000, AA, 354, 787Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, AA, 128, 212McDowell, M. R. C. 1961, The Observatory, 81, 240McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, AA, 550, A36McKee, C. F., & Krumholz, M. R. 2010, ApJ, 709, 308Mebold, U., D¨usterberg, C., Dickey, J. M., Staveley-Smith, L., & Kalberla,P. 1997, ApJ, 490, L65Meyer, D. M., Jura, M., & Cardelli, J. A. 1997, ApJ, 493, 222Miyake, S., Stancil, P. C., Sadeghpour, H. R., et al. 2010, AA, 709, L168Murray, C. E., Stanimirovi´c, S., Goss, W. M., et al. 2018, ApJS, 238, 14Norman, C. A., & Spaans, M. 1997, ApJ, 480, 145Nussbaumer, H., & Storey, P. J. 1969, AA, 126, 75Omukai, K. 2000, ApJ, 534, 809Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41Parravano, A. 1988, AA, 205, 71—. 1989, ApJ, 347, 812Patra, N. N., Kanekar, N., Chengalur, J. N., & Roy, N. 2018, MNRASLetters, 479, L7Peebles, P. J. E., & Dicke, R. H. 1968, ApJ, 154, 891Pequignot, D. 1996, AA, 313, 1026R´emy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2013, AA, 563A, 31Richings, A. J., & Schaye, J. 2016, MNRAS, 458, 270R¨ollig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & Sternberg, A. 2006,AA, 451, 917Roueff, E., & Flower, D. R. 1999, MNRAS, 305, 353Roueff, E., & Le Bourlot, J. 1990, AA, 236, 515Roueff, E., & Zeippen, C. J. 1999, AA, 343, 1005Roy, N., Kanekar, N., Braun, R., & Chengalur, J. N. 2013, MNRAS, 436,2352Saslaw, W. C., & Zipoy, D. 1967, Nature, 216, 976Saury, E., Miville-Deschˆenes, M.-A., Hennebelle, P., Audit, E., & Schmidt,W. 2014, AA, 567, A16Schaye, J. 2004, ApJ, 609, 667Schneider, R., Omukai, K., Limongi, M., et al. 2012, MNRAS: Letters, 423,L60Schroder, K., Staemmler, V., Smith, M. D., Flower, D. R., & Jaquet, R. 1991,Journal of Physics B: Atomic, Molecular and Optical Physics, 24, 2487Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32Sofia, U. J., Cardelli, J. A., Guerin, K. P., & Meyer, D. M. 1997, ApJ, 482,L105Spitzer, L. 1978, in Physical processes in the interstellar medium. New YorkWiley-Interscience, 1978, ed. L. Spitzer (Weinheim, Germany)Staemmler, V., & Flower, D. R. 1991, Journal of Physics B: Atomic,Molecular and Optical Physics, 24, 2343Stanimirovi´c, S., Murray, C. E., Lee, M.-Y., Heiles, C., & Miller, J. 2014,ApJ, 793, 132Sternberg, A., & Dalgarno, A. 1995, ApJS, 99, 565Sternberg, A., Petit, F. L., Roueff, E., & Bourlot, J. L. 2014, ApJS, 790, 10STayal, S. S. 2008, AA, 486, 629Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1Tielens, A. G. G. M. 2013, Reviews of Modern Physics, 85, 1021Valdivia, V., Hennebelle, P., Gerin, M., & Lesaffre, P. 2015, AA, 587van Dishoeck, E. F., & Black, J. H. 1986, ApJS, 62, 109 hermal Phases – from Solar Metallicity to Primordial Gas 21hermal Phases – from Solar Metallicity to Primordial Gas 21