Molecular outﬂows in local galaxies: Method comparison and a role of intermittent AGN driving

We report new detections and limits from a NOEMA and ALMA CO(1-0) search for molecular outﬂows in 13 local galaxies with high far-infrared surface brightness, and combine these with local universe CO outﬂow results from the literature. The CO line ratios and spatial outﬂow structure of our targets provide some constraints on the conversion steps from observables to physical quantities such as molecular mass outﬂow rates. Where available, ratios between outﬂow emission in higher J CO transitions and in CO(1-0) are typically consistent with excitation R i 1 (cid:46) 1. However, for IRAS13120 − 5453, R 31 = 2 . 10 ± 0 . 29 indicates optically thin CO in the outﬂow. Like much of the outﬂow literature, we use α CO(1 − 0) = 0.8, and we present arguments for using C = 1 in deriving molecular mass outﬂow rates ˙ M out = C M out v out R out . We compare the two main methods for molecular outﬂow detection: CO millimeter interferometry and Herschel OH-based spectroscopic outﬂow searches. For 26 sources studied with both methods, we ﬁnd an 80% agreement in detecting v out (cid:38) 150kms − 1 outﬂows, and non-matches can be plausibly ascribed to outﬂow geometry and signal-to-noise ratio. For a published sample of 12 bright ultraluminous infrared galaxies with detailed OH-based outﬂow modeling, CO outﬂows are detected in all but one. Outﬂow masses, velocities, and sizes for these 11 sources agree well between the two methods, and modest remaining di ﬀ erences may relate to the di ﬀ erent but overlapping regions sampled by CO emission and OH absorption. Outﬂow properties correlate better with active galactic nucleus (AGN) luminosity and with bolometric luminosity than with far-infrared surface brightness. The most massive outﬂows are found for systems with current AGN activity, but signiﬁcant outﬂows in nonAGN systems must relate to star formation or to AGN activity in the recent past. We report scaling relations for the increase of outﬂow mass, rate, momentum rate, and kinetic power with bolometric luminosity. Short ﬂow times of ∼ 10 6 yr and some sources with resolved multiple outﬂow episodes support a role of intermittent driving, likely by AGNs.


Introduction
Gas outflows triggered by intense star formation and/or active galactic nuclei (AGNs) are important agents in the evolution of galaxies, in regulating the fraction of baryons inside a dark matter halo that are converted to stars, and for the metal enrichment of the circumgalactic and intergalactic medium (see, e.g., reviews by Veilleux et al. 2005;Fabian 2012). A significant step forward was recently made in characterizing the molecular phase of outflows from local galaxies via detection of P-Cygni profiles of far-infrared OH lines (e.g., Fischer et al. 2010;Sturm et al. 2011;Veilleux et al. 2013;Spoon et al. 2013;González-Alfonso et al. 2017) and detection of broad line wings to the millimeter (mm) emission lines of CO (e.g., Feruglio et al. 2010;Alatalo et al. 2011;Cicone et al. 2014;Pereira-Santaella et al. 2018;Fluetsch et al. 2019, and references therein). Molecular outflows have been detected preferentially but not exclusively in dusty (ultra)luminous infrared galaxies ((U)LIRGs). Outflow ratesṀ out up to several 100 M yr −1 and mass loading factorsṀ out /SFR around and significantly above 1 have been observed. In extreme cases and if velocities are high enough to escape the galaxy, extrapolation of the current outflow rate would clear out most of the gas content of a galaxy within 10 7 yr. While both star formation and AGNs may contribute to driving outflows, the OH outflow velocity in ULIRGs as well as the CO-based outflow rates seem to be best correlated with AGN luminosity Cicone et al. 2014;Fiore et al. 2017;Fluetsch et al. 2019). Ionized, atomic, and molecular gas phases may all significantly contribute to the outflow rates (Contursi et al. 2013;Rupke & Veilleux 2013;Janssen et al. 2016;Rupke et al. 2017), in this paper we focus on the usually dominant molecular phase. There is a continuum of properties from the fast and energetic outflows reported in the studies mentioned above, towards smaller-scale and typically slower and less energetic molecular outflows. These are increasingly reported from high spatial-and spectral-resolution ALMA and NOEMA interferometric studies of nearby galaxies (e.g., Combes et al. 2014Combes et al. , 2019Querejeta et al. 2016;Salak et al. 2016;Slater et al. 2019;Ramakrishnan et al. 2019;Alonso-Herrero et al. 2019), along with indications of inflow, warps, and noncircular motions in complex, for example barred, potentials. Here, we focus on fast and energetic outflows with potential nonlocal impact.
Thanks to efficient Herschel spectroscopy, the most statistically robust studies currently are those concerning searches for OH P-Cygni absorptions. For those, the step from measuring outflow kinematics (e.g., Veilleux et al. 2013;Stone et al. 2016) to deriving outflow masses and outflow rates involves complex modeling (e.g., Sturm et al. 2011;González-Alfonso et al. 2017;Stone et al. 2018). For the CO-based searches, the obvious and highly desirable next steps are to extend the local galaxy samples under study, and to broaden the range of galaxy parameters covered. Additional value can come from truly spatially resolved studies of very nearby outflow sources.
Three major sources of uncertainty currently affect the derivation of molecular outflow rates from interferometric CO data: (i) the emission from outflowing gas has to be separated from that of the host; (ii) a CO conversion factor has to be adopted to convert from luminosity to gas mass, but may not necessarily be the same as for typical Galactic molecular clouds; and (iii) if only an outflow velocity and overall size are known, the conversion from outflow mass to outflow rate is uncertain by at least a factor of three due to different possible geometrical assumptions (see discussions for example in Cicone et al. 2014;Veilleux et al. 2017, and Sect. 4.2). Poorly known inclination of the outflow can add to these geometrical uncertainties.
Even if separated reliably from gas orbiting in the host galaxy, the outflowing gas will face a different fate depending on its velocity and the depth of the potential. There is a continuum from "fountains" that will re-join the host in much less than a gigayear to gas staying in the circumgalactic medium for longer periods until being re-accreted, and to gas able to escape the galaxy's dark matter halo. In Sect. 3 below we conservatively identify outflowing gas with line wings only. This emphasizes the second and third paths, but at typical escape velocities of the order of 600 km s −1 the fraction of gas truly escaping a galaxy's potential is typically small if ballistic motion is assumed (e.g., Fluetsch et al. 2019). Focus on fast wings outside the line core of the galaxy also avoids ambiguities between more modest velocity outflow and other noncircular motions such as inflow or flows in the potential of a stellar bar.
Using Herschel 70 µm data, Lutz et al. (2016Lutz et al. ( , 2018 derive sizes of the far-infrared emission region for large samples of local galaxies, in order to derive scaling relations of far-infrared size and surface brightness with salient galaxy properties, and establish that AGN hosts preferentially have more compact star formation than inactive galaxies. This sample can be used for an objective selection of compact far-infrared sources ( Table 1 in Lutz et al. 2016), based on observations in the Herschel archive that are not affected by bright neighbors. The 18 sources listed in Table 1 of Lutz et al. (2016) include several for which CO outflows have been recently detected, as one might expect for their high far-infrared surface brightness Σ FIR 10 11.75 L kpc −2 . If star formation powers the far-infrared, this is well beyond thresholds in star formation rate surface density above which star formation-driven outflows from galaxy disks tend to be observed in the local universe (Heckman 2002) or at high redshift (Newman et al. 2012;Davies et al. 2019). Alternatively, if an AGN significantly contributes to the observed far-infrared emission, it must be luminous and could drive an outflow. The likelihood of hosting outflows, and the proximity of the sample (z = 0.0072−0.0676 with median 0.027) make it an obvious basis for a search for local outflow laboratories, well suited for detailed study.
We here report on a NOEMA and ALMA search for CO(1-0) outflows in galaxies from the Lutz et al. (2016) sample of compact far-infrared galaxies that lack previous deep CO data, and discuss the results in conjunction with literature results for the remainder of the sample and outside the sample's definition criteria. We adopt an Ω m = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 cosmology, redshift-independent distances from the NASA Extragalactic Database (NED), if available, for z < 0.01 galaxies, a Chabrier (2003) initial mass function (IMF), a conversion for far-infrared luminosity to star formation rate (SFR) SFR = 1.9 × 10 −10 L FIR = 40−120 µm as appropriate for the Kennicutt (1998) conversion scaled to Chabrier IMF, and L IR (8-1000 µm)/L FIR (40-120 µm) = 1.9.

NOrthern Extended Millimeter Array observations
Observations for the Dec > 10 • targets were taken with the IRAM NOrthern Extended Millimeter Array (NOEMA) on Plateau de Bure between May 2016 and February 2018 (Projects s16bh, s17ar, w17cp). The CO(1-0) line was observed in the 3 mm band, mostly using the WideX correlator (3.6 GHz bandwidth) and typically 6 or 7 antennae in configurations between C and D. These low-spatial-resolution (2-4 ) configurations were chosen to initially provide maximum detection sensitivity, and not over-resolve faint and possibly extended high-velocity outflow components. Table 1 lists the baselines used, achieved RMS sensitivity per beam for a 20 km s −1 channel, and beam size. For III Zw 035, a strong outflow was detected and already spatially resolved in the low-resolution data. This source was followed up in the higherresolution A configuration using the new PolyFix correlator and new receivers. The combined data yield the smaller beam listed for this source in Table 1.
Calibration and mapping followed standard procedures in the GILDAS 1 environment (Guilloteau & Lucas 2000). For   Table 1 of Lutz et al. (2016). In reducing the ALMA data, we adopted the same channel widths as for the NOEMA data. Self calibration was applied for IRAS F09111−1007W, NGC 4418, IRAS F17207−0014, IRAS F20551−4250 but not the other sources for which we found no noticeable improvement.

Outflow properties of the targets
First, we determine whether an outflow is present, and then measure its CO flux, velocity, and size. Outflow identification can be relatively straightforward in cases where an extremely broad outflow line component is superposed on a core of velocity width that is typical for undisturbed gas moving in the gravitational potential of a galaxy. This distinction is more difficult if the velocity range with observed "broad" CO emission does not strongly exceed plausible nonoutflow rotational velocities within the galaxy's potential. Different radial distributions of molecular gas disks rotating in the potential, misaligned disk components, bar flows, or gas moving in the potential of a merger all can lead to relatively complex total line profiles over a velocity range that is similar to local circular velocities. For outflows from spatially resolved disks in regular rotation, the kinematic separation of outflowing and nonoutflowing gas can be improved by removing the regular velocity field (e.g., Newman et al. 2012), but this approach is less suited for a sample of compact nuclei that also includes mergers. Needing two Gaussian components for a satisfactory fit of the line profile is alone not sufficient to ascribe the broader component to an outflow.
We hence continue with a brief discussion of individual targets on the basis of total CO maps and selected channel maps, spectra, positions in different velocity channels, and positionvelocity diagrams, and visualize in Figs. A.1-A.13 such basic CO(1-0) properties. As appropriate for the observed spatial extent of an outflow, we extract nuclear or aperture spectra and fit multiple Gaussians to the spectrum, which is the sum of host emission and (if present) broad outflow. Visual inspection is used to choose the number of Gaussians needed to fit the host emission. Figure 1 illustrates the adopted conventions for outflow velocity and flux. For intrinsically somewhat asymmetric outflows and little obscuration, the shift ∆v of the broad component can take positive as well as negative values, and this is indeed observed (Table 3). Unless mentioned otherwise, we adopt as outflow velocity which is computed from the shift of the centroid of the broad Gaussian line component for the outflow with respect to the systemic velocity, and the full width at a tenth of the maximum of this broad outflow component. With this definition, v out will be in the line profile wings that are dominated by the outflow.
Because outflow components at velocities that also show host CO emission are hard to quantify, we quote below in Table 3 two measured CO fluxes for the outflow: S CO (gauss) is the integral of the full Gaussian broad line component, while S CO (wings) integrates this Gaussian only over the range(s) where the outflow dominates emission (i.e., contributes more than 50% of the total fitted profile). In further analysis, we adopt the more conservative S CO (wings). This is similar to the integration of wings only (e.g., Cicone et al. 2014) or to removal of host flux on the basis of a disk model (e.g. Veilleux et al. 2017;Herrera-Camus et al. 2019), but different from adopting the full broad Gaussian flux as in several other references in the literature. For objects without outflow detection, we assign limits for the outflow CO flux: S CO,wings < 10 × rms 200 × 200, based on the rms noise (Jy beam −1 ) over a 200 km s −1 channel. Using a factor ten here rather than three considers that undetected outflows may be spatially extended or spectrally broad. The outflow radius is directly taken from the maps for a few well-resolved outflows, otherwise we derive a radius Here, |∆R| is the spatial shift between line position and continuum position and FWHM refers to the line emission. Both line position and FWHM are derived from a Gaussian spatial fit in the UV plane using a velocity range that is dominated by outflow, that is, by the broad component of the spectral multi-Gaussian fit. R out is effectively a beam-deconvolved intrinsic radius. For symmetric bipolar outflow we use average properties of the two lobes. We do not attempt to correct the measured R out or v out , or the literature results adopted in Sect. 3.2, for the effects of outflow opening angle and inclination. This emphasizes homogeneity of results with respect to more complete information that may be available for a few of the best spatially resolved and strongly collimated outflows. Observed outflow quantities or limits are listed in Table 3.  Notes. See Sect. 3 for definitions and methods to measure these quantities.

Properties of individual targets
III Zw 035 ( Fig. A.1). Extended CO emission is centered on the brighter northern lenticular-shaped component of this galaxy pair (for an optical image see, e.g., Kim et al. 2013). The second fainter optical component, offset roughly along the major axis of the brighter galaxy by ∼7 toward PA ∼200 • does not stand out in CO emission or kinematics. On top of this extended disk-like CO emission with a weak SW-NE velocity gradient that does not strictly follow the major axis of the galaxy, we detect a strong and asymmetric spatially resolved bipolar CO outflow up to velocities of ±500−600 km s −1 . The outflow has PA ∼−45 • and we adopt as R out the 4 extent of its stronger NW (redshifted) side.
The line profile in an r = 3 aperture, offset by 1 each to N and W, includes both sides of the outflow and is well fit by the sum of three Gaussians. One is narrow and the other two are broad and (on average) redshifted, with FWHM 270 km s −1 and 633 km s −1 . We assign both broad components to the outflow. The sum of the two broad Gaussians dominates the flux within this aperture at all velocities. Hence, it is not possible to follow our standard procedure and separate the outflow line wings from the line core at a velocity where the outflow starts to dominate the line profile. Instead, guided by the PV diagrams, we adopt the velocity range from −120 km s −1 to 120 km s −1 as line core and include into our "wings" CO-flux-only emission from both broad components that is outside that range. The adopted ∆v, FW 10% , and v out are the flux-weighted means of the respective values for the two broad components.
IRAS F05189−2524 (Fig. A.2). Extended CO emission centered on the nucleus shows only little velocity change on 1-2 scales, consistent with a mostly face-on view of this Type 1 AGN host. A gradient from approaching (W) to receding (E) velocities is seen in the inner 0.3 . Beyond the ∼±200 km s −1 covered by these host structures, a clear red wing reaches out to ∼500 km s −1 and traces an outflow. This emission is extended and we use an r = 1 aperture to extract its flux (green component in the multi-Gaussian fit shown top right in Fig. A.2). There may be a hint of near-nuclear blueshifted emission down to ∼−1000 km s −1 , but this never reaches >3σ, and hence is not discussed further. We adopt v out = 446 km s −1 from the Gaussian fit to the r < 1 outflow spectrum, and R out = 0.41 + 0.5 × 1.66 from a Gaussian spatial fit to the 210-510 km s −1 visibilities.
NGC 2623 (Fig. A.3). The inner 2 show a roughly E-W velocity gradient from gas rotating around the nucleus of this merger. At |v| > 300 km s −1 , line wings trace a bipolar outflow in approximately N-S (polar) direction. We adopt v out = 587 km s −1 from the fit to the r < 4 outflow spectrum and R out = 1.62 from averaging the properties of the blue and red outflows at |v| = [300, 500] km s −1 .
In addition to the bipolar outflow with its blue component shifted to the N, the [−500, −300] km s −1 channel shows emission at ∼14 from the nucleus, at PA ∼130 • . This emission peaks at 4.2σ and has total flux 0.37 Jy km s −1 . Given the relatively distinct morphology, we do not assign this emission to the main bipolar outflow but interpret it as the result of an earlier outflow episode with different orientation (see also Sect. 5.3). We list this outflow as "NGC 2623 fossil". IRAS F09111−1007W (Fig. A.4). The inner 2 show a roughly E-W velocity gradient from gas rotating around the nucleus of this galaxy, which is part of a wide pair. At |v| 200 km s −1 , line wings trace a bipolar outflow in approximately N-S (polar) direction. We adopt v out = 331 km s −1 from the fit to the nuclear outflow spectrum and R out = 0.58 from averaging the properties of the blue and red outflows at |v| = [210, 510] km s −1 .
IRAS F10173+0828 (Fig. A.5). A rotating disk with kinematic major axis at PA ∼45 • is seen out to a radius of ∼8 . The bright and compact nuclear CO emission does not obviously connect to this rotating disk and has wings out to ±300 km s −1 , beyond the ±100 km s −1 LOS velocities of the disk. While this could indicate a modest velocity outflow, we cannot firmly exclude bound motion in the near-nuclear potential and hence assign an upper limit to the emission of a faster molecular outflow.
IRAS F12224−0624 (Fig. A.6). Weak blue-and redshifted outflow components are 4σ detected in addition to bright spatially compact CO that has a velocity gradient along PA ∼−30 • . We adopt v out = 566 km s −1 from the fit to the nuclear outflow spectrum. The spatial width of the weak outflow is very poorly constrained, we adopt a fiducial R out = 0.5 .
NGC 4418 (Fig. A.7). CO emission shows a compact core that is on average redshifted by ≈15 km s −1 relative to the extended emission. The redshift in Table 1 refers to extended emission at radii of 3-8 . At first glance, the nuclear spectrum and PV diagram suggest a strong red and a weak blue wing to CO(1-0). But the compact, warm, and dense nuclear region of NGC 4418 is known to have a rich chemistry and a complex mm or submm spectrum (e.g., Sakamoto et al. 2013;Costagliola et al. 2015), meaning that sensitive ALMA data detect rare molecules in addition to CO and CN, which are detected for all sources of the sample within the frequency range of our data. Specifically, there is emission near ∼−740 km s −1 relative to CO, which we identify with emission from the nitrogen sulfide radical NS, J = 5/2−3/2 at ∼115.55 GHz rest frequency (see also line surveys such as Aladro et al. 2015;Meier et al. 2015). This includes six hyperfine components, and Λ-doubling means that NS emission of a similar strength is superposed on the red wing of CO at ∼290 km s −1 , with another six hyperfine components at ∼115.15 GHz rest. We use the "blue" complex and local thermal equilibrium line ratios for the 6 + 6 lines in the blue and red complex (M. Drosdovskaya, priv. comm.) to infer and subtract contamination by NS for the red wing of CO. The adopted intensities are for T = 150 K but the relative strength of blue and red complex varies by less than 1% between 100 and 200 K, the range of central dust temperatures suggested for NGC 4418 by Sakamoto et al. (2013). We fix relative fluxes of the NS line components and their velocities relative to CO, but vary the overall flux normalization and line width to fit the blue nitrogen sulfide complex, and then subtract both complexes from the spectrum (see Fig. A.7).
Indeed much of the wing redwards of the CO(1-0) line is contributed by nitrogen sulfide, but CO line wings out to ±300 km s −1 remain on both the blue and red sides after subtraction of that contamination. We proceed with the NSsubtracted spectrum and adopt v out = 324 km s −1 from a three host plus one outflow component multi-Gaussian fit. As geometrical size we use R out = 0.23 from data that are averaging v = [−330, −190] km s −1 in the blue wing only, because the red wing is contaminated by nitrogen sulfide. Fluetsch et al. (2019) present CO(2-1) data for NGC 4418 with a line profile that is very similar to our NS-subtracted CO(1-0) profile (see also Fig. 2 below), but adopt a simple two-Gaussian decomposition that assigns a larger fraction of the total line flux to an outflow that dominates at |v| 100 km s −1 . We maintain our more conservative decomposition but note that none of the available data exclude that such a large fraction of low-velocity emission is due to slow outflow-like motion, rather than gravitationally bound motion.
IC 860 (Fig. A.8). CO emission shows a prominent compact core with an intrinsic diameter of about 1.5 . We find no significant emission at |v| > 250 km s −1 . The |v| 200 km s −1 line wings near the nucleus are broader than the velocity of larger-scale rotation and the positions of these channels fold back to near the nucleus rather than extending the PA ∼20-30 • rotation pattern. We cannot discriminate some modest outflow motion (unable to escape the galaxy) from bound motion in the near-nuclear potential. We assign a limit of 0.58 Jy km s −1 to a fast outflow.
IRAS 13120−5453 (Fig. A.9). This merger shows prominent CO emission (see also Sliwa et al. 2017). Our data indicate a clear rotating disk extending to at least 4 , with a projected rotation velocity of ∼150 km s −1 . Channel maps near systemic velocity show a depression at the center of this disk. The central line profile is double-horned with peak velocities similar to the larger-scale rotation, and shows additional outflow wings in the blue and red. Wings suggesting outflow are also seen in HCN(4-3) ) and CO(3-2) (Fluetsch et al. 2019, see also Sect. 4.1). Spectral decomposition using four host plus two outflow Gaussian components results in the outflow Gaussians contributing roughly half of the r < 1 CO flux (15% of the total flux), if integrating over the full profile. This type of profile decomposition is backed by the need for a strong optically thin central outflow component in order to explain the relation of CO(3-2) and (1-0) emission in this source (see discussion in Sect. 4.1 below, including Figs. 2 and 3). We adopt v out = 301 km s −1 from this fit to the r = 1 outflow spectrum and R out = 0.63 from averaging the properties of the blue and red outflows at v = [−690, −270] km s −1 and v = [270, 690] km s −1 . The adopted ∆v, FW 10% and v out are the flux-weighted means of the respective values for the two broad components.
IRAS F14378−3651 (Fig. A.10). The data show CO emission extended over a radius of ∼2 , with a velocity gradient along PA ∼20 • . The nuclear line profile is well fitted with two Gaussians. While there may be a slight flux excess over this fit at |v| ∼ 250 km s −1 , we do not consider this clear evidence for an outflow and assign a flux limit of 0.43 Jy km s −1 to a faster outflow.
IRAS F17207−0014 ( Fig. A.11). The inner arcsecond of this merger shows a disk-like velocity gradient along PA ∼−50 • but larger scales do not follow this trend. The spatial and velocity structure of CO in this merger is complex, including a tidal arm that is detected in CO out to at least 15 . The line profile of the galaxy is extremely broad and reaches out to ∼±500 km s −1 . Beyond that, an outflow is detected on the blue side only, out to v −1000 km s −1 . An outflow covering a similar velocity range has been reported for CO(2-1) (García-Burillo et al. 2015, see also Sect. 4.1). Our lower-spatial-resolution data are consistent with their PA ∼160 • for the blue outflow. We adopt v out = 1116 km s −1 from the fit to the nuclear outflow spectrum and an uncertain R out = 0.45 from the properties of the blue outflow at v = [−1000, −500] km s −1 .
IRAS F20551−4250 ( Fig. A.12). The velocity field of this merger is complex with a general NE-SW velocity gradient in the inner arcsecond. A clear red outflow is detected out to v ∼ 500 km s −1 , spatially offset towards the opposite direction as the red "disk" emission. A weaker blue counter-outflow is indicated at >4σ in the contours for the [−510, −210] km s −1 range as shown in Fig. A.12, but only ∼3σ for ranges ending at −270 km s −1 . We base quantitative analysis on the red outflow only, but remain aware of the tentative blue outflow when comparing to other outflow indicators. We adopt v out = 490 km s −1 from the spectral fit at the position of the red outflow, and R out = 0.27 .
NGC 7479 (Fig. A.13). There is strong CO(1-0) emission near the center and along the bar of this galaxy. The spectrum inside r = 3 shows no evidence for fast outflow at |v| > 300 km s −1 . The centroid positions of different channel maps as well as the PV-diagram indicate a pronounced bipolar CO structure at |v| ∼ 200 km s −1 , at PA ∼30 • . We do not ascribe this structure to outflow, but follow previous works (e.g., Sempere et al. 1995;Laine et al. 1999;Baker 2000), which on the basis of shallower data ascribed CO kinematics in the inner region of NGC 7479 to gravitationally bound gas flows within the barred potential of this galaxy. We assign a limit of 0.73 Jy km s −1 to a fast outflow.

Outflow data from the literature
We have searched the literature for local galaxies with interferometric low-J CO outflow detections or limits, starting with A134, page 6 of 35  2015) and Fluetsch et al. (2019). Spectra are divided by J 2 up so that they will overlap for excitation R i j = 1, i.e., same brightness temperature. the sources from Table 1 of Lutz et al. (2016) that we did not observe, but also including galaxies outside that parameter range. As noted in the introduction, our selection is not complete towards some of the smallest-scale and least energetic outflows that were recently reported from very-high-resolution interferometric studies. These would escape detection with "outflow wing" methods similar to ours, and be of lower impact beyond the local (nuclear) region of the host. Since different references follow a range of approaches for measuring outflow properties; we use the published results to reconstruct "observed" S CO (wings), v out , R out that are homogenized to our adopted conventions and consider differences in adopted distance and CO conversion factor where applicable. Where full parameters of multiple Gaussian fits are reported in the literature, our approach as described in Sect. 3 can be directly followed for deriving S CO (wings) and v out . For references that only quote outflow fluxes for a velocity range outside the host-dominated line core (e.g., Cicone et al. 2014), we adopt those fluxes and the average velocity of the outflow bin. If only the total S CO (gauss) for a broad Gaussian component is available, we assume S CO (wings) to be one-third of S CO (gauss), where the factor of one-third is adopted from the typical value measured for our own data. Again, we do not attempt to correct for outflow opening angle and inclination. This provides a homogeneous sample but does not make full use of all the information at hand for some of the best studies cases (e.g., Barcos-Muñoz et al. 2018;Aalto et al. 2016). In a few cases we use published plotted spectra to estimate v out .
For 4C 12.50 we adopt the Dasyra et al. (2014) limit for CO outflow in emission but we add their outflow velocity from CO(3-2) absorption and their size estimate. For the single dish outflow detection of IRAS F17020+4544 and the low-resolution interferometry of IRAS 05083+7936 we adopt a fiducial 1 kpc for R out . Our adopted homogenized literature data are summarized in Table 4.
From this point on, we treat our observational results and the adopted literature values jointly and consistently, as described in Sect. 4, and refer to this as the "combined" sample. For literature sources observed in CO transitions other than CO(1-0), we adopt the brightness temperature ratios in Table 4 of Bothwell et al. (2013) for conversion to CO(1-0).

Outflow structure and derived properties
Deriving outflow masses M out , mass outflow ratesṀ out , and energy and momentum flows from the direct observed quantities (Sect. 3, Tables 3 and 4) is a multi-step procedure. We describe here our assumptions, and to what extent our data constrain the uncertainties in the conversions involved.

Outflow masses and CO conversion factor
Outflowing molecular gas may not necessarily be in selfgravitating molecular clouds for which the standard "Galactic" conversion factor from CO luminosity to molecular gas mass is A134, page 7 of 35 A&A 633, A134 (2020) Fig. 3. Brightness temperature ratio R 31 of CO(3-2) and CO(1-0) emission for the ±200 km s −1 line core of IRAS 13120−5453. Data have been matched to the same beam. The nuclear region exceeds the maximum of 1 for optically thick thermal emission.
originally derived (see, e.g., Bolatto et al. 2013a). A range of possible conversion factors have been discussed for molecular outflows, ranging from a Galactic α CO(1-0) = 4.36 M /(K km s −1 pc 2 ) to α CO(1-0) = 0.80 M /(K km s −1 pc 2 ) as often used for ULIRGs (adopted partly due to the first outflow detections being in ULIRGs), and finally still lower values for α CO(1-0) that may arise in gas with optically thin CO emission. In that case, α CO(1-0) ≈ 0.34 M /(K km s −1 pc 2 ) for an excitation temperature T ex = 30 K (Bolatto et al. 2013a), and the ratio of a rotational transition i and a lower transition j (in brightness temperature units) may exceed the value R i j 1 for thermalized optically thick gas. Ratios may reach higher values up to R 21 ≈ 2.8 and R 32 ≈ 1.3, for T ex = 30 K (e.g., Leroy et al. 2015a).
Indeed, high CO line ratios and optically thin CO have been suggested for the outflow of IC 5063 (Dasyra et al. 2016). For four sources in our sample we can compare our CO(1-0) outflow results to published spectra for higher-J CO transitions, which we have re-extracted for consistent apertures (Fig. 2). Values consistent with R i1 1 are found for three targets: García-Burillo et al. 2015), and R 31 = 1.34 ± 0.23 for NGC 4418 (CO(2-1) from Fluetsch et al. 2019). We derived these brightness temperature ratios by fitting multiple Gaussians to the higher-J transition where mean velocity and width of the outflow component(s) are fixed to the CO(1-0) derived values, and only their normalization is varied. The error estimates combine statistical errors and an assumed 10% calibration uncertainty for each line. These results do not point toward the low conversion factor of the optically thin case. This is consistent with similar ratios reported for other outflows (Weiß et al. 2005;Cicone et al. 2012Cicone et al. , 2018Feruglio et al. 2015;Zschaechner et al. 2018), with α CO(1-0) = 2.1 ± 1.2 M /(K km s −1 pc 2 ) suggested by Cicone et al. (2018) on the basis of CO and [CI] in the N6240 outflow, and with an attempt to constrain the outflow conversion factor for M 82, leading to a α CO(1-0) that is two to four times below the Galactic value (Leroy et al. 2015a).
The clear exception in our sample is IRAS 13120−5453. The CO(3-2) spectrum (see also Fluetsch et al. 2019), divided by nine to get consistent brightness temperature scale with CO(1-0) clearly has wings above the CO(1-0) spectrum (Fig. 2). We derive R 31 = 2.10 ± 0.29, above the maximum value of 1 that is possible for thermalized optically thick CO emission, and consistent with a large contribution of optically thin gas that may have a moderate T ex ∼ 30 K. This high excitation is not only observed in the line wings. Figure 3 shows the beam-matched ratio map of CO(3-2) and CO(1-0) flux (in brightness temperature units R i j ), integrated over the line core within ±200 km s −1 . Values corresponding to R 31 ∼ 0.6 are observed in the disk, but R 31 > 1 is found in the nuclear region. Combined with the general CO morphology of the source (see moment-zero map and PV diagrams in Fig. A.9), the source is best described as a combination of a rotating disk with a radius of at least 4 and normal "molecular cloud" CO properties, and a largely optically thin CO outflow that is dominating the nuclear region over the full line width.
Three of the four targets studied here and the vast majority of literature results are consistent with optically thick CO emission in the outflows. However, the clear detection of a high ratio (optically thin CO) for IRAS 13120−5453 is a strong warning that this is not true for all outflows, and that associated uncertainties on conversion factor and on outflow masses cannot be fully removed with the present data. Below, we follow much of the local universe molecular outflow literature and adopt α CO(1-0) = 0.8 M /(K km s −1 pc 2 ) but note that our CO data alone do not determine a preference for this "ULIRG" value compared to a Galactic one. For a Galactic α CO(1-0) = 4.36 M /(K km s −1 pc 2 ), recently found to be consistent also with ULIRG gas-to-dust comparisons (e.g., Tacconi et al. 2018), derived outflow masses and rates would increase more than five-fold. It is also important to bear in mind that outflows may be complex and multi-phase. For detection and current constraints from high-density tracers in the Mrk 231 and NGC 253 outflows, see Aalto et al. (2012aAalto et al. ( , 2015, and Walter et al. (2017).
In the absence of a reliable source-specific conversion factor, we also apply α CO(1-0) = 0.8 M /(K km s −1 pc 2 ) to IRAS 13120−5453, noting that the associated overestimate of outflow mass for the gas that is optically thin in CO and the underestimate due to including only line wings may partly compensate. Concerning the full sample, the good agreement that we report below with OH-based outflow masses (Sect. 5.1) is sensitive to our adopted CO conversion factor.

Outflow rates
Given an observed outflow mass M out , radius R out and velocity v out , mass-outflow rates are commonly derived by assuming a spherical or (multi-)cone geometry with constant velocity within that volume, and with the simplifying assumption of no temporal or spatial acceleration or deceleration of the flow. This is a nontrivial assumption given that a multitude of mechanisms such as ballistic deceleration in the potential of the galaxy, acceleration by radiation pressure, or interaction with an entraining nonmolecular wind can all affect velocities, but it is appropriate for interpretation of the current body of observations which is typically not highly resolved. The outflow rate is then derived via  Notes. These parameters may differ from those published in the original references. They are reconstructed from the published information in order to be consistent with our adopted separation between line core and line wings, and our definitions for v out and R out . For consistency, they are not corrected for outflow inclination, even if this is constrained in some original references. See also Sect. 3.2 for details.
The value of C depends on the adopted outflow history (see also illustration in Fig. 4). Assuming that the outflow has started at a point in the past at −t flow = −R out /v out and has continued with a constant mass outflow rate until now directly leads to C = 1. This is what we adopt, and agrees with the "time averaged thin shell" approach often used to derive outflow rates from absorption line data (e.g., Rupke et al. 2005;González-Alfonso et al. 2017). In this case, the average volume density of the outflowing gas decreases outwards ρ ∝ R −2 but local cloud or filament densities may not follow this drop.
Several recent works have derived molecular outflow rates assuming C = 3 (e.g., Cicone et al. 2014;Fiore et al. 2017). This is obtained when one assumes that the outflowing gas fills the spherical or (multi-)cone volume with constant average volume density. This corresponds to a mass outflow rate that starts with its maximum in the past at −t flow , when the material was emitted that is now reaching R out , and for whichṀ out is quoted. After the start at −t flow the mass outflow rate drops in time ∝(t/t flow ) 2 , now reaching zero.
Adopting a fixed CO conversion factor means CO light traces molecular gas mass. This is plausible if the outflow is seen as A134, page 9 of 35 A&A 633, A134 (2020)Ṁ an ensemble of flowing clouds or filaments with conversion factor reflecting local cloud conditions rather than the larger-scale average density of the flow. This is an important assumption for both cases, but in particular for the first case (C = 1) with average density decreasing outwards.
Both scenarios are strongly simplified, but for resolved outflows that are strongly inclined with respect to the line of sight and have moderate opening angle, the radial profiles can provide insights into the outflow history. Figure 4 schematically shows outflow history, radial profile of the average density, and projected surface brightness for these two scenarios and for a thin shell case. Most CO outflows in our sample and in the current literature are barely spatially resolved, preventing any related conclusions. An exception is III Zw 035 (Fig. A.1), the data for which suggest an inclined outflow of modest opening angle, lacking the brightening at large radius that is expected for the constant volume density or thin-shell cases (Fig. 4). Similar arguments can be made for the strongly inclined outflows of M 82 (Walter et al. 2002;Leroy et al. 2015a), NGC 253 (Bolatto et al. 2013b), and NGC 3256S (Sakamoto et al. 2014). For ESO 320−G030, substructure is observed, but again no brightening at large radii (Pereira-Santaella et al. 2016). We therefore proceed by computing outflow rates with C = 1 in Eq. (3).

Consistency of Herschel OH-based and interferometric CO outflow detections
Samples of dozens of galaxies each have now been studied using the two complementary techniques for detecting molecular outflows in local galaxies: P-Cygni profiles of far-infrared OH lines enabled by the Herschel mission and CO line wings enabled by wide-band receivers at mm interferometers. Comparing results from these two approaches is demanded in the first place by their basic difference. Hydroxyl absorption probes columns along the lines of sight towards the far-infrared-emitting region, inducing a need for model derivations of covering factor and radial location under the assumption of overall spherical symmetry. Carbon monoxide emission traces all emitting gas, but only to the extent that it can be kinematically and/or spatially separated from the host CO emission. Both methods invoke conversions from observables to total outflowing molecular gas mass that include explicit conversion factors -OH abundance and α CO(1-0) , respectively. We compare the two approaches in two ways. First, we use our combined CO sample and a large set of OH data from the literature Spoon et al. 2013;Falstad et al. 2015;Stone et al. 2016) to check for consistency of detections and of the measured outflow velocities. Second, we compare the A134, page 10 of 35 outflow properties that González-Alfonso et al. (2017) derive for a dozen well-characterized OH outflow sources with CO resultsinterferometric data are available for their entire sample from our observations and from the literature. Figure 5 compares CO-and OH-derived outflow velocities. v out (CO) is derived as described is Sect. 4. Our v out (CO), which is not corrected for outflow inclination, is compared with the OH absorption in the line of sight. For OH, several different conventions are in use. We adopt v 84 (abs) as defined by Veilleux et al. (2013), since it is more robust at limited S/N than the alternative v max , and better matched to our v out convention for CO which samples the outflow wings but not the maximum detected velocity. In one case where the original OH reference quotes only v max (IRAS F20100−4156, Spoon et al. 2013), we estimate v 84 (abs) ∼ 800 km s −1 from their Fig. 10.
The agreement between the two approaches is good: Of 26 sources with data from both methods, 21 show both CO outflow and absorption dominated or P-Cygni OH spectra with available v 84 (abs). The CO-and OH-based velocities are in broad agreement, but there is considerable ∼0.36 dex scatter around the 1:1 line. Correlations of outflow velocities with galaxy properties will be subject to related scatter, likely for either of the two methods in use.
Some insights can be gained also from the few galaxies for which CO and OH disagree in reported outflows. NGC 7479 shows a clear fast OH outflow with v 84 (abs) ∼ −650 km s −1 (Stone et al. 2016) but no detectable CO outflow at such velocities (Fig. A.13). The line of sight sampling the OH outflow against the far-infrared continuum must have conditions that cannot be extrapolated to a larger scale and large-coveringfactor outflow. Conversely, NGC 1068 exhibits no OH outflow (Stone et al. 2016) but molecular outflows have been suggested in CO studies. This discrepancy may be due to its moderatevelocity CO outflow occurring mostly in the disk plane (e.g., García-Burillo et al. 2014), and not affecting our line of sight. The faster outflow reported by Gallimore et al. (2016) might be closer to the line of sight, but restricted to a much smaller volume and covering only a small fraction of the far-infrared continuum. For NGC 4418, González-Alfonso et al. (2012) and Veilleux et al. (2013) report OH absorption but the absorption is slightly redshifted, suggesting molecular gas in slow, ≈100 km s −1 inflow. Inflow is also suggested for atomic gas by an inverse P-Cygni profile in [Oi] 63 µm (González-Alfonso et al. 2012). Any CO emission from this slow inflow component would be hard to separate from the line core in Fig. A.7. Finally, the modest OH outflow velocity for IRAS F17207−0014 clearly differs from the fast (but low flux density) outflow seen both in CO(1-0) (this paper) and CO(2-1) (García-Burillo et al. 2015), indicating that the OH column at these fast velocities must be low towards our line of sight.
As in the explanation of these mismatches in OH and CO outflow detection, we believe that for the comparison of OH and CO velocities of the outflows in Fig. 5 that are detected in both tracers, systematic factors and geometry are main contributors to the scatter. Similar arguments apply to the comparisons of outflow masses, velocities, and rates below in Fig. 6. Concerning CO outflow velocities, the nominal errors of typically 10% for the fitted width of the broad CO outflow component do not consider the effects of possibly non-Gaussian true profiles, nor the impact of varying flux and width of host CO emission on the splitting off of the outflow wings that are adopted in our analysis (Fig. 1). Concerning the OH v 84 (abs), Veilleux et al. (2013) estimate a typical uncertainty of 50 km s −1 that is dominated by continuum placement uncertainties. Both estimates do not consider the different responses of the CO emission and OH absorption methods to variations in alignment of complex outflow geometries with respect to the line of sight, among other sources of systematic error such as the properties of the background farinfrared source against which OH is observed.
For the vast majority (≈80%) of the objects with both OH and CO diagnostics in hand, both OH absorption and CO outflow are detected, with plausible, mainly geometrical explanations for the remaining mismatches. For a sample of obscured nuclei characterized by HCN mm emission from vibrationally excited levels, Falstad et al. (2019) discuss in more detail geometrical and evolutionary factors that might lead to absence of far-infrared OH outflow signatures. The different sampling of the outflow geometry by emission and by absorption, as well as the modest S/N of many detections certainly also contribute to the scatter around the 1:1 line in the velocity comparison of Fig. 5. This good agreement would be too optimistic if dominated by CO followup of OH outflows, with preferential reporting of detections in the literature. However, CO detections and limits for our high-far-infrared surface brightness sample, and the analysis of the González-Alfonso et al. (2017) sample in the following paragraph, argue that such effects are not dominant.
A deeper comparison of CO-based outflow masses, rates, and energetics is possible for the 12 ULIRGs for which González-Alfonso et al. (2017) were able to derive reliable spherically symmetric outflow models that are constrained by both ground-state and excited level transitions of OH, observed with Herschel-PACS. All twelve sources are part of our combined CO sample. Figure 6 compares basic outflow mass, radius, and velocity as well as derived mass, momentum, and energy outflow rates. The OH models include up to three components per source; we follow González-Alfonso et al. (2017)   Of the 12 ULIRGs, only IRAS F14378−3651 lacks a CO outflow detection and our assigned upper limit for CO-based outflow mass is within a factor of two of the OH result. In fact, Fig. A.10 indicates some excess at v ∼ 250 km s −1 above the fit with two Gaussians, but we have conservatively opted in Sect. 3 to not assign this excess at modest velocity to an outflow. Also, any weak CO flux excess at the velocities −500 to −1000 km s −1 which are reached by the OH outflow is not significant. For the remaining 11 ULIRGs, outflow masses agree very well (top left of Fig. 6), with mean log(M out,OH /M out,CO ) = 0.06 dex and dispersion 0.29 dex. This agreement is noteworthy given the independent derivation of these masses, including in particular the conversion factor α CO = 0.8 M /(K km s −1 pc 2 ) on the CO side, and an interstellar medium OH abundance X OH = 2.5 × 10 −6 on the OH side. A consistent lower limit on X OH was obtained by Stone et al. (2018). The good agreement can be regarded as support for the adopted α CO(1-0) , but see also the discussion below in the context of outflow sizes.
Outflow velocities are in excellent agreement for these bright and well-studied galaxies (top center), with mean log(v out,OH /v out,CO ) = −0.03 dex and dispersion 0.08 dex.
Larger but still modest differences are found for the outflow radii. Outflow radii from the OH models are around 0.3 kpc for all ULIRGs, while the CO observations vary more widely and detect radii up to several kiloparsecs. A mean log(R out,OH /R out,CO ) = −0.54 dex is found with a dispersion of 0.39 dex. Similar ≈0.5 dex differences are found for the outflow rates for mass, momentum, and mechanical energy (bottom panels of Fig. 6). While the OH radii are based on spherically symmetric models and the CO radii are measured, it is improbable that all differences stem from modeling uncertainties, especially when considering the observed OH absorptions in the 84 µm and/or 65 µm transitions that are radiatively excited by intense far-infrared radiation fields close to an optically thick and necessarily small far-infrared source (González-Alfonso et al. 2017). Breakup into several individually optically thick sources could enlarge the sizes but this is not supported by the usually compact structure of high-resolution (sub)mm images (e.g., Barcos-Muñoz et al. 2017).
Assuming that much of the size difference between OH and CO outflow (top right of Fig. 6) is real, one is led to a scenario where a more compact and embedded molecular outflow region traced by the OH absorptions connects to a larger and older component that can still be traced by the collisionally excited CO emission. In contrast, the covering by clumps absorbing OH will decrease with outward motion. In that case, the good agreement of outflow masses (top left of Fig. 6) will be partly fortuitous or erroneous since the two molecules trace different though overlapping regions. A further factor in this comparison arises from the definition of outflow masses. Hydroxyl-based outflow masses from González-Alfonso et al. (2017) reflect all directions of a spherically symmetric outflow. Emission line components of such a flow that are moving close to the plane of sky blend with the bright core of the line and are excluded by our definition of S CO (wings); see also Fig. 1. Toy-model calculations for a spherically symmetric case suggest S CO (wings) would measure a mass two to three times smaller than the 4π outflow, but this correction is obviously dependent on outflow velocity and the width of the line core, and A134, page 12 of 35 may be an upper limit if replacing the spherically symmetric flow with a wide-angle flow with a large line-of-sight component.
A non-unique scenario that further reduces the small tensions in Fig. 6 between OH-and CO-based results can be described as follows: (i) An outflow that is constant in rate and velocity over the flow time derived from CO data is sampled in its entirety by CO emission, and its inner half or somewhat less is sampled by the OH absorptions. The molecular outflow masses traced by CO and by OH should reflect the different fractions of the flow that are sampled -the CO-based molecular outflow mass should be larger than the OH-based one. (ii) The similar derived outflow masses (Fig. 6) must therefore be due to an additional correction. Depending on actual outflow geometry, this could be achieved by considering the fraction of an approximately spherically symmetric CO outflow that is missed by the emission line wings (see above). Alternatively, increased α CO(1−0) or increased OH abundance (decreased molecular hydrogen mass traced by OH) could help in achieving an appropriate ratio of masses traced by CO and OH. (iii) Either of the corrections described in (ii) would also bring the outflow rates in the bottom panels of Fig. 6 closer to agreement, consistent with the assumed constant outflow rate.
In summary, the comparisons provide a very successful cross-validation of the OH P-Cygni and CO interferometric methods of characterizing molecular outflows. There is no indication for a need to modify the methods and conversions adopted for CO outflows in Sects. 3 and 4, nor for modifying the OH modeling methods of González-Alfonso et al. (2017). At least for ULIRG-like systems, which constitute the bulk of current outflow samples and of the comparisons above, the comparison also supports the large molecular outflow covering factor reported in OH outflow studies. This finding is originally based on the high incidence of P-Cygni absorptions in ULIRGs (e.g., Sturm et al. 2011;Veilleux et al. 2013), and is supported by the model results of González-Alfonso et al. (2017). Current knowledge of outflow geometries and histories for large samples is however insufficient to clarify whether remaining modest differences between OH-and CO-based results can find a purely geometrical explanation. Use of the time-efficient OH spectroscopy for high-z outflow searches is encouraged, be it with ALMA (Spilker et al. 2018;Herrera-Camus et al. 2020a) or with planned future far-infrared/submm space facilities (e.g., SPICA, Roelfsema et al. 2018).

Relation of outflow properties to far-infrared surface brightness, AGN luminosity, and bolometric luminosity
Several recent papers have studied scaling relations of molecular outflow masses, mass outflow rates, momentum flows, and energy flows with star formation rate, AGN luminosity, and bolometric luminosity of the host galaxies (Cicone et al. 2014;Fiore et al. 2017;González-Alfonso et al. 2017;Fluetsch et al. 2019), and to what extent the flows may be momentum or energy driven. Here, we focus on the additional comparison with farinfrared surface brightness, and briefly revisit the relations to AGN and bolometric luminosity. Table C.1 summarizes these and related properties for the combined sample. As far as other quantities are concerned, our sample concurs with observational results from these references with modest modifications due to different assumptions (e.g., C = 1 or 3 in Eq. (3) and outflow line fluxes corresponding only to wings or to entire broad components). One should note that all these studies have overlapping samples, since they are enlarging their original samples from the literature (Sect. 3.2 for this paper), or are entirely literature based. Figure 7 shows outflow velocity, outflow molecular mass, mass outflow rate, and the ratio of outflow and host molecular gas mass as a function of far-infrared surface brightness for 39 sources in our combined sample, which have measurements or limits for far-infrared surface brightness (Lutz et al. 2016(Lutz et al. , 2018. None of these quantities shows a significant trend with far-infrared surface brightness, but see the note below on selection effects at low surface brightness. If we first focus on the high-surface-brightness end with log(Σ FIR ) 11.75 L kpc −2 , we indeed find the expected very high incidence of molecular outflows that motivated our observations. Of the 18 sources with log(Σ FIR ) ≥ 11.75 L kpc −2 and/or close to optically thick far-infrared emission that are listed in Table 1 of Lutz et al. (2016), all but ESO 173−G015 have deep CO observations probing for outflows, from our own data or from the literature in the combined sample. Of these 17 objects, CO outflows have been detected in 13 sources while upper limits have been derived for IRAS F10173, IC 860, IRAS F14378, and NGC 7479. In several cases, these limits are conservative and a contribution of an outflow to the observed moderate velocity line wings remains a possibility. As expected, galaxies with high-far-infrared surface brightness are excellent hunting grounds for molecular outflows. This is reflected in the high detection rate: at least 13 out of 18 sources. On the other hand, the cases with sensitive upper limits demonstrate that far-infrared surface brightness alone is not sufficient to produce strong and fast outflows in all cases, despite strongly exceeding SFR thresholds above which outflows from galaxy disks are prevalent (Heckman 2002;Newman et al. 2012;Davies et al. 2019).
Given its selection above an already high surface-brightness threshold and the modest sample size, it is not surprising that no trends of outflow properties with surface brightness are detected within this high-surface-brightness sample. More noteworthy is the absence of trends in Fig. 7 over the full range of the combined sample, which extends to surface brightness that is more than an order of magnitude lower. A high incidence of molecular outflow is also found among the low-surface-brightness members of this sample. However, severe and poorly known selection effects may be at work here. Much of the molecular outflow literature is in papers studying single or very few sources, and selection criteria are diverse. Upper limits may often remain unpublished. We certainly cannot exclude the abundant existence of low-surface-brightness galaxies that have only weak or no molecular outflows. An obvious speculation is that the molecular outflow detection rate is smaller in this regime than for the high-surface-brightness sample. However, unbiased studies of larger samples are needed to establish such a difference, or to detect any trends in outflow properties over a large range of farinfrared surface brightness. Figures 8 and 9 show the same outflow quantities for the combined sample, now in relation to AGN luminosity and to bolometric luminosity, as also listed in Table C.1. Consistent with previous results that are based on samples that are partly overlapping with our combined sample (Sturm et al. 2011;Veilleux et al. 2013;Spoon et al. 2013;Cicone et al. 2014;Fiore et al. 2017;González-Alfonso et al. 2017;Fluetsch et al. 2019), we find higher outflow velocities, outflow masses, and outflow rates for sources that host a luminous AGN. However, to be clearly noted in our larger sample (Fig. 8) is the presence of some quite important outflows also in sources with low values or upper limits for L AGN (e.g., III Zw 035) that were not included in early work. There is no obvious link of outflow properties to L AGN at low L AGN 10 10 L ; the observed trends are driven by the high L AGN sources only. Such increased scatter in the relation A134, page 13 of 35 A&A 633, A134 (2020) Fig. 7. Relation of outflow properties to far-infrared surface brightness. The AGN contribution to the bolometric luminosity is color-coded on a square-root scale between zero (dark blue) and one (dark red). Black symbols in the left panels, plotted at an arbitrary low value, refer to sources lacking an outflow detection (upper limits in the right panels). High far-infrared surface brightness sources from Table 1 of Lutz et al. (2016) are at log(Σ FIR ) 11.75 L kpc −2 (vertical dotted lines). Carbon monoxide observations are close to complete here (17/18), unlike the incomplete literature results for lower surface brightness.
to AGN luminosity was also noted by Fluetsch et al. (2019). A further contribution to variation in outflow properties for a given L AGN will arise from different geometries affecting the coupling of AGN wind to the larger-scale interstellar medium. In combination, these findings suggest an important though not exclusive role of driving that is linked to present AGN activity.
In fact, quite clear relations of molecular outflow mass and mass outflow rate with bolometric luminosity are observed (Fig. 9), which we quantify as: log(M out /M ) = (7.76 ± 0.09) + (0.78 ± 0.12) Here we used linmix_err (Kelly 2007) and adopted errors of 0.15 dex for the bolometric luminosities and 0.3 dex for the outflow quantities. The value of 0.15 dex is conservatively larger than uncertainties in measured L IR or optical AGN continuum, in order to also consider uncertainties in the extrapolation from these to bolometric luminosity, and possible AGN variability. The fit for outflow mass is derived including both detections and upper limits, while the fit for mass-outflow rates includes only the outflow detections since R out and v out are needed to derive outflow rate or limit. We have verified that adopting fiducial R out = 0.7 kpc and v out = 450 km s −1 to tentatively include the limits on outflow mass also for outflow rate leaves the slope of Eq. (5) effectively unchanged, and reduces the normalization by ∼0.1 dex. While these relations can be helpful in assessing the likely importance of a molecular outflow from a given object, or for estimating total outflow rates for ensembles, we caution again A134, page 14 of 35 Fig. 8. Relation of outflow properties to AGN luminosity. The AGN contribution to the bolometric luminosity is color-coded on a square-root scale between zero (dark blue) and one (dark red). Black symbols in the left panels, plotted at an arbitrary low value, refer to sources lacking an outflow detection (upper limits in the right panels).
about possible selection effects in the combined sample that is significantly literature based. In addition, bolometric luminosity, AGN luminosity, and SFR are all correlated for this sample, causing a link of outflow properties to all three, including the SFR that we do not discuss in detail here.
There are two obvious ways to reconcile the preference of strong molecular outflows for AGN hosts with the detection of outflows also in galaxies that are lacking an AGN, and with the increased scatter at low L AGN in Fig. 8. First, intense star formation is a significant contributor to driving molecular outflows, which becomes dominant in the absence of a powerful AGN. Second, AGN-driven outflows may not be recognized as such because they persist for some time even if rapid accretionrate variations lower the AGN luminosity. A third possibility is driving by highly obscured AGN activity, but we note that for many cases the AGN luminosities in Table C.1 use infrared spectroscopy which is sensitive to AGN-heated dust even for obscured cases. We believe that both star formation driving and intermittent AGN activity are at work here, and address the second option in some more detail in Sect. 5.3. First however we discuss constraints from kinetic power and momentum rate balance. Figure 10 compares the outflow kinetic power P kin = 0.5 M out v 2 out with AGN luminosity and bolometric luminosity. Most AGN hosts stay below a wind kinetic power of 5% of the AGN luminosity, which is often considered the limit for coupling AGN radiative power into outflows via "energy-conserving" winds (e.g., King & Pounds 2015). Three AGN hosts (NGC 1266, Mrk 876, IRAS F17020+4544) hover around or slightly above that line, still consistent with plausible uncertainties of wind parameters or with modest AGN variability. Notable is that limits on AGN luminosity place many of the sources with nondetected AGNs well above that line. Also, their outflow momentum rate would sometimes require "momentum boosts" A134, page 15 of 35 A&A 633, A134 (2020) Fig. 9. Relation of outflow properties to bolometric luminosity. The AGN contribution to the bolometric luminosity is color-coded on a square-root scale between zero (dark blue) and one (dark red). Black diamonds in the left panels, plotted at an arbitrary low value, refer to sources lacking an outflow detection (upper limits in the right panels). Dashed lines reflect the fitting relations from Eqs. (4) and (5).
far above the maximum AGN radiative momentum that corresponds to the limit on L AGN . This boost reaches up to more than a factor 100 (Fig. 11), beyond plausibility even for energy-conserving flows (Faucher-Giguère & Quataert 2012;Zubovas & King 2012). These two findings suggest that some of the observed molecular outflows cannot be driven by the currently observed level of AGN activity.
The right panels of Figs. 10 and 11 compare outflow power and momentum rate to bolometric luminosity and associated momentum rate, which for the nonAGN systems are both dominated by star formation. Comparing outflow kinetic power to the total mechanical luminosity delivered by supernovae P kin,SN = 7 × 10 41 (SFR/M yr −1 ) erg s −1 , a relatively high conversion efficiency 10% would be needed in some cases. Similarly, the momentum rate of some of the stronger outflows from nonAGN hosts exceeds the radiative momentum rate. While this is in principle possible via multiple scattering in heavily obscured regions (e.g., Murray et al. 2005), we consider it more likely that we are facing a mix of driving by star formation and by past AGN activity.

Evidence for a role of intermittent AGN driving
We now present arguments that the contribution of AGNs to the driving of outflows must have a short-timescale and intermittent component in order to match the observed outflow properties. Figure 12 shows that typical flow times t flow = R out /v out for the combined sample are short, just above a million years (median 1.5 Myr). For an outflow which is so massive that it needs to be driven by not only a single young cluster or A134, page 16 of 35 Fig. 10. Relation of outflow kinetic power to AGN luminosity (left) and bolometric luminosity (right). The AGN contribution to the bolometric luminosity is color-coded on a square-root scale between zero (dark blue) and one (dark red). The continuous line indicates coupling of 5% of the luminosity into outflow power, while the dotted line indicates coupling of 10% of the supernova kinetic power into outflow power (assuming the bolometric luminosity is SFR dominated). The dashed line in the right panel reflects the fit from Eq. (6). Fig. 11. Relation of outflow momentum rate to AGN luminosity momentum rate (left) and bolometric luminosity momentum rate (right). A 1:1 relation is over-plotted. The dashed line in the right panel reflects the fit from Eq. (7). The AGN contribution to the bolometric luminosity is color-coded on a square-root scale between zero (dark blue) and one (dark red). star formation region but rather by a large part of a galaxy's total (SFR+AGN) energetics and/or momentum, such a short timescale is difficult to achieve with star formation events. Given the star formation history and the smoothing effect of massive star lifetimes, star formation events would naturally lead to a 10 7 yr timescale, even in merger-induced starbursts with peaked star formation histories (e.g., Mihos & Hernquist 1996;Hopkins et al. 2008Hopkins et al. , 2013Hayward et al. 2014). In a suite of 112 hydrodynamic major merger simulations from Cox et al. (2006) and Wuyts et al. (2010), the star formation rate stays above half of its peak value for a median duration of ∼70 Myr, with a minimum of 5 Myr and only three cases with duration below 10 Myr. Short outflow events of ∼1 Myr in massive galaxies are therefore much easier to explain with the highly and rapidly variable accretion onto an AGN, even if considering that observed outflows may average to some extent over rapid repetitive AGN "flickering" (Zubovas & King 2016).
This argument assumes that the R out , v out values obtained from CO interferometry and used to derive the flow time are representative for the full duration of the outflow event, have A134, page 17 of 35 constant velocity, and start at the center of the galaxy. This assumption would be incorrect if molecular gas from earlier phases of this outflow event were to remain undetected due to observational limitations, or due to being converted to atomic or ionized gas phases. What might happen to such molecular gas? If it is found at low surface density and CO surface brightness, observations with overly small beams might miss an important fraction of the outflow. While such effects surely exist, we infer from the lack of a trend in the inset of Fig. 12 that beam size effects are not dominating the observed flow-time distribution. Also, accumulated fast molecular gas from long outflow episodes has not been reported in large-aperture single-dish data that are sensitive to low-surface-brightness CO. In a situation where a large fraction of the outflowing gas will not finally escape the halo (e.g., Fluetsch et al. 2019), flow times as derived here are also limited by deceleration. Such effects will however work on 10 7 yr timescales related to the dynamical timescale, longer than the short typical t flow observed here.
Several studies have compared molecular outflows to outflows in the atomic phase and the ionized ∼10 4 K phase and concluded that, while the mass in those phases may sometimes approach the outflowing molecular mass, it does not exceed it by large factors (Contursi et al. 2013;Rupke & Veilleux 2013;Leroy et al. 2015a;Janssen et al. 2016;Rupke et al. 2017;Fiore et al. 2017;Fluetsch et al. 2019). Specifically for M 82, Leroy et al. (2015a) suggest a molecular fraction in the outflow that is decreasing with radius, and falling below half at ∼1 kpc, but with overall similar molecular and atomic outflowing mass. Less observationally constrained is a path to very hot 10 7 K dilute gas that is not easily detected in common X-ray data (Strickland & Heckman 2009). However, only for the fastest observed molecular outflows could the post-shock temperature T = 3 16 µv 2 s k reach beyond 10 7 K, and potentially high local densities may favor re-cooling. Finally, some outflowing molecular gas could be converted to stars (Zubovas & King 2014;El-Badry et al. 2016;Wang & Loeb 2018), but observational confirmation is still scarce (Maiolino et al. 2017;Gallagher et al. 2019;Rodríguez del Pino et al. 2019), not yet having reached unambiguous detection of young stars that are both moving at the outflow speed and trapping a significant fraction of the outflow rate. Very large masses of young stars above the molecular outflow mass would build up, and should be little obscured, since the outflows typically extend beyond the dusty galaxy centers. In summary, there seems to be no obvious storage of molecular outflow in other phases that would help to drastically shorten the observed molecular gas flow times.
A second argument for time-dependent AGN driving rests on the energy and momentum budget. For about 10% of their sample, Fluetsch et al. (2019) argue that outflow kinetic power and momentum cannot be driven by the combination of the currently available AGN luminosity and SFR -a past episode of higher AGN luminosity must have been present.
Finally and most directly, in a few spatially well resolved cases, distinct outflow episodes can be inferred from the morphology. NGC 2623 shows both a clear bipolar outflow with the blue lobe north of the nucleus, and a second "fossil" blue component at a greater offset to the southeast (see Sect. 3, Fig. A.3). Flow times for these two components differ by more than an order of magnitude: 1 Myr versus 13 Myr. Rapidly evolving near-nuclear conditions of this advanced merger or driving by different nuclei (but note the single component appearance in K-band, Rossa et al. 2007) might be responsible for the different orientation of these two flows. For IRAS 08572+3915, deep IRAM-NOEMA data detect a second spatially separated outflow cloud that is offset by ∼6 kpc (Janssen 2016;Herrera-Camus et al. 2020b) in addition to the previously known outflow (Cicone et al. 2014). On a smaller scale, spatially resolved ALMA data of ESO 302−G030 (Pereira-Santaella et al. 2016, their Fig. 6) show several clumps over a total scale of ∼0.8 kpc, with a rather well defined blue-red symmetry. Driving by individual supernovae may just be energetically possible in this case, but intermittent AGN emission appears equally plausible.
Accretion rate variations on many timescales are a characteristic of the AGN phenomenon. Beyond the timescales of direct flux monitoring programs or archival light curve studies, their imprint on local AGN phenomena has been demonstrated in various ways. Light echoes of past activity can be seen in the form of (extended) narrow line regions that are too powerful for the currently observed AGN luminosity (e.g., Lintott et al. 2009;Keel et al. 2012;Schirmer et al. 2013; see also Sartori et al. 2018). Conversely, AGNs with deficient NLR compared to their X-ray luminosity may trace the onset of activity (Schawinski et al. 2015). Perhaps most directly, X-ray light echoes trace accretion variations of the (much less luminous) activity in our Galactic Center (e.g., Clavel et al. 2013). Recent simulations start to reproduce the observed "flickering" AGN accretion variations (e.g., Novak et al. 2011;Gabor & Bournaud 2013). At high redshift, the smoothing of accretion rate variations provided by the observed outflow can explain the higher incidence of AGN-driven nuclear outflows compared to direct AGN indicators such as X-rays (Genzel et al. 2014;Förster Schreiber et al. 2019).
Molecular outflow statistics in a complete mass-selected sample could provide further insights into the timescales of AGN driving. Limitations arise from the additional role of starformation-driven outflows that can be hard to discern from fossil AGN-driven outflows, and from the diverse target selection of outflow studies in the current literature. Tentatively, both the identification of multiple episodes in a few targets and the large outflow incidence (13 out of 17) for the Lutz et al. (2016) high-surface-brightness sample argue that the ∼1 Myr outflow events may be episodes in a longer outflow history. However, the latter sample is a special population that favors a compact dense circumnuclear ISM, and the low incidence of OH outflows in the Seyfert sample of Stone et al. (2016) clearly shows that its statistics cannot be transferred to AGNs in general. In any case, "depletion times" needed to expel a the entire molecular gas mass of a galaxy have to be viewed with great caution if one assumes that the current outflow rate persists unchanged over a very long time.

Conclusions
We have presented new CO(1-0) interferometry from ALMA and NOEMA, probing molecular outflows in 13 local galaxies with high far-infrared surface brightness, and combined our results with literature data on local universe molecular outflows. Our main findings are: (1) Galaxies with high far-infrared surface brightness log(Σ FIR ) 11.75 L kpc −2 show the expected high incidence of molecular outflows. Suitable CO data is now available for 17 of the 18 sources in the Lutz et al. (2016) sample, and 13 of these exhibit molecular outflows. Several of these nearby targets, in particular III Zw 035, are well-suited to spatially resolved studies.
(2) Partly constrained by observed line ratios and by spatial structure of the best cases, we adopt the following constants for the conversion of CO observables to outflow physical properties: α CO(1−0) = 0.8 M /(K km s −1 pc 2 ) and C = 1 iṅ M out = C M out v out R out . The high R 31 = 2.1 in IRAS 13120−5453 clearly indicates that this outflow is optically thin in CO. This suggests that α CO(1−0) may vary between outflows, but the good agreement with independent OH-based results globally supports the adopted CO conversion factor for outflows. We adopt conservative definitions of outflow flux and outflow velocity that focus on the line wings. While this may miss slow outflow components, which are likely to fall back, it improves robustness against confusion of outflow with undisturbed host gas.
(3) We compare our results with outflow properties that are based on Herschel far-infrared spectra of OH absorption, with special focus on the 12 ULIRGs with detailed OH models by González-Alfonso et al. (2017). The good agreement in derived properties is an important cross-validation of the two main methods for characterizing molecular outflows. Modest remaining differences may relate to different but overlapping regions sampled by the two methods: smaller for the OH absorption versus larger for CO emission. (4) Outflow properties correlate better with AGN luminosity and with bolometric luminosity than with far-infrared surface brightness. The outflows with the greatest mass and mass outflow rate prefer systems with current luminous AGN activity, suggesting an important role of AGN driving, but significant outflows in some nonAGN systems must relate to star formation or to AGN activity in the recent past. (5) We report scaling relations of outflow properties with bolometric luminosity (Eqs. (4)- (7)), but note that these might still be influenced by selection effects. (6) Short ∼10 6 yr flow times and some sources with spatially resolved multiple outflow episodes support a role of intermittent or episodic driving, likely by AGNs.   In the absence of a clear outflow signature or of a clear velocity gradient in the line core, we simply plot two orthogonal directions of which PA +45 • follows a velocity gradient on scales of several arcseconds. We conservatively adopt an upper limit for outflow in this source but cannot exclude that part of the line wings to ±300 km s −1 are caused by outflow.