STAR FORMATION RELATIONS AND CO SPECTRAL LINE ENERGY DISTRIBUTIONS ACROSS THE J-LADDER AND REDSHIFT

The pan-chromatic (FUV/optical to radio) spectral energy distributions (SEDs) of our sample galaxies were modeled using CIGALE (Code Investigating GALaxy Emission; Burgarella et al. 2005; Noll et al. 2009). CIGALE employs dust-attenuated stellar population models to fit the FUV/optical SED, while at the same time ensuring that the dust-absorbed UV photons are re-emitted in the FIR, thus maintaining energy balance between the FUV and FIR. The FIR/submm continuum is modeled using the templates by Dale & Helou (2002) and Chary & Elbaz (2001). The Salpeter initial mass function was used for the stellar emission population synthesis models from Maraston (2005), and for the reddening we used attenuation curves from Calzetti et al. (1994) with a wide range of V-band attenuation values for young stellar populations. Despite having carefully checked our samples against AGNs, we allowed for the possibility of additional dust emission from deeply buried AGNs by including the 32 AGN models from the Fritz et al. (2006) library in our SED fits. Reassuringly, the AGN fraction never exceeded 20% of the total IR luminosity. Excellent fits were obtained for all the local galaxies due to their well-sampled SEDs. For the high-z galaxies, only sources with data points longward and shortward of (or near) the dust peak (λ_rest ∼ 100 μm) were included in the final analysis, which was a total of 35 out of the original 74 DSFGs. All SED fits used in this paper can be found at http://demogas.astro.noa.gr, and will also be presented in a forthcoming paper (E. M. Xilouris et al., in preparation).

From the SED fits we derived the IR (L_IR, from 8 μm to 1000 μm rest-frame), and the FIR (L_FIR, from 50 μm to 300 μm rest-frame), luminosities of our sample galaxies (Tables 1 and 2). We use the latter for our analysis in order to minimize the effects of AGNs, which are strongest in the mid-IR regime (i.e., ∼8–40 μm). In addition, the mid-IR is rich in PAH emission/absorption features, which could affect L_IR estimates. For the uncertainty on our IR/FIR luminosity estimates, we adopted the 1σ dispersion of the luminosity distributions obtained through bootstrapping the photometry errors 1000 times. Typical uncertainties, δL_FIR, were ∼20% and ∼40% for the local and high-z samples, respectively, and were adopted across the board for the two samples. We stress that the above FIR luminosities are total luminosities, that is, derived from aperture fluxes that encompass the full extent of the galaxies, and thus match the CO measurements.

Figure 2 shows the separate $\log L_{\rm FIR}\hbox{--}\log L^{\prime }_{\rm CO}$ relations (where log  is for base 10) for each CO transition (from CO J = 1–0 to J = 13–12) for the galaxy samples analyzed here. Highly significant correlations are seen in all transitions, as given by their near unity linear correlation coefficients (r; see Figure 2). Even for the highest transition (J = 13–12), where the dynamical range spanned in luminosities is relatively small, we see a statistically significant correlation. To ensure that the observed correlations are not simply due to both L_IR and $L^{\prime }_{\rm CO}$ being $\propto D_{\rm L}^2$ (the luminosity distance squared), we calculated for each correlation the partial Kendall τ-statistic (Akritas & Siebert 1996) with $D_{\rm L}^2$ as the test variable. In all cases (up to J = 13–12), we find probabilities P < 10⁻⁶ that the observed FIR−CO correlations are falsely induced by the fact that ${luminosity} \propto D_L^2$.

Zoom In Zoom Out Reset image size

A function of the form $\log L_{\rm FIR} = \alpha \log L^{\prime }_{\rm CO} + \beta$ was adopted to model the correlations, and the optimal values of the model parameters (α and β) were fitted (to this end we used the IDL routine linmix_err; Kelly 2007). The slopes (α) and intersection points (β) inferred from fits to the combined low- and high-z samples are given in Table 3, along with the scatter (s) of the data points around the fitted relations. The fits are shown as dashed lines in Figure 2. The best-fit (α, β)-values obtained by using L_IR instead of L_FIR are also listed in Table 3. Within the errors, the fitted parameters are seen to be robust against the adopted choice of FIR or IR luminosity. Also, our results did not change in any significant way when omitting the lensed DSFGs from the analysis. Often lensing amplification factors are uncertain, and strong lensing cannot only skew the selection of sources toward more compact (and thus more likely warm) sources, but for a given source it may also boost the high-J CO lines relative to the lower lines (this is discussed further in Section 6.2).

Figures 3 and 4 show the slopes and normalizations, respectively, of the $\log L_{\rm FIR}\hbox{--}\log L^{\prime }_{\rm CO}$ relations derived above as a function of the critical densities probed by the various CO transitions. The critical densities are calculated as n_crit = A_ul/∑_{i ≠ u}C_ui, where A_ul is the Einstein coefficient for spontaneous decay, and ∑_{i ≠ u}C_ui is the sum over all collisional coefficients (with H₂ as the collisional partner) out of the level u, upward and downward (see Table 4 where, as a reference, we also list n_crit-values for a number of HCN and CS transitions). Although it is the first three levels up or down from the u-level (i.e., |u − i| < 3) that dominate the sum, often in the literature molecular line critical densities are only calculated for a two-level system (i.e., |u − i| = 1) or for the downward transitions—both practices can significantly overestimate the true n_crit for a given transition. The collision rates were adopted from the Leiden Atomic and Molecular Database (LAMDA; Schöier 2005) for T_k = 40 K, which is within the range of typical dust and gas temperatures encountered in local (U)LIRGs and high-z DSFGs (e.g., Kovács et al. 2006). We do not correct for optical depth effects (i.e., line-trapping) as these are subject to the prevailing average ISM conditions, but we note that large optical depths (especially for low-J CO and HCN lines) can significantly lower the effective critical density to $n^{(\beta)} _{\rm crit}=\beta _{ul} n_{\rm crit}$, where β_ul is the average line escape probability (=[1 − exp (− τ_ul)]/τ_ul for spherical geometries). The collisional excitation of CO to higher rotational states (E_J) is set not only by the gas density but also by its kinetic temperature. The minimum temperature (T_min) required for significant collisional excitation of a given rotational state is approximately given by ∼E_J/k_B = B_rotJ(J + 1)/k_B, where B_rot is the rotational constant of CO, and k_B is the Boltzmann constant. As a rule of thumb, high kinetic temperatures are needed in order to excite the high-J CO lines (see Table 4), although due to the n–T_k degeneracy, this can also be achieved for very dense, low-temperature gas.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two trends regarding the $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ relations become apparent from Figures 3 and 4 (see also Table 3). First, the slopes are linear for J = 1–0 to J = 5–4, but then become increasingly sublinear the higher the J level. Second, the normalization parameter β remains roughly constant (∼2) up to J = 4–3, 5–4, but then increases with higher J level, reaching β ∼ 8 for J = 13–12, which for a given CO luminosity translates into ∼6 orders of magnitude higher L_FIR. We stress that although the $L_{\rm FIR}-L^{\prime }_{\rm CO}$ relations are linear, and β is roughly constant up to J = 5–4, it does not generally imply that the CO lines are thermalized (i.e., $L^{\prime }_{\rm CO_{\rm J,J-1}}/L^{\prime }_{\rm CO_{1,0}}\simeq 1$) up to this transition. There is significant scatter within the samples, and while a few sources do have nearly thermalized J = 2–1, 3–2, and/or 4–3 lines, in general, $L^{\prime }_{\rm CO_{\rm J,J-1}}/L^{\prime }_{\rm CO_{1,0}} \lesssim 1$. In fact, rewriting the $L_{\rm FIR}-L^{\prime }_{\rm CO}$ relations as

$$\begin{equation} L^{\prime }_{\rm CO_{\rm J,J-1}}/L^{\prime }_{\rm CO_{\rm 1,0}} = L_{\rm FIR}^{\alpha _{\rm J,J-1}^{-1} - \alpha _{\rm 1,0}^{-1}} \times 10^{\beta _{\rm 1,0} - \beta _{\rm J,J-1}}, \end{equation} \tag{ 1 }$$

and inserting the fitted values from Table 3 yields $L^{\prime }_{\rm CO_{\rm J,J-1}}/L^{\prime }_{\rm CO_{\rm 1,0}} < 1$ over the range L_FIR = 10^{9 − 14} L_☉.

In the following sections, we discuss these empirical relations in the context of existing theoretical models and the new observational studies of such relations using dense gas tracers like HCN and CS.

Before comparing our derived FIR−CO slopes with those from the literature, we must add two cautionary notes: (1) many studies examine the $L^{\prime }_{\rm CO}\hbox{--}L_{\rm FIR}$ relation rather than $L_{\rm FIR}-L^{\prime }_{\rm CO}$, and one cannot compare the two simply by inferring the inverse relation; (2) often only the errors in one variable (typically $L^{\prime }_{\rm CO}$) are taken into account when fitting such relations, when in fact the uncertainties in both L_FIR and $L^{\prime }_{\rm CO}$ must be considered (see Mao et al. 2010 for a further discussion). Failing to do so can result in erroneous estimates of the slope.

For these reasons, we have refitted the data from a number of studies (see below) using the method described in Section 3, that is, with errors in both L_FIR and $L^{\prime }_{\rm CO}$ and including only sources with L_FIR ≳ 10¹¹ L_☉ in our analysis. Finally, not all studies use the FIR definition used here to infer L_FIR, and other studies use the full (8–1000 μm) luminosity. These differences can result in a different overall normalization (i.e., β), but are not expected to affect the determination of α (see Table 3 where there is little change in α when switching between L_{FIR[50-300 μm]} and L_{IR[8-1000 μm]}).

From Figure 3 we note the overall good agreement between the FIR−CO slopes derived here and values from the literature. For CO(1–0), however, one set of measurements found super-linear slopes ($\alpha _{{\rm CO}_{1,0}} \sim 1.3\hbox{--}1.4$; Juneau et al. 2009; Bayet et al. 2010), while most others favor a slope of unity (Gao & Solomon 2004b; Mao et al. 2010; Ivison et al. 2011; this work). Note that our reanalysis of the Yao et al. (2003) and Baan et al. (2008) data revised their slopes from super-linear to linear: $\alpha _{{\rm CO}_{1,0}} = 0.94\pm 0.07$ and 1.08 ± 0.06, respectively (a similar result was found by Mao et al. 2010). Gao & Solomon (2004b) finds a super-linear FIR−CO slope (α = 1.3–1.4) from their entire sample (combining L_IR ∼ 10¹⁰ L_☉ objects with LIRGs and ULIRGs), yet when including only the LIRGs and ULIRGs in the analysis, we obtain a linear slope ($\alpha _{{\rm CO}_{1,0}} = 0.91\pm 0.22$). Refitting the data presented in Juneau et al. (2009) and Bayet et al. (2010) we reproduce their super-linear slopes.

For CO(2–1) our slope of unity is consistent within the errors with Bayet et al. (2010), who find a slightly super-linear slope based on 17 sources. Our reanalysis of the CO(2–1) data by Baan et al. (2008) yields a sublinear slope of α = 0.82 ± 0.10. However, as pointed out by the authors, a nonnegligible fraction of the total CO(2–1) emission is likely to have been missed due to the smaller telescope beam at higher frequencies, thus biasing the FIR–CO relation to a shallower value of α.

In the case of CO(3–2), the existing slope-determinations (Yao et al. 2003; Narayanan et al. 2005; Iono et al. 2009; Bayet et al. 2010; Mao et al. 2010), including our own, are in agreement and favor a value of unity within the errors. This includes a reanalysis of the Yao et al. (2003) data, which yielded α = 1.00 (see also Mao et al. 2010).

For CO(4–3) to CO(7–6) there is agreement within the errors between our results and the slopes found by Bayet et al. (2010), which were determined using observed and model-extrapolated CO luminosities of 7 low-z and 10 high-z sources. A departure from linear toward sublinear is also found by the latter study, albeit we find this turnover to occur at J = 6–5, rather than at J = 4–3 as deduced by Bayet et al. (2010). Here we must note, however, that the use of models to extrapolate to high-J CO luminosities is not safe, and artificial turnovers can be introduced because of the inability of such models—in the absence of appropriate line data—to reliably account for the existence of warmer and denser gas components. This further underscores the value of our observed $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ relations from J = 1–0 to J = 13–12 in safely determining such departures from linearity and/or in normalization before proceeding toward any interpretation based on ISM/SF physics.

The few super-linear slopes of FIR−CO luminosity relations for low-J CO lines that survive careful reanalysis (Juneau et al. 2009; Bayet et al. 2010) could be a byproduct of dense molecular gas being the direct SF fuel in all galaxies (with a constant SFE) and a f_{dense, X} = M_dense/M_X (where X could be the total H₂ gas mass traced by CO J = 1–0, J = 2–1 lines), which varies within the galaxy sample with d(f_dense)/dL_FIR > 0. Variations of this simple scenario have been suggested throughout the literature (Wong & Blitz 2002; Gao & Solomon 2004b), and unlike more sophisticated interpretations of such super-linear slopes offered by the two theoretical treatises available on this matter (Krumholz & Thompson 2007; Narayanan et al. 2008), the only assumption here is that d(f_dense)/dL_FIR > 0. The latter is a well-documented fact as starbursts/mergers, which dominate the high-L_FIR end, are observed to have larger dense/total gas-mass fractions than lower-L_FIR isolated disks (e.g., Gao & Solomon 2004b; García-Burillo 2012).

Within this scheme the further the gas phase X is from the dense, star forming phase in terms of physical conditions and relevance to the star formation, the higher the d(f_dense)/dL_FIR value is and the deduced super-linear slope of the $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ (or corresponding S-K) relation. On the other hand, for galaxy samples with a smaller range of IR–luminosities that have a nearly constant f_{dense, X}—such as the (U)LIRG+DSFG sample considered in this paper (see Section 5.2)—a linear slope of an $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ relation can still be recovered, for example, CO J = 1–0, which traces all metal-rich molecular gas mass rather than only the dense SF one. We return to this point later in our discussion.

As already mentioned in the Introduction, the S-K relation, and especially the important one involving the dense gas component in galaxies, is not easily accessible observationally. So one falls back to the much more observationally accessible proxy $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ relations. Their index can then be linked to that of an assumed underlying S-K relation in galaxies using two theoretical models (Krumholz & Thompson 2007; Narayanan et al. 2008). Both models posit the same intrinsic S-K relation of $\rho _{\rm SFR}\propto \rho ^{1.5}_{\rm gas}$ (or $\Sigma _{\rm SFR}\propto \Sigma ^{1.5}_{\rm gas}$ (for disks of near-constant scale-height), justified under the assumption of a constant gas fraction transformed into stars per free fall time (t_ff∝(Gρ_gas)^−1/2). The same S-K relation emerges if the SF timescale is instead set by the dynamical timescale of a marginally Toomre-stable galactic disk with SF converting a fixed fraction of gas into stars over such a timescale (Elmegreen 2002).

Both models give expected values of α versus n_crit(X) for $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ luminosity relations over a large range of line critical densities, but both models are applicable only for lines that require low temperatures to excite E_J/k_B ≲ 30 K (for the Krumholz & Thompson 2007 model this limit is ∼10 K). The reasons behind this limitation are explicit assumptions about isothermal gas states at a set temperature (Krumholz & Thompson 2007), or the tracking of such states over a small range of T_k ∼ 10–30 K (Narayanan et al. 2008). The low-J lines of heavy rotor molecular line (e.g. HCN and HCO⁺) data sets found in the literature and our low-J CO lines are certainly within the range of applicability of these models. The indices of the corresponding $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ power law relations can thus be compared to these theoretical predictions.

Gao & Solomon (2004a, 2004b) found a linear FIR−HCN correlation for the J = 1–0 transition (see also Baan et al. 2008), which extends from local ULIRGs/LIRGs (L_IR ∼ 10¹¹ to 10¹² L_☉) to normal star forming galaxies (L_IR ∼ 10⁹ to 10¹⁰ L_☉), down to individual Galactic molecular clouds with L_IR ≳ 10^4.5 L_☉ (Wu et al. 2005, 2010). A weakly super-linear FIR−HCN(1–0) slope (α ≃ 1.2) was found by Graciá-Carpio et al. (2008a) and García-Burillo (2012) over a combined sample of local normal galaxies and LIRG/ULIRGs. However, careful analysis by García-Burillo (2012) demonstrated that a bimodal fit (i.e., a different normalization parameter β) is better, with each galaxy sample well-fit by a linear relation. Finally, a weakly super-linear $L_{\rm FIR}\hbox{--}L_{\rm HCN_{1,0}}$ appears when high-z observations of the most IR-luminous starburst galaxies and QSOs are included in the locally established relation (Gao et al. 2007; Riechers et al. 2007). This could be bimodal instead, but with an otherwise linear $L_{\rm FIR}\hbox{--}L_{\rm HCN_{1,0}}$ relation (and with insufficient high-L_FIR objects to decide the issue). A physical reason for such bimodalities is discussed in Section 5.2.

Extending HCN observations to include many more objects in the crucial L_IR > 10¹² L_☉ regime is necessary for deciding such issues. Even then one must eventually obtain the underlying SFR−M_dense relation before arriving at secure conclusions about a varying SFE = SFR/M_dense of the dense gas in (U)LIRGs. The latter is the crucial physical quantity underlying the normalization of such $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm HCN}$ relations and, for example, a rising or bimodal $X_{\rm HCN}=M_{\rm dense}/L_{\rm HCN_{1,0}}$ factor toward high-L_IR systems can easily erase purported SFE trends obtained by using single-line proxies of dense gas.

A sublinear FIR−HCN slope (∼0.7–0.8) for the J = 3–2 transition has been reported (Bussmann et al. 2008; Juneau et al. 2009), but is very likely biased low due to not having performed any beam correction (see Section 2) for some of their very nearby extended objects, where the HCN beam does not cover the entire IR emitting region. For these sources, the HCN measurements do not match the IR luminosities, and since they reside at the lower end of the HCN luminosity distribution, the net effect will be to bias the relation toward shallower values. For this reason we have chosen to ignore the sublinear FIR−HCN slopes from Bussmann et al. (2008) and Juneau et al. (2009). A recent survey of HCN J = 4–3 and CS J = 7–6 (Zhang et al. 2014), as well as CS J = 1–0, 2–1, 3–2, and 5–4 (Z.-Y. Zhang et al., in preparation), toward nearby star forming galaxies (L_IR ∼ 10⁹–10¹² L_☉), where such effects have been adequately accounted for, establishes a slope α ∼ 1 for these transitions (see also Wu et al. 2010 and Wang et al. 2011). Many of these CS transitions have higher critical densities than HCN, and CS is furthermore less prone to IR pumping effects than HCN (CS is pumped at 7.9 μm compared to 14 μm for HCN). Pumping of HCN (and also HNC), however, typically only becomes important at dust temperatures ≳ 50 K (Aalto et al. 2007), and would typically require even higher temperatures for CS. This is important because pumping could affect the CS/HCN luminosities (especially at high-J), and thus result in linear IR–CS/HCN relations. In Figure 3 we summarize all the observationally determined $L_{\rm IR}\hbox{--}L^{\prime }_{\rm HCN}$ and $L_{\rm IR}\hbox{--}L^{\prime }_{\rm CS}$ slopes from the literature, along with those derived from our CO lines. Overall, the data suggest α ∼ 1 for our $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ relations from J = 1–0 up to J = 5–4, 6–5, and for the HCN and CS lines. The latter cover a range of n_crit ∼ 10⁴–10⁷ cm⁻³ (i.e., reaching up well into the high-density regime of the star forming gas phase).

From Figure 3 it becomes clear that most observations are incompatible with current model predictions (shown as the gray shaded area) both for the heavy rotor and the low-J CO lines (where such models remain applicable). Super-linear slopes do appear for some CO J = 1–0 data sets but then, unlike model predictions, the slopes remain linear for lines with much higher critical densities, including those of mid-J CO lines J = 3–2, 4–3, 5–4 (the FIR–CO luminosity relation for J = 6–5 is also compatible with a linear one).

In summary, we conclude that the global L_FIR–$L^{\prime }_{\rm X}$ relations in (U)LIRGs and DSFGs, as parameterized by $\log L_{\rm FIR} = \alpha \log L^{\prime }_{\rm X} + \beta$, are linear for heavy rotor lines, and our CO line data set up to J = 6–5, at which point the relations become increasingly sublinear for higher J. Moreover, the normalization factor β shows a similar behavior by being nearly constant β ∼ 2 up to J = 5–4/6–5, but then starting to increase systematically with increasing J, reaching β ∼ 8 for J = 13–12.

The simplest scenario outlined in Section 4.2 seems to work both for the low-J CO and the HCN and CS lines with much higher critical densities. Given that our sample solely consists of (U)LIRGs (i.e., L_{IR[8-1000 μm]} ⩾ 10¹¹ L_☉), for which the dense gas fraction (i.e., $f_{\rm dense} = L^{\prime }_{\rm HCN_{\rm low-J}}/L^{\prime }_{\rm CO_{\rm low-J}}$) is nearly constant (see discussion in Section 6), linear slopes are to be expected for low-J FIR–CO relations. For other samples in the literature (e.g., Bayet et al. 2010) that reach lower IR luminosities (∼10¹⁰ L_☉) and thus span a wider range in L_IR (over which f_dense changes appreciably) the super-linear slopes of their low-J FIR–CO relations seen in Figure 3 are also expected. In this simple picture, neither the occasional super-linear nor the linear slope of the $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ low-J relations carry any profound ISM physics other than that more dense gas mass corresponds to proportionally higher SFRs (as also suggested by Gao & Solomon 2004b and Wu et al. 2005).

For dense gas tracer lines, this shows itself directly with $L_{\rm IR}\hbox{--}L^{\prime }_{\rm X}$ relations that always have linear slopes. Here there is actually more ISM physics to be found in exploring what sets the value of the normalization parameter β, rather than the slope of $L_{\rm IR}\hbox{--}L^{\prime }_{\rm X}$ relations.

The low- and high-z (U)LIRGs studied here are highly dust-obscured galaxies, and the radiation pressure exerted by the strong absorption and scattering of FUV light by dust grains could be an important feedback mechanism, possibly setting the value of the normalization (and ultimately regulating the SF). The maximum attainable L_IR/M_dense ratio of a star forming region before radiation pressure halts higher accretion rates is ultimately set by the Eddington limit, giving, L_IR/M_dense ∼ 500 L_☉ M_☉⁻¹ (Scoville & Polletta 2001). Andrews & Thompson (2011) expressed the expected $L_{\rm IR}\hbox{--}L^{\prime }_{\rm CO}$ and $L_{\rm IR}\hbox{--}L^{\prime }_{\rm HCN}$ relations in the case of Eddington-limited SFRs, and found that for CO luminosity tracing only the active star forming gas the maximal possible luminosity is given by $L_{\rm Edd}=4\pi G c \kappa ^{-1} X_{\rm CO} L^{\prime }_{\rm CO}$, where κ is the Rosseland-mean opacity, and X_CO is the $L^{\prime }_{\rm CO}$-to-$M_{\rm H_2}$ conversion factor. A similar expression holds for HCN, albeit with different κ and X values (see Andrews & Thompson 2011 for details). Although the exact normalization of this relation for each molecular line depends on poorly constrained quantities like κ and X, the Eddington limit set by the strong FUV/optical radiation from embedded SF sites acting on the accreted dust and dense gas can naturally provide the normalization of the observed $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ and $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm HCN}$ relations. In fact, adopting κ = 5–30 cm² g⁻¹ and X_CO = 0.8 K km s⁻¹ pc², which are perfectly reasonable values for (U)LIRGs (Thompson 2005; Solomon et al. 1997), we find β = log (4πGcκ⁻¹X_CO) = 2.5–3.3. We note that the high κ-value (30 cm² g⁻¹), which corresponds to the three-fold increase in the dust-to-gas mass ratio for the Rosseland-mean opacity (see Andrews & Thompson 2011 for details) that might be expected in (U)LIRGs, yields β = 2.5, which is close to the observed normalization values obtained (β ≃ 2) for the low-J CO lines in Section 3 (Table 3; see also Figure 4).

Given that the Eddington limit is ultimately set within individual SF sites embedded deep inside molecular clouds, it will operate on all galaxies, not just (U)LIRGs. For ordinary star forming spirals, the global β normalization value of the $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ relations for low-J CO lines will be lower than its (Eddington limit) set value by a factor approximately equal to the logarithm of its dense gas fraction, that is, log (f_dense) = log (M_dense/M_tot). By the same token, the offset in the $\log (L_{\rm FIR}) \hbox{--} \log (L^{\prime }_{\rm CO})$ plane between two populations with significantly different dense gas fractions (e.g., f_{dense, 1} and f_{dense, 2}) can be shown to be Δβ ∼ log (f_{dense, 2}/f_{dense, 1}). Thus, an increasing f_dense(L_IR) function can cause the super-linear FIR−CO(low-J) relations seen in some galaxy samples (which in reality are a varying β(L_IR) rather than a superlinear α).

Local (U)LIRGs and high-z DSFGs, on the other hand, form stars closer to the Eddington limit on a global scale (Andrews & Thompson 2011). Thus, the linear FIR–CO (low-J) and FIR−HCN/CS relations observed for local (U)LIRGs and high-z DSFGs is consistent with the notion that radiation pressure is an important physical mechanism that underlies the observed star formation laws in highly dust-obscured galaxies. In effect, the extreme merger/starbursts that dominate the (U)LIRGs and high-z DSFG populations resemble dramatically scaled-up versions of dense gas cores hosting SF deep inside Giant Molecular Clouds (GMCs), with the balance between radiation pressure and self-gravity setting their equilibrium during the IR–luminous phase.

In this framework, the failure of the available theoretical models to account for the observed $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm CO}$ relations of low-J CO and heavy rotor molecular lines might be attributed to the role that radiation pressure feedback plays in ultimately determining such relations. This has not been taken into account in all current theoretical considerations that either seek to explain the S-K relation in galaxies (e.g., Elmegreen 2002), or use an (S-K) relation of ρ_SFR∝(ρ_gas)^1.5 to determine the emergent $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm line}$ relations for molecular gas (Krumholz & Thompson 2007; Narayanan et al. 2008). These use self-gravity and the associated timescale of free-fall time, along with models on how the SF efficiency (gas mass fraction converted into stars per freefall time) varies per phase in turbulent gas as the main ingredients toward a complete understanding of SF, S-K relations, and the proxy $L_{\rm FIR}\hbox{--}L^{\prime }_{\rm X}$ relations. A nongravitational force like that exerted by radiation pressure on accreted gas and dust near SF sites can greatly modify such a picture by reducing or eliminating the dependence on the free fall time, especially for the high-density gas (which, presumably, is the one closest to active SF sites). Alternatively, it has been suggested that in lower luminosity systems the star formation may be regulated by feedback-driven turbulence (kinetic momentum feedback) rather than by radiation pressure (Ostriker & Shetty 2011; Shetty & Ostriker 2012; Kim et al. 2013). Assuming a continuum optical depth at FIR wavelengths (τ_FIR) of the order of unity for our sample of starburst/merger (U)LIRGs and typical dust temperatures of ∼50 K, we can make a rough estimate of the expected radiation pressure, namely $P_{\rm rad} \sim \tau _{\rm FIR} \sigma T_{\rm d}^4/c \sim 1.2\times 10^{-8}\,{\rm erg\,cm^{-3}}$ (where σ is Stefan–Boltzmann's constant and c the speed of light). This is comparable to the turbulent pressure $P_{\rm turb}\sim \rho \sigma _{\rm v}^2/3 \sim 1.4\times 10^{-8}\,{\rm erg\,cm^{-3}}$, obtained assuming a turbulent velocity dispersion of σ_v ∼ 5 km s⁻¹ and an average gas mass density of $\rho \sim 2\mu n_{\rm H_2}$ corresponding to $n_{\rm H_2}\sim 10^5{\rm cm^{-3}}$. Both of these greatly exceed the expected thermal pressure—$P_{\rm th}\sim n_{\rm H_{\rm 2}} k_{\rm B} T_{\rm k}\sim 1.4\times 10^{-9}\,{\rm erg\,cm^{-3}}$ (for $n_{\rm H_2}\sim 10^5\,{\rm cm^{-3}}$ and T_k ∼ 100 K)—thus highlighting the point made previously that a complete physical model of the S-K relations has to incorporate the effects of radiation pressure and/or turbulence.

A more direct indication of significant amounts of warm and dense gas in our (U)LIRG-dominated sample, and to what extent its thermal state is likely to be maintained by the SFR-powered average FUV radiation fields, is provided by the CO SLEDs. In Figure 6 we show the FIR- and CO(1–0)-normalized CO SLEDs (top and middle panels, respectively), as well as the raw CO luminosities (bottom panel). The first version allows for an assessment of the CO SLEDs for the full samples (not all of our sources have CO J = 1–0 measurements), and shows the cooling power of the CO lines with respect to the continuum. The CO J = 1–0 normalized representation of the CO SLEDs makes for a direct comparison with observed CO line ratios in the literature, and is furthermore what is usually used to constrain the excitation conditions of the gas.

Zoom In Zoom Out Reset image size

A detailed analysis of the CO SLEDs, in conjunction with the multi-J HCN, CS, and HCO⁺ line data setsavailable for many of the (U)LIRGs in Figure 6, is needed for a full understanding of the heating and cooling mechanisms of the molecular gas and for quantifying the relative mass-fractions of the gas phases. Nevertheless, the marked contrast between their CO SLEDs and that of the Milky Way disk (where most of bulk of the molecular gas is warmed by photoelectric heating induced by the ambient FUV radiation field), already indicates the presence of a different heating source. While intense X-ray radiation fields (1–5 keV) generated by AGN can penetrate and heat gas up to ≳ 100 K at column densities of 10²²–10²⁴ cm⁻² (and without dissociating all the CO), it is unlikely to be the case here because great care has been taken in removing AGNs from our sample. The integrated power emitted in the CO J = 7–6 to J = 13–12 transitions for all the (U)LIRGs in our sample constitutes the bulk (about 60%) of the total energy output of all the CO lines. Exploring the effects of different CR and mechanical heating rates on the thermal structure of clouds, Meijerink et al. (2011) found that even for extreme CR fluxes (∼10²–10⁴ × the Milky Way value) it is difficult to maintain temperatures ≳ 100 K, and the effect on the high-J CO lines appears to be minor. Mechanical heating, such as supernova driven turbulence and shocks, however, was found to heat the gas more efficiently, and we favor this as the most likely explanation for hot gas and the boosted high-J CO lines observed in our (U)LIRGs.

Highly excited CO SLEDs have been found for the merger/starburst NGC 6240, where a recent analysis found FUV photons (and the resulting photoelectric heating) to be inadequate as the main heating source for the high temperatures of its dense gas (Meijerink et al. 2013; Papadopoulos et al., 2014), and for Mrk 231 where X-rays from the AGN are thought to heat the dense molecular gas reservoir (van der Werf et al. 2010; although Mrk 231 has recently been shown to be much less X-ray luminous than previously thought; see Teng et al. 2014). While such CO SLEDs, and the need for alternative heating mechanisms than FUV–photons to explain them, might be linked to the unusually high CO line-to-continuum ratios of both of these two sources, a similar conclusion was reached for M 82 and NGC 253 based on analyses of their full CO SLEDs (Panuzzo et al. 2010; Rosenberg et al. 2014a). Our work is the first to demonstrate such highly excited high-J CO SLEDs as a near generic characteristic of merger/starbursts (the galaxies that dominate the sample shown in Figure 6).

As mentioned in Section 2, incorrect FIR–CO relations may be inferred if the FIR and CO measurements cover different regions within galaxies. This can be a serious problem for local extended sources where single-dish CO beams can be smaller than the extent of the IR emission (see discussion in Zhang et al. 2014). We are confident that this is not an issue for our HerCULES sample, where all SPIRE-FTS CO line fluxes were scaled to a common 42'' angular resolution, which is beyond the extent of the IR emission in these sources (as traced by LABOCA 870 μm maps). However, if this was not the case, and the SPIRE-FTS measurements did not capture all the CO emission, it would imply that the derived FIR–CO slopes are biased high (since the CO luminosity will be underestimated relative to the total FIR luminosity, see Figure 2). In short, the sublinear FIR–CO slopes at high-J transition found here are robust against the (unlikely) possibility that some (small) fraction of the CO emission is unaccounted for.

The effects that the presence of strong AGNs would have on the SF relations are two-sided. On the one hand, it could lead to an overestimate of the IR luminosity attributed to star formation, thus biasing the FIR–CO slopes high. On the other hand, AGN-dominated environments, where penetrating X-rays may be dominating the gas heating, tend to have boosted high-J CO lines compared to star forming regions (e.g., Meijerink et al. 2007). If the AGN was deeply buried it would not be easily detectable in X-rays and would be optically thick in the IR, possibly down to mm wavelengths. The effect could be what is observed here—a change of slope in the correlation and the high-J CO lines reflecting a hot, deeply embedded AGN. As mentioned in Section 2, we used L_{FIR[50-300 μm]} instead of L_{IR[8-1000 μm]} in the FIR–CO relations in order to minimize the effects of AGN. More importantly, we only included local (U)LIRGs for which the bolometric AGN contribution was deemed to be <20%, as estimated from several MIR diagnostics (Section 2). In the case of the high-z sample, systems that were obviously AGN-dominated were immediately discarded (Section 2). Furthermore, we note that deep X-ray observations and MIR spectroscopy of millimeter and submillimeter selected DSFGs (which show no obvious signs of harboring an AGN) have shown that any AGN that might be present typically contribute ≲ 20 % to the total IR luminosity (Alexander et al. 2005; Mendez-Delmestre et al. 2007). Thus, we feel confident that neither the FIR nor the CO luminosities are biased high due to AGNs, and therefore our findings are not systematically affected by AGNs.

Galaxies that are gravitationally lensed are prone to differential magnification, an effect in which regions within a galaxy are magnified by different amounts due to variations in their location within the galaxy, and/or spatial extent (Blain 1999). This can significantly skew the observed relative contributions from hot versus cold dust to the IR luminosity, as well as low- versus high-J CO line luminosity ratios (Serjeant 2012). Furthermore, a flux-limited sample of strongly lensed sources will tend to preferentially select compact sources (Hezaveh et al. 2012), which may be more likely to have extreme CO excitation conditions. From Figure 6(a), however, there is nothing to suggest that the lensed DSFGs have markedly different $L^{\prime }_{{\rm CO}_{J,J-1}}/L_{\rm FIR}$ values than the nonlensed and local (U)LIRGs. In Figure 6(b), however, we do see a few lensed DSFGs that have markedly higher $L^{\prime }_{{\rm CO}_{J,J-1}}/L^{\prime }_{{\rm CO_{1,0}}}$ ratios at high-J than the other samples. This is exactly what we would expect to see if these sources were differentially lensed (and the high-J lines tracing more compact regions than the J = 1–0 line), or if the lensing preferentially selects compact sources (which would tend to have more extreme excitation conditions). A strong argument against our analysis being affected by differential magnification effects is the fact that when fitting the FIR-CO relations without the lensed DSFGs we obtain slopes nearly identical to the ones given in Table 3, and fully consistent within the errors. Thus, we conclude that our findings are unlikely to be affected in any significant way by differential magnification effects.