© The Authors 2023

1. Introduction

In the local (z ≲ 0.2) Universe, ultra-luminous infrared galaxies (ULIRGs, where L_IR(8 − 1000 μm)≥10¹² L_⊙) and luminous infrared galaxies (LIRGs, where L_IR(8  −  1000 μm)≥10¹¹ L_⊙) pinpoint gas-rich galaxy mergers undergoing intense starbursts (SBs) and supermassive black hole accretion (Sanders & Mirabel 1996; Genzel et al. 1998; Lonsdale et al. 2006; Pérez-Torres et al. 2021; U 2022). These processes together deeply modify the physical and dynamical properties of the interstellar medium (ISM), likely leading to a permanent morphological transformation and quenching (e.g., Hopkins et al. 2008).

(Sub)millimeter interferometric observations of CO lines (e.g., Downes & Solomon 1998; Wilson et al. 2008; Ueda et al. 2014) and dense H₂ gas tracers such as HCN and HCO⁺ (Aalto et al. 2012; Imanishi & Nakanishi 2014; Imanishi et al. 2019; Ledger et al. 2021) show that the extreme star formation rates (SFRs) of (U)LIRGs are fueled by massive (M_H2 > 10⁹ M_⊙), dense H₂ gas reservoirs, characterized by a high surface brightness in the central kiloparsec-scale region. However, feedback mechanisms and tidal forces can disperse ISM material outside of the nuclear regions (Springel et al. 2005; Narayanan et al. 2006, 2008; Duc & Renaud 2013). Hence, we can expect a portion of the ISM of (U)LIRGs to reside in diffuse, low-surface-brightness structures, possibly missed by high-resolution interferometric observations.

Galactic outflows have been ubiquitously observed in (U)LIRGs for decades, in the ionized (Westmoquette et al. 2012; Arribas et al. 2014) and atomic (Rupke et al. 2005; Martin 2005; Cazzoli et al. 2014) gas phases, as expected in sources affected by strong radiative feedback from SBs and active galactic nuclei (AGNs; e.g., Costa et al. 2018; Biernacki & Teyssier 2018). More recent is the discovery that the outflows of (U)LIRGs can embed large amounts of molecular gas, traveling at speeds of up to v ∼ 1000 km s⁻¹. Such molecular outflows have been detected unambiguously by Herschel, via observations of P-Cygni profiles and/or blueshifted absorption components of far-infrared (FIR) OH, H₂O, and OH⁺ transitions (Fischer et al. 2010; Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Stone et al. 2016; González-Alfonso et al. 2017, 2018), as well as through the investigation of broad and/or high-velocity components of CO (Feruglio et al. 2010, 2015; Cicone et al. 2012, 2014; Pereira-Santaella et al. 2018; Lutz et al. 2020; Fluetsch et al. 2019; Lamperti et al. 2022), HCN, HCO⁺ (Aalto et al. 2012, 2015; Barcos-Muñoz et al. 2018), and CN (Cicone et al. 2020) emission lines. These sensitive observations have also shown that the molecular ISM of (U)LIRGs, and especially the low-surface-brightness outflow components (e.g., Feruglio et al. 2013; Cicone et al. 2018; Herrera-Camus et al. 2020), can extend to several kiloparsecs, up to the edge of the field of view of single-pointing interferometric data.

Obtaining robust H₂ mass measurements of the total ISM reservoirs as well as of the gas embedded in outflows is crucial for understanding the impact of gas-rich galaxy mergers – and of the collateral powerful SB and AGN feedback mechanisms – on galaxy evolution. Low-J CO lines such as CO(1–0) and CO(2–1) can be used to estimate H₂ masses through a CO-to-H₂ conversion factor (hereafter α_CO); it is, however, highly dependent on the physical state of the gas. The α_CO parameter can vary by up to a factor of ∼10 in different ISM environments, depending on the CO optical depth, on the metallicity of the medium, and on the exposure to far-UV radiation and cosmic rays (CRs) that can destroy CO more than H₂ (see, e.g., Bisbas et al. 2015; Glover et al. 2015; Offner et al. 2014). For the molecular ISM of disk galaxies, the conventional α_CO factor is 4.3 M_⊙ (K km s⁻¹ pc⁻²)⁻¹ (Strong & Mattox 1996; Abdo et al. 2010; Bolatto et al. 2013), while for more perturbed galaxies, such as gas-rich mergers and SBs, a lower α_CO factor of ∼0.8 − 1.0 M_⊙ (K km s⁻¹ pc⁻²)⁻¹ is often preferred (Downes & Solomon 1998).

Combining multiple molecular transitions can help constrain the physical properties of molecular gas and so derive a better estimate of the α_CO factor. A particularly valuable H₂ gas tracer is the forbidden ³P₁−³P₀ fine structure line of atomic carbon, hereafter [CI](1–0), which has an excitation temperature of T_ex = 23.6 K and a critical density similar to that of CO(1–0) (i.e., n_crit, [CI] ∼ 1000 cm⁻³). The [CI](1–0) line is optically thin and has a simple three-level partition function, which makes it easier to interpret than low-J CO lines (see the discussion in Papadopoulos et al. 2022). Early models of photodissociation regions (PDRs) predicted [CI] to exist in a thin transition layer between the central region of molecular clouds (molecular gas – CO) and its envelope (ionized gas – [CII]; for the standard PDR view, see Tielens & Hollenbach 1985). However, observations have shown [CI] to be well mixed with CO widely throughout the cloud (Valentino et al. 2018; Saito et al. 2020; Papadopoulos et al. 2022), suggesting that both species might trace the bulk of molecular gas mass (Ojha et al. 2001; Papadopoulos et al. 2004; Kramer et al. 2008; Salak et al. 2019; Izumi et al. 2020; Jiao et al. 2017). Moreover, theoretical models show that CO (but not H₂) may be destroyed in environments dominated by CRs, shocks, or intense radiation fields, leaving behind CO-dark or CO-poor reservoirs; this provides support to the idea of [CI] being an alternative H₂ gas mass tracer.

Observing the atomic carbon emission from local galaxies requires a sensitive submillimeter telescope located at a very high and dry site. Indeed, the [CI](1–0) transition, at a rest frequency of ${\nu }_{[\mathrm{CI}](1-0)}^{\mathrm{rest}}=492.161$ GHz (609.135 μm), in the absence of a significant redshift, can be observed only if the atmospheric opacity is low (precipitable water vapor PWV < 1 mm). For this reason, sensitive observations of [CI](1–0) in the local Universe are still very sparse, even for bright (U)LIRGs. Cicone et al. (2018) used the Atacama Large Millimeter/sub-millimeter Array (ALMA) and the Morita Array (also known as the Atacama Compact Array, ACA) to obtain high S/N [CI](1–0) observations of NGC 6240. These data, combined with archival CO(1–0) and CO(2–1) observations, were used to study the $r_{21}\equiv L_{\mathrm{CO}(2-1)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }$ and $r_{\mathrm{CICO}}\equiv L_{[\mathrm{CI}](1-0)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }$ line ratios (the latter can be used to estimate α_CO) in the massive molecular outflow of NGC 6240. From a spatially resolved analysis, Cicone et al. (2018) find that the outflowing ISM in NGC 6240 is robustly characterized by an α_CO that is lower than the non-outflowing H₂ medium, and that r₂₁ is higher for high-σ_v outflow components, especially at large distances from the nuclei. The Cicone et al. (2018) analysis suggests that: (i) despite its obvious limitations, a multicomponent decomposition of galaxy-integrated spectra, performed simultaneously to multiple transitions, which enables an investigation of line ratios separately for spectral components with different widths and central velocities, can deliver results that are consistent with a proper spatial decomposition of the ISM for the outflowing and disk material; and (ii) the outflowing H₂ gas may be characterized by different physical properties from the non-outflowing ISM of NGC 6240, and in particular by a lower CO optical depth and a higher CO excitation. These results have been obtained on a single, extreme source, and further statistics are required.

Our study builds upon these previous results and aims to expand the analysis of Cicone et al. (2018) to a sample of 36 local ULIRGs and 4 additional LIRGs with low-J CO (up to J = 3) and [CI](1–0) line observations. In designing our survey, we paid particular attention to capturing the total flux from these sources, including possible extended low-surface-brightness components that may be dominated by outflows and tidal tails and may be missed by high-resolution, low S/N interferometric data. The final survey contains proprietary and archival data from the Atacama Pathfinder EXperiment (APEX) telescope, ALMA, ACA, and the IRAM Plateau de Bure Interferometer (PdBI, now the NOrthern Extended Millimeter Array, NOEMA), which we have re-reduced and reanalyzed in a consistent and uniform way. Therefore, we can rely both on a consistent data analysis and on high-quality spectra, all taken with receivers whose large instantaneous intermediate frequency (IF) bandwidth can properly sample the extremely broad emission lines of (U)LIRGs.

This paper is organized as follows. In Sect. 2 we describe the sample selection. In Sect. 3 we describe the observing strategy, the observations, and the data reduction. In Sect. 4 we explain the methodology used for the spectral fitting and the data analysis. Our results are presented in Sect. 5, where they are also contextualized through a comparison with relevant results from the literature. A more general discussion is reported in Sect. 6. Finally, Sect. 7 summarizes the main results and presents the conclusions of our work. Throughout this work, we adopt a Λ cold dark matter cosmology, with H₀ = 67.8 km s⁻¹Mpc⁻¹, Ω_M = 0.307, and Ω_Λ = 0.693 (Planck Collaboration XVI 2014).

2. The sample

In absence of additional spatial information, the only unambiguous method for assessing the presence of molecular outflows trough spectroscopy, is the detection of P-Cygni profiles or blueshifted absorption components in molecular transitions, such as the OH119 μm line observed by Herschel (Fischer et al. 2010; Sturm et al. 2011). However, the absence of these features does not necessarily rule out the presence of outflows, as in the case of NGC 6240, studied in Cicone et al. (2018). In this source, OH is detected only in emission despite the presence of an extreme molecular outflow detected in multiple tracers. For this reason, molecular emission line observations can provide valuable and unique information on the presence and properties of galactic outflows, complementary to OH data.

The Herschel OH targets studied by Veilleux et al. (2013) and Spoon et al. (2013) represent the only conspicuous sample of (U)LIRGs that has uniform and unambiguous prior information about their molecular outflows, hence providing a robust comparison data set for our investigation based on emission lines. Moreover, the southern targets in this sample also have plenty of ancillary data from ALMA and APEX, allowing us to capitalize on public archives, which is a main focus of this work. For these reasons, the targets in our sample are selected from the Spoon et al. (2013) and Veilleux et al. (2013) samples, regardless of the detection of a molecular outflow in OH.

From the parent Herschel samples, we have included all sources with declination δ < 15 deg, except IRAS 12265+0219 and IRAS 00397-1312 for which we did not have any data available¹. NGC 6240 satisfies our selection criteria but is excluded from our work because it was the target of the pilot study by Cicone et al. (2018). Our sample includes 36 ULIRGs, whose physical properties such as redshifts, L_IR, SFRs, and AGN fractions (α_AGN ≡ L_AGN/L_bol) are reported in Table 1 with their corresponding references. The 4 additional LIRGs reported in Appendix A have been reduced and analyzed consistently with the rest of the sample, but they have been excluded from the main body of the paper to avoid biasing the relations given the low statistics for these low-L_IR sources. Hereafter, we use the term “(U)LIRGs” when referring to the entire sample, and “ULIRGs” when we exclude the 4 LIRGs.

The parent Herschel samples from which our targets were selected have a redshift upper limit of z < 0.2. As a result, our study investigates local (U)LIRGs with redshifts ranging from z = 0.00708 (IRAS F12243–0036) to z = 0.1935 (IRAS F05024–1941). The sample is by definition composed by high-L_IR galaxies with L_IR ranging from 10¹¹ L_⊙ to 10^12.8 L_⊙, covering luminosities within the (U)LIRG regime. The sources span a wide range in α_AGN values from 0.0 up to 0.92, with ∼50% of the sources having α_AGN ≥ 0.5. The SFRs are taken from the parent sample papers, with the exception of galaxies selected in Veilleux et al. (2013) for which there are no SFRs reported. In those cases, we followed the method by Sturm et al. (2011; also used in Spoon et al. 2013) to obtain the SFRs using: SFR = (1 − α_AGN)×10⁻¹⁰L_IR, so that all values are computed uniformly. The SFRs range from a couple of solar masses a year up to ∼300 M_⊙ yr⁻¹. We checked that our sample of ULIRGs is representative, in terms of physical properties (e.g., SFR, L_AGN, α_AGN) of the local ULIRG population by comparing it with the QUEST (Quasar/ULIRG Evolutionary Study) sample at z < 0.2 (see Veilleux et al. 2009a). Most previous works studying the molecular gas in local (U)LIRGs have included both LIRGs and ULIRGs. We compare some of our results with the works of Herrero-Illana et al. (2019) and Jiao et al. (2017), whose samples are however heavily dominated by LIRGs as opposed to ours. As the distinction between LIRGs and ULIRGs is based on an arbitrary L_IR cut, several galaxies that are officially LIRGs (such as NGC 6240) belong to the same population as the more infrared-luminous ULIRGs.

Table 1 lists the OH outflow velocity values for the sources (29 out of the whole sample of 40) that show an OH outflow detection according to Veilleux et al. (2013) and Spoon et al. (2013). We also report the OH equivalent widths for the whole sample.

Our sample, which focuses on southern (U)LIRGs targeted by previous Herschel OH observations, has naturally a large overlap with the APEX and ALMA/ACA public archives. Indeed, many of these sources have been observed in previous projects targeting different molecular tracers. In this work, we make the most out of such archives, focusing on the low-J CO transitions and [CI] atomic carbon line, and we complement them with our own new proprietary high-sensitivity single-dish observations with APEX. The observations and the data reduction process are described in detail in Sect. 3.

3. Observations

3.1. Observing strategy and data reduction

We want to study simultaneously the total integrated line emission from the three lowest-J transitions of CO and from [CI](1–0) in our sample of 40 (U)LIRGs. To do so, we combine proprietary and archival single-dish (APEX) and interferometric (ACA, ALMA, and IRAM PdBI) observations. The final, reduced spectra employed in our analysis are all shown in Figs. B.1–B.6.

In Table B.1 we report all the data sets considered in this paper with their respective project IDs. In those cases where multiple spectra are available for the same source and transition, we report at the top of the corresponding row in Table B.1 the data set that was used for our main analysis, followed by the one(s) that are not employed in the analysis. In such cases of duplication, we assign higher priority to data sets with the highest sensitivity to large-scale structures, namely: (1) APEX PI data, (2) APEX archival data, (3) ACA archival data, and, lastly, (4) ALMA archival data. In this way, we prioritize single-dish data that better trace the total flux, including possible extended emission that can be filtered out by interferometric observations. If, for a given transition, single-dish data exist but are of poor quality (i.e., have a low S/N or are affected by instrumental issues), we prefer the ACA or ALMA data when available for the same transition, after carefully checking that there is no significant flux loss. The ALMA/ACA archival data used here are not tailored to the aim of our study, and therefore, the angular resolutions and maximum recoverable scales (MAS) of the interferometric observations vary over a wide range, from a fraction of arcsec (in the most extended ALMA antenna configurations), up to ∼40 arcsec (in the ACA antenna configurations) for the angular resolution, and a few arcsec up to ∼90 arcsec for the MAS. We do, however, pay special attention to define an aperture for extracting the total flux that is equal to or smaller than the MAS of the observation. The duplicated spectra that were discarded from our main analysis, have been nevertheless reduced and are shown in Figs. C.1–C.3.

As summarized in Table B.1, we have CO(1–0) line spectral data for 22 galaxies (20 ULIRGs and 2 LIRGs), where 21 data sets are obtained from the ALMA/ACA data archive, and one from the IRAM PdBI (analyzed by Cicone et al. 2014). CO(2–1) line observations are available for all 40 sources of the sample (36 ULIRGs and 4 LIRGs); of these, 32 galaxies were observed with APEX through PI observations, 7 have archival APEX data, and 6 have ALMA/ACA archival data. As many as 31 galaxies have CO(3–2) line coverage (30 ULIRGs and only one LIRG): 18 sources have APEX PI data, 16 have APEX archival data and 12 have ALMA/ACA archival data. Lastly, we have APEX PI observations of the [CI](1–0) line for 17 galaxies of the sample (14 ULIRGs and 3 LIRGs), one of which resulted in a non-detection (the LIRG PG2130+099). For 7 of these sources there are also ACA archival [CI](1–0) data. Summarizing, we cover all three CO transitions for 18 galaxies (45% of the sample), two CO transitions for 17 galaxies (42% of the sample), and a single CO transition (CO(2–1)) for the remaining 5 galaxies (13% of the sample). Additionally, we probe the [CI](1–0) emission line for 16 galaxies (40% of the sample), and have an [CI](1–0) upper limit for 1 additional target. In Appendix C we discuss specific instances where additional archival data were available but have been discarded in our analysis because of poor quality and unreliable fluxes.

In the following, we describe the data reduction and analysis procedure in more detail, separately for the single-dish and interferometric data.

3.2. APEX

3.2.1. Observations

The APEX PI CO(2–1) observations for 32 sources of our sample (see Table B.1) were conducted between August and December 2019 (project ID E-0104.B-0672, PI: C. Cicone). Our observing strategy was to reach a line peak-to-rms ratio of S/N > 5 on the expected CO(2–1) peak flux density in velocity channels δv ∼ 5 − 50 km s⁻¹. The observations of the CO(2–1) line (${\nu }_{\mathrm{CO}(2-1)}^{\mathit{rest}}=230.538$ GHz) were performed with the receivers SEPIA180 and PI230 (similar frequency coverage as the ALMA Band 5 and 6 receivers), depending on the target’s redshift. Both PI230 and SEPIA180 are frontend heterodyne with dual-polarization sideband-separating (2SB) receivers. The instrument PI230 can be tuned within a frequency range of 195 − 270 GHz with an IF coverage of 8 GHz per sideband and with 8 GHz gap between the sidebands. The backends are fourth-generation fast Fourier transform spectrometers (FFTS4Gs) that consist of two sidebands, upper (USB) and lower (LSB), of 4 GHz (2x4 GHz bandwidth), which leads to the total bandwidth of Δν = 8GHz². The instrument SEPIA180 covers a frequency range from 159 − 211 GHz. For this instrument, the backends are the eXtended bandwidth fast Fourier transform Spectrometers (XFFTSs), which also consist of two sidebands, upper and lower, each covering 4–8 GHz, for a total of Δν = 16GHz IF bandwidth. Both receivers have an average noise temperature (T_rx) of ∼55 K (Belitsky et al. 2018).

Each sideband spectral window covered 4 GHz and was divided into 65 536 (64k) channels, resulting in a resolution of ∼61kHz, which corresponds to ∼80 − 95ms⁻¹ in velocity units at the range of redshifts covered by our sample. The CO(2–1) emission line was placed in the LSB and the telescope was tuned to the expected CO(2–1) observed frequency for each source computed by using previously known optical redshifts (see Table 1). All our PI observations were performed in the wobbler-switching symmetric mode with 60″ chopping amplitude and a chopping rate of R = 0.5Hz. The data were calibrated using standard methods. The on-source integration times (without overheads) varied from source to source and were calculated using the APEX Observing Time Calculator tool³, and are reported in Table B.2. During the observing runs, the PWV varied from 0.8 < PWV[mm] < 3.

APEX PI CO(3–2) observations were obtained between October and December 2020 (project ID E-0106.B-0674, PI: I. Montoya Arroyave). In this case, our observing strategy was to reach a line peak-to-rms ratio of S/N ∼ 7 on the expected CO(3–2) peak flux density in velocity channels δv  ∼  50 − 100 km s⁻¹. These CO(3–2) observations (${\nu }_{\mathrm{CO}(3-2)}^{\mathrm{rest}}=345.796$ GHz) were carried out with the SEPIA345 receiver, which has a frequency coverage similar to ALMA Band 7. The instrument SEPIA345, similar to SEPIA180, is a frontend heterodyne with dual-polarization 2SB receiver and works with an XFFTS backend. It can be tuned within a frequency range of 272 − 376 GHz and it has two IF outputs per polarization (two sidebands: USB and LSB), each covering 4–12 GHz, leading to a total of up to Δν = 32 GHz IF bandwidth (Meledin et al. 2022). Each sideband spectral window covered 8 GHz and was divided into 4096 channels. Initially, the requested resolution was for 65536 (64k) channels (∼122 kHz per channel), corresponding to ∼107 − 126ms⁻¹; however, due to the necessity of performing remote operations during the pandemic and to the limited band for transferring data, we applied a spectral binning at the acquisition stage so that the data could be transferred quickly to Europe after acquisition. This however did not affect the scientific output, since these extragalactic targets are characterized by broad emission lines and the new spectral resolution of ∼1953kHz (∼1.7 − 2kms⁻¹) was still very high for our science goals. The CO(3–2) emission line was placed in the LSB, and the tuning frequency was the expected CO(3–2) observed frequency plus 2 GHz: by doing so, we centered the line at IF = 8 GHz (center of sideband), rather than at IF = 6 GHz (center of backend unit) in order to have better sampling of the baselines on both sides of the line. Observations were performed in the wobbler-switching symmetric mode with 100″ chopping amplitude and a chopping rate of R = 0.5Hz, and we adopted standard calibration. The on-source integration times (calculated similarly as for the CO(2–1) observations) are reported in Table B.2. During the observing runs, the PWV varied from 0.7 < PWV[mm] < 2.5.

The [CI](1–0) APEX PI observations were obtained between October 2020 and June 2021 (project ID E-0104.B-0672, PI: C. Cicone). Our observing strategy for the atomic carbon line was to reach a S/N peak-to-rms of ∼5 on the expected [CI](1–0) peak flux density in velocity channels δv ∼ 25 − 100 km s⁻¹. The expected [CI](1–0) line flux was conservatively estimated by assuming $L_{[\mathrm{CI}](1-0)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }=0.2$, that is, the lowest value observed in NGC 6240 by Cicone et al. (2018), which turned out to be a reasonable assumption. The [CI](1–0) line observations (${\nu }_{[\mathrm{CI}](1-0)}^{\mathrm{rest}}=492.161$ GHz) were carried out with nFLASH460, which covers a similar frequency range as the ALMA Band 8 receiver. The instrument nFLASH460 is a frontend heterodyne with dual-polarization 2SB receiver with instantaneous coverage in two bands (USB and LSB) of 4 GHz each, where the separation between the center of the two sidebands is 12 GHz. It covers the frequency window between 378 and 508 GHz, and works with a fast Fourier transform spectrometer backend in each sideband. For our observations, each sideband spectral window covered 4 GHz and were divided into 65 536 (64k) channels (∼61 kHz per channel), corresponding to a resolution of ∼37 − 44 m s⁻¹ in velocity units at the range of redshifts covered by our sample. The telescope was tuned to the expected [CI](1–0) observed frequency for each source, with wobbler-switching symmetric mode with 60″ chopping amplitude and a chopping rate of R = 0.5Hz. The data were calibrated using standard procedures. The observing times, computed similarly as for CO(2–1) and CO(3–2), are reported in Table B.2. During the observing runs, the PWV varied from 0.3 < PWV[mm] < 1.0.

Additionally, we used APEX archival CO(2–1) and CO(3–2) data for part of the sample, from different projects with observing dates ranging from 2010 to 2017, using the SHeFI and nFLASH receivers. All project codes of the archival data sets used throughout this work are also reported in Table B.1. The APEX archival CO(2–1) observations used in this paper have an average S/N ∼ 5, while APEX CO(3–2) archival observations reach an average S/N ∼ 3 − 5 (computed peak-to-rms with δv ∼ 50 km s⁻¹channels).

For all single-dish data, PI and archival, we adopted the same reduction and analysis steps, which are described below in Sect. 3.2.2.

3.2.2. Data reduction

For the reduction of the APEX data sets we used the GILDAS/CLASS software package⁴, and applied the following four steps to all science targets.

First, after collecting all spectral scans of interest for a given target and transition (which could have different observing dates), we checked all scans individually, and discarded those affected by baseline instabilities and instrumental features following a similar procedure to Cicone et al. (2017). At the time of our APEX PI CO(2–1) observations, we found that the receiver SEPIA180 had slightly more stable baselines than PI230. For PI230, we verified that one of the polarization windows was heavily affected by standing waves, which led to discarding an average of 40% of sub-scans in that polarization. In the worst cases, this polarization window had to be discarded completely, therefore effectively cutting the integration time by half. For the SEPIA180 instrument, the average percentage of discarded scans was ∼20%. For the PI SEPIA345 CO(3–2) observations, the average percentage of discarded scans was 25% (an additional flagging was performed at the edge of the window to compensate for platforming issues). The nFLASH460 [CI](1–0) observations were heavily affected by instrumental features and/or sky-lines (more common at these high frequencies) and therefore the average fraction of discarded scans was ∼50% for most sources.

Second, we collected all the selected sub-scans corresponding to a given source and transition, and fitted and subtracted a linear baseline from each sub-scan after masking the central v ∈ ( − 500, 500) km s⁻¹ in order to avoid the expected line emission. We then averaged together all baseline-subtracted scans to produce a high S/N spectrum for each source.

Third, we smoothed the combined spectrum to a common velocity bin of δv ≈ 50 km s⁻¹. We fitted and subtracted a final linear baseline by using the same masking of the central v ∈ ( − 500, 500) km s⁻¹. We fitted a single Gaussian function to the spectrum and, based on the result, we refined the central masking adjusting it to the width of the detected line emission. For the non-detections, we kept the initial mask (v ∈ ( − 500, 500) km s⁻¹). These spectra were binned to δv ≈ 50 km s⁻¹resolution to compute the rms values reported in Table B.2.

Lastly, we produced the final spectrum to be used in the spectral analysis. For all the galaxies, the spectrum was extracted with a resolution of δv ≈ 5 km s⁻¹, in order to allow for further smoothing, if needed, in the spectral fitting stage. The final spectrum was imported into Python⁵, where we performed the remaining analysis.

The output signal from the APEX telescope corresponds to the antenna temperature corrected for atmospheric losses, $T_A^{\ast }$ in [K], and it must be multiplied by a calibration factor (or telescope efficiency) in order to obtain the flux density in units of Jansky [Jy]. Both for PI230 and SEPIA180 the average Kelvin to Jansky conversion factor measured during our observation period is ∼36 ± 3 Jy K⁻¹, for SEPIA345 it is ∼37.5 ± 3 Jy K⁻¹ and for nFLASH460 it is ∼58 ± 5 Jy K⁻¹. For the archival data, the calibration factor used for CO(2–1) observations was ∼39 ± 5 Jy K⁻¹ and ∼46.5 ± 7 Jy K⁻¹ for CO(3–2) observations.

3.3. ALMA and ACA

The archival data used here have observing dates ranging from 2012 to 2018. They use different antenna configurations of the 12 m (ALMA) and 7 m (ACA) arrays, corresponding to different angular resolutions and different MAS. All the project IDs used in this work are reported in Table B.1. We performed calibration and imaging using the Common Astronomy Software Applications (CASA) package⁶. The calibrated measurement sets (MSs) of data sets older than 2018 were provided by the European Southern Observatory (ESO) ALMA helpdesk. For the newer data obtained after 2018, we retrieved the MS by running the CASA pipeline (version 5.6.1) and executing the calibration script provided with each corresponding data set.

We then analyzed the MS within CASA 5.6.1 and separated the spectral windows including the lines of interest (i.e., CO(1–0), CO(2–1), CO(3–2) or [CI](1–0), depending on the data set), using the task split. A first deconvolution and cleaning were performed in interactive mode with the task tclean, by adapting the mask to the source size. This first clean provided us with an initial datacube that we used to identify the line-free channels for continuum subtraction, and to optimize the parameters for the final clean. We ran the cleaning process until reaching uniform residuals, and produced an image of the source in which we measured the noise level. We then performed the continuum subtraction using the task uvcontsub, by fitting a first-order polynomial and estimating the continuum emission in the line-free frequency ranges previously identified. We produced the final data cube using once again the task tclean on the continuum-subtracted MS file. We constructed all image cubes with the highest spectral resolution available, ranging from Δv∼1 to ∼10 km s⁻¹, depending on the data set. Final cubes are obtained using Briggs weighting with robust parameter equal to 0.5 and primary-beam corrected. We extracted the final spectrum from a circular aperture that is size-matched to maximize the recovered flux. The apertures used for spectra extraction are reported in Table B.2. For the sources where only high-resolution ALMA observations were available (IRAS F01572+0009 and IRAS F12072−0444), we applied an uv tapering to enhance the sensitivity to extended structures. This however does not overcome the possible issue of missing flux from faint extended structures due to poor sampling of short uv baselines.

The spectra are exported from CASA in flux density [Jy] units, extracted in suitable format and imported into Python for further analysis. The quoted errors refer to the systematic errors on the absolute flux calibration of ALMA/ACA data (estimated to be 5% for Band 3 data and 10% for Bands 6–8 data, and typically the dominant source of error in the data used in this study), added in quadrature to the statistical rms of the spectra.

4. Methodology

4.1. Spectral line fitting

Our analysis is aimed at deriving source-averaged line ratios for different kinematic components of the molecular and atomic ISM in local (U)LIRGs. We also want to investigate possible statistical trends between the molecular (CO and [CI]) line ratios as a function of the central velocity v and line width σ_v of the different components, to understand whether broader and/or higher-v H₂ gas components are the origins for the extremely high global CO excitation previously found in (U)LIRGs (e.g., Papadopoulos et al. 2012), which Cicone et al. (2018) suggested based on their pilot study on NGC 6240.

We base the analysis reported in this paper exclusively on total (i.e., galaxy-integrated) molecular line spectra, and on the results of a multi-Gaussian spectral fitting of the CO and [CI](1–0) lines. We acknowledge that, without spatially resolved information, it is not possible to link in a straightforward way the different (e.g., broad or narrow) spectral line components to outflowing or non-outflowing gas. However, performing such a classification is not needed in our case, because we can rely on a large, statistically significant sample. Indeed, we are interested in studying in a statistical sense any trends observed between molecular line ratios and the central velocity (v) and/or velocity dispersion (σ_v) of the different spectral components. Once the presence of statistical correlations is assessed (independently of any arbitrary classification of such components in terms of outflow or disk), we can interpret the results by assuming that, in this sample of (U)LIRGs, the line luminosities of high-σ_v and/or high-v spectral components are more likely to be dominated by molecular gas embedded in outflows compared to the low-σ_v and low-v components. Such an assumption would be supported by (i) the results of the pilot study on NGC 6240 (see Cicone et al. 2018), (ii) the unambiguous detection of OH119 μm outflows in most targets (Sturm et al. 2011; Veilleux et al. 2013; Spoon et al. 2013), and (iii) the widespread evidence for high-velocity molecular outflows in local (U)LIRGs reported in the recent literature (see the review by Veilleux et al. 2020).

After having clarified our strategy, we now describe our spectral fitting procedure. We used mpfit for Python⁷. This tool uses the Levenberg-Marquardt technique to solve the least-squares problem in order to fit a user-supplied function (the model) to the user-supplied data points (the data). Spectral lines were modeled using single or multiple Gaussian profiles characterized by an amplitude, a peak position (v_cen), and velocity dispersion (σ_v, or, equivalently, the full width at half maximum, FWHM). Having multiple CO transitions allows us to partially break any degeneracy in the spectral line decomposition with Gaussian functions. We therefore fitted simultaneously all CO transitions available for each source, by constraining the v_cen and σ_v of the Gaussian components to be equal in all CO transitions, allowing only their amplitudes to vary freely in the fit. We allowed the fit to use up to a maximum of three Gaussian functions to reproduce the observed global line profiles, so that line asymmetries and broad wings are properly captured with separate spectral components when the S/N is high enough (see for example IRAS 13120-5453 in Fig. B.2). For all the fits, we verified using a reduced χ² criterion that a fourth Gaussian component was not required for any of the sources, with the data at hand. In those cases where the statistical criterion (reduced χ²) does not indicate a clear preference between a fit performed with one, two, or three Gaussians, we performed a visual inspection of the fit. For the low S/N spectra (e.g., S/N ∼ 3) without clear asymmetries, we used only one Gaussian component, in order not to over-fit the data.

The fitting procedure for the [CI](1–0) line spectra is carried out separately from the CO lines due to their overall lower S/N. Furthermore, by doing so, we can avoid assuming a priori that the atomic carbon and CO lines trace the same gas clouds and share the same kinematics. This assumption will be tested and discussed in Sect. 5.1.2.

The final fits for all sources and transitions are shown in Figs. B.1–B.6, where the spectral transitions are color-coded. Based on the fit results, we computed velocity-integrated line fluxes for both the individual Gaussian components and the entire line profiles, and the latter are reported in Table D.1. As a sanity check, we verified that the total line fluxes measured through the fit (by adding up the individual Gaussians) are consistent with the total line fluxes calculated by directly integrating the spectra within v ∈ ( − 1000, 1000) km s⁻¹, after setting a threshold of > 2σ for each channel. We find the values to be consistent within the errors, and therefore using either value will not affect the analysis performed throughout our work.

4.2. Line luminosities and ratios

We calculated the CO and [CI] line luminosities from the integrated line fluxes following the definition from Solomon et al. (1997):

$$\begin{array}{c}\hfill L_{\mathrm{line}}^{\prime }[\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^2]=3.25\times {10}^7\frac{D_{\mathrm{L}}^2}{{\nu }_{\mathrm{obs}}^2{(1+z)}^3}{\displaystyle \int }{S}_v\mathrm{d}v,\end{array}$$(1)

where D_L is the luminosity distance measured in [Mpc], ν_obs is the corresponding observed frequency in [GHz], and ∫S_v dv is the total integrated line flux in [Jy km s⁻¹]. In Table D.1 we report the total integrated fluxes and respective luminosities calculated for the different transitions. The CO line ratios are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill r_{21}& =L_{\mathrm{CO}(2-1)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }\hfill \\ \hfill r_{31}& =L_{\mathrm{CO}(3-2)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }\hfill \\ \hfill r_{32}& =L_{\mathrm{CO}(3-2)}^{\prime }/L_{\mathrm{CO}(2-1)}^{\prime }.\hfill \end{array}\end{array}$$(2)

In our analysis we use both the global CO line luminosity ratios as well as those computed for individual Gaussian components. Additionally, we calculate global [CI](1–0)/CO(1–0) line luminosity ratios, as

$$\begin{array}{c}\hfill r_{\mathrm{CICO}}=L_{[\mathrm{CI}](1-0)}^{\prime }/L_{\mathrm{CO}(1-0)}^{\prime }.\end{array}$$(3)

Line ratios were computed for all combinations of lines and sources where the individual line luminosity measurements pass a loose criterion of S/N > 1, in order not to penalize cases where one of the two lines is constrained at very high significance. As a result, some line ratios have very large error bars, which are taken into account in our analysis.

5. Results

5.1. Atomic carbon as an alternative gas tracer

5.1.1. A tight relation between CO(1–0) and [CI](1–0) luminosities

In Fig. 1 we plotted the measured total CO(1–0) and CO(2–1) line luminosities as a function of [CI](1–0) line luminosity, for the ULIRGs in our sample. The relation with CO(3–2) was not studied as this line starts to trace denser and more excited H₂ gas rather than the global molecular reservoir, while we are interested in exploring the potential of [CI](1–0) to probe similar regions as the CO J_up = 2, 1 transitions.

Fig. 1. CO(1–0) vs. [CI](1–0) luminosity (top) and CO(2–1) vs. [CI](1–0) luminosity (bottom) for the ULIRGs in our sample. The best-fit relations are shown as dashed orange lines. The best-fit parameters are reported at the bottom-right corner of the plots. We also display the Pearson correlation coefficients (ρ) and their associated p-values. The solid black lines in both panels represent the corresponding relations reported by Jiao et al. (2017) for a sample of 71 (U)LIRGs, and the dotted lines represent the relations of Jiao et al. (2019) for a sample of 15 nearby spiral galaxies, between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$.

The plots in Fig. 1 show that $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{\text{CO}(\text{2-1})}^{\prime }$ are both tightly correlated with $L_{[\text{CI}](\text{1-0})}^{\prime }$, showing Pearson correlation coefficients (ρ) equal to 0.77 and 0.71, and p-values of 9.1 × 10⁻³ and 4.7 × 10⁻³, respectively, for $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{\text{CO}(\text{2-1})}^{\prime }$. We performed a fit to the relation between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$ using least squares, which gives

$$\begin{array}{c}\hfill {\log }_{10}{L}_{\mathrm{CO}(1-0)}^{\prime }=(4.03\pm 2.36)+(0.64\pm 0.26){\log }_{10}{L}_{[\mathrm{CI}](1-0)}^{\prime }.\end{array}$$(4)

We performed a similar fit with the CO(2–1) line, whose results are reported on the corresponding plot. We find that $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{\text{CO}(\text{2-1})}^{\prime }$ follow very similar relations as a function of $L_{[\text{CI}](\text{1-0})}^{\prime }$, with variations in the best-fit parameters within one standard deviation.

In Fig. 1 we also report the relation found by Jiao et al. (2017) in a study of unresolved neutral carbon emission in a sample of 71 (U)LIRGs based on Herschel observations, for which they derive a best-fit relation of ${\text{log}}_{\text{10}}{\text{L}}_{\text{CO}(\text{1-0})}^{\prime }\text{=}(\text{-0}.\text{19}\pm \text{1}.\text{26})\text{+}(\text{1}.\text{09}\pm \text{0}.\text{15})$ ${\text{log}}_{\text{10}}{\text{L}}_{[\text{CI}](\text{1-0})}^{\prime }$. In a similar study performed on 15 nearby spiral galaxies with spatially resolved Herschel data, Jiao et al. (2019) obtain: log${}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(0.74\pm 0.12)+(1.04\pm 0.02)$ ${\text{log}}_{\text{10}}{\text{L}}_{[\text{CI}](\text{1-0})}^{\prime }$. Our best-fit $L_{[\text{CI}](\text{1-0})}^{\prime }$ versus $L_{\text{CO}(\text{1-0})}^{\prime }$ relation has a flatter slope than the one found by Jiao et al. (2017), likely due to our sample only covering a narrower dynamic range in luminosities than the sample in Jiao et al. (2017), while our sources are exclusively ULIRGs, the Jiao et al. (2017) sample is heavily dominated by LIRGs (62 LIRGs and only 9 ULIRGs). When we explore the same relations in our extended sample (including the LIRGs), shown in Fig. A.1, we find that the best-fit relation between $L_{[\text{CI}](\text{1-0})}^{\prime }$ and $L_{\text{CO}(\text{1-0})}^{\prime }$ (${\text{log}}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(-0.45\pm 0.73)+(1.13\pm 0.08){\text{log}}_{10}{L}_{[\text{CI}](\text{1-0})}^{\prime }$) is almost linear and well in agreement with that obtained by Jiao et al. (2017) and clearly shifted to higher $L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$ ratios with respect to the Jiao et al. (2019) fit performed on non-infrared luminous local galaxies. Such a difference in $L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$ ratios between (U)LIRGs and other galaxies is further explored in Sect. 5.6.

The tight correlations in Fig. 1 suggest that the CO(1–0) and CO(2–1) lines arise from similar regions as the [CI](1–0) emission, at least when averaged over galactic scales, and strengthen the hypothesis (see, e.g., Papadopoulos et al. 2004) that the [CI](1–0) line is an excellent molecular gas tracer, and a valid alternative to low-J CO line emission.

5.1.2. Comparison between CO and [CI] line widths

In the previous section we find a tight relation between the [CI](1–0) and CO total line luminosities for ULIRGs (which becomes almost linear when expanding the dynamic range in luminosity values by including the LIRGs). Here we test whether the two tracers share the same kinematics. We compared the line widths using the CO(2–1) and the [CI](1–0) lines. We preferred CO(2–1) over CO(1–0) to maximize the sample size, since CO(2–1) spectra are available for all 16 sources with a [CI](1–0) detection.

We used two different approaches to study the line widths of the two tracers. Firstly, we performed a dedicated, single Gaussian spectral fit, run independently for each line. Secondly, we computed the 16–84 (v₈₄ − v₁₆) and 2.3–97.7 (v_97.7 − v_2.3) percentile velocity intervals, derived from the analytical form of the overall best-fit line profile obtained by a multi-Gaussian fit. The latter approach is likely a more robust method when dealing with complex line profiles, as is the case for some of the sources in our sample, for example IRAS 19254–7245.

The results are shown in Fig. 2. The plot of σ_v, [CI] versus σ_v, CO obtained through a single Gaussian fit is plotted on the top panel of Fig. 2, while the percentile velocity plots obtained using the second approach are shown in the bottom panels of Fig. 2. We find that the sources characterized by broader profiles sit preferentially below the 1:1 relation. The best-fit relations obtained through a least squares regression analysis have a slope below unity.

Fig. 2. [CI](1–0) line width as a function of the CO(2–1) line width. Top panel: Velocity dispersion (σv) obtained via a single Gaussian spectral fit to the [CI](1–0) and CO(2–1) emission lines. Middle and bottom panels: 16–84 (v84 − v16) and 2.3–97.7 (v97.7 − v2.3) percentile velocity intervals for CO and [CI], respectively, derived from the best-fit function obtained from the multi-Gaussian fit. Orange data points correspond to the ULIRGs sample analyzed in this work. The purple pentagon and the pink dot represent the values obtained for NGC 6240, respectively for the total region and for the inner ∼2″ (see Cicone et al. 2018, and the text in Sect. 5.1.2 for a more detailed explanation of the apertures used for spectra extraction). The solid black line indicates the 1:1 relation. The dashed orange line shows the best-fit relation obtained using only our sample, while the dashed purple line is the best-fit relation obtained for our sample plus NGC6240 (total). The shaded gray areas corresponds to the 1σ confidence interval of the two fits. At the bottom-right corner we report the Pearson correlation coefficient (ρ) and the p-values, as well as the best-fit coefficients of the purple fit.

As an additional test, we over-plotted in Fig. 2 the values obtained for NGC 6240, which is a source characterized by an extremely turbulent ISM, strongly affected by outflows. The purple pentagon represents the total measurement available for NGC 6240, computed from spectra extracted from a 12″ × 6″ rectangular aperture encompassing the nuclei and the molecular outflow, while the pink dot represents the central 2″ × 2″ region (for a more in depth explanation of how the apertures are defined, see Cicone et al. 2018). Quite strikingly, both NGC 6240 data points sit on the best-fit relation obtained for our sample when probing the core of the lines via σ_v and v₈₄ − v₁₆. Instead, as we probe more toward the high-velocity wings of the line by using, e.g., the v_97.7 − v_2.3 values, the nuclear spectrum results to be more consistent with the 1:1 relation between [CI] and CO line width, while the total spectrum of NGC 6240, including the extended outflows, sits on the best-fit relation obtained from the analysis of the other sources. The fact that the total spectrum includes more of the extended outflow than the nuclear one (see Cicone et al. 2018), and is also the one that departs more from the 1:1 relation, may indicate that this deviation (i.e., a narrower width of [CI] with respect to CO) is accentuated by the inclusion of diffuse outflowing gas.

The fit that includes the total emission from NGC 6240 displayed in the top panel of Fig. 2 which is consistent with the one obtained from our sample alone, is

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{v},[\mathrm{CI}]}=(30\pm 9)+(0.64\pm 0.07){\sigma }_{\mathrm{v},\mathrm{CO}},\end{array}$$(5)

and the corresponding fit for the v_97.7 − v_2.3 velocity percentiles (shown in the bottom panel) is

$$\begin{array}{c}\hfill {\mathrm{v}}_{97.7}-{\mathrm{v}}_{2.3}[\mathrm{CI}]=(74\pm 66)+(0.70\pm 0.12){\mathrm{v}}_{97.7}-{\mathrm{v}}_{2.3}\mathrm{CO}.\end{array}$$(6)

Therefore, our data indicate that the [CI](1–0) line is narrower than CO(2–1). The average line width ratio is ⟨σ_v, [CI]/σ_v, CO⟩=⟨r_σ⟩ = 0.91 ± 0.07, computed using all sources. The ratios below unity are driven by targets with σ_v, CO > 150 km s⁻¹, while those with σ_v, CO ≲ 150 km s⁻¹, which represents the majority of our sample, are consistent with the 1:1 relation.

Few comparisons of CO and [CI] line widths can be found in the literature, and most of these previous studies report a 1:1 correspondence between the line widths. Michiyama et al. (2021), by comparing ACA CO(4–3) and [CI](1–0) observations of a sample of 36 local (U)LIRGs, found a 1:1 relation between the FWHMs of the two transitions. However, their analysis excludes sources with complex profiles (e.g., double peak emission), which we did not do. Similarly, Bothwell et al. (2017) analyzed ALMA [CI](1–0) emission line observations in a sample of strongly lensed dusty star-forming galaxies spanning a wide redshift range of 2 < z < 5, and for 11 of such sources they compared [CI] and CO(2–1) line widths, using literature CO data. Their results are consistent with a 1:1 relation. In Sect. 6 we discuss possible explanations for the difference in CO and [CI] line widths observed in our sample and specifically in the high-σ_v (U)LIRGs.

5.2. Molecular gas mass estimates and the CO-to-H₂ factor

Building upon Sect. 5.1.1, we use the [CI](1–0) and CO(1–0) emission lines to derive independent estimates of the molecular gas mass (M_mol) of the ULIRGs of our sample. We also use the [CI]-based M_mol to derive an average value for the α_CO factor, similarly to Cicone et al. (2018).

Both tracers rely on calibration factors in order to compute M_mol. For [CI](1–0)-based estimates, we need to assume the optically thin condition (which applies to most extragalactic environments), a value for the [CI] abundance with respect to H₂ (X_CI = [C/H₂]), and a value for the parameter Q₁₀ (i.e., the [CI] excitation factor). The molecular hydrogen gas mass can then be computed, following Dunne et al. (2021), as

$$\begin{array}{c}\hfill M_{\mathrm{mol},[\mathrm{CI}]}[M_{\odot }]=1.36\frac{(9.51\times {10}^{-5})}{{\mathrm{X}}_{\mathit{CI}}{\mathrm{Q}}_{10}}{L}_{[\mathrm{CI}](1-0)}^{\prime }.\end{array}$$(7)

For CO-based mass measurements, we need to assume an α_CO factor:

$$\begin{array}{c}\hfill M_{\mathrm{mol},\mathrm{CO}}[M_{\odot }]={\alpha }_{\mathrm{CO}}{L}_{\mathrm{CO}(1-0)}^{\prime }.\end{array}$$(8)

Both Eqs. (7) and (8) include the Helium contribution to the molecular gas mass through a multiplicative factor of 1.36.

All tracers of H₂ are affected by uncertainties. Using CO and an arbitrary α_CO value can introduce large errors in the computed M_mol due to its sensitivity to metallicity in a nonlinear fashion, and to the turbulence and kinematics of the CO-emitting clouds that affect the global optical depth of galaxy-averaged CO measurements. In the case of [CI], different combinations of X_CI and Q₁₀ may yield different results; X_CI may be the easiest parameter to model in terms of the ISM conditions if in fact the C⁰ abundance is determined by CRs (e.g., Bisbas et al. 2015; Dunne et al. 2021). In the following Sects. 5.2.1 and 5.2.2, we discuss separately the [CI]-based and CO-based M_mol estimates.

5.2.1. [CI]-based M_mol estimates

We compute M_mol, [CI] using the [CI](1–0) luminosities reported in Table D.1 and Eq. (7). We adopt a carbon abundance of X_CI = (3.0 ± 1.5)×10⁻⁵, which is an appropriate value for local star-forming galaxies and has been used by several previous studies (e.g., Weiß et al. 2005; Papadopoulos et al. 2004; Walter et al. 2011; Jiao et al. 2017; Cicone et al. 2018). For the [CI] excitation factor Q₁₀, we adopt a value of 0.48 (with < 16% variation), following the prescriptions by Papadopoulos et al. (2022), who find that, for the most expected average ISM conditions in galaxies (n_H2 = [300 − 10⁴] cm⁻³ and T_kin = [25 − 80] K), the [CI] lines are globally sub-thermally excited.

By defining a parameter α_[CI] analogous to α_CO to represent a [CI](1–0)-to-H₂ conversion factor, we can rewrite Eq. (7) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M_{\mathrm{mol},[\mathrm{CI}]}=& {\alpha }_{[\mathrm{CI}]}{L}_{[\mathrm{CI}](1-0)}^{\prime }\hfill \\ \hfill & \mathrm{with}{\alpha }_{[\mathrm{CI}]}=1.293\times {10}^{-4}{({\mathrm{X}}_{\mathrm{CI}}{\mathrm{Q}}_{10})}^{-1}.\hfill \end{array}\end{array}$$(9)

And thus, plugging in our assumptions for the X_CI and Q₁₀ values, we obtain α_[CI] = 9.0 M_⊙ (K km s⁻¹ pc²)⁻¹.

The [CI]-based M_mol values can then be used to infer α_CO as follows:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{CO}}=M_{\mathrm{mol},[\mathrm{CI}]}/L_{\mathrm{CO}(1-0)}^{\prime }.\end{array}$$(10)

The ratio $M_{\text{mol},[\text{CI}]}/L_{\text{CO}(\text{1-0})}^{\prime }$, which in practice represents the α_CO factor required to force agreement between [CI]- and CO-based H₂ mass estimates, is plotted in Fig. 3 against the infrared luminosity, for the 10 ULIRGs with available CO(1–0) and [CI](1–0) detections. The distribution of the resulting α_CO values is also shown on the right of Fig. 3.

Fig. 3. Ratio between [CI]-based molecular gas mass estimates and CO(1–0) line luminosity, computed for the ULIRGs in our sample that have both lines available, and plotted as a function of the infrared luminosity. This ratio is interpreted as the αCO factor (see Eq. (10)). The right part of the plot shows the resulting distribution of [CI]-based αCO values. The dotted line indicates the CO-to-H2 conversion factor for the Milky Way, and the dashed line corresponds to the value commonly used in the literature for (U)LIRGs (Downes & Solomon 1998). The mean and median values obtained for our sample (reported at the top-right corner) are αCO = 1.9 ± 0.4 M⊙ and 1.7 ± 0.5 M⊙ (K km s−1 pc2)−1, respectively.

Figure 3 demonstrates that 9 out of 10 targets require an α_CO value higher than the one commonly assumed for (U)LIRGs of 0.8 M_⊙ (K km s⁻¹pc²)⁻¹, see Downes & Solomon (1998). The mean value measured for our sample is 1.9 ± 0.4 M_⊙, and the median value is 1.7 M_⊙ (K km s⁻¹ pc²)⁻¹, with 16–84 velocity percentile equal to 1.2 – 2.8 (K km s⁻¹ pc²)⁻¹. Our results are consistent with the dust-based α_CO estimate equal to $1.8_{-0.8}^{+1.3}{M}_{\odot }$ (K km s⁻¹ pc²)⁻¹) derived by Herrero-Illana et al. (2019) for 55 (U)LIRGs from the Great Observatories All-sky LIRG survey (GOALS). These authors estimate first the dust mass (M_dust) from a FIR spectral energy distribution (SED) fit using Herschel data, following the strategy proposed by Scoville et al. (2016) of fixing T_dust = 25 K for every source, and then estimate α_CO by requiring that the gas-to-dust mass ratio of (U)LIRGs matches the one of local star-forming spirals. Similarly, a study performed by Kawana et al. (2022) on a nearby LIRG (NGC 3110) estimates an α_CO value of 1.7 ± 0.5 M_⊙ (K km s⁻¹ pc²)⁻¹ based on thermal dust continuum emission and assuming that the dust and the rotational temperature of the CO molecule are equal.

In Appendix A we extend the calculations to include the LIRGs in our sample with available [CI](1–0) and CO(1–0) data (i.e., two additional sources; see Fig. A.2). The resulting mean and median values are perfectly consistent with the ones obtained in Fig. 3.

5.2.2. CO-based M_mol estimates

We then use the lowest-J CO transition available for each source to estimate a CO-based total molecular gas mass. For sources without CO(1–0) data, we first estimate $L_{\text{CO}(\text{1-0})}^{\prime }$ based on the $L_{\text{CO}(\text{2-1})}^{\prime }$ transition, which is available for the whole sample. We assume a line ratio of r₂₁ = 1.1 ± 0.4, which is the median value computed in this work based on the data available in our sample (see Sect. 5.4 and Fig. 6). We then proceed to compute M_mol, CO using Eq. (8). We compute M_mol, CO using three different α_CO values, namely, 4.3, 0.8 and 1.7 M_⊙ (K km s⁻¹pc²)⁻¹, corresponding to the value measured for the Milky Way galaxy (Bolatto et al. 2013), the value commonly employed in the past for (U)LIRGs (Downes & Solomon 1998), and the median value obtained for our sample using the [CI]-based method described previously (see Sect. 5.2.1).

In Fig. 4 we show the comparison between the [CI]-based M_mol values, and the CO-based M_mol estimates obtained using different conversion factors. These plots report the same result as Fig. 3 (i.e., the α_CO value), but visualized in a different way and including the uncertainties. For each relation we perform two linear fits, one with a free-varying slope and another with a slope fixed to unity. The 1:1 relation is also over-plotted using a dotted black line. The plots in Fig. 4 show clearly that an α_CO = 4.3 M_⊙ (K km s⁻¹pc²)⁻¹ overestimates the M_mol, while an α_CO = 0.8 M_⊙ (K km s⁻¹pc²)⁻¹ underestimates it, for all ULIRGs of the sample. Instead, as is obvious from the definition, the value of α_CO = 1.7 M_⊙ (K km s⁻¹pc²)⁻¹, which is the median of the individual [CI]-based α_CO values estimated in the previous section, brings all data points closer to the 1:1 relation. We note, however, that this is an average value for the CO-to-H₂ conversion factor, and it is most likely to vary from galaxy to galaxy (as shown by a few outliers visible in Fig. 4), as well as a function of aperture size and spatial scales probed. This variability averages out on galaxy-averaged measurements leading to a typical value corresponding to the dominant radiation-emitting regions within the galaxy. We note that large uncertainties on empirically estimated α_CO values are still expected, as there are many factors that may impact on its value, for example density, temperature, metallicity, and optical depth. These results are however reassuring and indicate that the adoption of α_CO  ≃  1.7 M_⊙ (K km s⁻¹pc²)⁻¹ for local (U)LIRGs is reasonable.

Fig. 4. [CI]-based Mmol vs. CO-based Mmol values obtained with different αCO assumptions (0.8, 4.3, and 1.7 M⊙ (K km s−1pc2)−1). The best-fit linear relation with a free-varying slope is shown as a dashed orange line, and the best-fit parameters are reported at the bottom-right corner of each plot. The solid orange line shows the best-fit relation obtained by fixing the slope to unity, and the 1:1 relation is reported as a dotted black line.

5.3. Total line luminosities as a function of galaxy properties

Before investigating molecular line ratios and their dependencies on galaxy properties, it is worth exploring first the trends involving the line luminosities that are used to compute such ratios. We recall that our sample is, by construction, biased toward high L_IR. This could cause underlying galaxy scaling relations not to be properly captured by our targets, or to be detected with different slopes compared to the global star-forming galaxy population (as already seen, e.g., in Fig. 1), due to the limited range of intrinsic properties (e.g., SFR) probed by local (U)LIRGs (see the discussion on scaling relations in, e.g., Cicone et al. 2017). It is therefore important to identify the portion of the $L_{\mathrm{line}}^{\prime }-L_{\mathrm{IR}}$ (or $L_{\mathrm{line}}^{\prime }$- SFR, $L_{\mathrm{line}}^{\prime }-L_{\mathrm{AGN}}$) parameter space occupied by our sources in order to place our results into perspective. To this aim, Fig. 5 shows the CO(1–0), CO(2–1), CO(3–2), [CI](1–0) line luminosities as a function of L_IR, SFR, and L_AGN, for all targets with corresponding line measurements available (see Table D.1).

Fig. 5. Global CO(1–0), CO(2–1), CO(3–2) (top panels) and [CI](1–0) (bottom panels) line luminosity plotted as a function of LIR (left), SFR (middle), and LAGN(right) for our sample of ULIRGs. In each plot, the dashed lines are the best-fit relations obtained from a least squares regression analysis conducted for each transition separately (color-coded according to the transition; see the legend in the top-left panel). In the bottom panels, the shaded gray areas correspond to the 1σ confidence interval of the fit. The bottom panels also report the Pearson correlation coefficients (ρ) and their associated p-values. The top-left panel also shows the $L_{\text{CO}(\text{1-0})}^{\prime }$ − LIR relations found by Sargent et al. (2014) for local main-sequence (dot-dashed black line) and SB galaxies (dotted black line), and the $L_{\text{CO}(\text{1-0})}^{\prime }$ − LTIR relation obtained by Herrero-Illana et al. (2019) using IRAM 30m CO(1–0) observations of 55 local sources in the GOALS sample, which comprises mostly LIRGs with LIR < 1012 L⊙ and with a < 20% AGN contribution to LIR. In the top-middle panel, for better visualization, we over-plot with darker colors the data in bins of SFRs. In this plot, the solid purple line represents the $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR best fit obtained by Cicone et al. (2017) for the ALLSMOG and COLDGASS samples of local star-forming galaxies, described by ${\text{log}}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(8.16\pm 0.04)+(1.34\pm 0.07){\text{log}}_{10}\text{SFR}$.

For the $L_{\mathrm{CO}}^{\prime }$ versus L_IR relations (top-left panel of Fig. 5), we measure correlation coefficients of ${\rho }_{{\mathrm{L}}_{\mathrm{CO}(1-0)}^{\prime }-{\mathrm{L}}_{\mathrm{IR}}}=0.41$, ${\rho }_{{\mathrm{L}}_{\mathrm{CO}(2-1)}^{\prime }-{\mathrm{L}}_{\mathrm{IR}}}=0.45$, ${\rho }_{{\mathrm{L}}_{\mathrm{CO}(3-2)}^{\prime }-{\mathrm{L}}_{\mathrm{IR}}}=0.38$, with p-values = 7.6 × 10⁻², 6.4 × 10⁻³, 4.1 × 10⁻², respectively). These relations trace, essentially, the Schmidt-Kennicutt (S-K) law (Schmidt 1959; Kennicutt 1998); the correlation coefficients, although hinting of positive relations, do not show significant relation, most likely caused by the narrow range in L_IR probed by our sample. Such an hypothesis is strengthened as the ρ coefficients (and their respective p-values) yield much tighter relations between the quantities when the LIRGs in our sample are included in the analysis (see Fig. A.3). Our sample lies close to the $L_{\text{CO}(\text{1-0})}^{\prime }$ − L_IR relation obtained by Sargent et al. (2014) for local SBs (see the top-left panel of Fig. 5), which is offset by ≃0.46 dex from the main-sequence galaxies’ relation. We note that our sample is instead significantly offset with respect to the Herrero-Illana et al. (2019) relation, we ascribe this discrepancy to a combination of two factors: (i) the extrapolation of a relation that is based on a lower-L_IR sample than ours; and (ii) their use of a different method for computing the L_IR based on an SED fitting, which delivers lower L_IR (by ∼0.5 dex) per given $L_{\text{CO}(\text{1-0})}^{\prime }$. We verified that for the four targets in common with Herrero-Illana et al. (2019), the CO fluxes are consistent, but the L_IR computed through their SED fitting (reported in their Table 5) are 0.45–1.0 dex lower than our L_IR values, which are instead consistent with the L_IR values reported by Armus et al. (2009) for the same sources. The best-fit $L_{\text{CO}(\text{1-0})}^{\prime }$ − L_IR relation obtained by running a least square regression analysis on our data is ${\log }_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(0.8\pm 2.2)+(0.74\pm 0.18){\log }_{10}{L}_{\text{IR}}$, with $L_{\text{CO}(\text{1-0})}^{\prime }/L_{\text{IR}}$ ratios similar to the SB sample of Sargent et al. (2014). This is not surprising since all of the sources in our sample of ULIRGs show enhanced star formation (see Table 1).

The bottom-left panel of Fig. 5 reports $L_{[\text{CI}](\text{1-0})}^{\prime }$ as a function of L_IR. The [CI](1–0) luminosities span a range ${log}_{10}{L}_{[\text{CI}](\text{1-0})}^{\prime }=8.7-9.8$ [K km s⁻¹pc²], consistently lower than the CO(1–0) luminosities (${log}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=9.4-10.4$ [K km s⁻¹pc²]). Figure 1 shows a tight relation between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$, and hence a similar relation to that found between $L_{\text{CO}(\text{1-0})}^{\prime }$ (or $L_{\text{CO}(\text{2-1})}^{\prime }$) and L_IR is expected. Indeed, we measure a slightly higher Pearson correlation coefficient of ${\rho }_{L_{[\text{CI}](\text{1-0})}^{\prime }-L_{\text{IR}}}=0.56$ (p-value = 0.04).

The middle panels of Fig. 5 display the CO and [CI] line luminosities as a function of SFR. These relations are not much dissimilar from those with L_IR (left panels), as expected since much of the L_IR in (U)LIRGs is powered by star formation. Since the SFRs have been computed by removing the contribution to L_IR estimated to arise from AGN-heated dust, the middle panels of Fig. 5 should more truthfully trace the S–K relation. However, we struggle to retrieve a tight S–K law for this sample, probably because of selection biases due to their narrow distribution in SFRs, combined with the inevitably large uncertainties on the AGN fraction (Veilleux et al. 2009b). Although we measure slightly higher Pearson correlation coefficients for the $L_{\text{CO}}^{\prime }$ versus SFR relations (${\rho }_{L_{\mathrm{CO}(1-0)}^{\prime }-\mathrm{SFR}}=0.42$, ${\rho }_{L_{\mathrm{CO}(2-1)}^{\prime }-\mathrm{SFR}}=0.46$, ${\rho }_{L_{\mathrm{CO}(3-2)}^{\prime }-\mathrm{SFR}}=0.50$, with p-values = 6.8 × 10⁻², 5 × 10⁻³ and 5.8 × 10⁻³, respectively) than for the $L_{\text{CO}}^{\prime }$ versus L_IR relations, the values are still marginal at best for their correlations. When using [CI](1–0) as a H₂ tracer, we obtain ${\rho }_{L_{[\mathrm{CI}](1-0)}^{\prime }-\mathrm{SFR}}=0.47$ and p-value = 8.7 × 10⁻². The best-fit $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR relation (${log}_{10}{L}_{\mathrm{CO}(1-0)}^{\prime }=(9.04\pm 0.22)+(0.44\pm 0.11){log}_{10}\mathrm{SFR}$) has a flatter slope than the $L_{\text{CO}(\text{1-0})}^{\prime }$-L_IR one. Similarly, the $L_{[\text{CI}](\text{1-0})}^{\prime }$-SFR relation (reported on the plot), also shows a considerably shallower slope than the $L_{[\text{CI}](\text{1-0})}^{\prime }$-L_IR relation, with a value of 0.6 ± 0.3. In the top-middle panels of Fig. 5, we plot the best-fit $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR relation obtained for a much more unbiased sample of local star-forming main-sequence galaxies (drawn from the COLDGASS and ALLSMOG surveys; see Cicone et al. 2017)⁸. Our ULIRGs are characterized by significantly lower $L_{\text{CO}(\text{1-0})}^{\prime }$/SFR ratios than main-sequence galaxies, especially in the high-SFR regime. The result of lower $L_{\text{CO}(\text{1-0})}^{\prime }$/SFR ratios in (U)LIRGs than in normal main-sequence galaxies is well known, and due to a combination of (i) a higher star formation efficiency (SFE; i.e., higher SFE=SFR/M_mol, or equivalently lower τ_dep = M_mol/SFR) and (ii) a lower α_CO compared to normal disk galaxies, as also confirmed by our analysis in Sect. 5.2.

The right panels of Fig. 5 show the CO and [CI] line luminosities as a function of AGN luminosity. Neither of the two gas tracers, in any of the transitions, show any correlation with L_AGN. The measured Pearson correlation coefficients and p-values for the CO emission lines are ${\rho }_{L_{\text{CO}(\text{1-0})}^{\prime }-L_{\text{AGN}}}=-0.05$ (p-value = 0.8), ${\rho }_{L_{\text{CO}(\text{2-1})}^{\prime }-L_{\text{AGN}}}=-0.06$ (p-value = 0.7) and ${\rho }_{L_{\text{CO}(\text{3-2})}^{\prime }-L_{\text{AGN}}}=-0.2$ (p-value = 0.3). These results are opposed to what has been found in the literature by, for example, Husemann et al. (2017) and Shangguan et al. (2020a), who find signs of correlation between the CO(1–0) (or CO(2–1)) line luminosity and L_AGN on samples of tens of local quasars, although the interpretation is not entirely clear. Indeed, while it is generally accepted that the S-K law traces a fundamental causal relation between the availability of fuel and the resulting star formation activity, the relation between $L_{\mathrm{CO}}^{\prime }$ and L_AGN can be only investigated globally, because it relates two processes that affect very different scales in galaxies. We note, however, the change in the relation when the LIRGs in our sample are included (see Fig. A.3), for which we obtain Pearson correlation coefficients and p-values of ${\rho }_{L_{\text{CO}(\text{1-0})}^{\prime }-L_{\text{AGN}}}=0.53$ (p-value = 0.01), ${\rho }_{L_{\text{CO}(\text{2-1})}^{\prime }-L_{\text{AGN}}}=0.37$ (p-value = 0.02), and ${\rho }_{L_{\text{CO}(\text{3-2})}^{\prime }-L_{\text{AGN}}}=0.32$ (p-value = 0.09). Hence, the inclusion of LIRGs, with the consequent widening of the dynamic range in properties probed, produces positive trends between molecular line luminosities and L_AGN, in particular, for CO(1–0) and CO(2–1) (although weak); possibly originating from an underlying scaling of both quantities with SFR and/or M_* (see Cicone et al. 2017). It is, however, necessary to confirm these results with a sample including more sources in the L_IR < 12 L_⊙ regime.

The relation between $L_{[\text{CI}](\text{1-0})}^{\prime }$ and L_AGN in our sample of ULIRGs (bottom-right panel of Fig. 5) is equally weak as the relations for CO lines, with a measured correlation coefficient of ${\rho }_{L_{[\text{CI}](\text{1-0})}^{\prime }-L_{\text{AGN}}}=0.27$ (p-value = 0.3). Interestingly, when we include the LIRGs in the analysis (see Fig. A.3), [CI](1–0) shows a stronger correlation with L_AGN than CO(1–0) or CO(2–1), we measure a correlation parameter of ${\rho }_{L_{[\text{CI}](\text{1-0})}^{\prime }-L_{\text{AGN}}}=0.68$ (p-value = 3.7 × 10⁻³), with a slope equal to 0.74 ± 0.26, a similarly tight relation to the one found for $L_{[\text{CI}](\text{1-0})}^{\prime }$ and SFR for the entire sample. Indeed, our sample poorly populates the L_IR < 10¹² L_⊙ regime, and these results are solely driven by the galaxy IRAS F12243-0036 (also known as NGC 4418), a source that is an outlier in many respects. This galaxy has also the lowest redshift (z = 0.00708) in our sample, and was discovered by Sakamoto et al. (2010) to host an extremely compact obscured nucleus (CON). Here, ∼ 10⁸ M_⊙ of molecular gas are concentrated in a region with < 20 pc size, pointing to extremely high column densities of N_H > 10²⁵cm⁻² (see Falstad et al. 2021, and references therein). If such a correlation between $L_{[\text{CI}](\text{1-0})}^{\prime }$ and L_AGN is confirmed with larger statistics, this result could suggest different behaviors of CO and [CI] as H₂ gas tracers in AGN hosts.

5.4. Galaxy-integrated CO line ratios

The distributions of galaxy-integrated CO line ratios ($r_{21}\equiv L_{\text{CO}(\text{2-1})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$, $r_{31}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$, $r_{32}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{2-1})}^{\prime }$) for our local ULIRGs are shown in Fig. 6. Our sample spans a very wide range of values: 0.6 ≲ r₂₁ ≲ 1.9, 0.4 ≲ r₃₁ ≲ 1.3, and 0.4 ≲ r₃₂ ≲ 1.7. For comparison, we report in Fig. 6 also the measurements obtained by Leroy et al. (2022) across tens of local disk galaxies (N ∼ 40, 30, 20 objects with r₂₁, r₃₂, and r₃₁ measurements, respectively), which are based on resolved CO maps and are thus unaffected by beam mismatch. We note that the Leroy et al. (2022) sample is only representative of local massive (log₁₀M_*[M_⊙] ∼ 10.2 − 10.8), star-forming galaxies on the main sequence (− 10.5 ≲ log₁₀sSFR [yr⁻¹] ≲ −9.9).

Fig. 6. Distribution of galaxy-integrated CO line ratios obtained for our sample of ULIRGs, with mean and median values reported at the top right of each plot. Left: $r_{21}\equiv L_{\text{CO}(\text{2-1})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$. Middle: $r_{31}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$. Right: $r_{32}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{2-1})}^{\prime }$. The mean CO line ratios measured by Leroy et al. (2022) in local massive main-sequence disk galaxies are shown using vertical dotted lines, with the gray shaded regions representing their 16-84 velocity percentile ranges.

Our mean and median r₂₁ values, respectively 1.12 and 1.09, are significantly higher than the average ⟨r₂₁⟩ = 0.79 ± 0.03 measured by Saintonge et al. (2017) for the xCOLD GASS sample⁹. Although dominated by massive main-sequence galaxies, the xCOLD GASS sample, being only M_*-selected, includes targets above the main-sequence as well as AGNs; hence, it is not necessarily representative of purely star-forming disks. Indeed, the molecular ISM of main-sequence disks generally presents even lower values of r₂₁ ≈ 0.6 − 0.7 at z ∼ 0 (den Brok et al. 2021; Leroy et al. 2022), and perhaps also at higher redshift (see the results at z ∼ 2 by Aravena et al. 2014). Whether the r₂₁ ratio is higher or not in AGNs is not clear, indeed the mean r₂₁ measured by Husemann et al. (2017) and Shangguan et al. (2020b) in nearby unobscured quasars are ⟨r₂₁⟩ = 0.9 ± 0.3 and $⟨r_{21}⟩=0.6_{-0.13}^{+0.15}$, respectively. It is however quite well established that the central regions of star-forming galaxies, and in general regions with higher Σ_SFR, present systematically higher ratios than the outskirts (den Brok et al. 2021; Leroy et al. 2022).

Despite slight variations in different samples, global r₂₁ luminosity ratios above unity are extremely rare in the Universe, yet they are predominant in our sample of local (U)LIRGs. This finding is in agreement with Papadopoulos et al. (2012) who, through their analysis of 70 (U)LIRGs with heterogeneous single-dish multi-J CO observations, measured ⟨r₂₁⟩ = 0.91. This is slightly lower than our mean value, probably because of their larger statistics. It is also probable that the smaller apertures of the James Clerk Maxwell Telescope (JCMT) and IRAM 30m telescopes used for the CO(2–1) observations, combined with the narrower bandwidths of the heterodyne receivers used at that time, have contributed to some CO(2–1) flux loss for the most extreme (U)LIRGs of the Papadopoulos et al. (2012) sample. In any case, the r₂₁ value distribution obtained by these authors is similarly broad to the one we obtain (Fig. 6) with many sources globally characterized by r₂₁ > 1 values.

The offset in global CO ratios from the local galaxy population is even more extreme in the r₃₁ values, whose mean and median of 0.81 and 0.76 measured in our sample of ULIRGs are significantly higher than ⟨r₃₁⟩ = 0.31, which is the mean galaxy-integrated value measured by Leroy et al. (2022). Strikingly, Fig. 6 shows that there is no overlap between our ULIRGs sample and normal massive star-forming galaxies in the r₃₁ values. The results by Leroy et al. (2022) are consistent with literature data focusing on local massive main-sequence star-forming galaxies, overall confirming that the low-J CO line emission in these sources is dominated by optically thick clouds with moderately sub-thermally excited CO(2–1) and CO(3–2) transitions. As expected due to the diversity of their sample, Lamperti et al. (2020) measure higher values of r₃₁ compared to Leroy et al. (2022). Specifically, Lamperti et al. (2020) obtain ⟨r₃₁⟩ = 0.55 ± 0.05 for 25 xCOLD GASS objects with JCMT CO(3–2) data, which probe a massive (M_* > 10¹⁰ M_⊙), highly star-forming (with log₁₀sSFR [yr⁻¹] > −10.5 and with most sources being above the main sequence) portion of the parent xCOLD GASS sample presented by Saintonge et al. (2017). The same study by Lamperti et al. (2020) includes also 36 hard X-ray-selected AGNs from the BASS survey with JCMT CO(3–2) and CO(2–1) observations, and, by extrapolating the CO(1–0) luminosities from the CO(2–1) data, they report consistent r₃₁ values compared to a non-AGN sample with matched specific SFR. The average LIRG value obtained by Papadopoulos et al. (2012) is ⟨r₃₁⟩ = 0.67¹⁰, with a distribution displaying a significant tail including many sources with r₃₁ ≳ 1.

For completeness, we show in Fig. 6 also our r₃₂ measurements, with average and median respectively r₃₂ = 0.80 and r₃₂ = 0.76. These are also quite offset from normal galaxy disks, for which Leroy et al. (2022) report ⟨r₃₂⟩ = 0.50.

5.5. Global CO line ratios as a function of galaxy properties

In Fig. 7 we plotted the total CO line ratios as a function of different galaxy properties, namely: L_IR (top left), SFR (top middle), L_AGN (top right), M_mol (bottom left), molecular gas depletion timescale due to star formation, defined as τ_dep ≡ M_mol/SFR, which is the inverse of the SFE (SFE ≡ SFR/M_mol) (bottom middle), and the AGN fraction (α_AGN ≡ L_AGN/L_bol). The values of L_IR, SFR, L_AGN, and α_AGN are taken from the literature and are listed in Table 1 together with their corresponding references. Instead, the M_mol and τ_dep values have been computed in this work using, for each source, the lowest-J CO transition available and adopting the median α_CO factor computed in this work, that is, 1.7 M_⊙ (K km s⁻¹pc²)⁻¹ (see Sect. 5.2).

Fig. 7. CO line ratios as a function of galaxy properties: LIR (top left), SFR (top middle), molecular gas depletion timescale due to star formation τdep ≡ Mmol/SFR (top right), molecular gas mass Mmol (computed by using the average αCO derived in this work; bottom right), LAGN (bottom middle), and AGN fraction (LAGN/Lbol; bottom right). We show the binned values using darker colors, overplotted on the individual data points. Dashed lines indicate the best-fit relations obtained from a least squares regression analysis conducted for each transition line separately. The color coding refers to different line ratios, as indicated in the legend at the top-right corner of each plot.

With the aim of better highlighting any underlying trends, we divided the ULIRG sample into bins according to the quantity on the x-axis and calculated the mean values of the ratios in these bins, which are shown using darker symbols in Fig. 7. The dashed lines show the best-fit relations, with color coding indicating different CO line ratios. An inspection of Fig. 7 shows the absence of strong correlations between the global low-J CO line ratios and the other quantities investigated here. This is not completely surprising given the narrow dynamic range in galaxy properties spanned by our sample, which is representative of the most extremely infrared-bright galaxies of the local Universe (log₁₀L_IR[L_⊙] ∈ [12.0, 12.8], making it a very special sample in its galaxy and ISM properties, already evident in Fig. 6, where almost no overlap with the normal population of star-forming galaxies is seen. However, we retrieve (although still weak) positive correlations between some of the low-J CO line ratios and quantities related to the strength of the star formation activity (namely, L_IR, SFR and SFE$={\tau }_{\mathrm{dep}}^{-1}$).

In the top-left panel of Fig. 7 we investigate the relations between the line ratios and L_IR. With measured Pearson correlation coefficients and p-values of ρ_{r21 − LIR} = 0.50 (p-value = 0.02), ρ_{r31 − LIR} = 0.62 (p-value = 0.008) and ρ_{r32 − LIR} = −0.02 (p-value = 0.93), we find that r₂₁ and r₃₁ show a positive relation with the infrared luminosity, while r₃₂ does not. In the relations with SFR (middle top panel), we find the ratio r₃₁ to be the only one showing a positive correlation (ρ_{r31 − SFR} = 0.69 and p-value = 0.002), while no significant correlation is found for r₂₁ (ρ_{r21 − SFR} = 0.40 and p-value = 0.09), and r₃₂ (ρ_{r32 − SFR} = 0.10 and p-value = 0.61). Equivalently to the results seen in Fig. 5, the relations with L_IR and SFR are expected to be similar since much of the infrared luminosity in these galaxies is due to star formation. The correlations found with L_IR are in agreement with the study performed by Rosenberg et al. (2015) on the HERCULES sample¹¹, where a correlation was found between high-J CO ratios and L_FIR (as well as with dust color). Lamperti et al. (2020) and Leroy et al. (2022) found, for massive main-sequence galaxies with or without an AGN, positive correlations between CO line ratios and SFR, in particular for the ratio r₃₁. If the CO(1–0) emission traces the total H₂ reservoir including the more diffuse components, and the CO(3–2) emission traces the somewhat already denser gas, then the ratio r₃₁ can be interpreted as a measure of the fraction of molecular gas that is in the denser star-forming regions. This would lead to higher values of r₃₁ for ULIRGs and SB galaxies, as these sources are expected to have higher fractions of dense gas, leading also to an increased SFE (and decreased τ_dep). Indeed, in our sample of ULIRGs, we compute Pearson coefficient for the CO ratios versus τ_dep of ρ_{r21 − τdep} = −0.27 (p-value = 0.27), ρ_{r31 − τdep} = −0.58 (p-value = 0.02) and ρ_{r32 − τdep} = −0.26 (p-value = 0.17), where only r₃₁ shows a (negative) correlation with τ_dep, consistent with the (opposite) trend seen by Lamperti et al. (2020) between r₃₁ and SFE.

We do not find any trend between these global ratios and molecular gas mass estimates (ρ = 0.19, ρ = 0.25 and ρ = −0.24, for r₂₁, r₃₁, and r₃₂, all having p-values > 0.05). This result is in agreement with Yao et al. (2003), while Leroy et al. (2022) find a weak anticorrelation between r₂₁ or r₃₁ and the lowest-J CO transition available ($L_{\mathrm{CO},\mathrm{low}}^{\prime }$). We also do not find any significant trend between CO line ratios and L_AGN (ρ = 0.29, ρ = 0.10 and ρ = −0.10, for r₂₁, r₃₁, and r₃₂, all having p-values > 0.05), consistent with Lamperti et al. (2020) and Yao et al. (2003). The relations between CO line ratios and the AGN fraction α_AGN (lower-right panel of Fig. 7), may suggest a negative trend, in particular for r₃₁, though still not statistically significant as measured by ρ_{r31 − αAGN}=-0.43 (p-value = 0.09), while r₂₁ and r₃₂ show no correlation at all with α_AGN, having ρ_{r21 − αAGN} = −0.03 and ρ_{r32 − αAGN} = −0.2 and p-values > 0.05.

5.6. Global [CI](1–0)/CO(1–0) line ratios

The distribution of integrated [CI](1–0)/CO(1–0) line luminosity ratios ($r_{\text{CICO}}\equiv L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$) obtained for our sample is reported in Fig. 8. The measurements span the range 0.08 ≲ r_CICO ≲ 0.4, with a median of ⟨r_CICO⟩^median = 0.18 and a mean of ⟨r_CICO⟩^mean = 0.21.

Fig. 8. Distribution of galaxy-integrated $r_{\text{CICO}}\equiv L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$ line ratios obtained for our sample of ULIRGs, with mean and median values reported on the plot.

Due to the paucity of [CI](1–0) detections available in the literature, it is unfortunately not possible to compare our results with statistically significant measurements performed on less extreme samples of local star-forming galaxies. Furthermore, often local targets observed in [CI](1–0) do not have adequate, aperture-matched CO(1–0) line data that can be used to compute the [CI](1–0)/CO(1–0) ratio, so there are only a few works to which we can compare our results. Jiao et al. (2019) analyze Herschel SPIRE [CI](1–0) maps (with 1 kpc resolution) of 15 nearby galaxies, and combined them with single-dish CO(1–0) observations, convolved to the same beam, obtained with the single-dish Nobeyama 45m telescope. The 15 sources of Jiao et al. (2019) are very famous nearby galaxies: a few SBs (such as M 82, NGC 253, and M 83), one HII galaxy (NGC 891), six low-ionization nuclear emission-line region galaxies (LINERs) and three Seyferts (including NGC 1068). Their sample is small but it is diverse and covers almost three orders of magnitude in SFR surface density. Therefore, their median r_CICO of 0.11 (mean is 0.12), based on multiple measurements per galaxy (one per resolution element), may be considered the closest we can get to a value that is representative of typical local star-forming galaxies. Michiyama et al. (2021) perform simultaneous CO(4–3) and [CI](1–0) observations with ACA Band 8 in 36 local (U)LIRGs, and report single-dish CO(1–0) line luminosity measurements for 25 of their targets that have also an [CI](1–0) line detection. We used these values to compute the average and median r_CICO ratio for the Michiyama et al. (2021) sample with both CO(1–0) and [CI](1–0) measurements, both coincidentally being 0.13 with a standard deviation of 0.07. Their median value of 0.13 is lower than what we find for our sample of ULIRGs, despite our samples overlapping by seven sources (of which only five have a CO(1–0) value in Michiyama et al. 2021) and span a similar range of SFRs. We note that only in two cases we decided to employ in our analysis the ACA Band 8 data collected by Michiyama et al. (2021; project ID 2018.1.00994.S; see Table B.1). For the other five overlapping targets we preferred our own APEX PI observations over the ACA archival data, because of their better quality and/or higher flux recovered (see the duplicated [CI](1–0) data shown in Fig. C.1). This is consistent with the considerations by Michiyama et al. (2021), who estimate an [CI](1–0) line flux loss of ∼30 − 40 % for their ACA data.

We then investigated the r_CICO ratio as a function of the following galaxy properties: L_IR, SFR, L_AGN, M_mol, τ_dep, and α_AGN. According to the Pearson correlation coefficients, none of these relations show a significant trend. Hence, in the main body of the paper, we show only the r_CICO versus L_IR relation (Fig. 9) and report the other plots in Fig. E.1. As explained earlier, it is hard to make a comparison between our results and previous studies of [CI] and CO in galaxies, because these sparse previous analyses have mostly considered very heterogeneous samples and data sets, with large aperture and sensitivity variations between different data. With these caveats in mind, we can discuss Fig. 9 in relation to some previous findings. A study performed by Israel et al. (2015) in a sample of 76 local galaxies with Herschel [CI] line data and ground-based JCMT ¹³CO(2–1) line spectra, shows a correlation between the [CI](1–0)/¹³CO(2–1) flux ratio and L_FIR. Though with a Pearson parameter indicating a non-significant correlation, Fig. 9 shows a hint of a positive trend between the [CI](1–0)/¹²CO(1–0) line ratio and infrared luminosity in our sample of local ULIRGs. There are however two important things to consider. Firstly, the [CI](1–0)/¹³CO(2–1) ratio investigated by Israel et al. (2015) may behave differently from the [CI](1–0)/¹²CO(1–0) line ratio that we explore here. Indeed the different isotopologue and rotational level of the CO transition used to compute the ratio implies a different excitation temperature, abundance, and opacity. Secondly, the study by Israel et al. (2015) includes sources spanning L_IR ∼ 10⁹ L_⊙ to ≳10¹² L_⊙, while all the data points studied in Fig. 9 are in the L_IR > 10¹² L_⊙ regime. Performing a visual inspection of the Israel et al. (2015) [CI](1–0)/¹³CO(2–1) – L_FIR relation, considering only the data points at L_FIR > 10¹² L_⊙, we observe no clear trend. Taking these considerations into account, our results are actually in agreement with these previous findings: our average r_CICO is clearly higher than that computed in less infrared-luminous galaxies by Jiao et al. (2019), as discussed above for Fig. 8. Another study we can examine is Valentino et al. (2018). These authors investigate the [CI](1–0)/CO(2–1) luminosity ratio in a heterogeneous sample including star-forming galaxies on the z ∼ 1.2 main-sequence observed with ALMA, and ∼30 local SBs (mostly LIRGs and AGNs) with archival Herschel [CI] and ground-based single-dish low-J CO data, respectively from Liu et al. (2015) and Kamenetzky et al. (2016). Their total sample spans 10¹⁰ ≲ L_IR [L_⊙] ≲ 10^13.5; they find no [CI]/CO ratio variations as a function of L_IR, with an average value of $L_{[\mathrm{CI}](1-0)}^{\prime }/L_{\mathrm{CO}(2-1)}^{\prime }\simeq 0.20\pm 0.02$, consistent with our results.

Fig. 9. rCICO as a function of LIR for our sample of ULIRGs. Similar relations with other galaxy properties are reported in Fig. E.1. Sources with direct measurements of both the [CI](1–0) and CO(1–0) lines are plotted using purple crosses. The pink diamonds represent sources observed in [CI](1–0) but without CO(1–0) line data, for which we inferred $L_{\text{CO}(\text{1-0})}^{\prime }$ from the CO(2–1) line luminosity by assuming r21 = 1.1, which is the mean value computed for our sample (see Fig. 6). The dashed purple line represents the best fit using a least squares regression analysis, and the shaded gray area corresponds to the 1σ confidence interval of the fit. The Pearson correlation coefficients (ρ) are reported at the upper-left corner of the panel with their associated p-values.

Similarly, for SFR and L_AGN, we retrieve Pearson correlation coefficients of ρ_{rCICO − SFR} = 0.31 and ρ_{rCICO − LAGN} = 0.39 (with p-values of 0.28 and 0.17, respectively).

6. Discussion

The discussion is divided into three parts. First, in Sect. 6.1, we discuss the possible origin of the high CO line ratios measured in our sample, by exploring also possible dependences on gas kinematics (ISM turbulence, and presence of outflows). Then, in Sect. 6.2 we discuss the origin of the high r_CICO ratios, and investigate their relationship with prominent molecular outflows. Finally, in Sect. 6.3 we discuss the role of AGN feedback in setting the molecular gas properties of (U)LIRGs.

6.1. Physical drivers of low-J CO line ratios in (U)LIRGs

6.1.1. A weak dependence on ISM excitation

Our sample is characterized by extremely high global r₂₁, r₃₁, and r₃₂ values, which show little or no overlap with those measured in the local main-sequence galaxy population (see Fig. 6). Only marginal correlations were found, mainly for r₃₁ with L_IR, SFR, and τ_dep, while no statistical trends were found with L_AGN, or α_AGN (see Fig. 7). This suggests that, in these (U)LIRGs, galaxy-averaged low-J CO line ratios are unable to trace effects originated by AGN activity, while ratios depending on CO(3–2) line emission may start to reflect effects of star formation activity and gas density (seen in the relations with SFR and τ_dep). Overall, our results agree with the conclusion of Papadopoulos et al. (2012), who stated that the low-excitation J_up ≤ 3 portion of CO spectral line energy distributions (SLEDs) are “highly degenerate to the average state of the ISM,” as we verified here for these (U)LIRGs that are not characterized by cold CO SLEDs.

We highlight the observed large variations of CO line ratios within our sample, much larger than those measured in main-sequence galaxies. Does Fig. 7 inform us about the drivers of such strong differences? The largest r_{J, J − 1} variations appear to be typical at the high-L_IR end of the sample, in particular for the global r₂₁ values, which become consistently r₂₁ > 1 at log₁₀L_IR [L_⊙] ≃ 12.3. However, Fig. 7 shows only a weak a link between large CO ratio variations and SFR, but no link to AGNs. This may seem puzzling since AGNs are strong sources of far-UV photons and CRs, which can excite CO, but it is consistent with the hypothesis that low-J CO line ratios in these (U)LIRGs are actually degenerate for CO excitation effects. Higher-J CO transitions are needed to better constrain the effects of AGN on the excitation of the ISM (Gallerani et al. 2014; Rosenberg et al. 2015; Jarvis et al. 2020).

Ratios > 1.0 can only be obtained through highly excited gas coupled with low optical depths. More specifically, CO lines can become partially transparent (i.e., reduced line opacities) in the presence of a diffuse, warm, turbulent component characterized by a large velocity gradient. Cicone et al. (2018), based on the analysis of NGC 6240, suggest that such a turbulent (and/or envelope in) H₂ phase may be typical of massive molecular outflows, which are widespread in (U)LIRGs. In Sect. 6.1.2 we explore the hypothesis of a connection between high CO line ratios and the gas kinematics.

6.1.2. Effects of high-velocity and turbulent gas on the CO line ratios

The simultaneous multi-Gaussian component spectral fitting procedure described in Sect. 4.1 (see the spectral fits in Figs. B.1–B.6) allowed us to study CO line ratios as a function of the kinematic properties (v_cen and σ_v) of the CO-emitting gas. The results of this investigation are shown in Fig. 10, where we plot the r₂₁, r₃₁, and r₃₂ obtained for each Gaussian component employed in the fit as a function of their kinematic parameters, namely: σ_v, |v_cen|, and the combined parameter σ_v + |v_cen|. We also studied the dependence on v_cen, and the results were similar to the ones shown in Fig. 10, which is why we chose to only show the absolute value.

Fig. 10. Individual spectral component analysis of the CO line ratios as a function of the gas kinematics for our sample of ULIRGs. We set a detection threshold of 2σ for the flux of each Gaussian component in order to compute the ratios. The top, central, and bottom rows show, respectively, the r21, r31, and r32 plots. Left column: Line ratios as a function of the velocity dispersion, σv (line width of the spectral Gaussian component). Middle column: Line ratios as a function of the module of the central velocity of the spectral Gaussian component (|vcen|), measured with respect to the galaxy’s redshift, reported in Table 1. Right column: Line ratios as a function of the σv + |vcen| parameter. The dashed lines correspond to the best-fit relations obtained from a least squares regression analysis, and the shaded gray areas correspond to the 1σ confidence interval of the fit. The Pearson correlation coefficients and associated p-values are reported at the top-right corner of each plot, although they are less informative than the best fit because they do not take the (large) error bars into account.

Concerning the gas kinematics, the σ_v values range from ∼10 km s⁻¹ up to ∼230 km s⁻¹, and v_cen from ∼ − 300 km s⁻¹ up to ∼315 km s⁻¹. Most Gaussian components display |v_cen|≲100 km s⁻¹, which are most likely tracing “disk” material with systemic velocities, though we observe components reaching up to |v_cen|∼315 km s⁻¹ tracing blueshifted or redshifted velocity gas, with motions clearly deviating from ordered rotation. We caution against an over-interpretation of the specific best-fit kinematic properties for the individual targets, because the possibility to fit multiple components depends mainly on the S/N of the lines and on the presence of spectral substructures. As a result, a source with generally broad CO lines may be fitted either with one single broad component or with multiple narrow components at different velocities. We therefore stress again that the plots presented in Fig. 10 should be regarded for their statistical value.

The CO line ratios measured for the individual Gaussian components span similar ranges as the global galaxy-integrated ratios (see Fig. 6), namely: 0.7 ≤ r₂₁ ≤ 1.8, 0.3 ≤ r₃₁ ≤ 1.8, and 0.5 ≤ r₃₂ ≤ 1.6. The Pearson correlation coefficients (ρ) and associated p-values (see Fig. 10) do not evidence any significant trend among the nine relations explored. There are clearly instances where spectral components characterized by higher σ_v and/or higher central velocity show higher CO line ratios, but these trends are not statistically significant if considering the whole sample.

However, we note that the correlation coefficients do not take into account the error bars, which are quite large in this case. We therefore proceed with a regression analysis, whose results are over-plotted in Fig. 10. The r₂₁ versus σ_v plot (top-left panel) shows a statistical positive trend, described by the relation

$$\begin{array}{c}\hfill r_{21}=(0.76\pm 0.08)+(2.0\pm 0.7)\times {10}^{-3}{\sigma }_v.\end{array}$$(11)

A similar one is found between r₂₁ and σ_v + |v_cen| (r₂₁ = (0.79 ± 0.07)+(1.4 ± 0.5)×10⁻³(σ_v + |v_cen|)). Since no clear trend is detected between r₂₁ and |v_cen|, the relation found between r₂₁ and σ_v + |v_cen| is consistent with being entirely driven by the trend with σ_v.

Among the CO ratios explored here, r₂₁ is the least sensitive to ISM gas excitation, and so we can hypothesize that the scatter observed in r₂₁ values is dominated by optical depth effects (see also Zschaechner et al. 2018). In particular, high r₂₁ ≳ 1 can be explained by a low opacity, especially in the CO(1–0) transition. The positive correlation between r₂₁ and σ_v reported in Eq. (11) suggests that the low CO opacity is driven by large velocity gradients, which can in turn be due to turbulence in the ISM and/or bulk motions within molecular outflows. A similar positive trend between r₂₁ and σ_v was observed in NGC 6240, and in that case the high-r₂₁ and high-σ_v gas was predominantly linked to the massive molecular outflow (Cicone et al. 2018).

We note that neither the r₃₁ nor the r₃₂ values show significant trends as a function of gas kinematics. It is possible that ratios involving the CO(3–2) line transition are less sensitive to optical depth effects and more affected by gas excitation, and hence the absence of trends with these ratios may support an interpretation of the r₂₁ versus σ_v relation in terms of opacity of CO lines, rather than gas excitation effects.

As explained in Sect. 2, our sources were selected from the Herschel OH parent sample, and we made use of the literature OH119 μm data to study possible dependences between the line ratios measured in this paper and the available OH outflow parameters, mainly the OH maximal outflow velocity (OH_max) for sources that show absorption components and the equivalent width of the OH line (OH_EQW). The latter parameter is a negative value for sources where the absorption component is stronger than the emission component of the OH119 μm doublet. The OH parameters are also tracers of molecular gas kinematics in these galaxies and we can test whether there are correlations with the global CO line ratios measured in our sample. In Fig. 11 we show the global r₂₁, r₃₁ and r₃₂ line ratios as a function of the OH_max and OH_EQW. As we did also in Sect. 5.5, we divided the total sample into bins according to the quantity on the x-axis and calculated the mean values of the ratios in these bins, which are shown using darker symbols. The dashed lines show the best-fit relations, with color coding indicating different CO line ratios. We do not observe any statistically significant trend between the quantities, with Pearson correlation coefficients of ρ_{r21 − OHmax} = −0.33 (p-value = 0.23), ρ_{r31 − OHmax} = −0.13 (p-value = 0.68), ρ_{r32 − OHmax} = 0.10 (p-value = 0.64), ρ_{r21 − OHEQW} = −0.26 (p-value = 0.31), ρ_{r31 − OHEQW} = −0.32 (p-value = 0.24), ρ_{r32 − OHEQW} = 0.06 (p-value = 0.78). Additionally, the regression analysis performed on each luminosity ratio, which takes the error bars into account, reveal best-fit parameters that are consistent with flat trends for all the variables.

Fig. 11. Global CO line ratios as a function of OHmax (left) and OHEQW (right) for the sources with OH119 μm data as reported in the literature (see Table 1). The binned values (computed in the same way as in Sect. 5.4) are shown using darker colors and overplotted on the individual data points. Dashed lines indicate the best-fit relations obtained from a least squares regression analysis conducted for each transition line separately. The color coding refers to different line ratios, as indicated in the legend at the top-right corner of each plot.

Summarizing, we study any possible trends between the CO line ratios and the kinematics of the gas as traced by (i) the velocity of the individual spectral components of the CO multi-Gaussian fit, and (ii) the OH outflow properties. The former allows us to study the ratios for each spectral component, for which we retrieve a statistically significant positive trend (via regression analysis) between r₂₁ and σ_v, likely due to a lower CO opacity in the high-σ_v gas, which could be tracing outflows. On the other hand, we do not find significant trends between the global r₂₁ ratios and OH119 μm outflow properties, showing that stronger OH outflows are not statistically associated with globally enhanced r₂₁ values. These two results are only in apparent contradiction. Indeed, the OH119 μm outflows are detected in absorption against the FIR continuum and so trace mainly the inner components of the galaxy-scale molecular outflows, while the low-J CO lines, observed in emission and integrated over the whole extent of these galaxies, probe the totality of the molecular gas including extended and diffuse components far from the central engines. Indeed, it would be remarkable if we were to find a trend between the global r₂₁ values and the OH outflow properties. Higher-J CO lines, tracing the medium exposed to the excitation sources, may show tighter relations with the OH outflows. This hypothesis needs to be tested with J_up≥4 CO line observations. In conclusion, our analysis of CO line ratios as a function of gas kinematics, suggests that molecular outflows are not statistically characterized by higher CO line ratios in ULIRGs, at least for low-J CO lines up to J = 3. In other words, the low-J CO line ratios do not appear sensitive to the presence or strength of massive molecular outflows.

6.1.3. Absence of a relation between r₂₁ and α_CO

It has been proposed that, although the CO(2–1)/CO(1–0) luminosity ratio alone cannot be used to constrain the α_CO in galaxies (which require either independent H₂ tracers or a much more accurate multi-J analysis of CO lines), variations in r₂₁ should be correlated to those observed in α_CO, because both these parameters are affected by the optical depth of the gas. This holds as long as other factors, such as the CO abundance, are fixed. In particular, Gong et al. (2020) predict a close-to-linear anticorrelation between r₂₁ and α_CO (r₂₁ ∝ (1/α_CO)^1.2), based on their models tailored to the ISM conditions typical of galaxy disks. Such a trend is also strengthen by the findings in Leroy et al. (2022): these authors, by fitting radiative transfer models to measured line ratios, find an anticorrelation between r₂₁ and α_CO. However, such a dependence is not found in our sample. As shown by Fig. 12, r₂₁ and α_CO values are not correlated (ρ = 0.16, p-value = 0.66). The best-fit relation, which takes the large uncertainties on α_CO into account, is consistent with a flat trend.

Fig. 12. Global r21 ratios plotted as a function of [CI]-derived αCO values. The purple dashed line shows the best-fit relation, with its corresponding 1σ confidence interval represented by the gray shaded area. Neither the Pearson correlation coefficient (ρ) nor the regression analysis suggest any correlation between these quantities.

This result further strengthens the hypothesis that the molecular ISM of (U)LIRGs, as traced by low-J CO line ratios, is extremely different from that of normal galaxy disks. The global r₂₁ > 1 measured in our sample are not even allowed by the models proposed by Gong et al. (2020). Furthermore, even if such high r₂₁ values strongly suggest optically thin CO emission, the corresponding α_CO estimates of ∼ 1 − 4 M_⊙ (K km s⁻¹pc²)⁻¹ measured in those same sources that display r₂₁ > 1 (see Fig. 12), are significantly higher than the optically thin α_CO values measured in, for example, disk galaxies (see, e.g., Bolatto et al. 2013). This indicates an apparent contradiction that was already observed in NGC 6240 (Cicone et al. 2018). A possible explanation is the coexistence of a diffuse, turbulent molecular ISM phase, which is characterized by large velocity gradients (such as those typical of outflows) and dominates the low-J CO luminosity, with a denser phase, which dominates the gas mass and drives up the [CI]-based α_CO measurements. The latter phase may be preferentially traced by denser molecular gas tracers, such as CN, HCN, CS, and higher-J CO transitions.

6.2. Effects of high-velocity and turbulent gas on the [CI]/CO line ratio

In Sect. 5 we show that, in our targets, $L_{\mathrm{CO}(1-0)}^{\prime }$ (as well as $L_{\mathrm{CO}(2-1)}^{\prime }$) is tightly correlated with $L_{[\mathrm{CI}](1-0)}^{\prime }$, and that the [CI](1–0)/CO(1–0) luminosity ratios, spanning 0.07 ≲ r_CICO ≲ 0.4, are significantly higher than those measured in normal star-forming galaxies. Such high r_CICO values appear to be typical of infrared-bright galaxies, as suggested by previous works (Israel et al. 2015; Valentino et al. 2018). Within our sample, we do not find any significant relation between r_CICO and intrinsic galaxy properties such as SFR, M_mol, τ_dep, or L_AGN. In this section we test the hypothesis of a link between r_CICO and gas kinematics, and in particular with molecular outflows.

As described in Sect. 4.1, the [CI](1–0) line spectra were fitted independently from the CO lines, due to their generally lower S/N, and to the possibility that CO and [CI] lines may not have identical profiles (see also results in Sect. 5.1.2). Therefore, it is not possible to investigate r_CICO ratios as a function of kinematic parameters for individual Gaussian components as we did in Fig. 10 for the CO line ratios. However, we can still investigate the global r_CICO values as a function of velocity dispersion of the gas, as traced by, for example, σ_v (from the single Gaussian fit to [CI](1–0) and CO(2–1); see Sect. 5.1.2), the v_97.7 − v_2.3 interval for CO and [CI], and as a function of OH properties (i.e., maximal outflow velocity and equivalent width). In Fig. 13, we show the relations between r_CICO and σ_v, CO (top-left panel), v_97.7 − v_2.3 for CO (top-right panel), OH outflow maximal velocity (bottom-left panel) and OH equivalent width (bottom-right panel). All Pearson correlation coefficients and the respective p-values are reported in the figures, evidencing no indication of a correlation between any of the quantities. The regression analyses, which take the error bars into account, are consistent with flat trends in all cases. Similarly, we do not find any trend between r_CICO and [CI] line widths, as probed by σ_v, [CI] and its v_97.7 − v_2.3 velocity interval. Therefore, our results do not evidence a significant impact of the gas kinematics on the global [CI]/CO luminosity ratio in (U)LIRGs.

Fig. 13. Global rCICO as a function of the velocity dispersion of the gas. Top row: [CI](1–0)/CO(1–0) luminosity ratio as a function of σv, CO from the single Gaussian fit (see Sect. 5.1.2, left), and as a function of v97.7 − v2.3 velocity interval for CO (right). Bottom row: rCICO line ratio as a function of OHmax (left) and OHEQW (right). Sources with direct measurements of both the [CI](1–0) and CO(1–0) lines are plotted using purple crosses. For the other sources (pink diamonds), the CO(1–0) line luminosity was estimated assuming the average CO(2–1)/CO(1–0) luminosity ratio of r21 = 1.1 computed for our sample of ULIRGs. The dashed purple line represents the best fit using a least squares regression analysis, and the shaded gray area corresponds to the 1σ confidence interval of the fit. The Pearson correlation coefficients (ρ) are reported at the upper left corner of the panel with their associated p-values.

In light of this independence of r_CICO from gas kinematics, it remains difficult to explain the detection, shown in Fig. 2, of a marginal deviation of the [CI] and CO line widths from a 1:1 relation, where [CI] lines appear narrower than CO lines in sources with broad CO lines. One hypothesis is that a slight depletion of [CI] with respect to CO in the broad wings of the lines conspires with an enhancement of [CI]/CO in the narrow core of the same lines, so that σ_v, [CI] < σ_v, CO while the total luminosity ratio, r_CICO, remains unchanged. The reason why [CI] may be fainter in the high-velocity wings of the molecular lines could be due to its optically thin nature. This makes it challenging to detect this line in the presence of more diffuse and extended components of the medium, where the column density is much lower than in the disks. On the contrary, low-J CO lines such as CO(1–0) and CO(2–1) remain bright and relatively easy to detect in such components. The apparent bi-modality of the data observed in the two top panels of Fig. 13, where a cluster of sources with high r_CICO and narrower CO lines populate the upper-left regions of the plots, while the remaining data points show low r_CICO and high CO line widths, could strengthen that hypothesis. We warn the reader, however, that this needs to be confirmed with larger statistics.

Moreover, the premise of the different behaviors for the two gas tracers studied here is not supported by any theoretical prediction. Actually, such an effect would be at odds with the theoretical prediction that massive molecular outflows contain a significant CO-dark H₂ component, due to leakage of far-UV photons and CRs (Papadopoulos et al. 2004, 2018), which instead seems to be a promising avenue to explain the observational results by Saito et al. (2022) and Ueda et al. (2022). In particular, Saito et al. (2022) use ALMA to explore the [CI](1–0) and CO(1–0) emission within the central 1 kpc of the Seyfert 2 galaxy NGC 1068. They measure the highest r_CICO values (median 0.72, 16th–84th range of 0.33–1.86) in a region offset by ∼300 pc from the AGN position, beyond the circumnuclear disk, which appears to correspond to the shocked environment between the molecular outflow (see García-Burillo et al. 2014) and the galactic disk. These authors interpret the high [CI]/CO ratios as arising from the shocked gas, due to efficient shock dissociation of H₂ gas by the shock wave, and the consequent dissociation of CO via atomic hydrogen endothermically reacting with said molecule (see Hollenbach & McKee 1980). On the contrary, the SB ring of NGC 1068, which is unaffected by the outflow, displays a much lower median r_CICO = 0.15 (16th–84th ranges of 0.08–0.23). Given the limited field of view of the [CI](1–0) observation by Saito et al. (2022), an r_CICO map beyond the central 1 kpc, which could probe the outflow in regions not overlapping with the galactic disk (therefore not shocked), is not available.

Another example is the study of Arp 220 performed by Ueda et al. (2022), who obtain an r_CICO map of the central 5″ (1.9 kpc) of the galaxy. These authors, however, report that the archival ALMA data used in their analysis suffer from missing flux, whose fraction they estimate to be 34% for the [CI](1–0) line through a comparison between their data and JCMT observations (Papadopoulos et al. 2004). Therefore, the results presented by Ueda et al. (2022) apply only for components smaller than the maximum recoverable scale (MRS) of their data sets, 1.75 kpc. Ueda et al. (2022) measure $r_{\mathrm{CICO}}^{redOF}=0.9\pm 0.3$ in the collimated molecular outflow hosted by the western nucleus of Arp 220, clearly higher than the galaxy-integrated value of r_CICO = 0.22 ± 0.04. They interpret such a high r_CICO value to be arising from regions with elevated [CI]/CO abundance ratio, which may be caused either by CRs or by shocks. Compared to the results obtained by Saito et al. (2022), the regions with high r_CICO in Arp 220 are well within its central region and, as such, could also be driven by shock waves, leading to CO dissociation, as suggested for NGC 1068. However, disentangling between CRs, high radiation fields or shocks, is challenging. Additionally, Ueda et al. (2022) also study how r_CICO relates to r₃₁ and find a positive trend between r_CICO > 0.3 measurements and high r₃₁ > 1 values. Hence, regions with high [CI]/CO luminosity show very high, optically thin CO(3–2)/CO(1–0) ratios. In our sample, although most of the sources with r₃₁ ≳ 0.9 have also r_CICO > 0.2, the correlation found between the two quantities is not statistically significant (ρ = 0.52, p-value = 0.15, best-fit relation: r_CICO = ( − 0.05 ± 0.14)+(0.24 ± 0.20) r₃₁). However, any kiloparsec-scale effects of the type observed by Saito et al. (2022) and Ueda et al. (2022), possibly triggered by collimated outflows and localized shocks, would probably be washed out in the galaxy-integrated measurements used in this work.

6.3. Role of AGN feedback

The CO and [CI]/CO line ratios investigated in this work do not show any significant correlation with L_AGN or α_AGN, suggesting that AGN radiation has little effect on these ratios, in this sample of (U)LIRGs.

AGN feedback could manifest in different ways in galaxies. One manifestation is through powerful outflows, which expel gas from the central regions and may deplete the galaxy of its M_mol reservoir (Cicone et al. 2014). Alternatively, AGN feedback could inject turbulence in the gas (through, e.g., compact jets or outflows), making star formation inefficient even in the presence of cold and dense H₂ gas (see, e.g., Alatalo et al. 2015).

We exploited our refined total molecular gas mass estimates, which are based on higher quality data as well as on an independent check on the α_CO factor, to investigate any links between the AGN luminosity and signatures of suppressed star formation in these (U)LIRGs. Figure 14 shows the molecular gas depletion timescale as a function of L_AGN. We note that this is the depletion timescale due to the consumption of gas by the star formation activity, without including the gas expelled via outflows (see, e.g., Cicone et al. 2014 for an investigation including AGN feedback).

Fig. 14. Depletion timescale (τdep ≡ Mmol/SFR) plotted as a function of LAGN. The dashed line shows the best-fit relation and its 1σ confidence interval plotted as a shaded gray area. The best-fit parameters and the Pearson correlation coefficient are reported at the upper left corner of the panel.

We find a weak though still statistically significant correlation (ρ = 0.33, p-value = 0.05), implying that the ULIRGs with more luminous AGNs consume their gas slower than the ones with less luminous AGNs, or, in other words, have a less efficient star formation. We note, however, that the L_AGN values, taken from Veilleux et al. (2013) and Spoon et al. (2013), were computed using method 6 by Veilleux et al. (2009b), which is based on the 15 to 30 μm continuum ratio (f₃₀/f₁₅). This method carries uncertainties on the order of ∼20% on average for the reported L_AGN. Moreover, we also explored the τ_dep versus OH outflow velocity relation, and we did not find any significant trend.

7. Summary and conclusions

We have performed a combined analysis of high S/N CO(1–0), CO(2–1), CO(3–2), and [CI](1–0) emission line spectra of 36 ULIRGs and an additional 4 LIRGs at z < 0.2, selected because they had OH119 μm observations from Herschel. These data sets are particularly sensitive to the total line emission from the sources, including faint, low-surface-brightness components that may be missed by high-resolution interferometric observations. All data, both PI and archival, were reduced and analyzed in a uniform way, with the goal of minimizing the aperture biases that usually affect the ratios computed between lines observed at different frequencies. The LIRGs, which are a minority in our sample, were discussed separately from the ULIRGs to avoid the low statistics at L_IR < 10¹² L_⊙ biasing our relations. Our findings can be summarized as follows:

• We find a tight relation between CO and [CI] luminosity, ${log}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(4\pm 2)+(0.64\pm 0.26){log}_{10}{L}_{[\text{CI}](\text{1-0})}^{\prime }$, and a very similar one for CO(2–1). This strengthens the hypothesis that atomic carbon is a valid alternative to CO for tracing the bulk of H₂ gas in galaxies.

• We compared the line widths of [CI] and CO, and derived a best-fit relation of σ_v, [CI] = (30 ± 9)+(0.64 ± 0.07) σ_v, CO when comparing the velocity dispersion obtained via a single Gaussian fit. We also explored the line widths computed at the wings of the lines, using the difference in percentile velocities. We obtain v_97.7 − v_2.3(CO) = (74 ± 66)+(0.70 ± 0.12)v_97.7 − v_2.3([CI]). Hence, in our sample, [CI] lines are on average narrower than CO lines, especially for sources that show broader line profiles.

• We combined [CI](1–0) and CO(1–0) observations to derive a galaxy-averaged α_CO factor. Based on the ten sources with data available for both emission lines, we computed a median value of ⟨α_CO⟩ = 1.7 ± 0.5 M_⊙ (K km s⁻¹pc²)⁻¹. This value is higher than the α_CO commonly used in the literature for (U)LIRGs and is consistent with recent estimates by Herrero-Illana et al. (2019) and Kawana et al. (2022) based on the thermal dust continuum.

• In the ULIRGs sample, the total CO and [CI] line luminosities show weak correlations with L_IR and the SFR. However, when including the LIRGs, the correlations become much stronger, though they remain weaker for SFRs than for L_IR. Our sample of (U)LIRGs is characterized by $L_{\text{CO}}^{\prime }/L_{\text{IR}}$ ratios that are similar to or lower than the SBs’ best fit derived by Sargent et al. (2014), and so it is clearly not representative of the general local star-forming population, but rather only its most extreme (merger-driven) SBs.

• We measure extremely high galaxy-integrated CO line ratios, significantly offset from those measured in local massive main-sequence galaxies. For example, more than 50% of our targets have global r₂₁ > 1 (⟨r₂₁⟩^median = 1.09). The distribution of r₃₁ values in our (U)LIRGs, with a median of ⟨r₃₁⟩^median = 0.76, does not even overlap with that of massive local galaxies. The extremely high low-J CO ratios of local (U)LIRGs were noted by several previous studies (e.g., Papadopoulos et al. 2012), and here we confirm this trend with better-quality data.

• We investigated possible drivers of such high CO line ratios. The global r₂₁ and r₃₁ values show a positive relation with L_IR, and r₃₁ also shows a correlation with the SFR. We studied links between CO line ratios and gas kinematics, as traced by the central velocity and σ_v of the spectral line components simultaneously fit to all three CO transitions. We find a positive trend between r₂₁ and the velocity dispersion of the spectral components, probably tracing CO opacity effects, which is described by the relation r₂₁ = (0.76 ± 0.08)+(2.0 ± 0.7)×10⁻³σ_v. From these results we infer that, in local (U)LIRGs, low-J CO line ratios are generally poor tracers of CO excitation, contrary to what has been found for local massive main-sequence galaxies (e.g., Leroy et al. 2022; Lamperti et al. 2020). We have also explored the CO line ratios as a function of OH outflow properties, as retrieved from the literature (Veilleux et al. 2013; Spoon et al. 2013), and conclude that these ratios are not sensitive to the presence or strength of molecular outflows.

• We measure $r_{\text{CICO}}=L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$ values between 0.07 and 0.4, with a median of ⟨r_CICO⟩^median = 0.18. The values we obtained, based on APEX observations, are significantly higher than previous (sparse) measurements in local galaxies, including those obtained from interferometric observations of a partially overlapping sample of ULIRGs. This highlights the importance of sensitive, single-dish observations that can capture the total line flux, especially at high frequencies, where the interferometric flux loss becomes more severe. The r_CICO ratios do not show any significant trend with the galaxy properties explored here.

• We do not observe any significant correlation between r_CICO and the velocity of molecular outflows (as probed by OH absorption) or gas kinematics (as probed by the CO or [CI] velocity dispersion, σ_v, and the v_97.7 − v_2.3 velocity intervals), and hence we rule out gas kinematics or turbulence having, statistically, any significant effect on the global [CI]/CO luminosity ratio. This result makes the discovery that [CI](1–0) lines are narrower than CO lines even more puzzling, because it implies that any effect causing [CI] lines to be narrower (in high-σ_v sources) conspires with some other effect in such a way that the total luminosity ratio remains unchanged.

• Finally, we used our refined total M_mol estimates to investigate links between L_AGN and signatures of suppressed star formation. We find a statistically significant correlation between τ_dep ≡ M_mol/SFR and L_AGN, implying that (U)LIRGs with more luminous AGNs are statistically less efficient at forming stars. This could result from injected turbulence by radio jets or AGN winds (see, e.g., Alatalo et al. 2015).

Overall, our analysis confirms the exceptional molecular gas properties of local (U)LIRGs, which are completely different from those of normal star-forming galaxies in the local Universe. This study highlights the necessity of gathering large samples of high S/N observations of alternative H₂ gas tracers in local galaxies, especially [CI] lines, high-J CO lines, and high-density gas tracers (e.g., HCN, CN, CS, and SiO), to understand the extreme mechanisms at work during accelerated phases of galaxy evolution. Due to the faintness of most of these tracers and their potentially different spatial extents in host galaxies, such high S/N multi-tracer analyses need to be conducted with a sensitive large-aperture single-dish telescope with a large instantaneous field of view. Furthermore, observations of [CI] lines and high-J CO lines require a high atmospheric transparency that can only be obtained in a high, dry site. Finally, the large line widths (σ_v) of these local galaxy mergers require stable spectral baselines. The new concept for the Atacama Large Aperture Submillimeter Telescope (AtLAST¹²) satisfies all of these conditions, and so it can be a transformational facility for this field.

Acknowledgments

We are very grateful to the APEX and ESO staff for their unceasing efforts at preserving APEX operations and completing our projects despite the Covid19 pandemic. We thank the referee for their constructive and thorough reports, which significantly helped us improve the paper. We thank our friend and collaborator Padelis P. Papadopoulos who provided, as usual, very useful feedback and expert input on this work, and who always warns us against superficial, mainstream interpretations. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 951815 (AtLAST). This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under programme IDs 0104.B-0672, 0106.B-0674, 086.F-9321, 090.B-0404, 092.F-9325, 099.F-9709, 077.F-9300, and 084.F-9306. APEX is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00177.S, ADS/JAO.ALMA#2015.1.01147.S, ADS/JAO.ALMA#2013.1.00535.S, ADS/JAO.ALMA#2016.1.00140.S, ADS/JAO.ALMA#2016.1.00177.S, ADS/JAO.ALMA#2015.1.00287.S, ADS/JAO.ALMA#2013.1.00180.S, ADS/JAO.ALMA#2018.1.00503.S, ADS/JAO.ALMA#2017.1.01235.S, ADS/JAO.ALMA#2016.2.00006.S, ADS/JAO.ALMA#2017.1.01398.S, ADS/JAO.ALMA#2013.1.00659.S, ADS/JAO.ALMA#2018.1.00699.S, ADS/JAO.ALMA#2012.1.00377.S, ADS/JAO.ALMA#2017.1.00297.S, ADS/JAO.ALMA#2015.1.00102.S, ADS/JAO.ALMA#2012.1.00611.S, ADS/JAO.ALMA#2016.2.00042.S, ADS/JAO.ALMA#2013.1.00180.S, ADS/JAO.ALMA#2018.1.00888.S, ADS/JAO.ALMA#2018.1.00994.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work is based on observations carried out with the IRAM Plateau de Bure Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). M.A. acknowledges support from FONDECYT grant 1211951, ANID+PCI+INSTITUTO MAX PLANCK DE ASTRONOMIA MPG 190030, ANID+PCI+REDES 190194 and ANID BASAL project FB210003. B.B. and S.S. acknowledge the support from the Research Council of Norway through NFR Young Research Talents Grant 276043. P.S. and C.C. acknowledge financial contributions from Bando Ricerca Fondamentale INAF 2022 Large Grant “Dual and binary supermassive black holes in the multi-messenger era: from galaxy mergers to gravitational waves” and from the agreement ASI-INAF n.2017-14-H.O.

References

Appendix A: LIRGs

This appendix shows several plots presented in the main body of the text, now including the 4 LIRGs in our sample. These sources populate poorly the 11.0 < log₁₀L_IR L_⊙ < 12.0 regime, with IRAS F12243-0036 (also known as NGC 4418), with the lowest infrared luminosity and redshift in our sample, being the main driver of the relations found for our extended sample of (U)LIRGs. In Fig. A.1 we plot the [CI](1–0) with respect to CO(1–0) and CO(2–1) line luminosities. The results are similar to the ones obtained in Fig. 1, although the slope is now steeper, with a close-to-linear relation between the quantities.

Fig. A.1. Same as Fig. 1, now including the two LIRGs (in purple data points) for which we have [CI](1–0) observations. The best-fit relations (shown as dashed purple lines) now include the LIRGs data points. The best-fit parameters are reported at the bottom-right corner of the plots. We also display the Pearson correlation coefficients (ρ) and their associated p-values. The solid black lines in both panels represent the corresponding relations reported by Jiao et al. (2017) for a sample of 71 (U)LIRGs, and the dotted lines represent the relations of Jiao et al. (2019) for a sample of 15 nearby spiral galaxies, between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$.

In Fig. A.2 we explore how the CO-to-H₂ conversion factor changes when taking into account the two LIRGs in our sample for which we have both [CI](1–0) and CO(1–0) data (i.e., IRAS F12243-0036 and IRAS F00509+1225). These two sources show average α_CO values of ∼2.5 M_⊙ (K km s⁻¹ pc²)⁻¹ and ∼0.7 M_⊙ (K km s⁻¹ pc²)⁻¹, respectively. The overall mean and median of the α_CO factor, remains the same for the extended sample of (U)LIRGs.

Fig. A.2. Same as Fig. 3, but now including the two LIRGs for which we have [CI](1–0) and CO(1–0) observations. The right part of the plot shows the resulting distribution of [CI]-based αCO values. The dotted line indicates the CO-to-H2 conversion factor for the Milky Way galaxy, and the dashed line corresponds to the value commonly used in the literature for (U)LIRGs (Downes & Solomon 1998). The mean and median values for the sample of (U)LIRGs remain the same as the ones obtained solely for the ULIRGs: αCO = 1.9 ± 0.3 M⊙ (K km s−1 pc2)−1 and 1.7 ± 0.4 M⊙ (K km s−1 pc2)−1 for the mean and median, respectively.

Figure A.3 shows the relations between the line luminosities and the general galaxy properties for the extended sample of (U)LIRGs. By including sources with L_IR < 12.0 L_⊙, we now retrieve tighter relations between all CO and [CI] line luminosities and infrared luminosity (top-left and bottom-left panels). On the contrary, the relations with SFR (middle panels) become weaker for all CO lines, and the [CI] line shows now a stronger correlation to SFR. Lastly, the relations with L_AGN become stronger for all line luminosities, with statistically significant relations for CO(1–0) (${\rho }_{L_{\text{CO}(\text{1-0})}^{\prime }-L_{\text{AGN}}}=0.53$, and p-value = 0.01), CO(2–1) (${\rho }_{L_{\text{CO}(\text{2-1})}^{\prime }-L_{\text{AGN}}}=0.37$ and p-value = 0.02), and interestingly, the highest correlation is measured for [CI](1–0) (${\rho }_{L_{[\text{CI}](\text{1-0})}^{\prime }-L_{\text{AGN}}}=0.68$ and p-value = 3.7 × 10⁻³). If this result were to be confirmed with larger statistics, this could hint of different behaviors for the two gas tracers in the presence of AGNs.

Fig. A.3. Same as Fig. 5, but now including the four LIRGs in the extended sample. In each plot, the dashed lines are the best-fit relations obtained from a least squares regression analysis conducted for each transition separately (color-coded according to the transition; see the legend in the top-left panel). In the bottom panels, the shaded gray areas correspond to the 1σ confidence interval of the fit. The bottom panels report also the Pearson correlation coefficients (ρ) and their associated p-values. The top-left panel shows also the $L_{\text{CO}(\text{1-0})}^{\prime }$-LIR relations found by Sargent et al. (2014) for local main-sequence (dot-dashed black line) and SB galaxies (dotted black line), and the $L_{\text{CO}(\text{1-0})}^{\prime }$-LTIR relation obtained by Herrero-Illana et al. (2019) using IRAM 30m CO(1–0) observations of 55 local sources in the GOALS sample including mostly LIRGs with LIR < 1012 L⊙ and with a < 20% AGN contribution to LIR. In the top-middle panel, for better visualization, we over-plotted with darker colors the data in bins of SFRs. In this plot, the solid purple line represents the $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR best fit obtained by Cicone et al. (2017) for the ALLSMOG and COLD GASS samples of local star-forming galaxies, described by ${\text{log}}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(8.16\pm 0.04)+(1.34\pm 0.07){\text{log}}_{10}\text{SFR}$.

Appendix B: The CO and [CI] line data sets: Tables and final reduced spectra

This appendix collects plots and tables describing the CO(1–0), CO(2–1), CO(3–2) and [CI](1–0) spectral line data sets used in this work.

Figures B.1 to B.6 report the final reduced CO(1–0), CO(2–1), CO(3–2), and [CI](1–0) line spectra of our sample of (U)LIRGs. All data, including archival ones, were re-reduced and reanalyzed by us in a consistent and uniform way, following the description provided in Sect. 3. The spectra are presented together with their best-fit multicomponent spectral models, computed following Sect. 4.1.

Fig. B.1. Integrated, continuum-subtracted CO(1–0), CO(2–1), CO(3–2) and [CI](1–0) spectra, including only the data employed in the analysis. We show in red the CO(1–0) transition, in green the CO(2–1) transition, in blue the CO(3–2) transition, and in magenta the [CI](1–0) transition. The source name, spectral binning, and telescope are reported on the top-left corner of each plot. This figure shows sources: IRAS 00188-0856, IRAS 01003-2238, IRAS F01572+0009, IRAS 03521+0028, IRAS F05024-1941, IRAS F05189-2524, and IRAS 06035-7102.

Fig. B.2. Continued from Fig. B.1. This figure shows sources: IRAS 06206-6315, IRAS 07251-0248, IRAS 08311-2459, IRAS 09022-3615, IRAS 10378+1109, IRAS 11095-0238, and IRAS F12074-0444.

Fig. B.3. Continued from Fig. B.2. This figure shows sources: IRAS F12112+0305, IRAS 13120-5453, IRAS F13305-1739, IRAS F13451+1232, IRAS F14348-1447, IRAS F14378-3651, and IRAS F15462-0450.

Fig. B.4. Continued from Fig. B.3. This figure shows sources: IRAS 16090-0139, IRAS 17208-0014, IRAS 19254-7245, IRAS F19297-0406, IRAS 19542+1110, IRAS 20087-0308, and IRAS 20100-4156.

Fig. B.5. Continued from Fig. B.4. This figure shows sources: IRAS 20414-1651, IRAS F20551-4250, IRAS F22491-1808, IRAS F23060+0505, IRAS F23128-5919, IRAS 23230-6926, and IRAS 23253-5415.

Fig. B.6. Continued from Fig. B.5. This figure shows sources: IRAS F23389+0300, IRAS F00509+1225, PG1126-041, IRAS F12243-0036, and PG2130+099.

Table B.1 collects the telescope and project number information for all data used in this work, including the duplicated data not employed in the main analysis (see Appendix C).

Table B.2 lists the values of spectral extraction aperture, the observing time, and the RMS noise values for all the spectral line data sets employed in the analysis.

Appendix C: Duplicated data sets not employed in the main analysis

Here we report additional spectra reduced for the sources in our sample, which were not included in the main analysis because of the availability of better quality and/or higher priority data, as explained in Sect. 3.1. Figures C.1, C.2, and C.3 show the [CI](1–0), CO(2–1), and CO(3–2) duplicated spectra, respectively.

Fig. C.1. Duplicated [CI](1–0) spectra not used in the analysis.

Fig. C.2. Duplicated CO(2–1) spectra not used in the analysis.

Fig. C.3. Duplicated CO(3–2) spectra not used in the analysis.

Additional details concerning individual sources are described below.

IRAS F15462-0450 The APEX archival CO(2–1) data set for the source IRAS F15462-0450 was discarded from the analysis because an inspection of the header signals a suspected pointing issue during the observation. This may be the cause of significant flux loss that led to a non detection. The ALMA archival CO(1–0) data for this source show a total flux of 43 ± 16 [Jy km s⁻¹], which would lead to detectable CO(2–1) transition by APEX, which supports our hypothesis of pointing issues. We checked the literature for previous CO(2–1) observations of this source and found an IRAM 30m telescope data set from 2008 published by Xia et al. (2012), where they computed a total integrated CO(2–1) flux of 21.04 ± 1.58 [Jy km s⁻¹] for IRAS F15462-0450.

IRAS F12112+0305 The ALMA CO(1–0) line observations of IRAS F12112+0305 show a low flux value of 34.6 ± 1.1 [Jy km s⁻¹], hinting at a significant missing flux. We proceeded to check the literature for previous CO(1–0) observations of this particular source, and found the work by Chung et al. (2009) reporting a Five College Radio Astronomy Observatory (FCRAO) 14 m telescope observation carried out between 2007 and 2008, with a total CO(1–0) line flux of 64.3 ± 21.04 [Jy km s⁻¹].

Appendix D: Tables listing line fluxes, luminosities, and luminosity ratios

Table D.1 lists the integrated line fluxes and luminosities, obtained through the spectral line fitting, computed using the relations in Sect. 4.2.

Appendix E: Additional relations explored

In this appendix, we show the relations between the r_CICO and different galaxy properties: SFR, L_AGN, M_mol, τ_dep, and α_AGN (Fig. E.1). In the main body of the text we plot the relation with L_IR (see Fig. 9). We do not retrieve any significant relations for any of the quantities studied here, as can be seen from the Pearson correlation coefficients and p-values reported in each graph.

Fig. E.1. rCICO plotted as a function of different galaxy properties, namely: SFR, LAGN, Mmol, τdep, and αAGN. Symbols and notation are as in Fig. 9 for our sample of ULIRGs.

All Tables

All Figures

Fig. 1. CO(1–0) vs. [CI](1–0) luminosity (top) and CO(2–1) vs. [CI](1–0) luminosity (bottom) for the ULIRGs in our sample. The best-fit relations are shown as dashed orange lines. The best-fit parameters are reported at the bottom-right corner of the plots. We also display the Pearson correlation coefficients (ρ) and their associated p-values. The solid black lines in both panels represent the corresponding relations reported by Jiao et al. (2017) for a sample of 71 (U)LIRGs, and the dotted lines represent the relations of Jiao et al. (2019) for a sample of 15 nearby spiral galaxies, between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$.

In the text

Fig. 2. [CI](1–0) line width as a function of the CO(2–1) line width. Top panel: Velocity dispersion (σv) obtained via a single Gaussian spectral fit to the [CI](1–0) and CO(2–1) emission lines. Middle and bottom panels: 16–84 (v84 − v16) and 2.3–97.7 (v97.7 − v2.3) percentile velocity intervals for CO and [CI], respectively, derived from the best-fit function obtained from the multi-Gaussian fit. Orange data points correspond to the ULIRGs sample analyzed in this work. The purple pentagon and the pink dot represent the values obtained for NGC 6240, respectively for the total region and for the inner ∼2″ (see Cicone et al. 2018, and the text in Sect. 5.1.2 for a more detailed explanation of the apertures used for spectra extraction). The solid black line indicates the 1:1 relation. The dashed orange line shows the best-fit relation obtained using only our sample, while the dashed purple line is the best-fit relation obtained for our sample plus NGC6240 (total). The shaded gray areas corresponds to the 1σ confidence interval of the two fits. At the bottom-right corner we report the Pearson correlation coefficient (ρ) and the p-values, as well as the best-fit coefficients of the purple fit.

In the text

Fig. 3. Ratio between [CI]-based molecular gas mass estimates and CO(1–0) line luminosity, computed for the ULIRGs in our sample that have both lines available, and plotted as a function of the infrared luminosity. This ratio is interpreted as the αCO factor (see Eq. (10)). The right part of the plot shows the resulting distribution of [CI]-based αCO values. The dotted line indicates the CO-to-H2 conversion factor for the Milky Way, and the dashed line corresponds to the value commonly used in the literature for (U)LIRGs (Downes & Solomon 1998). The mean and median values obtained for our sample (reported at the top-right corner) are αCO = 1.9 ± 0.4 M⊙ and 1.7 ± 0.5 M⊙ (K km s−1 pc2)−1, respectively.

In the text

Fig. 4. [CI]-based Mmol vs. CO-based Mmol values obtained with different αCO assumptions (0.8, 4.3, and 1.7 M⊙ (K km s−1pc2)−1). The best-fit linear relation with a free-varying slope is shown as a dashed orange line, and the best-fit parameters are reported at the bottom-right corner of each plot. The solid orange line shows the best-fit relation obtained by fixing the slope to unity, and the 1:1 relation is reported as a dotted black line.

In the text

Fig. 5. Global CO(1–0), CO(2–1), CO(3–2) (top panels) and [CI](1–0) (bottom panels) line luminosity plotted as a function of LIR (left), SFR (middle), and LAGN(right) for our sample of ULIRGs. In each plot, the dashed lines are the best-fit relations obtained from a least squares regression analysis conducted for each transition separately (color-coded according to the transition; see the legend in the top-left panel). In the bottom panels, the shaded gray areas correspond to the 1σ confidence interval of the fit. The bottom panels also report the Pearson correlation coefficients (ρ) and their associated p-values. The top-left panel also shows the $L_{\text{CO}(\text{1-0})}^{\prime }$ − LIR relations found by Sargent et al. (2014) for local main-sequence (dot-dashed black line) and SB galaxies (dotted black line), and the $L_{\text{CO}(\text{1-0})}^{\prime }$ − LTIR relation obtained by Herrero-Illana et al. (2019) using IRAM 30m CO(1–0) observations of 55 local sources in the GOALS sample, which comprises mostly LIRGs with LIR < 1012 L⊙ and with a < 20% AGN contribution to LIR. In the top-middle panel, for better visualization, we over-plot with darker colors the data in bins of SFRs. In this plot, the solid purple line represents the $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR best fit obtained by Cicone et al. (2017) for the ALLSMOG and COLDGASS samples of local star-forming galaxies, described by ${\text{log}}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(8.16\pm 0.04)+(1.34\pm 0.07){\text{log}}_{10}\text{SFR}$.

In the text

Fig. 6. Distribution of galaxy-integrated CO line ratios obtained for our sample of ULIRGs, with mean and median values reported at the top right of each plot. Left: $r_{21}\equiv L_{\text{CO}(\text{2-1})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$. Middle: $r_{31}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$. Right: $r_{32}\equiv L_{\text{CO}(\text{3-2})}^{\prime }/L_{\text{CO}(\text{2-1})}^{\prime }$. The mean CO line ratios measured by Leroy et al. (2022) in local massive main-sequence disk galaxies are shown using vertical dotted lines, with the gray shaded regions representing their 16-84 velocity percentile ranges.

In the text

Fig. 7. CO line ratios as a function of galaxy properties: LIR (top left), SFR (top middle), molecular gas depletion timescale due to star formation τdep ≡ Mmol/SFR (top right), molecular gas mass Mmol (computed by using the average αCO derived in this work; bottom right), LAGN (bottom middle), and AGN fraction (LAGN/Lbol; bottom right). We show the binned values using darker colors, overplotted on the individual data points. Dashed lines indicate the best-fit relations obtained from a least squares regression analysis conducted for each transition line separately. The color coding refers to different line ratios, as indicated in the legend at the top-right corner of each plot.

In the text

	Fig. 8. Distribution of galaxy-integrated $r_{\text{CICO}}\equiv L_{[\text{CI}](\text{1-0})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }$ line ratios obtained for our sample of ULIRGs, with mean and median values reported on the plot.
In the text

Fig. 9. rCICO as a function of LIR for our sample of ULIRGs. Similar relations with other galaxy properties are reported in Fig. E.1. Sources with direct measurements of both the [CI](1–0) and CO(1–0) lines are plotted using purple crosses. The pink diamonds represent sources observed in [CI](1–0) but without CO(1–0) line data, for which we inferred $L_{\text{CO}(\text{1-0})}^{\prime }$ from the CO(2–1) line luminosity by assuming r21 = 1.1, which is the mean value computed for our sample (see Fig. 6). The dashed purple line represents the best fit using a least squares regression analysis, and the shaded gray area corresponds to the 1σ confidence interval of the fit. The Pearson correlation coefficients (ρ) are reported at the upper-left corner of the panel with their associated p-values.

In the text

Fig. 10. Individual spectral component analysis of the CO line ratios as a function of the gas kinematics for our sample of ULIRGs. We set a detection threshold of 2σ for the flux of each Gaussian component in order to compute the ratios. The top, central, and bottom rows show, respectively, the r21, r31, and r32 plots. Left column: Line ratios as a function of the velocity dispersion, σv (line width of the spectral Gaussian component). Middle column: Line ratios as a function of the module of the central velocity of the spectral Gaussian component (|vcen|), measured with respect to the galaxy’s redshift, reported in Table 1. Right column: Line ratios as a function of the σv + |vcen| parameter. The dashed lines correspond to the best-fit relations obtained from a least squares regression analysis, and the shaded gray areas correspond to the 1σ confidence interval of the fit. The Pearson correlation coefficients and associated p-values are reported at the top-right corner of each plot, although they are less informative than the best fit because they do not take the (large) error bars into account.

In the text

Fig. 11. Global CO line ratios as a function of OHmax (left) and OHEQW (right) for the sources with OH119 μm data as reported in the literature (see Table 1). The binned values (computed in the same way as in Sect. 5.4) are shown using darker colors and overplotted on the individual data points. Dashed lines indicate the best-fit relations obtained from a least squares regression analysis conducted for each transition line separately. The color coding refers to different line ratios, as indicated in the legend at the top-right corner of each plot.

In the text

	Fig. 12. Global r21 ratios plotted as a function of [CI]-derived αCO values. The purple dashed line shows the best-fit relation, with its corresponding 1σ confidence interval represented by the gray shaded area. Neither the Pearson correlation coefficient (ρ) nor the regression analysis suggest any correlation between these quantities.
In the text

Fig. 13. Global rCICO as a function of the velocity dispersion of the gas. Top row: [CI](1–0)/CO(1–0) luminosity ratio as a function of σv, CO from the single Gaussian fit (see Sect. 5.1.2, left), and as a function of v97.7 − v2.3 velocity interval for CO (right). Bottom row: rCICO line ratio as a function of OHmax (left) and OHEQW (right). Sources with direct measurements of both the [CI](1–0) and CO(1–0) lines are plotted using purple crosses. For the other sources (pink diamonds), the CO(1–0) line luminosity was estimated assuming the average CO(2–1)/CO(1–0) luminosity ratio of r21 = 1.1 computed for our sample of ULIRGs. The dashed purple line represents the best fit using a least squares regression analysis, and the shaded gray area corresponds to the 1σ confidence interval of the fit. The Pearson correlation coefficients (ρ) are reported at the upper left corner of the panel with their associated p-values.

In the text

	Fig. 14. Depletion timescale (τdep ≡ Mmol/SFR) plotted as a function of LAGN. The dashed line shows the best-fit relation and its 1σ confidence interval plotted as a shaded gray area. The best-fit parameters and the Pearson correlation coefficient are reported at the upper left corner of the panel.
In the text

Fig. A.1. Same as Fig. 1, now including the two LIRGs (in purple data points) for which we have [CI](1–0) observations. The best-fit relations (shown as dashed purple lines) now include the LIRGs data points. The best-fit parameters are reported at the bottom-right corner of the plots. We also display the Pearson correlation coefficients (ρ) and their associated p-values. The solid black lines in both panels represent the corresponding relations reported by Jiao et al. (2017) for a sample of 71 (U)LIRGs, and the dotted lines represent the relations of Jiao et al. (2019) for a sample of 15 nearby spiral galaxies, between $L_{\text{CO}(\text{1-0})}^{\prime }$ and $L_{[\text{CI}](\text{1-0})}^{\prime }$.

In the text

Fig. A.2. Same as Fig. 3, but now including the two LIRGs for which we have [CI](1–0) and CO(1–0) observations. The right part of the plot shows the resulting distribution of [CI]-based αCO values. The dotted line indicates the CO-to-H2 conversion factor for the Milky Way galaxy, and the dashed line corresponds to the value commonly used in the literature for (U)LIRGs (Downes & Solomon 1998). The mean and median values for the sample of (U)LIRGs remain the same as the ones obtained solely for the ULIRGs: αCO = 1.9 ± 0.3 M⊙ (K km s−1 pc2)−1 and 1.7 ± 0.4 M⊙ (K km s−1 pc2)−1 for the mean and median, respectively.

In the text

Fig. A.3. Same as Fig. 5, but now including the four LIRGs in the extended sample. In each plot, the dashed lines are the best-fit relations obtained from a least squares regression analysis conducted for each transition separately (color-coded according to the transition; see the legend in the top-left panel). In the bottom panels, the shaded gray areas correspond to the 1σ confidence interval of the fit. The bottom panels report also the Pearson correlation coefficients (ρ) and their associated p-values. The top-left panel shows also the $L_{\text{CO}(\text{1-0})}^{\prime }$-LIR relations found by Sargent et al. (2014) for local main-sequence (dot-dashed black line) and SB galaxies (dotted black line), and the $L_{\text{CO}(\text{1-0})}^{\prime }$-LTIR relation obtained by Herrero-Illana et al. (2019) using IRAM 30m CO(1–0) observations of 55 local sources in the GOALS sample including mostly LIRGs with LIR < 1012 L⊙ and with a < 20% AGN contribution to LIR. In the top-middle panel, for better visualization, we over-plotted with darker colors the data in bins of SFRs. In this plot, the solid purple line represents the $L_{\text{CO}(\text{1-0})}^{\prime }$-SFR best fit obtained by Cicone et al. (2017) for the ALLSMOG and COLD GASS samples of local star-forming galaxies, described by ${\text{log}}_{10}{L}_{\text{CO}(\text{1-0})}^{\prime }=(8.16\pm 0.04)+(1.34\pm 0.07){\text{log}}_{10}\text{SFR}$.

In the text

Fig. B.1. Integrated, continuum-subtracted CO(1–0), CO(2–1), CO(3–2) and [CI](1–0) spectra, including only the data employed in the analysis. We show in red the CO(1–0) transition, in green the CO(2–1) transition, in blue the CO(3–2) transition, and in magenta the [CI](1–0) transition. The source name, spectral binning, and telescope are reported on the top-left corner of each plot. This figure shows sources: IRAS 00188-0856, IRAS 01003-2238, IRAS F01572+0009, IRAS 03521+0028, IRAS F05024-1941, IRAS F05189-2524, and IRAS 06035-7102.

In the text

	Fig. B.2. Continued from Fig. B.1. This figure shows sources: IRAS 06206-6315, IRAS 07251-0248, IRAS 08311-2459, IRAS 09022-3615, IRAS 10378+1109, IRAS 11095-0238, and IRAS F12074-0444.
In the text

	Fig. B.3. Continued from Fig. B.2. This figure shows sources: IRAS F12112+0305, IRAS 13120-5453, IRAS F13305-1739, IRAS F13451+1232, IRAS F14348-1447, IRAS F14378-3651, and IRAS F15462-0450.
In the text

	Fig. B.4. Continued from Fig. B.3. This figure shows sources: IRAS 16090-0139, IRAS 17208-0014, IRAS 19254-7245, IRAS F19297-0406, IRAS 19542+1110, IRAS 20087-0308, and IRAS 20100-4156.
In the text

	Fig. B.5. Continued from Fig. B.4. This figure shows sources: IRAS 20414-1651, IRAS F20551-4250, IRAS F22491-1808, IRAS F23060+0505, IRAS F23128-5919, IRAS 23230-6926, and IRAS 23253-5415.
In the text

	Fig. B.6. Continued from Fig. B.5. This figure shows sources: IRAS F23389+0300, IRAS F00509+1225, PG1126-041, IRAS F12243-0036, and PG2130+099.
In the text

	Fig. C.1. Duplicated [CI](1–0) spectra not used in the analysis.
In the text

	Fig. C.2. Duplicated CO(2–1) spectra not used in the analysis.
In the text

	Fig. C.3. Duplicated CO(3–2) spectra not used in the analysis.
In the text

	Fig. E.1. rCICO plotted as a function of different galaxy properties, namely: SFR, LAGN, Mmol, τdep, and αAGN. Symbols and notation are as in Fig. 9 for our sample of ULIRGs.
In the text