Molecular Gas Heating, Star Formation Rate Relations, and AGN Feedback in Infrared-Luminous Galaxy Mergers

We examine the origin of molecular gas heating in a sample of 42 infrared-luminous galaxies at z < 0.3 by combining two sets of archival data: first, integrated CO line luminosities in the 1–0 and 5–4 through 13–12 transitions; second, results from radiative transfer modelling that decompose their bolometric emission into starburst, AGN, and host galaxy components. We find that the CO 1–0 and 5–4 through 9–8 lines primarily arise via radiative heating in the starburst and the host galaxy. In contrast, the CO 10–9 through 13–12 lines may arise primarily in the starburst and AGN, with an increasing contribution from mechanical heating and shocks. For the sample as a whole, we find no evidence that AGN luminosity affects the heating of molecular gas by star formation. However, for starbursts with low initial optical depths, a more luminous AGN may reduce the efficiency of starburst heating of the CO 5–4 and above lines, consistent with negative AGN feedback.

The heating of molecular gas in infrared-luminous galaxies is thus an important diagnostic of both the conversion of molecular gas to stellar and SMBH mass, and of the impact of starburst and AGN activity on the interstellar medium (ISM) of their host galaxies [50][51][52][53][54]. While the bulk of the molecular gas in galaxies is in the form of 'cold' H 2 , this molecule's lack of a dipole means that cold H 2 is challenging to observe directly 1 . Instead, the molecular gas is usually traced via observations of rotational transitions of carbon monoxide ( 12 CO, CO hereafter), since CO has both a strong dipole moment and relatively low-lying energy transitions. The ground state transition (J = 1-0) and transitions up to about J = 4-3 of CO are good tracers of the total molecular gas reservoir in galaxies [56][57][58][59][60][61][62][63]. The higher J transitions trace the warmer, denser molecular gas associated with starbursts and AGN (e.g., [64][65][66][67][68]). As a result, observations of CO emission in high-redshift galaxies have produced important insights into galaxy assembly processes [39,69,70].
Significant uncertainties remain over how different transitions of CO are produced. In principle, CO emission can arise in star-forming regions, AGN, or in the ISM heated by passively evolving stellar populations; One way to address this is to model physical conditions in the line-emitting gas using multiple transitions of CO simultaneously [71,72]. Another approach is to compare CO line luminosities to the total infrared luminosity of the galaxy [65]. However, ambiguities in the power source behind the infrared emission persist in this approach.
In this paper, we take an alternative approach to determining the origin of CO heating in infrared-luminous galaxies. We take a sample of 42 gas-rich mergers with infrared luminosities in excess of 10 12 L at z < 0.3, for which we have completed radiative transfer modelling which constrains their starburst, AGN, and host galaxy luminosities. We then assemble archival observations of the 1-0 and 5-4 through 13-12 transitions of CO for this sample, and compare them to the component luminosities from the radiative transfer modelling to infer the origin of the CO gas heating as a function of rotational quantum number. Throughout, we assume H 0 = 70 km s −1 Mpc −1 , Ω = 1, and Ω Λ = 0.7. We convert all literature data to this cosmology where necessary.

Sample Selection
The sample selection is described in full in [2,73], and is summarized here. The sample includes ultraluminous infrared galaxies (ULIRGs, systems with rest-frame 1-1000 µm luminosities in excess of 10 12 L ) with IRAS 60 µm fluxes greater than ∼2Jy that were assembled for the HERschel ULIRG Reference Survey [16,68,74,75]. This sample is representative of ULIRGs in the low-redshift Universe as it comprises nearly all known systems with L IR < 10 12 L at z < 0.27. The sample spans approximately one order of magnitude in total infrared luminosity, includes merger stages from early to late stage, and all four primary optical spectral classes (HII, LINER, Sy2, Sy1). The host galaxy properties are typical of massive, gas-rich galaxies in the local Universe. These data and diagnostic plots are given in [2,73].

Molecular Gas Data
We use the CO 1-0 transition, which traces cold molecular gas, and the CO 5-4 through 13-12 transitions, which trace hotter molecular gas. We explored the possibility of using the 2-1 through 4-3 lines, and in lines above 13-12. There are, however, insufficient observations of our sample in these lines to enable a statistical analysis. For CO 1-0 we use the compilation presented in Farrah et al. [2]. These data are mostly taken from Kamenetzky et al. [76], and/or a wider set of archival observations [56,67,[77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94]. For the 5-4 through 13-12 lines, we assemble data from Pearson et al. [68], who present a homogeneous analysis of these lines as observed by the SPIRE instrument onboard Herschel. Data for a few objects are also taken from Kamenetzky et al. [76]. In all cases, the CO line observations are spatially unresolved. Spatially resolved, interferometric observations of about a quarter of our sample exist, but we use spatially unresolved observations in all cases to give a consistent analysis.

Infrared Emission Data
We adopt the results from Efstathiou et al. [73]. This study decomposes the emission from our sample into contributions from a starburst, an AGN, and the host galaxy. The results from this study are consistent with previous SED modelling works, and provide derived physical parameters for each component. The complete set of parameters for each model, and their adopted ranges, are given in Table 2 of [73].
The starburst model [95] assumes an exponentially decaying burst of star formation. The model includes a range of possible initial optical depths of the star-forming clouds, efolding times of the starburst, and a range in starburst ages, thus tracking multiple possible evolutionary paths for the molecular clouds hosting the star formation.
The AGN model [96][97][98] assumes a smooth, tapered disk in which the half-opening angle of the torus, the inclination angle of the torus relative to the observer, the ratio of inner to outer radius, and equatorial optical depth, can all vary. The AGN model also includes luminosities derived by integrating the observed emission over 4π steradians, as well as luminosities that include a correction for the anisotropic structure of the obscurer.
The host model assumes a Sérsic profile with n = 4. The stars and dust are mixed, rather than the extinction applied as a foreground screen. The three parameters of the host model are the e-folding time of the star formation, the optical depth of the host, and the intensity of starlight.
For our study, the key extracted physical quantities from the radiative transfer models are: • L o tot : the total rest-frame infrared (1-1000 µm) luminosity, uncorrected for anisotropic AGN emission. • L Sb ,Ṁ Sb : the rest-frame infrared luminosity of the starburst, and the SFR in the starburst averaged over the starbursts age, excluding the host. • L h ,Ṁ h : the rest-frame infrared luminosity and SFR of the host galaxy, excluding the starburst. • L o AGN , L c AGN : Observed and anisotropy-corrected rest-frame infrared luminosity of the AGN. • τ V : the initial optical depth of giant molecular clouds in the starburst, with an allowed range of 50 < τ V < 250.
A description of these quantities is given in Efstathiou et al. [73].

Correlations with Total Infrared Luminosity
We first examine the relations between the CO luminosities and the total observed infrared luminosities (L o tot ). Using the Kendall-τ test, we find significant correlations (τ > 0.38, p < 0.01) between all CO line luminosities and L o tot . The p-values are plotted as the black line in Figure 1. We next derive relations between CO luminosity and L o tot , of the form: We used the Orthogonal Distance Regression algorithm [99] as implemented within the SciPy Python library using a random seed of 2001. This yields: log(L 10-9 ) = 1.39 ± 0.27 log(L o Tot ) − 9.50 ± 3.30 (8) log(L 11-10 ) = 0.25 ± 0.15 log(L o Tot ) + 4.31 ± 1.86 (9) log(L 12-11 ) = 0.04 ± 0.15 log(L o Tot ) + 6.79 ± 1.88 (10) log(L 13-12 ) = 0.22 ± 0.17 log(L o Tot ) + 4.56 ± 2.01 (11) (see also, e.g., Greve et al. [65]). The slopes as a function of J up are plotted in Figure 2. The slopes are consistent with a linear relation between L CO and L o tot for J = 10-9 and below, but a sub-linear, possibly flat relation for J = 11-10 and above. Some example fits are presented in Figure 3. Host Infrared Luminosity (L ) Observed AGN Infrared Luminosity (L )  (Table 1), normalized by the relevant mean CO luminosity L CO , and total, starburst, host, and observed AGN infrared luminosities (Sections 3.1 and 3.2). The best-fit relations are plotted, if there is evidence for a correlation. No strong dependence on optical spectral type is observed.

Component Luminosity Correlations
We next evaluate the Kendall-τ coefficients for each CO luminosity against L Sb , L o AGN , L c AGN , and L h (Figure 2). All of the CO luminosities correlate with L Sb . For L h , L 1-0 and L 5-4 through L 9-8 show correlations, though they may be weaker than those with L Sb . However, L 10-9 and above do not show correlations with L h . The starburst and host SFRs exhibit similar behavior to L Sb and L h . For L o AGN and L c AGN the L 10-9 and above luminosities show a correlation, but L 9-8 and below do not (see also, e.g., Esposito et al. [100]). The L o AGN and L c AGN luminosities show identical behaviour in these correlations, consistent with the inclination of the AGN obscurer relative to the line of sight not playing a dominant role in determining the observed CO line luminosities.

Multi-Component Fits
To facilitate comparisons with galaxy evolution models, we present two-component fits of the form L CO = αL A + βL B + γ, where the luminosities L A and L B are chosen based on the results of the Kendall-τ tests; starburst and host of CO 9-8 and below, and starburst and AGN for CO 10-9 and above. This yields: These fits are consistent with the log-linear fits, though they suggest the AGN may be the dominant contributor to the CO 11-10 and above emission. The fits are given in Figure 5. We do not present fits of the form: L CO = αL Sb + βL AGN + γL h as they proved challenging to constrain with the number of objects in our sample.

Star Formation Rate-Molecular Gas Scaling Relations
As an alternative diagnostic of the relation between the cold ISM and star formation (see, e.g., [103,104] for reviews) in our sample, we consider the relation between SFR and cold molecular hydrogen mass. To do so, we use the cold H 2 masses for our sample as summarized in Farrah et al. [2]. Since the total SFR is dominated in most cases by the starburst, we first fit a relation of the form in Equation (1). We obtain: log(M H 2 ) = 1.11 ± 0.19 log(Ṁ Sb ) + 6.92 ± 0.49 (52) This is consistent with the proposed universal scaling relation of Lada et al. [105]. Instead using a two-component model yields: We note though that Equation (53) does not give substantially reduced scatter relative to Equation (52). These results are summarized in Figure 6.

Starburst Efficiency
We next examine whether or not AGN luminosity can affect the heating of CO by the starburst. To do so, we consider the L Sb /L CO ratio as a function of L c AGN , via the relation For the whole sample, we find a mildly negative relation between L Sb /L CO and L c AGN in most CO lines (Figure 7). Combined with the results from Section 3.1, this can be interpreted as a significant contribution to heating the CO by the AGN in the CO 9-8 and above lines. We do not find that the AGN affects the starburst heating of the molecular gas, at least for the sample as a whole (see also, e.g., Shangguan et al. [106], Zhuang et al. [107], Valentino et al. [108] but also McKinney et al. [109]). We note though that the factors that affect fueling of star formation in infrared-luminous mergers are subtle, and may vary by location within individual sources (e.g., [110]). Evidence for a relation between L Sb /L CO and L c AGN emerges when dividing the sample into two sub-samples at τ V ∼ 170 and fitting relations of the form in Equation (54). The τ V > 170 sub-sample yields the same or slightly more negative slopes than the full sample-that is, similar trends in CO luminosity per unit starburst luminosity, as AGN luminosity increases. Conversely, the τ V < 170 sub-sample shows marginally positive slopes, suggesting less CO luminosity per unit starburst luminosity, as AGN luminosity increases. However, though the formal uncertainties imply this difference is highly sig-nificant, simulations that divide the sample in two equal sub-samples at random give an estimated confidence of 3σ.
To interpret this result, we assume that the direct AGN contribution to producing CO emission is insignificant in all lines below CO 10-9, and sub-dominant for CO 10-9 and above. We also assume that the CO emission is not significantly absorbed. Our result is then consistent with more luminous AGN reducing the ability of star formation to produce warm and hot CO in starbursts with lower initial optical depths (see also, e.g., Song et al. [111]). This effect can be interpreted as negative AGN feedback [89,[112][113][114][115][116][117][118], though the effect is subtle, and unlikely to be a viable channel to globally quench star formation in infrared-luminous galaxies. We speculate that the origin of the effect is that higher optical depth starbursts have higher ISM densities, which are harder to disrupt by a luminous AGN (see also, e.g., [119,120]).
Other explanations seem less plausible. We do not see a dependence on the depth of the 9.8 µm silicate feature, or on the ratio of polycyclic aromatic hydrocarbon equivalent widths, suggesting that this effect does not arise in compact obscured nuclei [121]. More extincted starbursts could just be less luminous, but we see no such trend between L Sb and τ V . Alternatively, more extincted starbursts could be intrinsically more efficient at heating molecular gas. If this were true though, then no trends with AGN luminosity would be expected. Linear fit slope α 10 11 10 12 10 13 Anisotropy-corrected AGN Infrared Luminosity (L )

Conclusions
We have studied the origin of CO emission in a sample of 42 infrared-luminous galaxy mergers at z < 0.27 with infrared luminosities in excess of 10 12 L . To do so, we take results from a recent radiative transfer modelling study of this sample, which decomposes their infrared luminosities into contributions from star formation, AGN activity, and the host galaxy. We then compare these component luminosities to archival measures of the integrated CO luminosities in the 1-0, and 5-4 through 13-12 transitions. Our conclusions are: 1.
The luminosities of the CO 1-0 and 5-4 through 13-12 lines are consistent with an origin in one or both of ongoing star formation, and heating by starlight in the host galaxy. We do not find evidence for a significant contribution to these lines from AGN activity.

2.
We present log-linear relations between the CO luminosities in the 1-0 and 5-4 through 13-12 transitions, and with the starburst and host galaxy luminosities. The slopes of these relations suggest an origin of the CO emission in gas heating.

3.
The luminosities of the CO 10-9 through 13-12 lines are consistent with an origin in one or both of the starburst and the AGN. We do not find evidence for a significant contribution from the host galaxy. The slopes of the relations between these CO line luminosities and the starburst and AGN luminosities may be sub-linear, which is consistent with a contribution to these lines from mechanical heating and shocks.

4.
We also present two-component linear model fits between each CO line luminosity and infrared component luminosities. For CO 9-8 and below we fit against starburst and host luminosity. For CO 10-9 and above we fit against starburst and AGN luminosity. These fits are consistent with the log-linear model fits, and give straightforward conversions that can be used in galaxy simulations.

5.
For the sample as a whole, we find no evidence for a relation between AGN luminosity and the L CO /L Sb , in any CO line. This suggests that a more luminous AGN does not reduce the efficiency by which the starburst heats the molecular gas. There is, however, evidence for a dependence on starburst initial optical depth (τ V ). Starbursts with low τ V may heat the CO 5-4 transitions and above less efficiently. This is consistent with mild negative AGN feedback. Funding: This research received no external funding.

Data Availability Statement:
All the data used in this paper are publicly available within this paper, or within the references cited.