© ESO, 2014

1. Introduction

Molecular gas is an important phase of the interstellar medium (ISM). This phase contains a significant fraction of the total mass, and stars form in it. But the study of molecular gas presents some complications. First, the lower energy levels of H₂, the main component of the ISM phase, have energies >500 K, thus in cold molecular gas (T< 100 K) most of the H₂ is in the fundamental state and no H₂ emission lines are produced. And second, only the near infrared (IR) ro-vibrational H₂ transitions, with E_up> 6000 K, are observable from ground telescopes, so only very high-temperature (T> 1500 K) molecular gas can be detected.

To overcome the first caveat, other abundant molecules with observable transitions in the millimeter range (like CO, HCN, etc.) are used as tracers of molecular gas. In particular, the lowest rotational transitions of CO, the second most abundant molecule, are commonly used to study the molecular gas content of galaxies. However, these low-J CO transitions mainly originate in the coldest molecular gas. Thus, ground observations are limited to the study of either the warmest or the coldest molecular gas.

Just recently, thanks to IR and sub-millimeter space observatories like the Infrared Space Observatory (ISO; Kessler et al. 1996), the Spitzer Space Telescope (Werner et al. 2004), and the Herschel Space Observatory (Pilbratt et al. 2010), the rotational H₂ transitions as well as the intermediate-J CO transitions became accessible for a large number of local galaxies (e.g., Rigopoulou et al. 2002; Roussel et al. 2007; van der Werf et al. 2010; Pereira-Santaella et al. 2013). Therefore, now for the first time, it is possible to obtain a complete snapshot of molecular gas emission and study its physical properties (temperature, density, column density, etc.) and the excitation mechanisms (ultraviolet (UV) radiation, shocks, and X-ray and cosmic rays).

In this work, we present new data obtained by the Fourier transform spectrometer (FTS) module of the Spectral and Photometric Imaging Receiver (SPIRE) instrument on-board Herschel (Griffin et al. 2010; Naylor et al. 2010; Swinyard et al. 2010) for six local active luminous IR galaxies. These SPIRE/FTS data cover the 210–670  μm (450–1440 GHz) spectral range, so the mid-J CO lines (J = 4–3 to J = 13 − 12) are observed. We completed the CO spectral line energy distributions (SLEDs) with ground-based observations of the three lowest J CO transitions.

In addition, we complemented the CO SLEDs with the H₂ SLEDs obtained from near- and mid-IR observations of these galaxies. We used the available mid-IR spectroscopy obtained by the Spitzer IR spectrograph (IRS; Houck et al. 2004) to measure the lowest rotational H₂ transitions (e.g., Wu et al. 2009; Pereira-Santaella et al. 2010), and near-IR integral field spectroscopy obtained by the Spectrograph for INtegral Field Observations in the Near-Infrared (SINFONI; Eisenhauer et al. 2003) on the Very Large Telescope (VLT) for the ro-vibrational H₂ transitions. For the first time, we have performed a radiation transfer analysis of the whole set of molecular lines together (i.e., CO rotational and H₂ rotational and ro-vibrational) in local IR bright galaxies. In total, the compiled CO and H₂ SLEDs contain information for 26 transitions with upper level energies between 5 and 15 000 K, thus the emission from most of the molecular gas is included.

The paper is organized as follows. In Sect. 2, we present the sample and the data reduction. Sections 3 and 4 describe the radiative transfer models used, and the fitting of the SLEDs. The cold-to-warm molecular gas ratio and the heating mechanisms are discussed in Sects. 5 and 6, respectively. We summarize the main results in Sect. 7.

2. Observations and data reduction

2.1. Sample

Our sample contains six local (50–180 Mpc) IR bright Seyfert galaxies observed by SPIRE/FTS through two programs (PIs: P. van der Werf and L. Spinoglio). Their log L_IR/L_⊙ ranges between 11.4 and 12.2 (see Table 1). All of them host an active nucleus (AGN), although the AGN dominates the energy output of the galaxy only for IRASF 05189–2524. For the remaining galaxies, the AGN contributes between 2 and 35% of the total luminosity. Three of the galaxies (NGC 34, IRASF 05189–2524, and UGC 05101) are advanced major mergers. The other three galaxies (NGC 5135, NGC 7130, and NGC 7469) are spirals, although NGC 7130 and NGC 7469 have peculiar morphologies suggesting recent minor interactions (Genzel et al. 1995; Bellocchi et al. 2012).

2.2. Herschel SPIRE/FTS spectroscopy

We obtained SPIRE/FTS high-resolution (1.45 GHz) spectroscopic observations of six nearby Seyfert luminous IR galaxies. Integration times varied between 5 and 34 ks corresponding to a 3σ line detection limit ~0.5–1 × 10^-15 erg cm^-2 s^-1.

The SPIRE/FTS consists of two bolometer arrays, the spectrometer short wavelength (SSW; 925–1545 GHz) and the spectrometer long wavelength (SLW; 446–948 GHz). They sparsely cover a field of view (FoV) of ~2′ with the central bolometer centered at the nuclei of the galaxies. The full width at half maximum (FWHM) of the SSW beam is 18′′ almost constant with frequency. For the SLW bolometers the beam FWHM varies between ~30 and ~42′′ with a complicated dependence on frequency (Swinyard et al. 2010). The SSW beams size correspond to 5–15 kpc at the distance of our galaxies, therefore these galaxies are almost point like sources for both the SSW and SLW bolometers (see also Fig. 1 of Pereira-Santaella et al. 2013).

The data was reduced as described by Pereira-Santaella et al. (2013), but using the more recent Herschel interactive pipeline environment software (HIPE) version 11 (Ott 2010). In brief, the pipeline first creates the interferograms from the bolometer timelines. After the interferogram phase errors are corrected and the baselines removed a Fourier transform is applied to obtain the spectra. These are dominated by the thermal emission of the telescope that is later removed. Residual background emission is subtracted by averaging the off-source bolometers. Finally, the point-source calibration is applied to the spectra extracted from the central bolometers (SLWC3 and SSWD4). The final spectra are plotted in Fig. 1.

To measure the line fluxes, we fitted a sinc function to the line profiles¹. The local continuum was estimated using a linear fit around the line frequency. When two lines were close in frequency we fitted them simultaneously. The fluxes of the ¹²CO, HF, H₂O, and fine structure atomic transitions are listed in Table 2. Transitions of CH⁺, H₂O⁺, and OH⁺ detected in some galaxies are reported in Table 3.

The SPIRE/FTS data of two of these galaxies, UGC 05101 and NGC 7130, were already presented by Pereira-Santaella et al. (2013). To take advantage of the newest pipeline and calibration they were reprocessed with the rest of the sample. For these two galaxies, the line fluxes in Tables 2 and 3 agree with those in Pereira-Santaella et al. (2013) within the 1σ uncertainties.

Fig. 1 Observed SPIRE/FTS spectra of our sample. The black and gray lines are the SSW and SLW spectra, respectively. Notice the overlap between the two spectra in the 960 and 990 GHz spectral range. The dashed red lines mark the position of the detected lines. The error bars indicate the median 1σ uncertainty of each spectrum.

Fig. 1 continued.

2.3. Spitzer/IRS spectroscopy

Mid-IR Spitzer/IRS spectroscopic observations of these galaxies were available on the Spitzer data archive. Most of the data is already published in several papers (e.g., Armus et al. 2004; Wu et al. 2009; Pereira-Santaella et al. 2010; Tommasin et al. 2010; Esquej et al. 2012), but the fluxes (and upper limits) of the H₂ transitions are not always reported. For this reason, we retrieved the data from the archive and reprocessed it in a uniform way.

These observations include low- and high-resolution spectroscopy (R ~ 60–120 and 600, respectively) covering the spectral range 5.5–14  μm (short-low resolution modules) and 10–36  μm(short- and long-high-resolution modules).

The observations obtained in the staring mode were processed using the standard pipeline version S18.18 included in the Spitzer IRS Custom Extraction software (SPICE). We assumed the point source flux calibration for these galaxies. To process the spectral mapping observations, we used CUBISM (Smith et al. 2007a). This software combines the individual slit observations to create the spectral data cubes. From the cubes, we extracted the spectra using a 15′′ × 15′′ square aperture and then we applied an aperture correction (see Alonso-Herrero et al. 2012). The spectra are shown in Fig. 2.

We measured the fluxes of the H₂ rotational transitions S(0), S(1), and S(2) in the high-resolution spectra² fitting a Gaussian profile and a linear function for the local continuum (Table 4). The H₂ S(3) and S(5) transitions were measured in low-resolution spectra. Because of the large number of spectral features in the mid-IR range, it is not trivial to determine the continuum for the low-resolution spectra. Therefore, we used PAHFIT (Smith et al. 2007b) to fit the complete spectrum and estimate the continuum and line fluxes. It is not possible to obtain a reliable measurement of the H₂ S(4) transition at 8.03  μm because it is blended with the 7.7  μm PAH broad feature. Similarly, the H₂ S(5) transition at 6.91  μm is blended with the [Ar iii]6.99  μm line, which in our galaxies is always stronger than the H₂ line. Therefore, the H₂ S(5) fluxes should be considered with caution.

The strength of the silicate feature at 9.7  μm (S_Si) was measured in the low-resolution spectra as described by Pereira-Santaella et al. (2010). It is defined as S_Si = ln(f_obs/f_cont), thus negative S_Si implies that the feature is seen in absorption (Table 5). We used the strength of the silicate feature to estimate extinction in our galaxies (see Sect. 2.6).

Fig. 2 Rest-frame Spitzer/IRS spectra. The short-low (5.5–14  μm), short-high (10–19  μm), and long-high (19–36  μm) spectra are plotted for each galaxy. For clarity, the spectra are multiplied by the following factors 1, 2.5, 14, 100, 220, and 2000, from bottom to top. The wavelengths of the H2 rotational transitions are indicated.

2.4. SINFONI integral field spectroscopy

We made use of K-band (1.95–2.45  μm) seeing-limited and adaptive optics assisted near-IR SINFONI observations for several objects of the sample. The observations where carried out during different periods, i.e., 60A (NGC 7469), 77B (NGC 5135 and NGC 7130), 80B (IRASF 05189–2524), and 82B (NGC 34), with a spectral resolution of R ~ 4000, using different plate scales that yield different FoVs. For further details of the NGC 7469 data, see Hicks et al. (2009), and for NGC 5135 and NGC 7130, please refer to Piqueras López et al. (2012).

We made use of archive data for IRASF 05189–2524 and NGC 34, observed on December 2007 and October 2008, respectively. The observations were made using the 0.̋05 × 0.̋10 pixel^-1 setup, which yields a FoV of ~3′′ × 3′′, and split into individual exposures of 900 s following a jittering O-S-O pattern for sky and on-source frames. The total on-source integration times for each galaxy were 5400 s for IRASF 05189–2524 and 1800 s for NGC 34. To perform the flux calibration and to correct for the instrument response and atmospheric absorption, we used three spectrophotometric standard stars, Hip032193 and Hip052202 for IRASF 05189–2524, and Hip001115 for NGC 34, which were observed with their respective sky frames.

The reduction and calibration of the raw data were performed following the same procedure outlined in Piqueras López et al. (2012), i.e., we used the standard ESO pipeline ESOREX (version 2.0.5) to perform the standard corrections of dark subtraction, flat fielding, detector linearity, geometrical distortion, and wavelength calibration. After these corrections were applied to each frame, the sky emission was subtracted from each individual on-source exposure. We constructed individual cubes from each sky-corrected frame that were then flux-calibrated separately, and combined into a single final cube taking the relative shifts in the jittering pattern into account. The absolute flux calibration to translate the counts into physical units is based on the K-band magnitudes from the 2MASS catalog (Skrutskie et al. 2006).

To obtain the integrated K-band spectra, first we fitted the brightest H₂ line (1–0 S(1) at 2.12  μm) in every spaxel to produce the intensity and velocity map of the H₂ emission. Then we combined the spectra of all the spaxels with signal to noise ratio >3 in the 1–0 S(1) line, correcting for the relative velocity of each spaxel. The integrated spectra are shown in Fig. 3. Finally, we measured the line fluxes and upper limits of the H₂ transitions in the integrated spectra by fitting a Gaussian profile with fixed FWHM and position that were first determined from the 1–0 S(1) fit. Fluxes and upper limits are presented in Table 4.

2.5. Aperture matching

We combined data from instruments with different FoVs and angular resolutions. The SPIRE/FTS beam size varies between 20 and 40′′, which correspond to 5–35 kpc at the distance of our galaxies, so they appear as point-like sources at this resolution. By contrast, the Spitzer/IRS slit width varies between 3.6 and 11.1′′ (1–10 kpc), depending on the IRS module. For the most distant galaxies in our sample, IRASF 05189-2524 and UGC 05101, these apertures encompass their mid-IR emission. Likewise, most of the mid-IR emission of NGC 34 and NGC 7469 is produced in the compact nuclear region (few kpc, Esquej et al. 2012; Díaz-Santos et al. 2010), thus no aperture matching corrections are needed. For the more nearby NGC 5135 and NGC 7130, Spitzer/IRS spectral mapping observations are available and it is possible to extract their mid-IR spectra from a 15 × 15′′ aperture (5 × 5 kpc) comparable to the SPIRE/FTS beam size and large enough to include most of their bright nuclear mid-IR emission.

The SINFONI FoV for the IRASF 05189-2524 observations is 3 × 3′′ (2.5 × 2.5 kpc). This source is unresolved in the SINFONI data (FWHM ~ 0.3′′), and therefore we expect that most of the flux is included in our measurements. For NGC 7469, the FoV of SINFONI data is 0.8 × 0.8′′ (0.2 × 0.2 kpc), so only the nuclear emission is included. NGC 7469 has a kpc scale star-forming ring, thus to account for the missing near-IR H₂ emission we added the emission measured in the star-forming ring by Díaz-Santos et al. (2007) to the nuclear fluxes (Table 4). This correction increases the fluxes by a factor of ~2.

Finally, the SINFONI FoVs are 8 × 8′′  and 8 × 16′′, respectively, for NGC 5135 and NGC 7130. This corresponds to 2.5–5 kpc which is similar to the source sizes (Arribas et al. 2012), so no correction is applied to the measured integrated SINFONI fluxes.

2.6. Extinction correction

Near- and mid-IR emission in luminous and ultra-luminous IR galaxies (LIRGs and ULIRGs) is affected by extinction due to the large amount of dust present in these galaxies. To estimate the extinction, we measured the S_Si (see Sect. 2.3), S_Si and A_K are related according to the extinction law. In particular, using the Chiar & Tielens (2006) local ISM extinction law we obtained that A_K/S_Si = −2.20. The S_Si values in our galaxies range from −1.50 to −0.13 (Table 5), which correspond to A_K between 0.3 and 3.3 mag.

Regions producing H₂ and ~10  μm continuum emissions (where the S_Si is measured) may have different extinction levels. We cannot use the near-IR lines to test if we can use the estimated A_K values to correct for the obscuration effects in the H₂ transitions because the differential extinction between the different transitions is small (A_λ/A_K = 0.93–1.16, see Table 4). On the contrary, the H₂ 0–0 S(3) transition at 9.67  μm lies close to the peak of the silicate absorption, so its flux is more affected by extinction than the rest of the mid-IR H₂ lines. This is clear for UGC 05101, whose 0–0 S(1)/0–0 S(3) ratios are 2.3 ± 0.4 and 0.8 ± 0.1 before and after the extinction correction. Actually, all the galaxies (except NGC 34, see below) have very similar 0−0 S(1)/0–0 S(3) ratios, 0.76 ± 0.06, after the extinction correction based on S_Si.

NGC 34 has a strong silicate absorption (S_Si = −1.10) and an observed 0–0 S(1)/0–0 S(3) = 0.77 ± 0.07, which is compatible with the average extinction corrected ratio of the other galaxies. Therefore, it is possible that the regions producing the H₂ and ~10  μm continuum emissions in NGC 34 suffer from very different extinction levels. Moreover, the NGC 34 uncorrected H₂ SLED does not have any indication of attenuation of the 0–0 S(3) transition (see Fig. 5). This is at odds with the UGC 05101 (S_Si = −1.50) SLED that shows a clear deficit of 0–0 S(3) emission. As a consequence, we do not apply any correction to NGC 34 since we do not have an accurate measurement of extinction toward the H₂ emitting regions.

3. Radiative transfer models

Previous works have analyzed the Herschel mid-J CO SLEDs of galaxies by fitting non-LTE radiative transfer models (e.g., Rangwala et al. 2011; Spinoglio et al. 2012; Pereira-Santaella et al. 2013; Rigopoulou et al. 2013). In general, a model with two components at different temperatures is used. In these models, the “cold” component (usually T_kin< 100 K) reproduces the lowest J CO lines, whereas the “warm” component (T_kin several hundred K) accounts for the higher J emission.

Similarly, when the mid-IR rotational H₂ transitions are analyzed, models with two components are used (e.g., Rigopoulou et al. 2002; Roussel et al. 2007). In these studies, the derived T_kin are always higher than ~100 K since at lower temperatures the rotational H₂ emission is negligible. In addition, it is assumed that H₂ is in LTE because the critical density of these lines is low (n_H2< 10⁴ cm^-3, see Fig. 3 of Roussel et al. 2007).

For the near-IR ro-vibrational H₂ lines, the molecular gas is usually assumed LTE too, although the critical densities of these transitions are higher (n_crit> 10⁷ cm^-3). The derived temperatures from these lines are around 2000 K. In general, it is concluded that they are collisionally excited and UV excitation has relatively low importance for the excitation of the H₂ν = 1–0 transitions (e.g., Davies et al. 2003; Bedregal et al. 2009).

In this work, we analyze the SLED of low- and mid-J CO transitions as well as pure rotational and ro-vibrational H₂ transitions. Therefore, to explain the emission of all these lines, our model would need at least four components with different temperatures ranging from ~50 K to 2000 K. Instead, we use an alternative approach assuming a continuous distribution of temperatures. This is advantageous since we avoid using a somewhat arbitrary number of model components to fit the data.

It is common to find that the observed SLEDs have a positive curvature, that is, the rotational temperature derived from a pair of consecutive transitions increases with the energy of their upper levels (e.g., Neufeld 2012). Therefore, emission could come from gas with a continuous distribution of temperatures instead of several discrete temperatures. To model this distribution of temperatures we assumed that the column density follows a power-law (dN ~ T^{− β}dT). A similar method has been used to model the ro-vibrational and rotational H₂ emission in shocked gas (e.g., Brand et al. 1988; Neufeld & Yuan 2008; Shinn et al. 2009), or the CO emission (e.g., Goicoechea et al. 2013).

In the following, we describe our model in detail. First, we explain how the models for single temperature and density are obtained, and then how they are combined to reproduce a power-law temperature distribution.

Fig. 3 Rest-frame SINFONI integrated spectra. For clarity the spectra are shifted by the following factors 0, 0, 3.5, 4.0, and 4.0 × 10-11 erg cm-2 s-1 from bottom to top. The wavelengths of the H and H2 transitions are indicated.

3.1. Single temperature non-LTE model

To solve the radiative transfer equations for the molecular emission we used the code RADEX (van der Tak et al. 2007). This code uses the escape probability approximation to obtain the molecular level populations and the intensities of the emission lines.

We constructed two sets of models, one for CO and other for H₂, covering a wide range of kinetic temperature (T_kin = 10–2800 K) and H₂ density (n_H2 = 10²−10⁹ cm^-3). For the H₂ grid, we also consider collisions with atomic H with n_H/n_H2 ratios between 1 and 10^-5. In these grids, we assume the LTE H₂ ortho-to-para ratio, which varies between 0 for low temperatures and 3 for T_kin> 200 K (see Fig. 4 of Burton et al. 1992).

The H₂ transitions are optically thin for N_H2/ Δv lower than ~10²⁴ cm^-2 (km s^-1)^-1, thus in our models we assume optically thin H₂ emission³. Similarly, we assume optically thin mid-J CO emission, this condition is fulfilled when N_CO/ Δv is lower than ~10¹⁵ cm^-2 (km s^-1)^-1. In general, this is the case for the mid-J CO observations of galaxies (e.g., Kamenetzky et al. 2012; Pereira-Santaella et al. 2013).

For the CO grid, we used the collisional rate coefficients of Yang et al. (2010) expanded by Neufeld (2012) for high-J levels (J = 41 to 80). This is important for the high temperature models as a non negligible fraction of the CO molecules are in energy levels with J> 40 that otherwise are ignored.

We created two sub-grids of H₂ models. One for high-temperature (T_kin> 100 K) and another for low temperature. For the high-temperature sub-grid we used the H₂-H₂ collisional rate coefficients of Le Bourlot et al. (1999) and the H₂-H coefficients of Wrathmall et al. (2007) that cover kinetic temperatures between 100 and 6000 K. For the lower temperature sub-grid (T_kin< 100 K) we used the H₂-H₂ collisional coefficients of Lee et al. (2008). They only include the lowest nine rotational energy levels of H₂ (E_upper< 5900 K), but at T_kin< 100 K, no molecules are expected at higher rotational energy levels. Collisions with atomic H are mostly important for vibrational excited H₂ levels, so they are neglected for the low-temperature sub-grid. The grids for ortho-H₂ and para-H₂ were constructed separately, although they were later combined assuming the LTE H₂ ortho-to-para ratio.

3.2. Power-law temperature distribution

We want to obtain the level populations for gas following a power-law temperature distribution, dN ~ T^{− β}dT. From the grids described in Sect. 3.1, we can obtain the fractional population of the molecular level i, n_i(T,n_H2), for a given kinetic temperature and H₂ density. Therefore, the total column density of molecules in the energy level i for a kinetic temperature distribution following a power-law between T₀ and T₁ is (1)where is a normalization constant and N is total molecular gas column density⁴. In our models, we use T₀ = 20 K and T₁ = 3500 K. This temperature range is required to model the lower-J CO and the near-IR H₂ transitions, respectively.

In Sect. 4, we show that it is not possible to fit simultaneously the SLEDs of both CO and H₂ using a single power-law temperature distribution. We find that the temperature distribution is much steeper for H₂ (β ~ 5) than for CO (β ~ 3). Therefore, for the SLED fitting we used a broken power-law model with two exponents (β₁ and β₂). This can be modeled with the following equation: (2)where T_b is the break temperature, and .

4. CO and H₂ SLED fitting

To fit the SLEDs, we converted the observed fluxes into beam averaged column densities using the following relation (3)where F_ij is the flux of the transition between levels i and j, ν_ij and A_ij are its frequency and Einstein coefficient, respectively, h is the Planck constant, and Ω is the beam solid angle. Since our sources are barely resolved by Herschel we assume a beam FWHM of 18′′, which is the SSW beam size. This beam size is also used for the rest of measurements as they should include most of the emission in this area (see Sect. 2.5 for details). Using this relation, we assume that the emission is optically thin, which is reasonable for the H₂ transitions as well as for the mid-J CO lines (see Sect. 3).

First, we attempted to fit the SLEDs using a power-law temperature distribution. The free parameters of this model are the CO and H₂ column densities, n_H2 and n_H densities, and β. However, it was not possible to obtain a good fit to both the CO and H₂ emissions. Nevertheless, the CO and H₂ SLEDs individually are well reproduced by a power-law model, although each SLED has different β values.

In these models, there is a strong degeneracy between β and n_H2, β decreases for decreasing n_H2 when densities are below the critical density of the transitions (n_H2< 10⁶ cm^-3 for CO and n_H2< 10⁴ cm^-3 for H₂). This is because changes in n_H2 and β produce similar variations in the SLED shape when the gas is not thermalized. In Fig. 4, we plot the χ² values for one of the galaxies where this degeneracy is evident. In addition, it is clear that the CO and H₂ SLEDs are not simultaneously reproduced by any pair of β and n_H2 parameters. For our galaxies, we find that β is always higher for H₂ (steeper temperature distribution) for any given n_H2. The β values range between 2.5 and 3.5 for CO and between 4.5 and 5.5 for H₂ for n_H2 higher than the critical densities of the transitions.

Fig. 4 χ2 values for the power-law temperature distribution model as a function of β and nH2 for NGC 7130. The red (green) contours are the 1, 2, and 3σ confidence levels for the CO (H2) SLED.

As anticipated in Sect. 3.2, to account for the difference on the β values for the CO and H₂ SLEDs, we used an alternative model consisting of a broken power-law, which has two exponent values, β₁ and β₂, for temperatures below and above the threshold temperature T_b, respectively. With this model, it is possible to fit all the available CO and H₂ transitions for our galaxies, with the exception of the CO  transition, which is always underpredicted by the model. The later can be caused by the presence of low temperature molecular gas that would not follow the power-law distribution. Consequently, we included a cold gas component at 15 K assuming optically thin emission for the CO  transition. In general, CO  emission is not optically thin, so the obtained column densities are lower limits.

The value of the breaking temperature (T_b) is not well constrained because the upper level energies of the CO and H₂ transitions do not overlap. When T_b is left as a free parameter it varies between 100 and 500 K with 1σ uncertainties ~100 K, and its variation affects the β values slightly. The average T_b is 200 ± 130 K. It is not possible to accurately determine T_b with our data, therefore, we fixed it to the average value.

As discussed above, β₁ and n_H2 are strongly correlated, specially for CO, when n_H2 is below the critical densities. The CO column density (N_CO) is relatively constant when β₁ and n_H2 vary. But, on the contrary, the H₂ column density variation is large. This is because we only measure the warm H₂ column density directly, which represents a small fraction of the total H₂, and the total H₂ column density is extrapolated using the β₁ exponent. Therefore, in our model, lower n_H2 also implies lower H₂ column densities (and higher CO abundances) because of the smaller value of β₁.

We used this to estimate a lower limit for n_H2 given that the CO abundance should be lower or equal to that of C. Assuming the solar atomic C abundance (2.4 × 10^-4; Lodders 2003), the lower limit for n_H2 is ~10^4.0±0.5 cm^-3, depending on the galaxy. In the case of sub-solar abundances in these galaxies, the n_H2 upper limit would be slightly higher. Similarly, we can obtain a lower limit for the CO abundance by setting n_H2 to a value above the CO critical density (e.g., 10⁶ cm^-3). The minimum CO abundance ranges between 10^-5 and 10^-6. In Table 7, we list the best-fit parameters for n_H2 = 10^4.5 cm^-3 and n_H2 = 10⁶ cm^-3. We plot the best-fit models in Fig. 5.

The n_H value depends mostly on the ratio between the intensities of the rotational and near-IR ro-vibrational H₂ lines. For our galaxies we obtain n_H values between 10³–10⁴ cm^-3. The critical densities for collisions with H of the near-IR H₂ transitions are n_H ~ 10⁴ − 10⁵ cm^-3 at 1500 K, much lower than the critical densities for collisions with H₂. Therefore, a relatively small abundance of atomic H can enhance the luminosity of the near-IR transitions. For the assumed n_H2 = 10^4.5 and 10⁶ cm^-3 the ro-vibrational near-IR emission would be one or two orders of magnitude lower than the observed if we did not include collisions with H in our model.

Fig. 5 CO and H2 rotational diagrams. The red (blue) line is the best-fit model for the warm (cold) component. The number close to the data points is the J quantum number of the level. For H2, circles correspond to the vibrational level ν = 0, diamonds to ν = 1, and squares to ν = 2.

5. Cold-to-warm molecular gas ratio

To estimate the cold molecular gas mass, we used the CO to H₂ conversion factor from Leroy et al. (2011), N_c(H₂)(cm^-2) = 2 × 10²⁰I_{COJ = 1 − 0} (K km s^-1). We compare this cold molecular gas column density with the warm molecular gas column density obtained from the fit of the CO and H₂ SLEDs (Table 7). The warm molecular gas fraction ranges between 1 and 20% for the models with n_H2 = 10^4.5 cm^-3 and between 15 and >100% for the models with n_H2 = 10⁶ cm^-3. This suggests that for the two galaxies where the warm column density exceeds the cold column density (UGC 05101 and NGC 5135) the model with n_H2 = 10⁶ cm^-3 overestimates the N_w(H₂). Thus, in the following we use the n_H2 = 10^4.5 cm^-3 model for these two galaxies. Moreover, these are the only two sources where the χ² is lower for the n_H2 = 10^4.5 cm^-3 model.

For our sample of IR bright galaxies it might be adequate to use the lower CO to H₂ conversion factor N_H2 (cm^-2) = 0.5 × 10²⁰I_{COJ = 1 − 0} (K km s^-1) measured by Downes & Solomon (1998) in active galaxies. Therefore, the cold-to-warm ratio could be lower by a factor of 4.

6. Molecular gas heating

In the previous sections, we showed that the CO and H₂ SLEDs are compatible with a broken power-law temperature distribution, however, we did not discuss how the molecular gas is heated. Three main mechanisms are invoked to explain the heating of the molecular gas: UV radiation, shocks, and X-ray and cosmic rays. In this section, we compare the molecular gas temperature distribution expected in photodissociation region (PDRs; UV radiation) and shocks regions with the observed temperature distribution.

X-ray and cosmic ray heating can be important for AGNs, however, in this work we analyze active galaxies with strong star-formation (see Table 1). The AGN contribution to the mid-J CO emission seems to be small compared to that of SF (e.g., Pereira-Santaella et al. 2013). Likewise, SF dominates the H₂ emission in LIRGs (e.g., Roussel et al. 2007; Pereira-Santaella et al. 2010). For these reasons, we consider that X-ray heating should not contribute significantly in the integrated spectra of these galaxies.

Fig. 6 Left panel: chemical structure and temperature of a PDR with G0 = 104 and nH = 104 cm-3. The thick black line is the gas temperature (left axis). The other lines represent the abundances (left axis) of H (black), H2 (orange), C+ (blue), C i (violet), and CO (green). Right panel: H2 (orange) and CO (green) column densities as a function of the temperature. The dashed lines are the best power-law fits to the distribution (see Sect. 6.1).

6.1. UV radiation: PDRs

We used the radiative transfer code Cloudy version 13.02 (Ferland et al. 2013) to estimate the gas temperature distribution of PDRs. We used a constant density slab model illuminated by the interstellar radiation field of Black (1987). We created a grid of models with n_H between 10⁴ and 10⁶ cm^-3 and radiation field intensity between 10¹ and 10⁴G₀, where G₀ = 1.6×10^-3 erg cm^-2 s^-1 is the intensity integrated between 6 and 13.6 eV. Since we are interested only in the PDR region we extinguished the input spectra to remove the ionizing radiation.

We extracted the CO and H₂ abundances and kinetic temperatures as a function of the column density from the output of the models, which we used to compute the temperature distribution of CO and H₂. Then we fitted these distributions with a power-law function. In the fit we only included gas with T> 20 K and T> 150 K, for CO and H₂, respectively, because lower temperature gas does not produce significant emission of the observed transitions.

In Fig. 6, we show an example Cloudy PDR model and the fit to the temperature distribution. We find that, in general, H₂ distributions are steeper than those of CO. The power-law index⁵ for CO is ~3–10 for n_H2 = 10^4.5 cm^-3, and ~5–7 for n_H2 = 10⁶ cm^-3, and for H₂β ranges between ~10 and 15. That is, these β values are higher than the values derived from our fits to the observed SLEDs (Sect. 4).

6.2. Shocks

To calculate the molecular gas temperature distribution in shocks we used a simple model. We assumed that molecular gas is heated at a constant rate up to a certain temperature, which depends on the shock velocity and is high enough to produce H₂ emission, and after that it is cooled radiatively by the dominant cooling molecule.

The main coolant of warm molecular gas is H₂, CO, or H₂O depending on the density and temperature of the gas (Neufeld & Kaufman 1993). For our galaxies, the H₂O transitions are weaker than the CO transitions for those with similar E_up values. That is, H₂O cooling does not seem to dominate the cooling for the average conditions in these galaxies, and therefore, we do not consider H₂O in our model.

First, we calculated the non-LTE cooling functions of CO and H₂ as a function of the gas density and kinetic temperature using the grid of models described in Sect. 3.1. For the H₂ cooling function we do not include the effect of collisional excitation by H since it would only affect the cooling of high-temperature gas through ro-vibrational emission of H₂.

The temperature dependence of these cooling functions can be approximated with a power-law (e.g., Draine & McKee 1993). Thus we fitted the calculated cooling functions using a power-law in the temperature range between 30 and 150 K for CO, and 150 and 1000 K for H₂. These are the temperature ranges where CO and H₂ dominate the gas cooling (Neufeld & Kaufman 1993).

Fig. 7 Power-law indices of the temperature dependence of the CO (top) and H2 (bottom) non-LTE cooling functions (solid line; see Sect. 6.2). The fits are valid for a temperature range between 30 and 150 K for CO and 150 and 1000 K for H2. The fitted β values form the models with nH2 = 104.5 cm-3 and nH2 = 106 cm-3 are also plotted for each galaxy. For visualization the abscissa of the data points are slightly offset.

Then, assuming that the molecular gas cools radiatively after being heated according to these fitted cooling functions (i.e., dU/ dt ~ − T^γ, where U = 5/2 kT is the internal kinetic energy), the observed temperature distribution of the molecular gas is dN/ dT ~ T^{− γ}. That is, the β we fitted for the SLEDs should be compared with γ. The power-law indices, γ, for CO and H₂ as a function of the density are shown in Fig. 7.

6.3. Comparison with the observations

In our fits, we obtained that the power-law index varies between 1.2 and 3.5 for CO and between 4.3 and 4.9 for H₂. In general, these values are lower than those expected in PDR, but higher than those expected from shocks (see Table 7 and Fig. 7). Therefore, the combination of PDRs and shocks are needed to explain the intermediate power-law indices.

The CO emission is compatible with that predicted by shock models only for NGC 5135⁶, whose best-fit model has n_H2 = 10^4.5 cm^-3. This result, however, strongly depends on the adopted density. For instance, for the model with n_H2 = 10⁶ cm^-3, which provides only a slightly worse fit to the CO SLED of NGC 5135, the CO emission would need a PDR component.

If we combine these PDR and shock models, the PDR SLEDs would dominate the lower-J CO transitions because of their steeper power-law index, while the contribution from shocks would be dominant for higher-J CO transitions. Moreover, only gas with T> 100 K produces H₂ emission so, depending on the G₀ values of the PDRs and the shock velocities, the contributions of PDRs and shocks may differ in the CO and H₂ SLEDs.

Our sample includes three spiral galaxies (NGC 5135, NGC 7130, and NGC 7469) and three major mergers (NGC 34, IRASF 05189–2524, and UGC 05101). The latter have lower β₁ values for their CO SLEDs (Fig. 7) suggesting that the relative contribution from shocks to the CO emission in these mergers is larger than in the spirals. For the H₂ emission only IRASF 05189–2524 has a β₂ value lower than the rest of galaxies, so in this source, shocks able to heat the molecular gas up to T> 100 K could be important. Also, this galaxy is the most luminous in our sample and the galaxy with the higher AGN contribution (Table 1), so the bright AGN might affect its CO and H₂ emissions.

6.4. Comparison with two-component models

The CO SLEDs of two of our galaxies (UGC 05101 and NGC 7130) were presented by Pereira-Santaella et al. (2013). In that paper, a model consisting of two components was used to fit the data: a cold component that contributed mainly to the CO emission, and a warm component with variable H₂ density and kinetic temperature.

Once corrected for the different beam sizes used, the beam-averaged CO column densities derived by both methods agree within the uncertainties, although, this could be a particular case due to the specific SLEDs of these two galaxies.

The interpretation of the n_H2 and T_kin of the two-component model (TCM) is not so straightforward. The high kinetic temperature needed by the TCM reveals the presence of warm molecular gas, but the exact temperature value is not necessarily a representative temperature of the molecular gas. The obtained n_H2 and T_kin are biased toward the conditions in the most luminous regions. Actually, the luminosity of each CO transition, for a given molecular mass, depends on n_H2 and T_kin. So, the derived n_H2 and T_kin will depend on the unknown distribution of densities and temperatures of the molecular gas. This problem is partially solved assuming a functional form for the dN/ dT, as we adopted in our models.

In addition, n_H2 and T_kin are correlated (lower n_H2 and higher T_kin, and higher n_H2 and lower T_kin produce similar CO SLEDs) and this adds an extra uncertainty to the derived values. This is equivalent to the uncertainty due to the β-n_H2 correlation for our models (see Sect. 4).

Although the TCM provides useful information, when the analyzed molecular transitions comprise a wide range of excitation temperatures, the need for additional components for the SLED fitting can make the modeling contrived. Instead, the power-law models used here could better represent the integrated temperature distribution of the molecular gas in a galaxy.

7. Conclusions

We have studied the integrated CO and H₂ emission of six local IR bright galaxies using non-LTE models. Assuming a broken power-law distribution for the molecular gas temperatures, our model reproduces both the CO SLED (from J_up = 1 to J_up = 13) and the H₂ SLED (J_up ≤ 7 for the lowest three vibrational levels) in our sample of galaxies.

The main findings of this work are summarized in the following:

• 1.

  With a single power-law temperature distribution it is notpossible to fit simultaneously the CO andH₂ SLEDs. The H₂ SLEDs have amuch steeper power-law index(β₂ ~ 4–5) than the CO SLED (β₁ ~ 1–3). This is the expected behavior for the temperature distributions in PDRs and shocks.

• 2.

  We found that for most of the galaxies, the models with n_H2 = 10⁶ cm^-3 provide the best fit to the observed data, thus the majority of the transitions are close to LTE. The minimum acceptable density for the warm gas is ~n_H2 = 10^{4 ± 0.5} cm^-3, lower densities would imply CO abundances higher than the atomic C abundance. Likewise, we obtained that the minimum CO abundance in the warm gas is x_CO ~ 10^-6 − 10^-5 assuming LTE conditions.

• 3.

  The column densities of the warm molecular gas (T> 20 K) represents between 10 and 100% of the molecular gas traced by the CO transition.

• 4.

  We used PDR and shock models to determine the excitation mechanism of the molecular gas. Our models show that the temperature distributions are steeper for PDRs than for shocks. We also found that the temperature distribution of the warmest gas (T> 100 K) emitting in H₂ is steeper than that of coldest gas (T> 30 K), which produces the mid-J CO emission for both PDR and shocks models. This is because of the different main coolant of the warm and cold molecular gas (H₂, and CO, respectively).

• 5.

  Neither PDR nor shocks alone can explain the derived temperature distributions. A combination of both is needed to reproduce the observed SLEDs. In this case, the lower-J CO transitions would be dominated by PDR emission, whereas the higher-J CO transitions would be dominated by shocks.

• 6.

  The three major mergers among our targets (NGC 34, IRASF 05189–2524, and UGC 05101) have shallower temperature distributions for CO than the other three spirals. This suggests that the relative contribution of shocks to the heating of warm molecular gas (T< 100 K) in these major mergers is higher than in the other three spirals.

• 7.

  For only one of the mergers, IRASF 05189–2524, the shallower H₂ temperature distribution (hot molecular gas) suggests that the relative importance of shocks is high. This galaxy also has a bright AGN that dominates the bolometric luminosity, which can contribute to the molecular gas heating. For the other two mergers, the H₂ temperature distribution is similar to that of the spiral galaxies. Therefore the shocks producing the extra contribution to the CO emission in these mergers are not able to heat the molecular gas to temperatures higher than 100 K, which would be necessary to see the H₂ emission.

Acknowledgments

We thank anonymous referee for comments that improved the paper. This work was funded by the Agenzia Spaziale Italiana (ASI) under contract I/005/11/0. J.P.L. acknowledge support from the Spanish Plan Nacional de Astronomía y Astrofísica AYA2010-21161-C02-01 and AYA2012-32295. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

All Tables

All Figures

	Fig. 1 Observed SPIRE/FTS spectra of our sample. The black and gray lines are the SSW and SLW spectra, respectively. Notice the overlap between the two spectra in the 960 and 990 GHz spectral range. The dashed red lines mark the position of the detected lines. The error bars indicate the median 1σ uncertainty of each spectrum.
In the text

	Fig. 1 continued.
In the text

	Fig. 2 Rest-frame Spitzer/IRS spectra. The short-low (5.5–14  μm), short-high (10–19  μm), and long-high (19–36  μm) spectra are plotted for each galaxy. For clarity, the spectra are multiplied by the following factors 1, 2.5, 14, 100, 220, and 2000, from bottom to top. The wavelengths of the H2 rotational transitions are indicated.
In the text

	Fig. 3 Rest-frame SINFONI integrated spectra. For clarity the spectra are shifted by the following factors 0, 0, 3.5, 4.0, and 4.0 × 10-11 erg cm-2 s-1 from bottom to top. The wavelengths of the H and H2 transitions are indicated.
In the text

	Fig. 4 χ2 values for the power-law temperature distribution model as a function of β and nH2 for NGC 7130. The red (green) contours are the 1, 2, and 3σ confidence levels for the CO (H2) SLED.
In the text

	Fig. 5 CO and H2 rotational diagrams. The red (blue) line is the best-fit model for the warm (cold) component. The number close to the data points is the J quantum number of the level. For H2, circles correspond to the vibrational level ν = 0, diamonds to ν = 1, and squares to ν = 2.
In the text

Fig. 6 Left panel: chemical structure and temperature of a PDR with G0 = 104 and nH = 104 cm-3. The thick black line is the gas temperature (left axis). The other lines represent the abundances (left axis) of H (black), H2 (orange), C+ (blue), C i (violet), and CO (green). Right panel: H2 (orange) and CO (green) column densities as a function of the temperature. The dashed lines are the best power-law fits to the distribution (see Sect. 6.1).

In the text

Fig. 7 Power-law indices of the temperature dependence of the CO (top) and H2 (bottom) non-LTE cooling functions (solid line; see Sect. 6.2). The fits are valid for a temperature range between 30 and 150 K for CO and 150 and 1000 K for H2. The fitted β values form the models with nH2 = 104.5 cm-3 and nH2 = 106 cm-3 are also plotted for each galaxy. For visualization the abscissa of the data points are slightly offset.

In the text