Physics of ULIRGs with MUSE and ALMA: The PUMA project. IV. No tight relation between cold molecular outflow rates and AGN luminosities

We study molecular outﬂows in a sample of 25 nearby ( z < 0 . 17, d < 750Mpc) ultra-luminous infrared galaxy (ULIRG) systems (38 individual nuclei) as part of the Physics of ULIRGs with MUSE and ALMA (PUMA) survey, using ∼ 400pc (0.1–1.0


Sample
The PUMA sample consists of 25 nearby (z < 0.165, d < 750 Mpc) ULIRG systems (38 individual nuclei) with IR luminosities (8−1000 µm) in the range log L IR /L = 11.9−12.7.The individual nuclei have luminosities from log L IR /L < 10.5 to log L IR /L = 12.7, based on the relative contribution of the nuclei to the total ALMA continuum fluxes (see Pereira-Santaella et al. 2021).The sample has been selected to cover a range of interacting phases: 12 systems are classified as advanced mergers (with distance between the nuclei d nuclei < 1 kpc) and 13 systems are classifies as interacting (with d nuclei = 1.8−8.3kpc).The nuclei of the ULIRGs can be classified as AGN-dominated or starburst-(SB-) dominated based on the AGN contribution in the MIR (α AGN , see Perna et al. 2021): eight ULIRG systems are dominated by AGN (α AGN ≥ 0.5) and 17 by a SB (α AGN < 0.5).Based on the optical classification, our sample includes the following: nine Seyfert galaxies (two Seyfert 1 and seven Seyfert 2), eight HII and eight low ionisation nuclear emission regions (LINERs, Perna et al. 2021).In this paper we adopt a combined classification, in which we consider as AGN all nuclei classified as AGN either based on the MIR criterion (α AGN ≥ 0.5) or based on the optical classification (Seyfert 1 or Seyfert 2).According to this combined classification, 11/25 systems (14/38 individual nuclei) are classified as AGN, while the others are classified as SBs.
An overview of the sample properties is shown in Table 1. Figure 1 shows the IR luminosities and redshift of the sample.
Depending on the redshift of the targets, the CO(2-1) transition falls in Band 5 or Band 6.The synthesised beam full-width at half-maximum (FWHM) has been selected for each target so that it corresponds to a similar physical spatial resolution for all targets (∼400−500 pc).As a result, the synthesised beam FWHM are in the range 0.13−1.05 .The maximum recoverable scales are ∼10 times the beam FWHM.Table 2 presents the details of the observations: the beam FWHM, line sensitivity and total CO(2-1) fluxes.The data-reduction is described in detail in the PUMA II paper (Pereira-Santaella et al. 2021).We note that two targets (12071−0444 and 13451+1232) are not presented in that work because the ALMA observations were not available at the time of publication, but they are included in this work.The ALMA observations of 12071−0444 have been reduced following the same procedure used for the other sources.The second target, 13451+1232 (4C+12.50), is a radio AGN with a strong 230 GHz continuum dominated by synchrotron radiation (Pereira-Santaella et al. 2021).Given that for this source the continuum is strong enough, we decide to apply self-calibration to these data in order to increase the signal-to-noise.We apply five rounds of phase self-calibration and one round of amplitude self-calibration.Due to the steep slope of its continuum, we use a first-order polynomial to subtract the continuum from the CO(2-1) spectral window in the uv plane.For the rest of the data-reduction, we follow the same procedure as presented in Pereira-Santaella et al. (2021).
(3) Name of the nucleus.( 4) and (5) Coordinates of the nuclei derived from the ALMA 1.4 mm continuum.If the continuum is not detected, the coordinates are measured from the NIR or optical HST (see Perna et al. 2021).( 6) Infrared luminosity of the nuclei.For the interacting systems, the proportion of IR luminosity in each nucleus is based on their relative ALMA continuum fluxes (see Pereira-Santaella et al. 2021).( 7) Morphology of the system.I: interacting system with nuclear separation >1 kpc; M: advanced merger with nuclear separation <1 kpc.(8) Nuclear separation.(9) Nuclear activity classification based on optical spectroscopy (see Perna et al. 2021).
(10) Fraction of AGN contribution to L bol derived from the 30 µm to 15 µm flux ratio (see Perna et al. 2021).( 11) Combined nuclear activity classification used in this paper.Nuclei are classified as AGN either if α AGN ≥ 0.5 or if their optical classification is Seyfert (Sy 1 or Sy 2).( 12) Objects classified as compact obscured nuclei (CONs) using ratios of the equivalent widths of different polycyclic aromatic hydrocarbons (PAHs) based on the method by García-Bernete et al. ( 2022). ( ): For 14348−1447, there is evidence that the AGN is located in the SW nucleus based on high-angular resolution mid-IR imaging (Alonso-Herrero et al. 2016) . ( ): Even though this target (Arp 220) is classified as a merger (d nuclei = 0.37 kpc), we report the position of the two nuclei that have been identified thanks to the low redshift of this source.For consistency, we treat this target as a single nucleus in the rest of the paper.

CO(2-1) moment maps
We produced the maps of the CO(2-1) intensity (moment 0), velocity (moment 1), and velocity dispersion (moment 2) for our sample.Before producing the moment maps, we masked pixels with low signal-to-noise ratio (S/N) in each velocity channel, where the noise was estimated using the median absolute deviation (MAD).Specifically, for each velocity channel, we masked spaxels with S /N < 3.5−5 depending on the overall S/N of the observations.Moreover, we masked individual spaxels which show spurious emission applying the following procedure.For each spectral channel, we applied a Gaussian kernel with a size of three pixels to smooth the channel image and then we masked individual pixels with S /N < 0.4 in the smoothed image.In this way, for each velocity channel, we remove isolated pixels with spurious emission, which could bias the calculation of the moment maps.The moment 0, moment 1, and moment 2 maps were produced using the cubes obtained with this masking process, while the CO(2-1) peak map was obtained from the original data cubes.We define the zero velocity (and consequently the redshift) based on the moment 1 velocity at the position of the continuum peak.The CO(2-1) redshifts (reported in Table 2) are in good agreement with previously reported redshifts.Throughout this paper, we use the radio velocity definition.
In 09022−3615, we detect an increase in velocity dispersion south of the nucleus (distance ∼1 , equivalent to ∼1 kpc), which corresponds to the most blue-shifted velocities in the moment 1 map.The spectrum in this location shows two peaks.This could be due to a blue-shifted outflow (or inflow) in this location, or a cloud pushed by an outflow.An alternative explanation is the presence of a second very obscured nucleus, which is not detected in the ALMA millimetre continuum.As we do not see evidence of rotation at the position of the putative second nucleus, we consider more likely the former explanations.
The CO(2-1) peak maps of 10/38 nuclei show a dip in the centre, corresponding roughly to the position of the continuum peak (e.g., 00188−0856, 12112+0305 NE, 13120−5453, 13451+1232 N).This drop in the centre is compatible with an extreme central optical depth.In a future work, we will present ongoing ALMA observations of the optically thin 13 CO isotopologue to investigate this.

Spectro-astrometry
We perform a spectro-astrometry analysis to identify highvelocity molecular gas that is decoupled from the main rotation of the galaxy disk.This analysis consists in determining how the centroid position of the CO(2-1) emission changes as a function of velocity.We follow a similar methodology to the one used by García-Burillo et al. (2015) and Pereira-Santaella et al. (2018). 1 The channel maps for the rest of the sample can be found at https: //doi.org/10.5281/zenodo.7022665.
We binned together the velocity channels in order to achieve a minimum S/N of five, necessary to reliably determine the position of the peak of the emission.To determine the peak position in each binned channel, we first identify the spaxel with the highest flux.Then, we consider a region of 5 × 5 spaxels centred on the maximum spaxel and we perform a 2D Gaussian fit.In this way, we can determine the position at sub-pixel scales.In some targets, the CO(2-1) peak map presents a dip in the centre (see previous section) in the central velocity-channels.Thus, the position determined from the pixel with the maximum flux is not representative of the centroid of the emission.In these cases (07251−0248, 10190+1322, 13120−5453, 13451+1232, 19297−0406, 19542+1110, and 20414−1651), we determine the centroid emission by fitting a 2D Gaussian.The uncertainties on the centroids are calculated as ∆x = FWHM beam /(2 × S /N), where FWHM beam is the beam size and S/N is the signal-tonoise ratio of the binned channel (Condon 1997).The spectroastrometry analysis of 13120-5453 and 20100−4156 SE are shown in Figs. 3 and 4 as examples.
We perform a linear bisector fit of the centroids of the low-velocity channels, to determine the orientation of the kinematic major axis (see dashed line in Figs. 3, 4, and E.1).As low-velocity channels, we consider absolute velocities |v| < 300 km s −1 .If there are channels at |v| < 300 km s −1 whose position deviates significantly from the direction described by the other low-velocity channels (e.g., 00091−0738 N), we do not include them in the fit of the kinematic major axis.To determine the direction of the high-velocity gas, we perform a bisector fit to the high-velocity centroids (see dotted line in Figs. 3, 4, and E.1).In general, we consider as high-velocity the channels with |v| > 300 km s −1 , and we exclude channels that do not agree with the direction of the highest-velocity centroids or that follow the directions of the kinematic major axis.In some cases, the centroids of the blue-shifted and red-shifted high-velocity channels occupy a very similar region (e.g., 13451+1232 W or 14378−3651), which prevents us from determining the direction of the outflow axis.It is possible that the emission of the highvelocity gas is compact and unresolved, or that the outflow is pointing towards our line of sight, so that the blue-and red-sides overlap.In some cases, we observe some gas deviating from the kinematic major axis, however, because of its low relative velocity, ∼200−300 km s −1 , it is more likely due to a tidal tail than to an outflow (e.g., 00091−0738 N, 01572+0009).

Comparison of position angles of the molecular gas, stellar and ionised gas disks
From the fit of the spectro-astrometry low-velocity channels, we derive the position angles (PAs) of the kinematic major axis of the molecular gas disk (listed in Table 2).We obtain the PA of the major axis on the receding half of the galaxy, measured east of north (anticlockwise).We compare these PAs with the ones of the stellar and ionised gas (traced by Hα) disks presented in Perna et al. (2022).Figure 5 shows the absolute differences between PA(CO) and the PA derived from the stellar and ionised gas kinematics.
Overall, there is a general agreement between these measurements, with differences within ∼20 • , hence it is consistent with what observed in non-interacting galaxies (see Perna et al. 2022).We identify only six outliers with PA differences 20 • (01572+0009, 07251−0248 E, 09022−3615, 14348−1447 SW, 16090−0139, and 17208−0014).The PAs of CO are measured on smaller scales (∼1 kpc), compared to the scales used to derive the PAs of the stellar disk and ionised gas (∼5-10 kpc).Thus, A45, page 4 of 49   6) CO(2-1) redshift, calculated based on the velocity of the moment 1 map at the ALMA continuum peak position.( 7) CO(2-1) integrated flux.For the interacting systems for which it is not possible to separate the flux belonging to the two nuclei (07251-0248, 11095-0238, and 19297-0406), we report the total flux of the system, as well as the total CO(2-1) luminosity (Col.7) and FWHM (Col.8). ( )For the undetected nucleus 16156+014546 SE, we provide a 3σ upper limit, calculated based on the rms and a typical line FWHM of 300 km s −1 .( 8) CO(2-1) luminosity, calculated as , where ν rest is the line rest-frequency, D L the luminosity distance and z the redshift (Solomon et al. 1997).( 9) Molecular gas mass calculated as (Bolatto et al. 2013b) and α CO = 0.78 M /(K km s −1 pc 2 ) is the ULIRGs-like CO-to-H 2 conversion factor.(10) FWHM of the integrated CO(2-1) line, corrected for instrumental resolution.(11) Position angle (PA) of the kinematic major axis measured from the spectro-astrometry of the low-velocity CO channels.The PA is measured east of north (anticlockwise) for the receding half of the galaxy.
we expect to see some differences between the PAs, especially considering the fact that our targets are mergers or interacting systems, many of which do not show regular rotating disk (only 27% and <50% of nuclei in the PUMA sample show rotating disks in the ionised gas and stars, respectively, Perna et al. 2022).In summary, we find that in most of the targets the molecular gas disk has an orientation (PA) similar to the stellar and ionised gas disk.

Outflow velocity range definition
One of the main goals of this work is to identify and characterise high-velocity (>300 km s −1 ) outflows of cold molecular gas, which produce broad wings in the line profile.Moreover, outflowing gas does not have to follow the disk rotation, thus it can be identified as high velocity gas that does not follow the main rotation pattern of the galaxy.For each nucleus, we define the velocity range where a potential outflow is located, using the spectro-astrometry analysis and the integrated CO(2-1) spectrum.We define the minimum (v min ) and maximum (v max ) velocities of the outflow (separately for the blue and red part of the outflow) and consider the emitting gas in the [v min , v max ] velocity range as part of the outflow.
To define v min and v max for each nucleus, we perform the following procedure, separately for the blue-shifted and red-shifted emission.We used the spectro-astrometry plot to select the minimum velocity at which the centroid position of the gas starts to deviate from the direction of the major axis of rotation.In A45, page 5 of 49  1).The magenta ellipse shows the FWHM and position angle of the ALMA beam .The contours in the maps are: continuum map: [0.3, 0.4, 0.6, 0. 8, 0.9] of the maximum; peak map: 1.5 × σ (where σ is the rms) and [0.1, 0. 2, 0.4, 0.6, 0.8] of the maximum; moment 0: [3,6,25,50,75] × σ (where σ is the rms), moment 1: every 50 km s −1 (every 25 km s −1 if the maximum value <100 km s −1 ), moment 2: every 25 km s −1 (every 15 km s −1 if the maximum value <150 km s −1 ).
particular, we look for velocity channels whose centroid position deviates significantly from the direction of the kinematic major axis, or that does not follow the rotation pattern (from blue to red).The v min values are in the range |v min | = 180−450 km s −1 .
To define v max , we started from the channel corresponding to v min and we continued adding velocity channels to create the map of the high-velocity emission, until the peak S/N of the integrated map starts to decrease.We checked by looking at the integrated CO(2-1) spectrum that we are not missing significant emission at v > v max due to a particular low S/N channel.The v max values A45, page 6 of 49 One caveat of our analysis is that since we are observing projected velocities, we are not sensitive to high velocity gas in the plane of the sky.Additionally, with our method we are not considering outflowing gas with low projected velocities, that is with velocities v such that v min (blue) < v < v min (red), because it overlaps with the velocities of the rotating disk.In order to identify this gas, we would need to model the rotation of the system to identify non-rotating gas (e.g., Brusa et al. 2018;Gao et al. 2021;Ramos Almeida et al. 2022).We plan to investigate this in a future work.

Measurements of the outflow properties and energetics
In this section, we describe the method we use to measure the main outflow parameters: outflow radius (R out ), outflow velocity (v out ), and molecular gas mass in the outflow (M out ).We use these parameters to derive the outflow energetics: mass outflow rate ( Ṁout ), mass-loading factor (η = Ṁout /SFR), outflow momentum rate ( Ṗout ), and kinetic luminosity (L out ).In the following, we explain how we measure the 'projected' R out and v out .We discuss the inclination corrections in Sect.4.5.1.
Different methods have been used in the literature to separate the outflow and rotating disk emission.A possible method (used for example by Pereira-Santaella et al. 2018) consists in subtracting the flux belonging to the rotating disk, by fitting the central velocity channels with one or two Gaussians and then considering only the flux in the residuals as part of the outflow.However, this method may underestimate the outflow mass as outflow flux with low (projected) velocities is generally assigned to the disk regardless of the actual position of the emission.An alternative method used in the literature consists in fitting the line profile using a systemic Gaussian component and a broader Gaussian component for the outflow (e.g., Fluetsch et al. 2019).This method may be overestimating the flux of the outflowing gas, since it considers that a large portion of the outflowing gas is at low velocities.
Since the resolution of the observations allows us to determine the position of the gas and to identify the gas that is not following the galaxy rotation, we prefer to consider only the gas with high velocities as part of the outflow.Moreover, most of the A45, page 7 of 49 line profiles of our targets cannot be well fitted using a simple model with only one or two Gaussians (see Figs. 3, 4 and E.1).The line profiles are asymmetric and show multiple peaks, which could also include self-absorption (see Sect. 4.1).Thus, to measure the outflow gas mass, we consider the total flux in the highvelocity channels, highlighted in the blue and red shaded regions on the spectra (see Fig. 3, previous section), without subtracting the low-velocity Gaussian fit.In Sect.4.5.2,we discuss the systematic effects affecting the derived outflow quantities depending on the different methods.A comparison of our outflow parameters with the ones reported by Fluetsch et al. (2019) and Lutz et al. (2020) is shown in Sect.4.6.
Outflow size: to measure the outflow radius, we fit a 2D Gaussian model to the high-velocity maps, obtained integrating the flux over the high-velocity channels (see bottom left panel in Fig. 3), separately for the blue and red part.To take into account the beam size and obtain the 'intrinsic radius', we convolve our model with a 2D Gaussian with the shape of the ALMA beam, and we fit this 'convolved model' to the maps.The radius of the blue (red) part of the outflow is defined as: where d b c is the centroid distance from the nucleus and FWHM b is the average size (deconvolved from the beam) of the two axes of the 2D Gaussian fit.
The outflow radius R out is the mean of the radii of the blueand red-shifted wings: The outflow radii are shown in Fig. 3 as dashed circles (bottom left panel).This method is analogous to the method used by Lutz et al. (2020), although applied here to higher spatial resolution data (400 pc vs. 700 pc), which is enough to resolve the outflow structure.Due to the limited S/N of the single channels, it is not possible to measure the radius in each channel, which would give a more accurate measurement to derive the mass outflow rate.We also measure the maximum extent of the outflow, by taking the maximum distance from the continuum position reached by the 3σ contour.We subtract the beam size in quadrature to obtain the 'intrinsic' distance.Then, we take the mean between the radius of the red-and blue-shifted channels as the representative maximum radius of the outflow (R 3σ ).
The outflow radii R out measured from the 2D Gaussian fit are in the range 0.18−0.94kpc (0.1−1.5 ).The maximum outflow radii R 3σ are in the range 0.1−2.1 kpc (0.05−2.5 ).The ratio of the observed R 3σ /R out is in the range 1-3, with a mean of 1.45.We note that in some cases the R 3σ values reported in Table 3 are smaller than R out .This is due to the different way used to deconvolve the beam and to the uncertainty on the Gaussian fit.We calculate that >50% of the outflow flux is within R out , with a median of 76%.Given that most of the outflow flux is located within R out , we decide to use R out to compute the mass outflow rate and the energetics.If we were to use the outflow flux A45, page 8 of 49  Notes.
(2) Name of the nucleus.( 3) and (4) Velocity range considered to measure the blue-and red-shifted wings of the CO(2-1) profile with respect to the systemic velocity.( 5) and ( 6) CO(2-1) flux in the blue-and red-shifted wings.(7) Maximum extent of the outflow estimated from the emission above 3σ (see Sect. 4.5).(8) Outflow radius (of the blue-and red-shifted wings) deconvolved from the beam, but not corrected for inclination.(9) Flux-weighted outflow velocity (see Eq. ( 5)), not corrected for inclination.(10) Outflow molecular gas mass, calculated assuming a ULIRG-like conversion factor α CO of 0.78 M /(K km s −1 pc −2 ) −1 and r 21 ratio of 0.91 (Bolatto et al. 2013b).( 11) Mass outflow rate calculated using Eq. ( 6). ( 12) Angle between the outflow axis and the kinematic major axis, derived from the spectro-astrometry, for the cases where it could be determined. (* ) Average outflow properties, excluding upper limits.The average of R out and v out have been corrected for inclination assuming the average corrections: within R 3σ instead of within R out to calculate the outflow mass, M out would increase on average by a factor of 1.3 (0.11 dex).
The mass outflow rate is proportional to M out /R out , thus the larger molecular mass included within R 3σ is counterbalanced by the larger radius.If we were to use R 3σ and the outflow flux within this radius, we would obtain very similar values (less than 8% smaller, −0.04 dex) compared to the Ṁout estimated using R out .
We do not attempt to correct the radii of the single targets for inclination, since information about the inclination is not available for the full sample (see discussion in Sect.4.5.1).
Outflow mass: to derive the outflow mass, we extract the central spectrum from a radius equal to the observed R out (not deconvolved from the beam) and we integrate the flux in the high-velocity channels between v min and v max , separately for the blue and the red part of the outflow.The central spectra are shown in Figs. 3, 4 and E.1 (bottom row, middle panel).Then, we sum the flux of the blue and red part of the outflow to obtain the total outflow flux (F out ).
We estimate the uncertainties on F out by extracting a spectrum from a region away from the source and measuring the standard deviation of the flux density in the high-velocity channels (σ in units of mJy) .The uncertainty on F out is : where ∆v ch is the width of a velocity channel in km s −1 , and N ch is the number of velocity channels in the high-velocity windows.We transform the CO(2-1) flux into luminosity (in units of K km s −1 pc 2 ) using the formula: where S CO is the velocity-integrated CO line flux in Jy km s −1 , ν rest is the line rest-frequency in GHz, D L is the luminosity distance in Mpc, and z is the redshift (Solomon et al. 1997).We A45, page 9 of 49 convert the CO(2-1) luminosity to CO(1-0) luminosity using r 21 = L CO(2-1) /L CO(1-0) = 0.91 (Bolatto et al. 2013b).Then we multiply it by the ULIRGs-like CO-to-H 2 conversion factor α = 0.78 M /(K km s −1 pc 2 ), to obtain the outflow molecular (H 2 ) gas mass M out .Although, we note that the cold molecular gas conditions in the outflow likely differ from those in the disk and, therefore, the CO-to-H 2 outflow conversion factor is uncertain (see e.g., Pereira-Santaella et al. 2020).
Mean outflow velocity: we calculate the mean velocity of the outflow separately for the blue-and red-shifted high-velocity wings, by taking the flux-weighted mean of the velocity in the channels identified as part of the blue-shifted (or red-shifted) emission (see Sect. 4.4, Fig. 3, middle panel of the bottom row): where v i is the velocity of channel i and F i is the CO(2-1) flux density in that channel.Different methods have been used in the literature to estimate the outflow velocity (see Sect. 4.5.2).We decide to use this 'flux-weighted velocity' to calculate the mass outflow rate, because it is independent from the modelling of the emission line profile and it gives more weight to the velocities at which most of the emission takes place.We note that the outflow velocities measured with this method are sensitive to the choice of the velocity window defined as 'high-velocity gas'.In Sects.4.5.2 and 4.6 we explore this possible bias.
Mass outflow rate: for the red and blue part of the outflow separately, we calculate the mass outflow rate (in units of [M yr −1 ]) using the formula: where R out is the outflow radius, M i is the H 2 gas mass in the channel with velocity v i , and the sum is over the velocity channels identified as part of the blue-shifted (or red-shifted) emission (see Fig. 3).The total mass outflow rate is the sum of the blue and red Ṁout .We note that this formula corresponds to the assumption that the outflow has started at a point in the past (−t = −R out /v out ) and has continued with a constant Ṁout (Rupke et al. 2005b;Veilleux et al. 2005;Lutz et al. 2020).Under this assumption, the average volume density of the outflowing gas decreases with radius (∝R −2 ).Assuming that the outflowing gas fills a spherical or multiconical volume with a constant average volume density, would increase Ṁout by a factor of three (e.g., Maiolino et al. 2012;Cicone et al. 2014;Lutz et al. 2020).
In the cases where the outflow is not detected, we estimate 3σ upper limits on Ṁout as: where v out = 390 km s −1 and R out = 0.52 kpc are the median outflow velocity and radius of our sample.M out,err is calculated based on F out,err (Eq.( 3)), measured from the spectrum extracted from a radius R out .The mass outflow rates are in the range ∼5 to ∼300 M yr −1 .We compare the outflow properties (v out , R out and M out ) measured from the blue and red-shifted wings.The v out , R out and M out measured from the red and blue parts are similar within a factor of 1.6, 2.4 and 1.2 respectively.The mean and corresponding standard deviation of the ratio of the red-and blue-shifted outflow quantities are: 0.97 ± 0.16, 1.05 ± 0.51, 0.99 ± 0.04, respectively.

Inclination corrections
We do not attempt to correct the outflow radii and velocities of the single targets for the inclination, since this information is not available for all objects.In order to apply an inclination correction to v out and R out , we need to know the inclination of the outflow with respect to our line of sight.Even if we knew the inclination of the molecular gas disk, in order to use this information we would need to know the inclination of the outflow with respect to the disk.In Sect.5.2 we discuss the projected angles between the outflow axis and the major kinematics axis.For a handful of targets ( 6), there is evidence that the outflow may be perpendicular to the major kinematics axis.However, only for one of them we have a measurement of the ionised gas (or stellar) disk inclination from Perna et al. (2022).Thus, we decide not to attempt to correct for inclination v out and R out of the single targets, and all the quantities reported here are the 'projected' ones.
However, we derive an average inclination correction that we apply to the average outflow properties of the sample.To convert the observed (projected) mean outflow velocity to intrinsic velocity, we need to divide v out by sin(i), where i is the inclination.Analogously, the observed outflow radius needs to be divided by cos(i) to recover the intrinsic value.Following Law et al. (2009), who considered the average for a collection of objects oriented isotropically in space, the average correction for the velocities is 1/ sin(i) = 1/0.79= 1.27.Analogously, we calculated the average correction for the radius: 1/ cos(i) = 2.We use these values to correct the mean outflow velocity and mean outflow radius reported in Table 3.
It is not possible to derive the average inclination correction for the mass outflow rate Ṁout in a similar way, since the calculation of the integral over the entire solid angle gives infinity.However, the average inclination correction for the dynamical time (t dyn = R out /v out ∝ M out / Ṁout ) is unity (see Cicone et al. 2015).

Caveats: Methods to measure outflow parameters
In this section, we discuss how the outflow parameters (in particular v out , M out and Ṁout ) change depending on the method adopted for the measurements.Readers who are less interested in the details about the methodology and the comparison with previous works may wish to go directly to Sect.4.7.
1) Low projected velocity (v < 300 km s −1 ) gas in the outflow: since we are not considering low projected velocities (|v| 200 − 300 km s −1 ), it is possible that we are missing part of the outflow flux.If we were to consider also this low velocities, the outflow mass M out would increase and the flux-weighted outflow velocity v out would decrease.Overall, we expect Ṁout to increase, but by a lower factor than M out , because the increase in M out is counterbalanced by the decrease of v out .
To test how much this effect could affect our measurements of the outflow properties, we consider 3D models of biconical outflows based on the models presented by Bae & Woo (2016).The outflow is a bicone with a half opening angle of 40 • .We set the maximum velocity of the outflow to 750 km s −1 .We consider several outflow radial velocity profiles, motivated by A45, page 10 of 49 previous works in the literature (e.g., Förster Schreiber et al. 2014;Venturi et al. 2018;Meena et al. 2021), and different outflow inclinations with respect to the line of sight.More details on the simulations can be found in the Appendix A. We measure v out , M out and Ṁout from the total simulated profile (out-flow+systemic component) considering only |v| > 300 km s −1 , to mimic the method we are using with our data.Then, we measure the outflow quantities (v out , M out , and Ṁout ) from the simulated outflow emission profile (without the systemic component) considering the full velocity range.The M out measured from |v| > 300 km s −1 are underestimated compared to the values measured from the entire velocity range by a factor of 0.2−1 (average 0.5), while the v out are overestimated by up to a factor of 2.2 (average of 1.6).Consequently, Ṁout would be underestimated by up to ∼0.45 dex (65%) for outflow inclinations close to the plane of the sky (90 • ).For an inclination of 10 • , the measured Ṁout would be underestimated by up to 0.1 dex (see Fig. 6).
For targets with a wide CO(2-1) line core (large FWHM), the flux of the gas in the rotating disk can overlap with the outflow flux up to larger velocities.Indeed, there is a positive correlation between the velocity at which we start to consider the flux to be dominated by the outflow (v min ) and the FWHM of the line (Spearmann rank correlation coefficient r = 0.81, p-value < 0.1).We did a test to estimate how much flux we may be missing in our measurements of the outflow flux for the targets with large FWHM CO(2-1) line profiles.In particular, we consider the 13 targets with |v min | > 300 km s −1 (either in the blue or red side).To estimate the amount of possible flux belonging to the outflow in the velocity range between |v| = 300 km s −1 and |v min |, we assume the outflow flux density in this range to be equal to the value at v min .Using this assumption, we estimate the outflow parameters (M out , v out , and Ṁout ) starting from v min = 300 km s −1 .We find that the value of M out increases by 0.28 dex on average (maximum 0.67 dex), while v out decreases by a factor of −0.07 dex on average (minimum −0.11 dex).So, the Ṁout estimates increase by 0.2 dex on average (maximum 0.6 dex).
2) Possible overestimation of Ṁout due to the rotating disk contribution at |v| = 250−300 km s −1 : since we do not model and subtract the disk rotation, it is possible that at low velocities (250−350 km s −1 ) we are including in the outflow flux some flux emitted by the gas in the rotating galaxy disk.To test how large this contribution could be, we subtract from the spectra the flux due to rotation estimated by modelling the core of the emission profile (with absolute velocities smaller than ∼300 km s −1 ) with one, two or three Gaussians (e.g., Pereira-Santaella et al. 2018).Then, we compute the outflow parameters (M out , v out , Ṁout ) from the residuals, considering the velocity range between v max and the velocity (v min ) at which the residuals approach zero.The outflow masses vary in the range −1.0 dex to 0.5 dex (−0.12 dex on average).In some cases, the measured M out increases because we can extend the outflow velocity range to smaller velocities, since there is no risk of including flux from the rotation.The outflow velocities vary by less than ±50 km s −1 .The Ṁout vary between −0.82 dex and +0.45 dex (−0.11 dex on average).
3) Different methods to estimate the outflow velocity: in this work, we adopted the 'flux-weighted' outflow velocity definition to compute the mass outflow rate.Other works instead have used different definitions of 'maximum outflow velocity' (e.g., Fluetsch et al. 2019;Lutz et al. 2020).If we assume a biconical outflow with constant gas velocity within the outflow, the range of velocities observed in the broad wings of the CO profile would be solely due to different orientation angles of the outflow gas clouds with respect to our line of sight.The part of the outflow closer to our line of sight would have the highest observed velocity.Thus, one could assume that this maximum velocity is the closest to the true/intrinsic velocity.We do not know the true radial profile of the velocity or the geometry of the outflowing gas for our targets.However, we can estimate how much our 'flux-weighted outflow velocity' differ from the 'maximum outflow velocity' assumption in our sample.We consider two definition of the maximum velocity.Both definition require the fit of the line profile with multiple Gaussian components: a broad component for the outflow, and other components for the systemic emission.The first definition that we consider is the prescription from Rupke et al. (2005a), used for example by Fluetsch et al. (2019): where |∆v| is the shift of the broad Gaussian outflow component with respect to systemic velocity and FWHM broad its full width at a half maximum.The second definition, adopted for example by Lutz et al. (2020), is: where FW 10% is the full width of the broad component at a tenth of its peak.We measure v max,FWHM and v max,FW10 for our sample of galaxies with outflow detection.We use up to three Gaussian components to model the systemic emission.We stress that in most of the cases, the parameters of the Gaussian components are highly degenerate.Thus, the derived maximum velocity can vary depending on the assumptions used in the fit.In Left: ratio between the maximum outflow velocity v max,FWHM (Eq.( 8)) and the flux-weighted outflow velocity (v out ) in our sample.
Fig. 7 we compare the flux-weighted v out with the maximum velocity v max,FWHM and v max,FW10 for our sample.The ratio v max,FWHM /v out varies between 0.5 and 2.0, with a median 0.9.
The ratio v max,FW10 /v out instead is almost always larger than one, with a maximum of 3.2 and a median of 1.6.Assuming that v max is the closest measure of the 'true' outflow velocity, it does not need to be corrected for inclination, while the observed (projected) v out would need to be corrected by an average factor of 1.3 (see Sect. 4.5.1).This factor can account for most of the average difference between v max,FW10 and v out .
We decide to avoid the methods based on the fit of the emission line profile with multiple Gaussian components, because the components are degenerate and this will introduce further uncertainties in the measurements.
An additional caveat is related to the definition of v max , which will impact the v out estimates.In selecting v max based on the S/N of the integrated high-velocity maps, it is possible that we are excluding diffuse high-velocity flux which is below the 3σ level.In this way, we may underestimate v out , especially for cases in which the outflow velocity increases radially.Unfortunately, with our current data it is not possible to estimate the impact of this possible additional component.Higher sensitivity observations are needed for this purpose.
In Table 4, we summarise the average biases in the outflow properties (v out , M out and Ṁout ) due to the effects described in 1), 2) and 3).Since the two effects described in 1) and 2) go in opposite directions, they tend to compensate each other.Taking into account the two effects, we may be underestimating M out ∼ 0.2 dex and Ṁout by ∼0.04 dex on average.
Based on the variations of outflow quantities from the tests using different methods, we estimate the typical uncertainties on the outflow quantities.For the outflow velocity, we estimate a typical uncertainty of 0.1 dex, while for M out and Ṁout we adopt a typical uncertainty of 0.3 dex.For R out , we estimate a typical uncertainty of 0.2 dex, based on the difference between R out and R 3σ .
The uncertainties on Ṁout are dominated by the uncertainties on the α CO conversion factor, which can be up to 0.7 dex (Papadopoulos et al. 2012;Bolatto et al. 2013b;Pereira-Santaella et al. 2020), and on the outflow geometry.If the α CO is more similar to the Galactic value (α CO = 4.3 M /(K km s −1 pc 2 )), it would imply that all our M out and Ṁout measurements are underestimated by a factor of ∼5 (0.7 dex), while if the optically thin case applies (α CO =

Comparison with previous works
In this Section we compare the derived outflow parameters with previous works that study properties of molecular outflow in nearby ULIRGs using CO observations.

Comparison with Lutz et al. (2020)
Lutz et al. ( 2020) study the outflow properties in a sample of 54 nearby (z < 0.2) galaxies with CO(1-0), CO(2-1) or CO(3-2) observations, with a range of spatial resolutions (30 pc to 5 kpc; 0.5−5 arcsec beam FWHM).They collect 41 nearby galaxies with molecular outflow detections or upper limits from the literature.They also present new NOEMA and ALMA observations, with a spatial resolution of ∼700 pc (0.5-3.7 arcsec beam FWHM), for 13 compact far-infrared galaxies from the Lutz et al. (2016) sample.To derive the outflow velocity and flux, they fitted the line profile with two Gaussian components (systemic and broad (outflow) component).They defined the outflow velocity using Eq. ( 9).To derive the outflow flux, they integrated the broad component only over the velocity ranges for which it contributes at least 50% of the total flux density of the line profile.The outflow radius was defined similarly to our method: R out = |∆R| + FWHM/2, where |∆R| is the distance of the outflow emission centroid from the continuum position and FWHM was derived from a Gaussian spatial fit in the uv-plane using a velocity range that is dominated by outflow, although, the spatial resolution is a factor of ∼2 lower than in our work.For the other 41 targets, they collect information about the outflow parameters (v out , R out , M out ) from the data published in the literature, trying to be consistent with their adopted definition of these parameters and their adopted methodology to separate the flux of the line core and high velocity wings.
There are 12 objects in common with our sample, with data published by Cicone et al. (2014) Barcos-Muñoz et al. (2018), Gowardhan et al. (2018), Pereira-Santaella et al. (2018), and Lutz et al. (2020).We compare our measurements of the outflow parameters (v out , R out , M out , Ṁout ) with these works in Fig. 8.There are some discrepancy between our measurements and the literature values.In particular for two galaxies (17208−0014 and 20100−4156) the literature v out is larger than 900 km s −1 , while we measure v out < 500 km s −1 .For 20100−4156, it is possible that the spectral range of our observations (v = [−1200, 1200] km s −1 ) is not wide enough to detect the emission at very high velocities (|v out | > 1000 km s −1 ).Additionally, this high-velocity outflow is detected in CO(1-0), while in CO(3-2) no outflow is detected (Gowardhan et al. 2018).Thus, it is possible that this difference also depends on the J-level observed.For 17208−0014, the S/N of the outflow component presented in Lutz et al. (2020) is not very high, thus the uncertainty on v out should be large (even though it is not reported in the paper).For four galaxies, R out values from the literature are considerably higher than our measurements (∼2−6 times higher).For three of these galaxies, the difference is due to the different definition of R out as the maximum radius at which the outflow is detected (see R max definition in Pereira-Santaella et al. 2018).
For the other target (20100−4156), the difference could be due to the different spatial resolution (1.5 arcsec vs. 0.2 arcsec).For Arp 220 instead R out from the literature is smaller than our measurement; also in this case the difference could be due to the different spatial resolution (0.09 arcsec vs. 1.0 arcsec).The M out values agree within a factor of 2.5, with some of our values being smaller and other larger than the values from Lutz et al. (2020).Since Ṁout is proportional to v out and R −1 out , the differences in R out and v out tend to balance each other and lead to similar Ṁout (within a factor of 3).

Comparison with Fluetsch et al. (2019)
We also compared the M out and Ṁout of our sample with the values reported by Fluetsch et al. (2019) for a sample of 45 local (z < 0.2) star-forming galaxies with CO(1-0), CO(2-1), or CO(3-2) ALMA observations.Fluetsch et al. ( 2019) used a different method to estimate the outflow mass: they measured the outflow flux by fitting the line profile with two Gaussians, one for the core of the line and one broad component for the outflow, and they considered the total flux of the broad component as the outflow flux.We expect that the outflow fluxes measured in this way will be higher than the ones measured with our method, since in addition to the flux in the wings, they are also considering the low-velocity emission of the broad component as part of the outflow.They measured the outflow velocity using Eq. ( 8).They fitted a 2D-Gaussian profile to the wing maps and used the beam-deconvolved major axis (FWHM) divided by two as the radius of the outflow.Compared to our methodology, they did not include the distance between the centroids of the blueshifted and red-shifted emission in the calculation of R out , thus their R out estimates are expected to be smaller.
Figure 9 shows the comparison of the outflow parameters derived by Fluetsch et al. (2019) with our results for the five sources in common between the two samples.Their v out tend to be higher than ours (maximum by factor of 1.7) , while their R out tend to be smaller (maximum by a factor of 1/4) .The largest difference is in the outflow masses, that are larger by up to a factor of 16 (1.2dex).This difference can be attributed to the different method used to estimate the flux belonging to the outflowing gas.This difference in M out propagates to the Ṁout .
Since the overlap between the two samples is small, we decide to also compare the outflow properties of the two samples at the same IR luminosity.Figure 10 shows M out and Ṁout as a function of IR luminosity for the Fluetsch et al. (2019) sample and our sample.The diamond symbols with red borders show the average values of the two samples for different bins of IR luminosity (with bin width of 0.5 dex).For the same IR luminosity (in the range log L IR /L = 12.0−12.5),their average M out and Ṁout are ∼0.6−0.8 dex higher than our measurements (factor of ×4−6).
Given the different methodology used to measure these parameters, the difference is not surprising.As discussed in Sect.4.5.2, this comparison highlights the large effect that the choice of methodology can have on the measured outflow A45, page 13 of 49  We extract the spectra from the central 9.4 × 9.4 spaxel and apply the point source aperture correction.We fit the OH spectra to derive the velocity of the outflow.In particular, we want to compare the velocity ranges where the OH outflow is detected with the ones of the CO outflow.Outflow parameters were derived by Veilleux et al. (2013) for 14 of the ULIRGs in our sample, but we repeat the analysis in order to obtain consistent parameters for the full sample.We set the zero velocity of the OH spectra based on the redshift of CO (see Table 1), so that we can directly compare the outflow velocities of the two tracers.
We fit the line profile following the method used by Veilleux et al. (2013), and we check that our derived parameters are consistent with their results.We first perform a linear fit of the continuum around the OH line and normalise the spectra to the continuum level.The OH doublet can appear in absorption (11/22), emission (4/22) or as P-Cygni profile (7/22).The P-Cygni profile is considered a clear indication of the presence of an outflow (e.g., Prochaska et al. 2011;Fischer et al. 2010;Sturm et al. 2011;González-Alfonso et al. 2012, 2013;Veilleux et al. 2013).In our sample, emission and P-Cygni profiles are more common in AGN (73%) than in SB-dominated nuclei, which is consistent with previous findings (Veilleux et al. 2013;Stone et al. 2016;Runco et al. 2020).
For the absorption profiles, we fit a model with two Gaussians (one for the systemic and one for the outflow) for each line of the OH doublet.The separation between the two lines of the Λ-doublet is fixed to 0.208 µm (in rest-frame) and the amplitude and width of the two lines were tied to be the same in each component.We convolve our model with the Herschel/PACS instrumental resolution (FWHM ∼ 270 km s −1 , Veilleux et al. 2013), to recover the intrinsic shape of the line.For profiles in emission, we fit a model with two Gaussians (one for the systemic and one for the outflow component) in emission for each OH line.For the P-Cygni profiles, we consider one component in absorption and one in emission for each line.Adding more components is not possible due to parameter degeneracy (Veilleux et al. 2013).We add as an additional constraint that the model absorption can not be deeper than the observed absorption (similar to Veilleux et al. 2013), to avoid a fitting result with unrealistically large amplitude of the emission and absorption components that cancel each other out.
We use the best-fit model to derive the characteristic velocity of the emission and absorption profiles, separately.For the absorption components, v 50 , v 84 and v 98 are the velocities corresponding to the 50th, 84th and 98th percentile of the absorption line profile, i.e. the velocities above which 50%, 84%, and 98% of the absorption takes place.Similarly, for the emission components, v 50 , v 84 and v 98 are the velocities corresponding to the 50th, 84th and 98th percentile of the emission line profile.Veilleux et al. (2013) consider v 84 to be a more robust estimate of the outflow velocity compared to the 'maximum outflow velocity'.We use v 98 as an estimate of the maximum outflow velocity, keeping in mind that it may be more susceptible to noise variation than v 84 .
We apply a Monte Carlo (MC) approach to estimate the uncertainty on the derived parameters.We estimate the noise level on a region of the continuum away from the line (excluding also the region of the CH + (3-2) and 18 OH 120 µm lines), then A45, page 14 of 49 we add random Gaussian noise, proportional to the noise level, to the best-fit model and we run the line fitting on this artificial spectrum.We repeat this procedure 50 times and we use the 16th and 84th percentiles of the distribution of best-fit parameters to estimate the 1σ uncertainties on the derived parameters.The 50 MC realisations are shown as orange curves on Fig. 11.The velocities derived from the OH profiles are shown in Table 5.The properties of molecular outflows derived from OH will be discussed and compared with those derived from the CO tracer in Sect.5.4.

Mean outflow properties
Here we summarise the ranges of outflow properties of our sample and their average values, corrected for inclination as described in Sect.4.5.1.We measure projected outflow velocities of ∼260−490 km s −1 , with a mean outflow velocity (corrected for inclination) of 485 ± 16 km s −1 .Outflow radii are in the range 0.17-0.94kpc and the mean inclination-corrected outflow radius is 1.1 ± 0.1 kpc.The outflow masses are between 1−35 × 10 7 M , with an average of (10 ± 2) × 10 7 M .These outflow masses corresponds to 0.2-6.5% of the total molecular gas masses.The mass outflow rates are in the range 6−302 M yr −1 , with an average outflow rate of 78 ± 16 M yr −1 .The ranges and averages of the outflow parameters measured in this work are summarised in Table 6.(Feruglio et al. 2010;Cicone et al. 2014).We note that, as discussed in Sect.4.6, the method used to derive the mass outflow rates can introduce systematic differences (up to factor of ∼10) between different samples.
We also compare our measurements with the properties of the ionised outflows measured by Arribas et al. (2014) for a sample of nearby U/LIRGs.For ULIRGs, they find maximum outflow velocities in the range 100-1000 km s −1 , with a mean 393 ± 38 km s −1 and mass outflow rates in the range 1−100 M yr −1 (from Fig. 13 in the paper), with an average 38 M yr −1 .The average Ṁout of the ionised gas is about a factor of 2 smaller than the one of the molecular phase.For U/LIRGs, they find outflow masses in the range 0.14−28 × 10 7 M , with a average of 6.7 × 10 7 M , similar to the masses of the molecular phase (Table 6).
We find molecular outflow dynamical times (t dyn = R out /v out ) in the range 0.45−2.77Myr.The mean t dyn , based on the mean observed R out and v out , is 1.37 Myr, assuming a average inclination correction of unity (Cicone et al. 2015).This is similar to the outflow dynamical times 0.63-2.51Myr reported in Pereira-Santaella et al. (2018) 00091−0738 Notes.
The top panel of Fig. 12 shows the percentage of nuclei in the sample belonging to each category: AGN, SB, merger (M), interacting (I), and the mixed categories (merger AGN, merger SB, interacting AGN, and interacting SB).For the interacting systems, we consider only the nuclei with log L IR /L > 11.8 in this statistical analysis.As shown in Pereira-Santaella et al. (2021), in most of our interacting systems the IR luminosity is dominated by one nucleus, thus by applying this luminosity threshold, we discard the fainter nuclei (with log L IR /L = 10.3−11.7)which would not be classified as ULIRGs.In this way, we avoid that the outflow detection rate in interacting nuclei is artificially lower only because of the 'secondary' faint nuclei.
The second panel of Fig. 12 presents the percentage of outflow detections in each category.We detect an outflow, defined as high-velocity (|v| > 300 km s −1 ) CO(2-1) emission which deviates from the main rotation, in 20/26 of the nuclei with log L IR /L > 11.8 (77±7%2 ).The nuclei with outflow detections are equally divided between mergers and interacting systems.The percentage of detections in SB (93 ± 4%, 14/15) is higher than in AGN (55 ± 14%, 6/11).A possible explanation for this difference is related to the outflow inclination: if AGN outflows are in the plane of the disk (contrary to SB outflows that tend to be perpendicular to the disk), they are more difficult to detect with our method, but by modelling and subtracting the disk rotation, it may be possible to identify them.Additionally, outflows in the plane of the disk may be braked and prevented from reaching high velocities, while outflows perpendicular to the disk can escape more freely, reaching higher velocities and being more easily detectable.Another possibility is that AGN outflows could be more collimated, and therefore more difficult to detect for some orientations (i.e. on the plane of the sky).A third possible explanation could be the different amount of molecular gas in the nucleus.AGN in our sample have on average lower nuclear molecular gas masses (log M(H 2 )/M = 9.2 ± 0.1), than SB nuclei (log M(H 2 )/M = 9.4 ± 0.1).The lower amount of material close to the nucleus may also be the reason that allows us to identify them as AGN, contrary to deeply buried nuclei.This is in agreement with an evolutionary scenario in which in a first merger phase the nuclei are more obscured and produced outflows, which expel gas and dust from the nuclear region; in a second phase, after some of the material has been removed, the nuclei are visible as AGN (e.g., Hopkins et al. 2008).Similarly, Stone et al. (2016) find a lower OH outflow detection fraction in X-ray selected AGN (24%) than in ULIRGs (∼70%) and suggest that outflow detection in ULIRGs may be easier due to their higher gas fraction.
With the spectro-astrometry analysis (see Sect. 4.2), we determine the outflow direction in 16/20 (80%) of the nuclei with outflow detection, while in the remaining four nuclei the outflow direction is not clear.This could be due to the fact that the outflow is unresolved, or to the fact that the outflow is pointing towards our line of sight, so that the blue and red-shifted sides of the outflow overlap.For the nuclei with a well determined outflow direction, we measure the angle between the outflow and the major axis of rotation (see Sect. 4.2).We qualitatively compare the direction of the high-velocity gas with the integrated maps of the blue and red-shifted high-velocity channels (see lower left panel in Figs. 3, 4 and E.1) to check that the direction is consistent with the position of the gas in the high-velocity channels.The bottom panel of Fig. 12 shows, for each category, the fraction of outflows with 1) direction perpendicular to the rotating disk (angle 90 ± 20 • ), 2) direction non-perpendicular, or 3) unclear direction.The percentages of outflow with determined direction in mergers (70 ± 12%, 7/10) is smaller than in interacting nuclei (90 ± 6%, 9/10).We could determine the outflow direction in 86±7% (12/14) of SB nuclei with outflow detection, but only in 67 ± 15% (4/6) of AGN.If in AGN the outflow is oriented within the plane of the disk, it is possible that it can not travel very far, thus it appear very compact and unresolved.This could explain why we could not determine the outflow direction in many AGN.
We find that in six nuclei the outflow direction is nearly perpendicular (angle 90 ± 20 • ) to the kinematic major axis.These nuclei are 12112+0305 NE, 12112+0305 SW, 14348−1447 NE (already presented in Pereira-Santaella et al. 2018), 10190+1322 E, 19542+1110 and 20100−4156 SE.All these nuclei are SB dominated.This supports the idea that outflow powered by SB tend to be perpendicular to the disk, while AGN outflow can have any orientation (e.g., Pjanka et al. 2017).However, there are also many SB nuclei (50±13%) for which the outflow direction is not perpendicular.This could be due to a hidden (deeply buried) AGN that is powering the outflow, even though it is not detected in the MIR or optical (see Pereira-Santaella et al. 2021).An alternative explanation is that in these SB the molecular gas is still strongly disturbed and has not yet settled into a galactic plane, and consequently also the path of least resistance is not well defined.We note that we measure the angle between outflow and major kinematic axis projected in the plane of the sky.Thus, there is the possibility that we measured a projected angle of 90 • for an outflow that is not perpendicular to the plane if we observe it from a particular orientation.However, for SB it is more likely that the outflow escape from the path of least resistance (perpendicular to the disk) than from any other direction, where the outflow will encounter more material.The majority of the nuclei with perpendicular outflow (5/6) are interacting systems.Although the number statistics are small, this could be due to the poorly defined disk rotation axis in some of the more advanced mergers.
Our sample was selected to have a similar number of objects in each category (AGN/SB, interacting/mergers), thus it does not have the L IR distribution or AGN fraction of the general population of ULIRGs.We estimate the outflow detection fraction that would be measured in a sample of ULIRGs with the same AGN fraction as the local ULIRGs population.Veilleux et al. (2009) study the AGN contribution in a representative sample of 74 ULIRGs from the IRAS 1 Jy Survey (Kim & Sanders 1998) and find that 24% of their sample has an AGN contribution >50% in the MIR.We use this AGN fraction to correct our outflow detection statistics and we find that the expected outflow detection fraction in local ULIRGs is 84% (compared to 77% measured in our sample).

Outflow quantities
In this Section, we investigate whether there are differences in the outflow properties depending on these categories.We calculate the mean outflow parameters (v out , R out , M out , Ṁout ) for each category and we compare it with the mean of the total sample.We do not include upper limits in this comparison.Even though there are small differences, the mean quantities in all categories are consistent (within 2σ) with the sample mean (see Fig. 13).bottom).The solid lines show the CDFs for the outflow detections only, while the dashed lines show the CDFs including upper limits (only for M out , Ṁout and η) .The shaded area mark the area between the CDFs with and without upper limits.The square symbols show the average values for the outflow detections in the two categories, the points with arrows show the upper limits (with arbitrary values on the y-axis).According to a survival analysis Two Sample test, the differences between AGN and SB and between interacting and mergers, are not statistically significant.Additionally, we compare the distribution of outflow properties for AGN vs. SB and mergers vs. interacting systems, taking into account also the upper limits on the non-detections.Figure 13 shows the cumulative distribution functions (CDF) of the outflow properties (v out , R out , log M out , log Ṁout , and log η = log( Ṁout /SFR)) for the AGN and SB categories (upper) and mergers and interacting (bottom) categories.To test whether two samples have different distributions, we perform a Two Sample test using the survival analysis package ASURV (Feigelson & Nelson 1985) which allows us to take into account upper limits.We find that the distributions of the outflow properties of AGN and SB are not significantly different according to the Gehan's, Logrank and Peto-Prentice's Two Sample Tests (p-value = 0.2-0.7).Similarly, we do not find significant differences between mergers and interacting systems (p-value = 0.2-0.9).
We also consider the mean outflow rate normalised by the total L IR (see Fig. 14), in order to remove the effect of the correlation between Ṁout and L IR (see Sect. 5.3).Also in this case there are no significant differences in the average Ṁout /L IR of the different categories.The average for the total sample is Ṁout /L IR = (3.8± 0.6) × 10 −11 M yr −1 L −1 .In the right panel, we also plot the average outflow rates normalised by AGN luminosity Ṁout /L IR,AGN for the AGN categories.The average Ṁout /L IR,AGN for AGN are ∼3−6 times higher than the Ṁout /L IR for starbursts.

Outflow properties and energetics
In this section, we investigate the relation between the outflow properties and the SFR, as well as with the AGN luminosity and AGN fraction.

Outflow and AGN properties
In this section, we look at trends between the outflow properties and the AGN luminosity and AGN fraction.We calculate the AGN luminosity (L  of η at large α AGN .On the other hand, our sample was selected without a prior knowledge of the presence of outflows.The difference could be also due to the method used to estimate the AGN fraction.We use α AGN estimated from the MIR 30 µm to 15 µm flux ratio.Cicone et al. (2014) estimate the AGN fraction using different methods: from the 5-8 µm spectral range using the method from Nardini et al. ( 2010), from a combination of MIR and FIR diagnostics as described in Veilleux et al. (2009), by taking the fraction of L AGN (estimated from [OIII] or from the X-ray luminosity) and L bol .If we were to use the method from Nardini et al. (2010), the α AGN for our sample would be lower (<0.7 for all sources).Cicone et al. ( 2014) also find a correlation between the molecular mass outflow rate and the AGN luminosity.A similar correlation was also reported in Fiore et al. (2017), who expanded the sample used by Cicone et al. (2014).Figure 16 shows Ṁout as a function of L AGN for our sample, together with the sample from Cicone et al. (2014), Lutz et al. (2020), the nearby (non-ULIRGs) AGN from Stuber et al. (2021), and the nearby (z < 0.13) type 2 AGN from Ramos Almeida et al. (2022).Ramos Almeida et al. (2022) find that their sample has Ṁout more than one order of magnitude below the relation derived by Cicone et al. (2014; dashed line on the Fig. 16).Our sample occupies the parameter space between the Ramos Almeida et al. ( 2022) objects and the Cicone et al. ( 2014) relation.Ramos Almeida et al. (2022) suggest that this scaling relation represents the upper boundary of the Ṁout versus AGN luminosity relation.Our sample confirms that for the same AGN luminosity, ULIRGs can have a wide range of Ṁout .We note that most of the objects in Cicone et al. (2014) were selected based on the presence of a molecular outflow (detected in CO or OH), while for our sample we do not apply this prior criterion.This may explain the fact that the sources in their sample have the highest Ṁout values for a given AGN luminosity.

Outflow energetics
In Fig. 17    To increase the dynamic range of SFR, we include in this figure also the 20 galaxies with candidate outflows from the Physics at High Angular resolution in Nearby GalaxieS (PHANGS 3 ) ALMA survey (Leroy et al. 2021;Stuber et al. 2021).The PHANGS-ALMA sample consists of 90 nearby (d < 24 Mpc) star-forming galaxies (log sSFR/yr −1 > −11) with stellar masses in the range 9.3 ≤ log M /M ≤ 11.1.Stuber et al. (2021) use ALMA CO(2-1) observations with a resolution of ∼100 pc to look for molecular outflows in their sample and they find 20 (22%) outflow candidates, of which 16 are classified as 'secure' candidates.Half of the candidates are classified as AGN hosts based on Véron-Cetty & Véron (2010).
The outflow parameters of this sample have been measured using a method comparable to ours.They integrated the total flux within the defined outflow velocity range to derive the outflow mass and measure the flux-weighted outflow velocity.They also measured flux-weighted outflow radii.Because of the limited S/N, we did not attempt to derive flux-weighted R out by measuring the outflow radii in each channel, but the R out derived from the 2D Gaussian fit should be representative of the radius where most of the flux is located.To convert the CO(2-1) outflow flux to outflow mass, Stuber et al. (2021) use the conversion factor α CO = 0.8 M /(K km s −1 pc 2 ) (Bolatto et al. 2013b), similar to our assumption of α CO = 0.78 M /(K km s −1 pc 2 ), and r 21 = 0.65 (Bolatto et al. 2013a;Leroy et al. 2013Leroy et al. , 2021;;den Brok et al. 2021) while we use the value of r 21 = 0.91 measured in ULIRGs (Papadopoulos et al. 2012).The SFR for this .These outflow parameters are lower than the ones we measure for PUMA, as it is expected for galaxies with lower SFR compared to ULIRGs.The left panel of Fig. 17 shows the relation between mass outflow rate and SFR.The dashed lines indicate constant massloading factors (η = Ṁout /SFR).The PUMA objects have η < 0.04−1.The PHANGS-ALMA galaxies have slightly higher η = 0.15−4.The dispersion in Ṁout for a given SFR is quite high (∼1 dex).Mass-loading factors in local starburst galaxies are typically lower than 2−3 (e.g., Bolatto et al. 2013b;Cicone et al. 2014;Salak et al. 2016), thus, the outflows of the PUMA sample are consistent with being powered by starburst.We note that the nuclei classified as AGN through MIR or optical diagnostics (orange symbols) do not show higher mass-loading factors compared with the rest of the sample.For the PUMA sample, the vertical errorbars are probably more extended towards higher values, given that with our method it is possible that we are underestimating Ṁout in some of the objects (see discussion in Sect.4.5.2).
In Fig. 17 we also show the outflow momentum rate ( Ṗout ) as a function of SFR.We follow Pereira-Santaella et al. ( 2018) to calculate the total momentum injected by supernova explosions.We assume that the momentum per supernova is 1.3×10 5 M km s −1 ×(n 0 /100 cm −3 ) −0.17 (Kim & Ostriker 2015), where n 0 = 100 cm −3 is the electron density; and that the supernova rate is 0.012×SFR(M yr −1 ) for the adopted IMF (Leitherer et al. 1999).Since all the PUMA galaxies have Ṗout smaller than the total momentum injected by supernovae (SNe, see dashed line in the middle panel of Fig. 17), their outflow momentum rate could be explained entirely by SNe.
The outflow kinetic luminosity is shown in the right panel of Fig. 17.The dashed lines show fractions of the energy produced by SNe (1%, 10%, 100%).The kinetic energy injected by supernova explosions is calculated as L SNe [erg s −1 ] = 9 × 10 41 × SFR [M yr −1 ] ( Leitherer et al. 1999), adapted for a Kroupa (2001) IMF.For the PUMA ULIRGs, the outflow kinetic luminosity is 0.2−5.9% of the energy produced by SNe, for the PHANGS-ALMA sample it is slightly lower (0.1−2.4%).The linear fit to the two samples gives: The slope is larger than one, meaning that, assuming that these outflows are driven by SNe, the coupling efficiency between ISM and SNe increases slightly with SFR.This slope is shallower than the value of 2 ± 0.2 found by Pereira-Santaella et al. ( 2018) using a smaller sample (15 objects).The difference could be partly explained by the different method used to derive the outflow parameters in Pereira-Santaella et al. ( 2018), which obtained higher mass outflow rates for ULIRGs (12-400 M yr −1 ).

Outflow launching mechanism
In this section we investigate the outflow launching mechanism.
In particular, we investigate whether the detected molecular outflows are momentum-driven or energy-driven.For momentumdriven outflows, the mass-loading factor (η = Ṁout /SFR) is proportional to v −1 out , while for energy driven outflows η ∝ v −2 out (e.g., Murray et al. 2005).It is possible to distinguish between the two scenarios by looking at the slope of the relation between log η and log v out : where α and b are the slope and intercept of the linear relation.In Fig. 18 we show log η as a function of log v out for the PUMA and PHANGS-ALMA samples.We fit a linear relation between these two quantities using the Monte Carlo Markov-Chain (MCMC) implementation PyStan4 .The lightblue shaded region shows the 1σ uncertainty on the fit, obtained by sampling the posterior distribution.We assume a systematic error  of 0.1 dex on v out and 0.3 dex on η.We also model the intrinsic scatter of the relation and we assume that the noise is normal distributed.We derive the best fit values and the corresponding 1σ uncertainties from the median and the 16th-84th percentiles, respectively, of the marginalised posterior distributions of the parameters (more details on the fitting procedure can be found in Appendix B).
The best linear bisector fit gives a slope α = −1.43 ± 0.28 (σ intr = 0.35).This slope is between the momentum-driven (α = −1) and energy-driven (α = −2) value.From Fig. 18, we can see that both slopes are included within the 1σ fit uncertainty (lightblue shaded area), and therefore we cannot distinguish between the two options on the basis of this figure.
In Fig. 19, we show the SFR as a function of v out .There is a strong positive correlation between these two quantities (r = 0.84, p-value < 0.01).The best linear fit gives: The slope of the relation is in agreement with the relation found for the ionised gas by Arribas et al. (2014, slope 0.24 ± 0.05) and for neutral gas by Rupke et al. (2005b, slope 0.21 ± 0.04).It is also in agreement with the theoretical prediction from Heckman et al. (2000) of the relation between velocity and starburst luminosity v max ∝ L 0.25 bol .This relation is derived under the assumption that starbursts have a maximum characteristic surface brightness (Lehnert & Heckman 1996;Meurer et al. 1997), and therefore their bolometric luminosity is proportional to their radius squared.We note that these relations consider v max , not the mean outflow velocity.We cannot use v max in our analysis, because it is not available for the PHANGS sample and we consider it less robust than v out .Nonetheless, for the PUMA sample we find a good correlation between v max and v out (r = 0.58, p-value < 0.01).Using this relation, we can substitute log SFR with 4 • log v out in the expression of the mass-loading factor: Substituting this in equation 11 gives: Thus, by looking at the relation between M out /R out and v out , we can derive the value of α.For the fit, we assume a systematic error of 0.1 dex on log v out and 0.3 dex on log(M out /R out ).
Figure 20 shows log(M out /R out ) as a function of log v out .
The best linear bisector fit relation has a slope of 2.61 ± 0.25, which implies α = −0.39.The α = −1 scenario (momentumdriven) is within the 1σ uncertainty (see shaded area in Fig. 20), while the α = −2 scenario does not agree with the fit.Therefore, our analysis favours the momentum-driven scenario as the primary launching mechanism for molecular outflows in ULIRGs.Previous works investigating the molecular outflows in ULIRGs also reach this conclusion.Cicone et al. (2014) find α = −1.0 ± 0.5 for a sample of five pure starburst galaxies, in agreement with the momentum-driven scenario.Pereira-Santaella et al. ( 2018) find α = −0.3± 0.2 combining a sample U/LIRGs and lower luminosity starbursts.We note that for this analysis we adopt a different r 21 value for the Stuber et al. (2021) sample (r 21 = 0.6) and for the PUMA sample (r 21 = 0.91).If we were to use the same value for both samples, the fit would agree even more with the momentum-driven scenario.

Escape fractions
We estimate the average fraction of the outflowing gas that can potentially escape from the gravitation potential of the galaxies.To estimate the escape fraction, we compare the outflow velocity with the escape velocity derived from a gravitational model of the host galaxy, which is assumed to be a truncated isothermal sphere (e.g., Veilleux et al. 2005).We use Eq. ( 7) in Arribas et al. (2014) to estimate the average escape velocity for our sample at the average outflow radius r = 1 kpc.Since estimates of the dynamical masses have been obtained by Perna et al. (2022) only for a subset of our targets (eight objects), we use the mean dynamical mass of the PUMA sample for our calculation M dyn = 3.9 × 10 10 M .The dynamical masses have been obtained from the modelling of the rotation of the ionised gas (see Perna et al. 2022, for more details).As truncating radius we assume two times the average effective radius R e = 2 kpc (Perna et al. 2022).We assume that the outflow velocity only has a radial component.With these assumptions, we estimate an average escape velocity of v esc = 486 ± 40 km s −1 .We apply an inclination correction to determine the average 'observed' v esc,obs = v esc /1.27 = 382 km s −1 (see Sect. 4.5).
Integrating the CO(2-1) emission at velocities higher than v esc,obs , we find that between 4−100% of the high-velocity gas will escape to the circumgalactic medium, with mean escape fraction f esc = 45 ± 6%.This calculation depends on the assumption on the truncating radius.If we were to assume a larger truncating radius, the escape velocity would increase, and the escape fraction would decrease.For example, assuming a truncating radius larger by a factor of two (r max = 4 × R e = 8 kpc), v esc would increase by 13%.The escape outflow rates are in the range 1-173 M yr −1 , with an average Ṁesc = 40 ± 10 M yr −1 .We find that interacting systems tend to have lower f esc than mergers (mean f esc = 33 ± 6% vs. f esc = 60 ± 10%).
We also compare the escape molecular gas mass (M(H 2 ) esc ) with the total M(H 2 ) of the galaxy (i.e.systemic and outflow).
The escape M(H 2 ) esc is <5% of the total M(H 2 ), with a mean of 1%.The mean depletion time based on the escape outflow rate (t dep = M(H 2 )/ Ṁesc ) is 158 Myr (range 23−3715 Myr).

OH vs. CO outflow properties
In this section we compare the molecular outflow properties derived from the CO(2-1) emission line with the ones derived from the OH 119 µm doublet, which is an alternative tracer used to identify molecular outflows.
We compare the OH and CO profile of each PUMA target in Figs. 3, 4, and E.1.The OH absorption is produced in front of the continuum emission, thus, we can assume that it is located in a region of the size of the continuum, which in most cases is equal or smaller than the beam size (Pereira-Santaella et al. 2021).Therefore, to have a fair comparison, we extract the nuclear CO spectrum from a region equal to the beam size.We also smooth the CO spectrum to the same spectral resolution of the OH spectrum.Qualitatively, it is evident from this comparison that there are some cases where a clear blue-shifted outflow is detected in OH but not in CO (e.g., 00188−0856), and vice-versa (e.g., 09022−3615, 17208−0014).In Fig. 21 we show four example cases: i) broad OH absorption reaching high negative velocities (v < −1000 km s −1 ), but no CO blue-shifted emission at v < −300 km s −1 (00188−0856); ii) CO outflow, but no sign of outflow in OH (22491−1808); iii) outflow detected in both OH and CO, but OH outflow maximum velocity (v 98 = −1540 km s −1 ) much larger than CO outflow maximum velocity v max = −800 km s −1 (20100−4156); iv) OH outflow maximum velocity (v 98 = −360 km s −1 ) smaller than CO outflow maximum velocity v max = −600 km s −1 (15327+2340).
In Fig. 22, we compare the OH and CO outflow velocities.In the upper row, we compare the OH v 98 velocity of the absorption profile with the 'weighted' outflow velocity (v out , left) and maximum outflow velocity (v max , right) of the blue-shifted wing of the CO profile.We show the objects with OH absorption or P-Cygni profile, for a total of 20 nuclei.For 12112+0305 and 14348−1447, we show both nuclei, since a CO outflow is detected around each nucleus and they have comparable continuum luminosity.For the other interacting systems, we consider only the nucleus with the brightest continuum, since the second nucleus does not have a CO outflow detection.
We identify three groups: i) targets with comparable OH and CO outflow velocity, ii) targets with OH outflow velocity larger than the ones from CO, and iii) targets with CO outflow velocity larger than the OH outflow velocity.In 13 targets, both |v 98 (abs)(OH)| and |v max (CO)| are >300 km s −1 .The |v 98 (abs)(OH)| values tend to be higher (on average by 170 km s −1 ) than the |v max (CO)| values, but this is not surprising given that v 98 and v max are measured using different methods.
For 3/20 targets, |v 98 (abs)(OH)| < 300 km s −1 , thus, there is no evidence of an outflow in the OH absorption profile, while a CO outflow is detected.This difference may be due to a collimated outflow in CO.In this case, the associated OH will absorb only a tiny fraction of the continuum behind the outflow and thus it will be undetectable.In Fig. 23 we compare the velocity difference |v 98 (abs)(OH) − v max (CO)| with the CO outflow radius R out .We find a weak anti-correlation between the R out and the OH-CO velocity difference (r = −0.4,or r = −0.3 if we exclude the point with v max (CO) = 0 kms −1 ), meaning that for more extended CO outflow, the CO velocity tend to be larger than the OH velocity.This is consistent with the scenario explained before, where for an outflow extending to a larger distance from the nucleus, a lower fraction of the outflowing gas may overlap with the  background continuum, making the OH absorption weaker (Veilleux et al. 2020).Additionally, if the velocity of the outflow increases with distance from the nucleus, the part with the highest velocity is more likely to not overlap with the continuum and thus, will not be detected in the OH absorption.Indeed, a radial velocity gradient of ∼1 km s −1 pc −1 has been detected in the molecular outflows of two galaxies: in the nearby starburst galaxy NGC 253 by Walter et al. (2017) and in the LIRG ESO 320-G030 by Pereira-Santaella et al. (2020).
For one target (00188−0856), the OH spectrum shows a clear outflow signature with v 98 (abs)(OH) ∼ −1200 km s −1 , but there is no outflow signature in the CO spectrum.The case of 00188−0856 can potentially be explained with extreme environments: if the ionisation fraction of the molecular gas is high, the abundance of CO would decrease while the abundance of OH remains high (González-Alfonso et al. 2018).This could be mostly associated with the highest-velocity gas.Alternatively, this difference may be explained by the geometry of the outflow.
Another factor that can explain part of the differences in outflow velocities is the difference in sensitivities of the CO and OH observations.Thus, the observation of one tracer may be more sensitive to weak outflow signatures at high velocities than the other.
We note that the AGN nuclei tend to have faster OH outflows than SB-dominated nuclei.Indeed, all AGN nuclei have |v 98 (OH)| > 300 km s −1 .There is a positive correlation (r = 0.51) between the outflow velocity difference |v 98 (OH) − v(CO)| and the AGN fraction (see Fig. 23).Veilleux et al. (2013) find that AGN-dominated nuclei (α AGN > 0.5) tend to have more negative OH outflow velocities than SB-dominated nuclei.Their interpretation of this result is that once a significant fraction of the material has been pushed away from the nucleus by the outflow, the AGN is more likely to be identifiable.At the same time, the central high-velocity part of the outflow can be seen more easily.
To investigate the differences between the OH and CO outflow velocities, we consider the possibility that the presence of a compact obscured nucleus (CON) could influence the OH outflow velocities.Falstad et al. (2019) find that CONs with the most luminous HCN-vib line lack signatures of high velocity outflows in the OH 119 µm absorption lines, even though some of them have molecular outflows detected in the CO or HCN emission lines.We apply the method described in García-Bernete et al. (2022) to identify CONs in our sample, based on the ratios of the equivalent widths of different polycyclic aromatic hydrocarbons (PAHs) features from the literature (see Table 1).We find that CONs do not occupy a specific region of the v 98 (abs)(OH) vs. v max (CO) diagram (see black circles in Fig. 22).We find that 6/7 CONs have |v 98 (abs)(OH)| > 300 km s −1 , a sign of molecular outflows.In Fig. 22, we also compared the velocities derived from the OH emission (v 98 (em)) with the velocities of the red-shifted CO outflow (v out (r) and v max (r)).Most targets have v 98 (em) > 300 km s −1 .For these targets, v 98 (em) generally agree, within the uncertainties, with the CO outflow velocities.The uncertainties on OH velocities derived from the emission profiles tend to be larger than the ones derived from the absorption, because the emission features are weaker.
One possible caveat of this analysis is the fact that the OH spectrum includes the nuclear region with a beam size ∼20 arcsec5 .For the interacting systems, these include the two nuclei.To detect the OH absorption, it is necessary to have the continuum emission in the background.In the majority of interacting cases, the ALMA continuum emission is dominated by one nucleus (>80% of the total emission for all but two cases, 11095−0238 (65%) and 14348−1447 (60%), see Table 7 in Pereira-Santaella et al. 2021), thus we can associate the OH absorption to the nucleus with the brightest continuum emission.

Discussion
In our sample, we find a lower molecular outflow detection rate in AGN (55 ± 14%) than in SB (93 ± 4%), for the nuclei with log L IR > 11.8 L (see Sect. 5.2).This finding is in contrast with the result of Stuber et al. (2021), who find higher molecular outflow detection rate in AGN (53%) compared to non-AGN (17%) in main-sequence galaxies from the PHANGS sample.A possible explanation is that the situation in ULIRGs is different than in main-sequence galaxies: in ULIRGs, SB-driven outflows are stronger and may be easier to detect than AGN-driven outflows.In main-sequence galaxies on the other hand, outflows driven by star-formation are weaker, and the AGN-outflows are more important.
Other than the detection rate, we do not find significant differences in the outflow properties (v out , R out , M out , Ṁout ) in AGN and SB in the PUMA sample.Thus, molecular outflow in AGN-dominated nuclei seem to be as powerful as SBdriven outflows, contrary to what was found in previous works (Cicone et al. 2014;Fluetsch et al. 2019).Cicone et al. (2014) and Fiore et al. (2017)  do not see this correlation, in agreement with the findings of Ramos Almeida et al. (2022).The previous relations maybe trace the most efficient coupling between AGN and molecular gas, so that the objects with the most extreme outflows follow those relations.This implies that there are other factors that set the mass outflow rate other than the AGN luminosity, as for example the geometry of the outflow and the coupling between the outflow and the CO disk, the efficiency of the energy/momentum transfer between the AGN and the molecular gas, or the distribution of the gas around the AGN.Additionally, starbursts could contribute to powering the molecular outflows even in AGNdominated nuclei.
In our sample, we find outflow dynamical times 0.5-2.8Myr (see Sect. 5.1).These dynamical times are significantly shorter than the expected age of the star-formation burst in ULIRGs (∼60−100 Myr, Rodríguez Zaurín et al. 2010), the depletion times of star-formation (10−100 Myr for our sample) or the outflow depletion times (∼10−700 Myr, Cicone et al. 2015;Pereira-Santaella et al. 2018).These short dynamical times of the molecular outflows may be due to the survival time of the molecular gas in the hot outflow environment (e.g., Decataldo et al. 2017).An alternative explanation may be related to the geometry of the outflow.If the outflow has a bi-conical geometry, the area of the outflow increases with the squared of the radial distance from the nucleus (r 2 ).Thus, the column density of the outflow decreases with radius, making the emission fainter and more difficult to detect at large radii (Pereira-Santaella et al. 2018).
In our sample, we find that between 4-100% of the outflow gas could escape from the gravitational potential of the galaxies into the circumgalactic medium (see Sect. 5.3.4).However, this represents only 1-5% of the total molecular gas reservoir of the galaxies.Thus, it is unlikely that these outflows are able to impact star-formation by removing gas from the molecular gas reservoir.

Summary and conclusions
In this work, we present ALMA CO(2-1) observations of 25 nearby (z < 0.165) ULIRG systems (38 individual nuclei) with a resolution of ∼400 pc.We have used these observations to study their molecular outflow properties.The main results of this work are as follows: 1. We detect molecular outflows in 20/26 (77%) nuclei with log L IR /L > 11.8.The molecular outflows have an average outflow velocity 485 ± 16 km s −1 , outflow masses 1−35 × 10 7 M , mass outflow rates Ṁout = 6−300 M yr −1 , and mass-loading factors η = Ṁout /SFR = 0.1−1.The majority of the outflows (18/20) are spatially resolved with radii of 0.2−0.9kpc and have short dynamical times (t dyn = R out /v out ) in the range 0.5−2.8Myr.We estimate the average escape velocity for our sample v esc = 486 ± 40 km s −1 .We find that, on average, 45 ± 6% of the high-velocity gas will escape to the circumgalactic medium, which represents less than 5% of the total molecular gas mass of the systems (see Sect. 4.5.) 2. We find that the outflow detection rate is higher in SBs (93±4%, 14/15) than in AGN (55±14%, 6/11) (see Sect. 5.2, Fig. 12).A possible interpretation is that SB outflows tend to be perpendicular to the disk, and thus they are easier to detect, while AGN outflow can have any inclination and/or can be more collimated.Indeed, we find that in 43% (6/14) of the SB nuclei with outflow detection, the projected position angle of the outflow is along the kinematic minor axis, which suggests that these outflows are perpendicular to the disk.This fraction is higher in early interacting SBs (50%) than in advanced merger SBs (17%).Outflows powered by AGN do not have a preferred orientation with respect to the disk.We do not find any significant difference in the mean outflow properties depending on the nuclear engine (AGN versus SBs) or on the merger stage (advanced mergers versus interacting systems, see Fig. 14).3. We find that our sample does not follow a tight Ṁout versus AGN luminosity relation, as reported in previous works (see Fig. 16).For the same AGN luminosity, we find Ṁout spanning up to 1.2 dex.This suggests that other factors may contribute in determining Ṁout in AGN ULIRGs, such as the efficiency of the energy or momentum transfer between the AGN and the molecular gas, the outflow geometry, or the distribution of the gas around the AGN. 4. We investigate whether the molecular outflows are momentum-driven or energy-driven, using also the PHANGS-ALMA sample (Stuber et al. 2021) to extend the analysis to galaxies with lower SFR.Based on the slope (α = −1.43 ± 028) of the mass-loading factor versus outflow-velocity relation (Fig. 18), we cannot distinguish between the momentum-driven (α = −1) and energy-driven (α = −1) scenarios.However, using the relation between the SFR and the outflow velocity (log v out ∝ 0.25 • log SFR), we could derive the value of α by fitting the relation between M out /R out and v out .From this analysis, we derived a slope α = −0.39 ± 0.25, which is more consistent with the momentum-driven than with the energy-driven scenario (see Sect. 5.3.3,Fig. 20). 5. We compare the outflow velocities derived from CO (v out ) with the ones derived from OH outflow (v 98 ; see Sect.5.4, Fig. 22).We find that for 13 targets, both the CO and OH velocities are above 300 km s −1 , which is a sign of outflow, while in three targets the OH velocities are considerably lower than v out , and therefore there is no evidence of outflows.This could be explained if the outflow is highly collimated, and thus the outflowing OH covers a small region of the background continuum.In one case (00188−0856), the OH outflow velocity is very high, but no outflow is detected in CO.This situation could be explained by an extreme environment, where the high ionisation would decrease the CO A45, page 25 of 49 abundance more than the OH abundance, by the outflow geometry or by the different sensitivities of the OH and CO observations.

Appendix A: Simulations of biconical outflows
In this section we provide some additional information about the simulations used to estimate the biases in the outflow properties measured with our method (see Sec. 4.5.2).We simulate a biconcal outflow with a half opening angle of 40 • .The maximum velocity of the outflow is 750 km s −1 and the maximum outflow extension is R out = 2 kpc.The outflow flux decreases with increasing radius (as it is expected in a scenario with constant mass outflow rate).We consider four distributions of the outflow velocity as a function of radius: 1) constant velocity, 2) velocity linearly decreasing with radius, 3) linear velocity law with an increasing velocity going from zero to v max at a given turnover radius R to = 0.5 × R out , followed by a linearly decreasing trend with v(R out ) = 0 km s −1 , 4) a radial velocity profile with turnover radius as in the previous case, but accounting for an initial velocity v(r = 0) = 400 km s −1 .We consider four values for the inclination of the outflow with respect to our line of sight: 10 • , 40 • , 70 • , and 90 • (where an inclination i = 90 • corresponds to the plane of the sky).The systemic component has a FWHM of 360 km s −1 and a flux ∼ 10 times the outflow flux.We measure v out , M out and Ṁout from the total simulated profile (out-flow+systemic component) considering only |v| > 300 km s −1 , to mimic the method we are using with our data.Then, we measure the outflow quantities (v out ,M out and Ṁout ) from the simulated outflow emission profile (without the systemic component) considering the full velocity range.Figure A.1 shows the ratios of the outflow quantities measured with our method (|v| > 300 km s −1 ) and from the total outflow profile.The differences between the two methods increase with increasing outflow inclination.A45, page 28 of 49

Appendix B: Bisector fit methodology
Here we explain in detail the methodology used for the bisector fit used in Figures 17,18,19,and 20.Considering two generic quantities x andy, we fit a linear relation between them using the Monte Carlo Markov-Chain (MCMC) implementation PyStan6 .We assume that the noise is normal distributed.We fit for the parameter α and β that minimises the quantity: where f (x i , α, β) = α • x i + β, y err,i is the error on the quantity y i , σ intr is the intrinsic scatter, and the sum is over the data points.The analogous expression can be defined with respect to the quantity x i and the error x err,i : We run the fitting algorithm to find the best parameters that minimise y and x , and then we perform a bisector fit.We run the algorithm for 2000 steps with 4 chains as 'burn-in' phase.We monitor the convergence by looking at the effective sample size (N e f f ), which is defined as the number of iterations divided by the integrated autocorrelation time N e f f = N iter /τ int .We check that the number of effective samples is > 10, which indicate that the algorithm has converged (Gelman et al. 2004).
We derive the best fit values and the corresponding 1σ uncertainties from the median and the 16th-84th percentiles, respectively, of the marginalised posterior distributions of the parameters.The intrisic scatter is obtained by adding in quadrature σ intr,x and σ intr,y .The lightblue shaded regions in the Figures shows the 1 σ uncertainty on the fit, obtained by sampling the posterior distribution.

Appendix C: CO outflow detection in absorption
We detect a blue-shifted absorption in the CO(2-1) spectrum of 07251-0248 E at velocities v = [−320, −500] km s −1 (see Fig. C.1).The absorption is compact, centred at the position of the continuum peak and nearly coincident with the position of the red-shifted outflow emission at v = [300, 550] km s −1 .This absorption can be interpreted as an outflow located in front of the continuum source and moving towards us.Alternatively, it could be due to an extremely compact outflow.The redshifted part of the outflow, identified in emission, appears compact and moves away from us along a direction compatible with our line of sight.This is the only target for which we could identify an absorption.Dasyra & Combes (2012) detected a blue-shifted absorption at −950 km s −1 in the CO(3-2) spectrum of 13451+1232 (4C+12.50),obtained with the IRAM 30m telescope.No absorption was detected in the CO(1-0) data observed with the IRAM Plateau de Bure Interferometer (PdBI) by Dasyra et al. (2014).We inspect the ALMA CO(2-1) spectrum of the west nucleus of 13451+1232 (shown in Fig. C.1), which is the nucleus dominating the FIR continuum emission and therefore, the nucleus where an absorption could be produced.We do not observe any significant absorption in the spectral channels on the blue-side of the CO(2-1) line in this target.Although, we note that the CO(3-2) absorption feature would be near the edge of the spectral window of our CO(2-1) spectrum.7.8 3.9 0 -3.9 -7.8

∆RA [arcsec]
1 kpc CO(2-1) moment 1 0.9 0 -0.9 -1.9The grey ellipse illustrates the ALMA beam FWHM.The grey contours on the moment 1 maps are every 50 km s −1 (every 25 km s −1 if the maximum value < 100 km s −1 ), and every 25 km s −1 (every 15 km s −1 if the maximum value < 150 km s −1 ) on the moment 2 map.In black are the CO(2-1) moment 0 contours ([3, 6, 25, 50, 75]×σ).d) Emission of the high-velocity channels(if detected), integrated over the velocity ranges indicated on the CO(2-1) spectrum (shown in panel e).Blue-and red-shifted channels are shown with blue and red contours, respectively (dashed lines indicate negative contour levels).The lowest contour corresponds to the 3σ level.The next contour levels are (0.5, 0.7, 0.9) of the peak of the emission, if these are above the 3σ level.The dashed circle shows the size of the outflow (R out ). e) CO(2-1) continuum-subtracted spectrum extracted from a circle with radius equivalent to the outflow size (R out ).The lower panel shows an y-axis zoom-in to highlight the emission in the wings.The horizontal dashed line show the 1σ noise level.The vertical dashed lines indicated the 'flux-weighted' velocity of the blue and red-shifted outflow (v out ). f) OH 119 µm spectrum (upper, if available) compared with the nuclear CO spectrum (bottom), convolved to the resolution of the OH spectrum (FWHM∼ 270 km s −1 ).The total fit to the OH lines is shown with a magenta line, while the Gaussian components of the fit are shown in lightblue and brown.The orange lines show the 50 Monte Carlo iterations used to estimate the uncertainties on the fit (see Sec.

Fig. 1 .
Fig. 1.Infrared luminosity (8-1000 µm) and redshift of the individual nuclei in the PUMA sample.AGN are shown in orange, while starburst (SB) dominated nuclei are in lightblue.Circles indicate interacting systems and squares indicate advanced mergers.The nuclei with log L IR /L < 12 are the fainter nuclei in interacting systems, where the IR luminosity is dominated by the other nucleus.

Fig. 3 .
Fig.3.Example of spectro-astrometry and outflow maps for one target (20100−4156 SE) with outflow direction perpendicular to the kinematic major axis.Panel a: spectro-astrometry of the CO(2-1) emission line, i.e. centroid position of the CO(2-1) emission in the different velocity channels.The points are colour-coded by the channel velocity.The pink diamond indicates the peak ALMA millimetre continuum position.The dashed line is a linear fit to the low-velocity points (kinematic major axis) and the dotted line is a fit to the high-velocity points (indicating the outflow direction, if present).Panels b,c: moment 1 and moment 2 maps, where the green square indicate the field of view of panel a.The grey ellipse illustrates the ALMA beam FWHM.The grey contours on the moment 1 maps are every 50 km s −1 (every 25 km s −1 if the maximum value <100 km s −1 ), and every 25 km s −1 (every 15 km s −1 if the maximum value <150 km s −1 ) on the moment 2 map.In black are the CO(2-1) moment 0 contours([3, 6, 25, 50, 75] × σ).Panel d: emission of the high-velocity channels, integrated over the velocity ranges indicated on the CO(2-1) spectrum (shown in panel e).Blue-and red-shifted channels are shown with blue and red contours, respectively (dashed lines indicate negative contour levels).The lowest contour corresponds to the 3σ level.The next contour levels are (0.5, 0.7, 0.9) of the peak of the emission, if these are above the 3σ level.The dashed circle shows the size of the outflow (R out ). Panel e: CO(2-1) continuum-subtracted spectrum extracted from a circle with radius equivalent to the outflow size (R out ).The lower panel shows an y-axis zoom-in to highlight the emission in the wings.The horizontal dashed line shows the 1σ noise level.The vertical dashed lines indicated the 'flux-weighted' velocity of the blue and red-shifted outflow (v out ). Panel f : OH 119 µm spectrum (upper) compared with the nuclear CO(2-1) spectrum (bottom), convolved to the resolution of the OH spectrum (FWHM ∼ 270 km s −1 ).The total fit to the OH lines is shown with a magenta line, while the Gaussian components of the fit are shown in lightblue and brown.The orange lines show the 50 Monte Carlo iterations used to estimate the uncertainties on the fit (see Sect. 4.7).The vertical dotted, dashed and dot-dashed lines show the v 50 , v 84 , and v 98 percentile velocities, respectively.The blue-shaded area in the upper panel shows the wavelength range between v 84 and v 98 .Figures for the rest of the sample are in the appendix.

Fig. 6 .
Fig.6.Simulated outflow profiles for outflow inclinations with respect to our line of sight of 10 • (upper) and 90 • (bottom).The outflow velocity in this particular simulation increases up to a turnover radius and then decreases.Up to four different velocity fields are considered in Appendix A. The magenta curve shows the outflow component, the grey curve the systemic component and the black curve the total profile.The vertical dashed lines show the flux-weighted v out measured from the total profile at |v| ≥ 300 km s −1 (coloured in blue and red).The vertical solid lines show the flux-weighted v out measured from the outflow component over the full velocity range.For high outflow inclination (close to the plane of the sky), the contribution to the outflow flux by the gas at low projected velocities (|v| < 300 km s −1 ) increases.

Fig. 7 .
Fig. 7. Comparison of different methods for measuring the outflow velocity.Left: ratio between the maximum outflow velocity v max,FWHM (Eq.(8)) and the flux-weighted outflow velocity (v out ) in our sample.Right: ratio between v max,FW10% (Eq.(9)) and v out .The horizontal lines show the 1-to-1 ratio.

Fig. 8 .Fig. 9 .
Fig. 8.Comparison of outflow parameters measured in this work with values reported inLutz et al. (2020).From left to right: outflow velocity, outflow radius, outflow mass and mass outflow rate.The dashed line marks the ratio of 1.Even though there are differences in outflow velocities and radii, the mass outflow rates agree within a factor of 3.

Fig. 12 .
Fig. 12. Outflow detection statistics for different categories: AGN, starburst (SB), mergers (M, d nuclei < 1 kpc), interacting systems (I, d nuclei > 1 kpc), and mixed categories (mergers AGN, mergers SB, interacting AGN, interacting SB).Upper panel: fraction of individual galaxy nuclei in each category with respect to the total sample.The scale on the right axis shows the number of objects.Middle panel: outflow detections fractions.The percentages have been calculated with respect to the total number of nuclei in each category.Lower panel: outflow orientation statistics divided in three groups: outflow projected orientation perpendicular to the kinematic major axis of the disk (angle 90 ± 20 • ), outflow projected orientation non-perpendicular, or orientation could not be determined.The percentages have been calculated with respect to the number of outflow detection in each category.Error bars in middle and lower panel show the 90% binomial confidence interval.
Fig. 13.Cumulative distribution function (CDF) of the outflow properties of AGN and starbursts (SB, top) and mergers (M) and interacting (I,bottom).The solid lines show the CDFs for the outflow detections only, while the dashed lines show the CDFs including upper limits (only for M out , Ṁout and η) .The shaded area mark the area between the CDFs with and without upper limits.The square symbols show the average values for the outflow detections in the two categories, the points with arrows show the upper limits (with arbitrary values on the y-axis).According to a survival analysis Two Sample test, the differences between AGN and SB and between interacting and mergers, are not statistically significant.

Fig. 14 .
Fig.14.Left: mean values of the ratio of Ṁout and total L IR for different categories: AGN, starbursts (SB), mergers (M), interacting (I), and the mix categories mergers AGN, mergers SB, interacting AGN and interacting SB.The errorbars show the uncertainty on the mean.The grey diamond and horizontal band shows the mean of the total sample and the corresponding uncertainty.Right: for the AGN categories, mean ratios of Ṁout and L IR,AGN .
we show the outflow rate, momentum rate ( Ṗout = Ṁout • v out ), and kinetic luminosity (L out = 1 2 Ṁout • v 2 out ) as a function of SFR.The SFR has been estimated from the IR luminosity, using the Kennicutt & Evans (2012) relation (withKroupa & (2021)  sample (stars), the AGN luminosities have been derived from the 14-195 keV X-ray luminosity.Symbols for the PUMA sample are as in Fig.1.The dashed and dotted lines show the relations presented byCicone et al. (2014) and byFiore et al. (2017), respectively, scaled down by a factor of three to match our outflow geometry definition.Weidner (2003) initial mass function (IMF)) after removing the fraction of luminosity associated with the AGN (α AGN ).

Fig. 17 .
Fig. 17.Mass outflow rate (left), outflow momentum rate (middle), and outflow kinetic luminosity (right) as a function of SFR.Lightblue and orange symbols indicate SB-and AGN-dominated nuclei, respectively.Circles indicate interacting systems and squares indicate mergers from the PUMA sample (including nuclei with log L IR /L < 11.8).Stars show the PHANGS-ALMA targets from Stuber et al. (2021).A representative errorbar is shown on the lower right.Left panel: the black dotted-dashed lines show lines of constant mass-loading factors (η = Ṁout /SFR) of 0.1, 1 and 10.Middle panel: the dashed line indicates the total momentum injected by SNe as a function of SFR.The dotted line shows the L SFR /c ratio, where L SFR is the IR luminosity for a given SFR derived using the Kennicutt & Evans (2012) relation.Right panel: the lines show the 1%, 10% and 100% of the energy produced by SNe.The lightblue lines in the three panels are the best linear fit to the data, with the shaded area indicating the 1σ uncertainty.

Fig. 18 .
Fig.18.Mass-loading factor versus outflow velocity for our sample (circles and square symbols) and the PHANGS-ALMA sample (stars).The lightblue line is the best linear bisector fit to the data and the shaded areas indicates the 1σ uncertainty on the fit with and without including the intrinsic scatter term (lighter and darker colour, respectively).The blue dashed line shows the predictions for a energy-driven outflow (α = −2) and the red dotted-dashed line for a momentum-driven outflow (α = −1).

Fig. 20 .
Fig. 20.Outflow mass divided by outflow radius (M out /R out ) versus outflow velocity.The lightblue line is the best linear fit to the data (slope = 2.61 ± 0.25).The shaded area indicates the 1σ uncertainty with and without including the intrinsic scatter term (lighter and darker colour, respectively).The blue dashed line shows the predictions for a energy-driven outflow (α = −2) and the red dotted-dashed line for a momentum-driven outflow (α = −1).Symbols as in Fig.17.

Fig. 22 .
Fig. 22.Comparison of the outflow velocities of OH and CO.Upper row: v 98 derived from the OH absorption profile compared with the mean (v out (b), left) and maximum (v max (b), right) velocities of the CO blue-shifted wing.The background colours highlight the three regions described in Sect.5.4: (i) targets with v 98 similar to v max (b), (ii) targets with v 98 > v max (b), and (iii) targets v 84 < v max (b).Lower row: v 98 derived from the OH emission profile compared with the mean (v out (r), left) and maximum (v max (r), right) velocities of the CO red-shifted wing.The blue lines mark the outflow detection threshold of |v| > 300 km s −1 .Black circles mark targets classified as compact obscured nuclei (CONs).Colour-code as in Fig. 1.

Fig. 23 .
Fig. 23.Comparison of the difference between the outflow velocities of OH (v 98 ) and CO (v max ) with the CO outflow radius R out (left) and the AGN fraction α AGN (right).Black circles mark targets classified as compact obscured nuclei (CONs).Colour-code as in Fig. 1.
Fig. A.1.Ratios of the outflow properties (v out , M out and Ṁout ) measured with our method from the total profile (systemic+outflow) considering only the |v| > 300 km s −1 range and from only the outflow profile considering the full velocity range.The ratios are shown as a function of the outflow inclination with respect to our line of sight (where i = 90 • is in the plane of the sky).The points are colour-coded according to outflow velocity radial distribution used in the simulation.The black lines show the one-to-one ratio.
Fig. C.1.Left: CO(2-1) spectrum (with continuum) of 07251-0248 E , showing an absorption at velocities v = [−320, −500] km s −1 .Middle:Emission of the blue-and red-shifted high-velocity channels of 07251-0248 E, shown with blue and red contours, respectively.The lowest contour corresponds to the 3σ level.The next contour levels are (0.5, 0.7, 0.9) of the peak of the emission.Dashed lines indicate negative contour levels ([-3, -4, -5, -6, -7]×σ).The dashed black circle shows the size of the outflow (R out ) and the grey cross the central position of the outflow.The black and green crosses shows the position of the continuum peaks of the E and W nuclei, respectively.Right: CO(2-1) spectrum (with continuum) of 13451+1232 W. No clear absorption is detected in this spectrum at the velocities of the an absorption detected in the CO(3-2) spectrum byDasyra & Combes (2012) at v = [−700, −1200] km s −1 (see orange shaded region).
Fig. E.1.Spectro-astrometry and outflow maps for 00091-0738 S and 00091-0738 N .a) Spectro-astrometry of the CO(2-1) emission line, i.e. centroid position of the CO(2-1) emission in the different velocity channels.The points are colour-coded by the channel velocity.The dotted line is a linear fit to the low-velocity points (kinematic major axis) and the dashed line is a fit to the high-velocity points (indicating the outflow direction, if present).Panel b) and c) show the moment 1 and moment 2 maps, where the green square indicate the field of view of panel (a).The grey ellipse illustrates the ALMA beam FWHM.The grey contours on the moment 1 maps are every 50 km s −1 (every 25 km s −1 if the maximum value < 100 km s −1 ), and every 25 km s −1 (every 15 km s −1 if the maximum value < 150 km s −1 ) on the moment 2 map.In black are the CO(2-1) moment 0 contours([3, 6, 25, 50, 75]×σ).d) Emission of the high-velocity channels(if detected), integrated over the velocity ranges indicated on the CO(2-1) spectrum (shown in panel e).Blue-and red-shifted channels are shown with blue and red contours, respectively (dashed lines indicate negative contour levels).The lowest contour corresponds to the 3σ level.The next contour levels are (0.5, 0.7, 0.9) of the peak of the emission, if these are above the 3σ level.The dashed circle shows the size of the outflow (R out ). e) CO(2-1) continuum-subtracted spectrum extracted from a circle with radius equivalent to the outflow size (R out ).The lower panel shows an y-axis zoom-in to highlight the emission in the wings.The horizontal dashed line show the 1σ noise level.The vertical dashed lines indicated the 'flux-weighted' velocity of the blue and red-shifted outflow (v out ). f) OH 119 µm spectrum (upper, if available) compared with the nuclear CO spectrum (bottom), convolved to the resolution of the OH spectrum (FWHM∼ 270 km s −1 ).The total fit to the OH lines is shown with a magenta line, while the Gaussian components of the fit are shown in lightblue and brown.The orange lines show the 50 Monte Carlo iterations used to estimate the uncertainties on the fit (see Sec. 4.7).The vertical dotted, dashed and dot-dashed lines show the v 50 , v 84 , and v 98 percentile velocities, respectively.The blue-shaded area in the upper panel shows the wavelength range between v 84 and v 98 .

Table 1 .
Properties of the sample.

Table 4 .
Possible biases in the outflow properties measured with our method, due to different issues.
Comparison of the outflow mass and mass outflow rate as a function of IR luminosity for our sample and the sample from Fluetsch et al. (2019, purple triangles).The red symbols show the average values (considering only the detections) in bins of log L IR for the Fluetsch et al. (2019) sample (triangles) and the PUMA sample (circles).The sources in common between the two samples are highlighted with black contours.Due to the different methods, the M out and Ṁout values estimated by Fluetsch et al. (2019) are higher by a factor of ×5−8 compared with our measurements at the same infrared luminosity.
Spoon et al. (2013), andGonzález-Alfonso et al. (2017)illeux et al. 2013;González- Alfonso et al. 2014cted veStone et al. 2016)mulations of biconical outflows (e.g.,Bae & Woo 2016), this scenario is possible when the outflow is oriented in a direction close to the plane of the sky, or if the outflow has low velocities.In general, this method would tend to over-estimate the outflow mass and it suffers from the degeneracy of the fit with multiple Guassian components.With our method on the other hand, we could be missing the outflow contribution at the low projected velocities, thus, we may be underestimating the outflow mass.4.7.Analysis of the OH 119 µm spectraOther lines that are often used as tracer of molecular outflows are the OH (hydroxyl) FIR lines (e.g.,Fischer et al. 2010;Sturm et al. 2011;Spoon et al. 2013;Veilleux et al. 2013;González- Alfonso et al. 2014, 2017;Stone et al. 2016).We compare the outflow parameters derived from CO(2-1) with the ones derived from the OH 119 µm doublet, which is the strongest of the OH ground-state lines (e.g.,González-Alfonso et al. 2017).We collect archival Herschel/PACS spectra of the OH 119 µm doublet transition (hereafter OH) for 22/25 of our ULIRGs systems (no data available for 00091−0738, 10190+1322, and 16156+0146).The majority of the spectra (20/22) were published inVeilleux et al. (2013),Spoon et al. (2013), andGonzález-Alfonso et al. (2017); 11095-0238 and 13451+1232 are not presented in these works but were found in the Herschel archive.
. The outflow depletion times (t dep = M(H 2 ) tot / Ṁout ) in our sample are in the range 15−644 Myr, with a median of 75 Myr.These are a bit longer than the values reported for ULIRGs inPereira-Santaella et al. (2018,  15−80 Myr)andCicone et al. (2015, 1.2−50 Myr).For comparison, the star-formation depletion times (t dep = M(H 2 ) tot /SFR) for targets with detected outflows in our sample are in the range 9-77 Myr (median 27 Myr).Example of the fit of the OH 119 µm doublet with absorption (upper), emission (middle) and P-Cygni (bottom) profile.The velocity is relative to the blue component of the doublet (at 119.233 µm) at the systemic velocity inferred from the CO(2-1) line.The best-fit model is shown in magenta, while the orange models are the results of the 50 Monte Carlo realisations that illustrate the uncertainty of the fit.Solid lines show the model convolved with the instrumental resolution, while dashed lines indicate the intrinsic (deconvolved) model.Blue and red mark the two components of the fit.The vertical dotted, dashed and dotdashed lines show the v 50 , v 84 and v 98 percentile velocities, respectively, derived from the absorption (blue) and/or emission (red) profiles.Grey dashed lines at ∼2000 km s −1 show the position of the 18 OH 120 µm doublet and the CH + (3-2) absorption line.
residuals Fig. 11.5.2.Outflow characteristics for AGN/SB and interacting/mergers5.2.1.Outflow detection rate and directionIn this Section, we present the statistics of the number of detected molecular outflows in our sample and we investigate A45, page 15 of 49

Table 5 .
Velocities derived from the OH 119 µm profiles.

Table 6 .
Mean cold molecular outflow properties in ULIRGs.
( )Not corrected for inclination.
Cicone et al. (2014)tion of the total infrared luminosity L IR due to AGN, based on the MIR AGN fraction (α AGN , seePerna et al. 2021).For the interacting systems, α AGN is determined only for the IR brightest nucleus, thus, for the second faint nucleus we cannot estimate L IR,AGN and L IR,SF .Cicone et al. (2014)reported a trend for the massloading factor (η = Ṁout /SFR) to increase with AGN fraction (α AGN = L AGN /L bol ), while we do not find a correlation between Cicone et al. (2014)hawinski et al. 2015)tion from star-formation to the mass outflow rate for small AGN fractions or ii) to the short timescale of AGN variability (10−10 5 yr,Gilli et al. 2000;Schawinski et al. 2015), compared to the outflow timescales (10 6 yr).SinceCicone et al. (2014)included in the study objects with known molecular outflows, it is possible that they sample the upper end of the distribution A45, page 18 of 49 Mass outflow rate versus AGN luminosity.Diamonds symbols show the sample from Cicone et al. (2014), for which Ṁout have been scaled down by a factor of three to match our outflow geometry definition.The sample of Ramos Almeida et al. (2022) is shown with triangles.Horizontal lines for the Ramos Almeida et al. (2022) sample indicate the position that they would occupy if we were to use the L AGN derived from SED fitting from Jarvis et al. (2019), instead of the one derived from the [O iii] luminosity (see Ramos Almeida et al. 2022, for details).The Lutz et al. (2020) sample is shown with red pentagons.For the Stuber et al.
Outflow velocity versus star-formation rate for our sample and the PHANGS-ALMA sample.The lightblue line is the best linear fit to the data (slope: 0.25 ± 0.01); the 1σ uncertainties on the fit (with and without including the intrinsic scatter term) are indicated with shaded areas (lighter and darker colour, respectively).Symbols as in Fig.17.