LIGHTING THE DARK MOLECULAR GAS: \textrm{H}_{2} AS A DIRECT TRACER

# Lighting the Dark Molecular Gas: H2 as a Direct Tracer

Aditya Togi11affiliation: Ritter Astrophysical Research Center, University of Toledo, 2825 West Bancroft Street, M. S. 113, Toledo, OH 43606
and J. D. T. Smith11affiliation: Ritter Astrophysical Research Center, University of Toledo, 2825 West Bancroft Street, M. S. 113, Toledo, OH 43606
###### Abstract

Robust knowledge of molecular gas mass is critical for understanding star formation in galaxies. The  molecule does not emit efficiently in the cold interstellar medium, hence the molecular gas content of galaxies is typically inferred using indirect tracers. At low metallicity and in other extreme environments, these tracers can be subject to substantial biases. We present a new method of estimating total molecular gas mass in galaxies directly from pure mid-infrared rotational  emission. By assuming a power-law distribution of  rotational temperatures, we can accurately model  excitation and reliably obtain warm ( K)  gas masses by varying only the power law’s slope. With sensitivities typical of Spitzer/IRS, we are able to directly probe the  content via rotational emission down to  K, accounting for of the total molecular gas mass in a galaxy. By extrapolating the fitted power law temperature distributions to a calibrated single lower cutoff temperature, the model also recovers the total molecular content within a factor of 2.2 in a diverse sample of galaxies, and a subset of broken power law models performs similarly well. In ULIRGs, the fraction of warm  gas rises with dust temperature, with some dependency on . In a sample of five low metallicity galaxies ranging down to , the model yields molecular masses up to larger than implied by CO, in good agreement with other methods based on dust mass and star formation depletion timescale. This technique offers real promise for assessing molecular content in the early universe where CO and dust-based methods may fail.

ISM:molecules—galaxies:ISM—infrared:ISM

## 1. Introduction

By a factor of more than ,  is the most abundant molecule in the universe, found in diverse environments ranging from planet atmospheres to quasars hosts (Hanel et al., 1979, 1981; Walter et al., 2003; Genzel et al., 2015). It was the first neutral molecule to form in the Universe and hence dominated the cooling of pristine gas at early times (Lepp et al., 2002). Stars form principally from molecular clouds, and most physical prescriptions for stellar formation rate (SFR) are therefore directly linked to the surface density of  gas (Kennicutt, 1998; Bigiel et al., 2008). To understand the star formation processes and how it varies and has evolved over the star-forming history of the Universe, it is necessary to accurately measure the mass and distribution of .

Despite its high abundance,  can be difficult to study directly. It possess no permanent dipole moment, which makes it a weak rotational emitter. In addition, the upper energy level for the lowest permitted rotational (quadrupole) transition is K above ground (Dabrowski, 1984; van Dishoeck & Black, 1986). Hence, the  gas that comprises the bulk of the molecular interstellar medium (ISM) is commonly believed to be too cold to be visible. And yet, despite these observational disadvantages, the advent of the Infrared Space Observatory (ISO) and in particular Spitzer’s Infrared Spectrograph (IRS, Houck et al., 2004) revealed pure rotational  emission from a rich variety of extragalactic sources, including normal star forming galaxies (Roussel et al., 2007), ultraluminous and luminous infrared galaxies (U/LIRGs Lutz et al., 2003; Pereira-Santaella, 2010; Veilleux et al., 2009; Stierwalt et al., 2014), galaxy mergers (Appleton et al., 2006), radio-loud AGN (Ocaña Flaquer et al., 2010), UV-selected galaxies (O’Dowd et al., 2009), quasar hosts (Evans et al., 2001), cooling-flow cluster systems (Egami et al., 2006), sources with extreme shock-dominated energetics (Ogle et al., 2010), and even star-forming sources at redshifts (Ogle et al., 2012).

In the absence of direct measurements of , the lower rotational line transitions of CO (the next most abundant molecule) are generally used as a molecular gas tracer. To convert the measured integrated intensities of CO to  column density a conversion factor, , is needed, typically calibrated against virial mass estimates in presumed self-gravitating clouds (Solomon et al., 1987; Scoville et al., 1987; Strong & Mattox, 1996; Abdo et al., 2010). However recent evidence indicates that the value of varies substantially, both within galaxies and at different epochs (Genzel et al., 2012; Sandstrom et al., 2013).

In molecular clouds in our Galaxy, the observed -ray flux arising from the cosmic ray interactions with  can be used to recover  gas mass accurately (Bhat et al., 1985; Bloemen et al., 1984, 1986), but this method requires information on the cosmic ray distribution not available in other galaxies. Another method to assess  content is to use dust as a tracer, with a presumed or recovered dust-to-gas (DGR) ratio (e.g. Sandstrom et al., 2013; Rémy-Ruyer et al., 2014). Among other biases, both the CO and dust-based methods lose reliability at low metallicity. Due to declining dust opacity and less effective self-shielding than , the CO abundance drops rapidly with decreasing metallicity (Wolfire et al., 2010; Bolatto et al., 2013). Variations in the radiation field can also lead to selective CO destruction (Hudson, 1971; Bolatto et al., 2013). At the lowest metallicities, DGR itself begins to scale non-linearly with the metal content of the gas (Herrera-Camus et al., 2012; Fisher et al., 2014; Rémy-Ruyer et al., 2014). Another method of inferring molecular gas content is to couple measured star formation rates with an assumed constant  depletion timescale (), equivalent to a constant star formation efficiency (SFE) in galaxies (Schruba et al., 2012). None of these methods provide direct tracers of the  reservoirs, and the relevant indirect tracers all rely on physical assumptions (e.g. constant or measurable values of , DGR, or ) which are unlikely to be valid in all environments.

Here we challenge the long-stated assumption that  rotational emission is a poor tracer of the total molecular content in galaxies. Our understanding of the structure of the molecular material in galaxies is undergoing significant revision. While CO traces the coldest component of the molecular ISM, half or more (and sometimes much more) of the molecular gas in galaxies is in a warmer, CO-dark state, where  persists (Field et al., 1966; Tielens, 2005; Draine, 2011; Wolfire et al., 2010; Pineda et al., 2013; Velusamy & Langer, 2014). What this means is that the common wisdom that all molecular gas is found at temperatures  K is incorrect. With molecular material existing in quantity at excitation temperatures 50 - 100 K, it becomes possible to use the rotational emission spectrum of the Universe’s dominant molecule to directly assess a substantial portion of the total molecular content in galaxies.

We introduce here a simple continuous temperature model which enables direct use of the  rotational emission from galaxies to recover their total molecular gas content. While several assumptions are required to make such a model possible, the resulting biases are expected to be completely distinct from those of the indirect gas tracers. The paper is laid out as follows. We describe the archival sample in § 2. A rotational temperature distribution model is presented in § 3, and in § 4 we describe the method and its calibration procedure. Section 5 presents results, discussions, applications and future prospects of our model. We summarize our conclusions in § 6.

## 2. Sample

We have selected a large sample of galaxies spanning a wide range of different physical properties with reliable  rotational emission detections. Our sample includes 14 Low Ionization Nuclear Emission Regions (LINERs), 18 galaxies with nuclei powered by star formation, 6 Seyfert galaxies, 5 dwarf galaxies, 11 radio galaxies, 19 Ultra Luminous Infrared galaxies (ULIRGs), 9 Luminous Infrared galaxies (LIRGs), and 1 Quasi Stellar Object (QSO) host galaxy. The galaxies were chosen to have at least three well-detected  pure-rotational lines with Spitzer/IRS (Houck et al., 2004), including either S(0) or S(1), occurring at 28.22 and 17.04 , respectively. We also required an available measurements of CO line emission (either J = 1–0 or 2–1) for an alternative estimate of the total molecular gas mass. See Table 1 for a complete sample list with relevant physical parameters and the references from which each galaxy was drawn.

The flux ratios are calculated using 70 and 160 intensities obtained from the Herschel-PACS instrument (Poglitsch et al., 2010). We convolved PACS 70 maps to the lower resolution of PACS 160 to calculate flux ratios for mapped regions, for which we have the  line flux. For ULIRGs and other unresolved galaxies, represents the global flux ratio. The PACS70/PACS160 ratios are listed in Table 1 for each galaxy.

### 2.1. Mir H2 rotational line fluxes

All MIR  rotational line fluxes were obtained from Spitzer-IRS observations at low (R 60 – 120) and high (R 600) resolution between 5 – 38 m. The Spitzer Infrared Nearby Galaxies Survey (SINGS, Kennicutt et al., 2003) is a diverse sample of nearby galaxies spanning a wide range of properties. The SINGS sample consists of LINERs, Seyferts, dwarfs, and galaxies dominated by star formation in their nuclei. We adopted the  rotational line fluxes for four lowest rotational lines (S(0) to S(3)) in SINGS galaxies from Roussel et al. (2007). The targets were observed in spectral mapping mode and the details of the observing strategy is described in Kennicutt et al. (2003) and Smith et al. (2004).

rotational line fluxes for 10 3C radio galaxies are from Ogle et al. (2010). This sample includes from S(0) to S(7). For the radio galaxy 3C 236, the fluxes of  rotational lines are obtained from Guillard et al. (2012).

The  rotational line fluxes for QSO PG 1440+356 were taken from the Spitzer Quasar and ULIRG Evolution Study (QUEST) sample of Veilleux et al. (2009), and those for the 19 ULIRGs were obtained from Higdon et al. (2006). The fluxes of  rotational lines for the NGC 6240 and other LIRGs in the sample are from Armus et al. (2006) and Pereira-Santaella (2010), respectively.  rotational lines are also observed in the molecular outflow region of NGC 1266 (Alatalo et al., 2011). Table 2 lists the flux of MIR  rotational lines for all our selected galaxies.

### 2.2. Cold molecular gas mass from CO line intensities

To test and calibrate a model which measures total molecular mass from the rotational lines of , an estimate of the “true” cold  gas mass in a subset of galaxies is required. For this purpose, we utilize the cold  gas masses for the SINGS sample compiled by Roussel et al. (2007), which are derived from aperture-matched CO velocity integrated intensities. The observed line intensities of the 2.6 mm CO(1–0) line within the measured CO beam size was chosen to match the aperture of the Spitzer-IRS spectroscopic observations. Roussel et al. (2007) assumed a CO intensity to molecular gas conversion factor () of 5.0 M(K km s pc) (equivalent to cm(K km s)). They applied aperture corrections by projecting the IRS beam and CO beam on the 8 Spitzer-IRAC map. The underlying assumption is that the spatial distribution of the aromatic band in emission and CO(1–0) line emission are similar at large spatial scales.

For the radio galaxies, molecular gas masses were compiled by Ogle et al. (2010), and the molecular gas mass was estimated from the CO flux density corrected and scaled to the standard Galactic CO conversion factor. For the radio galaxy 3C 236 Ocaña, Flaquer et al. (2010) estimated an upper limit to the cold  gas mass. The cold  gas masses for the LIRGs in our sample were obtained from Sanders et al. (1991). They used = 6.5 M(K km s pc) (equivalent to cm(K km s)). We derived aperture corrections to the CO intensities by projecting IRS and CO beam on the 8 IRAC map, since the CO beam size (55 aperture) is much larger than the IRS mapped region (). Required aperture corrections are typically less than a factor of 2. For the LIRGs NGC 7591, NGC 7130, and NGC 3256, the value of molecular gas mass was obtained from Lavezzi & Dickey (1998); Curran et al. (2000); Sakamoto et al. (2006), respectively. The global molecular gas masses for ULIRGs are from Rigopoulou et al. (1996); Sanders et al. (1991); Mirabel et al. (1989); Evans et al. (2002). For the quasar QSO PG 1440+356 Evans et al. (2001) derived the molecular gas mass.

Different studies have utilized substantially different values to calculate molecular gas masses. In the Milky Way disk, a modern value of M(K km s pc) (equivalent to cm(K km s)) has been found to be consistent with a wide variety of constraints (Bolatto et al., 2013). In this work, we scale all molecular gas mass estimates to correspond to . Since we require only the  gas masses in our analysis, the resulting molecular mass values were further reduced by a factor of 1.36 to remove the mass contribution of Helium and other heavy elements. Recent dust-based studies of a subset of the SINGS sample have revealed variations in both between and within these galaxies Sandstrom et al. (2013).

Table 2 lists the CO-based  gas mass estimated for each galaxy in the sample, aperture matched to the region of Spitzer/IRS coverage, along with  masses calculated using the conversion factor for central regions derived by Sandstrom et al. (2013) in parentheses, where available.

## 3. Model

The primary methods of measuring  mass all rely on a set of assumptions — that the factor is known and constant, that the dust-to-gas ratio is fixed or tied directly to metallicity, or that the star formation depletion time is unchanged from environment to environment. All of these methods also rely on indirect tracers: the CO molecule, which is 10,000 less abundant than , dust grains, with their complex formation and destruction histories and varying illumination conditions, or newly formed stars, which are presumed to be associated with molecular material with a consistent conversion efficiency. All of these methods suffer biases based on the physical assumptions made.

We propose a direct means of assessing total molecular gas mass using the MIR quadrapole rotational line emission of . Interstellar or direct ultraviolet radiation fields, shocks, and other mechanical heating processes are all potential sources of  excitation. The method we propose requires only the assumption that  rotational temperatures are, when averaged over the diversity of emitting environments and processes in galaxies, widely distributed, and smoothly varying. This is directly analogous to modeling the radiation field intensity which heats dust grain populations in galaxies with a smooth distribution (e.g. Draine et al., 2007; Galliano et al., 2003). We model the  temperature distribution as a smooth, truncated power law with fixed cutoff temperatures. Indeed, the  molecule is radiatively cooled by a cooling function which can be well approximated as a power law with slope 3.8 (Hollenbach & McKee, 1979; Draine & McKee, 1993). A similar analysis by Burton (1987) derived a cooling function with a power law index, n = 4.66, for  molecules in the temperature range 10–2000 K.

The temperature equivalents of  rotational levels are high, starting at T = 510 K (see Table 3). The Boltzmann distribution of energy levels, however, leads to substantial excitation even at more modest peak temperatures. The lowest upper energy levels at J = 2–4 are well populated at excitation temperatures of 50 to 150 K, whereas the higher-J levels are excited by temperatures from a few hundred to several hundred K. Typically, a small number of discrete temperature components are fitted to low-J rotational line fluxes, recovering warm  temperatures in this range (Spitzer et al., 1973; Spitzer & Cochran, 1973; Spitzer et al., 1974; Spitzer & Zweibel, 1974; Savage et al., 1977; Valentijn & van der Werf, 1999; Rachford et al., 2002; Browning et al., 2003; Snow & McCall, 2006; Ingalls et al., 2011).

Other studies have also adopted continuous temperature distributions for . Zakamska (2010) in her study of  emission in ULIRGs used a power law model to calculate the expected  rotational line ratios, where is the mass of  gas with excitation temperatures between and . In her sample, power law indices are required to reproduce the observed range of  line ratios. Similarly Pereira-Santaella et al. (2014) adopted a power-law distribution in  to model  rotational emission in six local infrared bright Seyfert galaxies, with power law indexes ranging from 4–5. In the shocked environment of supernova remnant IC 443, Neufeld & Yuan (2008) found a power law temperature distribution for  column density with power law index varying over the range 3–6 with an average value of 4.5.

The present practice for measuring warm (K)  gas mass is to use two or three distinct temperature components to model the  rotational line fluxes. Yet even on the scales of individual molecular clouds,  is present at a wide range of temperatures, calling into question whether discrete temperature components are physical. When used purely for assessing warm gas mass (T 100K), the power law method provides a more robust, unique and reproducible measure than discrete temperature fits, which are sensitive to the arbitrary choice of starting line pair. A continuous distribution also makes possible extrapolation to a suitable lower temperature to recover the total cold molecular gas mass (see § 4).

We fit the flux of  rotational lines using a continuous power law temperature distribution by assuming

 dN=mT−ndT, (1)

where is the column density of  gas between excitation temperature and , is the power law index, and is a constant. Integrating this distribution to recover the total column density, the scaling co-efficient is found to be

 m=Ntot(n−1)T1−nℓ−T1−nu, (2)

where and are the lower and upper temperatures of the distribution, respectively, and N is the column density of in the observed line of sight and temperature is presumed equal to the rotational temperature. For a continuous distribution of molecules with respect to temperature, the column density of molecules at upper energy level responsible for the transition line is

 Nj=∫TuTℓgjZ(T)×e−EjkTmT−ndT, (3)

where is the degeneracy value for the corresponding energy level . The degeneracy value for even and odd values of corresponding to para and ortho  is

 gj=2j+1, (4)

and

 gj=3(2j+1), (5)

respectively. The factor of 3 for odd values of is for ortho-hydrogen, which have parallel proton and electron spins, forming a triplet state. is the partition function at temperature . The model assumes local thermodynamic equilibrium (LTE) between ortho and para and the abundance ratio of ortho to para at any temperature, , in LTE is given by Burton et al. (1992). The partition functions for para and ortho  are

 Z(T)=Zp(T)=∞∑j=even(2j+1)e−EjkT (6)

and

 Z(T)=Zo(T)=∞∑j=odd3(2j+1)e−EjkT, (7)

respectively.

Table 3 list the wavelength of each  rotational line transition, the upper energy level of the corresponding transition in temperature units, and their radiative rate coefficients (Huber, 1979; Black & Dalgarno, 1976).

Assuming  emission to be optically thin, the flux observed in a given transition is

 Fj=hνANj+2Ω4π, (8)

where is the solid angle of the observation. Substituting the value of from equation 3 into Eq. 8 and then using equation 2, we obtain the total column density

 Ntot=4πFjλ(T1−nℓ−T1−nu)AhcΩ(n−1)∫TuTℓgj+2Z(T)e−Ej+2kTT−ndT (9)

From Eq. 8, we conclude that the column density is proportional to the flux and the corresponding transition wavelength , and inversely proportional to the spontaneous emission probability ,

 Nj+2∝FjλA (10)

To perform the fit, we develop an excitation diagram from the ratio of column densities,

 Nj+2N3=FjλjA1F1λ1Aj, (11)

where is the column density of the upper energy level of the transition S(1). We select , corresponding to the S(1) transition, since it is the brightest and most frequently detected line. Using equation 3, we find this column density ratio from the continuous temperature model to be

 Nj+2N3=∫TuTℓgj+2Z(T)×e−Ej+2kTT−ndT∫TuTℓg3Z(T)×e−E1kTT−ndT. (12)

We determined the parameters of our model by comparing the observed and the modeled column density ratios from equation 11 and 12, by varying the lower temperature and power law index (the upper temperature, , is fixed, see § 4). Knowing these parameters, we substitute their values in equation 9, to obtain the total column density. The total number of molecules is calculated using

 ntot=NtotΩd2, (13)

where is the total number of hydrogen molecules and is the distance to the object. We then calculate the total gas mass,

 MH2,tot=ntot×mH2, (14)

where is the mass of a hydrogen molecule.

## 4. Method & Calibration

Applying the model developed in § 3 permits direct assessment of the warm  column. The model-derived lower temperature is typically degenerate, so extrapolating the model to lower temperatures to recover the full molecular content requires a calibration procedure using a trusted independent estimate of  content. Here we describe both applications.

### 4.1. Warm H2

The power law temperature distribution of molecules in Eq. 1 has three primary parameters: upper temperature , power law index , and lower temperature . An excitation diagram formed from the  line flux ratios is the distribution of normalized level populations and constrains our model parameters. Excitation diagrams relate the column density of the upper level () of a particular transition, normalized by its statistical weight, , as a function of its energy level . In the following sections we describe how the model can be used to estimate the warm  gas mass. In practice, except in a few cases, the lower temperature cutoff is not directly constrained (see § 4.2.1). The model fit results themselves are given in Table 4.

#### 4.1.1 Upper temperature, Tu

The bond dissociation energy for  is 4.5 eV, corresponding to  K, hence  is not typically bound above this temperature.  present in photo-dissociation regions (PDRs) in the outer layers of the molecular cloud can reach temperatures of a few 100 K (Hollenbach & Tielens, 1997, 1999). For the SINGS galaxies, Roussel et al. (2007) concluded that the mass of  greater than 100 K contributes 1–30% of the total  gas mass. In the MOlecular Hydrogen Emission Galaxies (MOHEGs) of Ogle et al. (2010), K contributes only 0.01 of the total . K contributes less than 1 in ULIRGs (Higdon et al., 2006).

In our model, when is allowed to vary above 1000 K, we found negligible impact on the recovered total gas mass or the quality of the fit to the excitation diagram. We therefore fixed the upper temperature of the power law model distribution at K. ’s ro-vibrational excitation spectrum (including the 1-0 S(1) line at ) would provide sensitivity to even hotter gas, but the mass contribution of this very hot gas in the context of our model is insubstantial.

#### 4.1.2 Power law index, n

Keeping a fixed value of  K, the other two model parameters, and , are varied using a Levenberg-Marquardt optimization (Markwardt, 2009) to match the observed line flux ratios (suitably converted to column density ratios using Eq. 11). Figure 1 is an example excitation diagram for galaxy NGC 5033. The model fit converges to  K, , with a fixed value of  K. Model predicted ratio values for the unobserved lines S(4)–S(7) are indicated, as are the two discrete temperature fits chosen by Roussel et al. (2007).

The power law index in our sample ranges from 3.79–6.39, with an average value . Figure 2 is a frequency distribution of power law index required to fit the MIR  rotational line fluxes in the SINGS galaxies. The power law index range derived in our model is comparable to the indices required to fit the  rotational line fluxes for ULIRGs (2.5 n 5.0) and Seyfert galaxies (3.4–4.9) (Zakamska, 2010; Pereira-Santaella et al., 2014).

Galaxies with a steep power law index have low warm gas mass fractions, since a larger quantity of  in these systems are at lower temperatures. Molecular gas heated by shocks and other turbulent energetic phenomena have higher temperatures and hence, higher warm gas mass fractions than gas in photo-dissociation regions (PDRs) (Appleton et al. in prep). The power law index therefore gives information on the relative importance of gas heating by shocks, photoelectric heating, UV pumping, etc.

### 4.2. Total H2

With a continuous power law model well reproducing the rotational emission lines from warm , it is natural to consider whether the total molecular gas reservoir could be probed by suitable extrapolation of our model to lower temperatures. Typically, the entire reservoir of  cannot be probed through rotational emission, since the model loses sensitivity at temperatures far below the first rotational energy state. Recovering the total molecular content from rotational emission therefore typically requires an additional free parameter — an extrapolated model lower temperature, — which must be calibrated against known molecular mass estimates. We first explore the models sensitivity to low temperatures, and then describe the calibration procedure we adopt.

#### 4.2.1 Model Sensitivity Temperature, Ts

The upper energy level of the lowest rotational transition of  is 510 K (see Table 3) — substantially higher than typical kinetic temperatures in molecular regions. Yet, even at temperatures well below this value, molecules in the high temperature tail of the energy distribution can yield detectable levels of rotational emission. Below some limiting temperature, however, increasing the number of cold  molecules will increase the implied molecular gas mass, but result in no measureable changes to the  rotational line fluxes.

We define the sensitivity temperature, , as that temperature below which the  reservoir is too cold for changes in the amount of gas to lead to measurable changes in the excitation diagram. To evaluate this in the context of our model, we calculated the difference between the model-derived and observed column density ratios, . Note that with the adopted normalization to the line, a unique excitation curve exists for any given pair , (). We varied these model parameters and evaluated the quality of the model fit using

 χ2=m∑i=1(Ri,mod−Ri,obsσRi,obs)2, (15)

where , and and are, respectively, the modeled and observed flux ratios for the transition, with uncertainty , and the summation is over all independent line flux ratios. We map the space in and to values, following Avni (1976), via (where is the minimum value). Figure 3 shows example contours for the galaxy NGC 5033. The value of is determined at the maximum value of lower cutoff temperature along the contour. For the galaxy NGC 5033,  K.

The near-horizontal orientation of the contours indicates little correlation between and . Since for most of the sample, the contours do not close as you go towards lower temperatures , any chosen lower cutoff of the power law distribution of temperatures below remains consistent with the data.111This indicates that the precise form of the temperature distribution at is not well constrained; see § 5.2.

Generally, differences between the modeled and observed flux ratios decrease as decreases, and do not significantly change below 80 K. However, in some LINER and Seyferts systems (e.g. NGC 2798, NGC 3627, and NGC 4579), the contours close, and has a defined minimum. Figure 4 illustrates this phenomenon in the galaxy NGC 3627, evaluating along the ridge line of best-fitting indices for each , and contrasting this galaxy with the more typical case of an indefinite minimum. Evident in the fits to NGC 3627 is a single best-fitting lower cutoff temperature of  K. In such cases, a continuation of the power-law distribution to temperatures below this best-fitting value is counter-indicated, as it degrades the fit. This indicates that in these cases, the bulk of  is directly detected via rotational emission in a warm gas component excess.

To identify galaxies with warm molecular gas excess, we imposed the constraint , where is the value at 20 K. Figures 5 and 6 show the distribution of sensitivity temperature and the corresponding average confidence contour plot for all SINGS galaxies, excluding in total 8 galaxies, identified as having a warm gas excess (NGC 2798, NGC 3627, and NGC 4579) and few other galaxies (e.g. NGC 1266, NGC 1291, NGC 1316, NGC 4125, and NGC 5195) with low signal-to-noise (S/N) ratio in their S(0) line. The excluded galaxies are all LINERs or Seyferts with low S(0)/S(1) ratios (in the range 0.04–0.2, vs. median 0.33 in the SINGS sample, Roussel et al. (2007)). They exhibit evidence for shocks, unusually high dust temperatures, and merger morphologies — all processes that could result in substantial quantities of warm molecular gas.

As seen in Fig. 6, the average value of is 81 K. This implies that, with the quality of  rotational spectroscopy provided by the Spitzer/IRS instrument,  gas down to rotational temperatures of 80 K can be reliably detected.

#### 4.2.2 Model extrapolated lower temperature, Tℓ

Given our model sensitivity to  columns down to rotational temperatures 80 K, we evaluate the potential for extrapolating the power law distribution to lower temperatures to attempt to recover the the total molecular gas mass.

To calibrate an extrapolated lower temperature cutoff – — we require a set of galaxies with well-established molecular gas masses obtained from . We constructed a training sample from SINGS galaxies (except NGC 4725, where we have an upper mass limit of molecular gas from CO intensity), omitting galaxies with a central warm excess (see § 4.2.1). ULIRGs and radio galaxies are also avoided due to ambiguity in their values.

Recent studies have also shown in low metallicity galaxies can be substantially higher than the Galactic value (e.g. Bolatto et al., 2013; Schruba et al., 2012; Leroy et al., 2011) and hence we restricted our sample to galaxies with an oxygen abundance value of , using the average of their characteristic metallicities derived from a theoretical calibration and an empirical calibration as recommended by Moustakas et al. (2010).

To perform the calibration, the single power law model fitted to the rotational excitation diagrams of each of the 34 galaxy training sample is extrapolated to the temperature where the model mass equals the cold molecular gas mass, estimated from CO emission (), using (column 10 of Table 2).

Figure 7 shows the distribution of extrapolated model lower temperature, , required to fit the molecular gas mass in the SINGS sample (with and without warm excess sources). The full training sample is described by an average extrapolation temperature of  K. Excluding galaxies with warm gas excess narrows the distribution, but does not significantly change the average value of . Table 4 lists the value of for each galaxy, with the value in parentheses calculated when the conversion factor derived by Sandstrom et al. (2013) through dust emission is assumed.

It is important to note that the numerical value of the lower temperature cutoff of  K relies on the assumption of a single power law across all temperatures. See § 5.2 for an alternative viewpoint adopting broken power laws, which can change the value of without impacting the estimated total molecular gas mass.

#### 4.2.3 Mass distribution function

Since is a strong negative power law, the number of molecules decreases rapidly with increasing temperature — most of the  gas mass resides at the lowest temperatures. We can evaluate the fraction of the total molecular mass our model is sensitive to by comparing and

Figure 8 shows an example mass-temperature distribution for NGC 5033 as a function of molecular temperature. The total molecular gas mass of 1.9 M is distributed with an index in the temperature range  K to  K. The fraction of cold  gas mass ( K for NGC 5033) is .

Using the average power-law index and sensitivity temperature from the SINGS training sample (see Fig. 6) we find . Adopting instead the 2 contour, the value of becomes 97 K and the ratio M is . Even when considering a broken power law temperature distribution for H, the mass fraction remains similar, and we discuss the effects in §5.5. This demonstrates that, despite the many shortcoming of  as a rotational emitter, we can trace a substantial fraction of the  in galaxies directly.

## 5. Results, Discussions, & Applications

### 5.1. What are the typical molecular gas temperatures in galaxies?

UV pumping sets the level populations of upper level states, J 3, in particular at low column densities of . As a result, these states may result in rotational temperatures well in excess of their kinetic temperatures. For our purposes, whatever combination of collisional (including shocks) and UV pumping determines the level populations and hence average the  excitation, we assume a single power law distribution of rotational temperature describes the ensemble of molecules (although see § 5.2 for an alternative analysis employing broken power-laws). Since the mass is dominated by gas with the lowest rotational temperatures and hence lowest excitations, it is instructive to compare the derived lower temperature cutoffs (§4.2.2) with various theoretical models and measurements of temperature in molecular gas.

The average lower extrapolation cutoff,  K, is comparatively higher than what is typically assumed for molecular regions (10–30 K). And yet, this high lower cutoff is essential in describing the portion of molecular material following a single power-law distribution of temperatures. For example, if we forced the extrapolation temperature to a lower value of 20 K while retaining the average power law index (§ 4.1.2), the estimated molecular gas mass would be on average higher than the mass measured using .

The ratio of warm diffuse to cold dense molecular gas in galaxies is a key step in understanding the  temperature distribution. Assuming an incident radiation field G = 10 for a cloud with mass 10–10 M, Wolfire et al. (2010) estimated that molecules in the outer diffuse region () obtain temperatures in the range 50–80 K, with stronger incident radiation fields leading to even higher temperatures. Due to the differences in self shielding from dissociating radiation, molecular clouds have an outer layer of varying thickness which contains molecular hydrogen, but little or no CO. This gas cannot be traced through CO transitions and is hence called dark molecular gas (Wolfire et al., 2010). This dark gas can account for a significant fraction (24–55%) of the total molecular gas mass in galaxies (Pineda et al., 2013; Wolfire et al., 2010; Smith et al., 2014; Planck Collaboration et al., 2011).

CO-dark and diffuse molecular gas is heated to higher temperatures than dense CO-emitting gas in molecular cloud cores. Indeed, individual molecular cloud simulations find average mass-weighted  temperatures of 45 K for Milky-Way average cloud masses and radiation (Glover & Clark (2012a), Glover, priv comm.). Galaxy-scale hydrodynamic simulations with a full molecular chemistry network (Smith et al., 2014) also show that  gas in the CO-emitting regions is markedly biased to the low temperature end of the full temperature distribution (Glover, priv comm.). Also, in diffuse and translucent molecular clouds in the Galaxy, Ingalls et al. (2011) demonstrated that additional sources of heating are required to explain the observed  and atomic cooling-line power. Our extrapolated power law model traces both warm and cold molecular gas mass in galaxies. The common assumption that the bulk of molecular gas is found at rotational temperatures of 10–30 K may be true only in the CO emitting core regions of the molecular cloud. The average mass-weighted molecular gas temperature can be much higher when the full range of emitting environments are considered. It is therefore perhaps not surprising to find typical molecular gas temperatures of  K in galaxies, similar to our model derived values.

It must be noted, however, that to calculate the total molecular gas mass, we have extrapolated using a single power law index. The possibility of a broken or non-power-law temperature distribution for the 85% of  gas at temperatures lower than our sensitivity temperature cannot be excluded (and indeed in warm gas excess sources, will be required to explain ongoing star formation). In the following section, we consider the effect of broken power law models on the dominant  temperatures.

### 5.2. Broken power law models

In the previous subsection we discussed the results of applying a single power law model for the entire temperature distribution of , and found that a single lower cutoff temperature of = 49 K permits recovering the entire molecular mass. Yet as seen in § 4.2.1, with the typical performance of the Spitzer/IRS instrument, we are sensitive only to the 15% of  gas with temperature  K. While a single power law temperature distribution, varying only the slope , provides very high quality model fits to the  excitation diagram, at temperatures below , in principle any distribution form would provide an equally valid description. Chemo-hydrodynamic models with  formation and excitation do result in a broad distribution of  temperatures (Glover & Smith, in prep.). But shocks elevate  to higher rotational temperatures (e.g. Neufeld & Yuan, 2008), and are not well modeled, so as yet there exists no a priori expectation for the detailed 10–3000 K temperature distribution of  molecules. Any presumed distribution, however, must 1) reproduce the  excitation diagrams as well as a single power law does, and 2) permit recovering the CO-predicted  mass within a reasonable range of temperatures. A broken power-law provides a general framework for exploring changes in the low-temperature form of the distribution, and has been invoked by Pereira-Santaella et al. (2014) to jointly model  and CO excitation. Here we consider the changes resulting from adopting a broken power law model, with a fixed break temperature, , and slope below the break, , studying the resulting lower cutoff temperatures, as well as the ability of these models to satisfy these two criteria.

Holding the upper temperature constant as before, we evaluated a grid of broken power law models with fixed pre-break slope (from 0–3) and break temperature (from 70–290 K). For each such model, we fitted the observed MIR  rotational line ratios in excitation diagrams in the same training sample of 34 galaxies considered in § 4.2.2. As before, we then extrapolate the broken power law to recover the total  mass, as indicated by CO intensity, evaluating the associated extrapolation cutoff temperature . We define as the average value of across the training sample, analogous to for a single power law (§ 4.2.2).

Figure 9 shows the variation of the training sample’s average extrapolated lower cutoff temperature as a function of break temperature for differing values of below-the-break slope. Those combinations of (,) for which the excitation diagrams are reproduced as well as a single power law model ( relative to the single power law model fit; see § 4.2.1) for at least 50% of the training sample are shown in solid lines. In this allowable range of and , the dispersion of the ratio of  mass recovered from the broken power-law model using a single extrapolation , relative to that implied by CO ranges only from 0.28–0.34 dex. At higher temperatures the excitation diagram fit quality degrades rapidly (dotted lines). As the break temperature approaches the average slope for a single power law model,  K, the broken power law becomes equivalent to our single power law truncated at this temperature, independent of slope . At slopes approaching the single-power law average (see Fig. 2), lower extrapolated temperatures diverge less and less from the value obtained using a single power law, since the shape of the temperature distribution is then relatively unchanged. Low slopes require very low extrapolation temperatures, since the cumulative mass of  rises less rapidly at lower temperatures compared to the single power law. In fact, models with a flat slope () fail to recover the CO-derived  mass within a reasonable lower  temperature floor of  K.

Figure 9 demonstrates that the single lower cutoff temperature to which an  model must be extrapolated depends sensitively on the detailed shape of the distribution at low temperature. For reasonable slopes, the dominant mass-weighted temperature can range from 10–50 K. All models which produce satisfactory fits of the excitation diagrams and have physically reasonable extrapolation temperatures have nearly the same predictive capability for the total  mass. While optical depth effects introduce some difficulties (e.g. for low-J CO), ideally a family of profiles for the full temperature distribution of molecular gas in galaxies could be developed, based on tracers with sensitivity to different parts of the temperature range matched to chemo-hydrodynamic models of molecular gas across all densities. Since the single power law model provides excellent fits to  rotational emission line ratios with the minimum of free parameters, we retain this simplified model for the remainder of this work. It should, however, be kept in mind that broken power law models with flatter slopes below 100–120 K can produce similar quality fits and yield equivalently reliable model-based  masses, as long as the low temperature slope does not vary substantially from galaxy to galaxy.

### 5.3. Estimating total Mh2

By fitting a continuous temperature distribution to the MIR  rotational lines and calibrating an extrapolating temperature, this model can be used to calculate the total molecular mass in galaxies directly from  rotational emission. While this method does require a reliable source of known molecular masses for calibration (introducing secondary dependence on, e.g. ), once calibrated it is independent of any indirect tracer like CO, DGR, or assumptions about star formation depletion timescales. The biases inherent in this method (arising, for example, from the assumption of a single and smooth power law temperature distribution), are therefore expected to be relatively distinct from those of the aforementioned estimators. In this section we test the model’s capability to estimate molecular gas mass in different types of galaxies.

Adopting a fixed model lower extrapolation temperature K the total  gas mass is calculated by extrapolating the fitted model. The total  gas mass derived by our model is shown in Figure 10. It compares very well with along with based estimates of gas mass. The scatter in our model is about 0.31 dex (factor of 2) for the SINGS sample and increases to 0.34 dex (factor of 2.2) for the complete sample, including U/LIRGs and radio galaxies. Some of this scatter no doubt arises from ignorance of the true values in these systems.

Galaxies with warm gas excess and a high sensitivity temperature (§ 4.2.1), require similar values of extrapolated as the training sample. The warm molecular gas traced by MIR rotational lines may be completely isolated from the cooler gas in these galaxies. Extrapolating to a lower temperature of 49 K yields molecular gas masses which agree with CO-derived values, assuming . A possibility of enhanced CO excitation due to high molecular gas temperature with a corresponding low in such warm galaxies cannot be ruled out (see also § 5.4).

Taken together, a tight correlation between -derived and CO-based mass estimates is found, spanning seven orders of magnitude in mass scale and across a wide range of galaxy types. The mass calculated with a continuous temperature distribution model, extrapolated to a single fixed cutoff temperature of 49 K can provide an independent measurement of total molecular gas mass in galaxies, good to within a factor of 2.2 (0.34 dex). A dependence on does arise indirectly through our CO-based mass estimates in the training sample (§ 4.2.2), but the recovered dispersion is comparable to uncertainties in the conversion factor itself, even among normal galaxies (Bolatto et al., 2013), as well as methods using DGR (e.g. Sandstrom et al., 2013).

### 5.4. Model derived molecular gas mass in ULIRGs, LIRGs and radio galaxies

Many local ULIRGs are recent ongoing galaxy mergers. In the merging process a large amount of gas in the spiral disk is driven to the central nuclear region, increasing the gas temperature. The increase in temperature and turbulence increases the CO linewidth, resulting in a high value of for a given molecular gas mass. therefore gives an overestimate of  gas masses in these galaxies (Downes et al., 1993; Bryant & Scoville, 1999). Moreover, can yield molecular gas masses greater than the observed dynamical mass (Solomon et al., 1997). To avoid this, a lower value of conversion factor is suggested for ULIRGs and other merger systems — = 0.8 M(K km s pc, 5.5 lower than the standard Galactic value (Downes & Solomon, 1998). However, by considering the high-J CO ladder, some studies have suggested that even for ULIRGs values are possible (Papadopoulos et al., 2012). Some  emission may lie outside photo-dissociation and star forming regions in ULIRGs (Zakamska, 2010). This  may reside in CO-dark gas, so that applying a low value in ULIRGs may yield an underestimate  mass.

In radio galaxies molecular gas can be predominantly heated by shocks through powerful jets. The molecular gas clouds may be affected by turbulence, and not gravitationally bound, defying the use of standard . Ogle et al. (2014) using DGR in NGC 4258, a low luminosity AGN (LLAGN) harboring a jet along the disk, derived gas mass of about 10 M, an order of magnitude lower compared to the standard method of using . The molecular gas mass could be overestimated when used in radio galaxies, which harbor long collimated powerful jets.

In applying a power law model to the sample of ULIRGs, LIRGs and radio galaxies, the nominal model extrapolation temperature  K is used to calculate the total  gas masses, listed in column 3 of Table 5. The in column 4 of Table 5 is the required extrapolated temperature to match the cold molecular gas mass measured from the CO line intensity.

In radio galaxies 3c 424, 3c 433, cen A, and 3c 236 the estimated  gas mass using the power law model after extrapolation to = 49 K, is higher when compared to the CO luminosity derived values using . However, when accounted for the intrinsic variation in and the  gas masses are in agreement with each other except in 3c 424, where the difference in masses is more than 10.

Using = 0.8 M(K km s pc, a factor of 5.5 lower than the standard Galactic value, we derive a modified extrapolation temperature , which is required to match the lowered gas mass. Since the model mass rises rapidly to lower temperatures, . The values are listed in column 7 of Table 5. On average for ULIRGs,  K.

Figure 11 shows the temperature distribution of (adopting ) and (adopting ) in ULIRGs. The lower temperature cutoff in ULIRGs and radio galaxies (except 3c424) is very similar to the normal galaxy training sample when is adopted, but much higher with the reduced molecular mass of .

It is of interest that when is used with the nominal calibration cutoff temperature , the ULIRG sample in Fig. 10 does not exhibit any particular bias. Either and are similar to their normal Galactic values in these systems, or they have reduced and a higher  temperature floor. This may indicate that the same physical processes that lead to reduced in highly active system, including increased ISM pressure and radiation density, globally elevate the gas temperature.

Since ULIRGs could indeed have uniformly elevated molecular gas temperatures, this degeneracy between decrease and increase leads to a systematic uncertainty in the total  gas mass identical in form to that obtained directly from mass estimates based on . A suggested prescription which side-steps this ambiguity, in applying this model to systems with non-Galactic is to calculate the total gas mass using the the nominal  K, and scale it by , for the preferred . As can be seen in Fig. 10, which adopts a uniform , the molecular mass in ULIRGs is well recovered by this procedure.

### 5.5. Effect of dust temperature on the warm H2 fraction

Since typical sensitivity temperatures are  K (see § 4.2.1), we can directly calculate the warm  mass above  K without extrapoloation, and compare it to the total mass as recovered by the calibrated extrapolation of § 4.2.2. Given an estimate for the power law index, , and lower cut off temperature, , of the power law distribution, we can calculate the fraction of molecular gas mass at temperatures above 100 K as

 M(>100K)Mtotal=M(>100K)M(H2,CO)=∫Tu100KT−ndT∫TuTℓT−ndT, (16)

where M is the total molecular gas mass estimated from the CO line intensity. Assuming and we find

 M(>100K)Mtotal≈(100KTℓ)1−n. (17)

Table 4 lists the calculated mass fraction for each galaxy. Columns 6 and 8 of Table 5 are calculated using and as generally assumed for ULIRGs, respectively. At lower , the warm gas mass fraction is higher.

Figure 12 shows the fraction of warm molecular gas mass, , as a function of far infrared (FIR) dust color temperature, . The warm gas fraction obtained using ranges from 2–30%, and exhibits little correlation among different galaxy types with dust color temperature. The ULIRGs and normal star forming galaxies show similar warm gas fractions, though ULIRGs have warmer dust color temperatures, and LINERs and Seyferts have somewhat higher warm gas mass fractions than normal star forming galaxies on average.

In contrast, using the available dust-derived central estimates of Sandstrom et al. (2013) for normal galaxies and a reduced value for ULIRGs and QSO’s, however, leads to a strong correlation, with warmer dust implying an increasing warm  fraction. The average warm gas fraction for ULIRGs and QSO is then 45%, significantly above that of normal galaxies, and on the same increasing trend with dust temperature. Depending on which prescription for total gas mass is correct, this could indicate a dependence of on temperature.

When considering a broken power law model, the fraction M(100 K)/M is, as expected, very similar to that of a single power law for 100 K. For 120 K irrespective of low-temperature slope , the model yields a poor fit to the excitation diagrams (see § 5.2). At these intermediate break temperatures, the warm mass fraction M(100 K)/M changes by at most 10–20% compared to the case of a single power law.

### 5.6. Molecular gas in low metallicity galaxies

As traced by their CO emission, many low metallicity dwarfs appear to have vanishingly low molecular gas content, but retain high star formation rates. That is, they are strong outliers on the Schmidt-Kennicutt relation (Galametz et al., 2009; Schruba et al., 2012). This discrepancy implies either that in dwarf galaxies star formation efficiencies are higher compared to normal spirals, or they host large molecular gas reservoirs than is traced by the CO emission (Schruba et al., 2012; Glover & Clark, 2012b). It is possible that a significant fraction of  exists outside the CO region, where the carbon is in C (ionized) or C (neutral) states. Since  can self-shield from UV photons in regions where CO is photodissociated (Wolfire et al., 2010), at low metallicity, not only is the CO abundance reduced, but as dust opacity is reduced and ionized regions become hotter and more porous, can severely underestimate the molecular gas mass. Considerable effort has been invested in detecting and interpreting CO emission at metallicity 50 times lower than the solar metallicity, , to assess the molecular gas content (Leroy et al., 2011; Schruba et al., 2012; Cormier et al., 2014; Rémy-Ruyer et al., 2014). Dust emission can be used to estimate the molecular gas mass in ISM, assuming a constant DGR however, it is essential to know the change in DGR with metallicity. At very low metallicities, , DGR appears to scale non-linearly with metallicity (Herrera-Camus et al., 2012; Rémy-Ruyer et al., 2014).

Our direct detection of  gas mass through  rotational lines is independent of any indirect tracers, which are affected by changing metallicity and local radiation effects. Applying our model in low metallicity galaxies should yield an estimate of molecular gas mass without the same inherent biases introduced by these dust and CO abundance variations. In this section we estimate the  gas masses through our power law model in a low metallicity galaxy sample selected to have detected  rotational emission, faint CO detection, and (where available) estimates of dust mass. We then compare these -based estimates to models and other methods which attempt to control for the biases introduced at low metallicity.

The low metallicity galaxies were selected on availability of MIR  rotational lines to have atleast three rotational lines including S(0) or S(1) line fluxes along with CO derived molecular mass estimates.

#### 5.6.1 Metallicity estimation

To study the variation of  gas mass from CO derived measurements over the metallicity range, it is essential to estimate the metallicity of galaxies. The metallicities were determined applying the direct method. CGCG 007-025 is the lowest metallicity galaxy in our dwarf sample with the value of = 7.77 (Izotov et al., 2007). Guseva et al. (2012) estimated the value of  in the two H II regions, Haro 11B and Haro 11C as 8.1 and 8.33, respectively hence, we adopt the average value for Haro 11. The metallicity value, , for NGC 6822 is 8.2 (Israel, 1997). No literature value exist for the oxygen gas phase abundance for the specific region of Hubble V in NGC 6822, mapped by the IRS-Spitzer. Peimbert et al. (2005) estimated = 8.420.06 for Hubble V, which is inconsistent with the previous value of 8.2.

For selected SINGS galaxies, the metallicity values in the circumnuclear regions, which are approximately the size of our  line flux extracted regions, are estimated by averaging the theoretical (KK04) and an empirical metallicity calibration (PTO5) as recommended by Moustakas et al. (2010).

#### 5.6.2 Cold molecular gas from CO line emission

Although the CO abundance drops super-lineraly with decreasing metallicity, it is detected in the low metallicity sample, and as the most common molecular tracer, can be compared directly to our based method. We adopted the literature values for CO(1–0) line intensities for CGCG 007-025 and N66, while for Haro 11, UM 311 and Hubble V region CO(3–2) line intensities were scaled using the relation = 0.60 (in temperature units) due to unavailability of CO(1–0) line intensities. Calculating , (area integrated luminosity) the molecular gas masses are estimated using and were further scaled using the 8 map, to account for the difference in the extracted IRS spectrum and the CO beam regions for each galaxy. The molecular gas masses are listed in Table LIGHTING THE DARK MOLECULAR GAS:  AS A DIRECT TRACER.

#### 5.6.3 Molecular gas from dust emission

An alternative method for estimating molecular gas mass makes use of dust emission together with assumption of dust opacity and grain size distribution to calculate a total mass, scaling dust mass to the total gas mass using a presumed or modeled dust-to-gas ratio, and removing the measured atomic mass from the region.

Leroy et al. (2007) estimated  gas surface density using dust emission from FIR map in N66 region of SMC. They derived to be about 27. For Haro 11, UM 311 and Hubble V using metallicity-DGR relation from Sandstrom et al. (2013) and with the known dust mass, we estimated the total gas mass and after subtracting the atomic gas content subsequently calculated the molecular gas mass. However, we find a negative value for the molecular gas for UM311, suggesting the ISM mass is mainly dominated by the atomic gas. The  gas masses calculated from the dust emission are given in Table LIGHTING THE DARK MOLECULAR GAS:  AS A DIRECT TRACER.

#### 5.6.4 Molecular gas mass using our model

The  line flux extraction was performed for the similar region in SMC-N66, where the CO emission was measured by Rubio et al. (1996). The cubes were prepared using CUBISM (Smith et al., 2007b) (Jameson, K. et al. in prep) and the  line fluxes were estimated using the PAHFIT, a MIR spectral decomposition tool (Smith et al., 2007a). The flux of  rotational lines for Haro 11 are from Cormier et al. (2014), while for UM 311, and CGCG 007-025 Hunt et al. (2010) measured the  rotational line fluxes. For Hubble V region, Roussel et al. (2007) derived the flux of  rotational lines.

Figure 13 is an excitation diagram for low metallcity galaxy Haro 11, with the power law model fit. Our model prediction for the unobserved S(0) line is included, and is consistent with the estimated upper limit. The  total gas mass for each low metallicity galaxy is measured using the power law model with lower temperature extrapolation to = 49 K. Table LIGHTING THE DARK MOLECULAR GAS:  AS A DIRECT TRACER lists the value of metallicity with their distance and the measured value of  rotational line fluxes with the power law index and the  gas mass derived using CO, dust, and our model.

Figure 14 compiles  gas masses derived using the various indirect tracers, together with the results from our -only model. All values are shown relative to the  mass inferred using CO luminosities with , and are plotted as a function of metallicity. For the SINGS sample, a variation of about 2–3 times the -derived  gas mass is found, likely a consequence of the intrinsic variation in at high metallicity  . At intermediate metallicity, the molecular gas content from CO line emission for the Hubble V region in NGC 6822 compares well with our model-derived gas mass, which suggests a similar Galactic conversion factor ( = ), consistent with the results of de Rijcke et al. (2006). At lower metallicity, however, our model disagrees strongly with naive CO-based estimates, yielding up to the molecular gas mass inferred from CO emission.

The molecular masses we recover at low metallicity are in good agreement with other measures which attempt to account for the impact of reduced metal abundance. These include dust-derived measurements, where available (when the strong metallicity dependence of the dust-to-gass ratio is accounted for). They also agree well with the prescription for the power-law like variations with metallicity recovered from inverting star formation densities among all non-starburst galaxies in the HERACLES sample (Schruba et al., 2012). The theoretical model of varying by Wolfire et al. (2010), assuming  column density of 10 cm in a molecular cloud, agrees well with observational results based on dust-mass at low metallicity (Leroy et al., 2011; Sandstrom et al., 2013). Adopting a solar metallicity value of 8.66 and heating rate/H atom = -0.3 (from the studies of Malhotra et al. (2001)), the Wolfire et al. (2010) model also agrees well with our  rotational emission modeled masses.

The power law model reliably recovers the total molecular gas masses at metallicity as low as 10 of the Milky Way, where CO and other indirect tracers suffer strong and non linear biases.

### 5.7. Future prospects

At metallicities 0.25 Z, the direct power-law method recovers total molecular gas content as reliably as other tracers that account for or avoid the impact of reduced gas-phase metal content. At even lower metallicities, the CO abundance plummets, with essentially all of the molecular gas in a CO-dark phase. In the first epoch of the star formation in the universe, the extremely low abundance of heavy elements leaves  as a principal coolant (Lepp et al., 2002).

The Mid Infrared Instrument (MIRI) onboard JWST is sensitive enough to detect the S(1), S(2), S(3) and higher rotational lines of  in luminous galaxies (U/LIRGs) till redshift 0.6, 1.3, 1.9, and higher, respectively. Assuming average = 10 (Bonato et al., 2015), and = 310 L (typical for LIRGs), for S(1) at z = 0.5 to have signal to noise S/N = 5 will require 30 minutes of integration time with JWST-MIRI. It will be possible to measure pure  rotational lines at high redshifts of z 6–7, almost reaching the reionization era of universe, with the SPace Infrared telescope for Cosmology and Astrophysics (SPICA) and the Cryogenic Aperture Large Infrared Submillimeter Telescope Observatory (CALISTO), planned for the 2020 decade (Roelfsema et al., 2012; Bradford et al., 2015).

The above mentioned future projects for the next decade provides an opportunity to observe  rotational lines at high redshifts. The power law model can be an useful tool in estimating molecular gas mass and study its variation and consequences at different redshifts.

## 6. Summary

We present a new power law temperature distribution model capable of reliably estimating the total molecular gas mass in galaxies purely from as few as three mid-infrared pure-rotational  emission lines. Our model is independent of the biases affecting indirect tracers like CO, dust emission, or star formation prescriptions. It can hence be used even in environments where reliability of those indirect tracers is poor, such as at low metallicity. We calibrate the model on a sample of local star-forming galaxies with well-quantified CO-based molecular gas masses, and apply the model to local ULIRGs and low metallicity systems. Our key results are:

1. A model based on a continuous power law distribution of rotational temperatures reproduces well the  excitation from a wide range of galaxy types using only a single parameter (the power law slope ). This simple model can directly recover the warm  mass ( K) more reproducibly than arbitrary discrete temperature fits.

2. The power law index obtained for all SINGS galaxies ranges from 3.79 - 6.4, with an average value of 4.84 0.61.

3. With typical Spitzer detection sensitivities, the model can directly recover the  gas mass down to a limiting sensitivity temperature of  K (when the S(0) line is available), accounting for of the total molecular reservoir.

4. By calibrating the model using a subset of the SINGS sample with well determined CO-based molecular masses, we find that extrapolating the  temperature distribution to a single cutoff temperature of  K enables recovery of the total  gas mass within a factor of 2.2 (0.34 dex).

5. Evaluating a family of broken power law models, we find that the average lower cutoff temperature could range from 5–50 K, while recovering total  gas masses as well as a single power law model. However, low model slopes below the break () would require unphysically low cutoff temperatures, and for all power law indices , break temperatures 120 K result in unacceptable fits to the  MIR rotational lines.

6. When is used, the fraction of warm molecular gas mass (M 100 K) in this training sample ranges from 0.02 to 0.30. If a reduced is adopted for the ULIRGs, this fraction increases with increasing dust color temperature .

7. In ULIRGs, the total molecular gas mass obtained by extrapolating the model to  K is consistent with the molecular gas mass derived using . Alternatively, if a reduced is adopted, the model extrapolation temperature required rises to 80 K. Either the warm molecular gas fraction and lower temperature cutoff in ULIRGs is higher than in normal star forming galaxies, or is closer to the Galactic value than has been presumed.

8. At low metallicity (), where indirect tracers of molecular mass suffer increasingly strong and non-linear biases and the mass of the molecular reservoirs exceed their CO-derived estimates by or more, the direct power law model recovers the total molecular gas masses reliably, as assessed by other methods which account for metallicity.

9. With upcoming facilities including JWST, SPICA, and CALISTO, detection of the relevant  rotational lines in galaxies at intermediate to high redshifts becomes possible, opening a new window on the fueling history of star formation in the Universe.

Acknowledgements: This paper is dedicated to James R. Houck of Cornell University (1940-2015), the PI of the Spitzer IRS. This work is based on observations made with the Spitzer Space Telescope. We have made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology under NASA contract. Support for this research was provided by NASA through contract 1287374 issued by JPL/Caltech under contract 1407. We thank the anonymous referee for useful comments and suggestions. We are grateful to K. Jameson and A. Bolatto for sharing the N66-SMC region spectrum. We are also indebted to L. Armus, N. Scoville, B. Draine, A. Bolatto, S. R. Federman, K. Sandstrom, D. Neufeld, F. Walter, T. Diaz-Santos, E. Pellegrini for discussions which improved this work, and to S. Glover for providing model access and for helpful discussions of simulated  temperature in galaxies. JDS gratefully acknowledges visiting support from the Alexander von Humboldt Foundation and the Max Planck Institute für Astronomie, and support from a Cottrell Scholar Award from the Research Corporation for Science Advancement. AT acknowledges support for this work through the IPAC Visiting Graduate Fellowship program.â

## References

• Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 710, 133
• Alatalo et al. (2011) Alatalo, K., Blitz, L., Young, L. M., et al. 2011, ApJ, 735, 88
• Appleton et al. (2006) Appleton, P. N., Xu, K. C., Reach, W., et al. 2006, ApJ, 639, L51
• Armus et al. (2006) Armus, L., Bernard-Salas, J., Spoon, H. W. W., et al. 2006, ApJ, 640, 204
• Avni (1976) Avni, Y. 1976, ApJ, 210, 642
• Bhat et al. (1985) Bhat, C. L., Issa, M. R., Houston, B. P., Mayer, C. J., & Wolfendale, A. W. 1985, Nature, 314, 511
• Black & Dalgarno (1976) Black, J. H., & Dalgarno, A. 1976, ApJ, 203, 132
• Bloemen et al. (1984) Bloemen, J. B. G. M., Caraveo, P. A., Hermsen, W., et al. 1984, A&A, 139, 37
• Bloemen et al. (1986) Bloemen, J. B. G. M., Strong, A. W., Mayer-Hasselwander, H. A., et al. 1986, A&A, 154, 25
• Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846
• Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
• Bonato et al. (2015) Bonato, M., Negrello, M., Cai, Z.-Y., et al. 2015, MNRAS, 452, 356
• Bradford et al. (2015) Bradford, C. M., Goldsmith, P. F., Bolatto, A., et al. 2015, arXiv:1505.05551
• Browning et al. (2003) Browning, M. K., Tumlinson, J., & Shull, J. M. 2003, ApJ, 582, 810
• Bryant & Scoville (1999) Bryant, P. M., & Scoville, N. Z. 1999, AJ, 117, 2632
• Burton et al. (1992) Burton, M. G., Hollenbach, D. J., & Tielens, A.G.G.M. 1992, ApJ, 399, 563
• Burton (1987) Burton, M. G. 1987, Ph.D. Thesis
• Cormier et al. (2014) Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2014, A&A, 564, A121
• Curran et al. (2000) Curran, S. J., Aalto, S., & Booth, R. S. 2000, A&AS, 141, 193
• Dabrowski (1984) Dabrowski, I. 1984, Canadian Journal of Physics, 62, 1639
• de Rijcke et al. (2006) de Rijcke, S., Buyle, P., Cannon, J., et al. 2006, A&A, 454, L111
• van Dishoeck & Black (1986) van Dishoeck, E. F., & Black, J. H. 1986, Interstellar Processes: Abstracts of Contributed Papers, 149
• Downes et al. (1993) Downes, D., Solomon, P. M., & Radford, S. J. E. 1993, ApJ, 414, L13
• Downes & Solomon (1998) Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615
• Draine & McKee (1993) Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373
• Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium by Bruce T. Draine. Princeton University Press, 2011. ISBN: 978-0-691-12214-4,
• Draine et al. (2007) Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866
• Dufour (1984) Dufour, R. J. 1984, Structure and Evolution of the Magellanic Clouds, 108, 353
• Dufour et al. (1982) Dufour, R. J., Shields, G. A., & Talbot, R. J., Jr. 1982, ApJ, 252, 461
• Egami et al. (2006) Egami, E.; Rieke, G.H., Fadda,D. & Hines, D. C. 2006, ApJ, 652, L21
• Evans et al. (1999) Evans, A. S., Sanders, D. B., Surace, J. A., & Mazzarella, J. M. 1999, ApJ, 511, 730
• Evans et al. (2001) Evans, A. S., Frayer, D. T., Surace, J. A., & Sanders, D. B. 2001, AJ, 121, 1893
• Evans et al. (2002) Evans, A. S., Mazzarella, J. M., Surace, J. A., & Sanders, D. B. 2002, ApJ, 580, 749
• Evans et al. (2005) Evans, A. S., Mazzarella, J. M., Surace, J. A., et al. 2005, ApJS, 159, 197
• Field et al. (1966) Field, G. B., Somerville, W. B., & Dressler, K. 1966, ARA&A, 4, 207
• Fisher et al. (2014) Fisher, D. B., Bolatto, A. D., Herrera-Camus, R., et al. 2014, Nature, 505, 186
• Galametz et al. (2009) Galametz, M., Madden, S., Galliano, F., et al. 2009, A&A, 508, 645
• Galliano et al. (2003) Galliano, F., Madden, S. C., Jones, A. P., et al. 2003, A&A, 407, 159
• Genzel et al. (2012) Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69
• Genzel et al. (2015) Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20
• Glover & Clark (2012a) Glover, S. C. O., & Clark, P. C. 2012a, MNRAS, 421, 9
• Glover & Clark (2012b) Glover, S. C. O., & Clark, P. C. 2012b, MNRAS, 426, 377
• Guillard et al. (2012) Guillard, P., Ogle, P. M., Emonts, B. H. C., et al. 2012, ApJ, 747, 95
• Guseva et al. (2012) Guseva, N. G., Izotov, Y. I., Fricke, K. J., & Henkel, C. 2012, A&A, 541, A115
• Hanel et al. (1979) Hanel, R., Conrath, B., Flasar, M., et al. 1979, Science, 206, 952
• Hanel et al. (1981) Hanel, R., Conrath, B., Flasar, F. M., et al. 1981, Science, 212, 192
• Harnett et al. (1991) Harnett, J. I., Wielebinski, R., Bajaja, E., Reuter, H.-P., & Hummel, E. 1991, Proceedings of the Astronomical Society of Australia, 9, 258
• Herrera-Camus et al. (2012) Herrera-Camus, R., Fisher, D. B., Bolatto, A. D., et al. 2012, ApJ, 752, 112
• Higdon et al. (2006) Higdon, S. J. U., Armus, L., Higdon, J. L., Soifer, B. T., & Spoon, H. W. W. 2006, ApJ, 648, 323
• Houck et al. (2004) Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18
• Hollenbach & Tielens (1999) Hollenbach, D. J., & Tielens, A. G. G. M. 1999, Reviews of Modern Physics, 71, 173
• Hollenbach & Tielens (1997) Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 179
• Hollenbach & McKee (1979) Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555
• Huber (1979) Huber, K. P. & Herzberg, G. 1979, Constants of Diatomic Molecules (New York: Van Nostrand)
• Hudson (1971) Hudson, R. D. 1971, Reviews of Geophysics and Space Physics, 9, 305
• Hunt et al. (2010) Hunt, L. K., Thuan, T. X., Izotov, Y. I., & Sauvage, M. 2010, ApJ, 712, 164
• Hunt et al. (2015) Hunt, L. K., García-Burillo, S., Casasola, V., et al. 2015, A&A, 583, A114
• Ingalls et al. (2011) Ingalls, J. G., Bania, T. M., Boulanger, F., et al. 2011, ApJ, 743, 174
• Israel et al. (1996) Israel, F. P., Bontekoe, T. R., & Kester, D. J. M. 1996, A&A, 308, 723
• Israel (1997) Israel, F. P. 1997, A&A, 328, 471
• Israel et al. (1990) Israel, F. P., van Dishoeck, E. F., Baas, F., et al. 1990, A&A, 227, 342
• Israel et al. (2003) Israel, F. P., Baas, F., Rudy, R. J., Skillman, E. D., & Woodward, C. E. 2003, A&A, 397, 87
• Izotov & Thuan (1998) Izotov, Y. I., & Thuan, T. X. 1998, ApJ, 500, 188
• Izotov et al. (2007) Izotov, Y. I., Thuan, T. X., & Stasińska, G. 2007, ApJ, 662, 15
• Kennicutt (1998) Kennicutt, R. C., Jr. 1998, ApJ, 498, 541
• Kennicutt et al. (2003) Kennicutt, R. C., Jr., Armus, L., Bendo, G., et al. 2003, PASP, 115, 928
• Lavezzi & Dickey (1998) Lavezzi, T. E., & Dickey, J. M. 1998, AJ, 115, 405
• Lepp et al. (2002) Lepp, S., Stancil, P. C., & Dalgarno, A. 2002, Journal of Physics B Atomic Molecular Physics, 35, 57
• Leroy et al. (2007) Leroy, A., Bolatto, A., Stanimirovic, S., et al. 2007, ApJ, 658, 1027
• Leroy et al. (2011) Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12
• Lutz et al. (2003) Lutz, D., Sturm, E., Genzel, R., et al. 2003, A&A, 409, 867
• Malhotra et al. (2001) Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766
• Markwardt (2009) Markwardt, C. B. 2009, Astronomical Data Analysis Software and Systems XVIII, 411, 251
• Mirabel et al. (1989) Mirabel, I. F.; Sanders, D. B. & Kazes, I. 1989, ApJ, 340, L9
• Moustakas et al. (2010) Moustakas, J., Kennicutt, R. C., Jr., Tremonti, C. A., et al. 2010, ApJS, 190, 233
• Nesvadba et al. (2010) Nesvadba, N. P. H., Boulanger, F., Salomé, P., et al. 2010, A&A, 521, A65
• Neufeld & Yuan (2008) Neufeld, D. A., & Yuan, Y. 2008, ApJ, 678, 974
• Ocaña Flaquer et al. (2010) Ocaña Flaquer, B., Leon, S., Combes, F., & Lim, J. 2010, A&A, 518, A9
• O’Dowd et al. (2009) O’Dowd, M. J., Schiminovich, D., Johnson, B. D., et al. 2009, ApJ, 705, 885
• Ogle et al. (2010) Ogle, P., Boulanger, F., Guillard, P., et al. 2010, ApJ, 724, 1193
• Ogle et al. (2012) Ogle, P., Davies, J. E., Appleton, P. N., et al. 2012, ApJ, 751, 13
• Ogle et al. (2014) Ogle, P. M., Lanz, L., & Appleton, P. N. 2014, ApJ, 788, L33
• Okuda et al. (2005) Okuda, T., Kohno, K., Iguchi, S., & Nakanishi, K. 2005, ApJ, 620, 673
• Pereira-Santaella (2010) Pereira-Santaella, M., Alonso-Herrero, A., Rieke, G. H., et al. 2010, ApJS, 188, 447
• Pereira-Santaella et al. (2014) Pereira-Santaella, M., Spinoglio, L., van der Werf, P. P., & Piqueras López, J. 2014, A&A, 566, A49
• Papadopoulos et al. (2012) Padelis P. Papadopoulos, Paul van der Werf; Xilouris, E.; Isaak, K. G. & Yu Gao 2012, ApJ, 751, 10
• Peimbert et al. (2005) Peimbert, A., Peimbert, M., & Ruiz, M. T. 2005, ApJ, 634, 1056
• Pineda et al. (2013) Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103
• Planck Collaboration et al. (2011) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, AA19
• Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2
• Rachford et al. (2002) Rachford, B. L., Snow, T. P., Tumlinson, J., et al. 2002, ApJ, 577, 221
• Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31
• Rigopoulou et al. (1996) Rigopoulou, D.; Lawrence, A.; White, G. J., Rowan Robinson, M., & Church, S. E. 1996, A&A, 305, 747
• Roelfsema et al. (2012) Roelfsema, P., Giard, M., Najarro, F., et al. 2012, Proc. SPIE, 8442, 84420R
• Roussel et al. (2007) Roussel, H., Helou, G., Hollenbach, D. J., et al. 2007, ApJ, 669, 959
• Rubio et al. (1996) Rubio, M., Lequeux, J., Boulanger, F., et al. 1996, A&AS, 118, 263
• Sakamoto et al. (2006) Sakamoto, K., Ho, P. T. P., & Peck, A. B. 2006, ApJ, 644, 862
• Salomé & Combes (2003) Salomé, P., & Combes, F. 2003, A&A, 412, 657
• Sanders et al. (1991) Sanders, D. B., Scoville, N. Z., & Soifer, B. T. 1991, ApJ, 370, 158
• Sandstrom et al. (2013) Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5
• Saripalli & Mack (2007) Saripalli, L., & Mack, K.-H. 2007, MNRAS, 376, 1385
• Savage et al. (1977) Savage, B. D., Bohlin, R. C., Drake, J. F., & Budich, W. 1977, ApJ, 216, 291
• Schruba et al. (2012) Schruba, A., Leroy, A. K., Walter, F., et al. 2012, AJ, 143, 138
• Scoville et al. (1987) Scoville, N. Z., Yun, M. S., Sanders, D. B., Clemens, D. P., & Waller, W. H. 1987, ApJS, 63, 821
• Smith et al. (2004) Smith, J. D. T., Dale, D. A., Armus, L., et al. 2004, ApJS, 154, 199
• Smith et al. (2007a) Smith, J. D. T., Armus, L., Dale, D. A., et al. 2007a, PASP, 119, 1133
• Smith et al. (2007b) Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007b, ApJ, 656, 770
• Smith et al. (2014) Smith, R. J., Glover, S. C. O., Clark, P. C., Klessen, R. S., & Springel, V. 2014, MNRAS, 441, 1628
• Snow & McCall (2006) Snow, T. P., & McCall, B. J. 2006, ARA&A, 44, 367
• Solomon et al. (1987) Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730
• Solomon et al. (1997) Solomon, P. M.; Downes, D.; Radford, S. J. E. & Barrett, J. W. 1997, 478, 144
• Spitzer & Zweibel (1974) Spitzer, L., Jr., & Zweibel, E. G. 1974, ApJ, 191, L127
• Spitzer et al. (1974) Spitzer, L., Jr., Cochran, W. D., & Hirshfeld, A. 1974, ApJS, 28, 373
• Spitzer & Cochran (1973) Spitzer, L., Jr., & Cochran, W. D. 1973, ApJ, 186, L23
• Spitzer et al. (1973) Spitzer, L., Drake, J. F., Jenkins, E. B., et al. 1973, ApJ, 181, L116
• Stierwalt et al. (2014) Stierwalt, S., Armus, L., Charmandaris, V., et al. 2014, ApJ, 790, 124
• Strong & Mattox (1996) Strong, A. W., & Mattox, J. R. 1996, A&A, 308, L21
• Tielens (2005) Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium, by A. G. G. M. Tielens, pp. . ISBN 0521826349. Cambridge, UK: Cambridge University Press, 2005
• Valentijn & van der Werf (1999) Valentijn, E. A., & van der Werf, P. P. 1999, ApJ, 522, L29
• Veilleux et al. (2009) Veilleux, S., Rupke, D. S. N., Kim, D.-C., et al. 2009, ApJS, 182, 628
• Velusamy & Langer (2014) Velusamy, T., & Langer, W. D. 2014, A&A, 572, A45
• Walter et al. (2003) Walter, F., Bertoldi, F., Carilli, C., et al. 2003, Nature, 424, 406
• Wolfire et al. (2010) Wolfire, M. G., David Hollenbach., & Christopher, F. McKee. 2010, ApJ, 716, 1191
• Young et al. (2011) Young, L. M., Bureau, M., Davis, T. A., et al. 2011, MNRAS, 414, 940
• Zakamska (2010) Zakamska, N. L. 2010, Nature, 465, 60