# Molecular Outflows in Local ULIRGs: Energetics from Multitransition OH Analysis

E. González-Alfonso, J. Fischer, H. W. W. Spoon, K. P. Stewart, M. L. N.

Ashby, S. Veilleux, H. A. Smith, E. Sturm, D. Farrah, N. Falstad, et al.

► To cite this version:

E. González-Alfonso, J. Fischer, H. W. W. Spoon, K. P. Stewart, M. L. N. Ashby, et al. Molecular Outflows in Local ULIRGs: Energetics from Multitransition OH Analysis. The Astrophysical Journal, 2017, 836 (1), pp.11. 10.3847/1538-4357/836/1/11. [hal-02105252](https://hal.science/hal-02105252v1)

HAL Id: hal-02105252

https://hal.science/hal-02105252v1

Submitted on 28 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Molecular Outflows in Local ULIRGs: Energetics from Multitransition OH Analysis

E. González-Alfonso<sup>1,2</sup>, J. Fischer<sup>3</sup>, H. W. W. Spoon<sup>4</sup>, K. P. Stewart<sup>3</sup>, M. L. N. Ashby<sup>2</sup>, S. Veilleux<sup>5</sup>, H. A. Smith<sup>2</sup>, E. Sturm<sup>6</sup>, D. Farrah<sup>7</sup>, N. Falstad<sup>8</sup>, M. Meléndez<sup>5,9,10</sup>, J. Graciá-Carpio<sup>6</sup>, A. W. Janssen<sup>6</sup>, and V. Lebouteiller<sup>11</sup>

<sup>1</sup>Universidad de Alcalá, Departamento de Física y Matemáticas, Campus Universitario, E-28871 Alcalá de Henares, Madrid, Spain

<sup>2</sup>Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

<sup>3</sup>Naval Research Laboratory, Remote Sensing Division, 4555 Overlook Avenue SW, Washington, DC 20375, USA

<sup>4</sup>Cornell University, Cornell Center for Astrophysics and Planetary Science, Ithaca, NY 14853, USA

<sup>5</sup>Department of Astronomy and Joint Space-Science Institute, University of Maryland, College Park, MD 20742, USA

<sup>6</sup>Max-Planck-Institute for Extraterrestrial Physics (MPE), Giessenbachstraße 1, D-85748 Garching, Germany

<sup>7</sup>Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA

<sup>8</sup>Department of Earth and Space Sciences, Chalmers University of Technology, Onsala Space Observatory, Onsala, Sweden

<sup>9</sup>NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>10</sup>KBRwyle Science, Technology and Engineering Group, 1290 Hercules Avenue, Houston, TX 77058, USA

<sup>11</sup>Laboratoire AIM—CEA/Saclay, Orme des Merisiers, F-91191 Gif-sur-Yvette, France

Received 2016 September 5; revised 2016 December 7; accepted 2016 December 18; published 2017 February 3

## Abstract

We report on the energetics of molecular outflows in 14 local ultraluminous infrared galaxies (ULIRGs) that show unambiguous outflow signatures (P Cygni profiles or high-velocity absorption wings) in the far-infrared lines of OH measured with the Herschel/PACS spectrometer. All sample galaxies are gas-rich mergers at various stages of the merging process. Detection of both ground-state (at 119 and 79 µm) and one or more radiatively excited (at 65 and 84 µm) lines allows us to model the nuclear gas (≤300 pc) and the more extended components using spherically symmetric radiative transfer models. Reliable models and the corresponding energetics are found in 12 of the 14 sources. The highest molecular outflow velocities are found in buried sources, in which slower but massive expansion of the nuclear gas is also observed. With the exception of a few outliers, the outflows have momentum fluxes of (2–5) × $L_{IR}/c$ and mechanical luminosities of (0.1–0.3)% of $L_{IR}$. The moderate momentum boosts in these sources (≤3) suggest that the outflows are mostly momentum driven by the combined effects of active galactic nuclei (AGNs) and nuclear starbursts, as a result of radiation pressure, winds, and supernova remnants. In some sources (~20%), however, powerful ($10^{10.5–11} L_{\odot}$) AGN feedback and (partially) energy-conserving phases are required, with momentum boosts in the range of 3–20. These outflows appear to be stochastic, strong AGN feedback events that occur throughout the merging process. In a few sources, the outflow activity in the innermost regions has subsided in the past ~1 Myr. While OH traces the molecular outflows at subkiloparsec scales, comparison of the masses traced by OH with those previously inferred from tracers of more extended outflowing gas suggests that most mass is loaded (with loading factors of $\dot{M}/SFR = 1–10$) from the central galactic cores (a few × 100 pc), qualitatively consistent with an ongoing inside-out quenching of star formation. Outflow depletion timescales are <10<sup>8</sup> yr, shorter than the gas consumption timescales by factors of 1.1–15, and are anticorrelated with the AGN luminosity.

**Key words:** galaxies: ISM – infrared: galaxies – ISM: jets and outflows – line: formation – line: profiles – radiative transfer

## 1. Introduction

The correlations found between the masses of supermassive black holes (SMBHs) and the velocity dispersions, masses, luminosities, light concentrations, and Sérsic indices of the spheroidal components of their host galaxies (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001; Tremaine et al. 2002; Marconi & Hunt 2003; Ferrarese & Ford 2005; Graham & Driver 2007; Beifiori et al. 2012; Shankar et al. 2016) suggest a fundamental link between SMBH growth and stellar mass assembly. In addition, the color distribution of local galaxies (e.g., Strateva et al. 2001; Baldry et al. 2004; Schawinski et al. 2014), with the blue galaxies actively forming stars and red-and-dead early-type galaxies evolving passively, suggests that the color of red early-type galaxies must have evolved rapidly, with star formation terminated on short timescales (e.g., Hopkins et al. 2006b; Schawinski et al. 2014). Spatially resolved observations of z ~ 2.2 massive galaxies show an inside-out quenching of star formation, on timescales of <1 Gyr in the inner regions

(Tacchella et al. 2015). An appealing way to explain these observations is via a self-regulated feedback model involving the morphological transformation of late-type to early-type galaxies through mergers, which first funnel large amounts of gas into the circumnuclear regions<sup>12</sup> of the system, leading to both a circumnuclear starburst and the growth of an SMBH. Above an SMBH critical mass, the energy or momentum released by the SMBH limits efficiently the accretion onto the SMBH and quenches the starbursts through the expulsion of the interstellar gas from which stars are formed (negative feedback), ultimately yielding the SMBH–σ relationship (Silk & Rees 1998; Fabian 1999, 2012; King 2003, 2005; di Matteo et al. 2005; Murray et al. 2005; Springel et al. 2005; Hopkins et al. 2006a; King & Pounds 2015). Violent relaxation of the stellar component deeply changes the morphology from premerger disk galaxies to a coalescence system with a

<footnote n="12"><sup>12</sup> We use the terms “nuclear region” and “nuclear starburst” for spatial scales ≲300 pc, to differentiate from the term “circumnuclear” that is more generally used for scales ≲1–2 kpc.</footnote>
spheroidal component, which dominates over the reformed disk in the case of major mergers with limited gas fraction (e.g., Hopkins et al. 2009). Observational evidence for the quenching of star formation by active galactic nucleus (AGN) feedback has been reported recently (e.g., Farrah et al. 2012; Alatalo et al. 2015). There are alternative explanations of the $M_{BH}-M_{bulge}$ correlation, e.g., gas accretion onto the black hole (BH) from a viscous inner disk with limited replenishment by star formation in an outer disk (Burkert & Silk 2001), gravitational collapse of the inner regions of an isothermal bulge (Adams et al. 2001), and hierarchical assembly of BH and stellar mass through cycles of galaxy merging (Peng 2007; Jahnke & Macciò 2011). Other ways to account for the bimodality of galaxy colors and the quenching of star formation include suppression of cold inflows of gas (Dekel & Birnboim 2006) and gravitational heating of the intracluster medium in mergers (Khochfar & Ostriker 2008). While these environment mechanisms of quenching, as well as “strangulation” on long timescales (Peng et al. 2015), may dominate the secular evolution of late-type galaxies, a much more rapid morphological and quenching ($t_{quench} \lesssim 250$ Myr) evolution is favored for producing early-type galaxies (Schawinski et al. 2014).

While feedback processes such as superwinds in ultraluminous infrared galaxies (ULIRGs) and starbursts have been observed for decades in lines of ionized and neutral atomic gas (e.g., Heckman et al. 1990; Rupke et al. 2002, 2005a, 2005b, 2005c; Lípari et al. 2005, 2009; Veilleux et al. 2005; Spoon & Holt 2009), outflows observed in lines that trace the molecular medium are also key to understanding and quantifying this process, because the molecular phase may carry a significant or even dominant fraction of the momentum and mass outflow rates in buried stages. Far-IR spectroscopy with Herschel Space Observatory Photoconductor Array Camera and Spectrometer (PACS) (Pilbratt et al. 2010; Poglitsch et al. 2010) has indeed revealed powerful molecular outflows in ULIRGs traced by OH, with velocities exceeding $1000 \text{ km s}^{-1}$ in some sources and mass outflow rates of several hundred solar masses per year (Fischer et al. 2010; Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; González-Alfonso et al. 2014b, hereafter GA14). The high-velocity outflows discovered with Herschel were found to be ubiquitous and thus inferred to be wide angle in local ULIRGs (Veilleux et al. 2013; Stone et al. 2016). These investigations also revealed a correlation between the outflow velocity and the AGN luminosity. In the far-IR, the outflows are also traced by the line wing emission of the [C II] 158 µm transition (Janssen et al. 2016). High-velocity molecular outflows are detected at (sub)millimeter wavelengths in lines of CO, HCN, and HCO$^+$ (e.g., Feruglio et al. 2010, 2015; Aalto et al. 2012, 2015; Cicone et al. 2012, 2014; García-Burillo et al. 2015; Lindberg et al. 2016). Lower-velocity molecular outflows are also detected in millimeter lines of the above species and in CS (Sakamoto et al. 2009; Bolatto et al. 2013; Alatalo et al. 2015; Tunnard et al. 2015; Martín et al. 2016). The (sub)millimeter lines are now routinely observed with high angular resolution and in some sources trace the outflowing gas out to kiloparsec scales.

Far-IR molecular observations can provide key and unique insights into the outflow phenomenon: (i) the strength and optical depth of the far-IR continuum generate P Cygni line profiles in some lines, unambiguously indicating the presence of outflowing gas, discarding other alternatives such as high turbulence or noncircular rotation motions (e.g., Guillard

et al. 2015; Díaz-Santos et al. 2016); (ii) blueshifted absorption can be traced to low velocities, probing low-velocity outflows that may be missed from pure emission lines, due to confusion with the line core; (iii) despite the relatively poor spatial resolution of far-IR telescopes, multitransition observations including high-lying transitions provide a robust means to quantify the main outflow parameters (mass outflow rate, momentum flux, etc.).

OH in particular is an excellent tracer of these molecular outflows (Fischer et al. 2010; Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013, GA14), with high abundances in active regions, including photodissociated regions (PDRs), cosmic-ray-dominated regions (CRDRs), and X-ray-dominated regions (XDRs) (e.g., Goicoechea & Cernicharo 2002; Goicoechea et al. 2011; Meijerink et al. 2011; González-Alfonso et al. 2012, 2013). Owing to the spin-orbit interaction, its rotational level structure is characterized by two ladders ($^2\Pi_{3/2}$ and $^2\Pi_{1/2}$), with the cross-ladder transitions much more optically thin than intra-ladder ones, thus enabling the estimation of column densities. With high Einstein A-coefficients, the high-lying transitions along the $^2\Pi_{3/2}$ ladder, observed in absorption in galaxies, require strong far-IR radiation densities to be excited. Specifically, the high-lying OH doublet at 65 µm is an excellent tracer of the short-lived but most active and buried phase of (U)LIRG evolution (González-Alfonso et al. 2015, hereafter GA15), tracing the gas that lost its angular momentum during the merger, falling onto the ~100 pc nuclear region where it generates a nuclear starburst (Hopkins et al. 2009). It is just in such an obscured phase when most BH accretion is expected to occur, as indicated by the cosmic X-ray background (Fabian & Iwasawa 1999), and when outflows are expected to be most efficiently accelerated (Stern et al. 2016).

## 2. Guide to the Astrophysical Issues Addressed in This Paper

Here we present the results of our study of molecular outflows in a sample of local ULIRGs observed in OH with Herschel/PACS. Our multitransition analysis aims (i) to infer and describe the presence of outflowing gas in high-lying OH transitions, which probe the nuclear regions of ULIRGs that are first exposed to the feedback action of an AGN and a nuclear starburst, and compare and relate the line profiles of these high-lying lines with those of the ground-state lines; (ii) to analyze quantitatively the OH observations to infer the main physical parameters and energetics of the molecular outflows (outflowing mass, mass outflow rate, momentum flux, and energy flux), comparing our results with those obtained from other tracers and establishing depletion timescales; and (iii) to give insight into the accelerating mechanism of these outflows by assessing the relative roles of an AGN and a nuclear starburst, the relative contribution to the acceleration by radiation pressure on dust grains and winds, and whether the observed outflows are found in momentum- or energy-conserving phases.

The observations, data treatment, and sample selection criteria are described in Section 3, where we define a modeling sample of 14 ULIRGs observed in at least three far-infrared OH doublets. While all of them are gas-rich mergers or interacting systems, they show a high morphological diversity in the merging stage. The sample spans AGN luminosities of $(0.3-2) \times 10^{12} L_{\odot}$ and star-forming rates of 50–350 $M_{\odot}$ yr$^{-1}$.
In Section 4 we present the observational results of our multitransition profile analysis. We identify the presence of gas outflowing from the nuclear regions (<300 pc) traced by radiatively excited, high-lying OH transitions. We show that the highest molecular outflow velocities, traced by the OH doublet at 119 µm (Spoon et al. 2013; Veilleux et al. 2013), are found in buried sources that show strong absorption in the OH 65 µm transition at systemic velocities. The highly excited line profiles also show evidence for slower expansion of the nuclear regions of the ULIRGs that show P Cygni profiles in the ground-state OH doublet at 119 µm.

In Section 5 we describe our library of spherically symmetric radiative transfer models, our methods of achieving best-fit solutions to the line profiles, and our method of estimating the energetics of the outflows (outflowing mass, mass outflow rate, momentum flux, and energy flux).

The results of our quantitative analysis are presented in Section 6. We compare our results with the maximum estimated momentum and energy rates that the AGN and the starburst can provide. We find that while nuclear starbursts can provide a non-negligible contribution to the observed outflows, they are most likely unable to drive them alone. An AGN is required in most cases and clearly dominates in some ULIRGs. Outflow depletion timescales are <10⁸ yr, significantly shorter than the gas consumption timescales, and show an anticorrelation with $L_{AGN}$. We also find that the nuclear outflowing activity has recently subsided in a few sources of the sample. The outflowing masses inferred from OH, which probes the outflows at subkiloparsec scales, are similar to those obtained from other tracers of more extended outflowing gas.

The inferred outflow energetics are further interpreted in terms of simple dynamical models in Section 7. An important question discussed in the literature, and related to the growth of the SMBH and to the normalization of the $M_{BH}-\sigma$ correlation, is whether the outflows are momentum-conserving (e.g., driven by the ram pressure of winds and radiation pressure on dust grains) or energy-conserving (e.g., driven by the thermal pressure of a hot bubble that cannot cool on the timescale of its expansion). We find that the total momentum deposited into the interstellar medium (ISM) by the combined effects of the AGN and starburst is enough to explain the outflows in most sources, but significant momentum boosts are apparently required in some (~20%) ULIRGs. The inferred momentum boosts of the most powerful outflows depend on whether the gravitational potential well is assumed to be balanced by rotation and on the assumed geometrical dilution of the outflowing gas, and they are estimated to be in the range of 3–20. These are, in most cases, lower than the maximum momentum boosts predicted by theoretical studies on energy-conserving outflows. Partially energy-conserving phases, which we find associated with compact outflow components with high columns, appear to represent short stages with high AGN luminosities and are uncorrelated with the merger stage. Our analysis indicates that although radiation pressure on dust grains may (nearly) support the gas in the direction of the rotation axis, it cannot drive the outflows, and thus winds are required.

Our main conclusions are summarized in Section 8. We describe in more detail the individual sources and model fits in Appendix A. A comparison of the present results with those obtained from ionized lines in starburst galaxies (Heckman et al. 2015) is given in Appendix B. Radiation pressure support

is evaluated in Appendix C, and some discussion of the velocity fields in the outflows is given in Appendix D.

# 3. Observations and Data Analysis

## 3.1. *Selection of the ULIRG Outflow Modeling Sample*

We based the selection of our outflow modeling sample on all Herschel/PACS OH observations of local ULIRGs for which at least these three doublets were available: the $OH^2\Pi_{3/2}J = 5/2-3/2$ doublet at 119 µm (hereafter OH119), the cross-ladder $^2\Pi_{1/2} - ^2\Pi_{3/2}J = 1/2-3/2$ doublet at 79 µm (OH79), and the $OH^2\Pi_{3/2}J = 9/2-7/2$ doublet at 65 µm (OH65). In most sources, the $^2\Pi_{3/2}J = 7/2-5/2$ doublet at 84 µm (OH84) has also been observed. An energy level diagram of OH showing these transitions can be found in GA14. The OH119 and OH79 lines are ground-state transitions, with OH119 40× more optically thick than OH79 (Fischer et al. 2010). The OH84 and OH65 transitions have excited lower levels with energies of ≈120 and ≈300 K above the ground state, respectively.

Five Herschel observing programs included spectroscopic velocity-resolved observations of ULIRGs in the above rotational transitions of OH: the Herschel guaranteed time key program SHINING (PI: E. Sturm), the open time program HerMoLIRG (PI: E. González-Alfonso), a program that provided the full far-infrared spectra of two (U)LIRGs (PI: J. Fischer), the HERUS program (PI: D. Farrah; Farrah et al. 2013), and a directors discretionary program focused on the two outflow sources of the HERUS sample with the most prominent OH119 absorption at high velocities (PI: H. Spoon; Spoon et al. 2013). The first three programs observed three or more OH lines in all 20 ULIRGs in the Revised Bright Galaxy Sample (RBGS; Sanders et al. 2003). The last two programs observed the ground-state OH119 and OH79 transitions in an additional complete set of 24 more distant ULIRGs out to $z < 0.2$ and followed up with profiles of the excited OH84 and OH65 transitions in IRAS 03158+4227 and IRAS 20100–4156. In all, 22 ULIRGs at $z < 0.2$ were observed in three or more OH transitions.

To define the ULIRG outflow modeling sample discussed in Sections 5–7, we further constrained the sample described above to the ULIRGs in which (i) a P Cygni profile or a high-velocity blue wing in OH119 is detected (Veilleux et al. 2013; Spoon et al. 2013) and (ii) the excited OH84 and/or OH65 doublets are detected. Of all ULIRGs reported in the OH65 doublet (GA15), the former condition rules out the following sources: Arp 220, IRAS 15250+3609, IRAS F17207–0014, and IRAS F22491–1808, where the OH119 doublet is dominated by foreground absorption at redshifted velocities, and IRAS 07251–0248, where no OH119 redshifted emission feature is seen and the very strong and flat absorption feature does not clearly indicate outflowing gas. Modeling these sources requires a specific approach with more than the three-component scheme used in Section 5. We also do not model IRAS 12112+0305 and IRAS 19542+1110, both showing a P Cygni profile in OH119 (Veilleux et al. 2013), because of a mispointing in the three observed doublets (OH119, OH79, and OH65), resulting in the placement of the source near the edge of a spaxel, which causes a skewing of the instrumental profile. This results in a modeling sample of 14 local ULIRGs, which are listed in Table 1, together with some basic properties of the sources and the observation identification numbers (OBSIDs)
<table><caption>**Table 1**
Modeling Sample Galaxies and Herschel OBSIDs of the OH Doublets</caption><thead><tr><th>Galaxy Name<br>(1)</th><th>z<sub>CII</sub><br>(2)</th><th>D<br>(Mpc)<br>(3)</th><th>L<sub>IR</sub><br>(10<sup>12</sup> L<sub>☉</sub>)<br>(4)</th><th>α<sub>AGN</sub><br>(5)</th><th>OH119<br>OBSID<br>(6)</th><th>OH79<br>OBSID<br>(7)</th><th>OH84<br>OBSID<br>(8)</th><th>OH65<br>OBSID<br>(9)</th></tr></thead><tbody><tr><td>IRAS F03158+4227</td><td>0.13459</td><td>658</td><td>4.27</td><td>0.47</td><td>1342238963</td><td>1342263478</td><td>1342262940</td><td>1342263479</td></tr><tr><td>IRAS F05189-2524</td><td>0.04272</td><td>186</td><td>1.38</td><td>0.72</td><td>1342219441</td><td>1342219442</td><td>1342248556</td><td>1342219445</td></tr><tr><td>IRAS F08572+3915</td><td>0.05824</td><td>261</td><td>1.32</td><td>0.70</td><td>1342208956</td><td>1342184687<sup>a</sup></td><td>1342253600</td><td>1342208954</td></tr><tr><td>IRAS 09022-3615</td><td>0.05963</td><td>268</td><td>1.92</td><td>0.55</td><td>1342209402</td><td>1342209403</td><td>...</td><td>1342209406</td></tr><tr><td>IRAS F10565+2448</td><td>0.04309</td><td>193</td><td>1.14</td><td>0.47</td><td>1342207787</td><td>1342207788</td><td>1342254243</td><td>1342207790</td></tr><tr><td>Mrk 231</td><td>0.04218</td><td>186</td><td>3.37</td><td>0.67</td><td>1342186811</td><td>1342186811</td><td>1342253537</td><td>1342207782</td></tr><tr><td>IRAS 13120-5453</td><td>0.03107</td><td>136</td><td>1.86</td><td>0.33</td><td>1342214628</td><td>1342214629</td><td>1342248348</td><td>1342214630</td></tr><tr><td>Mrk 273</td><td>0.03780</td><td>166</td><td>1.45</td><td>0.34</td><td>1342207801</td><td>1342207802</td><td>1342257293</td><td>1342207803</td></tr><tr><td>IRAS F14348-1447</td><td>0.08257</td><td>376</td><td>2.09</td><td>0.17</td><td>1342224243</td><td>1342224242</td><td>...</td><td>1342224244</td></tr><tr><td>IRAS F14378-3651</td><td>0.06812</td><td>304</td><td>1.46</td><td>0.21</td><td>1342204337</td><td>1342204338</td><td>1342250130</td><td>1342204339</td></tr><tr><td>IRAS F19297-0406</td><td>0.08558</td><td>383</td><td>2.46</td><td>0.23</td><td>1342208890</td><td>1342208891</td><td>...</td><td>1342208893</td></tr><tr><td>IRAS F20100-4156</td><td>0.12971</td><td>632</td><td>4.68</td><td>0.27</td><td>1342216371</td><td>1342216371</td><td>1342267942</td><td>1342268108</td></tr><tr><td>IRAS F20551-4250</td><td>0.04295</td><td>185</td><td>1.02</td><td>0.57</td><td>1342208933</td><td>1342208934</td><td>1342253748</td><td>1342208936</td></tr><tr><td>IRAS F23365+3604</td><td>0.06449</td><td>281</td><td>1.43</td><td>0.45</td><td>1342212514</td><td>1342212515</td><td>1342257685</td><td>1342212517</td></tr></tbody></table>

**Note.** Columns: (1) Galaxy name. (2) Redshifts inferred from Gaussian fits to the [C II] 158 µm line. (3) Distance to the galaxy; adopting a flat universe with $H_0 = 71$ km s⁻¹ Mpc⁻¹ and $\Omega_M = 0.27$. (4) IR luminosity (8–1000 µm), estimated using the fluxes in the four IRAS bands (Sanders et al. 2003; Surace et al. 2004). (5) Estimated AGN contribution to the bolometric luminosity, as derived from the f15/f30 ratio (Veilleux et al. 2009; Spoon et al. 2013). This method yields $\alpha_{AGN} = 0.8$ in Mrk 231, but we have reduced the value in this source to 2/3 because the nuclear starburst has an estimated contribution to the bolometric luminosity of 25%–40% (Davies et al. 2004). Our adopted value, while basically consistent with the f15/f30 diagnostic, enables easy comparison with the works by Cicone et al. (2014) and Feruglio et al. (2015), who adopted just half that value ($\alpha_{AGN} \approx 1/3$, based on Nardini et al. 2010). We have adopted $L_{bol} = 1.15L_{IR}$ (Veilleux et al. 2009, 2013). The AGN and starburst luminosities are $L_{AGN} = \alpha_{AGN} L_{bol}$ and $L_* = (1 - \alpha_{AGN}) L_{bol}$; (5)–(8) OBSIDs of OH119, OH79, OH84, and OH65.

<sup>a</sup> The line was also targeted in instrument verification observations with OBSIDs 1342184688, 1342184689, 1342184690, 1342184691, and 1342184692, and the averaged spectrum is used.

of the OH observations. The excluded sources have similar luminosities to the modeled ones, but some have lower AGN contributions that may be lower limits due to high far-IR extinction (Arp 220, IRAS 07251–0248, and IRAS F17207–0014; Veilleux et al. 2013).

The modeling sample, individually described in Appendix A, is somewhat biased to the most prominent outflowing sources, but may still be considered a good representation of the diversity of the molecular outflow phenomenon in local ULIRGs. All sample sources are morphologically classified in the optical as mergers or interacting systems, albeit with a high diversity in evolving stages (Appendix A): from widely separated galaxies (~50 kpc, IRAS 03158+4227) to double nuclei systems with projected separation <10 kpc (IRAS 08572+3915, Mrk 273, IRAS 14348–1447, IRAS 20100–4156), multiple colliding systems (IRAS 10565+2448, IRAS 19297–0406), and advanced mergers showing a single nucleus with tidal tails (IRAS 05189–2524, IRAS 09022–3615, Mrk 231, IRAS 13120–5453, IRAS 14378–3651, IRAS 20551–4250, IRAS 23365+3604). They also have a diversity in the AGN contribution to the luminosity (Table 1), resulting in estimated AGN luminosities in the range of (0.3–2) × 10¹² L☉ (due to high extinction at far-IR wavelengths (GA15), some values of $L_{AGN}$ could be underestimated). Assuming Eddington luminosities, the implied SMBH masses are 10⁷–10⁸ M☉.

## 3.2. Data Analysis

For the basic analysis developed in Section 4 (Sections 4.3–4.8), we have used the OH observations of these ULIRGs together with those of luminous infrared galaxies (LIRGs, $L_{IR} > 10^{11}$ L☉), as well as the nearby bright infrared galaxy (BIRG) NGC 4945, also detected in excited OH transitions (GA15). In all these sources, we also used additional

observations of the [C II] 158 µm, [O I] 63 µm, and [O I] 145 µm lines from the programs listed above (Graciá-Carpio et al. 2011; Farrah et al. 2013).

All lines were observed in the highly sampled range mode of the PACS spectrometer, with the exception of NGC 4945, which was observed in line spectroscopy mode with Nyquist sampling. Most of our analysis was based on pipeline-processed spectra created by the Herschel Science Centre (HSC) using the Herschel Data Processing system; the Standard Product Generation version was HIPE 14.0.1, with calibration tree version 72. We downloaded the level-2 data products<fnref n="13" /> as needed from the Herschel Science Archive (HSA). Because both the molecular absorption lines and the continua are basically point-like in the sources studied in this work, we have used the point-source-calibrated spectra “c129,” which scale the emission from the central ≈9” × 9” spatial pixel to the total emission from the central 3 × 3 spaxels (“c9”), which is itself scaled according to the point-source correction. The absolute flux scale is robust to potential pointing jitter, with continuum flux reproducibility of ±15%. The PACS spectral resolution is 290, 160, 145, and 190 km s⁻¹ at 119, 79, 84, and 65 µm, respectively.

In some specific cases, the HSA spectra showed significant fluctuations in the continuum level close to or within the wavelength range where the line wings could be expected. We reprocessed those data by applying polynomial fits to the flat-fielding with masking windows around the observed absorption and emission features, which usually significantly improved the quality of the baselines. Nevertheless, the uncertainty in the equivalent width and flux in the line wings is usually dominated by the fitted baseline and in some cases may be up to ~25%.

<footnote n="13">¹³ http://herschel.esac.esa.int/Docs/PACS/html/pacs_om.html</footnote>
The observation of the OH119 doublet in IRAS 03158+4227 was mispointed by about half a spaxel (Spoon et al. 2013). We have examined spaxels 12 (central) and 13, where the OH119 absorption/emission is detected, and found a similar wing with velocities exceeding ~1500 km s⁻¹ in both spaxels, though the actual value is relatively uncertain, due to insufficient baseline on the blue side of the spectra. The redshifted emission feature is, however, slightly stronger in the central spaxel, so the central spaxel profile is used here for subsequent analysis.

Baseline fitting was subsequently performed with polynomials of order 1–2, and the fits were used to generate continuum-normalized spectra and to extract the far-IR continuum flux densities. Gaussians were fitted to absorption and emission features to obtain the peak velocities and line strengths, as discussed in Section 4.

# 4. Results

## 4.1. Description of the OH Spectra

Figure 1 shows the continuum-normalized OH spectra in all 14 sources of our modeling sample. The OH119 and some OH79 and higher-excitation OH profiles were previously presented by Fischer et al. (2010), Veilleux et al. (2013), Spoon et al. (2013), Sturm et al. (2011), and GA14 based on earlier HIPE pipeline reductions. The velocity scale in Figure 1 uses the redshifts obtained from Gaussian fits to the high signal-to-noise ratio [C II] 158 µm line, which is strongly dominated by the core of the line (Janssen et al. 2016). In all sources, the excited [O I] 145 µm line was also observed with Herschel/PACS, giving velocity shifts relative to the [C II] line of ≲40 km s⁻¹, and in most (nine) sources ≲20 km s⁻¹.

The blue arrows in Figure 1 indicate the location of possible contamination by lines of H₂O, dozens of which appear in the far-IR spectra of some (U)LIRGs (e.g., Falstad et al. 2015; González-Alfonso et al. 2012, hereafter GA12). The strongest of these lines is the $4_{23}-3_{12}$ ($E_{\text{lower}} = 250$ K) centered at ≈−1400 km s⁻¹ in the OH79 spectra. The line is most likely contaminating the OH79 blueshifted wing in IRAS 03158+4227 at velocities ≲−1400 km s⁻¹. It is clearly detected in other sources in which the OH79 wing does not blend with it. Potentially more problematic are instances of possible contamination of OH79 spectra by the H₂O $6_{15}-5_{24}$ at ≈−720 km s⁻¹. However, this line is much higher in energy ($E_{\text{lower}} \approx 600$ K), and it is a factor of 2–7 weaker than the lower-energy line in the extremely buried sources NGC 4418 and Arp 220, respectively (GA12). It is likely that this line significantly contaminates the OH79 spectra of IRAS 14348−1447 and IRAS 20100−4156. In some cases, the OH65 doublet may be slightly contaminated by the H₂O $6_{25}-5_{14}$ line ($E_{\text{lower}} \approx 580$ K) at +160 km s⁻¹, which could have the effect of shifting the velocity of peak absorption of the OH65 blue component to less negative velocities. There are other possible contaminating lines as well. In the OH119 profiles of IRAS 03158+4227, Mrk 231, Mrk 273, IRAS 14348−1447, IRAS 20100−4156, and IRAS 23365+3604, there is substantial absorption by CH⁺ $J = 3-2$ and $^{18}$OH at velocities >1000 km s⁻¹ (>500 km s⁻¹ relative to the red component of the OH119 doublet), which may significantly weaken the emission feature of the OH119 doublet.

The high-velocity absorption wings or P Cygni profiles observed in OH119 unambiguously indicate outflowing gas in

all sample galaxies. The maximum blueshifted velocities observed in OH119 range from the very moderate ~250 km s⁻¹ in IRAS 09022–3615 to ~1800 km s⁻¹ in IRAS 03158+4227, while most sources show maximum velocities in the range of 600–1300 km s⁻¹ (see Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013). With the exception of IRAS 19297–0406, all sources in Figure 1 also show P Cygni profiles or high-velocity blueshifted wings in the cross-ladder OH79 transition. Relative to the absorption strength, the redshifted emission feature is usually stronger in OH79 than in OH119. This is most likely because OH79 emits both through direct absorption and reemission of 79 µm continuum photons and via absorption of 53 and 35 µm photons and cascades down to the ground state through the Π₁/₂ $J = 1/2$ level (GA12), while the OH119 doublet only efficiently scatters through direct absorption and reemission of 119 µm continuum photons.

In contrast, P Cygni profiles in the excited OH84 and OH65 doublets are only observed in IRAS 03158+4227 (OH84), IRAS 08572+3915 (OH84 and OH65), and IRAS 23365+3604 (OH84). This can be attributed to the fact that the excited lines are formed closer to the optically thick far-IR source responsible for their excitation (GA14), and thus their emission from the far side is obscured. Nevertheless, OH84 blueshifted line wings are observed up to velocities of ~1000 km s⁻¹ in IRAS 03158+4227, IRAS 08572+3915, Mrk 231, IRAS 14378–3651, and IRAS 20100–4156, while significant wings with velocities <400 km s⁻¹ are also observed in IRAS 05189–2524, IRAS 13120–5453, Mrk 273, IRAS 20551–4250, and IRAS 23365+3604. The OH65 transition requires high columns to be excited (GA15); thus, high-velocity absorption wings in this line are only found in the most extreme sources: IRAS 03158+4227, IRAS 08572+3915, Mrk 231, and IRAS 20100–4156. We argue in Section 5 that a high-velocity wing in OH65 is the most reliable indicator of powerful AGN feedback.

It is also worth noting that in most sources in our sample, the main absorption features of the OH84 and OH65 doublet components, while peaking at around central velocities, are slightly blueshifted relative to the systemic velocities as measured in the [C II] 158 µm line. This systematic effect is further explored in Section 4.5, where we conclude that high columns of gas (if not most of the gas) are slowly expanding from the nuclear regions of ULIRGs, thus potentially changing the morphology of these regions and shifting the star formation to increasing radii.

## 4.2. Comparing Outflow Properties as Traced by Both OH and CO

It is instructive to compare some observational parameters of molecular outflows as derived independently from OH and CO observations. Studies of outflows in local ULIRGs based on CO observation have been carried out for five of our sources (Feruglio et al. 2010; Cicone et al. 2012, 2014; Feruglio et al. 2015): Mrk 231, Mrk 273, IRAS 08572+3915, IRAS 23365+3604, and IRAS 10565+2448. More recently, CO observations have also been obtained in two additional OH outflow sources, IRAS 03158+4227 and IRAS 20100–4156 (A. Gowardhan et al. 2017, in preparation). Given that the CO luminosity in the line wings is expected to be proportional to the gas mass of the molecular outflow, comparison between OH and CO sheds light on the physical properties associated
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Continuum-normalized spectra of OH in the 14 sources that are modeled in Section 5. The abscissa shows velocities referring to the blue component of each doublet, and the two vertical gray lines indicate the rest position of the two doublet components. Redshifts are obtained from Gaussian fits to the [C II] 158 µm line (Table 1). The blue arrows indicate the position of H₂O lines that could potentially contaminate the OH absorption (see the text for details). Note that the ordinate scale varies from line to line and from source to source, and that the plotted velocity range also varies with source.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the observed OH and CO outflow properties. The equivalent width of the (a) OH119 and (b) OH79 doublets, as a function of the CO flux in the blueshifted line wing multiplied by the square of distance. The OH $W_{eq}$ have been measured over the same velocity ranges as the CO fluxes. The dotted lines in panels (a) and (b) indicate chi-square fitting in the log-log plane, giving $\log_{10}(W_{eq} \text{OH119}) = (-3.56 \pm 0.53) + (1.10 \pm 0.10) \log_{10}(S_{\text{CO1-0}} D^2)$ (correlation coefficient of 0.79) and $\log_{10}(W_{eq} \text{OH79}) = (-3.17 \pm 0.47) + (0.91 \pm 0.09)\log_{10}(S_{\text{CO1-0}} D^2)$ (correlation coefficient of 0.93). (c) Comparison of the maximum blueshifted velocity of OH119 (corrected for the PACS spectral resolution at 119 µm) and CO (1–0). The dotted line indicates equal velocities. Abbreviated source names are indicated. The CO data are taken from Cicone et al. (2012, 2014) and A. Gowardhan et al. (2017, in preparation).</figcaption></figure>

with the OH doublets. In Figures 2(a) and 2(b) we plot the equivalent widths of OH119 and OH79 as a function of the CO (1–0) flux in the blueshifted line wing multiplied by the square of distance (a quantity proportional to the outflow mass). The OH equivalent widths have been calculated for the same velocity ranges as used for the CO blueshifted line wings given in Cicone et al. (2012, 2014) and A. Gowardhan et al. (2017, in preparation). There is a hint of correlation between $W_{eq}$ (OH119) and $S_{\text{CO1-0}} D^2$ (correlation coefficient of 0.79), though with notable dispersion. Since the OH119 absorption strength is basically measuring the covering factor of the 119 µm continuum by the outflow ($f_{119}$; see Section 4.3), the trend in Figure 2(a) generally indicates higher outflow mass with increasing covering factor, but the correlation is limited most likely because of the high optical depth of OH119 at most velocities. A better linear correlation is indeed found between $W_{eq}$ of the more optically thin OH79 absorption strength and $S_{\text{CO1-0}} D^2$ (slope of $0.91 \pm 0.09$ in the log–log plane, and correlation coefficient of 0.93). Therefore, the absorption strength of the OH79 doublet is sensitive to the mass of the outflow, from which we may expect to be able to infer reliable values for the outflow energetics. A trend was also found between the outflowing gas mass inferred from the [C II] 158 µm line wings and that from OH (discussed in Section 5; see Janssen et al. 2016).

Although OH119 might be expected to have moderate optical depths at the highest outflow velocities, the maximum outflow velocities observed on the blue sides of the OH119 and CO lines show intriguing differences (Figure 2(c)), with OH showing higher blueshifted velocities than CO in IRAS 03158+4227, Mrk 231, and IRAS 23365+3604 and lower velocities in Mrk 273. This is consistent with the fact that the two species may probe, at least at the highest outflowing velocities, somewhat different components of the galaxies. We may expect that OH generally traces more inner regions than CO: as the expanding gas breaks into clumps and moves away from the central region, the collisionally excited CO emission will still be observed if the density remains high enough within the clumps (i.e., if they are efficiently confined by the interclumpionized medium), but these clumps will cover a decreasing fraction of the far-IR continuum and will thus produce lower absorption and emission in the OH molecules that are

radiatively excited. On the other hand, CO is less sensitive to a compact outflowing shell with small radius and mass, even with enough column density to ensure detection in OH119.

### 4.3. The Covering Factor of the 119 µm Continuum

Given that the OH119 doublet is optically thick, its absorption strength constrains the covering factor of the 119 µm continuum by the outflowing OH. Figure 3 shows the peak absorption (panel (a)) and emission (panel (b)) strengths as a function of the OH119 peak absorption velocity ($V_{\text{p,abs}}^{\text{OH119}}$) for all the observed ULIRGs that show the peak absorption more blueshifted than $-80 \text{ km s}^{-1}$. We thus exclude here the sources that peak in OH119 at systemic velocities, because the strength of the peak absorption in these galaxies is determined by non-outflowing gas components. Figure 3(a) indicates that the strength of the peak absorption in OH119 attains a maximum of 44% in IRAS 03158+4227 and is as low as ≈10% in IRAS 05189–2524. Most sources have, however, peak absorption troughs of 20%–30% of the continuum (see also Figure 1), with hints of increasing absorption strength with increased blueshifted velocity.

The absorption strengths are significantly lower than those predicted from simple models of an outflowing spherical shell surrounding and fully covering a source of far-IR emission (GA14; see also Section 5). For instance, a shell with a velocity field varying linearly with radius, with gas velocities of $v_{\text{int}} = 400$ and $v_{\text{out}} = 200 \text{ km s}^{-1}$ at the inner and outer radius and velocity dispersion of $\Delta V = 100 \text{ km s}^{-1}$, generates in OH119 a P Cygni profile with peak absorption of ≈60% of the 119 µm continuum at $V_{\text{p,abs}}^{\text{OH119}} = -260 \text{ km s}^{-1}$. This is about twice the absorption observed in Mrk 231, IRAS 14378–3651, and IRAS 23365+3604, where the OH119 absorption peaks at similar velocities (Figure 3(a)). Since the OH119 doublet is optically thick, and thus insensitive to the OH column density at least at velocities close to the maximum absorption trough, the discrepancy indicates that OH only covers a fraction (half in the example above) of the observed 119 µm continuum. Hence, the covering factor of the galaxy 119 µm continuum by the outflowing OH can be estimated.

The observed values in Figure 3 are then compared with simple modeling results (colored lines and symbols) described
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Strength, relative to the continuum at 119 µm (C₁₁₉), of (a) the peak absorption feature ($I - S_{p,abs}^{OH119}/C_{119}$) and (b) the peak emission feature ($S_{p,ems}^{OH119}/C_{119} - 1$) of the OH119 doublet as a function of the velocity of the peak absorption feature ($V_{p,abs}^{OH119}$). Values are obtained from Gaussian fits to the emission and absorption features. Only sources with $V_{p,abs}^{OH119} < -80$ km s⁻¹ are included; abbreviated names are indicated. The green ($v_{out} = 200$ km s⁻¹) and blue ($v_{out} = 30$ km s⁻¹) curves show the results of single-component radiative transfer models for two covering factors of the 119 µm continuum by the outflowing OH, $f_{119} = 0.5$ (solid curves) and $f_{119} = 0.15$ (dashed curve; see the text for details).</figcaption></figure>

in Section 5. Rather than the multicomponent models generated in Section 5 to fit simultaneously all observed OH line profiles, we use here single models to explain the regularities observed in the data set of Figure 3, with no attempt to fit the OH119 line profile. We adopt a central continuum source with $T_{dust} = 55$ K and $\tau_{100} = 1$, surrounded by an extended shell of outflowing gas, with $v_{out} = 30$ and $200$ km s⁻¹ (blue and green curves, respectively), variable $v_{int} \le 900$ km s⁻¹, and velocity dispersion $\Delta V = 100$ km s⁻¹. From the emergent theoretical spectra, convolved with the PACS resolution, $V_{p,abs}^{OH119}$ and the strengths of the peak absorption and emission were calculated and multiplied by the appropriate factor $f_{119}$, the covering factor of the 119 µm continuum by the outflowing OH, to match the observations. Figure 3 indicates that most of the observed OH119 absorption and emission strengths can be explained with $f_{119}$ ranging from 0.15 (dashed curves) to 0.50 (solid curves).

P Cygni profiles occur when the absorbing and approaching gas in front of a continuum source and the emitting and receding gas behind it radiatively decouple from one another owing to the expansion. This explains the increase of the peak strengths with increasing outflow velocities for fixed $f_{119}$. In spite of the varying covering factors that may be expected to be present in our sample, the effect of this decoupling is apparent in the data. For high outflow velocities and an outflow significantly more extended than the continuum source, the absorption and emission features are expected to be similar, as

the OH merely redistributes the continuum in velocity space. The strength of the emission feature is decreased with decreasing velocities, due to extinction at 119 µm by the continuum source, which accounts for the fact that the absorption feature is stronger than the emission one. For low velocities, only the absorption feature is expected to be detectable. The modeling values in Figure 3 do not sensitively depend on fixed parameters such as $T_{dust}$, $\tau_{100}$, or $N_{OH}$, but depend on the gas velocity dispersion and of course on geometry, as our values apply to spherical symmetry.

The inferred range $f_{119} = 0.15-0.50$ is similar to the range obtained from the multicomponent analysis in Section 6 (0.17–0.67; the highest value is found in IRAS 03158+4227). These values are significantly lower than 1, indicating that the outflowing OH covers in most sources only a fraction of the observed 119 µm continuum. On the other hand, $f_{119} \sim 0.50$ appears to be too high to be produced by a collimated jet, thus favoring wide-angle coverage as recently found in the outflow of Mrk 231 from high angular resolution CO (2–1) observations (Feruglio et al. 2015). In IRAS 17208–0014, the outflow observed in CO (2–1) at the highest blueshifted velocities ($v < -450$ km s⁻¹) appears to be relatively collimated (García-Burillo et al. 2015), while the OH119 absorption at the same velocities is shallow (≤2%). This independently suggests that OH119 is primarily sensitive to wide-angle outflows. However, the sources in Figure 3(a) with relatively low absorption troughs, IRAS 05189–2524 and IRAS 20551–4250 (also Mrk 273; IRAS 09022–3615 is optically thin in the far-IR), may have a molecular outflow that is significantly collimated. In any case, the fractional absorption may be attributed to clumpiness of the outflowing gas, to biconical (two-lobed) structure, and to the fact that some far-IR-emitting regions of the host galaxy are not affected by the outflow.

## 4.4. Outflowing Gas and Buried Sources

GA15 showed that OH65 is a unique tracer of warm and optically thick cores that account for a significant fraction (if not most) of the (U)LIRG’s luminosity, confirming early indications that ULIRGs are optically thick in the far-IR (Downes et al. 1993; Solomon et al. 1997). The extreme properties inferred for these cores may imply that they represent the most deeply buried stages of AGN-starburst coevolution. The OH outflows are expected to emanate from these active central cores. For this reason we examine here the relationship between the OH65 absorption (at central velocities) and the velocity of the outflows as measured with the most sensitive ground-state OH doublet, OH119. $W_{eq}(OH65) \gtrsim 20$ km s⁻¹ indicates warm ($T_{dust} \gtrsim 60$ K) and optically thick (continuum optical depth at $100$ µm $\gtrsim 1$, or $N_H \gtrsim 10^{24}$ cm⁻²) cores. Results are shown in Figure 4, where the outflow velocity is characterized by $V_{84}(OH119)$, the velocity below which 84% of the absorption in the OH119 doublet takes place (Veilleux et al. 2013; Stone et al. 2016).

In the $V_{84}(OH119)-W_{eq}(OH65)$ plane, galaxies occupy three distinct regions, indicated in Figure 4 with shaded rectangles: galaxies with the highest $|V_{84}| \gtrsim 500$ km s⁻¹ are also strong in OH65 (region I), establishing the connection between the buried and warm cores and the high-velocity outflows. Since the maximum outflowing velocity in OH is correlated with the AGN luminosity (Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Stone et al. 2016), this relationship suggests that the warm material probed by the OH65 doublet
<figure><img src="image_6.png" /></figure>

**Figure 4.** Equivalent width of the high-lying OH65 doublet at central velocities (between −200 and +200 km s⁻¹, W<sub>eq</sub>(OH65); GA15) as a function of V<sub>84</sub>, the velocity below which 84% of the absorption in the OH119 doublet takes place (from Veilleux et al. 2013). We do not include NGC 6240 in this plot because of the presence of a blueshifted emission feature beyond the OH119 P-Cygni profile (J. Fischer et al., in preparation). Red, green, blue, cyan, and black colors indicate Ht, LINER, Seyfert 1, Seyfert 2, and unclassified optical spectral types, respectively (from Veilleux et al. 1995, 1999, 2009; Kim et al. 1998; Rupke et al. 2005b; García-Marín et al. 2006; Véron-Cetty & Véron 2006; or NED/SIMBAD). Circles and triangles indicate sources with fractional AGN contribution to the bolometric luminosity of α<sub>AGN</sub> ≥ 40% and <40%, respectively, as derived from f15/f30 (Veilleux et al. 2009). Three regions (I, II, and III) marked with shaded rectangles are identified in this plane (see the text). The highest outflow velocities, as measured from the sensitive and ground-state OH119 doublet, are found in sources with strong OH65 absorption (i.e., in buried and warm sources), but the reciprocal does not apply: sources with high OH65 absorption are also found with low |V<sub>84</sub>|. In addition, two sources with weak OH65 absorption have relatively high |V<sub>84</sub>| ~ 500 km s⁻¹.

ultimately provides the conditions for copious BH accretion and strong AGN feedback.

Nevertheless, the converse does not occur, because a number of sources with very strong OH65, including three very extreme sources (IRAS 15250 + 3609, IRAS 07251–0248, and the LIRG NGC 4418), have low |V<sub>84</sub>| (region III). Although most ULIRGs in region III are associated with a low attributed contribution of the AGN to the luminosity, the lack of high OH velocities in these sources does not fully preclude an energetically dominant AGN, because the velocities are expected to decrease with increasing column of gas that is accelerated, and the optically thick cocoons in these sources may cover a large solid angle.

In addition, at least two sources have relatively high |V<sub>84</sub>| ~ 500 km s<sup>-1</sup> but are weak in OH65 (region II; IRAS 10565 + 3609 and IRAS 13120-5453). The optical spectral types, denoted by symbol color in Figure 4, and the AGN fractions inferred from the mid-IR continua, denoted by symbol shape, do not account for the differences between sources in the three regions, although AGN-dominated sources are more clearly identified in regions I and II than in region III. Most region III sources have a redshifted OH119 and OH79 peak absorption and/or redshift absorption in the [O I] 63 µm line (IRAS 15250+3609, NGC 4418, IRAS F17207-0014, Arp 220, IRAS 22491-1808), indicating a complex velocity field with ground-state far-IR molecular line shapes dominated by non-outflowing foreground gas that, at least in some cases, appears to be inflowing (see Section 4.8). We thus focus our modeling in Section 5 on sources in regions I and II (see also Section 3.1).

4.5. Expansion of the Nuclear Regions

In this section we examine the evidence for low-velocity outflows detected in the excited OH65 and OH84 doublets.

Because OH65 is a key tracer of warm and optically thick nuclear regions of galaxies (GA15), outflowing signatures in OH65 may be expected to indicate expansion motions of the regions preferentially exposed to feedback from an SMBH or an extreme starburst. However, one crucial point when studying velocity shifts in OH65 is the redshift, or velocity center, of the nuclear region. Because OH65 is a very high lying transition, it is only excited in regions with a strong far-IR field and thus with high extinction in the far-IR, and hence the redshifted reemission feature from behind the continuum source is mostly obscured and not detected. Only one source, IRAS 08572+3915, shows direct evidence for a true P Cygni in OH65 (Figure 1). For all the other sources we must infer the nuclear molecular outflows by measuring velocity shifts and identifying absorption wings in OH65.

In Figure 5, we compare the profiles of the [C II] 158 µm line, the blue and red components of the OH65 doublet separately, and the OH79 doublet in six sources that show indications of nuclear outflowing gas. Velocity profiles for OH65 are then presented for both the blue and red components in separate panels, while the profiles for OH79 are presented for the blue component of the doublet (the zero velocity of the red component is also indicated by a dashed vertical line at ≈230 km s⁻¹). To infer the velocity of peak absorption in both the blue and red components of OH65 (V<sub>p,abs</sub><sup>OH65</sup>), “truncated” Gaussian fits have been applied to each of the OH65 components in such a way that contamination by the other *l*-doubling component and by high-velocity wings is minimized, and they are also shown in Figure 5 in blue and red, respectively. We note that V<sub>p,abs</sub><sup>OH65</sup> of the blue component may actually be lower (more negative) than the fitted value in some sources, due to contamination by the high-lying H<sub>2</sub>O 6<sub>25</sub>-5<sub>14</sub>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** [C II] 158 µm, the blue and red components of OH65, and the OH79 line shapes in six sources where evidence for nuclear outflowing gas is found. Redshifts used for all profiles have been taken in all cases from Gaussian fits to [C II] 158 µm (in green). The OH65 blue and red components have also been fitted with Gaussians (blue and red curves), showing significant velocity blueshifts or asymmetries relative to [C II] 158 µm in all these cases. The blue arrows indicate the position of the H₂O 6$_{25}$–5$_{14}$ line, which could be contaminating the blue component of some OH65 profiles. The cyan arrows in the lower panels indicate the approximate peaks of absorption and emission in the OH79 P Cygni line profiles, while the OH79 absorption features have their counterpart in the blueshifted OH65 doublet absorption features.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Velocity of peak absorption of the OH65 ($V_{p,abs}^{OH65}$, upper) and OH84 ($V_{p,abs}^{OH84}$, lower) doublets relative to the velocity of the [C II] 158 µm line, in all sources where the OH65 doublet is detected (GA15), including the 14 ULIRGs of the modeling sample. Blue and red symbols indicate peak absorption velocities of the blue and red components of the doublets, obtained from Gaussian fits to each doublet component. $V_{p,abs}^{OH65}$ of the blue component is systematically lower than that of the red component, most likely due to contamination of the blue component by the high-lying H$_2$O $6_{25}-5_{14}$ line. Error bars are $2\sigma$, and the error bars on the right-hand side indicate potential “skewness” error from pointing drifts, which we adopt as 0.01 µm (i.e., ≈45 km s$^{-1}$ for OH65; see Section 4.7.2 of the PACS Observer’s Manual, http://herschel.esac.esa.int/Docs/PACS/html/pacs.om.html). Abbreviated source names are indicated; those in magenta indicate that clear signatures of outflowing gas (P Cygni or wings) in the ground-state OH119 doublet are observed. Sources are sorted according to the component-average value of $V_{p,abs}^{OH65}$.</figcaption></figure>

line, the position of which is indicated with blue arrows. Therefore, we mostly rely on the velocities inferred from the OH65 red component to infer the detection of outflowing gas in OH65. The bottom panels in Figure 5 show the OH79 profiles in these sources, which all show P Cygni line shapes.

The OH65 profiles are broad, extending all the way to $-500 \text{ km s}^{-1}$ on the blue wing, and in some cases beyond the limits of the [C II] emission at zero intensity. (Note, however, that the [C II] line wings observed in Janssen et al. (2016) are hardly perceptible in Figure 5, due to the linear scale used for the ordinates.) In addition, clear velocity shifts of 50–200 km s$^{-1}$ between the OH65 peak absorption and the [C II] peak emission are seen in Figure 5. In IRAS 08572+3915, IRAS 14348–1447, Mrk 273, and IRAS 23365+3604, the OH65 feature is blueshifted relative to [C II]. In IRAS 20551–4250, OH65 peaks at the systemic velocity but shows an asymmetric blue wing in both OH65 components up to $-400 \text{ km s}^{-1}$.

Blueshifts in OH65 relative to [C II] can be attributed to (i) differences between the nuclear redshift(s) and that of the bulk of the host galaxy as traced by [C II], possibly due to nuclear motions associated with the merger; (ii) noncircular motions in the nuclear region due to bar-like or oval distortions and warps (leading to elliptical orbits; e.g., Sanders 1989); or (iii) OH65 absorption that is tracing expanding (radial) gas flows in the nuclear region. While the first two possibilities should be considered, particularly in ULIRGs where the rotating, merging nuclei might indeed generate such velocity shifts, comparison of the OH65 and OH79 profiles suggests that the third possibility is clearly favored in some sources. The OH79 line profiles show P Cygni profiles in all sources presented in Figure 5, and in IRAS 08572+3915, IRAS 14348–1447, and more tentatively in Mrk 273 the blueshifted absorption features in OH79 have counterparts both in

emission at positive velocities in OH79 (indicated by the cyan arrows) and in absorption at the same velocities in the OH65 doublet. It is this OH79–OH65 correspondence in the blue part of the profiles, together with the P Cygni in OH79 (as well as the preponderance of OH65 blueshifted absorption as discussed below), that enables us to distinguish between rotating and outflowing motions in the nuclear regions from the OH65 profiles. The inferred outflows do not preclude, however, the simultaneous incidence of a non-outflowing component in OH65.

Figure 6 shows $V_{p,abs}^{OH65}$ in all sources that are detected in OH65 (GA15, thus including more than the 14 sources of our modeling sample), as well as the corresponding quantity for OH84 in those sources where the doublet is available. With the exceptions of IRAS 07251–0248 and IRAS 17208–0014 (the latter showing evidence for outflowing gas in CO; García-Burillo et al. 2015), all other sources in Figure 6 with significant blueshifted OH65 absorption show a P Cygni profile in OH119, which is indicated in Figure 6 by magenta colored source names, while the incidence of a P Cygni OH119 line shape in sources with no OH65 blueshift is significantly lower. Three sources, IRAS 14378–3651, IRAS 20551–4250, and IRAS 12112+0305, show P Cygni profiles in OH119, but their OH65 doublets peak at central velocities. IRAS 20551–4250, however, has a blueshifted wing in OH65 that is also attributable to outflowing gas (Figure 5). The intermediate OH84 doublet confirms the velocity shifts seen in OH65 in most cases. Specifically, IRAS 10565+2448 is weak in OH65 (but detected; GA15), and its $V_{p,abs}^{OH65}$ value has large uncertainties. However, shifts of $\sim 130 \text{ km s}^{-1}$ are clearly seen in OH84. In IRAS 23365+3604, the velocity shifts in OH84 are larger than in OH65, as is also the case in Mrk 231 (GA14).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Detailed comparison of the blueshifted side of the OH119 (black histograms), OH79 (green), OH84 (blue), and OH65 (red) line profiles observed in the 14 sources that are modeled in Section 5. We only show the relevant velocity ranges of each doublet for clarity. The H~2~O *4~2~3--3~1~2* line is most likely producing absorption at velocities <~−1400 km s⁻¹ in the OH79 spectra of IRAS 03158+4227 and IRAS 20100−4156.</figcaption></figure>

Identifying outflowing gas in excited OH transitions may miss some molecular outflow sources, such as Arp 220 where clear indication of low-velocity outflowing molecular gas is seen in several species, including OH (Sakamoto et al. 2009; González-Alfonso et al. 2012; Tunnard et al. 2015; Martín et al. 2016).

A detailed comparison of the blueshifted absorption wings in OH119, OH79, OH84, and OH65 is shown in Figure 7 for the 14 sources modeled in Section 5. Besides Mrk 231, IRAS 03158+4227, IRAS 20100−4156, and IRAS 08572+3915, where the high-velocity outflowing gas in OH65 is evident, Mrk 273 and IRAS 14348−1447 show similar line shapes in OH65 and in the P Cygni OH119–OH79 doublets, indicating that OH65 is also tracing outflowing gas. Similarly, the OH65 line wing in IRAS 23365+3604 has its counterpart in OH119 and OH84, both showing P Cygni profiles (Figure 1). We conclude that, in addition to Mrk 231 (GA14), IRAS 03158+4227, and IRAS 20100−4156, the six sources in Figure 5 show evidence for outflowing gas in the high-lying OH65 transition. In addition, IRAS 09022−3615 and IRAS 10565+2448, although weak in OH65 (GA15), show hints of blueshifted OH65 absorption as well (Figure 1).

Could these velocity shifts observed in OH65 be associated with merging motions, rather than outflowing or inflowing gas? We compare the peak velocities of the OH65 doublet (in absorption) to those of the [O I] 63 µm (hereafter OI63) emission line in Figure 8(a). These velocities tend to be anticorrelated: when the OH65 peak is blueshifted, the OI63 peak tends to be redshifted and vice versa. Figure 8(b) shows $V_{\text{p,abs}}^{\text{OH65}}$ as a function of the OI63 red-blue asymmetry. Consistent with panel (a), panel (b) shows that, when OH65 is blueshifted, the OI63 line usually has a “red-type” profile (and vice versa), i.e., the profile is distorted by a reduction of the intensity on the blue side. This is fully consistent with the scenario of atomic oxygen also outflowing with the OH, and thus absorbing the continuum and line emission arising from behind.<fnref n="14" /> It also rules out that the OH65 shifts are the result of global nuclear motion relative to the more tenuous and

<footnote n="14">¹⁴ Although some continuum is absorbed, the OI63 line is not seen in absorption because it has an Einstein A-coefficient of ~10⁻⁴ s⁻¹ (much lower than OH65, ≈1.2 s⁻¹) and is thus easily excited through collisions in warm and dense nuclear regions. The OI63 line is only seen in absorption in sources where high columns of oxygen are found in extended, low-density regions in front of the 63 µm continuum source (GA12; Falstad et al. 2015).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Velocity of peak absorption of the OH65 doublet (blue component of the doublet) as a function of (a) the velocity centroid of the [O I] 63 µm emission line and (b) the [O I] 63 µm red–blue asymmetry, which is defined as $(F_+ - F_-)/(F_+ + F_-)$, where $F_+(F_-)$ is the flux measured at positive (negative) velocities. The dotted line in panel (a) indicates equal velocities. The error bars in the upper left corner of each panel indicate potential “skewness” uncertainties due to pointing drifts. The two sources in our sample with evidence for inflowing gas, IRAS 11506–3851 and IRAS 15250+3609 (see Section 4.8), are labeled in panel (b).</figcaption></figure>

extended region probed by the [C II] line, because in the latter case one would expect the OH65 and OI63 to peak at the same velocities, contrary to the observed trends.

To put these results on a first quantitative framework, the moderate velocities found here are similar to those predicted by King (2003) in his analytical modeling of outflows driven by a wind with momentum flux comparable to the Eddington-limited radiation field ($L_{edd}/c$), i.e., $v_m = (G L_{edd}/2f_g \sigma^2 c)^{1/2}$, giving $\approx 80 \text{ km s}^{-1}$ for $L_{edd} = 10^{12} L_\odot$, velocity dispersion $\sigma = 200 \text{ km s}^{-1}$, and gas fraction $f_g = 0.16$. The derivation of $v_m$ in King (2003) assumes an isothermal sphere of dark matter, involving an outflowing gas mass of $3 \times 10^8 \times (\sigma/200 \text{ km s}^{-1})^2 \times (r/100 \text{ pc}) M_\odot$ and a column density of $2 \times 10^{23} \times (\sigma/200 \text{ km s}^{-1})^2 \times (r/100 \text{ pc})^{-1}$. It ignores the potential well, implicitly assuming that the gas is supported by radiation pressure or rotation.

These results show that the nuclear regions of galaxies with high-velocity outflows (as seen in OH119) also show slower nuclear outflows, i.e., the far-IR “photospheres” where OH65 is formed are (partially or completely) expanding at low velocities. Since OH65 requires high columns to be detected, star formation may proceed at significant rates in this expanding material, with the consequence that star formation propagates outward.

## 4.6. Velocity Components

Figure 9 summarizes the relationship between $V_{p,abs}^{\text{OH65}}$ and the velocities inferred from OH79 and OH119. Some sources have $V_{p,abs}^{\text{OH65}} \approx V_{p,abs}^{\text{OH119}} \approx V_{p,abs}^{\text{OH79}} < 0$ (Figure 9), indicating the presence of nuclear outflowing gas with the three lines tracing basically the same gas. In other sources, OH65 peaks at different velocities than OH119 and OH79, indicating that different components are traced by the OH doublets. It is also

apparent from Figures 7 and 9 that the ground-state OH119 doublet usually traces gas at more blueshifted velocities than the other doublets in most sources. The peak absorption velocities of the four doublets are compared in Figure 10 for all sources where OH65 is detected, clearly showing a tendency for the high-lying OH65 and OH84 lines to trace gas at lower velocities than the ground-state OH119 and OH79 doublets. In the ground-state lines, absorption and reemission at systemic velocities take place over large volumes and tend to cancel each other (see Figure 4 in Spoon et al. 2013; Figure 8 in GA14), leading to clear P Cygni profiles due to absorption and reemission at higher velocities. This cancellation cannot happen in the high-lying lines, in which lower-velocity shifts indicate higher column densities of gas with lower outflowing velocities.

## 4.7. Compact and Extended Outflows

Despite the relatively low spatial resolution of Herschel/PACS, the OH excitation can nevertheless constrain the spatial extent of the outflows. In IRAS 08572+3915, Mrk 273, and IRAS 20551–4250, the absorption in the excited OH84 transition is almost as strong as (or even stronger than) the OH119 absorption (Figure 7). In these objects, essentially all outflowing gas traced by the ground-state and optically thick OH119 is significantly excited by the far-IR field, indicating relatively small distances to the nuclear region and thus a compact outflow. Mrk 273 in particular was found by Cicone et al. (2014) from high-resolution CO (1–0) observations to have the smallest size among all their detected ULIRGs. In the other sources, OH119 is much stronger than OH84 at least at some negative velocities, indicating the presence of a spatially extended outflow component.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Velocity of peak absorption of the OH65 doublet, $V_{p,abs}^{OH65}$, as a function of the velocity of peak absorption (a) of the OH119 doublet (blue transition) and (b) of the OH79 doublet (blue transition). Blue and red symbols indicate $V_{p,abs}^{OH65}$ of the blue and red transitions of the OH65 doublet. Abbreviated names are indicated for some sources. The dotted lines indicate equal velocities, $V_{p,abs}^{OH65} = V_{p,abs}^{OH119}$ and $V_{p,abs}^{OH65} = V_{p,abs}^{OH79}$. The error bars in the lower right corner of each panel indicate potential “skewness” uncertainties due to pointing drifts.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Velocity of peak absorption ($V_{p,abs}$) of the OH119, OH79, OH84, and OH65 doublets in all sources where OH65 and or OH84 are detected. [$V_{p,abs}$] decreases with increasing excitation energy. The error bar indicates potential “skewness” uncertainties from pointing drifts. Abbreviated source names are indicated.</figcaption></figure>

### 4.8. Inflows

Besides the prevalent outflowing motions traced by OH, inflow signatures are also seen in a few (U)LIRGs as observed with Herschel/PACS. In addition to the previously studied case of NGC 4418 (González-Alfonso et al. 2012), where the inflow has also been inferred from high angular resolution observations (Costagliola et al. 2013; Sakamoto et al. 2013), the case of Zw 049, where an inverse P Cygni profile is also observed in the [O I] 63 µm line (Falstad et al. 2015), the inverse P Cygni profile of OH119 in Circinus (Stone et al. 2016), and the inflow observed in OH119 around Arp 299a (Falstad et al. 2017), the clearest examples are the LIRG IRAS 11506–3851 and the ULIRG IRAS 15250+3609. Their OH spectra are shown in Figure 11,

where the redshifts are obtained from Gaussian fits to the [C II] 158 µm profiles. In IRAS 11506–3851, the OH119, OH79, and OH84 peak absorption velocities are redshifted relative to [C II], while the highest-lying OH65 velocity is more centered than the others but still slightly redshifted (see also Figure 10). The strongest evidence for inflow comes from OH79, which shows an inverse P Cygni profile with an emission feature at ~−100 km s⁻¹ and a concomitant asymmetry between the red and blue components of the doublet. The emission feature appears to have a weak counterpart in OH119 as well. In addition, the high-lying OH84 line shows hints of a blueshifted line wing in absorption as observed in CO (Pereira-Santaella et al. 2016), possibly indicating the coexistence of an outflow in the nuclear region. In IRAS 15250+3609, the observed OH119, OH79, and OH65 doublets are also redshifted relative to [C II], and OH119 shows a broad emission feature between ~−300 and −1000 km s⁻¹. Rather than a simple infall around a central warm source, the extreme velocities in this source rather suggest an approaching double-nucleus system, with the strongest source of far-IR radiation associated with the nucleus closer to the observer and illuminating the secondary nucleus responsible for the OH119 emission feature, which may be disrupted by tidal forces. Interestingly, this is the source with the highest equivalent width in OH65 (GA15; see also Figure 4), and a blueshifted line wing in this doublet also suggests the simultaneous occurrence of a nuclear outflow, as observed in HCO⁺ 3–2 (Imanishi et al. 2016). Both IRAS 11506–3851 and IRAS 15250+3609 show a blue-shift profile in the OI63 line, i.e., the line shows a reduction of the intensity on the red side (Figure 8(b)). These sources will be studied in more detail in future work.

## 5. Radiative Transfer Models

### 5.1. A Library of Model Components

In order to quantitatively estimate the energetics associated with the observed outflows, we have developed a library of model components with the radiative transfer code described in González-Alfonso & Cernicharo (1997, 1999). In short, the method calculates the statistical equilibrium populations in all shells of a spherically symmetric source, as well as the emergent continuum and velocity profiles of all lines after
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Observed OH line shapes in IRAS 11506-3851 and IRAS 15250+3609, where inflow signatures are detected. The vertical gray lines in all panels indicate the rest position of the two *l*-doubling components, with the velocity scale taken relative to the blue component of each doublet. The arrow in the OH65 panel indicates possible contamination by a high-lying H<sub>2</sub>O line.</figcaption></figure>

convolution with the PACS spectral resolution. The approach is thus non-LTE and nonlocal, takes into account the radiative pumping by both local and nonlocal dust (the dust is mixed with the gas), and accurately includes the effects of both line and continuum opacity effects and line overlap between the *l*-doubling components of OH. The gas velocity field can be composed of pure turbulence and radially outflowing motions with or without velocity gradients. The radial velocity is allowed to vary continuously within any shell. We developed two types of models to fit the continuum and line emission of our sample sources: “CORE” and “ENVELOPE” models (Figure 12). They are described below.

**CORE models.** These models are intended to describe the continuum emission and line absorption at systemic velocities produced in the nuclear regions of galaxies (GA15). They consist of a source of dust with temperature *T*<sub>dust</sub> and radial optical depth at 100 µm τ<sub>100</sub>. The OH, mixed with the dust, has turbulence ΔV, and its column is determined by τ<sub>100</sub> according to Equation (1) in González-Alfonso et al. (2014a), with an OH abundance relative to H nuclei of *X*<sub>OH</sub> = 2.5 × 10<sup>-6</sup>. The adopted abundance is, within a factor of ~3, consistent with the value inferred from multitransition observations of OH in the Galactic Sgr B2 and Orion KL outflow (Goicoechea & Cernicharo 2002; Goicoechea et al. 2006) and in buried galaxy nuclei (GA12; Falstad et al. 2015), and with chemical models

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Schematic representation of the two types of model components, CORE and ENVELOPE, and independent parameters (Section 5.1) used to fit the observed OH line profiles.</figcaption></figure>

of dense photodissociation regions (the peak value; Sternberg & Dalgarno 1995) and of cosmic-ray- and X-ray-dominated regions (Meijerink et al. 2011). The source structure includes expansion with constant velocity *v*<sub>out</sub> to describe the low-velocity blueshifts of the peak absorption of OH65 that indicate outflowing gas in the nuclear region (Section 4.5), and it is divided into shells to account for the varying excitation with
depth into the source. The free parameters of the CORE components that are varied from model to model are then $T_{\text{dust}}$, $\tau_{100}$, $\Delta V$, and $v_{\text{out}}$.

**ENVELOPE models.** These models also consist of a central source of only dust that is surrounded by an expanding shell of gas and dust (which is also divided into shells to account for the varying physical conditions with distance to the central source). The gas velocity is allowed to vary linearly with radius across the expanding shell (GA14), and thus the free parameters are $T_{\text{dust}}$ and $\tau_{100}$ of the central source and $R_{\text{out}}/R_{\text{int}}$, $v_{\text{int}}$, $v_{\text{out}}$, $\Delta V$, and $N_{\text{OH}}$ of the expanding shell. The density of dust across the shell is determined by the OH density, assuming also $X_{\text{OH}} = 2.5 \times 10^{-6}$ and a gas-to-dust ratio by mass of 100 (GA14). In all models we assume a constant mass outflow rate, so that $n_{\text{H}} v r^2$ is uniform between $R_{\text{int}}$ and $R_{\text{out}}$ (GA14).

## 5.2. $\chi^2$ Fitting

Our approach relies on the assumption that the OH profiles can be considered a superposition of the single-component models described above. This is questionable mostly at central velocities, where there may be radiative coupling between different gas components, but has the overall advantage of simplicity.

Up to $N_C = 3$ model components were allowed to be combined in order to simultaneously fit the continuum-normalized OH119, OH79, OH84 (when available), and OH65 spectra in a given source. For each combination of components, the modeled continuum-normalized value of line i and velocity channel j is

$$S_{ij} = 1 + \sum_{l=1}^{N_C} f_{il} s_{ijl}, \qquad (1)$$

where the sum extends to the $N_C$ components of the combination. Here $f_{il}$ is the ratio of the continuum flux density of model l to the total (observed) continuum, at the wavelength of line i:

$$f_{il} = C_{il}/C_{i,\text{OBS}}, \qquad (2)$$

where C stands for the continuum flux density and OBS refers to the observed values. Parameter $s_{ijl}$ is the continuum-normalized model prediction minus 1:

$$s_{ijl} = F_{ijl}/C_{il} - 1, \qquad (3)$$

where $F_{ijl}$ and $C_{il}$ are the predicted line flux density and continuum flux density, respectively, for model component l. Our minimization process calculates, for each combination of model components, the values of $f_{il}$ that best fit the observed spectra. For a given model component l, however, the values of $f_{il}$ for different lines are related to each other through

$$f_{il} = f_{\text{OH119},l} \frac{C_{il}/C_{\text{OH119},l}}{C_{i,\text{OBS}}/C_{\text{OH119,OBS}}}, \qquad (4)$$

and thus only the values of $f_{\text{OH119},l}$ are fitted. The reduced $\chi^2$ value for a given combination that we aim to minimize is

$$\chi^2 = \frac{1}{\sum_{l=1}^{N_C} N_v(i) - N_C} \sum_{i=1}^{N_l} \sum_{j=1}^{N_v(i)} \frac{1}{\sigma_i^2} [S_{ij} - E_{ij}]^2, \qquad (5)$$

where $N_l$ is the number of lines fitted in a given source, $N_v(i)$ and $\sigma_i$ are the number of velocity channels of and the rms noise around line i, respectively, and $E_{ij}$ is the observed continuum-normalized value of line i at velocity channel j. The $N_C$

parameters that are varied to minimize $\chi^2$ for each combination are $f_{\text{OH119},l}$, and this is done using standard methods.

In practice, we used the same “line-averaged” value for all $\sigma_i$ in a given source and increased the relative weight of channels in the line wings, in which we are mostly interested. Depending on the observed line shapes in a given source, we selected ($N_{\text{CORE}}$, $N_{\text{ENVELOPE}}$) as (1, 2), or (1, 1). The minimization procedure was done for all possible combinations allowed by the following restrictions: (i) The predicted spectral energy distribution (SED) for the components should be as warm as or warmer than the observed SED. As the OH expands, the molecules are subject to a progressively colder far-IR SED, but not colder than the average galaxy SED as seen by the observer. (ii) In some sources, a maximum $T_{\text{dust}}$ was imposed for the CORE models, to avoid overestimating the continuum at short wavelengths. More than 100 CORE models and more than 400 ENVELOPE models were generated, giving from $\sim 10^5$ to several $\times 10^6$ allowed combinations.

In addition to $\chi^2$, we calculated the corresponding value $\chi^2_{\text{BW}}$ that uses only the velocity channels at blueshifted velocities. A first selection was done including all combinations with a value of $\chi^2$ within 20% of the minimum $\chi^2$. Among all these combinations, the best fit was selected as the combination that gives the minimum value of $\chi^2_{\text{BW}}$. Figures 13–16 show our best-fit composite models for OH119, OH79, OH84, and OH65 in the 14 modeled sources (Section 3).

## 5.3. Approaching the Minimum $\chi^2$

Since the number of input parameters involved in the ENVELOPE models is high, the following three-step strategy to approach the minimum was generally adopted: (i) We first developed a coarse grid of models, common to all sources, by varying the velocity field, $N_{\text{OH}}$, and $T_{\text{dust}}$, for two values of $R_{\text{out}}/R_{\text{int}} = 1.3$ and $1.6$, and keeping $\tau_{100} = 1$ and $\Delta V = 100 \text{ km s}^{-1}$ fixed. A first search of the minimum was then performed, and the grid was additionally refined in some regions of parameter space if no satisfactory fit was found. (ii) The grid was additionally refined by varying the velocity field, $N_{\text{OH}}$, and $T_{\text{dust}}$ around the minimum and also varying to a greater extent $R_{\text{out}}/R_{\text{int}}$ and $\tau_{100}$. A second search of the minimum $\chi^2$ was carried out. (iii) The final step consisted in further refining the model grid around the minimum, varying also $\Delta V$ and with special consideration of $N_{\text{OH}}$. The region around the minimum was then sufficiently populated to determine the error bars associated with the energetics (see Section 5.5).

## 5.4. Continuum and Minimum Size

Our fitting procedure uses the observed continuum-normalized OH spectra and the corresponding theoretical continuum-normalized spectra, as well as the observed far-IR colors in the OH bands as inputs (Equation (4)). Once the value of $f_{\text{OH119},l}$ is fitted, we use the observed absolute value of the continuum flux density at 119 µm, $C_{\text{OH119,OBS}}$, to infer the absolute value of the 119 µm continuum flux density associated with model component l, $C_{\text{OH119},l} = f_{\text{OH119},l} \times C_{\text{OH119,OBS}}$. Since the normalized SED for any component is fully determined for given $T_{\text{dust}}$ and $\tau_{100}$, the inferred $C_{\text{OH119},l}$ value gives a specific prediction for the absolute SED of the underlying far-IR source associated with component l. The absolute calibration of the models is thus based on the far-IR continuum, which in turn is the source of the OH excitation.
<figure><img src="image_16.png" /><figcaption>Figure 13. Model fits to the observed OH119 doublets (shown as yellow-shaded histograms) in 14 sources. The models simultaneously fit the line profiles of all (OH119, OH79, OH84 when available, and OH65) doublets observed in each source through $\chi^2$ minimization using combinations of model components taken from a library. Dashed lines indicate different components (up to three in each source, in blue, green, magenta), and red is total. The two vertical lines in each panel indicate the rest position of the doublet transitions. The velocity scale is relative to the blue transition of the doublet. The fits to the other OH doublets are shown in Figures 14–16, and the predicted far-IR continua are compared with the observed SEDs in Figure 17.</figcaption></figure>

With the absolute SED known, the required solid angle $\Delta\Omega$ of component $l$ is obtained, and hence the minimum size of the component through $R_l = D\sqrt{\Delta\Omega/\pi}$, where $D$ is the distance to the source. This minimum radius implicitly assumes that the OH in each outflow component is not clumpy (we relax this assumption below when calculating the energetics), so that $f_{OH119,l} < 1$ is due to additional continuum from the galaxy not associated with the $l$-component (i.e., associated with other OH components or with far-IR emission unrelated to OH).

We compare the predictions for the galaxy continua with the observed SEDs for all modeled sources in Figure 17, where the red curves indicate the total continua (i.e., the sum of all components, including the non-outflowing one). Our superposition approach is consistent with observations as long as the modeled continuum does not exceed the observed one at any wavelength (it can underestimate the observed continuum if some of it is unassociated with OH). In two sources, IRAS 03158+4227 and IRAS 14348–1447, the modeled continuum slightly exceeds the observed SED at $\lambda < 50$ µm, indicating that the continua from the different components are not independent in these sources.

### 5.5. Estimating the Energetics

Using the thin-shell approximation, the mass of the outflowing gas associated with a given component $l$ is

$$M_{\text{out}} = 4\pi R_l^2 \mu m_H N_H, \quad (6)$$

where $\mu = 1.4$ is the mass per H nucleus and $N_H = N_{OH}/X_{OH}$. Both $R_l (=R_{\text{out}})$ for the CORE models, and $=(R_{\text{int}} + R_{\text{out}})/2$ for the ENVELOPE models) and $N_{OH}$ are directly inferred from the fitting procedure.

Since the continuum of each component underpredicts the observed SED, the observations can also be interpreted in terms of a clumpy distribution, where the source radius $R$ is actually larger than $R_l$, $R = f_c^{-1/2} R_l$, and the outflowing OH only covers a fraction $f_c$ of the underlying continuum source. We determine the minimum value of $f_c$ such that the modeled flux density does not overpredict the observed SED at any wavelength.<fnref n="15" /> The outflowing gas mass is, nevertheless, the same as in Equation (6), which we now express as

$$M_{\text{out}} = 4\pi f_c R^2 \mu m_H N_H, \quad (7)$$

but the mass outflow rate and momentum flux (see below) decrease by a factor of $f_c^{1/2}$. For the CORE models with $v_{\text{out}} > 0$, where the whole source is expanding, we limit $N_H$ to $5 \times 10^{23} \text{ cm}^{-2}$ if the column through the source exceeds that value because this is the column that yields a continuum optical depth of $\sim 1$ at 60 µm.

While we consider $R$ to be the most likely value of the outflow radius, and thus a clumpy distribution is adopted, the

<footnote n="15">¹⁵ Note that $f_c > f_{119}$ in warm components, because the predicted flux density at $\lambda < 100$ µm is closer to the observed one than that at 119 µm.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 14. Model fits to the OH79 doublets in 14 sources. Blue arrows indicate possible contamination by H₂O. See also caption of Figure 13.</figcaption></figure>

value of R₁ is relevant as it indicates the radius below which the OH is predicted to cover the full 4π sr. This value is used in our simple dynamical model in Section 7.

There are two limiting approaches to the estimation of the energetics associated with a given component: the local (maximum) or instantaneous values and the average (minimum) values (see also Rupke et al. 2005c). The local estimates use the inferred density n_H of the expanding gas:

$$ \dot{M}_{\text{loc}} = f_c 4\pi R^2 \mu m_{\text{H}} n_{\text{H}} v = \frac{M_{\text{out}} v}{\Delta R}, \quad (8) $$

$$ \dot{P}_{\text{loc}} = \dot{M}_{\text{loc}} v, \quad (9) $$

$$ \dot{E}_{\text{loc}} = \frac{1}{2} \dot{M}_{\text{loc}} v^2, \quad (10) $$

where $\Delta R = R_{\text{out}} - R_{\text{int}}$ is the thickness of the outflowing shell. These local estimates give the mass, momentum, and energy crossing instantaneously per unit time a spherical surface of radius R, and they can only be considered “time-averaged” if we assume that the outflowing gas extends to r = 0, although we only observe a fraction of it owing to far-IR opacity (extinction) and line opacity effects. This approach was used in Sturm et al. (2011) and GA14. The average values assume, by contrast, that there is no more outflowing gas inside the curtain of dust, so that

$$ \dot{M}_{\text{out}} = f_c 4\pi R^2 \mu m_{\text{H}} \frac{N_{\text{H}} v}{R} = \frac{M_{\text{out}} v}{R}, \quad (11) $$

$$ \dot{P}_{\text{out}} = \dot{M}_{\text{out}} v, \quad (12) $$

$$ \dot{E}_{\text{out}} = \frac{1}{2} \dot{M}_{\text{out}} v^2. \quad (13) $$

Evidently, the average values are a factor of $\Delta R/R$ lower than the local values. In this paper, we use the most conservative $\dot{M}_{\text{out}}$, $\dot{P}_{\text{out}}$, and $\dot{E}_{\text{out}}$ values to characterize the outflows. These are the “time-averaged thin-shell” values in Rupke et al. (2005c), which have also been used by a number of authors describing the energetics of the ionized and neutral phases of outflows (e.g., Arav et al. 2013; Borguet et al. 2013; Rupke & Veilleux 2013a; Heckman et al. 2015), and are most appropriate for comparison with outflow models (e.g., Faucher-Giguère & Quataert 2012; Stern et al. 2016; Thompson et al. 2015). In some studies, a factor of 3 higher values have been used because the emitting spherical (or multiconical) volume is assumed to be filled with uniform density (e.g., Feruglio et al. 2010, 2015; Maiolino et al. 2012; Rodríguez Zaurín et al. 2013; Cicone et al. 2014; Harrison et al. 2014; García-Burillo et al. 2015). For a steady flow with constant velocity, however, we would expect a density at the outer radius only 1/3 that of the average, thus also yielding the expression in Equation (11).

On the other hand, the assumed flow timescale ($R/v$) is overestimated in this expression if a fraction of the gas is loaded (and suddenly accelerated) at some distance from the center, which is relevant here given the compact sizes we infer below. However, if the gas is smoothly accelerated from the center, the flow timescale will be longer than $R/v$ by a factor of ~2. Our assumption of spherical symmetry (the 4π factor) may
<figure><img src="image_18.png" /><figcaption>Figure 15. Model fits to the OH84 doublets in 14 sources. See also caption of Figure 13.</figcaption></figure>

overestimate the energetics by a factor of ~2 in the case of bipolar emission with no gas flowing along the plane of sky. In systems where gas flows primarily along the plane of sky, the energetics will obviously be underestimated. In the above expressions involving v, we use the average of vint and vout.

The total values that characterize a given source are the sum of Mout, Mout, Pout, and Eout over all components (Rupke et al. 2005c),

$$ (M_{\text{tot}}, \dot{M}_{\text{tot}}, \dot{P}_{\text{tot}}, \dot{E}_{\text{tot}}) = \sum_{\text{comp}} (M_{\text{out}}, \dot{M}_{\text{out}}, \dot{P}_{\text{out}}, \dot{E}_{\text{out}}). \quad (14) $$

$\dot{M}_{\text{tot}}$, $\dot{P}_{\text{tot}}$, and $\dot{E}_{\text{tot}}$ are plotted in Figure 18 as a function of Lbol, LAGN, and L* (see also Figure 19).

To estimate the uncertainties of the total values in Equation (14), we calculate the values of Mout, Mout, Pout, and Eout, for all model combinations that yield $\chi^2$ and $\chi^2_{BW}$ to within 20% of their minimum values, and we conservatively use the maximum and minimum values of the energetics of all these combinations to define the error bars. Uncertainties usually range from 50% of the best-fit value to more than a factor of 2 in some cases, and do not include systematic uncertainties due to our adopted OH abundance.

We also calculate the mass-loading factor of each outflow component, $\eta = \dot{M}_{\text{out}}/\text{SFR}$, where SFR is the star formation rate (SFR = $10^{-10} L_{*} M_{\odot} \text{yr}^{-1}$, where L* has units of L*). The starburst and AGN contributions to LIR (see Table 1) are estimated from the flux densities at 15 and 30 $\mu$m according to Veilleux et al. (2009, 2013). The total value of $\eta$ is also obtained by summing the contributions from all outflowing components.

# 6. Modeling Results

## 6.1. Overall Model Fits

Our best fits to the OH119, OH79, OH84, and OH65 doublets are shown in Figures 13–16, respectively. Table 2 lists the parameters of the outflowing components for the best-fit combination, together with the derived energetics. The total energetics of the outflows are listed in Table 3.

In IRAS 09022–3615, where OH65 is not detected and OH84 is not available, we use only one model component to fit OH119 and OH79. In IRAS 13120–5453, where OH65 is also undetected, we use two model components. Similarly, we use only two components in IRAS 19297–0406, where OH79 is undetected. In the rest of the sources, three components are required to obtain a reasonable model fit. One of them is used to match the absorption in the excited OH84 and OH65 doublets at central velocities (CORE models), though the fit sometimes requires expansion velocities of up to 100 km s$^{-1}$. The other two components are needed to fit line absorption and emission at higher velocities ($\gtrsim 200$ km s$^{-1}$). With our limited number of components, some very high velocity wings seen in OH119 (IRAS 03158+4227, IRAS 13120–5453, IRAS 14348–1447, IRAS 14378–3651, and IRAS 23365+3604) are not fitted, because the associated energetics from only OH119 are relatively uncertain.

The overall goodness of the fits to the line shapes indicates that, in spite of the two-lobed structure generally observed for the outflows traced by CO (e.g., Alatalo et al. 2011; Cicone et al. 2014; Feruglio et al. 2015; García-Burillo et al. 2015), the
<figure><img src="image_19.png" /><figcaption>Figure 16. Model fits to the OH65 doublets in 14 sources. Blue arrows indicate possible contamination by H₂O. See also caption of Figure 13.</figcaption></figure>

spherical symmetry implicit in the model is, as a first approximation, a reliable approach to the outflow systems, in line with the wide-angle outflow geometry inferred from OH119 detection statistics (Veilleux et al. 2013; see also Section 4.3). However, some sources show indications of significant departures from sphericity, which is best seen in the fits to the OH119 doublet in Figure 13. In strict spherical symmetry, an OH119 redshifted emission feature that is weak relative to the absorption feature can only be attributed to far-IR extinction of the emission behind the continuum source, and detection of OH65 (a tracer of optically thick continuum at 119 µm) would be consistent with this explanation. In IRAS 13120–5453, however, there is essentially no reemission at redshifted velocities, but OH65 is undetected. The modeled OH119 profile indeed strongly overpredicts the redshifted emission, indicating that the gas responsible for the absorption is not accompanied by comparable amounts of OH on the other side of the continuum source. Significant overpredictions of the OH119 redshifted emission are also seen in IRAS 03158+4227 (absorption by CH⁺ may reduce the OH119 emission feature in this source; Section 4.1), IRAS 10565+2448, IRAS 19297–0406, and IRAS 23365+3604. On the other hand, the opposite effect is seen in IRAS 09022–3615, where the OH119 emission feature is stronger than the absorption feature. This can be explained by either collisional excitation of the OH119 line or a relative excess of outflowing gas behind the continuum source. We consider the energetics inferred from the models of IRAS 09022–3615 and IRAS 19297–0406 relatively

uncertain, and we do not include them in the analysis below. As argued in Section 4.3, the molecular outflows in IRAS 05189–2524, Mrk 273, and IRAS 20551–4250, in which the OH119 absorption trough is relatively weak, are most likely relatively collimated.

## 6.2. Overall Continua and Energetics

The predictions for the galaxy continua are compared with the observed SEDs in Figure 17, showing that the total predicted continuum in most sources (in red) is similar to the observed SED in the transition from mid- to far-IR wavelengths (i.e., from 25 to ~50 µm). This total predicted continuum also includes the warm components with non-outflowing (or slowly outflowing) OH, so that the good fits of the modeled continua indicate that OH is an excellent tracer of the warm regions in ULIRGs with and without outflowing gas. The close agreement we find in most sources, however, fails in IRAS 09022–3615 and IRAS 13120–5453, where OH65 is not detected. In most sources, the far-IR emission at $\gtrsim 60$ µm is underpredicted, indicating the presence of cold dust unrelated to OH.

The mass outflow rate and momentum and energy fluxes, $\dot{M}_{\text{out}}$, $\dot{P}_{\text{out}}$, and $\dot{E}_{\text{out}}$, associated with a given component depend on four parameters (Equations (11)–(13)): the size $R$, the column density $N_H$, the gas velocity $v$, and the covering factor $f_c$. The dependence of $M_{\text{out}}$ on all these parameters is shown in Figure 20 for all individual components involved in the best-fit models of our sample. The individual components can be classified according to their predicted spatial extent. Components with $R < 350$ pc are considered “compact” and are
<table><caption>**Table 2**

Physical Parameters of Individual Outflowing Components</caption><thead><tr><th>Galaxy Name<br>(1)</th><th>T<br>(2)</th><th>V<br>(km s<sup>-1</sup>)</th><th>R<sub>1</sub><br>(pc)</th><th>R<br>(pc)</th><th>N<sub>H</sub><br>(10<sup>21</sup> cm<sup>-2</sup>)</th><th>f<sub>c</sub><br>(7)</th><th>M<sub>out</sub><br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>M&#x0307;<sub>out</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>P&#x0307;<sub>put</sub><br>(10<sup>34</sup> dyn)</th><th>E&#x0307;<sub>out</sub><br>(10<sup>41</sup> erg s<sup>-1</sup>)</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="2">IRAS 03158+4227</td><td>E</td><td>800</td><td>160</td><td>270</td><td>120</td><td>0.36</td><td>420</td><td>1300</td><td>660</td><td>2700</td></tr><tr><td>E</td><td>1200</td><td>260</td><td>400</td><td>8</td><td>0.41</td><td>75</td><td>230</td><td>170</td><td>1000</td></tr><tr><td rowspan="2">IRAS 05189-2524</td><td>E</td><td>200</td><td>77</td><td>170</td><td>120</td><td>0.21</td><td>100</td><td>120</td><td>16</td><td>16</td></tr><tr><td>E</td><td>550</td><td>130</td><td>340</td><td>40</td><td>0.14</td><td>89</td><td>150</td><td>52</td><td>140</td></tr><tr><td rowspan="2">IRAS 08572+3915</td><td>E</td><td>500</td><td>62</td><td>110</td><td>160</td><td>0.33</td><td>88</td><td>410</td><td>130</td><td>330</td></tr><tr><td>E</td><td>950</td><td>50</td><td>110</td><td>80</td><td>0.20</td><td>28</td><td>240</td><td>150</td><td>700</td></tr><tr><td rowspan="2">IRAS 10565+2448</td><td>E</td><td>250</td><td>340</td><td>480</td><td>10</td><td>0.51</td><td>170</td><td>88</td><td>14</td><td>17</td></tr><tr><td>E</td><td>420</td><td>240</td><td>400</td><td>20</td><td>0.34</td><td>160</td><td>160</td><td>43</td><td>90</td></tr><tr><td rowspan="3">Mrk 231</td><td>C</td><td>100</td><td>100</td><td>160</td><td>500</td><td>0.39</td><td>730</td><td>460</td><td>29</td><td>14</td></tr><tr><td>E</td><td>550</td><td>300</td><td>600</td><td>20</td><td>0.26</td><td>260</td><td>240</td><td>84</td><td>230</td></tr><tr><td>E</td><td>700</td><td>74</td><td>170</td><td>120</td><td>0.18</td><td>91</td><td>380</td><td>170</td><td>590</td></tr><tr><td rowspan="2">IRAS 13120-5453</td><td>E</td><td>420</td><td>280</td><td>470</td><td>16</td><td>0.36</td><td>180</td><td>160</td><td>43</td><td>89</td></tr><tr><td>E</td><td>300</td><td>67</td><td>130</td><td>80</td><td>0.28</td><td>51</td><td>120</td><td>23</td><td>35</td></tr><tr><td rowspan="3">Mrk 273</td><td>E</td><td>700</td><td>79</td><td>340</td><td>120</td><td>0.05</td><td>110</td><td>220</td><td>99</td><td>350</td></tr><tr><td>C</td><td>100</td><td>190</td><td>240</td><td>500</td><td>0.63</td><td>2500</td><td>1100</td><td>69</td><td>34</td></tr><tr><td>E</td><td>600</td><td>160</td><td>370</td><td>20</td><td>0.17</td><td>68</td><td>110</td><td>43</td><td>130</td></tr><tr><td rowspan="3">IRAS 14348-1447</td><td>E</td><td>300</td><td>180</td><td>340</td><td>80</td><td>0.27</td><td>350</td><td>310</td><td>60</td><td>90</td></tr><tr><td>E</td><td>600</td><td>170</td><td>270</td><td>12</td><td>0.39</td><td>47</td><td>110</td><td>41</td><td>120</td></tr><tr><td>E</td><td>250</td><td>110</td><td>240</td><td>40</td><td>0.22</td><td>71</td><td>76</td><td>12</td><td>15</td></tr><tr><td rowspan="2">IRAS 20100-4156</td><td>E</td><td>500</td><td>170</td><td>300</td><td>120</td><td>0.35</td><td>510</td><td>880</td><td>280</td><td>690</td></tr><tr><td>E</td><td>1100</td><td>190</td><td>300</td><td>12</td><td>0.40</td><td>58</td><td>220</td><td>150</td><td>850</td></tr><tr><td rowspan="2">IRAS 20551-4250</td><td>E</td><td>300</td><td>86</td><td>160</td><td>80</td><td>0.27</td><td>83</td><td>160</td><td>30</td><td>44</td></tr><tr><td>E</td><td>600</td><td>67</td><td>190</td><td>20</td><td>0.13</td><td>13</td><td>41</td><td>16</td><td>47</td></tr><tr><td rowspan="3">IRAS 23365+3604</td><td>C</td><td>100</td><td>70</td><td>100</td><td>500</td><td>0.46</td><td>340</td><td>340</td><td>22</td><td>11</td></tr><tr><td>E</td><td>350</td><td>73</td><td>190</td><td>60</td><td>0.14</td><td>45</td><td>82</td><td>18</td><td>32</td></tr><tr><td>E</td><td>470</td><td>260</td><td>450</td><td>20</td><td>0.33</td><td>190</td><td>200</td><td>59</td><td>140</td></tr></tbody></table>

**Note.** Columns: (1) Galaxy name. (2) Type of model. C: CORE; E: ENVELOPE. (3) Average velocity. (4) Minimum radius. (5) Adopted radius, after correcting for the covering factor $f_c$. (6) Hydrogen column density, adopting an abundance $X(OH) = 2.5 \times 10^{-6}$ relative to H nuclei. (7) Covering factor of the far-IR continuum by the outflowing OH, as defined in Section 5.5. (8) Outflowing gas mass. (9) Mass outflow rate. (10) Momentum flux. (11) Energy flux.

<table><caption>**Table 3**

Total Outflow Energetics Inferred from Models of OH</caption><thead><tr><th>Galaxy Name<br>(1)</th><th>M<sub>tot</sub> (HV)<br>(10<sup>7</sup> M<sub>☉</sub>)</th><th>M<sub>tot</sub> (T)<br>(10<sup>7</sup> M<sub>☉</sub>)</th><th>M&#x0307;<sub>tot</sub> (HV)<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>M&#x0307;<sub>tot</sub> (T)<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>P&#x0307;<sub>tot</sub> (HV)<br>(10<sup>34</sup> dyn)</th><th>P&#x0307;<sub>tot</sub> (T)<br>(10<sup>34</sup> dyn)</th><th>E&#x0307;<sub>tot</sub> (T)<br>(10<sup>42</sup> erg s<sup>-1</sup>)</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>IRAS 03158+4227</td><td>50<sup>+9</sup><sub>-24</sub></td><td>50<sup>+120</sup><sub>-24</sub></td><td>1500<sup>+340</sup><sub>-280</sub></td><td>1500<sup>+830</sup><sub>-280</sub></td><td>840<sup>+290</sup><sub>-140</sub></td><td>840<sup>+290</sup><sub>-140</sub></td><td>370<sup>+170</sup><sub>-81</sub></td></tr><tr><td>IRAS 05189-2524</td><td>19<sup>+27</sup><sub>-11</sub></td><td>...</td><td>270<sup>+22</sup><sub>-130</sub></td><td>...</td><td>67<sup>+13</sup><sub>-29</sub></td><td>...</td><td>16<sup>+4</sup><sub>-7</sub></td></tr><tr><td>IRAS 08572+3915</td><td>12<sup>+5</sup><sub>-8</sub></td><td>12<sup>+33</sup><sub>-7</sub></td><td>650<sup>+230</sup><sub>-340</sub></td><td>650<sup>+390</sup><sub>-320</sub></td><td>280<sup>+170</sup><sub>-110</sub></td><td>280<sup>+170</sup><sub>-99</sub></td><td>100<sup>+100</sup><sub>-33</sub></td></tr><tr><td>IRAS 10565+2448</td><td>32<sup>+13</sup><sub>-10</sub></td><td>32<sup>+28</sup><sub>-10</sub></td><td>250<sup>+25</sup><sub>-81</sub></td><td>250<sup>+150</sup><sub>-81</sub></td><td>57<sup>+13</sup><sub>-21</sub></td><td>57<sup>+13</sup><sub>-20</sub></td><td>11<sup>+4</sup><sub>-4</sub></td></tr><tr><td>Mrk 231</td><td>35<sup>+1</sup><sub>-14</sub></td><td>110<sup>+20</sup><sub>-8</sub></td><td>620<sup>+25</sup><sub>-260</sub></td><td>1100<sup>+41</sup><sub>-150</sub></td><td>250<sup>+36</sup><sub>-110</sub></td><td>280<sup>+38</sup><sub>-97</sub></td><td>83<sup>+23</sup><sub>-26</sub></td></tr><tr><td>IRAS 13120-5453</td><td>18<sup>+4</sup><sub>-4</sub></td><td>...</td><td>160<sup>+34</sup><sub>-33</sub></td><td>...</td><td>43<sup>+9</sup><sub>-9</sub></td><td>...</td><td>9<sup>+2</sup><sub>-2</sub></td></tr><tr><td>Mrk 273</td><td>16<sup>+2</sup><sub>-8</sub></td><td>16<sup>+60</sup><sub>-8</sub></td><td>350<sup>+35</sup><sub>-150</sub></td><td>350<sup>+490</sup><sub>-130</sub></td><td>120<sup>+10</sup><sub>-63</sub></td><td>120<sup>+44</sup><sub>-63</sub></td><td>38<sup>+10</sup><sub>-26</sub></td></tr><tr><td>IRAS 14348-1447</td><td>42<sup>+71</sup><sub>-29</sub></td><td>290<sup>+97</sup><sub>-35</sub></td><td>430<sup>+180</sup><sub>-280</sub></td><td>1500<sup>+270</sup><sub>-250</sub></td><td>100<sup>+97</sup><sub>-63</sub></td><td>170<sup>+97</sup><sub>-48</sub></td><td>25<sup>+32</sup><sub>-12</sub></td></tr><tr><td>IRAS 14378-3651</td><td>12<sup>+13</sup><sub>-1</sub></td><td>...</td><td>180<sup>+180</sup><sub>-33</sub></td><td>...</td><td>53<sup>+37</sup><sub>-14</sub></td><td>...</td><td>14<sup>+6</sup><sub>-5</sub></td></tr><tr><td>IRAS 20100-4156</td><td>56<sup>+22</sup><sub>-41</sub></td><td>...</td><td>1100<sup>+830</sup><sub>-620</sub></td><td>...</td><td>430<sup>+620</sup><sub>-250</sub></td><td>...</td><td>150<sup>+320</sup><sub>-80</sub></td></tr><tr><td>IRAS 20551-4250</td><td>10<sup>+11</sup><sub>-5</sub></td><td>...</td><td>200<sup>+210</sup><sub>-100</sub></td><td>...</td><td>45<sup>+81</sup><sub>-21</sub></td><td>...</td><td>9<sup>+29</sup><sub>-4</sub></td></tr><tr><td>IRAS 23365+3604</td><td>24<sup>+26</sup><sub>-13</sub></td><td>58<sup>+73</sup><sub>-47</sub></td><td>280<sup>+260</sup><sub>-160</sub></td><td>620<sup>+490</sup><sub>-430</sub></td><td>77<sup>+90</sup><sub>-43</sub></td><td>99<sup>+110</sup><sub>-48</sub></td><td>18<sup>+45</sup><sub>-8</sub></td></tr></tbody></table>

**Note.** Columns: (1) Galaxy name. (2)–(3) Outflow gas mass. (4)–(5) Mass outflow rate. (6)–(7) Momentum flux. (8) Energy flux. HV refers to the values calculated for only the high-velocity (>200 km s⁻¹) component(s), and the T columns include the values of the low-velocity components (only tabulated when the value or error is different from that of the corresponding HV column). The contribution of the low-velocity components to the energy flux is negligible, and only the total value is given.

indicated with circles in Figures 20, while those with R > 350 pc are considered “extended” and shown with squares. With typical radii of ~100 pc, all CORE components are compact, as well as most ENVELOPE models. The small radii in these models are a direct consequence of the need to

match the outflowing gas observed in high-lying OH84 and OH65 doublets, which require high $T_{dust} \gtrsim 70$ K for the associated continuum source to pump the excited OH levels, and hence small sizes to avoid overestimating the continuum emission. Conversely, the extended components reproduce the
<figure><img src="image_20.png" /><figcaption>Figure 17. Predictions for the continuum of the models for OH shown in Figures 13–16, compared with the observed SEDs. Blue, green, and magenta lines indicate the continua associated with the OH components shown with the same colors in Figures 13–16, and red is total. In some sources, the continuum associated with the OH components closely matches the observed SED between 20 and 50 µm, but the continuum is usually underpredicted at longer wavelengths.</figcaption></figure>

observed absorption and emission in the ground OH119 and OH79 doublets at velocities where no contribution is seen from the high-lying lines, with effective $T_{dust} \sim 55$ K.

The highest mass outflow rates $\dot{M}_{out} \gtrsim 300 M_{\odot} \text{yr}^{-1}$ are mostly found in compact components and can exhibit both low and high gas velocities, but they consistently show high column densities ($N_{H} \gtrsim 10^{23} \text{cm}^{-2}$). From Equations (11)–(13), we express

$$ 
\begin{aligned}
\dot{M}_{\text{out}} &= 320 \times \frac{f_c}{0.3} \times \frac{R}{150 \text{ pc}} \times \frac{N_H}{10^{23} \text{ cm}^{-2}} \\
& \quad \times \frac{v}{500 \text{ km s}^{-1}}, \\
\dot{P}_{\text{out}} &= 10^{36} \times \frac{f_c}{0.3} \times \frac{R}{150 \text{ pc}} \times \frac{N_H}{10^{23} \text{ cm}^{-2}} \\
& \quad \times \left( \frac{v}{500 \text{ km s}^{-1}} \right)^2, \\
\dot{E}_{\text{out}} &= 2.5 \times 10^{43} \times \frac{f_c}{0.3} \times \frac{R}{150 \text{ pc}} \times \frac{N_H}{10^{23} \text{ cm}^{-2}} \\
& \quad \times \left( \frac{v}{500 \text{ km s}^{-1}} \right)^3,
\end{aligned}
\quad (15) $$

in $M_{\odot} \text{yr}^{-1}$, dyn, and erg s$^{-1}$, respectively. The observational requirement for high column densities, $N_{OH} \gtrsim 2 \times 10^{17} \text{cm}^{-2}$, comes from the absorption in the high OH84 and OH65 doublets.

Very high values of $\dot{M}_{out}$ are also associated with moderate expansion velocities of the CORE components, although they are uncertain owing to the dependence on the adopted redshift. The components showing the highest momentum fluxes ($\dot{P}_{out} > 10^{36}$ dyn; blue symbols in Figure 20) are found in IRAS 03158+4227, IRAS 08572+3915, Mrk 231, and IRAS 20100–4156 and are associated with either high columns or extreme velocities ($>1000 \text{ km s}^{-1}$). In spite of their higher radii and overall high outflow velocities, most extended components have moderate $\dot{M}_{out} < 300 M_{\odot} \text{yr}^{-1}$, because values of $N_{OH} \lesssim 5 \times 10^{16} \text{cm}^{-2}$ are enough to account for the ground-state OH doublets.

The total $\dot{M}_{tot}$ estimates (Figure 18) show hints of a positive correlation with $L_{IR}$. Most sources with $L_{IR} \lesssim 2 \times 10^{12} L_{\odot}$ have $\dot{M}_{tot} \sim 200 M_{\odot} \text{yr}^{-1}$, with $\dot{M}_{tot}$ increasing sharply for $L_{IR} > 2 \times 10^{12} L_{\odot}$. We note that IRAS 08572+3915 is not necessarily an outlier, because its luminosity may be significantly higher than the apparent one according to the model by Efstathiou et al. (2014). $\dot{M}_{tot}$ appears also to be weakly correlated with $L_{AGN}$ and $L_{*}$, the latter indicating loading factors of 1–10.

With the exception of IRAS 14348–1447, the low-velocity components transport relatively low momentum in comparison with the high-velocity components. Starburst99 models (Leitherer et al. 1999) predict that a 40 Myr old continuous starburst can supply a maximum momentum (see also Veilleux et al. 2005; Heckman et al. 2015) of $\dot{P}_{*} \sim 3.5L_{*}/c$ (including radiation pressure at a level of $L_{*}/c$), which is marked as a shaded rectangle in Figure 18. Our inferred total momentum
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Estimated energetics of 12 out of 14 sources modeled in this paper; modeling results for IRAS 09022–3615 and IRAS 19297–0406 are relatively uncertain and not considered. $M_{\text{tot}}$, $\dot{P}_{\text{tot}}$, and $\dot{E}_{\text{tot}}$ are the total values defined in Equation (14) (see Section 5.5). Filled circles indicate the values obtained by ignoring the (uncertain) low-velocity ($<200 \text{ km s}^{-1}$) components, while open circles include them. Sources are labeled as follows: (1) IRAS 03158+4227, (2) IRAS 05189–2524, (3) IRAS 08572+3915, (4) IRAS 10565+2448, (5) Mrk 231, (6) IRAS 13120–5453, (7) Mrk 273, (8) IRAS 14348–1447, (9) IRAS 14378–3651, (10) IRAS 20100–4156, (11) IRAS 20551–4250, (12) IRAS 23365+3604. The dashed lines in the upper right panel indicate mass-loading factors of $\eta = \dot{M}_{\text{tot}}/\text{SFR} = 2$, 4, and 8. Shaded rectangles mark the momentum and energy rates that can be supplied by an AGN and a starburst according to several models discussed in Section 6.2.</figcaption></figure>

flux measurements are $\dot{P}_{\text{tot}} \gtrsim 3L_{*}/c$ (Figure 18), indicating that starbursts alone are most likely unable to drive the outflows via momentum-conserving winds, but may play an important role in some sources. The momentum supplied by the AGN may approach $\dot{P}_{\text{AGN}} \sim 2L_{\text{AGN}}/c$ (treating both radiation pressure on dust grains and AGN inner winds as $L_{\text{AGN}}/c$), also below the inferred total values. Nevertheless, Figures 19 (a)–(c) show the absolute values of $\dot{P}_{\text{tot}}$ as a function of $\dot{P}_{*}$, $\dot{P}_{\text{AGN}}$, and $\dot{P}_{*} + \dot{P}_{\text{AGN}}$, indicating a remarkable similarity between $\dot{P}_{\text{tot}}$ and the combined momentum rates supplied by the starburst and the AGN in most galaxies. The outflows in IRAS 03158+4227, IRAS 08572+3915, Mrk 231, and IRAS 20100–4156 (Figures 19) seem to require significant momentum boosts. These rough estimates are refined in Section 7.1, where consideration of the covering factors, the extended nature of the starburst, the potential well, and the momentum boost due to radiation trapping is made, and in which we show that momentum-driven winds may indeed be responsible for some outflow components with the combined effect of the AGN (in most cases dominant) and the starburst.

It is estimated that supernovae and stellar winds can provide a mechanical luminosity of up to ~1.8% of the starburst luminosity (Leitherer et al. 1999; Veilleux et al. 2005; Harrison et al. 2014), and only a fraction of the mechanical power will go into bulk motion of the ISM (1/4 according to Weaver et al. 1977). Energy-conserving winds from the starburst are then probably unable to drive the molecular outflows in Mrk 231, IRAS 08572+3915, IRAS 03158+4227, and IRAS 20100–4156, where mechanical luminosities of >1% of $L_{*}$ are inferred (Figure 18). For the rest of the sources, an expanding energy-conserving bubble generated by supernovae could in principle account for the observed outflows, with coupling efficiencies of ~10%–50%. We also refine these estimates in Section 7.2, where we find that only some low-velocity components can be driven in this way. Energy-conserving bubbles generated by AGN winds are believed to be able to supply a power of up to ~5% of $L_{\text{AGN}}$ (e.g., King & Pounds 2015), from which a fraction of ~1/2 may go into bulk motion of the ISM (Faucher-Giguère & Quataert 2012). Figure 18 shows that the energy flux is <5% for all sources, including the strong AGNs Mrk 231, IRAS 08572+3915,
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Top: inferred total momentum fluxes ($\dot{P}_{\text{tot}}$, Equation (14)) as a function of (a) the estimated momentum rate supplied by the starburst ($\dot{P}_{*} \sim 3.5L_{*}/c$), (b) the estimated momentum rate supplied by the AGN ($\dot{P}_{\text{AGN}} \sim 2L_{\text{AGN}}/c$), and (c) the sum of both. The gray diagonal lines indicate ratios of $\dot{P}_{\text{tot}}/\dot{P}_{*,\text{AGN,TOTAL}} = 1, 2$, and 10. Bottom: estimated total energy fluxes ($\dot{E}_{\text{tot}}$) as a function of (d) the starburst luminosity, (e) the AGN luminosity, and (f) the sum of both. The gray diagonal lines indicate ratios $\dot{E}_{\text{tot}}/L_{*,\text{AGN,TOTAL}} = 0.2\%, 1\%, \text{ and } 5\%$. Filled circles indicate the values obtained by ignoring the (uncertain) low-velocity ($<200 \text{ km s}^{-1}$) components, while open circles include them. Galaxies are labeled as in Figure 18.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Mass outflow rate, $\dot{M}_{\text{out}}$, of the individual components involved in the best-fit models, shown as a function of the radius $R$, velocity $v$, column density $N_H$, and covering factor $f_c$. Red symbols indicate CORE components; others indicate ENVELOPE components (see Figure 12). Blue symbols indicate components with high $\dot{P}_{\text{out}} > 1.0 \times 10^{36} \text{ dyn.}$ (a) The dashed vertical line marks the value we use for the classification of compact ($R < 350 \text{ pc}$) and extended ($R > 350 \text{ pc}$) components. (b)–(d) Circles and squares indicate compact and extended components, respectively.</figcaption></figure>

IRAS 03158+4227, and IRAS 20100–4156, in general agreement with theoretical predictions. We show in Section 7.2 that *partially* energy-conserving flows driven by the AGN are generally required for the high-velocity components.

The most powerful outflows (Figures 19(d)–(f)) are found in IRAS 03158+4227 (belonging to a widely separated pair), IRAS 08572+3915 and IRAS 20100–4156 (mergers with projected separation of ≈5.5 kpc), and Mrk 231 (post-merger),
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Comparison between the total properties of the outflows, $M_{\text{tot}}$, $\dot{M}_{\text{tot}}$, and $\dot{P}_{\text{tot}}$, derived from OH and from CO (1–0) (Cicone et al. 2014) for the sources that have been analyzed in both species. For OH, only the high-velocity ($>200 \text{ km s}^{-1}$) components are considered. The values of $M_{\text{tot}}$ and $\dot{P}_{\text{tot}}$ in Cicone et al. (2014) have been divided by 3 to match our definitions in Equations (11) and (12).</figcaption></figure>

apparently unrelated to merging stage. They are mostly associated with high AGN luminosities (Figure 19(e)). They are best identified by high-velocity absorption in the high-lying OH65 doublet and represent ~20% of all local ULIRGs observed in OH65. This suggests that stochastic and episodic strong AGN feedback events occur throughout the merger process.

### 6.3. Comparison with Other Tracers

Our calibration of the energetics from OH observations is based on (i) the reliability of the present models, and in particular of the inferred outflow sizes that are not directly observed but inferred from the OH excitation and the observed and predicted continuum; and (ii) our adopted OH abundance relative to H nuclei, $X_{\text{OH}} = 2.5 \times 10^{-6}$ (GA14; GA15), which in turn is based on both observations and chemical models. These assumptions are worth checking by comparing our values of $M_{\text{tot}}$, $\dot{M}_{\text{tot}}$, and $\dot{P}_{\text{tot}}$ (Equations (14)) with the values inferred from totally independent measurements of CO for the five sources where both OH and CO have been observed and analyzed: Mrk 231, IRAS 08572+3915, Mrk 273, IRAS 10565+2448, and IRAS 23365+3604 (Cicone et al. 2012, 2014). Since CO is not sensitive to low-velocity outflows, due to contamination by the core of the line, the values we include here from OH exclude the low-velocity ($<200 \text{ km s}^{-1}$) components. On the other hand, we divide the values of $\dot{M}_{\text{OF}}$ and $\dot{P}_{\text{OF}}$ in Cicone et al. (2012, 2014) by a factor of 3 to match our definitions given in Equations (11) and (12).

The outflowing masses, $M_{\text{out}}$, are compared in Figure 21(a). The agreement is within a factor of 2 in all sources except IRAS 08572+3915, where $M_{\text{out}}(\text{OH})$ is a factor of ~4 lower than $M_{\text{out}}(\text{CO})$. Since $M_{\text{out}}(\text{OH}) \propto R^2$ but $M_{\text{out}}(\text{CO})$ is independent of $R$, the discrepancy is a consequence of the very compact outflow (~100 pc) inferred for OH: both components (moderate and high velocities) are very excited, with OH84 as strong as OH119, as well as absorption in OH65 (Figure 7). Since the observed absorption is dominated by the high amounts of outflowing nuclear OH, and the source of far-IR absorption is expected to be compact, OH is not sensitive to the more extended outflowing gas in this source. On the redshifted side of the OH spectrum, in addition, IRAS 08572+3915 shows an emission feature in OH84 (and marginally also in OH65) at positive velocities that is not fitted (see Figure 1) and does not appear to belong to the same components we model in Figure 15. Nevertheless, the values

of $\dot{M}_{\text{tot}}$ and $\dot{P}_{\text{tot}}$ inferred from both species (Figures 21(b), (c)) are similar in this source, because the flow timescale, $t_{\text{flow}} \sim R/v$, is correspondingly higher in CO.

In all other sample sources, $\dot{M}_{\text{tot}}$ and $\dot{P}_{\text{tot}}$ are somewhat higher for OH than for CO. The largest discrepancies are found for IRAS 10565+2448 and IRAS 23365+3604, where we infer higher values by factors of ~2 and ~4 for OH, respectively. These discrepancies are a consequence of the extended emission estimated in CO, $R > 1$ kpc (Cicone et al. 2014), while we infer $R \sim 0.4–0.5$ kpc and $R \sim 0.2–0.45$ kpc for OH in IRAS 10565+2448 and IRAS 23365+3604, respectively. High-resolution CO observations in these sources will check whether there is a decrease in $\dot{M}$ with increasing distance to the center, as found in Mrk 231 (GA14; Feruglio et al. 2015). Nevertheless, the difference in the OH and CO spatial distributions may be real, as the maximum velocity observed in CO is lower than in OH in both sources (Figure 2).

It is also interesting to compare the energetics inferred from OH with those derived from Na I D (Rupke & Veilleux 2013a) for the four sources that have been analyzed in both species. Rupke & Veilleux (2013a) used the same expressions for the mass, momentum, and energy flux as we use in this study (Equations (7), (11)–(13)). The comparisons are shown in Figure 22, where the comparable values of the outflow masses inferred from the two species is remarkable. This is surprising because $M_{\text{out}} \propto R^2 N_{\text{H}}$; while our inferred values of $R$ are <1 kpc, the radii observed by Rupke & Veilleux (2013a) are 2–5 kpc, but the much higher columns inferred from OH compensate for this, thus resulting in similar total masses. The mass, momentum, and energy fluxes are then much higher for the molecular gas.

The similar outflow masses found for OH and for the more spatially extended CO and Na I D suggest that a significant fraction of the total outflow mass is loaded at small radii (a few $\times$ 100 pc), and the outflows are less efficient in evacuating the gas initially located in the extended kiloparsec-scale disks.

The total momentum fluxes ($\dot{P}_{\text{tot}}$) inferred for ULIRGs are compared in Figure 23 with the analogous quantities inferred by Heckman et al. (2015) from observations of UV ionized lines in a sample of starburst galaxies. Most of the composite ULIRGs lie close to the line $\dot{P}_{\text{tot}} = \dot{P}_{*}$, which is also a good fit for a number of pure starbursts, but usually lie above it, probably due to the AGN contribution. The comparison appears to indicate that nuclear starbursts also contribute to
<figure><img src="image_25.png" /><figcaption>Figure 22. Comparison between the properties of the outflows, $M_{tot}$, $\dot{M}_{tot}$, $\dot{P}_{tot}$, and $\dot{E}_{tot}$, derived from OH and from Na I D (Rupke & Veilleux 2013a) for the sources that have been analyzed in both species. For OH, only the high-velocity (>200 km s⁻¹) components are considered.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 23. Total momentum fluxes ($\dot{P}_{tot}$) of ULIRGs (red symbols, inferred from OH) and of starburst galaxies (black symbols, inferred from ionized UV lines in Heckman et al. 2015) plotted as a function of the momentum rate supplied by the starburst ($\dot{P}_{*} = 3.5 \times L_{*}/c$). The blue line indicates equal values in both axes.</figcaption></figure>

the molecular feedback in most ULIRGs, but are not enough to drive the observed outflows alone in most sources. As a corollary, if our estimate of the relative contributions of the AGN and the starburst is sufficiently accurate, it appears that feedback from various combinations of the AGN and the nuclear starburst regulates the star formation and the growth of the SMBH in gas-rich mergers. However, the four outlier ULIRGs labeled in Figure 23 appear to require very strong AGN contribution, likely in the form of energy-conserving winds. Comparison of the inferred molecular outflow properties with those of the ionized outflow phase in starburst galaxies (Heckman et al. 2015) is further explored in Appendix B.

Overall, the comparison of the energetics inferred from OH and CO, especially the mass outflow rate and the momentum flux, and the comparison between the observed SEDs and the model predictions in Figure 17 are satisfactory. We now push onward to the analysis of the molecular outflows by considering the individual components of our fits in the next section, as different components in a given source have different associated energetics.

## 7. Discussion

In our initial analysis of the momentum deposition rates supplied by the AGN and the starburst in Section 6.2 (Figure 19), we ignore the clumpiness of the molecular outflows, the momentum boost due to trapping of infrared radiation, and the finite extension of the starbursts, as well as the gravitational potential well in the central regions of the galaxies. We attempt to overcome these limitations with a simple model that accounts for these aspects, in the framework of momentum-driven and energy-driven outflows.

### 7.1. Momentum-driven Outflows

We attempt to evaluate the impact of radiation pressure on dust grains from both the AGN and the compact, but not point-like, central starburst, as well as the role of momentum deposition by the AGN winds, stellar winds, and supernovae in launching the molecular outflows, using a simple analytical model of momentum-driven flows based on recent work in the literature (Heckman et al. 2015; Ishibashi & Fabian 2015; Murray et al. 2005; Thompson et al. 2015). As pointed out by Stern et al. (2016), the observed momentum fluxes are the result of the time-averaged forces acting on the outflowing gas. It is therefore necessary to integrate the equation of motion beginning from the time that the outflowing shell is close to the energy source (and presumably covering a large fraction of the $4\pi$ sr), taking into account the change of the covering factor with increasing radius. The dynamics of the gas are given by

$$ \frac{dP_{\text{mod}}}{dt} = -\frac{M_{\text{out}}v_{\text{cir}}^2}{r} + \dot{P}_{*}(r) + \dot{P}_{\text{AGN}}(r), \quad (16) $$

where $P_{\text{mod}}$ is the (modeled) momentum of the outflow, and $v_{\text{cir}}$ is the “circular” velocity characterizing the potential well of an isothermal sphere (Heckman et al. 2015). We follow the approach of Murray et al. (2005) and use $v_{\text{cir}}^2 = 2\sigma^2$ for the sources where $\sigma$, the stellar velocity dispersion tabulated in Dasyra et al. (2006) and Genzel et al. (2001), is available. Mrk 231 has an unusually low value of $\sigma$ (Tacconi et al. 2002), and we adopt $v_{\text{cir}} = 240 \text{ km s}^{-1}$ to account for the dynamical mass of $6.7 \times 10^9 M_{\odot}$ within $r = 500 \text{ pc}$ (Davies et al. 2004). For the rest of the sources we adopt $v_{\text{cir}} = 205, 150, 220, \text{ and } 260 \text{ km s}^{-1}$ in IRAS 08572+3915, IRAS 10565+2448, IRAS 03158+4227, and IRAS 20100-4156, respectively, based on our fits to the OH65 doublet, and $v_{\text{cir}} = 226 \text{ km s}^{-1}$ in IRAS 13120-5453 based on recent molecular observations
with ALMA (Privon et al. 2016). Our $v_{\text{cir}}$ values are not corrected for rotation and flattening and are thus lower limits; nevertheless, the high $M_{\text{out}}$ values in the nuclear regions of ULIRGs (Table 2) suggest gas that is emanating from the rotating structures where the bulk of the nuclear gas mass is stored, with significant rotational support.

$\dot{P}_*(r)$ and $\dot{P}_{\text{AGN}}(r)$ are the rates of momentum deposition onto the ISM at radial position $r$ due to the starburst and the AGN, respectively. For the starburst,

$$ \dot{P}_*(r) = (1 + \tau_{\text{IR}} + \tau_{*,w}) f_c(r) \frac{L_{*,\text{encl}}(r)}{c}, \quad (17) $$

where $\tau_{\text{IR}}$ is the infrared opacity of the outflowing shell and accounts for the momentum boost due to absorption of re-radiated IR emission (Thompson et al. 2015). Although $\tau_{\text{IR}}$ is expected to vary with $r$, we avoid overestimating its effect (see Roth et al. 2012) by simply taking a characteristic constant value of $\mu m_H k_{\text{IR}} N_{\text{H}}$, where an absorption-mass coefficient of $k_{\text{IR}} = 5$ cm$^2$ g$^{-1}$ of gas is adopted (Ishibashi & Fabian 2015; Thompson et al. 2015). $\tau_{*,w}$ describes the contribution of stellar winds and supernova remnants, for which we adopt $\tau_{*,w} = 2.5$ based on Starburst99 (Leitherer et al. 1999; Veilleux et al. 2005; Heckman et al. 2015). $f_c(r)$ is the covering factor of the starburst feedback by the outflowing OH, for which we adopt two approaches: (i) The first is a “mixed” shell-cloud approach, in which a shell with $f_c = 1$ is used for $r \le R_1$ (see Equation (6)), while $f_c = (R_1/r)^2$ for $r > R_1$ owing to shell breakup and geometrical dilution of the cloud ensemble. We thus neglect here the expansion of clouds included in Thompson et al. (2015) because the high-density molecular tracers HCN and HCO$^+$ are detected in the outflow of Mrk 231 at large distances from the center (Aalto et al. 2012, 2015; Lindberg et al. 2016). (ii) An “aligned-momenta” approach is also considered in which $f_c = 1$ at all radii, reflecting the possibility that the AGN and starburst momenta are preferentially released in the same directions as the clumpy, outflowing OH. $L_{*,\text{encl}}(r)$ is the starburst luminosity enclosed within radius $r$. While in Mrk 231 Davies et al. (2004) fitted the spatial distribution of the nuclear starburst intensity with an exponential $I(r) = I_0 \exp\{-r/r_d\}$, where $r_d = 150-200$ pc, we have approximated the corresponding enclosed luminosity as a linear function of $r$,

$$ L_{*,\text{encl}}(r) = \frac{r}{3r_d} L_*, \quad r \le 3 r_d, \quad (18) $$

where $L_*$ is the total starburst luminosity, and applied this expression to all galaxies with $r_d = 150$ pc. We thus assume a very compact starburst to avoid underestimating its momentum deposition.<fnref n="16" /> Since in Equation (16) we also neglect the inward force exerted by the stars formed outside $r$, the role of the nuclear starburst in launching the outflows is most likely not underestimated.

The momentum deposition due to AGN feedback is described as

$$ \dot{P}_{\text{AGN}}(r) = (1 + \tau_{\text{IR}} + \tau_{\text{AGN},w}) f_c(r) \frac{L_{\text{AGN}}}{c}, \quad (19) $$

with $\tau_{\text{AGN},w} = 1$ (e.g., Zubovas & King 2012). We implicitly assume that, after being scattered by electrons in a Compton-thick medium, most photons peak in the UV and are absorbed by dust, thus contributing at least twice (for $f_c = 1$) to the momentum deposition onto the ISM. Comparison of this simple momentum-driven model with observations enables us to evaluate the role of energy-driven flows (e.g., DeBuhr et al. 2012; Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Feruglio et al. 2015; Tombesi et al. 2015; Stern et al. 2016) in the OH outflows, which is considered below in Section 7.2.

In the following we ignore different terms in Equations (17) and (19) to isolate the relative roles of radiation pressure on dust grains and winds. The equation of motion is then integrated from the “launching” radius ($R_0$; see below) to the “observed” radius ($R$, Table 2), by assuming a constant outflowing gas mass, so that the first member in Equation (16) is $0.5 M_{\text{out}} dv^2/dr$. By defining the modeled momentum flux in the same way as for the observations, i.e., $\dot{P}_{\text{mod}} = M_{\text{out}} v^2/R$ (Equation (12)), the model predictions can be directly compared with our values for the observed momentum rates inferred for the individual components, $\dot{P}_{\text{out}}$ (Section 5.5, Table 2). We also rely on the AGN fractions ($\alpha_{\text{AGN}} = L_{\text{AGN}}/L_{\text{bol}}$; Table 1) given in Veilleux et al. (2013), some of them being uncertain owing to high extinction even at far-IR wavelengths (GA15).

### 7.1.1. Can Radiation Pressure on Dust Grains Drive the Outflows?

First, we ignore the terms proportional to $\tau_{*,w}$ and $\tau_{\text{AGN},w}$ in Equations (17) and (19), thus taking into account only radiation pressure from both the AGN and the starburst. By also neglecting the gravitational term (superscript NG) and assuming that the outflow is launched at very small radii (taking $R_0 = 0$), a strong upper limit is found for the momentum flux that can be attained as a result of radiation pressure in the shell-cloud approach:

$$ \dot{P}_{\text{mod,RP}}^{\text{NG}}(R) = \frac{R_1}{R} \dot{P}_*(R_1) \left[ 1 + 2 \ln \left( \frac{R}{R_1} \right) \right] + \frac{2R_1}{R} \dot{P}_{\text{AGN}}(R_1) \left( 2 - \frac{R_1}{R} \right), \quad (20) $$

where $\dot{P}_*(R_1) = (1 + \tau_{\text{IR}}) L_* R_1 / (3 c t_d)$ and $\dot{P}_{\text{AGN}}(R_1) = (1 + \tau_{\text{IR}}) L_{\text{AGN}} / c$ only include the radiative terms. We compare in Figure 24(a) the observed momentum fluxes $\dot{P}_{\text{out}}$ with the strong upper limits of Equation (20), showing that the observed values are above the limiting values for most high-velocity components.

By now including the gravitational term, we obtain the critical momentum flux $\dot{P}_{\text{crit}}(r) = M_{\text{out}} v_{\text{cir}}^2 / r$ that must be overcome at some position $r$ by $\dot{P}_*(r) + \dot{P}_{\text{AGN}}(r)$ in order to allow the outflow launching to proceed (Heckman et al. 2015; Thompson et al. 2015). In the shell-cloud approach, this condition should be met at least at $r = R_1$, where $(\dot{P}_*(r) + \dot{P}_{\text{AGN}}(r)) / \dot{P}_{\text{crit}}(r)$ is highest, and thus we compare in Figure 24(b) $\dot{P}_{\text{RP}}(R_1) \equiv \dot{P}_*(R_1) + \dot{P}_{\text{AGN}}(R_1)$ and $\dot{P}_{\text{crit}}(R_1)$. Interestingly, both values are found to be similar in most sources within the uncertainties, indicating that, together with rotation, radiation pressure (dominated by the AGN in most cases) is able to support the structures against gravity, particularly in the direction of the rotation axis (see also Murray et al. 2005; GA15). This support is further explored in Appendix C.

<footnote n="16">¹⁶ The implied stellar luminosity-to-mass ratio within $3r_d$ is $L_{*,\text{encl}}(r)/M_{*,\text{encl}}(r) = G L_*/(3r_d v_{\text{cir}}^2)$, which gives $240 L_{\odot}/M_{\odot}$ for $L_* = 10^{12} L_{\odot}$ and $v_{\text{cir}} = 200$ km s$^{-1}$.</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 24. Evaluation of the possibility that radiation pressure from the AGN and the starburst can drive the observed molecular outflows with the shell-cloud approach (Section 7.1.1). Circles and triangles indicate components with $\tau_{IR} \ge 0.5$ and $<0.5$, respectively, thus distinguishing components with significant momentum boost owing to radiation trapping. (a) Ratio of the observed momentum rates for the individual components, $P_{\text{out}}$, and the values predicted by Equation (20), which ignores the effect of the potential well. (b) Ratio of the rates of momentum deposition due to radiation pressure, $\dot{P}_{\text{RP}}$, to the critical momentum flux $P_{\text{crit}}$ that needs to be overcome in order to launch the outflow. Both values are evaluated at position $R_1$, where the ratio is expected to be maximum. (c) Ratio of the observed momentum rates to the values predicted by Equation (21) that include the gravitational term (only those components that satisfy $\dot{P}_{\text{RP}}(R_1) > \dot{P}_{\text{crit}}(R_1)$ have solutions). Galaxies are labeled as in Figure 18.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25. Evaluation of the possibility that momentum-driven winds and radiation pressure from the AGN and the starburst can drive the observed molecular outflows with the shell-cloud approach (Section 7.1.2). Circles and triangles indicate components with $\tau_{IR} \ge 0.5$ and $<0.5$, respectively, thus indicating the components with significant momentum boost due to radiation trapping. (a) Ratio of the observed momentum rates for the individual components, $P_{\text{out}}$, and the predicted values when the potential well is ignored. (b) Ratio of the rates of momentum deposition due to radiation pressure and winds, $\dot{P}_{\text{T}}$, and the critical momentum flux $P_{\text{crit}}$ that must be overcome in order to launch the outflow. Both values are evaluated at position $R_1$, where the ratio is expected to be maximum. (c) Ratio of the observed momentum rates and the values predicted by including the gravitational term (only those components that satisfy $\dot{P}_{\text{T}}(R_1) > \dot{P}_{\text{crit}}(R_1)$ have solutions). Galaxies are labeled as in Figure 18.</figcaption></figure>

For the model components that satisfy $\dot{P}_{\text{RP}}(R_1) > \dot{P}_{\text{crit}}(R_1)$, we calculate the expected momentum flux at the observed position R:

$$ \dot{P}_{\text{mod,RP}}(R) = -\frac{2R_1}{R}\dot{P}_{\text{crit}}(R_1)\ln\frac{R}{R_0} + \frac{R_1}{R}\dot{P}_{\text{AGN}}(R_1)\left[1 - \left(\frac{R_0}{R_1}\right)^2 + 2\ln\frac{R}{R_1}\right] + \frac{2R_1}{R}\dot{P}_{\text{AGN}}(R_1)\left[2 - \frac{R_0}{R_1} - \frac{R_1}{R}\right], \quad (21) $$

where $R_0$ is the launching radius ($\dot{P}_{\text{RP}}(R_0) = \dot{P}_{\text{crit}}(R_0)$). The resulting values are compared with the observed momentum rates in Figure 24(c), showing that no outflowing component can be explained with this model. We conclude that radiation

pressure alone can provide (partial) support against gravitation in the vertical direction (see also Thompson et al. 2005), but can hardly drive the observed outflows.

### 7.1.2. Inclusion of Winds

We now include the terms proportional to $\tau_{\text{AGN,w}}$ and $\tau_{\text{AGN,w}}$ in Equations (17) and (19) and follow the same three steps as in Section 7.1.1. In the shell-cloud approach, the relevant equations have the same forms as Equations (20) and (21). Ignoring the gravitational term, Figure 25(a) shows that the observed momentum fluxes are mostly lower than or similar to the predictions. This scenario may describe gas that is rotationally supported and radially pushed, flowing close to the equatorial plane of the disk. Because of this extra support, molecular outflows coplanar with the disk can be roughly
accounted for via momentum deposition by the AGN and the starburst, with no need for energy-conserving phases.

If rotational support is unimportant (e.g., if the outflow is driven in the polar direction), the total momentum deposited by the AGN and the starburst, $\dot{P}_T(R_1) \equiv \dot{P}_*(R_1) + \dot{P}_{AGN}(R_1)$, must be higher than the critical value $\dot{P}_{crit}(R_1)$ needed to launch the outflows. Figure 25(b) compares $\dot{P}_T(R_1)$ to $\dot{P}_{crit}(R_1)$, showing that the potential well can be overcome for most outflow components once winds are included. The calculated momentum fluxes at the effective observed position R are usually lower than those inferred from observations (Figure 25(c)), but at least one velocity component in IRAS 03158+4227, IRAS 05189–2524, IRAS 14348–1447, IRAS 14378–3651, Mrk 231, IRAS 20100–4156, IRAS 20551–4250, and IRAS 23365+3604 and possibly the two components in IRAS 10565+2448 appear to be outflows with little momentum boost (≤3).

The aligned-momenta approach, simulating the scenario in which winds from the AGN and starburst are released in the same directions as the molecular clumps, provides a much more efficient coupling between the released and outflow momenta. As shown in Figure 26, most modeled and observed momentum fluxes are similar in this scenario. The potential well is included in these calculations. Among the high-velocity components, only the components with high column densities in IRAS 03158+4227, IRAS 08572+3915, and IRAS 20100–4156 have momentum boosts higher than 3. The momentum fluxes of the 100 km s⁻¹, very massive components in IRAS 14348–1447 and IRAS 23365+3604 cannot be explained in this way; results for these components are very sensitive to $v_{cir}$ (see Figures 24(a) and 25(a)).

## 7.2. Energy-driven Outflows

Both starburst and AGN winds can create an adiabatic hot bubble that would drive a vigorous outflow. Since radiation pressure, together with rotation, can support the structures against gravity (Section 7.1.1), here we drop both the gravitational and radiation pressure terms and use for simplicity $R_0 = 0$. Starting with starburst winds (both stellar winds and supernovae), the energy injection in the outflowing ISM is

$$ \frac{dE_{\text{mod}}}{dt} = \beta_{\text{out,*}} f_c(r) K_w L_{*,\text{encl}}(r), \quad (22) $$

where $K_w L_*$ is the power of the winds and $\beta_{out,*}$ is the fraction of this energy that goes into bulk motion of the shocked ISM. We adopt $K_w = 0.02$ (Leitherer et al. 1999; Veilleux et al. 2005; Harrison et al. 2014) and $\beta_{out,*} = 1/4$ (Weaver et al. 1977). Integrating Equation (22) gives for the momentum flux in the shell-cloud approach

$$ \dot{P}_{\text{mod,*}}^{\text{E}} = \frac{\beta_{\text{out,*}} K_w L_* R_1}{2 r_d v} \times \frac{R_1}{R} \left( 1 + 2 \ln \frac{R}{R_1} \right), \quad (23) $$

where v is the velocity of the outflowing gas. Figure 27(a) compares the observed values of $\dot{P}_{out}$ with the values given by Equation (23) for all individual components. Some components could be explained by feedback from a compact starburst, including the ~100 km s⁻¹ outflows. However, most high-velocity components show momentum fluxes that exceed the calculated upper limit. The reason for this behavior is that the momentum boost, relative to the momentum-driven flows, is

<figure><img src="image_29.png" /><figcaption>Figure 26. Same as Figure 25(c), but for the aligned-momenta approach in which $f_c = 1$ is assumed at all radii. The potential well is included in these calculations. Galaxies are labeled as in Figure 18.</figcaption></figure>

given by $\beta_{out,*} \times 3000 \text{ km s}^{-1}/v$ (≤7.5 in all cases), and thus little boost is expected for the high-velocity components, due to the moderate velocities of the starburst winds (a few × 10³ km s⁻¹). This result suggests that bubbles generated by AGN winds are required for the high-velocity components.

For the AGN winds, we solve

$$ \frac{dE_{\text{mod}}}{dt} = \beta_{\text{out,AGN}} f_c(r) \dot{E}_{\text{AGN}}, \quad (24) $$

where $\beta_{out,AGN}$ is the fraction of the injected energy ($\dot{E}_{AGN}$) that goes into bulk motion of the shocked ISM. Following Faucher-Giguère & Quataert (2012) and Stern et al. (2016), we adopt $\beta_{out,AGN} = 1/2$ and integrate Equation (24) to give

$$ \dot{P}_{\text{mod,AGN}}^{\text{E}} = g(R) \frac{L_{\text{AGN}}}{c} \beta_{\text{out,AGN}} \frac{v_{\text{in}}}{v}, \quad (25) $$

where $g(R) = 3R_1/(2R) \times (2 - R_1/R)$ is of order unity, and the momentum of the inner wind is written as $L_{AGN}/c$ and has a velocity $v_{in}$. Figure 27(b) shows $\dot{P}_{out}/(g L_{AGN}/c)$, together with the solid lines indicative of the $\dot{P} \propto v^{-1}$ relationship expected for an energy-conserving outflow (e.g., Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Stern et al. 2016; Tombesi et al. 2015), with assumed nuclear wind speeds of $v_{in} = 0.1c$ and 0.033c. The required momentum boosts are found to be ≲10 in most sources and ≲20 in all cases, in general consistent with $v_{in} \lesssim 0.1c$. We also find that the expected trend of increasing momentum boost with decreasing velocity no longer holds in this regime of relatively low velocities (100–1000 km s⁻¹).

## 7.3. Interpretation of these Models

Stern et al. (2016) argue that the pressure of the hot shocked wind can be constrained by emission-line ratios measured in
<figure><img src="image_30.png" /><figcaption>**Figure 27.** Predictions for energy-driven outflows, which use a shell-cloud approach and assume that gravity and radiation pressure balance and thus are both ignored (Section 7.2). (a) Ratio of the observed momentum rates for the individual components, $\dot{P}_{\text{out}}$, and the predicted maximum values for starburst winds (Equation (23)). (b) Ratio of the observed momentum rates and $gL_{\text{AGN}}/c$, where g is of order unity (Equation (25)). In energy-driven flows, this ratio is expected to be $\beta_{\text{out,AGN}} v_{\text{in}}/v$, where $v_{\text{in}}$ is the velocity of the inner AGN wind and $\beta_{\text{out,AGN}}$ is the fraction of the injected energy that goes into bulk motion of the ISM gas. The two solid lines show the expected relationship for $v_{\text{in}} = 0.1c$ and $0.033c$, with $\beta_{\text{out,AGN}} = 1/2$. Filled symbols indicate components that, within a factor of 3, can be explained neither through momentum driven by the starburst and AGN nor through energy driven by the starburst, and thus require an energy-conserving phase driven by the AGN. Galaxies are labeled as in Figure 18.</figcaption></figure>

the colder line-emitting gas. They showed that in quasars radiation pressure is probably dominating at all scales, with an upper limit $P_{\text{hot}}/P_{\text{rad}} \le 6$ for the spatial domain we sample with OH. Stern et al. (2016) also indicated that high momentum boosts could only be obtained in buried quasars where the shocked wind was tightly confined, and the high momentum boosts reported for quasars in the literature could thus only be obtained early in their evolution, during the buried stage. This conclusion is consistent with Zubovas & Nayakshin (2014), who found that an energy-conserving bubble with escape channels for the hot gas cools adiabatically and hence is less efficient in driving the ISM, leading to small momentum boosts that approach the momentum-conserving scenario.

GA15 argued that OH65 is an excellent tool to probe this buried stage. All sources in the present sample except IRAS 09022–3615, IRAS 13120–5453, and IRAS 10565 +2448 show strong OH65 features. Figure 4 shows that the highest OH velocities are indeed found in sources with strong OH65 features, consistent with the model described by Zubovas & Nayakshin (2014) and Stern et al. (2016). However, some sources with very strong OH65 absorption only show low-velocity outflows, or even lack evidence of outflows (region III in Figure 4), suggesting that they are still in an early stage of this process, and the columns around some extremely obscured sources are so high that they cannot (yet) be accelerated to high velocities and dispersed. Therefore, perhaps it is not surprising that high-velocity molecular outflows are found in moderately (but not extremely) obscured sources, in which the shocked wind is moderately confined but finds paths with columns $\lesssim 10^{23} \text{ cm}^{-2}$ that can be efficiently accelerated. Once these paths are found and the wind emerges, its pressure will decrease and moderate momentum boosts are expected, as found here in most ULIRGs.

Our results for the momentum boost are, however, sensitive to the way in which the gravitational potential well and the geometrical dilution of the outflowing gas are treated (Figures 25(a), (c) and 26). If the AGN and starburst winds are not isotropic but restricted to the solid angles of the outflowing OH (the aligned momentum approach discussed in Section 7.1.2.), no geometrical loss of momentum is expected and most outflows approach the momentum-conserving scenario. On the other hand, rotationally supported structures (i.e., ignoring the potential well; Section 7.1.2) are obviously more easily launched and accelerated, and a momentum-conserving approach would be enough to account for most of the inferred momenta (Figure 25(a)). However, this only applies to the outflow components that are launched close to the equatorial plane of the nuclear disk or torus around the AGN. Because we require high columns to account for the molecular outflows, which are observed against strong sources of far-IR radiation (and hence against regions with high columns), at least some low-velocity components can indeed be associated with “equatorial” flows. This may not be the case for most high-velocity components, as recent CO interferometric observations show that the high-velocity gas is not coplanar with the rotating disk in NGC 1614 and IRAS 17208–0014 (García-Burillo et al. 2015) and in Mrk 231 (Feruglio et al. 2015). Nevertheless, the OH outflows are wide-angle, and the possibility that some low-velocity components flow close to the equatorial plane in a momentum-driven phase is again consistent with the theoretical approach by Zubovas & Nayakshin (2014), who found that in the directions of higher resistance, the ISM is mainly affected by the incident momentum of the inner outflow, while the outflow energy tends to escape along leaky paths. This makes the deposition of AGN energy in the ISM (feedback) less efficient, potentially explaining the large masses that SMBHs attain and, likewise,
<figure><img src="image_31.png" /><figcaption>**Figure 28.** (a) Outflow depletion time (filled circles) calculated from the total molecular mass ($M_{H2}$) of the galaxy and the total values of the mass outflow rates, $t_{dep} = M_{H2}/\dot{M}_{tot}$. Most $M_{H2}$ values have been inferred from the CO (1–0) emission, (with references in GA15 and A. Gowardhan et al. 2016, in preparation); in the case of IRAS 13120–5453, $M_{H2} = 7.5 \times 10^9 M_{\odot}$ in the central 1 kpc diameter region is estimated from the continuum associated with the H$_2$O submillimeter emission (Privon et al. 2016). The values of $\dot{M}_{tot}$ (Equation (14)) only include the high-velocity components ($>200$ km s$^{-1}$). Asterisks show the gas consumption timescales, $M_{H2}$/SFR, where SFR($M_{\odot}$ yr$^{-1}$) = $10^{-10}L_{*}(L_{\odot})$. (b) Flowing times for the individual outflowing components, calculated as $t_{flow} = r/v$. Open and filled squares indicate low-velocity ($<300$ km s$^{-1}$) and high-velocity ($>300$ km s$^{-1}$) components, respectively. (c) Fraction of the total molecular gas of the galaxy involved in the observed outflows. Open (filled) squares include (do not include) low-velocity ($<200$ km s$^{-1}$) components. Error bars in panels (a) and (c) only include uncertainties in $\dot{M}_{tot}$ and $M_{tot}$, respectively.</figcaption></figure>

the $M_{BH}-\sigma$ relation (Zubovas & Nayakshin 2014; King & Pounds 2015).

In some ULIRGs, energy-conserving phases with momentum boosts in the range of ~3–20 are most likely required to explain some outflowing components (Figures 25(c), 26, and 27(b)). The corresponding sources are IRAS 03158+4227, IRAS 08572+3915, IRAS 20100–4156, and possibly Mrk 231 and Mrk 273, although the momentum boost in the

<figure><img src="image_32.png" /><figcaption>**Figure 29.** Flow times for the individual outflow components (see Figure 28) as a function of the equivalent width of the OH65 doublet (between ~200 and +200 km s$^{-1}$; GA15). Open and filled squares indicate low-velocity ($<300$ km s$^{-1}$) and high-velocity ($>300$ km s$^{-1}$) components, respectively. Recent outflow events are found in active sources as diagnosed by the absorption in the high-lying OH65 doublet.</figcaption></figure>

last two objects is small if $f_c = 1$ is assumed (Figure 26). Most of these galaxies have strong absorption in OH65 at blueshifted velocities. Our model for Mrk 273, however, is relatively uncertain because of the weakness of the high-velocity (~700 km s$^{-1}$) spectral feature in OH119, OH79, and OH84, which suggests a jet-driven collimated outflow with energetics unrelated to the source luminosity. Momentum boosts of <20 are lower than the values obtained in theoretical approaches to energy-conserving outflows for an inner wind velocity of $v_{in} = 0.1c$ (e.g., Faucher-Giguère & Quataert 2012; King & Pounds 2015). Since this may be indicative of only partial covering and confinement of the hot shocked gas (Stern et al. 2016) even in optically thick ULIRGs, we denote this scenario as “partially energy-conserving outflows.”

## 7.4. Depletion Time and Flow Timescales

Figure 28(a) shows the depletion timescales ($t_{dep}$) calculated from the ratios of total molecular gas masses ($M_{H2}$) to total mass-loss rates $\dot{M}_{tot}$. The ratios lie mostly between $10^7$ and $10^8$ yr, except in the powerful IRAS 03158+4227, IRAS 08572+3915, and IRAS 20100–4156. These timescales are significantly shorter than the gas consumption timescales due to star formation ($t_{con}$; asterisks in Figure 28(a); see also Evans et al. 2002), and they are much shorter than the strangulation timescales of local galaxies (Peng et al. 2015). The values of $M_{H2}$ used to estimate both $t_{dep}$ and $t_{con}$, as inferred in most cases from simple scaling values of CO (1–0) luminosities (GA15), should be considered with caution because the fraction of the CO population in the lowest rotational levels is tiny; values of the mass-to-CO luminosity ratio a factor of ~2 higher are not excluded (Solomon et al. 1997). In addition, the mass stored in the neutral (atomic) and ionized gas components is not taken into account. The gas masses may be some factor higher than measured in CO, but an upper limit of $t_{dep} \lesssim 10^8$ yr can be established within uncertainties. The $t_{con}/t_{dep} = 1.1–15$ ratio is, however, independent of gas mass and indicates that the outflows mainly limit the amount of gas that will be converted into stars, i.e., the outflows are mainly responsible for the quenching of star formation. The flow timescales, defined as the time required for the outflows to arrive at the current (inferred)
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Depletion timescales (Figure 28(a)) as a function of $L_{IR}$, $L_{AGN}$, and $L_{*}$. The source with $t_{dep} \sim 2 \times 10^6$ yr is IRAS 08572+3915, for which $L_{AGN}$ (and hence $L_{IR}$) is probably underestimated (Efstathiou et al. 2014). The best anticorrelation is found between $t_{dep}$ and $L_{AGN}$ (correlation coefficient of -0.59, and -0.74 if IRAS 08572+3915 is excluded), indicating that the feedback from the AGN dominates the quenching of star formation at least for $L_{AGN} \gtrsim 10^{12} L_{\odot}$.</figcaption></figure>

radial position, are shown in Figure 28(b) for the individual components, indicating that they are much shorter than $t_{dep}$. Some low-velocity components ($<300$ km s$^{-1}$; open symbols) have $t_{flow}$ comparable to that of the high-velocity components in the same source, suggesting that they may be a consequence of the same AGN “event.” The fraction of the total molecular mass carried by the high-velocity outflows (Figure 28(c)) is $\lesssim$10%, so that the current “cycle” would have to be repeated a dozen times to deplete the galaxy’s molecular gas. We show in Figure 29, however, that sources with weak absorption in OH65 have $t_{flow} \gtrsim 1$ Myr, while shorter $t_{flow}$ are only found in sources with strong absorption in OH65, i.e., the regions that are expected to represent the most buried and active stage of nuclear starburst–AGN coevolution (GA15). The two sources with weak OH65 absorption, IRAS 13120–5453 and IRAS 10565+2448 (see also Figure 4), show outflowing signatures that may be a vestige of past strong nuclear activity, which has subsided over time as a result of the outflowing activity. This suggests that the molecular outflows are efficient at evacuating the nuclear regions of ULIRGs, but less efficient at expelling the gas located at larger distances.

## 7.5. Negative Feedback: The Relative Roles of the AGN and the Starburst

The outflow depletion timescale is plotted as a function of $L_{IR}$, $L_{AGN}$, and $L_{*}$ in Figure 30. Since $t_{dep}$ is sensitive to both the cumulative effect (the $M_{H2}$ reservoir decreases as the molecular gas is expelled) and the present effect (through $M_{tot}$) of the negative feedback, Figure 30 allows us to check the relative roles of the AGN and the starburst in quenching the star formation in the sources. The middle and right panels indicate that $t_{dep}$ is better correlated with $L_{AGN}$ than with $L_{*}$, indicating that the AGN is primarily responsible for the short $t_{dep} \lesssim 10^7$ yr found in IRAS 08572+3915, IRAS 03158+4227, IRAS 20100–4156, Mrk 231, IRAS 05189–2524, and most probably Mrk 273. The $t_{dep}-L_{AGN}$ correlation would be even tighter if IRAS 08572+3915 had an AGN luminosity well in excess of $10^{12} L_{\odot}$, as was argued by Efstathiou et al.

(2014). The $L_{AGN}$ threshold for AGN-dominated feedback, $\sim 10^{12} L_{\odot}$, is similar to that found by Veilleux et al. (2013), $10^{11.8} L_{\odot}$, from the analysis of the outflow velocities inferred from the OH119 doublet in a larger sample of (U)LIRGs.

For the rest of the sources, the starburst may have a significant role in the observed feedback, but the relative roles of the AGN and the starburst are more uncertain. According to the models developed in Sections 7.1 and 7.2, and taking at face value the $\alpha_{AGN}$ values in Table 1, the AGN appears to dominate in IRAS 20551–4250 and IRAS 23365+3604, and the starburst may have a comparable contribution in IRAS 10565+2448 and a dominant role in IRAS 13120–5453, IRAS 14348–1447, and IRAS 14378–3651. However, the main caveats for these assignments are as follows: (i) the extended outflow components may be the relic of an epoch when the AGN was very active, but has since subsided; (ii) $L_{AGN}$ may be underestimated in sources with strong OH65 absorption, as in IRAS 14348–1447; (iii) the presence of high-velocity wings in the OH119 doublet in all these sources (with velocities approaching or exceeding 800 km s$^{-1}$; see Figure 13) may be difficult to explain through supernova explosions (Martizzi et al. 2016), perhaps indicating the dominance of the AGN in driving at least this extremely high velocity molecular gas.

# 8. Conclusions

## 8.1. Observational Results

The main observational results of our analysis of the line intensities and profiles of four OH doublets (at 119, 79, 84, and 65 µm) in a sample of 14 (U)LIRGs observed with Herschel/PACS are as follows:

1. In ULIRGs in which a P Cygni profile or a high-velocity blueshifted absorption wing is detected in the ground-state OH119 transition, we confirm that the more optically thin OH79 transition also shows clear indication of outflowing gas in its line shape, though at lower velocities (see also Sturm et al. 2011).
2. OH119 is optically thick even in the line wings, and its strength provides a measure of the covering factor of the 119 µm continuum by the outflowing OH. The overall covering factor ranges from 0.1 to >0.5 (Section 4.3).

3. The highest outflowing velocities arise in sources with strong absorption in the high-lying OH65 doublet at central velocities (Section 4.4), indicating that the buried nuclear regions of ULIRGs are most efficient in launching and accelerating the molecular outflows.

4. From the velocity shifts observed in OH65 and their correspondence with those of OH79, we find evidence for low-velocity expansion of the nuclear regions in several ULIRGs, mostly in those sources that also present P Cygni line shapes in OH119. The peak absorption velocity in the high-lying lines is usually closer to systemic than in the ground-state OH119 and OH79 doublets (Section 4.5), indicating increasing column densities for lower outflow velocities.

5. Spatially compact and extended outflows are best distinguished by the OH84/OH119 ratio at blueshifted velocities (Section 4.7).

## 8.2. Modeling Results

We have modeled the observed OH doublet profiles in 14 ULIRGs through $\chi^2$ minimization of combinations of single-component models (Section 5.1). The models yield satisfactory fits to the observed OH doublets in 12 sources, while closely matching the observed SEDs in the transition from mid- to far-IR wavelengths (25–50 µm; Figures 13–17). The models enable us to estimate the energetics associated with the individual components and with each source globally (Section 5.5).

1. Excluding the low-velocity ($\sim$100 km s⁻¹) components, the total mass outflow rates range from $\sim$150 to $\sim$1500 $M_{\odot}$ yr⁻¹ (Table 3), the mass-loading factors are 1.5–8 (Figure 18), the momentum fluxes are (2–20) × $L_{IR}/c$, and the energy fluxes are (0.1–3)% of $L_{IR}$. The corresponding momentum boosts are <10 if the momentum rates due to radiation pressure, winds, and supernovae supplied by the AGN ($2L_{AGN}/c$) and the starburst (3.5$L_{*}/c$) are combined (Section 6.2, Figure 19). These values may be significantly higher in some sources (e.g., IRAS 14348–1447, IRAS 23365+3604, Mrk 231) if the low-velocity components are included.

2. The most powerful outflows, with energy fluxes of $\sim$10¹⁰.⁵–10¹¹ $L_{\odot}$ (Figure 19(f)), are found in IRAS 03158+4227 (belonging to a widely separated pair), IRAS 08572+3915 and IRAS 20100–4156 (mergers with projected separation of $\approx$5.5 kpc), and Mrk 231 (post-merger), apparently uncorrelated with the merging stage. They are best identified by high-velocity absorption in the high-lying OH65 doublet and represent $\sim$20% of all local ULIRGs observed in OH65. This suggests stochastic episodic strong AGN feedback events throughout the ULIRG phase of the merging process.

3. Overall, the energetics inferred from OH are consistent with those derived from CO. The similar outflow masses found for OH and for the more spatially extended CO and Na I suggest that a significant fraction of the total outflow mass is loaded at small radii (a few × 100 pc), and the

outflows are less efficient at evacuating the gas initially located in the extended kiloparsec-scale disks.

4. We use a simple analytical model for momentum-driven flows to check whether radiation pressure due to the AGN and starburst and momentum-conserving winds can account for the observed OH outflows (Section 7.1). Our primary adopted geometry is a combination of a shell for small radii and a collection of clouds at larger radii with decreasing covering factor relative to the released momenta, but we also check the results obtained with no geometrical loss of input momentum. The model also roughly simulates the extended nature of the starbursts and the momentum boost due to radiation trapping. We find that radiation pressure may partially support the gas against gravity (Section 7.1.1 and Appendix C), but is unable to drive the outflows. Inclusion of momentum-conserving winds (Section 7.1.2) could explain some of the outflowing components, in particular if the gravitational potential well is assumed to be balanced by rotation or the input momentum is not geometrically lost. However, partially energy-conserving phases are most likely required for some high-velocity components with high columns (IRAS 03158+4227, IRAS 08572+3915, IRAS 20100–4156, and possibly Mrk 231).

5. Energy-conserving flows are treated separately for the starburst and AGN (Section 7.2). Using limiting values for feedback from a nuclear starburst, some low-velocity components can be accounted for. Winds from the AGN with energy-conserving phases are required to explain some of the highest-velocity OH outflows, with momentum boosts relative to $L_{AGN}/c$ that, for the individual components, do not exceed a factor of $\sim$20 and are mostly $\lesssim$10. These values are mostly lower than the maximum momentum boosts predicted by theoretical models for energy-driven outflows. For the $<1000$ km s⁻¹ outflows probed by OH, the momentum flux carried by the highest-velocity components is usually higher than that carried by lower-velocity components. Our results indicate that energy-driven and momentum-driven outflows may coexist in the same sources.

6. Depletion timescales generally lie between 10⁷ and 10⁸ yr, except in the cases of IRAS 03158+4227, IRAS 08572+3915, and IRAS 20100–4156, where $t_{dep}$ is somewhat shorter (Figure 28). These estimates are comparable to but significantly shorter than the gas consumption timescales. The flow timescales ($\sim R/v$) are much shorter, $\sim$10⁵–10⁶ yr, and also favor recent or ongoing nuclear outflow activity in highly excited and buried sources as probed by the absorption in the high-lying OH65 doublet.

7. We find an anticorrelation between $t_{dep}$ and $L_{AGN}$, indicating that the AGN is primarily responsible for the short $t_{dep} \lesssim 10^7$ yr found in IRAS 08572+3915, IRAS 03158+4227, IRAS 20100–4156, Mrk 231, IRAS 05189–2524, and most probably Mrk 273. For the rest of the sources the relative roles of the AGN and the starburst are more uncertain, but the AGN either dominates or has a significant role in the observed negative feedback.

We thank Avani Gowardhan for providing us with the CO fluxes of IRAS 03158+4227 and IRAS 20100–4156 prior to
the publication. E.G.-A. acknowledges the hospitality of the Harvard-Smithsonian Center for Astrophysics, where part of the present study was carried out. We would also like to acknowledge the contribution made by Will Lewitus who carried out some early data analysis on these outflows while a student intern at NRL. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAFIFSI/OAA/OAP/ OAT, LENS, SISSA (Italy); and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). E.G.-A. is a Research Associate at the Harvard-Smithsonian Center for Astrophysics and thanks the Spanish Ministerio de Economía y Competitividad for support under projects FIS2012-39162-C06-01 and ESP2015-65597-C4-1-R and NASA for support through grant ADAP NNX15AE56G. Basic research in IR astronomy at NRL is funded by the US ONR; J.F. also acknowledges support from NHSC/JPL subcontracts 139807 and 1456609. H.S. acknowledges support from NHSC/JPL subcontract 1483848. S.V. thanks NASA for partial support of this research through NASA grant ADAP NNX16AF24G. This research has made use of NASA's Astrophysics Data System (ADS) and of GILDAS software (http://www.iram.fr/IRAMFR/GILDAS).

Facilities: Herschel Space Observatory (PACS).

## Appendix A
### Individual Sources

Mrk 231. Our three-component model fit is consistent with the more detailed, four-component model reported in GA14: we fit a “quiescent” component (QC) that accounts for the high-lying absorption at systemic velocities, a compact component with gas outflowing up to >1300 km s⁻¹ (HVC), and an extended component (EC) with intermediate velocities. The HVC and EC have associated $\dot{M} \sim 400$ and $\sim 250 M_\odot$ yr⁻¹, respectively. These values are roughly consistent with the inferred values from CO (2–1) by Feruglio et al. (2015) once their values are corrected to match our definition. Tentatively, we find in the current work that the QC is also expanding at $\sim 100$ km s⁻¹ and with $\sim 450 M_\odot$ yr⁻¹, matching the velocity of peak absorption in the OH84 doublet though not that of the OH65 doublet. On the other hand, we find that the EC could be more compact than adopted in GA14 (which was based on interferometric CO (1–0) emission in Cicone et al. 2012); a good match to the line profiles is also obtained with $R \sim 400$ pc.

IRAS 08572+3915. This is an AGN-dominated double-nucleus source with CO detected only in the NW nucleus (Evans et al. 2002), and with no trace of PAH emission at 3.3 µm (Imanishi et al. 2008). It shows strong absorption in the blue wings of OH119, OH79, and OH84 out to −1000 km s⁻¹, and somewhat more tentatively in OH65. At blueshifted velocities, OH84 is almost as strong as OH119 (Figure 7), indicating a compact outflow. There are also hints of the presence of a non-outflowing (i.e., peaking at systemic velocities) component, but it is spectrally diluted by the prominent wings. It is the only source in our sample that shows a P Cygni profile in the highly excited OH65 doublet. Our best-fit model for the blueshifted wing requires two components for the outflow, though both are very compact with similar $r \sim 110$ pc. The lower-velocity (200–800 km s⁻¹) component, with $N_H \sim 1.5 \times 10^{23}$ cm⁻², and the high-velocity

(800–1100 km s⁻¹) component, with $N_H \sim 8 \times 10^{22}$ cm⁻², both have similar $\dot{P}_{\text{out}} \sim (1.3–1.5) \times 10^{36}$ dyn, about 3 times the momentum supplied by the dominant AGN and the starburst ($\sim 5 \times 10^{35}$ dyn). However, Efstathiou et al. (2014) argue that this Seyfert 2 galaxy has a much larger AGN luminosity, but the radiation mainly escapes along the axis of the torus away from our line of sight. If this is the case and the AGN has a luminosity of $\sim 10^{13} L_\odot$ (Efstathiou et al. 2014), it would supply a momentum rate of $2.5 \times 10^{36}$ dyn, which is just the sum of the momentum fluxes estimated for the two OH outflowing components.

With $\sim 100$ pc resolution, Rupke & Veilleux (2013b) have resolved the outflow in H₂ rovibrational lines, showing a blueshifted wing similar to OH119 with extent up to 400 pc and centered at $\sim 200$ pc from the center. This is more extended than the radius inferred for OH, which is perhaps consistent with the high extinction of the nuclear source at far-IR wavelengths. Nevertheless, the OH appears to be missing the more extended component of the outflow. Cicone et al. (2014) indeed measured a radius of $\sim 0.8$ kpc based on CO emission, but, in contrast with other sources, OH does not show evidence for extended emission in the blueshifted wing. This may be due to the compactness of the far-IR continuum. OH84 shows a prominent P Cygni profile (see Figure 1), but the redshifted emission (at velocities between 1000 and 1600 km s⁻¹ relative to the blue component, very likely extended) does not appear to belong to the component responsible for the blue absorption wing. This redshifted feature appears to have its counterpart in CO (1–0) (Cicone et al. 2014) and has not been modeled.

Shirahata et al. (2013) have detected the outflow in the rovibrational lines of CO at $\sim 5$ µm, up to blueshifted velocities of $-400$ km s⁻¹, and derived a column of $N_H \sim 5 \times 10^{22}$ cm⁻². Since the 5 µm continuum emission that CO is absorbing is expected to arise from a very warm, small region close to the AGN (Shirahata et al. 2013), CO and OH probably probe different regions, but the two tracers indicate high columns of outflowing gas close to the central engine.

Mrk 273. This is a double system with the northern nucleus containing most of the molecular mass and probably also an AGN that accounts for roughly 1/3–1/2 of the total luminosity (Rupke & Veilleux 2013a). Mrk 273 is similar in OH79 to IRAS 05189–2524. In addition, both sources (and also IRAS 20551–4250) show an OH119 absorption trough of only $\sim 10\%$ of the continuum, and OH79 correspondingly shows weak absorption and emission in comparison with other ULIRGs with similar luminosity such as IRAS 23365+3604, but extending to higher velocities. This indicates a relatively low covering factor of the continuum by the outflow, perhaps associated with the biconical shape of the outflowing gas (Rupke & Veilleux 2013a). At systemic velocities, the high-lying OH84 and OH65 doublets show strong absorption, most likely associated with the very compact molecular core detected by Downes & Solomon (1998). Mrk 273 shows two well-defined velocity components on the blue side of the OH profiles, both detected in OH84, suggesting a compact molecular outflow. Cicone et al. (2014) indeed found the shortest radius in Mrk 273, $\sim 0.55$ kpc, among the ULIRG outflows observed in CO. The $\sim 300$ km s⁻¹ component, detected in both OH84 and OH65, appears to be extremely compact ($\sim 130$ pc) and has high columns ($N_H \sim 8 \times 10^{22}$ cm⁻²), yielding $\dot{P}_{\text{out}} \sim 2 \times 10^{35}$ dyn. If the northern nucleus harbors an embedded AGN with 1/3–1/2 of the total luminosity, the above momentum flux is 1–2 times that provided by the AGN. While suggesting a moderate momentum boost, overcoming the deep
potential well is the main drawback to explaining this outflow component through momentum-driven approach. For the high-velocity (~700 km s⁻¹) component, which is detected in OH84 with similar strength to OH119, we infer $N_H \sim 1 \times 10^{23}$ cm⁻² and $r \sim 350$ pc, giving $\dot{P}_{\text{out}} \sim 10^{36}$ dyn, i.e., nearly twice the momentum rate supplied by the starburst and the AGN combined. A high momentum boost is required (Figure 27), though we cautiously note the weakness of the OH119, OH79, and OH84 high-velocity wings, which may suggest a jet-driven collimated outflow component.

**IRAS 20551-4250.** This is an advanced merger with a single nucleus that shows evidence for both a deeply embedded AGN and a starburst from hard X-rays to mid-IR emission (Genzel et al. 1998; Franceschini et al. 2003; Haan et al. 2011; Sani & Nardini 2012). The relative contribution from the AGN and the starburst based on mid-IR line diagnostics may be misleading, as [Ne II] and [O IV] are not detected but the equivalent width of the PAH 6.2 µm feature is low, suggesting high extinction (Farrah et al. 2007). The strong OH65 absorption indeed indicates that the source is optically thick even at far-IR wavelengths (GA15). Sani & Nardini (2012) estimated a 26% contribution by the AGN, while Veilleux et al. (2013) estimated 57% from the $f_{15}/f_{30}$ flux density ratio. OH65 shows a blueshifted line wing up to -400 km s⁻¹, while the rest of the doublets show absorption up to -800 km s⁻¹. These are best matched with two outflowing components that provide a total momentum flux of ~4.5 × 10³⁵ dyn, 15% in excess of the combined momentum rate supplied by the AGN and the starburst.

**IRAS 23365+3604.** This LINER-type ULIRG shows no evidence of a double nucleus (Downes & Solomon 1998), and the tidal tails observed in the optical (Murphy et al. 1996) suggest a post-coalescence, late-stage merger (Martin & Soto 2016). About half of its luminosity is associated with an AGN (Veilleux et al. 2013). Besides OH119 and OH79, OH84 also shows a P Cygni profile in this source, and the OH65 peak absorption is blueshifted by ≈100 km s⁻¹ relative to [C II] 157 µm. OH119 shows blue absorption up to -1000 km s⁻¹ and possibly more, but the other doublets only show absorption up to -500 km s⁻¹. Our best-fit model involves three outflowing components, with the CORE component outflowing at ~100 km s⁻¹ ($\dot{M} \sim 340 M_\odot$ yr⁻¹). The other two components (green and magenta curves in Figures 13–16) are probably not independent. They combine to give a total $\dot{P}_{\text{out}} \sim 7.5 \times 10^{35}$ dyn, ~4× the momentum rate supplied by the AGN. If a compact starburst is considered, the outflow in IRAS 23365+3604 requires little momentum boost for both components (≤3).

**IRAS 14348-1447.** No evidence for broad-line emission or very high ionization lines in the optical and near-IR is found in this LINER-type ULIRG (Veilleux et al. 1997), suggesting either a weak AGN contribution to the luminosity or an extremely buried AGN. In the far-IR, the strong absorption found in OH65 (GA15, Figure 16) strengthens the latter possibility. It is a very massive merger, with projected nuclear separation of ~5 kpc, and with the two CO (1–0) components associated with the stellar nuclei (Evans et al. 2000). The southwest (SW) nucleus is the most luminous in both CO (Evans et al. 2000) and the near-IR (Scoville et al. 2000), and the OH profiles suggest that that it is the source of the outflow because its redshift, $z_{\text{SW}} = 0.0827$ (Evans et al. 2000), is in good agreement with the center of the OH119 P Cygni profile (Figure 13). We have used in Figures 13–16 $z_{\text{CII}} = 0.08257$ (see Table 1; but a better fit to OH119 and OH79 is obtained with $z_{\text{SW}}$). A very high velocity wing, in

excess of 1000 km s⁻¹, is observed in OH119. An interesting characteristic of this source is the blueshift of ≈100 km s⁻¹ observed in the OH65 absorption feature (Figures 5 and 16). Although the northeast (NE) nucleus has a redshift that coincides with that of the OH65 feature, the close correspondence between the OH65 and the OH79 absorptions strongly suggests that OH65 is also produced toward the SW nucleus, i.e., the nucleus is expanding. The associated momentum flux, ~7 × 10³⁵ dyn, is of the same order as the momentum rate that can be supplied by the dominant energy source, whether starburst or AGN, and can be explained through a momentum-driven approach if it is rotationally supported. Still, some non-negligible contribution of the NE nucleus may be present, and the above estimates are upper limits. The other, higher-velocity components together account for ~10³⁶ dyn.

**IRAS 05189-2524.** This is an AGN-dominated ULIRG, with $\alpha_{\text{AGN}}$ ranging from 0.6 to 1 using the six diagnostics in Veilleux et al. (2009). It has a single identifiable nucleus and is proposed to be a late-stage merger system, due to the observed tidal tails in the optical (Sanders et al. 1988). Its optical classification is Seyfert 2, but it shows broad (~2000 km s⁻¹) emission near-IR lines ($Pa\alpha$, $Pa\beta$, He) and is thus classified as Seyfert 1 in the near-IR (Severgnini & Risaliti 2000). Spoon & Holt (2009) found from blueshifted [Ne III] and [Ne V] mid-IR emission lines detected in this and other ULIRGs that the gas ionization increases with increasing blueshift. While the ground-state OH119 and OH79 profiles show well-defined P Cygni line shapes, the excited OH84 and OH65 profiles show strong absorption at systemic velocities (GA15). At moderate blueshifted velocities (~-200 km s⁻¹), OH79 has an outflowing counterpart in both OH84 and OH65 (Figure 5), indicating a compact component with $N_H \sim 10^{23}$ cm⁻² and momentum flux ($\dot{P}_{\text{out}} \sim 1.5 \times 10^{35}$ dyn) somewhat lower than that provided by the AGN (~3 × 10³⁵ dyn), suggesting a momentum-driven shell (Figures 25(c) and 27(b)). OH119 shows a flat absorption profile between -200 and -600 km s⁻¹ at ≈10% level, indicating that the high-velocity gas, not detected in OH84 or OH65, is more extended and with a low covering factor. With $r \sim 340$ pc, the high-velocity component can be roughly explained with $N_H \sim 4 \times 10^{22}$ cm⁻², yielding $\dot{P}_{\text{out}} \sim 5 \times 10^{35}$ dyn and thus requiring moderate momentum boost. Although the OH79 absorption at $v < -600$ km s⁻¹ is not quite reproduced, this is the AGN-dominated ULIRG with the lowest momentum boost (Figure 18).

**IRAS 14378-3651.** Duc et al. (1997) found no obvious optical companion to this distorted galaxy, which the authors suggest is a final-stage merger. This is consistent with the OH observations: the prominent outflow probed by the OH119, OH79, and OH84 doublets appears to have expelled a significant fraction of the nuclear molecular gas, leaving a relatively weak OH65 absorption (GA15, Figure 16). A moderate (250 km s⁻¹) velocity component, mostly responsible for the broad and blueshifted OH84 absorption, and a high-velocity component give a total momentum flux of ~5 × 10³⁵ dyn, which is of the same order as the momentum rate supplied by the (dominant) starburst. This suggests at most moderate momentum boost (≤3).

**IRAS 10565+2448.** In the optical (r band) this galaxy appears to be a triple system (Murphy et al. 1996), but CO is only detected in the western nucleus (Downes & Solomon 1998), which dominates the luminosity of the system (Shih & Rupke 2010; Rupke & Veilleux 2013a). The AGN luminosity
may be significantly lower than that of the starburst (Veilleux et al. 2009, 2013; Rupke & Veilleux 2013a) since the value of $\alpha_{AGN} = 0.47$ is higher than that obtained from other diagnostics in Veilleux et al. (2009). The source shows prominent OH119 and OH79 P Cygni profiles with blue absorption trough and red emission strength similar to or stronger than that of AGN-dominated sources, although with lower-velocity extent (Figures 13 and 14). The ground-state OH lines thus indicate outflow parameters for the extended component of the outflow that are similar to AGN-dominated ULIRGs. Our model fits two outflow components with $N_H \sim (1-2) \times 10^{22} \text{ cm}^{-2}$ and $r \sim 400-500 \text{ pc}$ each, yielding $\dot{R}_{\text{tot}} \sim 5.5 \times 10^{35} \text{ dyn}$. This is $\sim 1.2$ times the momentum supplied by the starburst and AGN combined, suggesting a momentum-conserving molecular outflow. The high total outflowing molecular mass of $\sim 3 \times 10^8 M_{\odot}$ still falls short of the outflowing atomic gas of $\sim 10^9 M_{\odot}$ inferred from Na I D by Rupke & Veilleux (2013a), which extends out to at least $\sim 4 \text{ kpc}$ from the western nucleus (Shih & Rupke 2010). In spite of the prominent extended outflow inferred from the ground-state OH119 and OH79 doublets, the source shows weak OH65 and OH84 absorption, indicative of relatively low columns or dust temperatures in the central region (GA15). This may be consistent with the model by Downes & Solomon (1998), who fitted an *unfilled* disk (i.e., with a central gap) to the CO (1–0) interferometric data. In addition, the weak OH84 appears to be globally shifted to the blue (Figure 15), indicating that a significant part of the excited, nuclear region is outflowing, and that the source recently rid itself of most gas in the nuclear region. Together, these clues appear to indicate that, in spite of the high luminosity of the source, the nuclear outflowing activity in IRAS 10565+2448 has subsided in the past $\sim 1$ Myr. This is consistent with the fact that the most obscured regions are found spatially coincident with the highest-velocity outflow, rather than with the western nucleus (Rupke & Veilleux 2013a).

**IRAS 13120–5453.** This ULIRG is classified as a post-merger with a single nucleus (Stierwalt et al. 2013). Its luminosity is dominated by the starburst, though it has a Compton-thick AGN (Iwasawa et al. 2011; Teng et al. 2015). The source is not detected in OH65 (GA15), indicating that the material responsible for the AGN obscuration is too compact to have any significant contribution to the far-IR. Our fit to the OH doublets consists of a central velocity component, best identified in OH84 (Figure 15), and an outflow component with moderate column density ($\sim 1.5 \times 10^{22} \text{ cm}^{-2}$) and probably extended. It gives an estimated momentum flux of $\sim 4 \times 10^{35} \text{ dyn}$, which is significantly lower than the predicted momentum rate supplied by the starburst and the AGN. This suggests a momentum-conserving flow. A weak but apparent blueshifted wing in OH119 up to velocities of $\sim -1500 \text{ km s}^{-1}$ (Figure 13) is not modeled and could reveal the feedback by the buried AGN.

**IRAS 19297–0406.** This is a multiple colliding system dominated by star formation. Though the source is undetected in OH79, it still shows a prominent absorption in OH65 at systemic velocities. The energetics of the outflow are very uncertain, but the model is consistent with an outflowing column of $N_H \sim 10^{22} \text{ cm}^{-2}$, $M_{\text{out}} \sim 4 \times 10^7 M_{\odot}$, $\dot{M} \sim 100 M_{\odot} \text{ yr}^{-1}$, and $\dot{P} \sim 3 \times 10^{35} \text{ dyn}$. This is significantly lower than the momentum rate that a galaxy with $L_{\text{IR}} \sim 2.5 \times 10^{12} L_{\odot}$ can

provide, suggesting that only a single component of the multiple system is involved in the outflow.

**IRAS 09022–3615.** Optically classified as a star-forming galaxy (Lee et al. 2011), the red near-IR continuum in the L-band suggests a buried AGN in this ULIRG (Lee et al. 2012). The OH65 absorption in this source is weak but detected and shows a slight blueshift relative to the systemic velocity. This is the only source in our sample where the emission feature in OH119 is stronger than the absorption feature, thus involving important geometrical effects or significant collisional excitation. Our model is thus relatively uncertain, though it involves $\dot{M} \sim 130 M_{\odot} \text{ yr}^{-1}$ and $\dot{P} \sim 1 \times 10^{35} \text{ dyn}$. As in the case of IRAS 19297–0406, this is significantly lower than the momentum rate that the galaxy could supply based on its luminosity. However, the energetics could be significantly underestimated if the weak absorption is due to low far-IR continuum brightness, rather than to low column densities.

**IRAS 03158+4227.** This galaxy has a companion (projected separation of $\sim 50 \text{ kpc}$) with no apparent sign of interaction in the *r* band (Murphy et al. 1996), but Meusinger et al. (2001) found a curved tail in the companion that matches the expected result of a binary interaction. The optical spectrum shows no indication of an AGN (Meusinger et al. 2001), and Risaliti et al. (2000) find no evidence for an AGN from hard X-ray activity. However, this is (together with IRAS 20100–4156) the most extreme source in our sample. At systemic velocities the OH65 absorption is very deep, with nearly flat absorption between the two doublet components (Figure 1). This indicates very high extinction in the far-IR, capable of obscuring any AGN signature at shorter wavelengths. On the other hand, the outflow is observed out to $\gtrsim 1500 \text{ km s}^{-1}$ in OH119, as emphasized by Spoon et al. (2013). In OH79 the observed P Cygni profile has a blueshifted wing apparently extending up to $-1700 \text{ km s}^{-1}$, though it may be contaminated by H$_2$O primarily at velocities more blueshifted than $-1300 \text{ km s}^{-1}$. The OH84 absorption trough is deep, 10% at $-800 \text{ km s}^{-1}$, and a very prominent OH65 absorption is seen up to $-1000 \text{ km s}^{-1}$. Our simple three-component model approach cannot account for all spectroscopic details in the wide velocity range depicted by the OH doublets. Even neglecting the high-velocity wing in OH79, the OH65 absorption at central velocities is underpredicted and OH84 is slightly overpredicted over the full velocity range. Nevertheless, the model captures the presence of a compact, massive ($4 \times 10^8 M_{\odot}$), and excited component flowing at high velocities (300–1300 km s$^{-1}$), with an estimated mass outflow rate of $1300 M_{\odot} \text{ yr}^{-1}$ and a momentum flux of $\sim 7 \times 10^{36} \text{ dyn}$. These astonishing values cannot be understood in terms of the momentum rate supplied by a galaxy with even such high luminosity as IRAS 03158+4227, and a very high momentum boost is required. This is probably an “exploding” quasar, with enormous masses of gas suddenly accelerated to high velocities, and a nuclear-energy-conserving phase within this extremely buried source is required to understand the involved energetics. In addition, another outflowing component is required to account for the extremely high velocity OH wings. Our model for the latter component is more uncertain but probably does not dominate the energetics of the source and can be explained through a momentum-conserving phase (Figure 25(c)).

**IRAS 20100–4156.** From optical and near-IR images, Duc et al. (1997) classify this galaxy as a close interacting pair of disk galaxies (projected separation of $\approx 5.4 \text{ kpc}$), with a starburst or
<figure><img src="image_34.png" /><figcaption>**Figure 31.** Comparison of the outflow properties of local starburst galaxies inferred by Heckman et al. (2015) (black symbols) with those of ULIRGs as inferred from OH (individual components, in red) in this work. The column density ($N_H$), outflowing mass ($M_{\text{out}}$), velocity ($v$), mass outflow rate ($\dot{M}_{\text{out}}$), and momentum flux ($\dot{P}_{\text{out}}$) are shown as a function of the estimated radius ($r$).</figcaption></figure>

LINER-type spectrum. A possible third component of the merger was identified by Bushouse et al. (2002) with NICMOS. The galaxy shows strong megamaser emission in OH at centimeter wavelengths (Harvey-Smith et al. 2016), with a plateau of OH emission whose blue wing traces the velocities of the absorption found in the far-IR OH doublets. Searches in X-rays have identified the possible presence of a buried AGN (Franceschini et al. 2003). The OH doublets indicate a galaxy with extreme properties. The OH65 absorption is very deep and broad, with little decrease of the absorption at wavelengths between the two doublet components. Absorption in OH119 is observed out to $-1300 \text{ km s}^{-1}$, and to $-1000 \text{ km s}^{-1}$ in the other three doublets. Similar to IRAS 03158+4227, our model involves a compact, massive ($5 \times 10^8 M_{\odot}$), and excited component flowing at velocities up to $800 \text{ km s}^{-1}$ with a momentum flux of $\sim 3 \times 10^{36} \text{ dyn}$. It most likely requires a high momentum boost. A second, even higher-speed component with one order of magnitude lower mass accounts for $\sim 1.5 \times 10^{36} \text{ dyn}$, and could be explained via a momentum-driven outflow. The OH119 absorption is underestimated by our model on the red side of the profile, suggesting additional, foreground absorption. The very high velocities in this source suggest that the AGN contribution to the luminosity is underestimated.

## Appendix B

### Comparison with Ionized Lines in Starburst Galaxies

We compare here the outflow properties of local starburst galaxies inferred by Heckman et al. (2015) and of ULIRGs as

derived in this work. Since these properties are derived in basically the same way (i.e., Equations (6), (11), and (12), while the starburst galaxies are analyzed in lines of ions, the comparison is relevant to contrast different phases of the outflow and to infer the main differences between pure starburst galaxies and ULIRGs where the AGN plays an important role. Figure 31 shows the outflow column density, mass, velocity, mass outflow rate, and momentum flux, as a function of the outflow radius, for the starburst galaxies (black symbols) and the ULIRGs (individual components, in red).

A quantitative difference between the two phases is found in their column densities, which are moderate in the ionized outflow phase of starburst galaxies, $\sim 10^{21} \text{ cm}^{-2}$, but are much larger, $\gtrsim 10^{22} \text{ cm}^{-2}$, for the molecular outflows in ULIRGs. The outflow size, which is adopted as twice the size of the starburst in Heckman et al. (2015), is also larger than the extents we have derived for the OH outflows. The outflowing gas mass (Equation (6)) is then found roughly similar in both types of sources and phases. The derived mass outflow rates and momentum fluxes are significantly higher for the OH outflows.

## Appendix C

### Radiation Pressure Support of the Highly Excited Nuclear Structures

Thompson et al. (2015) consider two types of structures that are subject to radiation pressure: a spherical shell and a cloud ensemble. We show here that the pressure exerted by the
<figure><img src="image_35.png" /><figcaption>**Figure 32.** Radiation pressure support in (U)LIRGs with measurable absorption in the high-lying OH65 doublet. (a) Planck-averaged mass-absorption coefficient ($k_{IR}$) as a function of the dust temperature ($T_{dust}$). (b) The gray curve shows $L_{EDD}/M_{gas}$ as predicted by Equation (26) with a gas fraction of $f_g = 0.15$ and the $k_{IR}$ profile of panel (a). The solid line (squares) shows the calculated $L_{IR}/M_{gas}$ for a continuum optical depth at 100 µm of unity (GA15). For these optically thick sources we expect radiation pressure support for $T_{dust} \sim 140$ K. (c) Infrared luminosities ($L_{IR}$) as a function of the velocity dispersion ($\sigma$) for all (U)LIRGs where OH65 is detected (GA15). $\sigma$ is calculated from Gaussian fits to the OH65 doublet. The dashed and dotted lines indicate the Eddington luminosities $L_{EDD}$ predicted by Equation (27) for $T_{dust} = 100$ and 150 K, respectively.</figcaption></figure>

radiation responsible for the observed molecular excitation in ULIRGs is close to that required for support against gravity, but probably not enough to overcome gravity and drive an outflow under typical physical conditions in the nuclear regions of ULIRGs. Starting with a shell at initial position $r_0$ and IR optical depth $\tau_{IR} \gg 1$, the ratio of Eddington luminosity to gas mass is (from Equation (8) in Thompson et al. 2015)

$$ \frac{L_{\text{EDD}}}{M_{\text{gas}}} = \frac{4\pi c G}{f_g k_{\text{IR}}} = 1.76 \times 10^4 \times \left(\frac{0.15}{f_g}\right) \times \left(\frac{5 \text{ cm}^2/\text{g}}{k_{\text{IR}}}\right) L_{\odot}/M_{\odot}, \quad (26) $$

where $M_{gas}$ denotes the gas mass within $r_0$, $f_g$ is the gas mass fraction in that region, and $k_{IR}$ (dependent on $T_{dust}$) is the IR opacity. Combes et al. (2013) estimate that the average $f_g$ in local galaxies is $\lesssim 0.10$, though with high dispersion, and Downes & Solomon (1998) derived $f_g \sim 0.17$ in the nuclear regions of ULIRGs. We adopt a reference value of $f_g = 0.15$.

Figure 32(a) shows the Planck-averaged opacity $k_{IR}$ as a function of $T_{dust}$, based on the mass-absorption coefficient of dust given in González-Alfonso et al. (2014a) and a nominal gas-to-dust ratio by mass of $g/d = 100$. The reference value given by Thompson et al. (2015), $k_{IR} = 5 \text{ cm}^{-2} \text{ g}^{-1}$ of gas, is obtained for $T_{dust} = 140$ K. Based on this $k_{IR}$ profile and on Equation (26), we plot in Figure 32(b) $L_{EDD}/M_{gas}$ (gray curve)

as a function of $T_{dust}$. We also show in Figure 32(b) the calculated values of $L_{IR}/M_{gas}$ (squares, solid line) for a spherical source with $\tau_{100} = 1$, the minimum value we inferred from radiative transfer models to the OH65 doublet (GA15). To calculate $M_{gas}$ from $\tau_{100}$, we adopt an absorption coefficient of $0.45 \text{ cm}^{-2} \text{ g}^{-1}$ of gas at 100 µm (based again on González-Alfonso et al. 2014a) and $g/d = 100$.

The figure indicates that $T_{dust} \approx 140$ K is required for radiation pressure support against gravity for $\tau_{100} = 1$. For $\tau_{100} > 1$, the $L_{IR}/M_{gas}$ values would be lower because $L_{IR}$ would remain nearly the same owing to the fact that the far-IR emission is optically thick; hence, higher $T_{dust}$ would be required.

In the far-IR, the values of $T_{dust}$ as directly inferred from the OH and H$_2$O excitation in buried nuclei (GA12; GA14; GA15; Falstad et al. 2015) are 100–150 K, suggesting that these nuclear regions are close to the regime of radiation pressure support against gravity in the vertical direction (Thompson et al. 2005), but may be insufficient to drive an outflow. The most excited source in our original sample is the LIRG NGC 4418, where $T_{dust} \sim 150$ K is inferred (GA12), but no molecular outflow is detected in this source. We emphasize that in these sources that are optically thick in the far-IR, higher $T_{dust}$ are expected to be prevalent within the cocoon of dust, due to radiation trapping, but the relevant values for radiation pressure support at the scales probed by the far-IR are the