Thermal Phases of the Neutral Atomic Interstellar Medium -- from Solar Metallicity to Primordial Gas

TL;DR

This paper investigates how the thermal phases of the neutral atomic interstellar medium change across different metallicities, revealing the disappearance of multiphase structure at very low metallicities and implications for galaxy star formation.

Contribution

It provides a comprehensive model of the thermal structure of the neutral atomic ISM over a wide metallicity range, including new analytic expressions for phase transition thresholds.

Findings

Multiphase ISM exists at solar metallicity with coexisting warm and cold phases.

At low metallicity, cosmic-ray heating dominates, shifting phase transition conditions.

H$_2$ cooling becomes significant at metallicities below 1 extperthousand, smoothing out multiphase phenomena.

Abstract

We study the thermal structure of the neutral atomic (H {\small I}) interstellar medium across a wide range of metallicities, from supersolar down to vanishing metallicity, and for varying UV intensities and cosmic-ray ionization rates. We calculate self-consistently the gas temperature and species abundances (with a special focus on the residual H$_2$), assuming thermal and chemical steady-state. For solar metallicity, $Z' \equiv 1$, we recover the known result that there exists a pressure range over which the gas is multiphased, with the warm ($\sim 10^4$ K, WNM) and cold ($\sim 100$ K, CNM) phases coexisting at the same pressure. At a metallicity $Z' \approx 0.1$, the CNM is colder (compared to $Z'=1$) due to the reduced efficiency of photoelectric heating. For $Z' \lesssim 0.1$, cosmic-ray ionization becomes the dominant heating mechanism and the WNM-to-CNM transition shifts to ever increasing pressure/density as the metallicity is reduced. For metallicities $Z' \lesssim 0.01$, H$_2$ cooling becomes important, lowering the temperature of the WNM (down to $\approx 600$ K), and smoothing out the multiphase phenomenon. At vanishing metallicities, H$_2$ 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 different regimes, and the critical metallicities for H$_2$ cooling and heating. We discuss potential implications on the star-formation rates of galaxies and self-regulation theories.

Tables

Table 1. Table 1 Parameters

notation	fiducial value / range	meaning	comments
$Z^{\prime }$	from 10-5 to 3	metallicity relative to Galactic	see §2.2.1, §4
$I_{\mathrm{UV}}$	0.1, 1, 10	the 6-13.6 eV UV intensity, relative to Galactic	see §2.2.3, §5.2
${\zeta }_{-16}=\frac{{\zeta }_{cr}}{{10}^{-16}{s}^{-1}}$	0.1, 1, 10	the total (primary+secondary) ionization rate per H atom	see §2.2.2, §5.2
$A_C$, $A_O$	1.4 and 3.2 $\times {10}^{-4}$	the carbon and oxygen solar elemental abundances	see §2.2.1 and Eq. (4)
${\delta }_C$, ${\delta }_O$	0.53 and 0.41	the carbon and oxygen dust depletion factors	see §2.2.1 and Eq. (4)
$Z_d^{\prime }$		the DGR relative to Galactic	see §2.2.1 and Eq. (5)
${\zeta }_p$	${\zeta }_{cr}/\left(1+{\varphi }_s\right)$	the primary ionization rate per H atom	see §2.2.2 and Eq. (A9)
${\varphi }_s$	$\approx 0.6$	the number of secondary electrons produced per ionization	see §2.2.2 and Eq. (A9)
$D_0$	$5.8\times {10}^{-11}$ s-1	the H2 photodissociation rate at $I_{\mathrm{UV}}=1$	see §2.2.3, §5.2
$D_{-}$	$5.6\times {10}^{-9}$ s-1	the H- photodetachment rate at $I_{\mathrm{UV}}=1$	see §2.2.3, §5.2
$R$	(cm3  s-1)	the H2 formation rate coefficient	see §2.2.4
$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-3  s-1)	the cooling and heating rates per unit volume	see §2.1 and A.2
${ℒ}_i\left(T\right)$	(erg  cm3  s-1)	the cooling efficiency of species $i$. For $n/n_{\mathrm{crit}}\ll 1$, $L_i=x_i{n}^2{ℒ}_i\left(T\right)$	see §A.2
${\Lambda }_i\left(T\right)$	(erg  s-1)	the cooling rate per particle. For $n/n_{\mathrm{crit}}\gg 1$, $L_i=x_{i}n{\Lambda }_i\left(T\right)$	see §A.2

Table 2. Table 2 Reduced Network

Reaction	Rate coefficient (subscript refers to the reaction no. in the text)
$H+H:\mathrm{dust}\to H_2$	$R_d=3\times {10}^{-17}{T}_2^{1/2}{Z}_d^{\prime }$ cm3 s-1
$H+e\to H^{-}+\nu$	$k_{\text{<a><span>10</span></a>}}=7.2\times {10}^{-16}{T}_3^{0.64}{e}^{-9.2K/T}$ cm3 s-1
$H^{-}+H\to H_2+e$	$k_{\text{<a><span>11</span></a>}}=2.6\times {10}^{-9}{T}_3^{-0.39}{e}^{-39.4K/T}$ cm3 s-1
$H^{-}+H^{+}\to H+H$	$k_{\text{<a><span>15</span></a>}}=4.1\times {10}^{-8}{T}_3^{-0.5}$ cm3 s-1
$e+H^{+}\to H+\nu$	${\alpha }_B=1.4\times {10}^{-12}{T}_3^{-0.75}$ cm3 s-1
$H_2+H\to H+H+H$	Formula from Martin et al. (1996)
Reaction	Photo/CR-rate
$H+\mathrm{CR}\to H^{+}+e$	$\zeta ={10}^{-16}{\zeta }_{-16}$ s-1
$H_2+\nu \to H+H$	$D_0{I}_{\mathrm{UV}}=5.8\times {10}^{-11}{I}_{\mathrm{UV}}$ s-1
$H^{-}+\nu \to H+e$	$D_{-}{I}_{\mathrm{UV}}=5.6\times {10}^{-9}{I}_{\mathrm{UV}}$ s-1
C+, O, He	assume all gas-phase carbon/oxygen/helium are in C+, O, He

Equations

$$L(n,T;x_i)=G(n,T;x_i),$$

$$\sum _{j>l}{x}_j{x}_{l}n{k}_{jl,i}+\sum _j{x}_j{q}_{j,i}=0.$$

$$X_m=\sum _i{\alpha }_{i,m}{x}_i,$$

$$X_{\mathrm{m}}=A_{\mathrm{m}}{Z}^{\prime }\left(1-{\delta }_{\mathrm{m}}\times (Z_d^{\prime }/Z^{\prime })\right),$$

$$Z_d^{\prime }=\{\begin{array}{ll}{Z}^{\prime }& \text{for }{Z}^{\prime }\ge Z_0^{\prime }\\ Z_0^{\prime }(Z^{\prime }/Z_0^{\prime }{)}^{\alpha }& \text{for }{Z}^{\prime }<Z_0^{\prime },\end{array}$$

$$\mathrm{H}+\mathrm{cr}$$

$${\mathrm{H}}_2+\mathrm{cr}$$

$$R_d=3\times 10^{-17}{\left(\frac{T}{100\mathrm{K}}\right)}^{1/2}{Z}_d^{\prime }{\mathrm{cm}}^3{\mathrm{s}}^{-1}.$$

$$\mathrm{H}+\mathrm{e}\to {\mathrm{H}}^{-}+\nu ,$$

$${\mathrm{H}}^{-}+\mathrm{H}\to {\mathrm{H}}_2+\mathrm{e}.$$

$$x_{\mathrm{e}}={\left(\frac{\zeta }{{\alpha }_{B}n}\right)}^{1/2}=2.7\times 10^{-3}{T}_3^{3/8}{\zeta }_{-16}^{1/2}{\left(\frac{n}{10\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1/2},$$

$$R_{-}\equiv x_{\mathrm{e}}{k}_{\text{ref}reac:H-formation}\eta \approx 1.9\times 10^{-18}{T}_3^{1.02}\eta {\zeta }_{-16}^{1/2}{\left(\frac{n}{10\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1/2}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}$$

$${\mathrm{H}}^{-}+\nu \to \mathrm{H}+\mathrm{e}.$$

$${\mathrm{H}}^{-}+{\mathrm{H}}^{+}\to \mathrm{H}+\mathrm{H}.$$

$$\eta ={\left(1+\frac{D_{-}{I}_{\mathrm{UV}}}{k_{\text{ref}eq:H2formgasreac}n}+\frac{k_{\text{ref}reac:mutualneutralization}{x}_{{\mathrm{H}}^{+}}}{k_{\text{ref}eq:H2formgasreac}}\right)}^{-1}$$

$$x_{{\mathrm{H}}_2}=\frac{x_{\mathrm{e}}{k}_{\text{ref}reac:H-formation}n}{D_0{I}_{\mathrm{UV}}}\eta =3.3\times 10^{-7}{T}_3^{1.01}\eta {\zeta }_{-16}^{1/2}{I}_{\mathrm{UV}}^{-1}{\left(\frac{n}{10\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{1/2}.$$

$$n_{\mathrm{cool},Z}=\frac{{\zeta }_{\mathrm{p}}{E}_{\mathrm{cr}}}{A_{\mathrm{O}}{\mathcal{L}}_Z(T_0)Z^{\prime }}\approx 4.0\left(\frac{10^{-2}}{Z^{\prime }}\right){\zeta }_{-16}\mathrm{c}{\mathrm{m}}^{-3},$$

$$P_{\max }\simeq 1.1n_{\mathrm{cool},Z}{T}_0=2.6\times 10^4\left(\frac{10^{-2}}{Z^{\prime }}\right){\zeta }_{-16}\mathrm{c}{\mathrm{m}}^{-3}\mathrm{K}.$$

$$n{\mathcal{L}}_{{\mathrm{H}}_2}{x}_{{\mathrm{H}}_2}={\zeta }_p{E}_{\mathrm{cr}}.$$

$$n_{\mathrm{cool},{\mathrm{H}}_2}={\left(\sqrt{\frac{{\alpha }_B}{\zeta }}\frac{D_0{I}_{\mathrm{UV}}}{k_{\text{ref}reac:H-formation}}\frac{E_{\mathrm{cr}}{\zeta }_p}{{\mathcal{L}}_{{\mathrm{H}}_2}}\right)}^{2/3}\approx 5.0{\zeta }_{-16}^{1/3}{I}_{\mathrm{UV}}^{2/3}\mathrm{c}{\mathrm{m}}^{-3}.$$

$$x_{{\mathrm{H}}_2,\mathrm{cool}}=1.4\times 10^{-6}{\left(\frac{{\zeta }_{-16}}{I_{\mathrm{UV}}}\right)}^{2/3}.$$

$$P_{\mathrm{cool},{\mathrm{H}}_2}\approx 1.1n_{\mathrm{cool},{\mathrm{H}}_2}{T}_0\approx 3.3\times 10^4{\zeta }_{-16}^{1/3}{I}_{\mathrm{UV}}^{2/3}\mathrm{c}{\mathrm{m}}^{-3}\mathrm{K}.$$

$$Z_{\mathrm{cool},{\mathrm{H}}_2}^{\prime }\equiv \frac{{\mathcal{L}}_{{\mathrm{H}}_2}{x}_{{\mathrm{H}}_2}}{{\mathcal{L}}_z{A}_{\mathrm{O}}}\approx 8.0\times 10^{-3}{\left(\frac{{\zeta }_{-16}}{I_{\mathrm{UV}}}\right)}^{2/3}.$$

$$n_{\mathrm{heat},{\mathrm{H}}_2}$$

$$=7.2\times 10^3{\zeta }_{-16}^{1/3}{T}_3^{-1.08}\mathrm{c}{\mathrm{m}}^{-3},$$

$$Z_{\mathrm{heat},{\mathrm{H}}_2}^{\prime }=1.2\times 10^{-5}{\zeta }_{-16}^{2/3}.$$

$$P_{\min }=8500I_{\mathrm{UV}}\frac{Z_d^{\prime }/(Z^{\prime }f)}{1+3.1(I_{\mathrm{UV}}{Z}_d^{\prime }/{\zeta }_{-16}{)}^{0.365}}\mathrm{c}{\mathrm{m}}^{-3}\mathrm{K},$$

$$x_{{\mathrm{H}}_2}=\frac{x_{\mathrm{e}}{k}_{\text{ref}reac:H-formation}{k}_{\text{ref}eq:H2formgasreac}{n}^2}{D_0{D}_{-}{I}_{\mathrm{UV}}^2}=3.3\times 10^{-8}{T}_3^{0.62}{\zeta }_{-16}^{1/2}{I}_{\mathrm{UV}}^{-2}{\left(\frac{n}{10\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{3/2}.$$

$$L_i=x_i{n}^2\stackrel{\mathcal{L}(T)}{\overbrace{\sum _j{x}_j{q}_{ij,\uparrow }(T)\Delta E}}\frac{1}{1+n/n_{\mathrm{crit}}},$$

$$L_i\simeq x_{i}n{n}_{\mathrm{crit}}(T)\mathcal{L}(T)\equiv x_{i}n{\Lambda }_{\mathrm{LTE}}(T),$$

$$G_{\mathrm{pe}}={\gamma }_{\mathrm{pe}}ϵZ_d^{\prime }{I}_{\mathrm{UV}}n,$$

Full text

1. Introduction

Emission and absorption line studies of the 21 cm hyperfine spin-flip transition of atomic hydrogen (HI) show that cold ($T<200$ K) dense clumps are commonly embedded within extended warm ($T\sim 8000$ K) diffuse atomic clouds. The cold neutral medium (CNM) clumps are observed as narrow absorption line features embedded within broad emission line envelopes tracing the warm neutral medium (WNM) component (e.g., in Galactic observations: Dickey et al. 1978; Kulkarni & Heiles 2011; Colgan et al. 1988; Garwood & Dickey 1989; Kolpak et al. 2002; Dickey et al. 2003; Heiles & Troland 2003; Roy et al. 2013; Stanimirović et al. 2014; Patra et al. 2018; Murray et al. 2018, and in external galaxies: Braun & Walterbos 1992; Dickey & Brinks 1993; Mebold et al. 1997; Dickey et al. 2000; Marx-Zimmer et al. 2000; Warren et al. 2012; Bolatto et al. 2013; Herrera-Camus et al. 2017). Remarkably, for the characteristic thermal pressures prevailing in the interstellar gas, theory predicts that a multiphase medium is possible in which the CNM and WNM can coexist (Field et al., 1969; Wolfire et al., 1995; Liszt, 2002; Wolfire et al., 2003). The multiphase nature of the H I gas is also clearly observed in numerous hydrodynamical simulations 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 interstellar HI, the WNM is cooled by electron impact Ly$\alpha$ emission (Gould & Thakur, 1970; Dalgarno & McCray, 1972), and the cold medium by metal line fine structure cooling, mainly neutral impact [CII] 158 $\mu$m emission (Dalgarno & Rudge, 1964). For a given heating rate and at sufficiently low densities metal line cooling becomes very inefficient and the gas heats up to temperatures ($\sim 8000$ K) at which the Ly$\alpha$ emission limits further heating. At high densities the metal line cooling becomes effective and the temperature drops to values characteristic of the fine-structure energy level spacings ($\sim 100$ K). Because of the particular shape of the combined Ly$\alpha$ plus fine-structure cooling function (we discuss this in detail), an intermediate range of densities and associated pressures exists at which the gas can be either warm or cold.

Hydrodynamical simulations show that CNM forms at sites where flows of WNM collide and compress the gas (e.g., Hennebelle & Pérault, 1999; Kritsuk & Norman, 2002; Audit & Hennebelle, 2005; Vázquez-Semadeni et al., 2006; Gazol & Kim, 2013; Saury et al., 2014; Gazol & Villagran, 2016). The formation of the dense CNM phase and the multiphase structure 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 low pressure prevailing at the outskirts of galaxy discs do not allow the formation of the cold phase resulting in a shut-off of star-formation. The multiphase phenomena where CNM and WNM phases coexist at the same pressure may provide a 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-formation rate drastically increases, the high heating rate induced by the stellar UV prevents the formation of the CNM phase, which suppresses star-formation.

Previous computations of the HI WNM/CNM thermal phase structures have generally assumed that metals dominate the CNM cooling. Then, as would be expected, as the metal abundance is reduced higher particle densities are required to compensate for the reduced cooling, and the multiphase behavior is predicted to occur at higher gas densities and thermal pressures (Wolfire et al., 1995; Liszt, 2002; Walch et al., 2011). The central question we address in this paper is how WNM/CNM multiphase behavior of atomic hydrogen (HI) gas is affected by H2 cooling and heating processes as the metallicity becomes vanishingly small.

H2 cooling in the context of primordial and/or low-metallicity gas and first star-formation has been considered by many 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 that in the absence of metals, just small amounts of H2 can enable cooling to a few 100 K (see for example the recent study by Glover & Clark 2014) thereby controlling the fragmentation and dynamical properties of collapsing star-forming clouds in the low-metallicity regime. Less attention has been given to the specifics of WNM to CNM transitions and possible multiphase behavior at very low-metallicities. For example, in their study of protogalactic disks Norman & Spaans (1997) presented pressure versus density thermal phase diagrams for atomic gas over a range of redshifts and associated metallicities. In their models multiphase behavior disappears at high redshifts (see their Fig. 5). However it remains unclear what are the individual roles of the metallicity and radiation field when H2 processes begin playing a role, and how these processes affect the phase structure. Similarly, Inoue & Omukai (2015) carried out numerical hydrodynamical simulations of isobarically contracting gas clouds and studied the effects of H2 cooling on thermal instabilities. They found that at sufficiently low metallicities H2 acts to smooth out thermal instabilities (see their Fig. 2). But why?

The goal of our paper is to study the competing microphysical processes that control the phase structures of atomic gas, into the regime of vanishing metallicity. We address the following questions. How are the WNM and CNM phases affected by the onset of H2 cooling (and heating)? Is multiphase gas possible when H2 replaces metal fine-structure line cooling? If not, why not? In this study we present numerical computations of the (equilibrium) thermal states, with UV radiation intensity, cosmic-ray (or X-ray) ionization rate, and metallicity, as distinct parameters, and with the addition of the H2 processes in a step by step fashion. This enables an analytic description of the behavior which we also present.

The structure of our paper is as follows. In §2 we present the governing equations, and discuss the chemical and thermal processes we include. We begin with a discussion of the solar metallicity model in §3. We discuss low metallicity behavior in §4 and 5, down to effectively pristine H2 dominated gas. We characterize the critical densities or pressures above which the gas may cool below $\sim 10^{4}$ K and condense into intermediate ($\sim 500$ K) or cold ($\sim 50$ K) phase (depending on metallicity), and whether two-phase media can exist, as function of metallicity, UV intensity and CR ionization rate. We summarize in §6.

2. Model Ingredients

We compute the temperatures and gas phase chemical compositions of the neutral atomic hydrogen (HI) components of the interstellar medium (ISM) of galaxies assuming thermal equilibrium and chemical steady state. We use the results to construct pressure versus density phase diagrams for the HI gas. We consider standard heating and cooling mechanisms, with the addition of molecular hydrogen (H2) processes that become important and eventually dominant in the limit of low metallicity. We investigate the role of the assumed far-UV radiation intensities and cosmic ray-ionization rates on the coupled thermal and chemical properties of the HI gas. We study the behavior of the thermal phase structures for standard ISM conditions into the limit of vanishing metallicity.

2.1. Basic Equations

In thermal equilibrium,

[TABLE]

where $G$ and $L$ are the total heating and cooling rates per unit volume (erg cm*-3* s*-1*). The rates depend on the total hydrogen volume density $n$ (cm*-3*), the gas temperature $T$ (K), the dust-to-gas abundance ratio, and the gas phase abundances of the chemical species, $x_{i}\equiv n_{i}/n$, especially for C+, O, electrons, and H2. For heating we include far-UV photoelectric (PE) ejection of electrons from dust grains and polycyclic aromatic hydrocarbons (PAHs), and cosmic-ray (CR) and/or X-ray ionization of the hydrogen atoms. For cooling, we include HI Ly$\alpha$ line-emission excited by electron impacts, and fine-structure [CII] (157.7$\mu$m), [OI] (63.2$\mu$m, 145.5$\mu$m), and [CI] (609.1$\mu$m, 370.4$\mu$m) line emissions excited primarily by collisions with the neutral hydrogen atoms. We also include energy losses due to electron recombination with dust grains and PAHs, and due to gas-dust collisions. As discussed by Wolfire et al. (2003, hereafter, WMHT) these are the dominant heating and cooling mechanisms for thermally stable and unstable warm/cold ($\sim 10^{4}$ to $\sim 100$ K) interstellar HI for standard ISM conditions 111Metastable line emission (e.g. [OI] $6302$ ${\rm\AA}$) also cools the WNM, but its contribution is subdominant (WMHT, see also Fig 30.1 in Draine 2011)..

Importantly, we also incorporate heating and cooling processes involving H2. This includes heating via the release of binding energy during H2 formation, kinetic energy imparted to the H-atoms produced by H2 photodissociation, and internal H2 energy released to the gas via collisional deexcitation of UV excited states (we refer to this latter process as ”UV pumping heating”). Cooling is mainly by H2 ro-vibrational line-emissions excited by collisions with electrons and neutrals. We also include cooling by isotopic HD line emissions, and due to dissociation of the H2 in collisions with thermal electrons and H-atoms, although these are minor. In the Appendix we present detailed discussions and expressions for the various heating and cooling processes we have included in our computations. The roles and effects of H2 heating and cooling processes have not been considered in previous analytic studies of the multiphase structure of the neutral ISM (Field et al. 1969; Wolfire et al. 1995; Liszt 2002; WMHT). We will show how the H2 heating and cooling processes influence the thermal structures as the heavy element abundances become small.

Eq. (1) is coupled to the equilibrium rate equations for the chemical abundances,

[TABLE]

This is a set of $N$ equations, equal to the number of species in the chemical network. The first sum is over the formation and destruction rates (s*-1*) of species $i$ via two-body reactions involving species $j$ and $l$, where $k_{jl,i}$ are the temperature-dependent reaction rate coefficients (cm3 s*-1*). The second sum 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*-1*) for these processes. For $j=i$, the coefficients $k_{jl,i}$ and $q_{j,i}$ are negative and represent the destruction of species $i$ in reactions with species $l$, and by UV and CR dissociations and ionizations. For $j\neq i$ the coefficients are positive and represent the formation channels for species $i$. Eqs. (2) are augmented by the (dependent) conservation equations for the total elemental abundances and electric charge. These are

[TABLE]

where $X_{m}$ is the total elemental abundance of element $m$, and $\alpha_{i,m}$ is the number of atoms of element $m$ in species $i$. Eq. (3) also describes the electric charge conservation equation, with $\alpha_{i,m}$ representing the number of positive charges in species $i$, and with $X_{m}=0$. In total, there are $M+1$ conservation equations, where $M$ is the number of elements in the chemical network. We discuss our chemical network in §2.3.

2.2. Parameters

The gas temperature $T$, and species abundances $x_{i}$, depend on the gas density, $n$, the far-UV radiation intensity, $I_{\rm UV}$, the CR ionization rate, $\zeta_{\rm cr}$, and the metallicity $Z^{\prime}$. We define these parameters below. Additional physical parameters affecting the 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 heavy elements, and the dust abundance, are crucial parameters for the thermal balance. The gas phase carbon and oxygen abundances determine the efficiency of metal-line cooling, and the dust-to-gas mass ratio determines the efficiency of PE heating and dust recombination cooling. We assume that the total gas-phase abundances of the carbon and oxygen species are given by

[TABLE]

where $m=$C, or O. In these expressions $Z^{\prime}$ is the normalized metallicity such that $Z^{\prime}=1$ corresponds to the Solar elemental abundances for which $A_{\rm O}=5.4\times 10^{-4}$ and $A_{\rm C}=3.0\times 10^{-4}$ (Asplund et al., 2009), $\delta_{\rm C}$ and $\delta_{\rm O}$ are dust depletion factors, and $Z^{\prime}_{d}$ is the normalized dust-to-gas ratio. We adopt

[TABLE]

with $Z_{0}^{\prime}=0.2$, and $\alpha=3$. This broken power-law relation is suggested by recent observations of low metallicity galaxies (Rémy-Ruyer et al., 2013). Here, $Z^{\prime}_{d}=1$ corresponds to a standard ISM dust-to-gas ratio of $6.2\times 10^{-3}$ (Zubko et al., 2003, Table 6, the BARE-GR-S model).

For the dust depletion factors, we adopt $\delta_{\rm C}=0.53$ and $\delta_{\rm O}=0.41$. With these values, we recover the characteristic interstellar Galactic carbon and oxygen gas-phase abundances $X_{\rm C}=1.4\times 10^{-4}$, $X_{\rm O}=3.2\times 10^{-4}$ at $Z^{\prime}=Z^{\prime}_{d}=1$ (Cardelli et al., 1996; Sofia et al., 1997; Meyer et al., 1997). On the other hand, at very low metallicities, when $Z^{\prime}_{d}/Z^{\prime}\ll 0.2$, the dust depletion terms vanish, and all the available carbon and oxygen remain in the gas-phase, with $X_{\rm C}=A_{\rm C}Z^{\prime}$, $X_{\rm O}=A_{\rm O}Z^{\prime}$.

2.2.2 Cosmic-ray and X-ray ionization

Cosmic rays (or X-rays) ionize atomic and molecular hydrogen,

[TABLE]

producing energetic non-thermal electrons. The electrons generated by these primary ionizations lose energy through coulomb scatterings that heat the gas, and through further atomic and molecular excitations and secondary ionizations (Dalgarno et al., 1999). The total (primary+secondary) ionization rate per hydrogen atom is $\zeta_{\rm cr}=\zeta_{p}\times(1+\phi_{s})$, where $\zeta_{p}$ is the primary ionization rate and $\phi_{s}\sim 0.7$ is the number of secondary electrons produced per primary ionization. The cosmic-ray (CR) ionization rate determines the CR heating rate and the electron abundance (see Eqs. 12 and A8, below). For typical Galactic values, $\zeta_{\rm cr}\approx 10^{-16}$ s*-1* (Indriolo & McCall, 2012; Tielens, 2013), and we define the normalized CR ionization rate, $\zeta_{-16}\equiv\zeta_{\rm cr}/(10^{-16}\ {\rm s^{-1}})$. Our fiducial value is $\zeta_{-16}=1$. We study the effects of varying $\zeta_{-16}$ in §5.2.

X-rays can play an important role in providing an additional source of ionization, affecting the chemistry and the heating rate. We do not explicitly include X-ray ionization in this study, however, these may be simulated by choosing an appropriate values for $\zeta_{\rm cr}$ and $\phi_{s}$.

2.2.3 Interstellar far-UV radiation field

The far-UV radiation intensity determines the efficiency of photoelectric heating, the H2 photodissociation rate, and the production rate of the C+ coolant. Draine (1978, 2011) estimated the local UV interstellar radiation field (ISRF), and obtained $u_{\rm UV,0}=0.056$ eV cm*-3*, for the energy density over the 6-13.6 eV band (=1.69 times the Habing 1968 estimate). We adopt the Draine spectral shape and allow for variations in the UV intensity by multiplying the overall spectrum with a normalized intensity parameter, $I_{\rm UV}$, such that $u_{\rm UV}=u_{\rm UV,0}I_{\rm UV}$. The UV intensity determines the PE heating rate through Eq. (A3). For a unit Draine ISRF, the free-space (unattenuated) H2 photodissociation rate is $D_{0}=5.8\times 10^{-11}$ s*-1* (Sternberg et al., 2014). At low metallicity, a harder spectrum may be more typical, since the first generations of stars (pop-III stars) are theorized to be very massive, with stellar effective temperatures as high as 105 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 in removing the H- anion via photodetachment (Miyake et al., 2010). H- is a critical intermediary for H2 formation at low metallicity (see §2.2.4, below). Mathis et al. (1983) estimated the Galactic ISRF within this band as a combination of diluted blackbodies at 3000, 4000 and 7500 K. For an ISRF combined of the Draine field over the 6-13.6 eV band, and with the Mathis field over the 0.75-6 eV band, the H- photodetachment rate is $D_{-}=2.7\times 10^{-7}$ s*-1*. In contrast, for a 105 K blackbody (BB) spectrum (normalized such that the photon density in the 6-13.6 eV band equals that of the Draine field), $D_{-}=2.8\times 10^{-9}$ s*-1* (Bialy & Sternberg 2015; see their Table 2). For our fiducial models we seek a compromise between these two extreme cases and extrapolate the Draine (1978) spectrum (see their Eq. 11) down to 0.75 eV. The resulting H- photodetachment rate is $D_{-}=5.6\times 10^{-9}$ s*-1*. This value is close to the 105 K BB value which represents the pristine metallicity gas, but is $\sim 2$ orders of magnitude smaller than the Draine+Mathis value, which represents the Galaxy and a high metallicity ISM. For moderately low metallicities $D_{-}$ may range between these two extremes, and will depend on the form of the initial mass function (IMF). We explore the effect of variations in $I_{\rm UV}$ and $D_{-}$ in §5.2.

The free-field photodissociation rate, $D_{0}$ may be strongly reduced in the interiors of a cloud due to dust absorption and H2 self-shielding, leading to partial or complete H I-to-H2 conversion (van Dishoeck & Black, 1986; Draine & Bertoldi, 1996; McKee & Krumholz, 2010; Sternberg et al., 2014; Bialy & Sternberg, 2016). In this paper we focus on the thermal properties of the global atomic ISM, rather than the structure of the H I-to-H2 transition or the interior shielded molecular gas. Thus, we always use the free-space UV field and do not include radiation attenuation in our models. For the WNM shielding is in any case not significant (except for extremely large shielding columns). For the CNM, our model applies to the cloud outskirts that are exposed to the free-field UV.

2.2.4 H2 Formation

At sufficiently high metallicity, the H2 is formed predominantly on dust grains. We adopt

[TABLE]

Sophisticated models of the H2 formation process on dust-grains (e.g., Cazaux & Tielens, 2002; Le Bourlot et al., 2012; Bron et al., 2014) derived more complicated dependences on gas temperature, as well as on other physical parameters (e.g., through the dust temperature which in turn depends on the UV radiation intensity). However for our purposes such a simple description is sufficient, because as we show below, H2 affects the thermal structure only at low metallicities where its formation is dominated through gas-phase processes.

At low metallicities ($Z^{\prime}\lesssim 0.1$), H2 is formed mainly through gas-phase reactions, with the most important one being the H- formation route (McDowell, 1961; Peebles & Dicke, 1968; Hirasawa, 1969; Galli & Palla, 1998). Although, the H2 formation rate is obtained through the numerical solution to the chemical network (as discussed in §2.1), it is useful to obtain approximate analytic equations for the effective gas-phase formation rate coefficient, and the resulting H2 fractional abundance.

This formation channel is initiated by radiative attachment

[TABLE]

followed by associative detachment

[TABLE]

Assuming that the electron abundance is set by a balance between CR ionization of the atomic hydrogen and proton-electron recombination, we have

[TABLE]

where $T_{3}\equiv T/({\rm 10^{3}\ K})$ and $\alpha_{B}$ is the case-B recombination rate coefficient. The effective H2 formation rate coefficient through this channel is

[TABLE]

where $k_{\ref{reac: H- formation}}$ is the radiative attachment rate coefficient222The rate coefficients are based on UMIST 2012 (McElroy et al., 2013). See Table 2 for a summary of the most important reactions and rates. See online material for the full network. (Eq. 10). The H2 formation sequence is moderated mainly by photodetachment,

[TABLE]

and also by mutual neutralization with protons,

[TABLE]

In Eq. (13) we define the branching ratio

[TABLE]

for H2 formation via associative detachment versus H- photodetachment and mutual-neutralizaion. For $D_{-}=5.6\times 10^{-9}$ s*-1* photodetachmnet becomes important when $n/I_{\rm UV}\lesssim 1$ cm*-3*, while mutual-neutralization is important when $n/\zeta_{-16}\lesssim 0.01$ cm*-3*. As we show below, H2 plays an important role in cooling and heating at densities $\gg 1$ cm*-3*. Thus, generally, $\eta\approx 1$ is a good approximation 333For high $D_{-}$ values, $\eta$ may become smaller than unity at the cooling point, reducing the H2 abundance. We discuss this in §5.2..

Assuming that H2 destruction is dominated by photodissociation, we obtain an analytic approximation for the H2 fractional abundance at low metallicities,

[TABLE]

We find Eq. (17) to be in good agreement with the numerical solution for the chemical network, over the gas densities and temperatures in our parameter space.

2.3. Numerical Method

Given a set of parameters, $I_{\rm UV}$, $\zeta_{\rm cr}$, and $Z^{\prime}$, we solve for the equilibrium temperature and species abundances as functions of the hydrogen density, $n$. We form a logarithmic grid, spanning $\log T=1$ to $4.2$ and $\log n=-2$ to 6, with $100\times 90$ sampling points. For each volume density in the grid, we calculate the species abundances coupled with the net cooling, $L-G$, as functions of the temperature, and find the temperature(s) for which 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*-3* K) by simply multiplying by the density, i.e., $P(n)=nT(n)$. For any 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 294 two-body gas-phase reactions among 34 species, composed of hydrogen, helium, carbon, and oxygen bearing atomic and molecular species. We adopt the rate coefficients suggested by the UMIST 2012 database (McElroy et al., 2013), except for the H2 collisional dissociation for which we adopt the Martin et al. (1996) rate (including both collisional dissociation and dissociative tunneling). We also include the formation of H2 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 in Bialy & Sternberg 2015, see their table 2.

In total we obtain a set of 34 coupled non-linear algebraic equations for the 34 species. The equations depend on $T$ and $n$, as well as on the additional physical parameters, $I_{\rm UV}$, $\zeta_{\rm cr}$, and $Z^{\prime}$. We replace five of these equations with the element conservation (H, He, C, and O), and charge conservation equations. For any density $n$ there is always a single unique solution for the temperature $T$ and abundances $x_{i}$. Rather than solving the equilibrium algebraic equations (e.g., via Newton-Raphson iteration), we instead integrate the time-dependent rate equations until a time, $t$, where $dx_{i}/dt=0$ and equilibrium is reached. We find this approach to be more robust.

To improve computational speed, we also consider a reduced chemical network, which allows us to calculate the steady-state temperature as a function of metallicity and volume densities (or pressures) over a very well sampled grid, as presented in Figs. 6 and 7. The reduced network is discussed in the Appendix.

3. Thermal Structure at Solar Metallicity

In this section we first present results for our solar metallicity reference models ($Z^{\prime}=Z^{\prime}_{d}=1$). The solid curves in Fig. 1 show the temperature (upper left panel), the thermal pressure (lower left panel), and the cooling and heating rates per hydrogen nucleus (right panels), as functions of the hydrogen nucleus gas density $n$, assuming $I_{\rm UV}=\zeta_{-16}=1$. The dashed curves in the left panels, located above and below the solid curve, are for $I_{\rm UV}=\zeta_{-16}=10$ and 0.1, respectively.

First we focus on the behavior at low to intermediate densities, $n\lesssim 10^{3}$ cm*-3*, where H2 cooling and heating processes are subdominant. In this regime, our results essentially reproduce those presented by Wolfire et al. (1995) and WMHT. The behavior is controlled by the temperature dependence of the Ly$\alpha$ and fine-structure metal-line cooling functions. We display these cooling functions in the upper panel of Fig. 3. The black curve shows the exponentially steep electron impact excitation Ly$\alpha$ cooling rate (erg s*-1*) per hydrogen nucleus. The blue curves show the much flatter metal cooling rates for $Z^{\prime}=1$, 10*-2*, and 10*-4*, assuming all of the gas-phase carbon is C+, and the oxygen is in neutral atomic form. The metal cooling is controlled by electron and neutral-hydrogen impact excitation of the fine-structure levels. The solid and dashed curves in Fig. 3 are for $x_{\rm e}=10^{-4}$ and 10*-2* respectively. For the high $x_{\rm e}$, and when $T\lesssim 500$ K, metal line cooling is dominated by electron impact excitation. The shapes of these curves are of the utmost importance.

Cooling becomes inefficient compared to heating as the density becomes small. For $n\lesssim 0.1$ cm*-3*, the gas becomes warm, and Ly$\alpha$ dominates the cooling. In this regime the temperature depends weakly on density and remains close to 104 K. This is because for $T\lesssim 10^{4}$ K Ly$\alpha$ cooling is exponentially sensitive to the temperature (see Fig. 3), and any increase in the cooling rate due to an increase in gas density requires just a small decrease in temperature to balance the heating rate. At $n\approx 1$ cm*-3*, the temperature decreases to $\sim 6\times 10^{3}$ K, at which point energy losses by the [CII] and [OI] emissions become comparable to Ly$\alpha$, as seen in Fig. 1 and Fig. 3. Above this density the metals dominate the cooling and the temperature falls abruptly. The sharp drop occurs because at these high temperatures (and as long as $T\gtrsim 300$ K, i.e. above the transition energies of the fine-structure lines) metal cooling is a weak function of temperature (Fig. 3). Any increase in density then requires a large decrease in temperature to balance the heating rate. We refer to the density at which metal cooling takes over as the “cooling point density”, and we provide analytic estimates in §4. Below a few 100 K, the metal line cooling becomes exponentially temperature-sensitive (Fig. 3), and the equilibrium $T(n)$ curve flattens again.

In Fig. 1, the low-density flat part of the $T(n)$ curve followed by the sharp temperature drop gives rise to a local maximum, $P_{\rm max}$, in the $P(n)$ curve. In general, a local maximum such as this occurs if and where the temperature starts decreasing more rapidly than $1/n$. For our fiducial model with $I_{\rm UV}=\zeta_{-16}=1$, we find $P_{\rm max}=7.6\times 10^{3}$ cm*-3* K occurring at a density $n_{\rm max}=1.1$ cm*-3*. For $n<n_{\rm max}$ the gas is warm with $T_{\rm WNM}=6000$ to $10^{4}$ K. This is the classical warm neutral medium (WNM), and $P_{\rm max}$ is the maximum pressure possible for the WNM. Starting from cold gas at high densities, the steep rise of the $T(n)$ curve as the density is *reduced

• gives rise to a local minimum, $P_{\rm min}$, in the pressure versus density curve. In general, a local minimum occurs if and where the temperature starts rising more rapidly than $1/n$ as the density is reduced. For our fiducial model $P_{\rm min}=2.6\times 10^{3}$ cm*-3* K at $n_{\rm min}=7.3$ cm*-3*. For $n>n_{\rm min}$, the gas remains cold with $T_{\rm CNM}\lesssim 300$ K. This is the cold neutral medium (CNM) and $P_{\rm min}$ is the minimum pressure possible for the CNM. We show $P_{\rm min}$ and $P_{\rm max}$ in lower (upper) panels of Fig. 1 as the pair of thin 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$\alpha$ to metal-line cooling the temperature falls more rapidly than $1/n$ over a significant density interval from $n_{\rm max}$ to $n_{\rm min}$. This enables multiple solutions for the gas temperature between $P_{\rm min}$ and $P_{\rm max}$ for densities ranging from $n_{\rm min}T_{\rm CNM}/T_{\rm WNM}$ to $n_{\rm max}T_{\rm WNM}/T_{\rm CNM}$. For densities between $n_{\rm max}$ and $n_{\rm min}$ the gas is thermally unstable and this is the UNM. For isobaric density perturbations $\delta n$ in the UNM, the gas either cools to the stable CNM branch (for positive perturbations), or heats to the stable WNM branch (for negative perturbations). For thermal pressures between $P_{\rm max}$ and $P_{\rm min}$ the gas may become multiphased with CNM condensations coexisting within an enveloping WNM at the same thermal pressure.

Fig. 1 shows that for $n\gtrsim 5\times 10^{3}$ cm*-3*, the temperature increases with density, as H2 heating processes start contributing. In this regime the H2 heating rate increases more steeply than $n^{2}$ (see Eq. A15 in the limit $n/n_{\rm crit}<1$) leading to the temperature increase with $n$. We stress that the H2 heating processes play this role so long as the gas is primarily atomic. If any self-shielding induces an HI-to-H2 transition, the H2 heating will be strongly suppressed. For solar metallicity a total (atomic plus molecular) column as small as $N\sim 10^{17}$ cm*-2* is sufficient to induce conversion to molecular form for $n\gtrsim 10^{4}$ cm*-3* (and $I_{\rm UV}=1$). However, at low metallicities, much large columns, $N\sim 10^{21}$ cm*-2*, are required (Bialy & Sternberg, 2016)

At still higher densities (beyond the considered parameter space) the hydrogen undergoes a conversion into H2 and the atomic carbon and oxygen convert into CO and H2O. These molecular formation processes become efficient as the temperature exceeds few hundred Kelvins, which triggers the neutral-neutral chemical network of H2O (Sternberg & Dalgarno, 1995; van Dishoeck et al., 2013; Bialy et al., 2015). As our focus is on the (predominantly) atomic ISM, we restrict ourselves to densities $n\lesssim 10^{5}$ cm*-3*.

The dashed curves in the left panel of Fig. 1 are the equilibria for $I_{\rm UV}=0.1$ and 10, assuming a constant ratio $\zeta_{-16}/I_{\rm UV}=1$. For all values of $I_{\rm UV}$ in Fig. 1 photoelectric heating dominates and is proportional to $I_{\rm UV}$. Thus, when $I_{\rm UV}$ is increased (decreased) by a factor of 10, the $T(n)$ curves shift to the right (left), by a factor of 10. As a result, the pressure range within which the gas may be multiphased shifts to higher or lower values in proportion to $I_{\rm UV}$. For example, for $I_{\rm UV}=10$, $(P_{\rm max},P_{\rm min})=(7.7,2.9)\times 10^{4}$ cm*-3* K. For $I_{\rm UV}=0.1$, $(P_{\rm max},P_{\rm min})=(7.5,2.7)\times 10^{2}$ cm*-3* K.

At solar metallicity, H2 cooling and heating become important only in the very dense CNM, but as we discuss below, at lower metallicities, H2 cooling starts to play a role also at low-to-intermediate densities, in the WNM, affecting the WNM-to-CNM transition.

4. Thermal Structure at Low Metallicity

In this section we investigate the thermal structures as the metallicity is lowered. We first consider the transition from photoelectric to cosmic-ray dominated heating as $Z^{\prime}$ is reduced, for a range of ratios $\zeta_{-16}/I_{\rm UV}$. We then solve the equation of thermal equilibrium, in two steps. First, we lower the metallicity and solve for $T$ with the exclusion of the H2 heating and cooling processes. Then we solve again with the H2 included. This procedure enables us to distinguish the effects of diminished metal line cooling from the appearance of the H2 heating and cooling processes that eventually dominate the behavior at very low metallicity.

4.1. PE versus CR

Near solar metallicity ($Z^{\prime}\sim 1$) photoelectric (PE) heating usually dominates, and the PE heating rate decreases in proportion to the dust-to-gas ratio $Z_{d}^{\prime}$. At sufficiently low $Z^{\prime}$ heating by cosmic-ray ionization becomes competitive and dominates. In Fig. 2 we show the boundary lines at which the PE and CR heating rates are just equal, in the $I_{\rm UV}/\zeta_{-16}$ versus $Z^{\prime}$ plane, for $n/\zeta_{-16}$ ratios equal to 0.1, 1, and 10 cm*-3*. For $\zeta_{-16}/I_{\rm UV}=1$ the transition from PE to CR heating occurs between $Z^{\prime}=0.3$ and 0.1, depending on density and temperature. With decreasing $\zeta_{-16}/I_{\rm UV}$, the transition metallicity is lowered, but not by much. This is because below $Z^{\prime}=0.2$ the dust-to-gas ratio falls super-linearly with metallicity (see Eq. 5). Conversely, for sufficiently large $\zeta_{-16}/I_{\rm UV}$ ($\gtrsim 100$), CR heating may dominate even at solar metallicity.

4.2. Step-1: Reducing the Metallicity

We now consider the thermal structures at reduced metallicity, with the exclusion of the H2 heating and cooling processes in Step-1. Results for $T(n)$ and $P(n)$ are shown as the dashed curves in Fig. 4 (top and middle panels) for $Z^{\prime}$ from 3 to $10^{-5}$. In these computations we have assumed $I_{\rm UV}=\zeta_{-16}=1$.

Fig. 4 shows that as $Z^{\prime}$ is lowered to $0.1$, ($P_{\rm max}$,$P_{\rm min}$) drop to $(4.6,0.43)\times 10^{3}$ cm*-3* K. These drops are due to the more rapid reduction in $Z^{\prime}_{d}$ and associated PE heating rate, compared to the diminished metal line cooling rate that decreases with $Z^{\prime}$. The reduced heating relative to cooling at $Z^{\prime}=0.1$, reduces the minimum CNM density by a factor of 2.6 (compared to $Z^{\prime}=1$), to $n_{\rm min}=2.8$ cm*-3*. The maximal density of the WNM is also lower, with $n_{\rm max}=0.71$ cm*-3*. As the metallicty is reduced further CR heating takes over, and the gas then remains warm to ever higher densities. This is expected given the reduction in the cooling efficiencies for a fixed CR heating rate that is independent of metallicity. Importantly, multiphase behavior is maintained, with $P_{\rm min}$ and $P_{\rm max}$ shifted to higher pressures.

The behavior in the CR heating regime may be analyzed straightforwardly as follows. As we discussed in §3, as the density is increased in the warm phase, metal line-emission energy losses start contributing at a temperature of $\sim 6\times 10^{3}$ K. The cooling point temperature is nearly independent of $Z^{\prime}$, and we denote the corresponding density as $n_{{\rm cool},Z}$ 444Here $n_{{\rm cool},Z}$ is the transition point for Ly$\alpha$-to-metal cooling. Later we also introduce $n_{\rm cool,H_{2}}$ for the transition point from Ly$\alpha$-to-H2 cooling.. For $n>n_{{\rm cool},Z}$, the temperature drops sharply as Ly$\alpha$ becomes inefficient and metal cooling dominates. To obtain an analytic formula for $n_{{\rm cool},Z}$ we equate the heating and metal line cooling rates at $T_{0}\equiv 6\times 10^{3}$ K. We assume heating by CR ionization. This gives

[TABLE]

where we have defined the Solar metallicity fine-structure line cooling efficiency (erg s*-1* cm3) per oxygen atom $\mathcal{L}_{Z}\equiv\mathcal{L}_{\rm[OI]}+(x_{\rm C^{+}}/x_{\rm O})\mathcal{L}_{\rm[CII]}$ 555As shown in Figs. 1 and 5, [CII] always dominates over [CI] cooling. . In Eq. (18) $\zeta_{p}$ is the primary ionization rate, and $E_{\rm cr}$ is the mean deposition energy per ionization (see Eq. A8). In the numerical evaluation we used $\mathcal{L}_{Z}(T_{0})=7.2\times 10^{-23}$  erg s*-1* cm3, $x_{{\rm C^{+}}}/x_{\rm O}=A_{\rm C}/A_{\rm O}=0.56$, $A_{\rm O}=5.4\times 10^{-4}$ (assuming no depletion), and $E_{\rm cr}=15.2$ eV, $\phi_{s}=0.57$. These values are appropriate for $x_{\rm e}=8.2\times 10^{-3}$, the equilibrium electron fraction at the transition point at $Z^{\prime}=0.01$. For $Z^{\prime}=10^{-2}$ Eq. (18) gives, $n_{{\rm cool},Z}=4.0$ cm*-3*, decreasing linearly with $Z^{\prime}$, in good agreement with the temperature drop points shown for $T(n)$ in Fig. 4 (dashed curves). For $Z^{\prime}\gtrsim 0.1$ PE heating contributes in addition to CR, and the cooling densities are larger than given by Eq. (18).

When the H2 processes are neglected (or negligible) the sharp drop in temperature induced by the transition from Ly$\alpha$ to metal line cooling gives rise to the multiphased $P(n)$ curves, and then $n_{{\rm cool},Z}$ is approximately equal to $n_{\rm max}$ the maximum density for which a WNM is possible. Thus,

[TABLE]

We assume a cosmological helium abundance, $A_{\rm He}=0.1$. Eq. (19) is consistent with our numerical results in Fig. 4.

The scalings given by Eqs. (18) and (19) are in agreement with previous models (Wolfire et al., 1995; Liszt, 2002, WMHT). As we now discuss, when H2 is included this simple picture is greatly altered.

4.3. Step-2: Inclusion of H2 Processes

The red curves in the upper panel of Fig. 3 show the H2 cooling function. In the lower panel we show the contribution of individual H2 rotational and vibrational levels to H2 cooling, for excitation by atomic hydrogen. As we discuss below, the fact that (unlike metal cooling) the H2 cooling function never saturates with temperature all the way until Ly$\alpha$ becomes operational, smooths out of the multiphase phenomena and eventually eliminates it entirely at low metallicity.

The solid curves in Fig. 4 show our results for the equilibrium $T(n)$ and $P(n)$ with the H2 processes included in our Step-2. Again, these are heating via H2 formation, photodissociation, and UV-pumping, and cooling via ro-vibrational line emissions, and collisional dissociation. The various heating and cooling rates for $Z^{\prime}=10^{-1,-2...,-5}$, are shown in Fig. 5.

The efficiencies of the molecular processes depend on the H2 abundance, $x_{\rm H_{2}}$, that we compute self-consistently with the temperature, for a given gas density and metallicity. In Fig. 4 we show $x_{\rm H_{2}}$ versus $n$, for $Z^{\prime}$ from 3 to 10*-5*. At all metallicities, H2 destruction by photodissociation is very efficient and the gas remains predominantly atomic. The H2 abundance is higher at high metallicities as the H2 is formed efficiently on dust-grains. For $Z^{\prime}\lesssim 0.1$, the H2 is formed much less efficiently in the gas phase. The resulting H2 abundances then range from $10^{-7}$ to $10^{-5}$, consistent with Eq. (17). Interestingly, it is in the low metallicity limit where $x_{\rm H_{2}}$ is very low, that H2 plays the most significant role in heating and cooling the gas. At high $Z^{\prime}$ molecular hydrogen is more abundant but metal line cooling is very efficient and induces a WNM-to-CNM transition at densities well below the densities at which H2 cooling starts playing a role. However, at sufficiently low $Z^{\prime}$, H2 cooling must eventually dominate over metals.

4.3.1 The onset of H2 cooling at low metallicity

When the H2 processes are included the $T(n)$ and $P(n)$ curves in Fig. 4 are altered in several ways. First, at sufficiently low metallicity it is H2 rather than the metals that compete with Ly$\alpha$ at the “cooling point” temperature of $\sim 6\times 10^{3}$ K below which Ly$\alpha$ becomes inefficient. This is also illustrated in the left panels of Fig. 5. We define $Z^{\prime}_{\rm cool,H_{2}}$ as the critical metallicity at which this switch from metals to H2 occurs at the cooling point. Numerically we find $Z^{\prime}_{\rm cool,H_{2}}\approx 10^{-2}$. The cooling point density is then shifted to a low value of $\approx 5$ cm*-3*, independent of metallicity (see Figs. 4 and 5). The H2 abundance at this point is $x_{\rm H_{2}}\approx 10^{-6}$.

For an analytic estimate we define the gas density for the H2 cooling point where

[TABLE]

We again assume heating by cosmic-ray ionization and $\mathcal{L}_{\rm H_{2}}$ is the H2 cooling efficiency (erg s*-1* cm3). With Eq. (17) for $x_{\rm H_{2}}$ with $\eta=1$, we obtain

[TABLE]

Plugging back into Eq. (17) we get the H2 abundance at the cooling point,

[TABLE]

In the numerical evaluations in Eqs. (21 and 22) we use $\mathcal{L}_{\rm H_{2}}=2.1\times 10^{-22}$ erg cm3 s*-1* at $T_{0}$, and $E_{\rm cr}$ and $\phi$, as in Eq. (18). The corresponding pressure is

[TABLE]

Unlike the Ly$\alpha$-to-metal transition which is accompanied by a maximum in the $P-n$ curve (Eq. 19), the Ly$\alpha$-to-H2 transition does not induce such a maximum but only a kink in the $P-n$ curve (see Fig. 4).

The critical metallicity below which H2 cooling becomes important is

[TABLE]

As we discuss in §5, this formula is in good agreement with our numerical results. Importantly, the temperature decreases more gradually for transitions from Ly$\alpha$ to H2 cooling. Because of the contributions from many lines in the ro-vibrational ladder (as illustrated in the bottom panel of Fig. 3), the H2 cooling function remains sensitive to temperature in the entire range $10^{2}<T<10^{4}$ K. This is illustrated in the upper panel of Fig. 3 where the red curves show the much steeper H2 cooling function compared to the flat metals function (blue curves). For this reason, the temperature drop at $n_{\rm cool,H_{2}}$ is moderated and does not produce a local maximum in the $P(n)$ curve. Most generally, the transition from Ly$\alpha$ to H2 cooling does not result in a multiphased behavior.

While H2 cooling does not induce a multiphase, Fig. 4 shows that multiphased behavior continues to be produced by metal-line cooling, although over a much narrower pressure range, down to a metallicity as low as $10^{-4}$. This multiphase is induced by the transition from H2 to metal cooling (rather than Ly$\alpha$ to metals) and occurs for a lower cooling point temperature of $T_{0}\approx 600$ K. So long as the heating is by CR ionization the cooling point density at which the sharp temperature drop occurs is still given by Eq. (18) but evaluated at the lowered $T_{0}$. Thus, in this regime $n_{{\rm cool},Z}$ continues to vary as $1/Z^{\prime}$, but the density prefactor in Eq. (18) is reduced. For example, for $Z^{\prime}=10^{-3}$, $n_{{\rm cool},Z}=40$ cm*-3*, leading to a pressure maximum of $P_{\rm max}\approx 1.1n_{{\rm cool},Z}T_{0}=7.6\times 10^{4}$ cm*-3* K for the WNM (see Fig. 4).

Because H2 cooling reduces the temperature of the WNM, the multiphase strip is now narrower than at solar metallicity, and continues to shrink with decreasing metallicity. For example, while at solar metallicity the multiphase range is wide, with $\delta_{n}\equiv(n_{\rm max}-n_{\rm min})/n_{\rm min}=5.4$ and $\delta_{P}\equiv(P_{\rm max}-P_{\rm min})/P_{\rm min}=1.9$, for $Z^{\prime}=10^{-3}$ we have $\delta_{n}=3.8$ and $\delta_{P}=0.88$. For $Z^{\prime}=10^{-4}$, the multiphase phenomena is already extremely suppressed, with $\delta_{n}=1.6$ and $\delta_{P}=0.05$.

4.3.2 The onset of H2 heating at extremely low metallicity

For still lower metallicities ($Z^{\prime}\leq 10^{-5}$) sharp gas cooling and associated multiphased structures disappear entirely. The density at which metal cooling can balance CR heating increases as $1/Z^{\prime}$, and at sufficiently low metallicities this density becomes so high that H2 heating becomes operational and dominant. In this limit metal cooling becomes ineffective. This is because H2 heating includes a dependence on the molecular formation rate (see Eq. A15) so that the heating rate increases more rapidly than $n^{2}$. However, metal cooling varies no faster than $n^{2}$, and so cannot compete. In this limit H2 is the dominant coolant, multiphase behavior does not occur, and the gas is driven to intermediate, thermally stable, temperatures of $\sim 600$ K. This behavior is illustrated in Figs. 4 and 5.

We obtain an analytic estimate for the critical metallicity, $Z^{\prime}_{\rm heat,H_{2}}$, below which metal line cooling is no longer effective as follows. We first estimate the density $n_{\rm heat,H_{2}}$ at which H2 and CR heating are equal. For $n>n_{\rm heat,H_{2}}$ we are in the H2 heating regime. The critical $Z^{\prime}_{\rm heat,H_{2}}$ is then the metallicity at which the metal cooling density, $n_{{\rm cool},Z}$, as given by Eq. (18), is at least as large as $n_{\rm heat,H_{2}}$. For $n_{{\rm cool},Z}\gtrsim n_{\rm heat,H_{2}}$ metal line cooling is ineffective because it cannot compete with H2 heating.

Equating the H2 and CR heating rates gives,

[TABLE]

where we used our H2 heating expression, Eq. (A15), assuming UV pumping heating dominates (but still in the $n<n_{\rm crit}$ limit), and used Eq. (13) for the effective gas-phase formation rate coefficient $R_{-}$. Next, equating the metal cooling density $n_{{\rm cool},Z}$ (Eq. 18) with $n_{\rm heat,H_{2}}$ above gives,

[TABLE]

In the numerical evaluation we adopted a metal cooling efficiency $\mathcal{L}_{z}=1.5\times 10^{-23}$ erg cm3 s*-1* at $n=10^{4}$ cm*-3*, $T=10^{3}$ K. For $Z^{\prime}\lesssim Z_{\rm heat,H_{2}}$ the gas never undergoes a second metal-line cooling phase, and there is no multiphase structure. This formula is in agreement with the numerical results shown in Fig. 4.

Eq. (27) is a strict requirement on the metallicity, expressing the metallicity below which a multiphase cannot possibly occur. In practice, as discussed above, much before the metallicity becomes that low, the multiphase region may already become very narrow, so that any thermal instability effects may become negligible.

5. Parameter Study

In this section, in Figs. 6 and 7, we present results of our thermal phase computations in the form of equilibrium temperature contour plots in the metallicity versus density ($Z^{\prime}-n$) plane, as well as temperature and density contour plots in the metallicity versus pressure ($Z^{\prime}-P$) plane. We investigate the effects of variations in $\zeta$ and $I_{\rm UV}$. We also discuss the time scales relevant for our assumptions of thermal and chemical equilibrium.

5.1. Contour plots

The top panels of Fig. 6 show the equilibrium gas temperature in the metallicity versus density ($Z^{\prime}-n$) plane, for $\zeta_{-16}=I_{\rm UV}=1$. The colors indicate the temperature, ranging from cold $\sim 10^{2}$ K (blue), to intermediate $\sim 10^{3}$ K (yellow), to warm $\sim 10^{4}$ K (red). In the upper-right panel we overplot contours of the cooling times. The cooling times range from $\sim 1-10$ Myrs in the WNM, depending on the metallicity, down to 0.1-0.4 Myrs in the CNM at low $Z^{\prime}$, and down to less then 0.01 Myrs at high $Z^{\prime}$. At higher metallicities the increase in the heating rate due to the onset of PE heating, is accompanied by appropriately higher cooling rates. For the CNM, the cooling rate naturally increases with increasing metallicity, as [CII] and [OI] dominates the cooling. In the WNM, where Ly$\alpha$ dominates, the cooling rate increases via a slight increase in the equilibrium temperature. Since at low $Z^{\prime}$ a large portion of the parameter space is WNM, or the intermediate $\approx 10^{3}$ K medium (at very low $Z^{\prime}$), with cooling times 1-10 Myrs, the gas may potentially be out of thermal equilibrium, depending on the magnitude of the local dynamical time (e.g., the sound crossing time, turbulent time, free-fall time, etc.). We further discuss the cooling and chemical time in the Appendix.

In the upper-left panel the various curves delineate the dominating heating and cooling processes that operate in different regions of parameter space. For heating, these are (a) photoelectric (PE) heating, above the dashed curve from the WNM into the CNM, (b) H2 heating, to the right of the dashed curve at high density in the CNM at moderately reduced metallicity, and again at high densities and very low metallicities in the intermediate temperature zone, and (c) cosmic-ray (CR) heating everywhere else. For cooling these are (a) Ly$\alpha$ to the left of the solid curve in the low density WNM, (b) metals to the right of the solid curve into the CNM, and (c) H2 cooling below the solid curve in the lower right portion of parameter space.

For any $Z^{\prime}$ the gas temperature decreases monotonically with increasing $n$, with sharp drops in $T$, from the WNM to CNM (red to blue), occurring near the Ly$\alpha$/metal or H2/metal cooling point boundaries. The WNM-to-CNM transition points depend on metallicity. The intermediate temperature zone appears at very low $Z^{\prime}$ and high density.

For $Z^{\prime}\gtrsim 0.2$ PE heating dominates over CR ionization in the WNM. Because PE heating and metal cooling both vary linearly with $Z^{\prime}$ (see Eq. 5), in this regime the Ly$\alpha$ to metal cooling point runs vertically (independent of $Z^{\prime}$) at an approximately constant density of $\approx 0.5$ cm*-3*. Below $Z^{\prime}\approx 0.2$ CR heating takes over, and the cooling point then moves diagonally, with $n_{{\rm cool},Z}\propto 1/Z^{\prime}$, consistent with Eq. (19) for a cooling point temperature $T_{0}=6\times 10^{3}$ K appropriate for the transition from Ly$\alpha$ to metal cooling. Below a metallicity of $8\times 10^{-3}$, H2 cooling (rather than metals) takes over from Ly$\alpha$ in the WNM. This is consistent with our analytic estimate (Eq. 24) for $Z^{\prime}_{\rm cool,H_{2}}$. The Ly$\alpha$ to H2 boundary line is vertical (independent of $Z^{\prime}$) and occurs well inside the WNM at a density $\approx 3$ cm*-3* consistent with our Eq. 21 for $n_{\rm cool,H_{2}}$. Below $Z^{\prime}_{\rm cool,H_{2}}\approx 8\times 10^{-3}$ the WNM-to-CNM boundary remains diagonal with $n_{{\rm cool},Z}\propto 1/Z^{\prime}$ as given by Eq. (19), but with a lower $T_{0}$ (i.e., $T_{0}=5000$ to $600$ K, for $\log Z^{\prime}=-2.3$ to $-4.7$), as appropriate for the transition from H2 to metal cooling.

For $Z^{\prime}<Z^{\prime}_{\rm heat,H_{2}}\approx 10^{-5}$ and $n>n_{\rm heat,H_{2}}\approx 10^{4}$ cm*-3*, H2 heating dominates over CR ionization, and therefore metal cooling becomes ineffective. The critical values $Z^{\prime}_{\rm heat,H_{2}}$ and $n_{\rm heat,H_{2}}$ occur at the intersection of the CR-to-H2 heating and H2-to-metal cooling boundary curves in Fig. 6, and are consistent with our analytic estimates given by Eqs. (25) and (27). In this portion of parameter space an intermediate equilibrium temperature of $\sim 10^{3}$ K is set by the balance of H2 heating and cooling.

The lower panels of Fig. 6 show the temperature (left) and density (right) in the pressure versus metallicity ($P-Z^{\prime}$) plane (again for $\zeta_{-16}=I_{\rm UV}=1$). The regions of pure CNM (cold/dense) and WNM (warm/diffuse) are indicated, as is the intermediate temperature H2 zone at very low metallicity and high density. The multiphased region where the CNM and WNM overlap occurs between the solid black curves, from $P_{\rm min}$ (to the left) to $P_{\rm max}$ (to the right).

The position of the multi-phased zone in Fig. 6 is determined primarily by (a) the switch from PE to CR heating, (b) the dependence of the metal cooling point on $Z^{\prime}$, and (c) the onset of H2 cooling and heating at low metallicity. For $Z^{\prime}$ down to $\sim 0.1$ PE heating dominates, and the multiphased zone widens and occurs at reduced values of $P_{\rm min}$ and $P_{\rm max}$. The width is maximal at $Z^{\prime}\approx 0.1$ where PE heating is least effective compared to metal cooling. For $Z^{\prime}\lesssim 0.1$ CR heating takes over (independent of metallicty) and $P_{\rm min}$ and $P_{\rm max}$ increase as the transition to metal cooling requires larger densities.

As the metallicity is lowered below $Z^{\prime}_{\rm cool,H_{2}}\approx 8\times 10^{-3}$ the multiphased zone become narrower due to H2 cooling in the WNM. At $Z^{\prime}\approx 10^{-4}$ there still exists a (small) temperature jump at the WNM to CNM boundary, but the pressure range of the multiphase zone is already negligible. At $Z^{\prime}<Z^{\prime}_{\rm heat,H_{2}}\approx 10^{-5}$ metal cooling becomes ineffective, the temperature then decreases smoothly with increasing pressure, from the $\sim 10^{4}$ K Ly$\alpha$ cooled WNM, to the $\sim 10^{3}$ K H2 cooled intermediate temperature zone.

5.2. Varying $\zeta$, $I_{\rm UV}$ and $D_{-}$

In Fig. 7 we explore the thermal structures for varying UV intensities and cosmic-ray ionization rates. The left column shows the cases $\zeta_{-16}=0.1$ and 10, assuming $I_{\rm UV}=1$. The middle column is for $I_{\rm UV}=0.1$ and 10, with $\zeta_{-16}=1$. The right column shows models in which $I_{\rm UV}$ and $\zeta_{-16}$ scale together by $\pm 1$dex relative to our fiducial case $\zeta_{-16}=I_{\rm UV}=1$. This last set reflects a physical scenario in which the massive stars that produce the energetic UV photons, are also the sources of the cosmic-rays originating in the remnants of core-collapse supernovae. For the $Z^{\prime}-n$ plots, the PE/CR/H2 and Ly$\alpha$/metals/H2 heating and cooling boundaries are shown as the solid and dashed curves, as in the upper left panel of Fig. 6. For the $Z^{\prime}-P$ plots, the various phases including the multiphased zones are indicated as in the lower panels of Fig. 6.

In the lower panels we also show, as dotted white curves, a comparison to the WMHT formula for $P_{\rm min}$ as a function of $I_{\rm UV}$, $Z^{\prime}$, $Z^{\prime}_{d}$, and $\zeta_{-16}$,

[TABLE]

(their Eq. 33) over its validity range (their Eqs. 35-38). We include the factor $f\equiv(1-\delta_{\rm C}Z_{d}^{\prime}/Z^{\prime})/(1-\delta_{\rm C})$ to correct666While WMHT use elemental abundances similar to ours for $Z^{\prime}=1$, they do not have a metallicity dependent depletion in their model, and $f$ corrects their formula for this effect. for the metallicity dependent depletion onto dust (i.e., our Eq. 4). For $Z_{d}^{\prime}$ we use our Eq. (5). The agreement is not perfect due to the different treatment of cosmic-ray and X-ray heatings, and PAH chemistry. However, overall, the general trends with $Z^{\prime}$, $\zeta$ and $I_{\rm UV}$ are recovered.

The left-panels of Fig. 7 show that as $\zeta$ increases, the WNM, CNM and the multiphase regions, are all shifted to higher densities and pressures. Within the CR-heating region (below the almost horizontal dotted line) higher densities are required for cooling to offset an increasing CR heating rate. Within the PE-heating region, as $\zeta$ increases, the electron abundance increases, the dust-grains are less positively charged, and the release of photoelectrons from the dust-grains into the gas is more efficient (i.e., see Fig. 1 in Wolfire et al. 1995). In the middle panels, where $I_{\rm UV}$ varies (with $\zeta$ kept constant), the shift is observed only at high metallicities ($Z^{\prime}\gtrsim 0.1$) where PE heating dominates, as expected.

The critical metallicity at which H2 cooling becomes important, $Z^{\prime}_{\rm cool,H_{2}}$, increases with $\zeta$ and decreases with $I_{\rm UV}$. When $\zeta$ and $I_{\rm UV}$ scale together, $Z^{\prime}_{\rm cool,H_{2}}$ remains unchanged. This is because the H2 abundance (a) increases with $\zeta$, as the gas phase formation is more efficient, and (b) decreases with $I_{\rm UV}$, as photodissociation is more efficient. The observed dependence of $Z^{\prime}_{\rm cool,H_{2}}$ on $\zeta$ and $I_{\rm UV}$ is in good agreement with Eq. (24). The location of the Ly$\alpha$-to-H2 transition (the vertical white contour) is in good agreement with Eq. (21). For example, looking at the left panels of Fig. 7 we see that as $\zeta_{-16}$ increases from 0.1 to 10 (with $I_{\rm UV}=1$ kept constant), the Ly$\alpha$-to-H2 transition shifts from $n\approx 2$ to 8 cm*-3*. On the other hand, when $I_{\rm UV}$ varies from 0.1 to 10 (and $\zeta_{-16}=1$ is kept constant, middle panels), the transition point shifts from $\approx 0.7$ to $\approx 20$ cm*-3*. Finally, when both $\zeta$ and $I_{\rm UV}$ vary together from 0.1 to 10, the transition point shifts by two orders of magnitude, from $\approx 0.4$ to 40 cm*-3*.

The critical metallicity, $Z_{\rm heat,H_{2}}$, at which H2 heating takes over, increases with $\zeta$, and is independent of $I_{\rm UV}$, as predicted by Eq. (27). As discussed in §4.3, this metallicity is set by the requirements that H2 heating equal CR heating, which in turn must equal metal cooling, all of which are independent of $I_{\rm UV}$. An increase in the CR intensity increases both the H2 formation rate and the CR heating rate, but the latter is affected more strongly, and thus $Z_{\rm heat,H_{2}}$ increases with $\zeta$.

As discussed in §2.2.3, the H- photodetachment rate at low metallicity is uncertain. We consider variations between $D_{-}=5.6\times 10^{-9}$ s*-1* for the extrapolated Draine field (our fiducial models), up to $D_{-}=2.7\times 10^{-7}$ s*-1* for a Draine+Mathis field. To investigate the effect of a varying $D_{-}$, we have also computed a set of models with the higher Draine+Mathis $D_{-}$ rate. The resulting temperature map is identical to that in Fig. 6, with the only difference that the H2 cooling line is shifted to a higher density, $n_{\rm cool,H_{2}}=15$ cm*-3*, i.e., an increase by a factor of 3.5. Since the $n_{\rm cool,Z}$ curve is unaffected by variations in $D_{-}$ (as it is set by metal cooling), the critical metallicity for H2 cooling is reduced, and $Z_{\rm cool,H_{2}}$ is reduced to $10^{-3}$. The shift to a higher cooling density is expected as the higher detachment rate reduces the H2 formation rate, requiring a higher gas density for the onset of H2 cooling. Unlike $n_{\rm cool,H_{2}}$ and $Z_{\rm cool,H_{2}}$, the heating points $n_{\rm heat,H_{2}}$ and $Z_{\rm heat,H_{2}}$ are unaffected by the value of $D_{-}$, as they occur at very high densities where H- photodetachment is subdominant.

Quantitatively, the factor 3.5 shift in $n_{\rm cool,H_{2}}$ may be understood as follows. For the higher $D_{-}$ rate, photodetachment dominates H- removal for densities $n<10^{2}I_{\rm UV}T_{3}^{0.39}$ cm*-3*. For such densities, $\eta\approx 0.1T_{3}^{-0.39}I_{\rm UV}^{-1}(n/10{\rm cm^{-3}})$ (see Eq. 16), and the H2 abundance is

[TABLE]

At $T_{3}=6$ and $n=10$ cm*-3* (i.e., near the Ly$\alpha$-to-H2 transition) this $x_{\rm H_{2}}$ is a factor of $\approx 20$ smaller than the corresponding value in Eq. (17), and has a steeper dependence on density. Plugging this back into Eq. (29), we see that $n$ will now enter with a power 5/2, thus the factor 20 lower $x_{\rm H_{2}}$, needs to be compensated by an increase in density by a factor of $\approx 20^{2/5}\approx 3.3$, in agreement with our numerical results.

Importantly, while $\zeta$, $I_{\rm UV}$ and $D_{-}$ affect the values of the critical densities and metallicities, the qualitative behaviour, the shift of the WNM-to-CNM transition to higher pressures with decreasing metallicity, the quenching of the multiphase zone due to H2 cooling at low $Z^{\prime}$, and the final disappearance of the multiphase phenomenon at very low $Z^{\prime}$, is robust.

6. Summary and Discussion

In this paper we have studied the thermal phase properties of atomic gas from high to vanishing metallicity, with special emphasis on understanding how multiphase behavior is altered when H2 heating and cooling processes become dominant. As discussed in §1, the critical role of H2 cooling at low metallicity has long been recognized in many studies. However, multiphased CNM/WNM behavior into the low metallicity and primordial regimes has received less attention.

For solar metallicity, we reproduce the known result that a CNM phase cooled by [CII] and [OI] line emission, ($T_{\rm CNM}\sim 100$ K) and a WNM phase cooled by Ly$\alpha$ emission ($T\sim 6000$ K) phases can coexist at thermal equilibrium at the same thermal pressure (Field et al. 1969; Wolfire et al. 1995; Liszt 2002, WMHT). The WNM-to-CNM transition and the multiphase region shifts to higher densities and pressures with increasing UV intensity, $I_{\rm UV}$, and/or increasing ionization rate, $\zeta$, as photoelectric (PE) heating becomes more efficient.

As the metallicity decreases, the metal cooling rate decreases. However, as long as $Z^{\prime}\gtrsim 0.1$, PE heating continues to dominate so that the heating rate also decreases. For sufficiently low metallicities, $Z_{d}^{\prime}$ is expected to scale superlinearly with $Z^{\prime}$ (in our model this occurs for $Z^{\prime}<0.2$), PE heating then falls faster than metal cooling, as $Z^{\prime}$ decreases, and the CNM is then colder and is less dense compared to solar metallicity models. The pressure range that allows a multiphase is then larger compared to solar metallicity models.

For metallicities $Z^{\prime}\lesssim 0.1$, the ionization of neutral gas by cosmic-rays (CR) becomes the dominant heating mechanism of the gas (the metallicity value depends on the CR ionization rate). In this limit as the metallicity decreases, the metal cooling rate decreases, while heating remains constant, shifting the WNM-to-CNM transition to higher densities and pressures. In this limit, and when H2 heating and cooling processes are ignored, we find that the Ly$\alpha$-to-metal transition density (and associated WNM-to-CNM transition) scales as $\zeta/Z^{\prime}$.

However, while previous analytic models of the HI phase structure often ignored H2 cooling, we find that H2 cooling dominates and modifies the phase structure at low metallicity (as we show H2 cooling and heating also play a role at solar metallicity gas but only at very high densities, $\sim 10^{6}$ cm*-3*, deep in the CNM). The critical metallicity below which H2 cooling becomes important is $Z^{\prime}_{\rm cool,H_{2}}=8.0\times 10^{-3}(\zeta_{-16}/I_{\rm UV})^{2/3}$. Below this metallicity, H2 cooling reduces the temperature of the WNM, to temperatures as low as 600 K (the exact value depending on metallicity). This lowers $P_{\rm max}$ and results in a narrower multiphase zone, both in pressure and density.

Finally, at extremely low metallicities $Z^{\prime}<Z^{\prime}_{\rm heat,H_{2}}=1.2\times 10^{-5}\zeta_{-16}^{2/3}$, the density at the WNM-to-CNM transition is already so high that H2 heating starts contributing. In this limit the thermal structure is determined only by HI (Ly$\alpha$) and H2 cooling. There is a single-phase solution at any pressure, as the temperature decreases smoothly with gas pressure, from the $\approx 10^{4}-6000$ K WNM, down to $T\approx 600$ K, as determined by the balance of H2 cooling and heating.

Gravitational collapse is favored in the cold-dense CNM, where the free-fall time is relatively short. Our result that the transition from warm to cold gas occurs at ever increasing pressures as the metallicity is reduced (over the range $10^{-5}\lesssim Z^{\prime}\lesssim 0.1$, for a constant $I_{\rm UV}$ and $\zeta$), suggests that galaxies with metal-poor ISM (i.e., galaxies in early evolutionary stages or low mass dwarf galaxies), should have maintained very high interstellar pressures in order to allow the formation of a CNM phase. For example, for our fiducial model of $I_{\rm UV}=\zeta_{-16}=1$, and for a metallicity $Z^{\prime}=10^{-3}$ a minimum pressure, $P/k\approx 10^{5}$ cm*-3* K is required. This requires abnormally large ISM column densities (for hydrostatic pressure equilibrium), deep gravitational potential wells of the dark-matter halos, or a highly pressurized hot ambient gas that may provide the needed pressure support.

Alternatively, a WNM-to-CNM transition and a multiphase ISM may still occur at normal ISM pressures at low metallicity if the heating rate is reduced. This requires a reduced $\zeta$ and $I_{\rm UV}$ values, suggesting lower star-formation (supernovae) rates in low metallicity galaxies. For example, for standard Galactic conditions ($Z^{\prime}=I_{\rm UV}=\zeta_{-16}=1$), the multiphase occurs at $P\approx 3000$ cm*-3* K (Fig. 6). For $Z^{\prime}=10^{-3}$ the multiphase may still occur at a similar pressure if $I_{\rm UV}$ and $\zeta$ are reduced by a factor of 10 (Fig. 7, rightmost column, 3rd row). When connected to a theory of a self-regulated star forming ISM (Ostriker et al., 2010), this would imply that the star-formation rate per unit gas column (i.e., the normalization of the Kennicutt Schmidt relation) should be lower at very low metallicities compared to solar (e.g., a factor of 10 lower for $Z^{\prime}=10^{-3}$ compared to $Z^{\prime}=1$). If neither of the above mechanisms are sufficient to maintain the high pressure needed, the ISM would remain warm, at $\approx 10^{4}$ K. In such a case, a different mode of star-formation would be expected, with only very massive clouds susceptible to gravitational collapse.

Previous studies have suggested that the atomic ISM of star-forming galaxies is naturally driven towards the multiphase pressure range, which provides a regulation mechanism for star-formation (Parravano, 1988, 1989; Ostriker et al., 2010). The idea is as follows. If the star-formation rate of a galaxy increases this leads to an increase in the UV intensity which in turn increase the PE heating rate. The $P-n$ curve is then pushed up and to the right (see Fig. 1, red-dashed curve), the gas is driven to the warm (WNM) phase, resulting in a decrease in the efficiency of star-formation, and vice versa (see Fig. 1 in Ostriker et al. 2010). In the limit where CR dominates over PE heating, the regulation is also expected to hold, since CRs are accelerated in supernovae remnants, and thus the 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 provides a prediction for the thermal pressure of the atomic ISM in galaxies, as a function of metallicity, CR ionization rate and UV intensity (i.e., see Figs. 6, 7). This pressure may be then related to the star-formation rate and galaxy gas and dark matter halo, through the requirement of hydrostatic equilibrium.

At extremely low metallicities, the multiphase structure disappears, and the multiphase-regulation mechanism cannot hold. This may imply a different (potentially a bursty) mode of star-formation for galaxies in this early evolutionary stage. In this limit the dense gas remains warm, at $\approx 600$ K, and is cooled by H2 cooling. This temperature is a factor of $\sim 10$ higher than the classic CNM at solar metallicity. The corresponding Jeans (1902) mass is a factor of $\approx 10^{3/2}\approx 30$ higher. Stellar structures forming in this low metallicity regime are expected to be more massive compared to their high metallicity counterparts.

Our finding that the cooling of the atomic ISM at low metallicity (below $\approx 0.01$) is often dominated by H2 line emission from the moderately warm phase, may open up a new observational window for the study of the low metallicity ISM, in local dwarfs and in galaxies of early evolutionary stages at high redshift. This is particularly timely with the upcoming launch of the James Webb Space Telescope (JWST) which will be able to observe IR H2 lines at extreme sensitivity.

We thank John Black, Emeric Bron, Franck Le Pettit, Evelyne Roueff, Chris McKee, and Eve Ostriker, for fruitful discussion and helpful suggestions. We thank the referee for constructive comments that improved our manuscript significantly. SB thanks the Center for Computational Astrophysics at the Flatiron Institute for hospitality and funding where some of this research was carried out with AS. This work was also supported by the German Science Foundation via DFG/DIP grant STE 1869/2-1 GE 625/17-1 at Tel Aviv University.

Appendix A A. Heating and Cooling Processes

A.1. Line-Emission Cooling

In general, gas cooling through line-emission occurs whenever a particle (ion/atom/molecule) is collisionaly excited to a higher energy state, and then relaxes back to a lower energy state by (spontaneously) emitting a photon. In this process, part of the kinetic energy 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-level system (e.g., the [CII] 158 $\mu$m transition), for which an analytic expression for the cooling rate may be obtained. Let $\Delta E$ be the transition energy, and $x_{i}$ and $x_{j}$ the abundances of the cooling species and the collisional partners, respectively. The cooling rate per unit volume (erg cm*-3* s*-1*) is

[TABLE]

where $n_{{\rm crit}}\equiv A/q_{i,\downarrow}$ is the critical density at which the collisional deexcitation and radiative decay rates are equal, $A$ is the Einstein coefficient for spontaneous decay, and $q_{ij,{\uparrow}}$ and $q_{ij,{\downarrow}}$ the collisional rate coefficients for the upward and downward transitions respectively. The latter are related through $(q_{ij,{\uparrow}}/q_{i,{\downarrow}})=(g_{u}/g_{d})\exp[-\Delta E/(k_{B}T)]$, where $g_{u}$, $g_{d}$, are the quantum degeneracies. The combination $\sum_{j}x_{j}q_{ij,{\uparrow}}(T)\Delta E\equiv\mathcal{L}_{i}(T)$ is the cooling coefficient. For densities $n\ll n_{{\rm crit}}$, the cooling rate $L\propto n^{2}$, as the density of the coolant and the rate of collisions (which populate the upper level) are both $\propto n$. At $n\gg n_{{\rm crit}}$ the energy levels are thermalized and $L\propto n$. In this limit,

[TABLE]

where $\Lambda_{{\rm LTE}}$ is cooling rate per atom/molecule (erg s*-1*) 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 [OI] and H2), with the appropriate choice of an effective $\Delta E$ and $n_{\rm crit}$.

We calculate the line emission cooling for Ly$\alpha$, [CII] 158 $\mu$m, [OI] 63 and 146 $\mu$m, and the ro-vibrational transitions of H2 and HD.

Ly$\alpha$ cooling: We compute the cooling via radiative decay from the 2p and 2s states of neutral hydrogen collisionaly excited by electrons. For the collisional de-excitation rate coefficient, we adopt $q_{\downarrow}=8.63\times 10^{-6}\Omega/(g_{u}\sqrt{T/{\rm K}})$ cm3 s*-1*, $g_{u}=2(2l_{u}+1)$ where $l_{u}=0$ and $1$ for 2s and 2p is the quantum angular momentum number, and $\Omega$ is the dimensionless collision strength (Gould & Thakur, 1970; Spitzer, 1978). We adopt a constant $\Omega=0.26$ for the 2s level, and $0.4$ for 2p. The resulting total (2s+2p) collisional de-excitation rate coefficient is $q_{\downarrow}=1.7\times 10^{-8}/\sqrt{T_{4}}$ cm3 s*-1*, 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<10^{4}$ K.

Metal fine-structure cooling: We calculate the cooling arising from the from the fine-structure transition of [CII] ${}^{2}P_{3/2}^{\rm o}-^{2}P_{1/2}^{\rm o}$ at $\lambda=157.7$ $\mu$m, and the fine structure transitions of [OI] ${}^{3}P_{i}-^{3}P_{j}$ where $[i,j]=([2,1],[2,0],[1,0])$, with the most important being the $[i,j]=[1,0]$ at $\lambda=63.2$ $\mu$m. We include collisional excitations with electron, hydrogen atoms and para-H2 and ortho-H2. 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 [CI] (although we find it to be a subdominant coolant), using the analytic fits from Draine (2011, see his Table 6), based on data from Roueff & Le Bourlot (1990); Abrahamsson et al. (2007); Schroder et al. (1991); Staemmler & Flower (1991).

H2 ro-vibrational cooling: For the H2 cooling function we use the polynomial fits provided by Glover & Abel (2008) with updates from Glover (2015, see their Appendix A1.2), which account for excitations of the H2 ro-vibrational ladder through collisions with H, He, para-H2, ortho-H2, H+, and e, assuming a constant ortho-to-para ratio o/p=3. We find that our results are not sensitive to the adopted o/p value. The cooling functions are based on collisional data from Ehrhardt et al. (1968); Crompton et al. (1969); Linder & Schmidt (1971); Draine et al. (1983); Gerlich (1990); Krstić & 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 have also calculated the H2 cooling efficiency assuming the updated H-H2 collision rates of (Lique, 2015). We find that adopting this alternative collisional dataset results in at most a negligible difference in the total H2 cooling function (which includes all the collisional partners).

HD ro-vibrational cooling: For HD cooling we use the cooling function of Lipovka et al. (2005) 777We spotted typos in the formulae in Lipovka et al. (2005) (confirmed by 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_{\rm HD}$ should be replaced with $W_{\rm HD}/n$., based on collisional data from Roueff & Flower (1999); Roueff & Zeippen (1999).

A.2. Heating and Cooling by Dust

Dust grains and PAHs, contribute both to gas heating, through the injection of energetic photoelectrons, and to cooling, though dust/PAH-assisted recombination. For the photoelectric (PE) heating we follow Bakes & Tielens (1994) and WMHT and adopt

[TABLE]

where the PE efficiency is

[TABLE]

and

[TABLE]

In Eq. A3 and A4, we assume that the PAH abundance scales with the dust-to-gas ratio, $Z_{d}^{\prime}$, and adopt $\gamma=2.2\times 10^{-24}$ erg s*-1* and $\phi_{\rm PAH}=0.5$ (see Eqs. 19-21 in WMHT 888The numerical values in Eq. A3 and A4 differ from those presented in WMHT because we used $I_{\rm UV}$ as the basic UV intensity unit, where WMHT use the Habing (1968) $G_{0}$ units. The conversion is, $1I_{\rm UV}=1.7G_{0}$.). For recombination cooling we adopt

[TABLE]

with $\gamma_{\rm rec}=4.65\times 10^{-30}$ erg cm3 s*-1*.

In our model, we do not solve for the PAH chemistry. At high metallicity, the PAHs may affect the electron fraction as they introduce an additional recombination channel for the H+ (see WMHT). At low metallicity, this effect becomes unimportant as the abundance of the PAHs falls rapidly with the vanishing DGR. Even at high metallicity, where the PAHs are importnat for recombination, we find that the net effect on the phase diagrams is mild, as shown in Fig. 7 where we compare to WMHT (see also 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

[TABLE]

where $T_{d}$ is the dust temperature, and $\gamma_{\rm g-d}=3.8\times 10^{-33}$ erg cm3 s*-1*, (Hollenbach & McKee, 1989; Glover & Jappsen, 2007). The dust is typically much colder than the gas, and is approximated by $T_{d}=16.4I_{\rm UV}^{1/6}$ K (Draine, 2011).

A.3. Cosmic-Ray Heating

Cosmic-ray or X-ray ionization liberate energetic electrons into the gas. These energetic electrons lose energy through coulomb interactions, and ionization, dissociation and excitations. The fraction of energy that is lost into coulomb interactions heat the gas. We adopt

[TABLE]

for the CR heating rate, where $\zeta_{\rm p}=\zeta/(1+\phi_{s})$ is the primary CR ionization rate,

[TABLE]

is the number of secondary ionizations per each primary electron, and

[TABLE]

is the characteristic deposition energy per ionization event (Dalgarno & McCray, 1972; Draine, 2011). The free electrons in the gas can also interact directly with the CR protons, giving an additional heating rate

[TABLE]

For CR proton energies of $\sim 50$ MeV, $E_{\rm cr}=287$ eV (Goldsmith et al., 1969). This direct heating channel dominates when the electron fraction is larger than $x_{\rm e}\sim 0.02$.

A.4. H2 Heating

When LW photons are absorbed by H2 molecules, they lead to population of the excited electronic (LW) states. $\sim 10$ % of the radiative decays lead to H2 photodissociation, whereas in the rest $90$ % of the cases, the molecule decays through ro-vibrational transitions back to the ground state. If the gas density is sufficiently high, such that collisions rather than spontaneous emission deexcite the H2, the energy of the excited level is transferred to the colliding particle (H, e, H2, etc.) leading to gas heating. Following Burton et al. (1990) and Röllig et al. (2006), the pumping heating rate may be approximated with

[TABLE]

where $9D_{0}$ is the H2 pumping rate (9 excitations occur per dissociation), and $E_{\rm pump}$ and $n_{\rm crit}$ are the effective energy and critical density of the pseudo two-level transition (see discussion above). The second equality holds for chemical equilibrium, for which $D_{0}I_{\rm UV}x_{\rm H_{2}}=Rnx_{\rm H}$ where $R$ is the total H2 formation rate coefficient. In this limit, the UV pumping heating rate is proportional to the H2 formation rate, $Rn^{2}$. At low densities, below $n_{\rm crit}$, another factor of $n/n_{\rm crit}$ enters, reflecting the fact that when $n<n_{\rm crit}$ spontaneous emission dominates deexcitation, not collisions. Following Röllig et al. (2006) we adopt $E_{\rm pump}=1.12$ eV, and $n_{\rm crit}=1.1\times 10^{5}/\sqrt{T_{3}}$ cm*-3* for the effective energy released per UV pumping event and for the critical density.

We also include heating through H2 photodissociation, and adopt

[TABLE]

where $E_{\rm pd}=0.4$ eV (Black & Dalgarno, 1977). Photodissociation heating dominates over pumping heating when $n/n_{\rm crit}<E_{\rm pd}/(9E_{\rm pump})$, which occurs at $n\lesssim 4.4\times 10^{3}/\sqrt{T_{3}}$ cm*-3*, for our adopted parameters.

When a new H2 molecule is formed, part of the binding energy (4.5 eV) is converted into heat, either by providing translational energy of the H2 molecules, or by exciting the H2 ro-vibrational levels which when followed by collisional de-excitation lead to gas heating. The latter is effective when the gas volume density is large, of order of or higher than the critical density of H2. The heating resulting from H2 may be written as

[TABLE]

where the summation is over the formation channels of H2, (i) formation on dust grains, and (ii) the gas-phase H- formation sequence, and $R_{i}$ are given by Equations (9) and (13), for the dust and gas formation routes, respectively. For each of the formation processes, the first term in brackets, $E_{{\rm form,1}}$, is the energy that goes into translational energy of the products (for dust-formation, this is the H2. For the gas phase H- formation route, most of the energy goes to the electron), and $E_{{\rm form,2}}$ is the energy that goes into internal molecule excitation, which can then heat the gas through collisional deexcitation. We adopt $E_{\rm form,1}=0.2$ eV and $E_{\rm form,2}=4.48$ eV for the dust catalysis (Hollenbach & McKee, 1979, see their Appendix VI(c)). For the H- formation we follow Čížek et al. (1998) who finds that on average $\approx 0.6$ eV is released as translational energy, and adopt $E_{\rm form,1}=0.6$ eV and $E_{\rm form,2}=3.13$ eV.

We may obtain a simple form for the total heating rate by all of the H2 heating processes discussed above. Assuming that always a single formation route dominates (either the dust catalysis or the gas phase H- route) and that the H2 system has reached chemical equilibrium, we get

[TABLE]

Here we defined $E_{\rm H_{2},1}=E_{\rm form,1}+E_{\rm pd}$ and $E_{\rm H_{2},2}=E_{\rm form,2}+9E_{\rm pump}$ and assumed that the effective critical densities are the same for H2 formation and UV pumping. In practice, since these critical densities are effective two-level system approximations, their values may differ for different processes. With our assumed parameters, $E_{\rm H_{2},1}=(0.6,1.0)$ eV and $E_{\rm H_{2},2}=(14.6,13.2)$ eV for H2 formation on dust and via H-, respectively. The second term in Eq. (A15) is dominated by H2 pumping heating, and becomes greater than the first term when $n\gtrsim(5.0,\ 8.4)\times 10^{3}/\sqrt{T_{3}}$ cm*-3*, for H2 formation on dust and via H-, respectively.

Appendix B B. Reduced Chemical Network

To improve computational speed, we consider a reduced chemical network consisting of the species: H, H-, H+, H2, 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 §2.2.4. In summary,

the H2 is formed via dust catalysis and by the H- gas-phase route, and is destroyed by photodissociation. The H- is formed via radiative association, and is destroyed by associative detachment, photodetachment, and mutual neutralization. For the electrons we include production by CR-ionization of hydrogen, and photoionization of carbon assuming all the carbon is in the form of C+. The electrons are removed by radiative recombinations the protons. We assume the oxygen and helium are in neutral 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 significantly the computation speed. The reduced network was used to produce Figs. 6 and 7. We have confirmed the accuracy of the reduced network by comparing the equilibrium $T(n)$, and $x_{\rm H_{2}}(n)$ curves for the reduced and full networks at the various metallicities considered in Fig. 4.

Appendix C C. Time-scales

Thermal and chemical equilibrium are good assumptions whenever the dynamical time of interest, $t_{\rm dyn}$, e.g., the turbulent crossing time, the free-fall time, the cloud lifetime, etc., is long compared to the cooling time, $t_{\rm cool}$, and the chemical time, $t_{\rm chem}$.

C.0.1 The cooling/heating time

The cooling time from an initial temperature $T_{i}$ to an equilibrium temperature $T$ is

[TABLE]

where we assumed isochoric cooling of predominantly atomic gas (an adiabatic index $\gamma=5/3$). If the cooling time is dominated by the cooling rate near the final equilibrium temperature, the cooling time simplifies to

[TABLE]

This is a good assumption when the cooling function increases superlinearly with $T$, as occurs when H2 or Ly$\alpha$ emission dominate the cooling, for $T\gtrsim 6\times 10^{3}$ K, and if $Z^{\prime}<Z^{\prime}_{\rm cool,H_{2}}$ down to $\sim 1000$ K (see discussion in §4, and Fig. 3 below), or at low temperatures, $T\lesssim 100$ K where metal cooling dominates but is a steep function of the temperature. On the other hand, When metals dominate the cooling at $T\gtrsim 100$ K, the cooling function depends weakly on $T$ (with a power-index $<1$), and the cooling time is then governed by the initial conditions. In this case the cooling time may still be approximated by Eq. (C2), but with $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 equilibrium temperature, where $G\approx L$. This latter form of $t_{\rm cool}$ is particularly useful when $G$ scales approximately linearly with the density while being weakly dependent on temperature. For example, for CR heating, $G=\zeta_{p}E_{\rm cr}n$ and we get

[TABLE]

In the numerical evaluation we used $E_{\rm cr}=15.2$ eV and $\phi_{s}=1.57$ (the same values as in Eq. 18). We see that for typical $\zeta_{-16}=1$, the cooling time in the WNM is long, with $t_{\rm cool}\sim 20$ Myrs, for $T=6000$ K. For the CNM, the cooling time is considerably shorter, with $t_{\rm 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$.

C.0.2 Chemical Times

Generally, the time for the H2 (and H) abundance to reach steady-state is

[TABLE]

where $D$ is the local photodissociation rate. In highly shielded cloud interiors, $D\ll 2Rn$, and the H2 time equal the H2 formation time, $1/(2Rn)$, which may become long. However, in our models we consider the state of the ambient atomic medium exposed to the free-space UV field, with $D=D_{0}$. In this regime the H2 time is the H2 (free-space) photodissociation time, $1/D_{0}=550\ I_{\rm UV}^{-1}\ {\rm yr}$, and is always very short.

Other relevant chemical timescales are the e, C+, and O, timescales. At low metallicity, the steady-state electron abundance is set by the ionization-recombination equilibrium, as given by Eq. (12). Under the assumption of $x_{\rm e}<1$, the time-scale for this process to achieve steady-state is

[TABLE]

where $n_{0}\equiv n/(10^{0}{\rm cm^{-3}})$. The electron time is already short ($\ll 1$ Myr) in the WNM, and becomes extremely short in the CNM, with $t_{\rm e}\approx 40T_{2}^{0.75}/n_{2}$ years. For the C+ time, consider a gas that recombines from a higher ionization state, ${\rm C^{+2}+e\rightarrow\rm C^{+}}$. The recombination time is $t_{\rm C^{+}}=1/(nx_{\rm e}\alpha_{\rm C^{+}})\approx 4(n_{0}\zeta_{-16})^{-1/2}$ Myr, where we used $\alpha_{\rm C^{+}}=6.5\times 10^{-12}$ cm s*-1* (Nussbaumer & Storey 1969, see Table 14.7 in Draine 2011), and used the equilibrium electron abundance (Eq. 12; we can assume equilibrium for $x_{\rm e}$ since $t_{\rm e}\ll t_{\rm C^{+}}$). The timescale for the chemical equilibrium of O may be longer than that of C+, as the recombination rate for O+ is longer than for C*+2*. 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.