© The Authors 2024

1. Introduction

In the local Universe (z ≤ 0.2), most (ultra-)luminous infrared galaxies ((U)LIRGs), defined by their high infrared luminosities (L_IR(8 − 1000 μm)≥10¹¹ L_⊙), correspond to galaxy mergers that can trigger accelerated star formation – starbursts (SBs) – or accretion of matter onto the supermassive black hole (SMBH) in the center of the galaxy – active galactic nuclei (AGNs; Sanders & Mirabel 1996; Genzel et al. 1998; Lonsdale et al. 2006; Pérez-Torres et al. 2021; U 2022). The stellar and AGN feedback in these sources generate the most powerful galactic outflows known, which can embed gas in different phases and affect the host galaxy from scales of a few parsecs in the interstellar medium (ISM) to several tens of kiloparsecs out to the circumgalactic medium (CGM; Cicone et al. 2015).

The molecular gas phase of the ISM is the raw material from which stars form, and so constraining its physical properties in different gas environments, and in particular diffuse and extended outflowing gas, is important to understand the evolution of galaxies. Being rich in dense gas and dust, (U)LIRGs are ideal targets for (sub)millimeter observations of the molecular gas phase. Indeed, in the past decade, thanks to Herschel and ground-based submillimeter observatories, massive molecular (H₂) outflows have been detected in the majority of local (U)LIRGs (see Veilleux et al. 2020 for a review). In recent years, observations of different gas tracers, such as CO, HCN, HCO+, among others, of these sources have also allowed to study for the implications of galactic outflows in the galaxy evolution scheme to be studied. Galactic outflows may be the culprits of either removing large fractions of gas from the host galaxy, and hence decreasing the star formation rate (SFR; negative feedback; e.g., Di Matteo et al. 2005; Tumlinson et al. 2017), or compressing the available gas, resulting in an increased SFR (positive feedback; e.g., Maiolino et al. 2017; Gallagher et al. 2019). Despite extensive theoretical and computational work aimed at understanding galaxy-scale outflows, the propagation of cold outflows to scales of ∼10 kpc or larger, and their interaction with the gas phases embedded in the circumgalactic medium (CGM), is particularly difficult to model because of high numerical resolution requirements in cosmological simulations (Hummels et al. 2019; Suresh et al. 2019; van de Voort et al. 2019; Schimek et al. 2024). Therefore, further observational constraints of cold gas outflows, including their most extended components that are often faint or missed in high-resolution interferometric observations, are strongly needed to inform theoretical models.

Given that near-IR and mid-IR H₂ lines are not sensitive to the cold (i.e., T < 100 K) molecular phase, we need proxy tracers. The most common molecular tracer is the CO molecule, as it is produced in similar conditions to molecular hydrogen and is collisionally thermalized by H₂ even at low densities (Bothwell et al. 2013). Its low-J rotational transitions are easily excited as they have low critical densities (∼300 cm⁻³ for CO(1−0)), and they are thus sensitive to the total gas reservoir. The higher-J transitions have larger critical densities (∼10⁴ cm⁻³ for CO(5−4)), and they arise from denser and/or warmer regions in the gas. As such, observations of the different rotational transitions of a molecule can be used as a probe of the physical conditions of the gas, as its excitation depends on a competition between radiative and collisional transitions, with the latter depending on the temperature and density of the gas. The shape of the rotational ladder of CO – also referred to as the spectral line energy distribution (SLED) –, which is a plot of the CO(J_up → J_up − 1) integrated line fluxes as a function of the upper rotational quantum number J_up, is therefore sensitive (albeit in a degenerate way) to the H₂ gas density and temperature.

It has been proposed that the ISM in ULIRGs (as well as other galaxy populations) is composed of a low-excitation cold component phase that traces most of the H₂ gas mass, embedded in a warm, turbulent, and diffuse component, traced by higher excitation CO emission lines (Imanishi et al. 2023; Aalto et al. 1995). In this sense, the CO SLEDs may be broken down into three different regions (e.g., Rosenberg et al. 2015; Decarli et al. 2020; Boogaard et al. 2020; Esposito et al. 2022): (i) the low-J transitions (J_up ≲ 3) that trace the bulk of the molecular gas mass with cold temperatures and low densities, (ii) the mid-J transitions (4 ≲ J_up ≲ 9) which typically trace warmer and denser gas related to star formation, the photon-dominated regions (PDRs), and (iii) the high-J transitions (J_up ≳ 10) tracing extreme conditions such as shocks, X-ray dominated regions (XDRs) or cosmic ray-dominated regions (CRDRs) often linked to AGN activity. Combining information from the different rotational transitions of the CO molecule leads to a better understanding of the heating mechanisms behind the observed excitation of the ISM in galaxies.

Because SBs and AGNs introduce intense mechanical and radiative feedback into the ISM, it is expected that these have strong effects on the gas excitation, and thereby, on CO SLEDs. The incidence of massive molecular outflows in ULIRGs is extremely high, > 70% (Veilleux et al. 2013; Lamperti et al. 2022), pointing to extremely powerful feedback mechanisms at work in these sources. As a matter of fact, studies performed on local (U)LIRGs over the last decade thanks to Herschel, have shown that the global CO SLEDs of these sources peak at J ≥ 6, and may even remain flat up to J ≃ 10 (Kamenetzky et al. 2014; Lu et al. 2014; Rosenberg et al. 2015; Greve et al. 2014; Liu et al. 2015; Mashian et al. 2015), while the CO ladder of normal disk galaxies has been observed to peak at lower J ≃ 3 (Dannerbauer et al. 2009; Crosthwaite & Turner 2007).

Given the high incidence (≳70%) of massive outflows in ULIRGs, their extreme CO SLEDs may be connected to the powerful feedback mechanisms at work in these galaxies. However, a clear connection between molecular outflows and CO gas excitation has never been established: on the one hand, CO excitation studies of outflows are sparse and only sample a few, low-J CO lines (see references and discussion in Montoya Arroyave et al. 2023); on the other hand, the Herschel spectroscopic observations of multiple CO transitions, which have been used to probe the CO SLEDs of ULIRGs up to very high J_up, are not adequate to study the connection between outflows and CO excitation due to their low spectral resolution, which prevent a reliable study of the CO line profile.

It is worth noting, however, the diversity in the observed CO SLEDs for a given galaxy population (Narayanan & Krumholz 2014). Indeed, while the CO fluxes measured in ULIRGs are characterized by higher excitation than those measured in normal SF disk galaxies, and trends have been found, for instance, between excitation and high IR luminosity and/or high dust temperature (traced by the 60 μm/100 μm color index) (see e.g., Lu et al. 2014; Rosenberg et al. 2015; Pearson et al. 2016), their SLEDs span a wide range of excitation for J_up ≳ 5 lines (see e.g., Papadopoulos et al. 2012; Mashian et al. 2015). Similarly, observed CO SLEDs of dusty star forming galaxies selected through their submillimeter fluxes (also known as submillimeter galaxies, SMGs), which are commonly considered ULIRG-type galaxies at high-z, also show a wide variety in their measured CO excitation (see e.g., Harrington et al. 2021; Yang et al. 2017; Spilker et al. 2014; Bothwell et al. 2013; Carilli & Walter 2013). Such galaxy selections based on arbitrary luminosity- or color-cuts can result in heterogeneous physical conditions. Consequently, no “template” SLED exists for a given population such as (U)LIRGs, and we must pay particular attention when deriving average gas conditions for galaxy samples.

While modeling individual CO SLEDs is necessary for understanding the different physical conditions that galaxies host, we can still, however, study statistical trends in galaxy populations. ULIRGs are known to host powerful galactic outflows, leading to large velocity gradients (LVG) which may be related to the high excitation CO SLEDs observed in these sources. Although spatially resolved observations are needed to determine and characterize outflowing and non-outflowing material, single-pointing single-dish observations that deliver the total integrated molecular line emission constitute an essential first step for recovering the total flux of the poorly explored high-J CO lines, whose luminosity is unknown especially for the broad line wings that can trace outflows. This was the aim of our APEX observing campaign targeting 40 local (U)LIRGs, complemented by archival APEX and ALMA data, whose first results based on the low-J (J ≤ 3) CO and [CI](1−0) lines have been reported in Montoya Arroyave et al. (2023).

One of the results presented in Montoya Arroyave et al. (2023) was the measurement of extremely high galaxy-integrated low-J CO ratios in local (U)LIRGs, significantly offset from those measured in local massive main-sequence (MS) galaxies. However, from a spectral decomposition of the CO lines into low-velocity/narrow and high-velocity/broad line components, which, in first approximation, are dominated respectively by disk and outflowing gas, we found no evidence of consistently higher CO excitation in the latter, based on CO transitions up to J_up = 3. These results highlighted the need for higher-J CO lines to better constrain the excitation in these sources, as CO(J_up < 3) lines have little diagnostic value concerning the molecular ISM excitation in local (U)LIRGs, contrary to what has been found for normal star-forming galaxies in the local Universe (e.g., Leroy et al. 2022; Lamperti et al. 2020).

In this follow-up paper, we extend the analysis performed in Montoya Arroyave et al. (2023) for a subsample of 17 local ULIRGs, with high-sensitivity single-pointing APEX observations of CO lines up to J_up = 7. This effort is the first of its kind as it can rely on sensitive, high frequency (up to ∼710 GHz), submillimetric ground-based single-dish data for a conspicuous sample of local ULIRGs, obtained with wide bandwidths and high spectral resolution, reduced and analyzed in a homogenous way, by paying particular attention to aperture effects and where the few scans with poor baselines have been removed from the spectra. The paper is organized as follows. In Sect. 2 we describe the sample, the observing strategy, and the data reduction and analysis. In Sect. 3 we construct the CO SLEDs and investigate the ISM excitation and physical properties, where they are also contextualized through a comparison with relevant results from the literature. In the same section we use radiative transfer models to extract the average density, temperature and column density of the global ISM observed in our sample, and comment on caveats of deriving global gas properties. In Sect. 4 we explore correlations between the higher CO excitation of ULIRGs and the appearance of broad and/or high-velocity spectral components, to investigate the links with massive molecular outflows. Finally, we summarize our results and conclusions in Sect. 5. 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. Observations and data

2.1. The sample

The parent sample of our study, that is, the 40 (U)LIRGs presented in Montoya Arroyave et al. (2023), is representative the local ULIRG population in terms of physical properties, such as star formation rate (SFR), AGN luminosity (L_AGN), and AGN fractions (α_AGN ≡ L_AGN/L_bol). We now extend the analysis to higher-J CO data (4 ≤ J_up ≤ 7) for 17 ULIRGs, which span the entire range of L_AGN and SFR of the parent sample, in the redshift range 0.04 < z < 0.18. Previous Herschel OH119 μm observations, which are discussed in more detail in Montoya Arroyave et al. (2023), have reliably detected molecular outflow signatures in 11 out of the 17 galaxies studied in this paper (65% detection rate, which is typical for ULIRGs). Table 1 lists the subsample analyzed throughout this paper.

2.2. Observing strategy

2.2.1. Low-J CO data

A full description of low-J (J_up ≤ 3) CO data was given in Montoya Arroyave et al. (2023), hence here we only summarize the main points of the data acquisition and processing. The low-J CO data combine proprietary and archival single-dish (APEX) and archival interferometric (ACA, ALMA and IRAM PdBI) observations. Since we aim to include in our analysis the most extended and diffuse molecular gas components of the ISM of these sources, we gave priority to observations with the highest sensitivity to large-scale structures. In this sense, when data from multiple instruments are available for the same source and transition, we set the following priorities: (1) our own APEX proprietary (PI) data, (2) archival APEX or IRAM PdBI data, (3) archival ACA data, and, lastly, (4) archival ALMA data. This means that, when available, we use our own APEX PI observations since they were designed expressly for our science goals, although this is not always possible (for example, APEX does not have a 100 GHz receiver providing CO(1−0) line coverage). For the subsample of 17 ULIRGs presented in this work, we have the following low-J CO data coverage: CO(1−0) spectral line data for 10 galaxies, CO(2−1) line observations for all 17 sources, and lastly, 16 galaxies with CO(3−2) line coverage.

The spectral setup and sensitivity goals of the APEX PI data are described in detail in Montoya Arroyave et al. (2023). The receivers used were SEPIA180, SEPIA3445 and nFLASH230 depending on the observed frequency of interest, and all data were reduced and analyzed in a consistent way. For the APEX archival data, a compendium of different project datasets were used, using mainly the SHeFI and nFLASH receivers. For the single-dish APEX data, the beam sizes are ∼27 arcsec for the CO(2−1) line, and ∼20 arcsec for the CO(3−2) line. The apertures used for extracting the spectra for the interferometric data are reported in Montoya Arroyave et al. (2023). For all single-dish data, PI and archival, we adopted the same reduction and analysis steps, similar to the one explained in Sect. 2.3 for the high-J CO data presented in this work.

The archival ALMA and ACA data were gathered from different projects with different array configurations, leading to different angular resolutions and maximum recoverable scales. We re-reduced all these data paying particular attention to recovering the extended and diffuse emission (which may be affected by the choice of the cleaning mask, for example), and extracted the final spectrum from an aperture that was size-matched to maximize the recovered flux (more details can be found in Montoya Arroyave et al. 2023). For the sources where only high-resolution ALMA observations were available, we applied uv-tapering to enhance the sensitivity to extended structures. In the subsample analyzed in this work, only the CO(1−0) line observation of the source IRAS F01572+0009 is a high resolution ALMA dataset that may be biased against the most extended and diffuse components. We are aware, however, that uv-tapering does not overcome the issue of missing flux from faint extended structures due to the poor sampling of short uv baselines inherent to interferometric observations.

2.2.2. Higher-J CO data

All mid-J CO rotational transitions (with J_up = 4,  5,  6,  7) were observed with the APEX telescope between October 2020 and June 2021 (project ID E-0106.B-0674, PI: I. Montoya Arroyave) for the sources presented in Table 1. Our observing strategy was to reach a line peak-to-rms ratio of S/N ∼ 5 on the expected peak flux density in velocity channels δv ∼ 50 − 100 km s⁻¹, assuming a line width equal to that measured in the CO(2−1) line observation for each source.

The CO(4−3) (${\nu }_{\mathrm{CO}(4-3)}^{\mathrm{rest}}$ = 461.041 GHz) and CO(5−4) (${\nu }_{\mathrm{CO}(5-4)}^{\mathrm{rest}}$ = 576.268 GHz) observations were carried out with the nFLASH460 receiver, and the CO(6−5) (${\nu }_{\mathrm{CO}(6-5)}^{\mathrm{rest}}$ = 691.473 GHz) and CO(7−6) (${\nu }_{\mathrm{CO}(7-6)}^{\mathrm{rest}}$ = 806.652 GHz) observations were performed with the SEPIA660 receiver. Both instruments are frontend heterodyne with dual-polarization sideband-separating (2SB) receivers, and each sideband (upper -USB-, and lower -LSB-) is covered by a fast Fourier transform spectrometer (FFTS) unit. The nFLASH460 receiver covers a frequency range between 378 and 508 GHz, and has a IF coverage between 4 and 8 GHz, with an instantaneous coverage of 4 GHz per sideband and a 12 GHz gap between the sidebands. The SEPIA660 receiver can be tuned in the frequency range 597 − 725 GHz, and it has two IF outputs per polarization, each covering 4 − 12 GHz separated by a 16 GHz gap, leading to a total of up to Δν = 32 GHz IF bandwidth (Meledin et al. 2022). The average noise temperatures of the receivers are (T_rx) of ∼150 K for nFLASH460, and between ∼75 − 125 K for SEPIA660 (Belitsky et al. 2018). The beam sizes for the mid-J CO observations are ∼14,  11,  9 and 8 arcsec for the CO(4−3), CO(5−4), CO(6−5) and CO(7−6) lines, respectively.

For our observations, each sideband spectral window was divided into 4096 channels (∼977 kHz per channel), corresponding to a resolution in velocity units at the range of redshifts covered by our sample of ∼0.6 km s⁻¹ for the nFLASH460 observations, and ∼0.4 km s⁻¹ for the SEPIA660 observations. The telescope was tuned to the expected observed frequency¹ for each transition, with wobbler-switching symmetric mode with 100″ chopping amplitude and a chopping rate of R = 0.5 Hz. The data were calibrated using standard procedures (see Sect. 2.3). The on-source integration times (without overheads) varied from source to source, from ∼5 min up to ∼2 h. During the observing runs, the PWV varied from 0.2 < PWV [mm] < 1.0.

Our final data coverage for the mid-J CO lines is as follows: CO(4−3) line observations for 16 sources, among which there is a non-detection due to a nearby atmospheric line (IRAS 06035−7102); CO(5−4) line detections in two sources, CO(6−5) line detections in 11 sources, and lastly, three sources with CO(7−6) line measurements. Therefore, in total, we detected 74 CO lines with J_up ≤ 7 out of the 75 observed lines. Our data set includes between three and six CO lines per source. We adopted the same reduction and analysis steps for all data, which are summarized below in Sect. 2.3.

2.3. Data reduction

The low-J CO lines data reduction is explained in detail in Montoya Arroyave et al. (2023), both for the proprietary APEX data, and the archival ALMA/ACA data. For the high-J CO APEX data, we follow a similar procedure as that of the low-J CO APEX data. Here, we briefly go through the steps of such data reduction process.

We reduced the data using the GILDAS/CLASS software package². We began by isolating all the spectral scans for each source and transition and inspected them individually in each polarization, we then discarded those scans that were affected by baseline instabilities and other instrumental issues. We subtracted the baselines in individual scans after masking the spectral range corresponding to v ∈ ( − 500, 500) km s⁻¹, in order to avoid the expected line emission, and followed by co-adding all scans and then subtracting a linear baseline to the averaged spectrum. We fit a single Gaussian function to the line profile, and based on the resulting fit parameters (line width and central velocity), we adjusted the central mask to the velocity channels of the detected line emission. Lastly, we produced the final spectra and exported them with a common spectral resolution of δv ≈ 5 km s⁻¹, in order to allow for further rebinning in the spectral fitting phase. The spectra were imported into Python³, where the data analysis was performed.

Since the output signal from the APEX telescope corresponds to the antenna temperature corrected for atmospheric losses ($T_{\mathrm{A}}^{\ast }$) in units of Kelvin [K], it must be multiplied by a calibration factor (or telescope efficiency) in order to obtain the flux density in units of Jansky [Jy]⁴. For the nFLASH460 receiver, the average Kelvin to Jansky conversion factor measured during our observation period is 63 ± 11 Jy K⁻¹, and for SEPIA660 it is 57 ± 5 Jy K⁻¹. The APEX observations that are presented for the first time in this work, covering the mid-J CO lines (4 ≤ J_up ≤ 7) in 17 sources, are shown in Figs. B.1–B.17.

2.4. Data analysis

The analysis performed in this work is based on galaxy-integrated molecular line spectra, and as such, it is aimed at estimating the average physical conditions of the molecular ISM in our sample of ULIRGs using source-averaged CO SLEDs.

Once exported in Python, all CO line spectra were fit using mpfit⁵. This code uses the Levenberg–Marquardt technique to solve the least-squares problem in order to fit data (the observed spectrum) with a model (a multi-Gaussian function). Each Gaussian component is characterized by an amplitude, a peak position (v_cen), and a velocity dispersion (σ_v, or, equivalently, the full width at half maximum, FWHM). 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 e.g., source IRAS 09022−3615, Fig. B.7). At this stage, we do not assign any physical meaning to the different Gaussian components fitted to each transition, nor force the fit to use the same Gaussians for all CO transitions of the same source, as our goal is to simply reproduce the total spectral line profiles.

We verified that the velocity-integrated line fluxes measured through the Gaussian fits are consistent with those calculated by directly integrating the spectra within v ∈ ( − 1000, 1000) km s⁻¹, after setting a signal-to-noise threshold of S/N > 1.5 per channel. We find that the two methods give consistent results within the statistical errors. The total line fluxes are reported in Table A.1.

We then converted the velocity-integrated fluxes to line luminosities following Solomon et al. (1997):

$$\begin{array}{c}\hfill {L^{\prime }}_{\mathrm{line}}[\mathrm{K}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}\mathrm{p}{\mathrm{c}}^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 A.2 we report the luminosities calculated for the different transitions.

From the line luminosities we can compute the CO line ratios with respect to the CO(1−0) line luminosity, where we use the following notation:

$$\begin{array}{c}\hfill r_{J1}\equiv {L^{\prime }}_{\mathrm{CO}(J\to J-1)}/{L^{\prime }}_{\mathrm{CO}(1-0)}\end{array}$$(2)

where J corresponds to the upper rotational quantum number, J_up. The average global ratios calculated for our sample are reported in Table 2.

3. CO excitation analysis

3.1. CO SLEDs

Our results are displayed in Fig. 1 in the form of observed CO SLEDs (using integrated flux density in Jy km s⁻¹ vs. J_up), for all 17 ULIRGs. This plot demonstrates the wide range of gas excitations in the ULIRG population, visible from the different shapes of the observed CO ladders. Already from a visual inspection, we can observe sources with SLEDs peaking at J_up ∼ 3, 4 (e.g., IRAS F05189−2524, IRAS 09022−3615, IRAS F12112+0305, IRAS F13451+1232, IRAS 17208−0014, IRAS 19254−7245, IRAS F20551−4250 and IRAS F23128−5919), others showing a plateau-like feature from J_up ≳ 3 (IRAS F01572+0009 and IRAS F23060+0505), and lastly, sources for which no plateau nor turn-over point is observed up to the highest detected CO transition (IRAS 03521+0028, IRAS 06035−7102, IRAS 06206−6315, IRAS 08311−2459, IRAS 10378+1109, IRAS 16090−0139 and IRAS F22491−1808). The peak of the CO SLED sheds light on the gas excitation conditions, since it directly depends on the gas kinetic temperature and density. For source IRAS 03521+0028, with only three detected CO (J_up ≤ 4) lines and no visible turn-over, we cannot constrain the H₂ gas properties.

Fig. 1. CO spectral line energy distributions of the ULIRGs in our sample. The blue circles with the error bars are the velocity-integrated line fluxes in units of Jy km s−1. The dotted gray lines show the thermalized, optically thick limit (in all lines), SLED profiles, normalized to the CO(2−1) transition (which is the lowest-J CO transition that is available for all sources). The source IDs are reported on the top-left corner of each panel.

We also show, indicated by a gray dotted curve in each panel of Fig. 1, the trend expected for thermalized CO emission with fluxes rising as the square of J_up (assuming all lines are in the Rayleigh-Jeans regime). Star forming ULIRGs are expected to have a larger fraction of dense, excited CO gas that may be thermalized up to J_up = 3 or beyond (e.g., Weiß et al. 2007; Papadopoulos et al. 2012), in contrast to the low-density, low-excitation gas dominating the CO emission from more quiescent galaxies (Braine et al. 1993; Crosthwaite & Turner 2007; Boogaard et al. 2020). Indeed, Fig. 1 shows that CO(3−2) is only weakly subthermally excited for most sources in our sample. The average r₃₁ line ratio is 0.77 ± 0.02 (see Table 2), consistent with moderately subthermally excited lines. Conversely, CO(J_up ≥ 4) lines are significantly subthermally excited in most sources. For the higher-J transitions, we derived average galaxy-integrated line ratios of r₄₁ = 0.54 ± 0.03, r₅₁ = 0.42 ± 0.07, r₆₁ = 0.17 ± 0.007 and r₇₁ = 0.58 ± 0.04. It is worth noting, however, the lower statistics for the J_up ≥ 5 ratios, and in particular J_up = 5, 7, for which we only have two measurements each; as such, the mean r₅₁ and r₇₁ values should not be considered representative of the ULIRG population.

The estimated line ratios are, overall, in agreement with those reported by Papadopoulos et al. (2012) for a sample of 70 local (U)LIRGs. Compared to other galaxy populations, such as high-z SMGs (see e.g., Bothwell et al. 2013; Spilker et al. 2014) and BzK galaxies (Daddi et al. 2015), our ULIRGs show higher CO ratios. For a sample of main sequence galaxies at z ∼ 1.25, Valentino et al. (2020) measured⁶ ⟨r₄₂⟩ = 0.36, ⟨r₅₂⟩ = 0.28 and ⟨r₇₂⟩ = 0.09, which are lower than the ones we measure for our sample of ULIRGs (r₄₂ = 0.50, r₅₂ = 0.38, and r₇₂ = 0.53). However, for galaxies above the MS (at z ∼ 1.25), classified as SBs and extreme SBs, they obtain higher values, closer to those measured in our work. We warn the reader not to compare these ratios directly, as the work by Valentino et al. (2020) has low statistics for CO(4−3) detections (four measurements), and our sample has low statistics for CO(5−4) and CO(7−6), and so, the values may be subject to biases. Even so, such comparison gives, at the very least, an idea of the extreme physical conditions of the molecular gas in local ULIRGs, with some presenting CO line ratios that resemble galaxies classified as SBs at high-z.

Sources IRAS 10378+1109 and IRAS 16090−0139 present some of the highest excitation SLEDs. Indeed, these two sources (out of three in our sample) with CO(7−6) detections, show a remarkable high flux in this transition and are the ones responsible for driving up the average r₇₁ ratio value. Both galaxies have interesting features in their observed spectra (Figs. B.8 and B.11), for which the CO(7−6) emission seems to be broader than their corresponding lower-J CO lines. In particular, the CO(7−6) emission from IRAS 10378+1109 appears to be red-shifted from the systemic velocity inferred from the low-J CO lines by δv ∼ 300 km s⁻¹, whereas no other source in the sample exhibits a similar behavior⁷. While there are galaxies for which we observe line wings becoming more prominent with increasing J_up, as can be seen, for instance, in IRAS 06035−7102 in Fig. B.4, in those cases the appearance of more prominent wings does not significantly affect the centroid of the line. The apparent decrease of CO(7−6) flux in the systemic velocity component of IRAS 10378+1109 is indeed puzzling and deserves further inspection, ideally with spatially resolved data that allow for a kinematical analysis of the gas distribution in this galaxy. Interestingly, the same source also shows a remarkable trend of nearly thermalized CO emission up to J_up = 7, indicating there is a large fraction of gas in this galaxy under extreme physical conditions, and a likely continuous injection of energy into the medium in order to significantly populate the high-J levels of the CO molecule. In Sect. 4 we discuss possible explanations for the observed high excitation in some of the ULIRGs in our sample.

The relative variations between the SLEDs of our ULIRGs sample are shown in the upper panel of Fig. 2, where we show all CO ladders in the same plot, normalized by the CO(2−1) velocity-integrated fluxes. By plotting all the sources together, the large variation in gas excitation within our sample, as traced by the CO(J > 3) lines, becomes more evident. In the bottom panel of Fig. 2 we show the average CO SLED of our sample (green circles), and compare it to other well-studied systems and galaxy populations. In particular, we show the average CO SLEDs measured by Papadopoulos et al. (2012) in local (U)LIRGs (blue stars in Fig. 2), and the one derived by Bothwell et al. (2013) for SMGs at z = 2.2 (blue crosses in Fig. 2). Our sample is about four times smaller than that studied by Papadopoulos et al. (2012), and hence our measured error bars are larger than those derived in their work. On the other hand, our data are more sensitive, can rely on broader bandwidths thanks to upgraded receivers, and are more uniform in terms of instrumental setups and analysis. Both Bothwell et al. (2013) and our survey can rely on only a few measurements of the CO(7−6) transition, with our ratio being dominated by two high excitation sources, and as such, affected by a large difference and uncertainty in the estimated mean value.

Fig. 2. Top panel: observed CO SLEDs for the ULIRGs in our sample. The color coding is shown in the legend on top of the figure. Bottom panel: we compute the average CO SLED of our sample, and compare it to other well-studied galaxies, the starburst galaxy M 82 (Mashian et al. 2015), the nearby ULIRG NGC 6240 (hosting a dual AGN and a SB; Papadopoulos et al. 2014), and the Milky Way inner disk (Fixsen et al. 1999). We also show the average CO SLEDs for a sample of local (U)LIRGs studied by Papadopoulos et al. (2012), and a sample of high-z SMGs (Bothwell et al. 2013). In both panels, the dotted gray line corresponds to a constant brightness temperature in the Rayleigh-Jeans regime (i.e., S​ν ∼ ν2), which is the trend expected for thermally excited CO emission, in the optically thick limit for all lines. All line fluxes are in units of Jy km s−1, normalized by the CO(2−1) flux measured by APEX.

While we are aware of the extreme nature of the ISM in ULIRGs, the difference between the CO SLEDs of our sample and the Milky Way inner disk (Fixsen et al. 1999) is still remarkable, as shown by the bottom plot of Fig. 2. In the Milky Way and in other disk galaxies, the H₂ gas excitation is dominated by much more quiescent gas, and so the CO SLED peaks at very low-J CO rotational transitions. The CO ladders of our sources are more in agreement with those measured in more extreme SBs, such as M 82 (the archetypal SB galaxy, see Mashian et al. 2015), and NGC 6240 (a known ULIRG to host both SB and AGN, see Papadopoulos et al. 2014). Nevertheless, while these two galaxies may show similar gas excitation up to J_up = 7, they differ strongly for higher rotational lines, as it was evidenced in the study performed by Mashian et al. (2015), analyzing Herschel CO data in the range J_up = 14 − 30. These authors found that the CO ladders of M 82 and NGC 6240, albeit similar up to J_up ≈ 8, differ strongly at higher-J, with the former peaking around J_up ∼ 8 and then quickly declining with increasing J_up, as opposed to NGC 6240 that shows strong emission lines even up to J_up = 28. This is a strong indicator of the need for CO(J_up > 8) lines to properly constrain the physical properties of the gas in SBs and (U)LIRGs, which will be evidenced also by our analysis.

3.2. Lower and higher excitation sources

To better visualize the different CO SLED shapes in our sample, we classify them into lower excitation and higher excitation⁸. We define lower-excitation sources those for which we observe a turn-over point in the CO SLED up to the observed rotational transition, and higher-excitation those for which we observe a plateau or increasing flux up to the observed CO transition. There is only one source for which such classification is not possible, i.e., IRAS 03521+0028, for which we only have data up to CO(4−3) and no detection of a clear peak, leaving the shape of the CO ladder highly uncertain. Figure 3 shows the CO SLEDs of lower-excitation (left panel) and higher-excitation (right panel) sources. The CO(6−5) transition appears to be the key emission line to identify the different trends of CO ladders. Indeed, the lower-excitation sources show an average CO SLED that peaks around J_up ∼ 3, 4, with flux ratios decreasing starting from J_up = 6. We note that these sources, albeit classified as “lower-excitation” in this work, have higher excitation CO SLEDs than those measured in Milky Way-like galaxies and typical star forming discs (Fixsen et al. 1999; Dannerbauer et al. 2009). The average ratios of the “lower-excitation” ULIRGs are ⟨r₃₁⟩_low ex = 1.01, ⟨r₄₁⟩_low ex = 0.45 and ⟨r₆₁⟩_low ex = 0.15.

Fig. 3. CO SLEDs of the ULIRGs classified as having lower excitation (left panel) and higher excitation (right panel) gas, following the classification criteria described in Sect. 3.2. Each source is shown in a different color, according to the legend reported on the top of each panel. The average CO SLED of each group of galaxies is also computed and plotted using black square symbols connected by a black solid line.

Conversely, the higher-excitation ULIRGs show peaks beyond the CO(4−3) transition. The average CO SLED derived for these sources shows a first plateau at J_up = 4, similar to the peak of the lower-excitation ones, and an unconstrained peak at higher J_up, which suggests the need of at least two components to reproduce the observed SLEDs. This apparent double-peaked CO ladder shape is driven, however, by the CO lines with lowest statistics in our sample, CO(5−4) and CO(7−6), for which we only have two and three measurements each, respectively, and so larger error bars; as such, the derived average CO SLED is subject to bias. It is, however, constructive to discuss the CO ladder trends observed in these sources, as it is not the first time that they have been detected in the ULIRG population. Indeed, it has been pointed out in previous studies that active SF galaxies likely host two (or even three) different molecular gas excitation components, one dominating the low-J CO line emission, tracing the quiescent gas in a more diffuse and cold phase, and a second, more excited component that traces more extreme environments than are associated to AGNs (XDRs) and/or shocks (Papadopoulos et al. 2014; Pereira-Santaella et al. 2014; Kirkpatrick et al. 2019). Such energetic processes are not uncommon in local ULIRGs, most of which are major galaxy mergers, where accelerated star formation and possibly feeding onto one or multiple SMBHs can be triggered by the gravitational interactions. In turn, the powerful radiative and mechanical feedback mechanisms resulting from the enhanced SMBH and SB activity can launch fast outflows embedding large amounts of molecular gas (up to ∼10¹⁰ M_⊙ in the most extreme cases such as NGC 6240, see Cicone et al. 2018), hence possibly affecting the global excitation conditions of the molecular ISM, which is the hypothesis that we aim to test in this work. The connection between CO excitation and massive molecular outflows will be further explored and discussed in Sect. 4.

3.3. LVG modeling

The focus of our work is on the relation between CO excitation and CO spectral broadening and asymmetries possibly tracing massive molecular outflows, which will be explored in Sect. 4. However, in the following we report for completeness an analysis of the global CO ladders performed with RADEX (see van der Tak et al. 2007), a non-LTE (local thermodynamic equilibrium) radiative transfer code, that uses an escape probability of β = (1−e^−τ)/τ (where τ is the optical depth of the line), derived from an expanding sphere geometry, to compute the molecular line intensities. This geometry, similar to LVG codes, works under the assumption that turbulent motions within star-forming clouds lead to large velocity gradients, and thus the photon escape probabilities are computed accordingly (Sobolev 1957; Ossenkopf 1997). These are commonly used to model molecular line intensities, hence gas excitation and SLEDs, in galaxies as they can predict the shape of the rotational ladder by computing the collisional excitation of molecules (CO in our case) for a range of physical conditions. In particular, we use the code developed by Yang et al. (2017), RADEX_EMCEE⁹, which combines pyradex¹⁰ and emcee¹¹ to model the CO fluxes, and find the best set of physical parameters that reproduce the SLEDs of our sample. The CO collisional data are taken from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005). Following Yang et al. (2010), we assume that H₂ is the only collisional partner of CO with an ortho-to-para ratio of 3. The input parameters of the code are: the H₂ volume density (n_H2), the kinetic temperature of the molecular gas (T_kin), the CO column density (N_CO), and the size of the emitting region (or solid angle of the source). The results for the relative populations in the different rotational levels (i.e., the excitation) will depend on N_CO/dv, while line fluxes will scale as ∼dv. The fluxes of all CO lines are equally scaled by the solid angle of the source, and as such the resulting shape of the CO SLED does not depend on this parameter and only acts as a normalization factor, hence only the first three parameters (n_H2, T_kin, N_CO/dv) are relevant.

Instead of using RADEX to generate a grid of line fluxes for a range of input parameters to then fit the observed CO SLEDs, the code developed by Yang et al. (2017) adopts a Bayesian approach for fitting the observed CO fluxes. A Markov chain Monte Carlo (MCMC) calculation is performed to explore the continuous parameter space and it calls pyradex in each iteration for computing the CO fluxes for a given combination of {n_H2, T_kin, N_CO/dv}, it then samples the posterior probability distribution function to obtain the next most probable set of parameters to use as inputs. This approach allows for a better sampling in the parameter space while avoiding unnecessary calculations, thus resulting in a faster convergence. We used a flat log-prior for all the parameters to set boundaries for the parameter space. The ranges are as follows: volume density n_H2 from 10^1.5 to 10^7.5 cm⁻³, kinetic temperature T_kin from the CMB temperature at the redshift of the source (T_CMB = 2.7315(1 + z)) to 10³ K, and CO column density per unit velocity gradient N_CO/dv from 10^14.5 to 10^19.5 cm⁻² km s⁻¹. Additionally, we also use a range for the velocity gradient dv/dr = 0.1 − 1000 km s⁻¹ pc⁻¹ (see Tunnard & Greve 2016), which sets a limit for the ratio between N_CO/dv and n_H2¹². The prior probabilities that fall outside of the boundaries are set to 0.

As discussed in Sect. 3.1, there is a wide range of CO excitation conditions within the ULIRG population, and while previous studies show that the CO emission is likely dominated by two excitation components, the lack of CO(J_up > 7) lines in our study prevents us from fitting two different components without being subject to significant degeneracy between the derived gas properties. For these reasons, in this work, we are fitting the observed CO SLEDs with a single-component LVG model, approximating all CO emission to be characterized by a single kinetic temperature, H₂ volume density, and column density per velocity gradient. The modeled CO SLEDs from the best-fitting estimates of (n_H2, T_kin, N_CO/dv) and their ±1σ uncertainties, are shown in Figs. B.1–B.17, together with our CO flux measurements.

The LVG modeling results are summarized in Table 3. We obtain best-fit values in the following ranges: n_H2 ≈ 10^2.0 − 10^5.8 cm⁻³, T_kin ≈ 10^1.2 − 10^2.9 K, and N_CO/dv ≈ 10^16.2 − 10^18.8 cm⁻² (km s⁻¹)⁻¹. These are overall consistent with the results of single-component LVG modeling of other samples of local ULIRGs (see e.g., Ao et al. 2008; Weiß et al. 2005; Rosenberg et al. 2015) and high-z SMGs (e.g., Cañameras et al. 2018; Yang et al. 2017; Spilker et al. 2014). However, for the higher-excitation sources, driven by the high mid-J CO line fluxes, we derive extreme conditions that are physically improbable to be representative of the average state of the ISM in these galaxies. For instance, the high fluxes measured in the CO(J_up = 4, 6) lines for IRAS 06206−6315 boost the average H₂ volume density of the galaxy to $log\left(n_{{\mathrm{H}}_2}[{\mathrm{cm}}^{-3}]\right)\approx 5.6_{-1.4}^{+1.1}$, and a similar effect is observed in IRAS F20551−4250, though in this latter case, the parameter boosted to higher value is the gas kinetic temperature to $log\left(T_{\mathrm{kin}}[\mathrm{K}]\right)\approx 2.8_{-0.2}^{+0.2}$. These values are indeed highly uncertain, as can be evidenced by the ±1σ uncertainties, or by the posterior probability distributions showed in Appendix B, where it becomes clear that the parameter space is highly biased to higher values, and so, overall poorly sampled. This effect is seen in IRAS F05189−2524, IRAS F12112+0305, IRAS 17208−0014, IRAS F20551−4250, and IRAS F23128−5919. The reason for this roots from the Ansatz of using a single LVG model. For instance, the CO ladder of IRAS 17208−0014 does not fall down steeply after its peak, but remains significantly excited. If we were to overplot a single LVG model with reasonable temperature (e.g., 50 K) that reproduces the data up to the peak at J_up = 4, it would drop much faster for transitions above peak compared to the measurement. This can easily be captured by adding a second component with a higher density that broadens the CO SLED. However, for our single component approach, the temperature has to reach very high values in order to fit the J_up = 6 line. A similar effect will be at play for those solutions with very high density (e.g., IRAS 06206−6315). In such cases, the CO SLED is not well sampled beyond the peak, thus the shape of the drop-off at higher-J CO lines remains unconstrained. The CO fluxes can then be reproduced with a low-temperature model, however, a correspondingly high density is needed to move the peak of the CO SLED to the highest J_up measurement. As a result, changing the flat log-priors in the MCMC to cover a wider range in the sampling of the parameter space, does not lead to better results, as the outcome still leads to unphysical high values of volume density and/or kinetic temperature in order to reproduce the high fluxes measured for higher-J CO lines. These results exemplify the degeneracy between the n_H2 and T_kin on the shape of the CO SLED, as an increase in either will lead to more excited gas. In order to break such degeneracy, we need complementary high resolution CO(J_up > 7) line data.

A single-component LVG fit is certainly an over simplification, and provides only a crude approximation to the physical properties that characterize the entirety of the ISM within these sources. The CO SLED plots reported in Appendix B confirm that the LVG single component modeling is not adequate for some of our sources. In particular, the sources mentioned above with poor sampling of the posterior probability distribution, as well as IRAS 06206−6315, IRAS 09022−3615, IRAS F22491−1808 and IRAS F23060+0505, are some of the galaxies that likely require a second excitation component. Many of these sources have already been studied with Herschel across the much broader range of CO transitions. Pearson et al. (2016) grouped the SLEDs of their ULIRG sample into three categories: flat, increasing fluxes (low- to high-J), and decreasing fluxes. These authors found that the SLED’s slopes correlate to the FIR colors as traced by the 60 μm/100 μm color index (lower values corresponding to colder sources), with warmer ULIRGs showing higher excitation CO SLEDs (a trend already seen in studies by Rosenberg et al. 2015; Lu et al. 2014). We also explored a possible correlation between the lower- and higher-excitation ULIRGs in our sample and their FIR color (using IRAS measured 60 μm and 100 μm fluxes), and we did not find a clear correlation of higher excitation for warmer ULIRGs. However, this may be due to the degeneracy caused from classifying sources into lower- and higher-excitation based on CO SLEDs that are sampled only up to CO(7−6) in our case. And indeed, some of the sources belonging to the lower-excitation group under our criteria, have been observed to have high-excitation CO SLEDs when probing CO(J > 7) lines (e.g., IRAS 09022−3615, IRAS F05189−2524, IRAS 17208−0014).

Pearson et al. (2016) showed that, while there are sources for which a single-component CO SLED fit is sufficient, the CO ladder of IRAS 09022−3615 (classified as a warm ULIRG in their study, and included also in our work) requires three excitation components in order to reproduce the observed CO fluxes up to J_up = 13: (i) a cold (47 K) component peaking at J_up ∼ 3, (ii) a warm (205 K) component peaking at J_up ∼ 6, and (iii) a third hot (415 K) component that dominates at high-J transitions, peaking around J_up ∼ 11. They propose that the hottest component traces the inner PDRs, shocks and/or AGN-heated gas (XDRs). In our analysis, the best-fit parameters derived for IRAS 09022−3615 from the single LVG component model are n_H2 ≈ 10^2.6 cm⁻³ and T_kin ≈ 250 K. This result is in agreement with the warm component (ii) of Pearson et al. (2016), which dominates the mid-J CO emission. Nevertheless, this source is classified as a lower-excitation ULIRG under the criteria used in our work. This is due to the lack of CO(J ≥ 7) lines that would evidence the high excitation nature of this source, demonstrating once more the need to observe higher-J CO lines to properly determine the physical conditions of galaxies.

Our sample overlaps also by two sources with that studied by Rosenberg et al. (2015). These authors classify the ULIRGs into three groups (similar to Pearson et al. 2016) depending on the drop-off slope of the CO ladder from the J = 5 − 4 transition and onwards. Their “Class I” objects have the steepest drop-offs, the “Class II” sources peak at around J = 6 − 5, and “Class III” objects are characterized by flat CO SLEDs. For Class II and III (the highly excited classes) sources, heating mechanisms besides UV heating are required to explain the high-J CO emission. Sources IRAS 17208−0014 and IRAS F05189−2524 are classified as Class III objects in their study, while in this work they both lie in our lower-excitation group. According to Pereira-Santaella et al. (2014), shocks play an important role in setting the higher-J CO excitation conditions of IRAS F05189−2524. Interestingly, in our analysis, the single-component LVG model does not provide a good fit to either of these two sources (see Figs. B.3 and B.12), with poorly sampled posterior probability distributions due to the bias caused by the high CO(J_up = 4, 6) fluxes, hence resulting in a derived kinetic temperature that is close to the upper boundary with misleading uncertainties.

Despite many of our results and several literature studies suggesting the need of at least two excitation components to properly reproduce the observed CO SLEDs of ULIRGs, and so derive accurate physical conditions of their ISM, we restrain ourselves from adding another component to our fitting procedure, as our ground-based CO data, lacking high-J CO transitions, do not allow us to resolve the ambiguity between the kinetic temperature and the H₂ volume density in the LVG models. Therefore, we caution against an over-interpretation of the specific best-fit parameters derived for individual sources reported in Table 3, specially those classified as higher-excitation sources, for which the single-component model proves to be inadequate. On the other hand, thanks to the high S/N and spectral resolution, our data are most adequate to perform an analysis of the links between spectral line profiles and global CO excitation conditions, which is not possible from the low spectral resolution Herschel spectra. Section 4 will focus on exploring such a link.

4. A link between CO excitation and molecular outflows

4.1. Analysis of line widths as a function of J_up of the CO transition

In the following, we explain how we tested the hypothesis of a relation between CO excitation and the appearance of broad and/or high-velocity CO spectral components, which may indicate that massive molecular outflows are responsible for setting the global CO excitation conditions of ULIRGs, an hypothesis already advanced in previous works but never confirmed or discarded (Cicone et al. 2018; Montoya Arroyave et al. 2023). We start by exploring the presence of any trends between the observed CO line widths and the upper rotational quantum number J_up.

Before showing the results, we discuss briefly what we expect from such analysis from a classical PDR perspective. Carbon monoxide is unique in being largely abundant and having low-J transitions that are easily excited under mild conditions (i.e., low-density and low-temperature environments), causing their emission to be commonly optically thick (Narayanan & Krumholz 2014). Conversely, mid- to high-J CO lines require the gas to be in more extreme conditions (either denser or warmer environments) to be excited. In the standard PDR scenario, such environments are typical of molecular cloud cores (i.e., high density gas) and/or active star-forming regions (high temperature and density) (Hollenbach & Tielens 1997). It follows that mid- and high-J CO lines may arise from regions with different kinematics: (i) dense cores of GMCs, where the gas is more shielded from processes that introduce turbulence to (or even disrupt) the gas cloud; and (ii) regions with active star formation, injecting both kinetic and thermal energy into the surrounding medium, leading to a layer of gas with turbulent motions emitting higher-J CO lines. Indeed, emission coming from the latter may broaden the spectral line profile due to turbulence; however, the amount of gas found in such regions is overall a small fraction of the total molecular gas budget in galaxies, and as such we would expect the contribution to the total integrated line profile to be minimal. In fact, in the classical PDR interpretation, high-J CO lines are expected to depart from thermalisation (Narayanan & Krumholz 2014), and are more likely to be optically thin, with line widths that are either narrower than the low-J CO lines, or at most consistent if we assume that all the gas shares the same galaxy-wide kinematics.

Since our APEX single-dish spectra enclose the majority of the ISM of these sources (beam values ranging from 8 − 27 arcsec, corresponding to ∼14 − 48 kpc), any observed CO spectral line broadening will be dominated by macroscopic processes. We rule out opacity (see e.g., Hacar et al. 2016) to produce any measurable line broadening effects on such scales, as demonstrated by the near 1:1 correspondence observed between the line widths of low-J CO lines and the optically thin atomic carbon lines measured in our sample, presented and discussed in Montoya Arroyave et al. (2023)¹³ Doppler broadening due to kpc-scale rotating molecular gas disks can produce line wings up to v ∼ 200 − 250 km s⁻¹ on each side of the line in galaxy-integrated spectra of ULIRGs. In the presence of massive molecular outflows, the total observed CO line spectra can extend up to absolute line-of-sight velocity values of ∼1000 − 1500 km s⁻¹ in the most extreme cases (see e.g., Cicone et al. 2014). In addition, inflows, tidal tails, and merging companions included in the beam can also contribute to the observed CO line broadening. However, due to the extremely high incidence of massive molecular outflows observed in ULIRGs (> 70%, see also OH outflow properties reported in Table 1 in Montoya Arroyave et al. 2023), it is reasonable to assume that high-velocity and broad CO line components in our sample contain also outflowing material.

The most direct measurement of CO line widths can be obtained through a single Gaussian fitting, by using the full width at half-maximum (FWHM). However, the total CO line profiles of local ULIRGs often present broad wings, asymmetries and/or double-peaked emission. In these cases, a single Gaussian component may be a poor fit. To account for this, we also estimate the 16th–84th percentile velocity interval (v₈₄ − v₁₆), derived from the analytical expression of the best-fit line profile obtained with the multi-Gaussian fit. In Fig. 4 we show the CO line widths, estimated with both methods (top and bottom panels), plotted as a function of J_up of the transition. The left panels show the whole sample, while the middle and right panels display separately the lower- and higher-excitation sources, respectively, as classified in Sect. 3.1. The sample shows a high variance, however, interestingly, when averaging the data (black solid line in the left panel of Fig. 4), we retrieve a mean increasing trend of CO line widths as a function of J_up, with a Pearson correlation coefficient of 0.91 and a p-value of 4.3 × 10⁻³ for the v₈₄ − v₁₆ line widths measurements (for the FWHM approach we retrieve 0.95 and 1.2 × 10⁻³, respectively), hence it is statistically significant at the 4σ level. The relations obtained with the two line width estimate methods show remarkably similar results. The middle panels of Fig. 4 show that for lower-excitation sources the line widths increase up to CO(J_up ∼ 3, 4), similar to the shape of the average CO SLED retrieved for these sources. For higher-excitation sources, the average trend shows increasing line widths up to the CO(6−5) transition, and a decrease for the CO(7−6) line (see e.g., Figs. B.5 and B.11 for sources with broader high-J CO lines). The behaviors for different sources indicate that the trend of broader line profiles for higher-J CO lines is consistent with being driven by the higher-excitation ULIRGs in our sample. To better visualize the broadening observed from low- to high-J CO transitions in the higher excitation ULIRGs, in Figs. B.1–B.17 we also show a plot overlaying the best-fit line profiles (normalized to their maximum value) obtained for all the available CO transitions for each source. The sources showing the most remarkable line broadening in higher-J CO transitions are: IRAS 06035−7102 (Fig. B.4), IRAS 06206−6315 (Fig. B.5), IRAS 08311−2459 (Fig. B.6), and IRAS 16090−0139 (Fig. B.11).

Fig. 4. Observed CO line widths (normalized by the CO(2−1) line width) as a function of the upper rotational quantum number Jup. Top row: CO line widths computed using the velocity percentile v84 − v16 interval, derived from the analytical expression of the best-fit line profile obtained with the multi-Gaussian fit. Bottom row: line widths obtained via a single Gaussian spectral fit to each CO emission line. Left column: all the sources in our sample; middle and right columns show the ULIRGs with low and high excitation CO SLEDs, respectively (see explanation in Sect. 3.1). The sources are color coded as presented in the upper legend, the mean trend is computed for each plot and shown with black square markers.

This result is remarkable considering that it is obtained from an analysis of single dish beam-integrated spectra, probing regions as large as ∼10 − 15 arcsec in the mid-J CO transitions, corresponding to tens of kpc. As such, the processes at the origin of the observed trends must be predominant on galactic scales and affect a significant portion of the molecular ISM in these sources, otherwise their signature would be highly diluted in the large beam. This result challenges classical nuclear (r < 1 kpc) PDR or XDR scenarios for the origin of the high CO excitation in these ULIRGs, as already, emerged from the results of the CO SLED analysis presented in Sect. 3.3.

Most previous galaxy-wide CO SLEDs surveys have been carried out with Herschel spectroscopic data (see Kamenetzky et al. 2014; Pearson et al. 2016; Liu et al. 2015; Greve et al. 2014; Lu et al. 2014; Pereira-Santaella et al. 2014; Mashian et al. 2015), and, while such studies provide a global view of the extreme physical conditions of the molecular gas in (U)LIRGs, they lack the sensitivity and spectral resolution needed to perform a spectral profile analysis. Consequently, only few literature studies have evaluated possible variations of the line widths of different CO rotational transitions. Most previous works found similar line widths throughout low- to high-J CO lines, and usually opted for deriving an average line width for all CO transitions for individual sources (e.g., Yang et al. 2017; Greve et al. 2005; Weiß et al. 2007). In the study of 40 SMGs at redshift z ∼ 1.2 − 4.1 performed by Bothwell et al. (2013), they found no significant differences between the line widths of low-J (i.e., J = 2 − 1) and higher-J (i.e., J_up ≥ 3) lines, with the mean values of the distributions agreeing to within their 1σ errors. From the FWHMs derived for each of their SMGs, Bothwell et al. (2013) found that “hotter dust” ULIRGs are kinematically similar to SMGs as traced by their line widths, with ULIRGs having on average narrower line profiles than SMGs. Boogaard et al. (2020) reached a similar conclusion based on a study of spectrally resolved CO(J_up ≤ 8) emission lines in a sample of SF galaxies at 0.5 < z < 3.6. Their measurements showed consistent widths for different CO rotational transitions. They argued that the gas distribution does not differ strongly in different transitions, leading to almost constant line widths at different J. The galaxies in their sample are, however, less extreme than our ULIRGs, with lower infrared luminosities and SFR surface densities than the SMGs and SF galaxies studied in, for example, Bothwell et al. (2013), Spilker et al. (2014), Daddi et al. (2015).

An interesting take on the diagnostic power of CO line widths is that of Rosenberg et al. (2015), who argued that the width of an emission line could indicate more mechanical energy available in the host galaxy, and, as such, more heating of the molecular gas. In their study, these authors combined ground-based low-J (J ≤ 3) CO observations¹⁴, with Herschel CO(4 ≥ J ≥ 13) data for 29 (U)LIRGs. They measured the FWHMs of the spectrally resolved CO(1−0) lines, and found that the lower-excitation sources display narrower CO(1−0) lines, while the higher-excitation galaxies span a wide range of CO(1−0) line widths. They propose that, by combining the information from the measured CO(1−0) line widths with the AGN luminosity, it is possible to shed light onto the origin of the high CO excitation. In particular, for the higher-excitation sources that cannot be explained only with UV heating from massive stars, these authors propose that broad lines may be indicative of mechanical heating, especially when the AGN contribution is low and so heating from X-rays and Cosmic Rays may not be important.

According to Bournaud et al. (2015), the high CO gas excitation measured in starburst mergers can be driven by strong turbulent compression motions associated with tidal interactions, which lead to high-density environments. This suggests that non-ordered motions in galaxy mergers cannot be neglected when studying molecular lines. However, the simulation used by Bournaud et al. (2015) does not include any prescription for AGN feedback, which is important in ULIRGs, and, more importantly, it does not account for cold molecular gas embedded in galactic outflows. The latter is now a well-established result, thanks to the significant observational efforts made in the past 10 years (see e.g., Veilleux et al. 2020), and it is even more relevant in ULIRGs where massive molecular outflows are almost ubiquitous and can be as massive as to embed a significant portion of their molecular ISM (Cicone et al. 2018).

4.2. A possible relation between CO excitation and molecular outflows

In this section we explore further the marginal trend of increasing line width in higher-J transitions presented in Sect. 4.1, with the aim of testing whether such relation can be linked more closely to massive molecular outflows. Since the positive relation observed in Fig. 4 is driven by the higher-excitation ULIRGs, here we focus only on these objects (8 out of the total sample of 17 ULIRGs).

Due to the turbulent nature of the ISM and the merging state of ULIRGs, it is not straightforward to link the broad/narrow components of galaxy-integrated CO line spectra to outflowing/non-outflowing gas, in absence of additional spatial information. Nevertheless, most of these ULIRGs are known to host powerful molecular outflows, which have been unambiguously detected through blue-shifted OH absorption lines (Spoon et al. 2013; Veilleux et al. 2013; Sturm et al. 2011). Therefore, as already discussed in Montoya Arroyave et al. (2023), where we performed a similar analysis limited to the three lowest-J CO transitions, we can assume that, in first approximation, the broad and/or high-v_cen components of the CO spectra are most likely dominated by molecular gas embedded in outflows, while the low-v_cen and low-σ_v components are most likely dominated by the gas in rotating disks and in other regions of the “non-outflowing” ISM. Following this line of reasoning, we performed a multi-Gaussian component simultaneous fit to all CO transitions available for each source, by tying v_cen and σ_v of each Gaussian to be equal in all CO transitions, allowing only their amplitudes to vary freely in the fit (for details see Montoya Arroyave et al. 2023). We allowed the fit to use up to three Gaussian functions for each line, in order to capture any possible asymmetries/wings in the observed spectra. Such simultaneous multi-Gaussian spectral fits, performed only for the higher-excitation sources, are displayed in Appendix B.

The next step is to classify the different Gaussian components. To do so, we plot them in the FWHM-v_cen parameters space, shown in Fig. 5, including all eight higher-excitation sources. The plot shows a clustering of spectral components with high S/N (small errorbars) at low-v_cen and low FWHM values, which are consistent with probing the dynamically quiescent (i.e., non-outflowing) gas. We then draw a region in the FWHM-v_cen space enclosed within −150 < v_cen [km s⁻¹]< 150 and FWHM < 400 km s⁻¹ (rectangle in Fig. 5), which delimits these components classified as quiescent. We classify the components lying outside this area as components containing outflows. We note that a similar approach was used by Cicone et al. (2018) to identify the outflow-dominated molecular line emission from NGC 6240. An important take-away point from Cicone et al. (2018) is that the broad wings are highly attenuated in the galaxy-integrated CO spectrum of NGC 6240, where the components at velocities above 500 km s⁻¹ show fluxes that are < 5% of the peak of the line, despite this source hosting one of the most massive and extended molecular outflows known. Indeed, the spatially resolved analysis of the CO data of NGC 6240, reported in the same paper, demonstrates that (i) along the eastern-western direction of expansion of the outflow, several kpc away from the known nuclear molecular disk(s) in NGC 6240 (Medling et al. 2019), the broad CO components can reach velocities beyond 1000 km s⁻¹, and at the same time still show a significant flux at low-to-zero line-of-sight velocity, and that (ii) as much as 60 ± 20% of the total H₂ reservoir in the central 6 kpc × 3 kpc region of NGC 6240 is embedded in the outflow. This means that the spectral signature of even the most extreme molecular outflows can be highly diluted in galaxy-integrated spectra of ULIRGs, where the detection of even faint (a few % of the peak) high-v/broad components can indicate a potentially extreme outflow event. Such aperture-dilution effects are one of the reasons why the CO-integrated galaxy spectra (almost) never show velocities as high as the maximal velocities seen in blue-shifted OH absorption lines detected against the FIR continuum, the other reason being that the OH lines generally probe the innermost parts of these outflows, close to the launching base (< 200 pc for the high-velocity components of OH, see e.g., González-Alfonso et al. 2017). Indeed, when exploring the relation between the OH outflow velocities (reported in Spoon et al. 2013; Veilleux et al. 2013) and gas excitation for our sample, we find no evidence of higher excitation for ULIRGs with measured larger OH outflow velocities. In fact, some sources classified as higher-excitation in our work, are only detected in absorption in OH, and thus lack outflow velocity measurements. This has been already discussed in Montoya Arroyave et al. (2023) for our parent sample of (U)LIRGs, where we argue that OH velocities are not representative of the strength of molecular outflows, as they depend on other factors (orientation of the outflow with respect to the Line of Sight, the presence and strength of the FIR continuum).

Fig. 5. FWHM as a function of the central velocity of all Gaussian components employed in the simultaneous CO spectral fits performed on the higher-excitation sources. The dashed blue rectangle encloses the region of the parameter space that we ascribe to the quiescent gas. The components falling outside this regions are classified as being dominated by outflowing gas.

Using the classification reported in Fig. 5, we computed velocity-integrated CO fluxes for each Gaussian component and, from these, we constructed CO SLEDs for all the quiescent and high-v/high-σ components, separately. The results of this analysis are shown in Fig. 6, where the top panel shows in blue color the CO SLEDs of the “quiescent/non-outflowing” components (i.e., data points lying within the rectangle in Fig. 6) and the bottom panel in green color shows the CO SLEDs of the high-v/high-σ components (i.e., data points lying outside of the rectangle in Fig. 6) that likely contain outflowing gas. The plots show a striking difference, where the CO ladders that are characteristic of the quiescent gas are significantly less excited than those produced by considering the high-v/high-σ gas components. This result, combined with the increase of line widths in higher-J CO transitions shown in Fig. 4, strongly suggests that the higher CO excitation seen in some ULIRGs is related to the appearance of prominent broad line wings, which likely contain massive molecular outflows in these sources. In other words, in those ULIRGs whose global CO SLEDs are characterized by a very high excitation, the source of such higher CO excitation is gas with a higher velocity shift and/or a higher velocity dispersion compared to disk gas, hence likely gas embedded in powerful, galaxy-wide outflows. Since our analysis is based on total (i.e., galaxy-integrated) line spectra, the trends shown in Fig. 6 are remarkable, as they signify that the outflows may be affecting at the same time the global gas kinematics and the excitation of the molecular ISM in these sources.

Fig. 6. Top panel: CO SLEDs for the spectral components that lie within the blue rectangle in Fig. 6 probing dynamically quiescent gas. Bottom panel: CO SLEDs for the spectral components that lie outside the blue rectangle in Fig. 6, likely containing molecular outflows in these sources. Only higher-excitation ULIRGs are included in this analysis (see Sect. 3.2). Quiescent components show significantly less excited CO SLEDs compared to the high-v/high-σ spectral components. The dotted gray lines show the thermalized optically thick limit (in all lines) SLED profiles, normalized to the CO(2−1) transition.

4.3. The origin of high CO excitation in ULIRGs

The source of the highly excited CO gas in (U)LIRGs has been a long standing issue. Previous studies of global CO SLEDs of ULIRGs used Herschel data, and so have been limited by the lack of spatial and/or spectral resolution, and by the use of heterogeneous datasets as the Herschel observations did not cover low-J CO transitions. Nonetheless, these works showed that a classical PDR scenario is unable to account for the highly excited mid-J CO SLEDs of ULIRGs (see Sect. 3 and references therein), and that other heating mechanisms are needed. While low-J CO lines are commonly optically thick and likely thermalized, populating high-J CO lines usually requires non-thermal excitation mechanisms such as shocks, cosmic rays, or X-ray heating (Esposito et al. 2022; Pereira-Santaella et al. 2013; Valentino et al. 2020), and these mechanisms have been suggested in those case where a single PDR component fails to reproduce the observed CO SLED (see e.g., Pereira-Santaella et al. 2014; Chen et al. 2023; Mingozzi et al. 2018; Pozzi et al. 2017). XDRs produced by the AGN are commonly invoked to account for higher-J CO fluxes in ULIRGs. While there have been studies showing that the chemistry and physical conditions of the gas can be strongly impacted by AGN, either due to X-rays or shocks (e.g., Aalto et al. 2011; Viti et al. 2014; García-Burillo et al. 2014; Usero et al. 2004), the consensus seems to be that emission associated with the AGN becomes significant only at J ≳ 9 CO lines (Brusa et al. 2018; Lu et al. 2017; Esposito et al. 2022; Vallini et al. 2019).

While we are not able to discern the exact origin of the CO excitation (e.g., radiative vs mechanical feedback) due to the lack of higher-J CO constraints, our results, shown in Figs. 4 and 6, challenge all scenarios for which the source of the high CO excitation in ULIRGs is a “classical” component of the ISM, such as PDRs formed around massive disk stars, or an XDR due to an AGN. Indeed, in such cases one would expect the highly excited gas to reside both in the disk and in the outflow, but possibly to affect more the former than the latter, as it is the disk gas that is first exposed to the radiation coming from massive stars and AGNs. Galactic-scale outflows develop from regions of massive star formation and/or from a central AGN, hence it is clear that radiation from these processes can leak through the outflows and so affect also the gas entrained at high velocities. However, this would not explain why the bulk of the kpc-scale outflowing material would be significantly more affected than the disk, which is the interpretation we favor from our analysis of galaxy-integrated spectra. Furthermore, because our analysis is based on galaxy-integrated single-dish spectroscopy, the observed trends cannot be driven by perturbed gas confined within the launching base of the outflow that is directly exposed to stellar and/or AGN radiation, as the emission from such component would be very compact, arising from scales of at most ∼100 pc, and so it would be diluted within our large beam. Even in the case of gas heated to extreme temperatures by, for instance, shocks, in such compact regions, the areas sampled by our beams of ∼10 kpc are typically 10 000 times larger than the areas defined by the central ∼100 pc of these sources. Hence, we infer that the culprit for the globally high CO excitation of these ULIRGs must be gas that is significantly dynamically disturbed, and at the same time occupies a substantial fraction of the beam of our observations, so that its luminosity contributes significantly to the observed spectra in all CO transitions, producing a detectable signature.

Previous studies by Meijerink et al. (2013) and Papadopoulos et al. (2014) focused on the galaxy NGC 6240 argue that the high CO excitation is arising from galaxy-wide shocks. While this is not exactly the same mechanism, it demonstrates that shocks can have indeed an impact on the shape of CO SLEDs and shocks are plausible to happen during the launch of molecular outflows. Thus, the main question is to determine how much of the gas is affected by such outflow-related processes, or, alternatively, how much of the source solid angle is occupied by these components, which would determine to what degree the global CO SLEDs will be impacted. To this end, it is then necessary to follow up with spatially resolved observations that allow for a proper characterization of the CO line luminosities in different galactic regions. In the context of previous results that point towards a mechanical feedback as the main culprit of high CO excitation in starburst galaxies, it is also worth mentioning the work of Harrington et al. (2021). Their analysis uses a more realistic approach to CO SLED modeling that, once applied to a sample of high-redshift lensed dusty star forming galaxies, reveals the need for mechanical energy to be injected onto scales larger than kpc to drive the observed CO ladders.

In conclusion, although we are not able to pinpoint the exact physical mechanism at the origin of such a link (e.g., heating due to outflow-driven shocks, large velocity gradients in the outflows, and radiation or cosmic rays propagating mainly through the outflows without affecting the quiescent gas) our analysis has reliably shown there is a large fraction of the molecular gas in these ULIRGs that is both highly excited and significantly dynamically disturbed, likely contained within galaxy-scale outflows, while the dynamically quiescent gas is significantly less excited. Although in our sample we did not detect a correlation between FIR colors and sources showing the most striking higher-J CO line broadening (as discussed in Sect. 3.3), our proposed link between high global CO excitation and massive molecular outflows is not in contradiction with previous studies linking the high CO excitation with warmer FIR colors such as Pearson et al. (2016). Indeed, one of the main mechanisms proposed for the acceleration of molecular outflows, especially in dusty galaxies hosting powerful AGN, is radiation pressure on dust (see e.g., Costa et al. 2018; Bieri et al. 2017; Ishibashi & Fabian 2023), which thus could explain both properties of ULIRGs.

We note that molecular outflows show different properties in different classes of galaxies for what concerns their CO excitation, morphology/extent and energetics, probably reflecting a combination of (i) different ISM conditions of their host galaxies, (ii) different evolutionary stages/ages of the outflows themselves, and (iii) different driving mechanisms. This field is still in its infancy, and we are only now starting to address these issues. The literature reports both cases of low-excitation (e.g., ESO 320-G020, see Pereira-Santaella et al. 2020) and high-excitation molecular outflows (e.g., IC 5063, see Dasyra et al. 2016). In this work, we are not proposing that all outflows (in ULIRGs or in other sources) have high excitation, but rather we are proposing a link between the exceptionally highly excited global CO conditions of ULIRGs (observed and noted since the early days of submillimeter astronomy) and the extremely high incidence (∼70%, see e.g., Lamperti et al. 2022; Veilleux et al. 2013) of massive molecular outflows in these sources. The latter is indeed much higher than the incidence observed in normal main sequence star forming galaxies (∼20 − 25%, see Stuber et al. 2021), which also show much lower global CO line excitation.

5. Summary and conclusions

We presented new, high S/N mid-J CO line observations (4 ≤ J_up ≤ 7) of 17 local ULIRGs obtained with the APEX telescope, using the nFLASH460 and SEPIA660 heterodyne receivers. We detected 32 out of the 33 targeted mid-J CO lines, adding up to 74 out of 75 CO(J ≤ 7) line detections. The only non-detection is due to an atmospheric line close to the tuning frequency, and most spectra have a high S/N that enable a detailed analysis of asymmetries and broad components of the spectral line profiles. For all these sources, we have very high S/N low-J CO line observations (CO(1−0), CO(2−1), CO(3−2)), obtained by combining proprietary and archival APEX and ACA/ALMA data, which were presented in a previous paper (Montoya Arroyave et al. 2023). Compared to previous multi-J CO studies of local ULIRGs that combined heterogeneous ground-based low-J CO data with high-frequency Herschel spectroscopy, our study can rely on a homogeneous, ground-based dataset, characterized by a high S/N and high spectral resolution. This survey was expressly designed to be sensitive to faint and diffuse components, such as those embedded in the massive molecular outflows that are almost ubiquitous in local ULIRGs, which may have been missed by previous interferometric and/or low S/N analyses. On the other hand, our study is limited in that it lacks high-J CO coverage, as there are no space-based facilities that can provide high-spectral resolution observations in the THz window. Our results can be summarized as follows:

• Our sample exhibits a variety of CO SLED shapes, with peaks or plateaus emerging at different J values, demonstrating that the ULIRG population can not be described by a single set of average physical conditions. Some sources show weakly subthermally excited lines up to CO(4−3), and in particular, galaxy IRAS 10378+1109 shows almost thermalized emission up to CO(7−6).

• However, in general, these ULIRGs are characterized by rather extreme CO excitation conditions, with average CO line ratios that are higher than those measured in high-z dusty star forming galaxies. We measure average ratios of ⟨r₄₁⟩ = 0.54 ± 0.03 (based on 9 sources) and ⟨r₆₁⟩ = 0.170 ± 0.007 (based on 6 sources), where $r_{J1}\equiv L^{\prime }{}_{\text{CO}(J\to J-1)}/L^{\prime }{}_{\text{CO}(\text{1-0})}$. The average r₅₁ and r₇₁ measurements, respectively ∼0.42 and ∼0.58, are biased towards higher values, since they are based only on two sources with high CO excitation.

• Based on the shape of their CO ladders, we classify our sample into lower-excitation sources, characterized by CO SLEDs that peak around J_up ∼ 3, 4, and higher-excitation sources, for which the CO SLEDs keep increasing or plateau up to the highest-J CO line detected.

• A single component LVG modeling of the CO SLEDs of our sample, performed following the method developed by Yang et al. (2017), delivers average volume densities n_H2 ≈ 10^2.0 − 10^5.8 cm⁻³, with the majority (14 out of 17) of sources having moderate densities with n_H2 < 10^3.3, and only three sources having high densities of n_H2 ≳ 10^4.5. Gas kinetic temperatures vary from T_kin ≈ 16 − 795 K, and CO column density per unit velocity gradient N_CO/dv ≈ 10^16.2 − 10^18.8 cm⁻² (km s⁻¹)⁻¹. These values are, however, likely to be overestimated, especially for the higher-excitation ULIRGs, for which extremely high (unphysical) values are derived for n_H2 and/or T_kin, in order to reproduce the observed CO fluxes, indicating the need of a second high excitation component. However, a two component LVG fit is highly degenerate with our data, because of the lack of CO(J_up > 7) lines, hence it is not presented in this work.

• We measured the line widths of each CO transition using both a Gaussian fit and a non parametric measurement, and analyzed their variations as a function of the CO transition (J_up). We found a positive trend, statistically significant at the 4σ level, of increasing CO line width with J_up of the transition, which appears to be driven by the higher-excitation sources in the sample.

• We then investigated further this trend, by focusing on the higher-excitation sources. Similar to the analysis already presented in Montoya Arroyave et al. (2023) for the low-J CO lines, we performed a simultaneous multi-Gaussian spectral fit, where all Gaussian components are constrained to share the same central velocity and FWHM in all CO transitions available for each source, while only their amplitude is allowed to vary. We then analyzed separately the CO SLEDs obtained from the spectral components displaying low central velocities and/or small FWHM values, and those characterized by disturbed gas kinematics, that is, with high velocity shifts and/or high FWHM values, which in these ULIRGs are likely dominated by massive molecular outflows. We found a clear difference in CO excitation between the CO ladders that are characteristic of dynamically quiescent gas are significantly less excited than those produced by considering the high velocity and/or high velocity dispersion gas.

• We find evidence, for the first time in a sample of local ULIRGs, that the higher CO excitation evidenced in some of these sources is related to the broadening of their CO lines. We suggest that the higher excitation gas is embedded in galaxy-scale molecular outflows, while it only contributes to a minor extent to the dynamically quiescent ISM of these sources.

Overall, we confirm the exceptional global molecular gas excitation of local ULIRGs, and we favor the interpretation that the origin of the higher excitation is associated to the large-scale outflows present in these sources, rather than the interpretation that the very hot components at the launching base of the outflows, or very hot turbulent gas confined in very small scales, are the drivers of the observed global physical conditions. Moreover, our results challenge scenarios for which the high excitation phase of the molecular medium in ULIRGs is related to a classical PDR or XDR component residing in the dynamically quiescent disk gas, while they suggest that such phase can be found in high-velocity gas, likely containing outflows, and filling the beam of single-dish observations.

This study demonstrates the need to perform extensive high S/N line surveys with sensitive large-aperture single-dish telescopes that can capture also the faint, most diffuse and extended components of the molecular ISM of galaxies including the outflowing gas, and so it provides a key scientific driver for the development of the future Atacama Large Aperture Submillimetre Telescope (AtLAST¹⁵, Klaassen et al. 2020; Ramasawmy et al. 2022; Mroczkowski et al. 2023) and its instrumentation. To perform the kind of analysis presented in this work, the latter needs to include high frequency (ν > 350 GHz), high spectral resolution instrumentation that can deliver stable baselines over a broad simultaneous bandwidth of at least ∼4 GHz, to enable the characterization of broad and high-velocity emission line components of extragalactic sources.

Acknowledgments

We want to give special thanks to the APEX and ESO staff for their efforts at continuing APEX operations and completing our projects despite the Covid19 pandemic. I.M.A. would like to thank Davide Decataldo and Alice Schimek for useful discussions. We thank the referee for valuable comments that helped improve the paper. 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#2013.1.00535.S, ADS/JAO.ALMA#2016.1.00140.S, ADS/JAO.ALMA#2016.1.00177.S, ADS/JAO.ALMA#2013.1.00180.S, ADS/JAO.ALMA#2018.1.00503.S, ADS/JAO.ALMA#2016.2.00006.S, ADS/JAO.ALMA#2017.1.01398.S, ADS/JAO.ALMA#2017.1.00297.S, ADS/JAO.ALMA#2013.1.00180.S, ADS/JAO.ALMA#2018.1.00888. 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).

References

Appendix A: Measured fluxes and luminosities

Tables A.1 and A.2 list the integrated line fluxes and luminosities, respectively, obtained from directly integrating the spectra within v ∈ ( − 1000, 1000) km s⁻¹, after setting a threshold of > 1.5σ for each channel (see Sect. 2.4).

Appendix B: Reduced spectra and CO SLEDs

Figures B.1 to B.17 show, for each source of our sample, the final reduced CO line spectra obtained with APEX, the observed CO SLEDs, and the results of the LVG modeling performed with the code RADEX_EMCEE developed by Yang et al. (2017), as explained in Section 3.3. Each spectrum is fitted a (multi)-Gaussian profile, following the description provided in Section 2.4. The total fit is plotted with a solid line, and the individual Gaussian components are shown in dashed lines, all in the corresponding color for that transition. For the higher-excitation sources we perform an additional simultaneous fit, as described in Section 4.2, and are also plotted on the spectra for the relevant sources in black lines, the total fit with a solid lines, and the individual components in dotted lines.

Fig. B.1. Top two rows: Integrated and continuum-subtracted CO emission lines for galaxy IRAS F01572+0009. The spectral binning and telescope are reported on the top-left corner of each spectrum. The multi-Gaussian fit performed individually for each spectral line, is shown in corresponding colors to each CO transition (in a solid line is the total, and the individual components are shown in dashed lines). For the higher-excitation ULIRGs, the simultaneous multi-Gaussian fit is also overplotted on each spectrum in dark gray color (in a dashed line is the total, and the individual components are shown in dotted lines). Bottom panels: On the top left are all the available spectral line profiles (fitted with multi-Gaussian profiles) normalized to their maximum value plotted together. This allows for better visualizing the variations in line widths for different CO rotational lines. On the bottom left is the CO SLED fitted with one gas component, where the maximum posterior probability is shown in a purple solid line, while the median value of the marginal probability distribution is shown with an orange dashed line. On the bottom right panel, and with the same colors, are the posterior probability distributions of the molecular gas density (nH2), gas temperature (Tkin) and CO column density per velocity (NCO/dv) parameters for the single and two-component fit.

Fig. B.2. Continued from Fig. B.1 but for source IRAS 03521+0028.

Fig. B.3. Continued from Fig. B.1 but for source IRAS F05189-2524.

Fig. B.4. Continued from Fig. B.1 but for source IRAS 06035-7102.

Fig. B.5. Continued from Fig. B.1 but for source IRAS 06206-6315.

Fig. B.6. Continued from Fig. B.1 but for source IRAS 08311-2459.

Fig. B.7. Continued from Fig. B.1 but for source IRAS 09022-3615.

Fig. B.8. Continued from Fig. B.1 but for source IRAS 10378+1109. The green data points in the bottom-left panel showing the CO SLED correspond to the fluxes of the systemic velocity components (with central velocity vcen ≈ 12 km s−1), resulting from the simultaneous fit. If considering only this systemic component, the CO SLED looks much less extreme, and in agreement with subthermal excitation conditions. Nonetheless, since the data analyzed in this work does not allow us to determine the origin of the velocity offset observed in the CO(7−6) spectrum, and for consistency with the other sources, we use the modeling results obtained by considering the CO SLED that uses the total velocity-integrated CO fluxes.

Fig. B.9. Continued from Fig. B.1 but for source IRAS F12112+0305.

Fig. B.10. Continued from Fig. B.1 but for source IRAS F13451+1232.

Fig. B.11. Continued from Fig. B.1 but for source IRAS 16090-0139.

Fig. B.12. Continued from Fig. B.1 but for source IRAS 17208-0014.

Fig. B.13. Continued from Fig. B.1 but for source IRAS 19254-7245.

Fig. B.14. Continued from Fig. B.1 but for source IRAS F20551-4250.

Fig. B.15. Continued from Fig. B.1 but for source IRAS F22491-1808.

Fig. B.16. Continued from Fig. B.1 but for source IRAS F23060+0505.

Fig. B.17. Continued from Fig. B.1 but for source IRAS F23128-5919.

In each CO SLED plot, the dark purple solid line corresponds to the CO SLED solution from the maximum posterior probability, while the dashed orange line shows the CO SLED model from the median value of the marginal probability distribution of each parameter. The fainter purple lines correspond to a random set of 100 models within the 1σ quartile. The posterior probability distributions of each parameter are shown by the density-contour plots and the purple histograms, with contours increasing in steps of 0.5σ.

All Tables

All Figures

Fig. 1. CO spectral line energy distributions of the ULIRGs in our sample. The blue circles with the error bars are the velocity-integrated line fluxes in units of Jy km s−1. The dotted gray lines show the thermalized, optically thick limit (in all lines), SLED profiles, normalized to the CO(2−1) transition (which is the lowest-J CO transition that is available for all sources). The source IDs are reported on the top-left corner of each panel.

In the text

Fig. 2. Top panel: observed CO SLEDs for the ULIRGs in our sample. The color coding is shown in the legend on top of the figure. Bottom panel: we compute the average CO SLED of our sample, and compare it to other well-studied galaxies, the starburst galaxy M 82 (Mashian et al. 2015), the nearby ULIRG NGC 6240 (hosting a dual AGN and a SB; Papadopoulos et al. 2014), and the Milky Way inner disk (Fixsen et al. 1999). We also show the average CO SLEDs for a sample of local (U)LIRGs studied by Papadopoulos et al. (2012), and a sample of high-z SMGs (Bothwell et al. 2013). In both panels, the dotted gray line corresponds to a constant brightness temperature in the Rayleigh-Jeans regime (i.e., S​ν ∼ ν2), which is the trend expected for thermally excited CO emission, in the optically thick limit for all lines. All line fluxes are in units of Jy km s−1, normalized by the CO(2−1) flux measured by APEX.

In the text

Fig. 3. CO SLEDs of the ULIRGs classified as having lower excitation (left panel) and higher excitation (right panel) gas, following the classification criteria described in Sect. 3.2. Each source is shown in a different color, according to the legend reported on the top of each panel. The average CO SLED of each group of galaxies is also computed and plotted using black square symbols connected by a black solid line.

In the text

Fig. 4. Observed CO line widths (normalized by the CO(2−1) line width) as a function of the upper rotational quantum number Jup. Top row: CO line widths computed using the velocity percentile v84 − v16 interval, derived from the analytical expression of the best-fit line profile obtained with the multi-Gaussian fit. Bottom row: line widths obtained via a single Gaussian spectral fit to each CO emission line. Left column: all the sources in our sample; middle and right columns show the ULIRGs with low and high excitation CO SLEDs, respectively (see explanation in Sect. 3.1). The sources are color coded as presented in the upper legend, the mean trend is computed for each plot and shown with black square markers.

In the text

	Fig. 5. FWHM as a function of the central velocity of all Gaussian components employed in the simultaneous CO spectral fits performed on the higher-excitation sources. The dashed blue rectangle encloses the region of the parameter space that we ascribe to the quiescent gas. The components falling outside this regions are classified as being dominated by outflowing gas.
In the text

Fig. 6. Top panel: CO SLEDs for the spectral components that lie within the blue rectangle in Fig. 6 probing dynamically quiescent gas. Bottom panel: CO SLEDs for the spectral components that lie outside the blue rectangle in Fig. 6, likely containing molecular outflows in these sources. Only higher-excitation ULIRGs are included in this analysis (see Sect. 3.2). Quiescent components show significantly less excited CO SLEDs compared to the high-v/high-σ spectral components. The dotted gray lines show the thermalized optically thick limit (in all lines) SLED profiles, normalized to the CO(2−1) transition.

In the text

Fig. B.1. Top two rows: Integrated and continuum-subtracted CO emission lines for galaxy IRAS F01572+0009. The spectral binning and telescope are reported on the top-left corner of each spectrum. The multi-Gaussian fit performed individually for each spectral line, is shown in corresponding colors to each CO transition (in a solid line is the total, and the individual components are shown in dashed lines). For the higher-excitation ULIRGs, the simultaneous multi-Gaussian fit is also overplotted on each spectrum in dark gray color (in a dashed line is the total, and the individual components are shown in dotted lines). Bottom panels: On the top left are all the available spectral line profiles (fitted with multi-Gaussian profiles) normalized to their maximum value plotted together. This allows for better visualizing the variations in line widths for different CO rotational lines. On the bottom left is the CO SLED fitted with one gas component, where the maximum posterior probability is shown in a purple solid line, while the median value of the marginal probability distribution is shown with an orange dashed line. On the bottom right panel, and with the same colors, are the posterior probability distributions of the molecular gas density (nH2), gas temperature (Tkin) and CO column density per velocity (NCO/dv) parameters for the single and two-component fit.

In the text

	Fig. B.2. Continued from Fig. B.1 but for source IRAS 03521+0028.
In the text

	Fig. B.3. Continued from Fig. B.1 but for source IRAS F05189-2524.
In the text

	Fig. B.4. Continued from Fig. B.1 but for source IRAS 06035-7102.
In the text

	Fig. B.5. Continued from Fig. B.1 but for source IRAS 06206-6315.
In the text

	Fig. B.6. Continued from Fig. B.1 but for source IRAS 08311-2459.
In the text

	Fig. B.7. Continued from Fig. B.1 but for source IRAS 09022-3615.
In the text

Fig. B.8. Continued from Fig. B.1 but for source IRAS 10378+1109. The green data points in the bottom-left panel showing the CO SLED correspond to the fluxes of the systemic velocity components (with central velocity vcen ≈ 12 km s−1), resulting from the simultaneous fit. If considering only this systemic component, the CO SLED looks much less extreme, and in agreement with subthermal excitation conditions. Nonetheless, since the data analyzed in this work does not allow us to determine the origin of the velocity offset observed in the CO(7−6) spectrum, and for consistency with the other sources, we use the modeling results obtained by considering the CO SLED that uses the total velocity-integrated CO fluxes.

In the text

	Fig. B.9. Continued from Fig. B.1 but for source IRAS F12112+0305.
In the text

	Fig. B.10. Continued from Fig. B.1 but for source IRAS F13451+1232.
In the text

	Fig. B.11. Continued from Fig. B.1 but for source IRAS 16090-0139.
In the text

	Fig. B.12. Continued from Fig. B.1 but for source IRAS 17208-0014.
In the text

	Fig. B.13. Continued from Fig. B.1 but for source IRAS 19254-7245.
In the text

	Fig. B.14. Continued from Fig. B.1 but for source IRAS F20551-4250.
In the text

	Fig. B.15. Continued from Fig. B.1 but for source IRAS F22491-1808.
In the text

	Fig. B.16. Continued from Fig. B.1 but for source IRAS F23060+0505.
In the text

	Fig. B.17. Continued from Fig. B.1 but for source IRAS F23128-5919.
In the text