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

PUMA IV. Spatially resolved cold molecular outflows in ULIRGs

I.  Lamperti 1E-mail<EMAIL_ADDRESS>M.  Pereira-Santaella1  M.  Perna1, 2  L.  Colina1 S.  Arribas1 S.  García-Burillo3 E.  González-Alfonso4

S.  Aalto 5 A.  Alonso-Herrero6 F.  Combes7 A.  Labiano8 J.  Piqueras-López1 D.  Rigopoulou9 P.  van der Werf10

I.  Lamperti et al.

Centro de Astrobiología (CAB),  CSIC-INTA,  Ctra. de Ajalvir Km. 4, 28850 Torrejón de Ardoz, Madrid, Spain
INAF - Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy
Observatorio Astronómico Nacional (OAN-IGN)-Observatorio de Madrid, Alfonso XII, 3, 28014, Madrid, Spain
Universidad de Alcalá,  Departamento de Física y Matemáticas,  Campus Universitario, 28871 Alcalá de Henares,  Madrid, Spain
Department of Space, Earth and Environment, Division of Astronomy and Plasma Physics, SE-412 96 Gothenburg, Sweden
Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, 28692, Villanueva de la Cañada, Madrid, Spain
LERMA,  Obs. de Paris,  PSL Univ.,  Collége de France, CNRS, Sorbonne Univ., Paris, France
Telespazio UK for ESA, ESAC, E-28692 Villanueva de la Cañada, Madrid, Spain

Department of Physics,  University of Oxford,  Keble Road, Oxford OX1 3RH, UK
Leiden Observatory,  Leiden University,  PO Box 9513, 2300, RA Leiden, The Netherlands

We study molecular outflows in a sample of 25 nearby ($z< 0.17$,  $d<750$ Mpc) ULIRG systems ( individual nuclei) as part of the  Physics of ULIRGs with  MUSE and ALMA (PUMA) survey,   using $\sim400$ pc (0.1-1.0" beam FWHM) resolution ALMA CO(2-1) observations.
We used a spectro-astrometry analysis to identify high-velocity ($> 300$ ) molecular gas disconnected from the galaxy rotation, which we attribute to outflows.
In 77% of the 26 nuclei with $\log L_{IR}/L_{\odot}>11.8$,  we identified molecular outflows with an average $v_{out}= 490$ , outflow masses $1-35 \times 10^7$  , mass outflow rates $=6-300$  yr$^{-1}$,  mass-loading factors $\eta = \dot{M}_{out}/SFR = 0.1-1$, and an average outflow mass escape fraction of $45\pm6\%$.
The majority of these outflows (18/20) are spatially resolved with radii of $0.2-0.9$ kpc and have short dynamical times ($t_{dyn}=R_{out}/v_{out}$)  in the range $0.5-2.8$ Myr.
The outflow detection rate is higher in nuclei dominated by starbursts (SBs,  $14/15=93\%$) than in active galactic nuclei (AGN,  $6/11=55\%$).
Outflows perpendicular to the kinematic major axis are mainly found in interacting SBs.
We also find that our sample does not follow the  versus  AGN luminosity relation reported in previous works.
In our analysis, we include a sample of nearby main-sequence galaxies  (SFR = $0.3-17$  yr$^{-1}$) with detected molecular outflows from the PHANGS-ALMA survey to increase the $L_{IR}$ dynamic range.
Using these two samples, we find a correlation between the outflow velocity and the star-formation rate (SFR),  as traced by $L_{IR}$ ($v_{out} \propto SFR^{0.25\pm0.01})$,  which is consistent with what was found for the atomic ionised and neutral phases.
Using this correlation,  and the relation between / and , we conclude that these  outflows are  likely  momentum-driven.
Finally, we compare the CO outflow velocities with the ones derived from the OH 119 doublet.
In 76% of the targets,  the outflow is detected in both CO and OH,  while in three  targets (18%)  the outflow is only detected in CO, and in one target  the outflow is detected in OH but not in CO.
The difference between the OH and CO outflow velocities could be due to the far-IR background source required by the OH  absorption which makes these observations more dependent on the specific outflow geometry.

§ INTRODUCTION

Local ultraluminous infrared galaxies (ULIRGs) are extreme objects with infrared (IR, 8-1000 $\mu$m) luminosities $ L_{IR}>10^{12} L_{\odot}$.
ULIRGs are mostly gas-rich major mergers  [Lonsdale et al., 2006] and represent an important stage in galaxy evolution.
The nuclei of ULIRGs host intense starbursts,  and active galactic nuclei (AGN) can also account for a significant fraction of their IR luminosity [Farrah et al., 2003, Nardini et al., 2010].

Nuclear outflows powered by starbursts and AGN  are thought to play a significant role in the evolution of galaxies.  They can influence star-formation by injecting energy into the interstellar medium (ISM),  heating  or expelling gas from the galaxy <cit.>.
In ULIRGs, outflows have been detected  in different gas phases: atomic neutral <cit.>,  atomic ionised <cit.>,
hot molecular <cit.>,
and cold molecular <cit.>.

In this work, we focus on the cold molecular phase, which is thought to account for most of the mass outflow rate <cit.>.
The most common tracers of the cold molecular phase are the low-J transitions of CO <cit.>, HCN    <cit.>,   and FIR OH lines <cit.>.

To  accurately measure the outflow properties, such as the  mass outflow rate, outflow energy and momentum rate,  it is necessary to know the distribution of the outflowing gas.  Spatially resolved studies of molecular outflows in ULIRGs  have only been  performed on a few sources <cit.>.
In most cases, the targeted  objects were selected based on the presence of an outflow,  which was detected through previous low-resolution observations. This could introduce a bias in the samples,  in which only ULIRGs with the most-extreme outflows are selected.

In this work, we present  high spatial resolution (0.1-1.0",  $\sim400$ pc) ALMA observations of the CO(2-1) line for an unbiased sample of 25 ULIRGs,  selected only based on their IR luminosity.

This paper is part of the Physics of ULIRGs with MUSE and ALMA (PUMA) project. The main goals of PUMA are the following: (i) to study the prevalence of outflows in different gas phases (ionised, neutral, and molecular) as a function of the galaxy properties, (ii) to determine the driving mechanisms of the outflows (e.g. distinguish between starburst- and AGN-powered outflows), and (iii) to  characterise the effects of outflow feedback on the host galaxies.
The PUMA survey combines VLT/MUSE optical integrated field spectroscopy and ALMA CO(2-1) and continuum observations to study the multi-phase (ionised,  neutral, and molecular) properties of outflows in ULIRGs. Perna et al., 2021 have presented the first results on the spatial distribution of the ionised gas and the resolved stellar kinematics derived from the MUSE data. A detailed analysis of the MUSE data for Arp 220 has been presented in Perna et al., 2020.  Perna et al., 2022 studied the properties and incidence of the ionised gas disks in the PUMA ULIRGs,  the associated velocity dispersion and its relation with the offset from the main sequence.
Pereira-Santaella et al., 2021  analysed the ALMA 220 GHz continuum and  provided evidence for the ubiquitous presence of deeply obscured AGN in ULIRGs, which could substantially contribute to their IR emission.

This paper is organised as follows.
In Section <ref>, we describe the sample selection criteria and the general properties of the PUMA targets.
Section <ref> describes the ALMA observations and the data reduction.
In Section <ref>,  we present the CO(2-1) moment maps (Sec. <ref>), the spectro-astrometry analysis  (Sec.  <ref>), and the method used to derive the properties of the molecular outflows  (Sec. <ref>- <ref>),  as well as the analysis of the OH 119  spectra (Sec.  <ref>).
The outflow properties,  detection rate,  energetics, and  launching mechanisms are discussed in Section <ref>-<ref>.
In Section <ref>, we compare the outflow properties derived from the CO(2-1) data with the ones derived from the OH 119 doublet.
The discussion is presented in Section <ref>.
In Section <ref>, we summarise the main results and our conclusions.
Throughout this work, we assume a cosmological model with $\Omega_{\lambda} = 0.7$, $\Omega_{\text{M}}= 0.3$ , and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$.

§ SAMPLE

The PUMA sample consists of 25 nearby ($z<0.165$,  $d<750$ Mpc) ULIRG systems ( individual nuclei) with IR luminosities ($8-1000~\mu$m) in the range $\log L_{IR}/L_{\odot} =11.9-12.7$.  The individual nuclei have luminosities from  $\log L_{IR}/L_{\odot} < 10.5$ to $\log L_{IR}/L_{\odot} =12.7$,  based on the relative contribution of the nuclei to the total ALMA continuum fluxes <cit.>.
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.3$ kpc).
The nuclei of the ULIRGs can be classified as AGN-dominated or starburst- (SB-) dominated based on the  AGN contribution in the MIR <cit.>: eight  ULIRG systems are dominated by AGN ($\alpha_{AGN}\geq0.5$) and 17  by a SB ($\alpha_{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 LINERs
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 ($\alpha_{AGN}\geq0.5$) or based on the optical classification (Seyfert 1 or Seyfert 2).  According to this combined classification,  11/25 systems  (14/ individual nuclei) are classified as AGN, while the others are classified as SBs.

An overview of  the sample properties is shown in Table <ref>.
Figure <ref> shows the IR luminosities and redshift of the sample.

Infrared luminosity (8-1000 $\mu$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_{\odot}<12$ are the fainter nuclei in interacting systems, where the IR luminosity is dominated by the other nucleus.

Properties of the sample.

IRAS name     other name     nucleus     RA      Dec      log $L_{IR}$    morph.       $d_{nuclei}$     opt. class.       $\alpha_{AGN}$     AGN/SB     CON

[deg]     [deg]       [log ]            [kpc]

(1)     (2)      (3)     (4)    (5)     (6)     (7)     (8)     (9)     (10)     (11)     (12)

00091-0738         S     2.930301     -7.368707     12.33     I     2.3     HII     0.46     SB     Y

N     2.930424     -7.368383     < 10.82

00188-0856          -     5.360472     -8.657220     12.42     M     < 0.3     Sy 2     0.50     AGN     N

00509+1225     I Zw 1     -     13.395560     12.693316     11.87     M     < 0.2     Sy 1     0.92     AGN     N

01572+0009     Mrk 1014     -     29.959380     0.394688     12.65     M     < 0.4     Sy 1     0.65     AGN     N

05189-2524          -     80.255834     -25.362582     12.20     M     < 0.1     Sy 2     0.72     AGN     N

07251-0248         E     111.906721     -2.915070     11.40     I     1.8     HII     0.30     SB     Y

W     111.906385     -2.915106     12.41

09022-3615          -     136.052940     -36.450537     12.33     M     < 0.1     HII     0.54     AGN     N

10190+1322         E     155.428143     13.115447     11.12     I     7.2     HII     0.17     SB     N

W     155.427056     13.114952     11.98

11095-0238         NE     168.014095     -2.906371     11.84     I     1.1     LINER     0.44     SB     Y

SW     168.013994     -2.906468     12.16

12071-0444         N     182.438000     -5.020377     12.41     I     2.3     Sy 2     0.75     AGN     N

S     182.437998     -5.020812     < 11.66

12112+0305         NE     183.441903     2.811542     11.27     I     4.2     LINER     0.17     SB     N

SW     183.441417     2.810866     12.28

13120-5453     WKK 2031     -     198.776347     -55.156339     12.27     M     < 0.1     Sy 2     0.33     AGN     N

13451+1232     4C+12.50    W     206.889004     12.290064     12.31     I     4.3     Sy 2     0.82     AGN     N

E     206.889583     12.289944     < 10.31

14348-1447         SW     219.409503     -15.006729     12.21     I     5.5     LINER     0.17     AGN$^{\star}$     Y

NE     219.409988     -15.005908     11.98                       SB$^{\star}$

14378-3651          -     220.245888     -37.075537     12.15     M     < 0.1     Sy 2     0.21     AGN     N

15327+2340     Arp220    E     233.738722     23.503135     11.65     M     0.4     LINER     0.06     SB     Y

W     233.738431     23.503177     12.04

16090-0139          -     242.918411     -1.785098     12.62     M     < 0.7     HII     0.41     SB     Y

16156+0146         NW     244.539016     1.656043     12.14     I     8.3     Sy 2     0.70     AGN     Y

SE     244.539754     1.655458     < 11.54

17208-0014          -     260.841481     -0.283582     12.43     M     < 0.1     LINER     0.00     SB     Y

19297-0406         N     293.092955     -4.000286     < 10.45     I     1.1     HII     0.23     SB     Y

S     293.092904     -4.000500     12.45

19542+1110          -     299.149103     11.318064     12.09     M     < 0.1     LINER     0.26     SB     N

20087-0308          -     302.849441     -2.997422     12.47     M     < 0.9     LINER     0.20     SB     N

20100-4156         SE     303.373149     -41.793113     12.65     I     6.5     LINER     0.26     SB     Y

NW     303.372813     -41.792383     < 10.96

20414-1651          -     311.075663     -16.671340     12.24     M     < 0.2     HII     0.00     SB     N

22491-1808         E     342.955620     -17.873368     12.22     I     2.7     HII     0.15     SB     Y

W     342.955158     -17.873239     < 10.53

(1) IRAS name.
(2) Alternative name.
(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  <cit.>.
(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 <cit.>.
(10) Fraction of AGN contribution to $L_{bol}$ derived from the 30to 15 flux ratio <cit.>.
(11) Combined nuclear activity classification used in this paper. Nuclei are classified as AGN either if $\alpha_{AGN}\geq 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.
$^{\star}$: 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] .
$^{\star \star}$: 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.

§ OBSERVATIONS AND DATA REDUCTION

In this work, we analyse ALMA 12-m array CO(2-1) 230.538 GHz observations of the 25 ULIRGs systems in the PUMA sample.
Most of these observations were carried out as part of our programmes 2015.1.00263.S, 2016.1.00170.S,  2018.1.00486.S, and 2018.1.00699.S  (PI: M.  Pereira-Santaella).  Additionally, we use archival data for  13120-5453 and 15327+2340 (Arp 220),  from programmes 2016.1.00777.S (PI: K.  Sliwa) and 2015.1.00113.S (PI: N.  Scoville),  respectively.
The observations were taken between 2015 and  2021.
We note that the analysis of the CO(2-1) observations of three of the PUMA  ULIRGs (12112+0305,  14348-1447,  and 22491-1808) has already  been presented in Pereira-Santaella et al., 2018,  but are included in this paper for completeness.

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 ($\sim400-500$ pc).
As a result, the synthesised beam FWHM  are in the range $0.13-1.05$".
The maximum recoverable scales are $\sim10$ times the beam FWHM.
Table <ref> 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.

Properties of the ALMA CO(2-1) observations.

IRAS name     nucleus     synthesised beam     sensitivity     int.  time     $z_{CO}$     $S_\text{CO} $     $\log L'_\text{CO}$     $\log M(H_2)$     FWHM     PA

(arcsec$\times$arcsec)      [$\mu$Jy beam$^{-1}$]     [min.]          [Jy ]       [K pc$^2$]     [$M_{\odot}$]      []      [deg.]

(1)     (2)      (3)     (4)    (5)     (6)     (7)     (8)     (9)     (10)     (11)

00091-0738     S     0.29$\times$0.24      552     100     0.117984     14.4$\pm$0.2     9.33     9.26     294$\pm$5     303

N                 0.118188     12.1$\pm$ 0.2     9.25     9.18     203$\pm$5     300

00188-0856          0.14$\times$0.12      265     69     0.128500     20.6$\pm$0.2     9.55     9.48     274$\pm$5     268

00509+1225          0.35$\times$0.31      332     48     0.061112     56.9$\pm$0.7     9.29     9.22     376$\pm$5     136

01572+0009          0.16$\times$0.13      279     69     0.163271     10.8$\pm$0.3     9.47     9.40     183$\pm$16     306

05189-2524          0.59$\times$0.49      440     34     0.042731     101.1$\pm$0.4     9.33     9.26     183$\pm$5     93

07251-0248     E     0.36$\times$0.31      332     41     0.087787                        55

W                 0.087817                       271

tot                      68.8$\pm$ 0.7     9.77     9.70     376$\pm$5

09022-3615          0.34$\times$0.30      348     37     0.059577     329.4$\pm$0.8     10.14     10.07     376$\pm$5     0

10190+1322     E     0.36$\times$0.31      330     49     0.075970     61.2$\pm$0.2     9.61     9.54     416$\pm$5     250

W                 0.076626     53.2$\pm$ 0.8     9.55     9.48     325$\pm$5     113

11095-0238     NE     0.31$\times$0.25      460     121     0.106535                        6

SW                 0.106302                        -

tot                      40.7$\pm$ 0.1     9.71     9.64     264$\pm$5

12071-0444     N     0.22$\times$0.16      367     62     0.128969     40.8$\pm$0.6     9.86     9.79     183$\pm$5     79

S                 0.128441     4.6$\pm$ 0.2     8.92     8.85     335$\pm$15     172

12112+0305     NE     0.34$\times$0.30      265     74     0.072717     122.3$\pm$0.5     9.88     9.81     335$\pm$5     80

SW                 0.073203     27.8$\pm$ 0.4     9.23     9.17     274$\pm$5     290

13120-5453          0.65$\times$0.65      1065     19     0.031114     650.7$\pm$2.1     9.87     9.81     335$\pm$5     91

13451+1232     W     0.23$\times$0.17      362     81     0.121680     41.5$\pm$0.6     9.82     9.76     386$\pm$15     242

E                 0.121291     2.8$\pm$ 0.4     8.65     8.58     223$\pm$36      -

14348-1447     SW     0.35$\times$0.30      345     79     0.082750     116.1$\pm$0.6     9.95     9.89     193$\pm$5     231

NE                 0.082351     64.9$\pm$ 1.0     9.70     9.63     213$\pm$5     196

14378-3651          0.41$\times$0.28      554     47     0.068097     63.7$\pm$0.6     9.52     9.45     183$\pm$5     211

15327+2340          1.28$\times$0.81      2145     2     0.018120     1693.7$\pm$3.5     9.84     9.77     477$\pm$5     41

16090-0139          0.20$\times$0.17      436     49     0.133690     75.5$\pm$1.0     10.15     10.09     345$\pm$10     181

16156+0146     NW     0.24$\times$0.16      498     76     0.132985     6.2$\pm$0.3     9.06     8.99     223$\pm$25     323

SE                  -     < 0.05$^\star$     < 6.96     < 6.89     -      -

17208-0014          0.47$\times$0.47      1043     7     0.042800     551.4$\pm$2.5     10.06     10.00     427$\pm$5     109

19297-0406     N     0.34$\times$0.31      282     89     0.085407                        230

S                 0.085785                        -

tot                      140.3$\pm$ 0.5     10.05     9.98     386$\pm$5

19542+1110          0.40$\times$0.34      368     44     0.062479     40.9$\pm$0.3     9.28     9.21     274$\pm$5     229

20087-0308          0.31$\times$0.25      287     192     0.105461     94.6$\pm$0.3     10.05     9.98     528$\pm$5     250

20100-4156     SE     0.20$\times$0.14      336     49     0.129848     31.3$\pm$0.3     9.74     9.68     376$\pm$5     228

NW                 0.129740     2.2$\pm$ 0.2     8.59     8.52     111$\pm$6     268

20414-1651          0.24$\times$0.18      375     47     0.086874     51.4$\pm$0.6     9.62     9.56     427$\pm$5     59

22491-1808     E     0.44$\times$0.33      342     58     0.077742     59.4$\pm$0.2     9.59     9.52     284$\pm$5     348

W                 0.077560     4.1$\pm$ 0.2     8.40     8.34     121$\pm$5      -

(1) IRAS name.
(2) Name of the nucleus.
(3) FWHM of the synthesised ALMA beam.
(4) 1$\sigma$ line rms sensitivity per 10 channel.
(5) Integration time on source.
(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  (column 7) and FWHM (column 8).
$^\star$For the undetected nucleus 16156+014546 SE, we provide a 3$\sigma$ upper limit, calculated based on the rms and a typical line FWHM of 300 .
(8) CO(2-1) luminosity, calculated as $L'_{CO}= 3.25\cdot 10^7 S_{CO} \nu_{rest}^{-2} D_L^2 (1+z)^{-1} $, where  $\nu_{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 $M(H_2) = 1/r_{21} \cdot \alpha_{CO} \cdot  L'_{CO(2-1)}$, where $r_{21}= L'_{CO(2-1)}/L'_{CO(1-0)}=0.91$ [Bolatto et al., 2013] and $\alpha_{CO}=0.78$ $M_\odot$/(K  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.

§ ANALYSIS

§.§ 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 <ref>)  are in good agreement with previously reported redshifts.
Throughout this paper,  we use the radio velocity definition.

Figure <ref> shows the CO(2-1) moment 0, 1, and 2 maps,  along with  the CO(2-1) peak map and continuum map  for four targets as an example.  The maps of the full sample are shown in the Appendix (Fig. <ref>).

The CO(2-1) emission is detected with S/N $ > 3$
in all individual nuclei but 16156+0146 SE.
In the majority of the nuclei (32/37), the kinematic major axis is clearly visible in the moment 1 maps. The moment 1 maps of five nuclei show a less ordered motion (01572+0009, 05189-2524,  12112+0305 NE,   14348-1447 SW,  22491-1808 E+W).
In 22/37 nuclei the moment 2 map shows an increase in the velocity dispersion close to the peak continuum position.

In  09022-3615, we detect an increase in velocity dispersion south of the nucleus (distance $\sim 1$",  equivalent to $\sim 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/  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.

In the appendix <ref>, we also show the maps of the CO(2-1) emission integrated in 50  channels for all targets.

Examples of the ALMA $\sim220-250$ GHz continuum and CO(2-1) moment maps (from left to right: peak map,  moment 0,  1,  and 2)  for a merger (20087-0308) and an interacting (10190+1322) starburst dominated system  and  for a merger (01572+0009) and an interacting (12071-0444) AGN dominated  system. The blue crosses mark the position of the nuclei (see Table  <ref>).
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$\times\sigma$ (where $\sigma$ is the rms) and [0.1, 0. 2,  0.4,  0.6,  0.8] of the maximum;
moment 0: [3, 6, 25, 50, 75]$\times\sigma$ (where $\sigma$ is the rms),
moment 1: every  50 (every 25  if the maximum value $< 100$ ),
moment 2: every 25  (every 15   if the maximum value $< 150$ ).

§.§ Spectro-astrometry

We perform a spectro-astrometry analysis  to identify high-velocity 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.

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\times5$ 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 $\Delta x= FWHM_{beam}/(2\times S/N)$, where $FWHM_{beam}$ is the beam size and S/N is the signal-to-noise ratio of the binned channel [Condon, 1997].
The spectro-astrometry analysis of  13120-5453 and 20100-4156 SE  are shown in Figure <ref> and  <ref>  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 Fig. <ref>, <ref>,  and <ref>).  As low-velocity channels,  we consider absolute velocities $|v|< 300$ .
If there are channels at $|v|< 300$  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 Fig. <ref>, <ref>,  and <ref>).   In general,  we consider as high-velocity the channels with  $|v|> 300$ , 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 high-velocity 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,  $\sim200-300$ , it is more likely due to a tidal tail than to an outflow  (e.g.  00091-0738 N, 01572+0009).

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).
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  (every 25  if the maximum value $< 100$ ),  and  every 25  (every 15   if the maximum value $< 150$ ) on the moment 2 map.  In black  are the CO(2-1) moment 0  contours ([3, 6, 25, 50, 75]$\times\sigma$).
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$\sigma$ level. The next contour levels are (0.5, 0.7, 0.9) of the peak of the  emission,  if these are above the 3$\sigma$ level.  The dashed circle shows the size of the outflow ().
Panel e) CO(2-1) continuum-subtracted spectrum  extracted from a circle with radius equivalent to the outflow size (). The lower panel shows an y-axis zoom-in to highlight the emission in the wings.  The horizontal dashed line shows the 1$\sigma$ noise level.
The vertical dashed lines indicated the `flux-weighted' velocity of the blue and red-shifted outflow ($v_{out}$).
Panel f) OH 119 spectrum (upper) compared with the nuclear CO(2-1) spectrum (bottom),  convolved to the resolution of the OH spectrum  (FWHM$\sim 270$ ).  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. <ref>). 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.

Same as Figure <ref>,  but for 20087-0308,  a target with outflow direction not perpendicular to the kinematic major axis.

§.§.§ 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 <ref>). 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$\alpha$) disks presented in Perna et al., 2022.
Figure <ref> 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 $\sim~20^{\circ}$,  hence it is  consistent with what observed in non-interacting galaxies  <cit.>.
We identify only six outliers with PA differences $\gg 20^{\circ}$ (01572+0009, 07251-0248 E, 09022-3615,14348-1447 SW,  16090-0139, and 17208-0014).
The PAs of CO are measured on smaller scales  ($\sim$1 kpc),   compared to the scales used to derive the PAs of the stellar disk and ionised gas ($\sim$5-10 kpc). Thus,  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 <cit.>.
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.

Absolute difference between the position angles (PAs) of the kinematic major axis derived from CO and the PAs derived from the stellar kinematics (left), and the ionised gas kinematics (right).  The PAs of the stars and ionised gas are taken from Perna et al., 2022.
The dashed lines show a difference of  $20^{\circ}$. The names of galaxies with large PA differences ($\gg 20^{\circ}$) are shown on the figure.  The colours and shapes of the symbols are as in Fig. <ref>.

§.§ Outflow velocity range definition

One of the main goals of this work is to identify and characterise high-velocity ($> 300$ ) 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 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$ .
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 are in the range $|v_{max}|=300-800$ .
Figure <ref>,  <ref>,  and <ref>  show the CO(2-1) spectra with the blue-shifted and red-shifted outflow velocity ranges highlighted with blue and red shaded areas.
The integrated maps of the blue- and red-shifted outflow channels are shown as contours in the bottom-left panel of the figures.

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 <cit.>.  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 (),  and molecular gas mass in the outflow ().  We use these parameters  to derive the outflow energetics: mass outflow rate (),  mass-loading factor ($\eta=\dot{M}_{out}/SFR$), outflow momentum rate ($\dot{P}_{out}$),  and kinetic luminosity ($L_{out}$).
In the following, we explain how we measure the `projected'  and  .  We discuss the inclination corrections in Sec. <ref>.

Different methods have been used in the literature to separate the outflow and rotating disk emission.
A possible method <cit.> 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 <cit.>. 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 line profiles of our targets  cannot be well fitted using a simple model with only one or two Gaussians (see Fig. <ref>,  <ref> and  <ref>).  The line profiles are asymmetric and show multiple peaks, which could also include self-absorption (see Sec. <ref>).
Thus, to measure the outflow gas mass,  we consider the total flux in the high-velocity channels,   highlighted  in the blue and red  shaded regions  on the spectra (see Fig.  <ref>, previous section), without subtracting the low-velocity Gaussian fit.

In Section <ref>),  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 Section <ref>.

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. <ref>),  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:

\begin{equation}
R_{out}^b = d_c^b + \frac{FWHM^b}{2},
\end{equation}

where $d_c^b$ 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 blue- and red-shifted wings:

\begin{equation}
R_{out} = 0.5 \cdot(R_{out}^{b}+R_{out}^{r}),
\end{equation}

The outflow radii are shown in Fig. <ref>  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$\sigma$ 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\sigma}$).

The outflow radii  $R_{out}$ measured from the 2D Gaussian fit are in the range $0.18-0.94$ kpc ($0.1-1.5$").
The maximum outflow radii  $R_{3\sigma}$  are in the range $ 0.1-2.1$ kpc ($0.05-2.5$").
The  ratio of the observed $R_{3\sigma}/R_{out}$ is in the range 1-3, with a mean of 1.45.  We note that in some cases the $R_{3\sigma}$ values reported in Table <ref> 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 within $R_{3\sigma}$ instead of within $R_{out}$ to calculate the outflow mass,   would increase on average by a factor of 1.3 (0.11 dex).
The mass outflow rate is proportional to /, thus the larger molecular mass included within $R_{3\sigma}$ is counterbalanced by the larger radius.
If we were to use $R_{3\sigma}$ and the outflow flux within this radius, we would obtain very similar values (less than 8% smaller,  $-0.04$ dex)  compared to the   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 Sec. <ref>).

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  and ,  separately for the blue and the red part of the outflow.  The central spectra are shown in Fig. <ref>,  <ref> and  <ref> (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 ($\sigma$ in units of mJy) .  The  uncertainty on $F_{out}$ is :

\begin{equation}
F_{out, err} = \sigma \cdot \Delta v_{ch} \cdot \sqrt{N_{ch}},
\label{eq:F_err}
\end{equation}

where $ \Delta v_{ch}$ is the width of a velocity channel in , 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  pc$^2$)  using the formula:

\begin{equation}
L'_{CO}= 3.25\cdot 10^7 S_{CO} \nu_{rest}^{-2} D_L^2 (1+z)^{-1},
\end{equation}

where $S_{CO}$  is the velocity-integrated CO line flux in Jy ,   $\nu_{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 convert the CO(2-1) luminosity to CO(1-0) luminosity using $r_{21} = L'_\text{CO(2-1)}/L'_\text{CO(1-0)} = 0.91$ [Bolatto et al., 2013]. Then we multiply it by the ULIRGs-like CO-to-H$_2$ conversion factor $\alpha$= 0.78 M$_\odot$/(K  pc$^2$),  to obtain the outflow  molecular (H$_2$) gas mass .
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 <cit.>.

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 Sec. <ref>,  Fig. <ref>,  middle panel of the bottom row):

\begin{equation}
v_{out} =  \frac{ \sum_{i} v_i \cdot F_{i}}{ \sum_{i}  F_{i}},
\label{eq:vout}
\end{equation}

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 Sec. <ref>).
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 Sections <ref> and <ref> 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 [  yr$^{-1}$]) using the formula:

\begin{equation}
\label{eq:M_rate}
\dot{M}_{out} =  \frac{|v_{out}| \cdot M_{out}}{ R_{out}}= \frac{ \sum_{i} |v_i| \cdot M_{i}}{ R_{out}},
\end{equation}

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. <ref>). The total mass outflow rate is the sum of the blue and red $\dot{M}_{out}$.
We note that this formula corresponds to the assumption that the outflow has started at a point in the past  ($-t =-$/) and has continued with a constant   [Rupke et al., 2005, Veilleux et al., 2005, Lutz et al., 2020].
Under this assumption,  the average volume density  of the outflowing gas decreases with radius ($\propto R^{-2}$).
Assuming that the outflowing gas fills a spherical or multi-conical volume with  a constant average volume density,  would increase $\dot{M}_{out}$ by a factor of three <cit.>.

In the cases where the outflow is not detected,  we estimate 3$\sigma$ upper limits on $\dot{M}_{out}$ as:

\begin{equation}
\dot{M}_{out, UL} =  3\cdot\frac{\langle v_{out}\rangle \cdot M_{out, err}}{\langle R_{out}\rangle},
\end{equation}

where $\langle v_{out}\rangle= 390$ and $\langle R_{out}\rangle=0.52$ kpc are the median outflow velocity and radius of our sample.
$M_{out, err}$ is calculated based on $F_{out, err}$ (eq. <ref>),  measured from the spectrum extracted from a radius $\langle R_{out}\rangle$.
The mass outflow rates are in the range $\sim5$ to $\sim300$  yr$^{-1}$.

We compare the outflow properties (,  and ) measured from the blue and red-shifted wings.  The ,  and  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\pm0.16$, $1.05\pm0.51$,  $0.99\pm0.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  and , 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 Sec. <ref> 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   and  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  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/\left<\sin(i)\right>=1/0.79=1.27$.  Analogously,  we calculated the average correction for the radius: $1/\left<\cos(i)\right>=2$.
We use these values to correct the mean outflow velocity and mean outflow radius reported in table <ref>.

It is not possible to derive the average inclination correction for the mass outflow rate  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} \propto M_{out}/\dot{M}_{out})$ is unity <cit.>.

CO(2-1) observed outflow properties.

IRAS name      n.      [$v_{min}, v_{max}$] (b)       [$v_{min}, v_{max}$] (r)      $S_{CO}$ (b)     $S_{CO}$ (r)    $R_{3\sigma}$     $R_{out}$      $v_{out}$     $\log M_{out}$       $\log \dot{M}_{out}$     angle$_{out}$

[]      []     [Jy ]     [Jy ]      [kpc]      [kpc]     []      []      [yr$^{-1}$]     [deg.]

(1)     (2)      (3)     (4)    (5)     (6)     (7)     (8)     (9)     (10)     (11)     (12)

00091-0738     S     -     -    -     -     -     -     -    <7.26     <1.15     -

N     -     -    -     -     -     -     -    <7.20     <1.08     -

00188-0856           -     -    -     -     -     -     -    <7.33     <1.22     -

00509+1225           -     -    -     -     -     -     -    <6.62     <0.50     -

01572+0009           -     -    -     -     -     -     -    <7.62     <1.51     -

05189-2524           [-200,-600]     [200,350]     1.27$\pm$0.05     1.27$\pm$0.04     1.39     0.83$\pm$0.01     294$\pm$16     7.66$\pm$0.01     1.22$\pm$0.01     -

07251-0248     E     [-320,-500]     [300,550]     -     0.32$\pm$0.01     0.64     0.17$\pm$0.01     377$\pm$12     7.37$\pm$0.01     1.72$\pm$0.01     -

W     -     -    -     -     -     -     -    <6.87     <0.75     -

09022-3615           [-400,-700]     [300,500]     6.87$\pm$0.05     2.05$\pm$0.04     1.68     0.94$\pm$0.01     423$\pm$8     8.50$\pm$0.01     2.17$\pm$0.01     12

10190+1322     E     [-290,-410]     [350,420]     0.06$\pm$0.02     0.09$\pm$0.02     0.60     0.62$\pm$0.02     366$\pm$143     6.96$\pm$0.08     0.74$\pm$0.08     75

W     -     -    -     -     -     -     -    <6.78     <0.66     -

11095-0238     NE     [-300,-500]     [300,500]     0.25$\pm$0.02     0.22$\pm$0.02     0.46     0.37$\pm$0.01     373$\pm$37     7.70$\pm$0.02     1.71$\pm$0.02     32

SW     -     -    -     -     -     -     -    <7.05     <0.94     -

12071-0444     N     [-250,-430]     [260,400]     0.22$\pm$0.03     0.22$\pm$0.03     0.36     0.40$\pm$0.01     385$\pm$65     7.83$\pm$0.04     1.82$\pm$0.04     37

S     -     -    -     -     -     -     -    <7.37     <1.26     -

12112+0305     NE     [-250,-700]     [250,700]     2.43$\pm$0.03     4.09$\pm$0.03     2.07     0.59$\pm$0.01     357$\pm$6     8.54$\pm$0.01     2.33$\pm$0.01     90

SW     [-200,-600]     [200,400]     0.51$\pm$0.04     0.29$\pm$0.03     0.73     0.63$\pm$0.01     292$\pm$32     7.63$\pm$0.03     1.31$\pm$0.03     84

13120-5453           [-300,-800]     [300,600]     5.38$\pm$0.15     2.48$\pm$0.13     0.70     0.37$\pm$0.01     432$\pm$25     7.89$\pm$0.01     1.96$\pm$0.01     29

13451+1232     W     -     -    -     -     -     -     -    <7.62     <1.50     -

E     -     -    -     -     -     -     -    <7.22     <1.10     -

14348-1447     SW     [-250,-700]     [250,700]     1.46$\pm$0.04     1.59$\pm$0.04     1.78     0.67$\pm$0.01     380$\pm$13     8.31$\pm$0.01     2.07$\pm$0.01     56

NE     [-300,-550]     [300,600]     0.57$\pm$0.03     0.46$\pm$0.03     1.36     0.69$\pm$0.01     393$\pm$35     7.83$\pm$0.02     1.60$\pm$0.02     82

14378-3651           [-200,-600]     [200,300]     1.16$\pm$0.05     0.32$\pm$0.02     0.69     0.51$\pm$0.01     311$\pm$24     7.82$\pm$0.02     1.61$\pm$0.02     -

15327+2340           [-450,-600]     [400,600]     4.06$\pm$0.15     2.08$\pm$0.18     0.96     0.56$\pm$0.01     474$\pm$43     7.33$\pm$0.02     1.26$\pm$0.02     64

16090-0139           [-400,-600]     [350,600]     0.33$\pm$0.04     0.70$\pm$0.04     0.69     0.50$\pm$0.01     454$\pm$58     8.22$\pm$0.02     2.19$\pm$0.02     -

16156+0146     NW     -     -    -     -     -     -     -    <7.32     <1.21     -

SE     -     -    -     -     -     -     -    <7.47     <1.35     -

17208-0014           [-400,-800]     [390,450]     3.18$\pm$0.17     0.18$\pm$0.06     0.49     0.51$\pm$0.01     487$\pm$148     7.78$\pm$0.03     1.77$\pm$0.03     3

19297-0406     N     [-350,-500]     [350,500]     0.28$\pm$0.02     0.27$\pm$0.02     1.01     0.64$\pm$0.01     409$\pm$42     7.58$\pm$0.02     1.40$\pm$0.02     49

S     -     -    -     -     -     -     -    <6.76     <0.64     -

19542+1110           [-250,-400]     [300,500]     0.19$\pm$0.02     0.16$\pm$0.02     0.48     0.38$\pm$0.01     346$\pm$52     7.14$\pm$0.03     1.11$\pm$0.03     72

20087-0308           [-350,-450]     [380,600]     0.22$\pm$0.01     0.39$\pm$0.02     1.12     0.67$\pm$0.01     425$\pm$32     7.79$\pm$0.01     1.60$\pm$0.01     40

20100-4156     SE     [-300,-800]     [250,800]     0.90$\pm$0.04     1.05$\pm$0.04     0.76     0.44$\pm$0.01     438$\pm$26     8.47$\pm$0.01     2.48$\pm$0.01     70

NW     -     -    -     -     -     -     -    <7.39     <1.28     -

20414-1651           [-350,-440]     [300,400]     0.18$\pm$0.03     0.25$\pm$0.03     0.44     0.52$\pm$0.01     364$\pm$84     7.48$\pm$0.05     1.33$\pm$0.05     25

22491-1808     E     [-200,-400]     [200,600]     1.69$\pm$0.01     0.68$\pm$0.01     0.92     0.26$\pm$0.01     263$\pm$5     8.12$\pm$0.01     2.14$\pm$0.01     57

W     -     -    -     -     -     -     -    <6.60     <0.49     -

average$^{*}$                         1.84     1.07     485     8.01     1.89

(1) IRAS name.
(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$\sigma$ (see Sec. <ref>).
(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. <ref>),  not corrected for inclination.
(10) Outflow molecular gas mass,  calculated  assuming a ULIRG-like conversion factor $\alpha_{CO}$ of 0.78 /(K  pc$^{-2}$)$^{-1}$ and $r_{21}$ ratio of 0.91 [Bolatto et al., 2013].
(11) Mass outflow rate calculated using equation <ref>.
(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  and   have been corrected for inclination assuming the average corrections: $ R_{out}= R_{out}(obs)/\left< \cos(i)\right> = 2\cdot R_{out}(obs)$,  $v_{out}= v_{out}(obs)/\left< \sin(i) \right> = 1.27\cdot v_{out}(obs)$ (see Section <ref>).

§.§.§ Caveats: Methods to measure outflow parameters

In this section, we discuss how the outflow parameters (in particular  ,  and ) 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 Section <ref>.

1) Low projected velocity ($ v < 300$ ) gas in the outflow: since we are not considering low projected velocities ($|v|\lesssim 200-300$ ), it is possible that we are missing part of the outflow flux.
If we were to consider also this low velocities, the outflow mass  would increase and the flux-weighted outflow velocity  would decrease.
Overall, we expect  to increase, but by a lower factor than , because the increase in  is counterbalanced by the decrease of .

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$^{\circ}$.
We set the maximum velocity of the outflow to 750 .
We consider several outflow radial velocity profiles, motivated by previous works in the literature <cit.>, and different outflow inclinations with respect to the line of sight.  More details on the simulations can be found in the appendix <ref>.
We measure ,  and  from the total simulated profile (outflow+systemic component) considering only $|v| > 300$ , to mimic the method we are using with our data.  Then, we measure the outflow quantities (, ,  and ) from the simulated outflow emission profile (without the systemic component) considering the full velocity range.

The   measured from $|v| > 300$  are underestimated compared to the values measured from the entire velocity range by a factor of $0.2-1$ (average 0.5), while the   are overestimated by up to a factor of 2.2 (average of 1.6).
Consequently,   would be underestimated by up to $\sim$0.45 dex (65%) for outflow inclinations close to the plane of the sky (90$^{\circ}$). For an inclination of 10$^{\circ}$,  the measured   would be underestimated  by up to 0.1 dex (see Fig.  <ref>).

Simulated outflow profiles for outflow  inclinations with respect to our line of sight of 10$^{\circ}$ (upper) and  90$^{\circ}$ (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 <ref>.
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  measured from the total profile at $|v|\geq300$   (coloured in blue and red).
The vertical solid lines show the flux-weighted  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$ )  increases.

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$ (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$  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  (, , and ) starting from $v_{min}=300$ .
We find that  the value of  increases by 0.28 dex on average (maximum 0.67 dex), while  decreases by a factor of $-0.07$ dex on average (minimum $-0.11$ dex).  So, the  estimates increase by $0.2$ dex on average (maximum 0.6 dex).

2) Possible overestimation of  due to the rotating disk contribution at $|v|=250-300$ : Since we do not model and  subtract the disk rotation,  it is possible that at low velocities ($250-350$ ) 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 $\sim300$ ) with one,  two or three Gaussians <cit.>.
Then, we compute the outflow parameters (, , ) 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  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 $\pm 50$ .
The   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'  <cit.>.
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., 2005, used for example by Fluetsch et al., 2019:

\begin{equation}
v_{max, FWHM} = | \Delta v| + \frac{FWHM_{broad}}{2},
\label{eq:vout_Fluetsch2019}
\end{equation}

where $|\Delta 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:

\begin{equation}
v_{max, FW10} = | \Delta v| +\frac{FW_{10\%}}{2},
\label{eq:vout_Lutz2020}
\end{equation}

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 Fig. <ref> we compare the flux-weighted  with the maximum velocity  $v_{max, FWHM}$  and $v_{max, FW10}$ for our sample.  The ratio $v_{max, FWHM}$/ varies between  0.5 and 2.0, with a median 0.9.  The ratio $v_{max, FW10}$/ 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)  would need to be corrected  by an average factor of 1.3 (see Sec. <ref>).  This factor can account for most of the average difference between $v_{max, FW10}$ and .

Comparison of different methods for measuring the outflow velocity.  Left: Ratio between the maximum outflow velocity $v_{max, FWHM}$ (eq.  <ref>) and the  flux-weighted outflow velocity ($v_{out}$) in our sample.
Right: Ratio between  $v_{max, FW10\%}$  (eq.  <ref>) and . The horizontal lines show the 1-to-1 ratio.

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$\sigma$ level.  In this way,  we may underestimate , 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 <ref>, we summarise the average biases  in the outflow properties (,  and ) 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  $\sim 0.2$ dex and  by $\sim 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  and  we adopt a typical uncertainty of 0.3 dex. For , we estimate a typical uncertainty of 0.2 dex, based on the difference between $R_{out}$ and $R_{3\sigma}$.

The uncertainties on $\dot{M}_{out}$ are dominated by the uncertainties on the  $\alpha_{CO}$ conversion factor,  which can be up to 0.7 dex  <cit.>,
and on the outflow geometry.
If the $\alpha_{CO}$ is more similar to the Galactic value ($\alpha_{CO}$=4.3 M$_\odot$/(K  pc$^2$)), it would imply that all our  and  measurements are underestimated by a factor of $\sim5$ (0.7 dex), while if the optically thin case applies ($\alpha_{CO}= 0.35$ M$_\odot$/(K  pc$^2$)), our measurements would be overestimated by a factor of $\sim2$ [Bolatto et al., 2013].

Possible biases in the outflow properties measured with our method, due to different issues.

3c biases in outflow properties

[dex]    [dex]      [dex]

4l 1) Not considering outflow flux at $|v|<300$  :

average:     $0.20$      $-0.30$      $-0.15$

range:     $(0, 0.34)$       $(-0.70, 0)$     $(-0.45, 0)$

4l2) Rotating disk contribution at $|v| = 250 -300$ :

average:     $-0.01$      $0.12$      $0.11$

range:     $(-0.08, 0, 08)$     $(-1.0, 0.5)$      $(-0.82,  0.45)$

4l3) Flux-weighted velocity instead of $v_{max}$:

average:     $-0.20$      $-$      $-0.20$

range:      $(-0.50, 0.30)$     $-$      $(-0.50, 0.30)$

§.§ 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 $\sim$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. <ref>.
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} = |\Delta R| + FWHM/2$, where $|\Delta 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 $\sim2$ lower than in our work.
For the other 41  targets,  they collect information about the outflow parameters (, , ) 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 (, , , ) with these works in Figure <ref>.
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 , while we measure  $v_{out}< 500$ .  For 20100-4156, it is possible that the spectral range of our observations ($v$ = [-1200, 1200] ) is not wide enough to detect the emission at very high velocities  ($|v_{out}|>1000$ ).  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,  values from the literature are considerably higher than our measurements ($\sim2-6$ times higher).   For three of these galaxies, the difference is due to the different definition of  as the maximum radius at which the outflow is detected <cit.>.  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  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  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   is proportional to   and $^{-1}$,  the differences in  and  tend to balance each other and lead to similar  (within a factor of 3).

Comparison of outflow parameters measured in this work with values reported in Lutz 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.

§.§.§ Comparison with Fluetsch et al., 2019

We also compared the  and   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. <ref>.
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 blue-shifted and red-shifted emission in the calculation of , thus their  estimates are expected to be smaller.

Figure <ref> 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  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.2 dex). This difference can be attributed to the different method used to estimate the flux belonging to the outflowing gas. This difference in  propagates to the .

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 <ref> shows   and  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_{\odot}= 12.0-12.5$),  their average   and   are $\sim0.6-0.8$ dex higher than our measurements (factor of $\times4-6$).

Given the different methodology used to measure these parameters, the difference is not surprising. As discussed in Section <ref>, this comparison highlights the large effect that the choice of methodology can have on the measured outflow parameters.
The approach used by Fluetsch et al., 2019 assumes that most of the outflow flux is at low projected velocity.   Based on simulations of biconical outflows <cit.>,  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.

Comparison of outflow parameters measured in this work with values reported in Fluetsch et al., 2019.From left to right: outflow velocity, outflow radius, outflow mass and mass outflow rate.  Differences in  and  can be up to a factor of 16.

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  and  values estimated by  Fluetsch et al., 2019 are higher by a factor of $\times5-8$ compared with our measurements at the same infrared luminosity.

§.§ Analysis of the OH 119 spectra

Other lines that are often used as tracer of molecular outflows are the OH (hydroxyl) FIR lines   <cit.>. We compare the outflow parameters derived from CO(2-1) with the ones derived from the  OH 119 doublet, which is the strongest of the OH ground-state lines <cit.>.
We collect archival /PACS  spectra of the OH 119  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 in Veilleux et al., 2013, Spoon et al., 2013, and González-Alfonso et al., 2017; 11095-0238 and  13451+1232 are not presented in these works but were found in the  archive.

We extract the spectra from the central $9.4''\times 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 <ref>), 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 <cit.>.
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 $\Lambda$-doublet is fixed to 0.208 (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 /PACS  instrumental resolution <cit.>,
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 <cit.>, 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 lines), then 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$\sigma$ uncertainties on the derived parameters. The 50 MC realisations are shown as orange curves on Figure <ref>.
The velocities derived from the OH profiles are shown in Table <ref>.
The properties of molecular outflows derived from OH will be discussed and compared with those derived from the CO tracer in Sec. <ref>.

Example of the fit of the OH 119 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 $\mu$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 dot-dashed 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 $\sim2000$  show the position of the $^{18}$OH 120 $\mu$m doublet and the CH$^{+}$(3-2) absorption line.

Velocities derived from the OH 119$\mu$m profiles.

IRAS name     $v_{50}$(abs)     $v_{84}$(abs)     $v_{98}$(abs)    $v_{50}$(em)      $v_{84}$(em)     $v_{98}$(em)

[]    []        []     []        []      []

(1)     (2)      (3)     (4)    (5)     (6)     (7)

00091-0738$^{\star}$       -     -     -      -     -     -

00188-0856     -330$\pm$29     -774$\pm$69     -1186$\pm$136      -     -     -

00509+1225      -     -     -     120$\pm$33     262$\pm$85     404$\pm$194

01572+0009      -     -     -     88$\pm$139     155$\pm$209     188$\pm$318

05189-2524     -306$\pm$84     -543$\pm$72     -781$\pm$105     145$\pm$23     316$\pm$22     489$\pm$33

07251-0248     -127$\pm$20     -347$\pm$22     -475$\pm$126      -     -     -

09022-3615     -139$\pm$36     -270$\pm$28     -398$\pm$72     214$\pm$11     342$\pm$14     472$\pm$22

10190+1322$^{\star}$       -     -     -      -     -     -

11095-0238      -     -     -     206$\pm$163     344$\pm$40     413$\pm$110

12071-0444     -249$\pm$17     -454$\pm$30     -694$\pm$43     194$\pm$8     296$\pm$12     398$\pm$17

12112+0305     21$\pm$12     -76$\pm$27     -207$\pm$48     341$\pm$263     503$\pm$327     664$\pm$405

13120-5453     -253$\pm$16     -586$\pm$18     -850$\pm$35      -     -     -

13451+1232      -     -     -     -51$\pm$35     155$\pm$123     360$\pm$274

14348-1447     -245$\pm$34     -560$\pm$40     -917$\pm$68     154$\pm$11     300$\pm$14     446$\pm$23

14378-3651     -296$\pm$27     -650$\pm$39     -1007$\pm$58     189$\pm$15     319$\pm$14     478$\pm$20

15327+2340     -28$\pm$0     -194$\pm$0     -359$\pm$13      -     -     -

16090-0139     -36$\pm$25     -410$\pm$43     -683$\pm$91      -     -     -

16156+0146$^{\star}$        -     -     -      -     -     -

17208-0014     50$\pm$10     -99$\pm$15     -466$\pm$94      -     -     -

19297-0406     -121$\pm$39     -500$\pm$84     -879$\pm$199      -     -     -

19542+1110     -302$\pm$50     -480$\pm$96     -622$\pm$157     267$\pm$72     302$\pm$114     373$\pm$168

20087-0308     48$\pm$19     -160$\pm$24     -403$\pm$52      -     -     -

20100-4156     -393$\pm$39     -985$\pm$73     -1539$\pm$145      -     -     -

20414-1651     17$\pm$0     -53$\pm$16     -88$\pm$85      -     -     -

22491-1808     88$\pm$20     25$\pm$55     -8$\pm$196      -     -     -

(1) IRAS name.
(2),  (3) and (4): 50th, 84th and 98th percentile velocities derived  from the absorption profile (when present).
(5),  (6) and (7): 50th, 84th and 98th percentile velocities derived  from the emission profile (when present).
$^{\star}$: OH 119$\mu$m spectra not available for 00091-0738, 10190+1322, and 16156+0146.

§ RESULTS

§.§ Mean outflow properties

Here we summarise the ranges of outflow properties of our sample and their average values,  corrected for inclination as described in Sec. <ref>.
We measure projected outflow velocities of $\sim260-490$ , with a mean outflow velocity (corrected for inclination) of $485\pm16$ . Outflow radii are in the range 0.17-0.94 kpc and the mean inclination-corrected outflow radius is $1.1\pm0.1$ kpc.
The outflow masses are between $1-35  \times 10^7$ , with an average of $(10\pm2)  \times 10^7$ . 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$  yr$^{-1}$, with an average outflow rate of $78\pm16$  yr$^{-1}$.

As a comparison,  Ramos Almeida et al., 2022 find molecular gas mass outflow rate $=8-16$  yr$^{-1}$ in a sample of nearby type 2 quasars ($\log L_{AGN}/[\text{erg s}^{-1}] = 45.7-46.3$).  Lower   $=0.3-5$  yr$^{-1}$ have been measured in lower luminosity ($\log L_{AGN}/[\text{erg s}^{-1}] = 43.2-44.2$) Seyfert galaxies [Alonso-Herrero et al., 2019, Domínguez-Fernández et al., 2020, García-Bernete et al., 2021].  Higher molecular gas mass outflow rates ($=60-400$  yr$^{-1}$) have been measured in ULIRGs hosting an AGN [Feruglio et al., 2010, Cicone et al., 2014]. We note that, as discussed in Sec. <ref>, the method used to derive the mass outflow rates can introduce systematic differences (up to factor of $\sim10$) 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 , with a mean $393\pm 38$  and mass outflow rates in the range  $1-100$  yr$^{-1}$ (from Fig. 13 in the paper), with an average 38  yr$^{-1}$.  The average  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  \times 10^7$ , with a average of  $6.7  \times 10^7$ ,  similar to the masses of the molecular phase.

We find molecular outflow dynamical times ($t_{dyn} = R_{out}/v_{out} $) in the range $0.45-2.77$ Myr.

The mean $t_{dyn},$ based on the mean observed  and , 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.51 Myr reported in Pereira-Santaella et al., 2018.
The outflow depletion times ($ t_{dep} = M(H_2)_{tot}/\dot{M}_{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 in Pereira-Santaella et al., 2018 ($15-80$ Myr) and Cicone 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).

Mean cold molecular outflow properties in ULIRGs.

[]      [kpc]      [10$^7$ ]       [ yr$^{-1}$]

mean:      $485\pm16 $       $1.07\pm0.08$        $10\pm2$	     $78\pm16$

range:     $ 260-490^{\star}$     $0.26-0.94^{\star}$      $ 1-35 $    $6-302$

$^{\star}$: Not corrected for inclination.

§.§ Outflow characteristics for AGN/SB and interacting/mergers

§.§.§ Outflow detection rate and direction

In this Section, we present the statistics of the number of detected molecular outflows in our sample and we investigate whether there is any dependency of the outflow detection rate on the nuclear classification (AGN or SB) or on the merger stage (advanced mergers or interacting systems).

The top panel of Figure <ref> 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_{\odot} > 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_{\odot} = 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. <ref> presents the percentage of outflow detections in each category.  We detect an outflow, defined as high-velocity ($|v|>300$ ) CO(2-1) emission which deviates from the main rotation,  in 20/26 of the nuclei with $\log L_{IR}/L_{\odot} > 11. 8$ ($77\pm7\%$ [The uncertainties on the percentages represent the 90% binomial confidence intervals. ]).

The nuclei with outflow detections are  equally divided between mergers  and interacting systems.
The percentage of detections in SB ($93\pm4$%,  14/15) is higher than in AGN ($55\pm14$%,  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_{\odot}= 9.2\pm 0.1$), than SB nuclei ($\log M(H_2)M_{\odot}= 9.4\pm 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  <cit.>.
Similarly,   Stone et al., 2016 find a lower OH outflow detection fraction in X-ray selected AGN (24%) than in ULIRGs ($\sim70\%$) and suggest that outflow detection in ULIRGs may be easier due to their higher gas fraction.

With the spectro-astrometry analysis (see Sec. <ref>),  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 Sec. <ref>).
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 Fig. <ref>, <ref> and <ref>) to check that the direction is consistent with the position of the gas in the high-velocity channels.
The bottom panel of Fig. <ref> shows, for each category,  the fraction of outflows with 1)  direction perpendicular to the rotating disk (angle  $90\pm20^{\circ}$),  2) direction non-perpendicular,  or 3) unclear direction.
The percentages of outflow with determined direction in mergers ($70\pm12 \%$,   7/10) is smaller than in interacting nuclei ($90\pm6\%$,  9/10).
We could determine the outflow direction in  $86\pm7\%$ (12/14) of  SB nuclei with outflow detection,  but only in $67\pm15\%$ (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\pm20^{\circ}$) to the kinematic major axis. These nuclei are 12112+0305 NE,  12112+0305 SW,  14348-1447 NE  <cit.>,   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 <cit.>.
However,  there are also many SB nuclei ($50\pm13$%) 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 <cit.>.
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^{\circ}$ 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 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 show 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\pm20^{\circ}$),  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.

§.§.§ 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 (, , , ) 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$\sigma$) with the  sample mean (see Fig. <ref>).

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 <ref> shows the cumulative distribution functions (CDF) of the outflow properties (, , $\log M_{out}$,   $\log \dot{M}_{out}$,  and $\log \eta=\log (\dot{M}_{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. <ref>), in order to remove the effect of the correlation between  and $L_{IR}$ (see Sec. <ref>).  Also in this case there are no significant differences in the average /$L_{IR}$  of the different categories.
The average for the total sample is /$L_{IR}=(3.8\pm 0.6 )\times10^{-11}$  $yr^{-1}$ $^{-1}$.  In the right panel, we also plot the average outflow rates normalised by AGN luminosity /$L_{IR, AGN}$ for the AGN categories.  The average /$L_{IR, AGN}$ for AGN are $\sim3-6$ times higher than the  /$L_{IR}$ for starbursts.

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  ,  and $\eta$) . 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.

Left: Mean values of the ratio of  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  and $L_{IR,AGN}$.

§.§ Outflow properties and energetics

In this section,  we investigate the relation between the outflow properties and the total infrared luminosity,  as well as with the SFR and AGN luminosity.

§.§.§ 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_{IR,AGN}$)  as the fraction of the total infrared luminosity $L_{IR}$  due to AGN,  based on the MIR AGN fraction  <cit.>.   For the interacting systems,  $\alpha_{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 mass-loading factor ($\eta=$/SFR) to increase with AGN fraction ($\alpha_{AGN}=L_{AGN}/L_{bol}$),  while we do not find a correlation between these two quantities in our sample (Spearmann rank correlation coefficient for the detections $r=0.27$, $p$-value=0.28),  as it is shown in Figure <ref>. The correlation found in Cicone et al., 2014 is driven mostly by the objects with high AGN fraction ($>0.8$), which have $\eta>10$ .
Fluetsch et al., 2019 expand the  Cicone et al., 2014 sample and  find that the correlation between $\eta$ and $\alpha_{AGN}$ is  only evident for $\alpha_{AGN}> 0.7$. They propose that the lack of correlation for $\alpha_{AGN}< 0.7$  could be due to i) the contribution from star-formation to the mass outflow rate for small AGN fractions  or ii) to the short timescale of AGN variability <cit.>, compared to the outflow timescales ($10^6$ yr).
Since Cicone et al., 2014  included in the study objects with known molecular outflows,  it is possible that they sample the upper end of the distribution of $\eta$ at large $\alpha_{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 $\alpha_{AGN}$ estimated  from the MIR 30 $\mu$m to 15 $\mu$m flux ratio.  Cicone et al., 2014 estimate the AGN fraction using different methods:  from the 5-8 $\mu 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 $\alpha_{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 <ref> shows    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  more than one order of magnitude below the relation derived by Cicone et al., 2014 (dashed line on the Fig. <ref>).
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  versus AGN luminosity relation.    Our sample confirms that for the same AGN luminosity,  ULIRGs can have a wide range of .
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  values for a given AGN luminosity.

Mass-loading factor ($\eta=\dot{M}_{out}/$SFR) versus AGN fraction ($\alpha_{AGN}$).  Star symbols show the Cicone et al., 2014 sample (scaled down by a factor of three to match our outflow geometry definition).  Symbols for the PUMA sample  are  as in Fig. <ref>.

Mass outflow rate versus AGN luminosity.  Diamonds symbols show the sample from Cicone et al., 2014, for which  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  luminosity <cit.>.    The Lutz et al., 2020 sample is shown with red pentagons.  For the Stuber et al., 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. <ref>.
The dashed and dotted lines show the relations presented by Cicone et al., 2014 and by Fiore et al., 2017, respectively,  scaled down by a factor of three to match our outflow geometry definition).

§.§.§ Outflow energetics

In Figure <ref> we show the outflow rate,  momentum rate ($\dot{P}_{out} =\dot{M}_{out} \cdot v_{out}$),  and kinetic luminosity ($L_{out} = \frac{1}{2}\dot{M}_{out} \cdot v_{out}^2$) as a function of  SFR.  The SFR has been estimated from the IR luminosity, using the Kennicutt & Evans, 2012 relation (with Kroupa & Weidner, 2003 initial mass function (IMF)) after removing the fraction of luminosity associated with the AGN ($\alpha_{AGN}$).

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[http://www.phangs.org]) 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 \text{sSFR/yr}^{-1} > -11$) with stellar masses in the range $9.3 \leq \log M_{\star}/M_\odot \leq 11.1$.
Stuber et al., 2021 use  ALMA CO(2-1) observations  with a resolution of $\sim100$ 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  by measuring the outflow radii in each channel,  but the  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 $\alpha_{CO}=0.8$ M$_\odot$/(K  pc$^2$) [Bolatto et al., 2013],  similar to our assumption of $\alpha_{CO}=0.78$ M$_\odot$/(K  pc$^2$),  and $r_{21} =0.65$ [Bolatto et al., 2013, Leroy et al., 2013, den Brok et al., 2021, Leroy et al., 2021] while we use the value of  $r_{21} =0.91$ measured in ULIRGs  [Papadopoulos et al., 2012].
The SFR for this sample are derived by Leroy et al., 2021  from the WISE band 4 photometry, with a calibration consistent  with Kennicutt & Evans, 2012,  the one used for the PUMA targets.
The molecular outflows detected in the PHANGS galaxies have weighted velocities $v_{out}= 65-238 $ , $M_{out}= 0.35-102\times 10^6$  and $\dot{M}_{out} =0.15-21$  yr$^{-1}$.  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  Figure <ref> shows the relation between mass outflow rate and SFR.  The dashed lines indicate constant mass-loading factors ($\eta = \dot{M}_{out}/$SFR). The PUMA objects have $\eta < 0.04-1$. The PHANGS-ALMA galaxies have slightly higher $\eta =0.15- 4$.    The dispersion in  for a given SFR is quite high ($\sim1$ dex).
Mass-loading factors in local starburst galaxies are typically lower than $2-3$ <cit.>,
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  in some of the objects (see discussion in Sec. <ref>).

In Figure <ref>  we also show the outflow momentum rate ($\dot{P}_{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\times 10^5$  $\times(n_0/100 \text{ cm}^{-3})^{-0.17}$ [Kim & Ostriker, 2015],  where $n_0=100 \text{ cm}^{-3}$  is the electron density; and that the supernova rate is $0.012\times$SFR ( yr$^{-1}$) for the adopted IMF [Leitherer et al., 1999].
Since all the PUMA galaxies have $\dot{P}_{out}$ smaller than the total momentum injected by supernovae (SNe,  see dashed line in the middle panel of Fig. <ref>), their outflow momentum rate could be explained entirely by SNe.

The outflow kinetic luminosity is shown in the right panel of Figure <ref>. 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\times 10^{41}\times$SFR [yr$^{-1}$] <cit.>, 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:

\begin{equation}
\log L_{out}\ \text{[erg s}^{-1}] = (39.4\pm0.08)+(1.3\pm0.05) \log SFR\ \text{[M}_{\odot}\ \text{yr}^{-1}].
\end{equation}

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\pm0.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  yr$^{-1}$).

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_{\odot}<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 ($\eta$=/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$\sigma$ uncertainty.

§.§.§ 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 momentum-driven outflows,  the mass-loading factor ($\eta = \dot{M}_{out}/SFR$) is proportional to $v_{out}^{-1}$, while for energy driven outflows $\eta \propto v_{out}^{-2}$ <cit.>.
It is possible to distinguish between the two scenarios by looking at the slope of the relation between $\log \eta$ and $\log v_{out}$:

\begin{equation}
\log \eta =\log \left( \frac{\dot{M}_{out}}{SFR} \right) = \alpha \cdot \log v_{out} + b,
\label{eq:mass_load}
\end{equation}

where $ \alpha$ and $b$ are the slope and intercept of the linear relation.
In Figure <ref> we show $\log \eta$ as a function of $\log v_{out}$ for the PUMA and ALMA-PHANGs samples.
We fit a linear relation between these two quantities using  the Monte Carlo Markov-Chain (MCMC) implementation PyStan[https://mc-stan.org/users/interfaces/pystan].
The lightblue shaded region shows the 1 $\sigma$ 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 $\eta$. 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$\sigma$ 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 <ref>).

The best linear bisector fit gives a slope $\alpha=-1.43\pm0.28$ ($\sigma_{intr}$ = 0.35). This slope is  between the momentum-driven ($\alpha=-1$) and energy-driven ($\alpha=-2$) value.  From Fig. <ref>, we can see that both slopes are included within the 1$\sigma$ fit uncertainty (lightblue shaded area),  and therefore we cannot distinguish between the two options on the basis of this figure.

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$\sigma$ 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 ($\alpha=-2$) and the red dotted-dashed line for a momentum-driven  outflow ($\alpha=-1$).

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\pm0.01$);  the 1$\sigma$ 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. <ref>.

In  Figure <ref>, 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:

\begin{equation}
\log v_{out} \propto (0.25\pm0.01) \cdot \log SFR.
\end{equation}

The slope of the relation is in agreement with the relation found for the ionised gas by Arribas et al., 2014 (slope $0.24\pm0.05$) and for neutral gas by Rupke et al., 2005 (slope $0.21\pm0.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}\propto L_{bol}^{0.25}$.
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 \cdot \log v_{out}$ in the expression of the mass-loading factor:

\begin{equation}
\begin{split}
& \log \eta =\log \left( \frac{\dot{M}_{out}}{SFR} \right)
= \log \dot{M}_{out}- \log SFR \\
& = \log \frac{v_{out}\cdot{M}_{out}}{R_{out}}- 4\cdot \log v_{out} \\
& = \log \frac{{M}_{out}}{R_{out}}- 3\cdot \log v_{out}.
\end{split}
\end{equation}

Substituting this in equation <ref> gives:

\begin{equation}
\begin{split}
\log \eta = \log \frac{{M}_{out}}{R_{out}}- 3\cdot \log v_{out}= \alpha \cdot \log v_{out} +b
\end{split}
\end{equation}

\begin{equation}
\begin{split}
\log \frac{{M}_{out}}{R_{out}}= (\alpha+3) \cdot \log v_{out} +b.
\end{split}
\end{equation}

Thus, by looking at the relation between ${M}_{out}/R_{out}$ and $v_{out}$, we can derive the value of $\alpha$.
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 <ref> 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\pm0.25$, which implies $\alpha=-0.39$. The $\alpha=-1$ scenario (momentum-driven) is within the 1$\sigma$ uncertainty (see shaded area in Fig. <ref>), while the $\alpha=-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 $\alpha =-1.0\pm 0.5$ for a sample of five pure starburst galaxies, in agreement with the momentum-driven scenario.
Pereira-Santaella et al., 2018 find $\alpha=-0.3\pm0.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.

Outflow mass divided by outflow radius (/) versus outflow velocity.  The lightblue line is the best linear fit to the data (slope=$2.61\pm0.25$). The shaded area indicates the 1$\sigma$ 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 ($\alpha=-2$) and the red dotted-dashed line for a momentum-driven  outflow ($\alpha=-1$).  Symbols as in Fig. <ref>.

§.§.§ 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 <cit.>.
We use equation (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 $\langle M_{dyn} \rangle = 3.9\cdot10^{10}$ . The dynamical masses have been obtained from  the modelling of the rotation of the ionised gas  <cit.>.
As truncating radius we assume two times the average effective radius $\langle R_e \rangle = 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\pm40$ .  We apply an inclination correction to determine the average `observed' $v_{esc, obs}=v_{esc}/1.27= 382$ (see Sec. <ref>).

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 $\langle f_{esc} \rangle= 45\pm6\%$.
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\times\langle R_e \rangle=8$ kpc),  $v_{esc}$ would increase by 13%.
The escape outflow rates are in the range 1-173  yr$^{-1}$,  with an average $\langle \dot{M}_{esc} \rangle = 40\pm10$  yr$^{-1}$.
We find that interacting systems tend to have lower  $f_{esc}$ than mergers (mean $f_{esc}=33\pm6$% vs.  $f_{esc}=60\pm10$%).

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)/\dot{M}_{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 doublet, which is an alternative tracer used to identify molecular outflows.

We compare the OH and CO profile of each PUMA target in Figure <ref>,  <ref>,  and <ref>. 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 Figure  <ref> we show four example cases:  i) broad OH absorption reaching  high negative velocities ($v< -1000$ ), but no CO blue-shifted emission at $v < -300$ (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$ ) much larger than CO outflow maximum velocity $v_{max}=-800 $ (20100-4156);
iv) OH outflow maximum velocity ($v_{98}=-360$ ) smaller than CO outflow maximum velocity $v_{max}= -600$ (15327+2340).

Comparison of the OH 119 and CO(2-1) line profile.  The CO(2-1) spectrum is extracted from a nuclear region corresponding to the ALMA beam size  and  is convolved to the resolution of the OH spectrum (FWHM$\sim270$ ). We identified four different cases: i) OH absorption reaching high velocity ($|v|> 1000$ ),  but no CO blue-shifted emission  (00188-0856);
ii) CO outflow, but no sign of OH outflow (22491-1808);
iii) outflow detected in both OH and CO,  but OH outflow maximum velocity ($v_{98}$) larger than CO outflow maximum velocity (20100-4156);
iv) OH outflow maximum velocity ($v_{98}$) smaller than CO outflow maximum velocity (15327+2340).

In Figure <ref>, 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 (, 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$ . The |$v_{98}$(abs)(OH)| values  tend to be higher (on average by 170 ) 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)\text{(OH)}|< 300$ ,  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 Figure <ref> we compare the velocity difference $|v_{98}(abs)\text{(OH)}-v_{max}\text{(CO)}|$ with the CO outflow radius . We find a weak anti-correlation  between the  and the OH-CO  velocity difference ($r=-0.4$,  or $r=-0.3$ if we exclude the point with $v_{max}$(CO)=0 ), 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 $\sim1$ 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)\text{(OH)}\sim -1200$ , 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}\text{(OH)}|> 300$ . There is a positive correlation ($r=0.51$) between the outflow velocity difference $|v_{98}\text{(OH)}$-$v\text{(CO)}|$ and the AGN fraction (see Fig. <ref>).
Veilleux et al., 2013 find that AGN-dominated nuclei ($\alpha_{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  absorption lines, even though some of them have molecular outflows detected in the CO or HCN emission lines.