FAST MOLECULAR OUTFLOWS IN LUMINOUS GALAXY MERGERS: EVIDENCE FOR QUASAR FEEDBACK FROM HERSCHEL

S. Veilleux; M. Meléndez; E. Sturm; J. Gracia-Carpio; J. Fischer; E. González-Alfonso; A. Contursi; D. Lutz; A. Poglitsch; R. Davies; R. Genzel; L. Tacconi; J. A. de Jong; A. Sternberg; H. Netzer; S. Hailey-Dunsheath; A. Verma; D. S. N. Rupke; R. Maiolino; S. H. Teng; E. Polisensky

doi:10.1088/0004-637X/776/1/27

1. INTRODUCTION

Gas-rich galaxy merging may trigger major starbursts, lead to the formation of elliptical galaxies, and account for the growth of supermassive black holes (BH; e.g., Sanders et al. 1988; Hopkins et al. 2009). This merger-driven evolutionary scenario starts with a completely obscured ultraluminous infrared galaxy (ULIRG). As the system evolves, the obscuring gas and dust is gradually dispersed, giving rise to dusty QSOs and finally to completely exposed QSOs. Powerful winds, driven by the central quasar or the surrounding starburst, have been invoked to stop the growth of both the BH and spheroidal component and explain the tight BH–spheroid mass relation (e.g., Fabian 1999; King 2003; Murray et al. 2005). These winds are purported to quench star formation in the merger remnants ("negative mechanical feedback"), creating a population of red gas-poor ellipticals and explaining the bimodal color distribution observed in galaxy surveys (e.g., Kauffmann et al. 2003). There is growing observational support for these winds: e.g., most galaxies with high star formation rate (SFR) densities show signatures of outflows, both locally and at high redshifts (e.g., Heckman 2002; Veilleux et al. 2005; Chen et al. 2010; Weiner et al. 2009; Rubin et al. 2010; Steidel et al. 2010; Kornei et al. 2012; Martin et al. 2012; Newman et al. 2012, and references therein). The local outflows are often spatially resolved, allowing determination of the mass outflow rates (∼0.1–5 × SFR) and kinetic energies (∼10^56–58 erg in ULIRGs; Rupke et al. 2002, 2005a, 2005b, 2005c; Martin 2005, 2006). These winds are primarily driven by the starburst rather than the active galactic nucleus (AGN), except in late-stage mergers with quasar signatures where velocities ≳1000 km s⁻¹ (cf. 100–400 km s⁻¹ in other systems) are sometimes observed (e.g., Rupke et al. 2005c; Krug et al. 2010; H. B. Krug et al. 2013, in preparation).

Studies conducted in the past three years are shedding new light on the wind phenomenon in the local universe: (1) Early results from our Herschel¹³ guaranteed time key program SHINING (PI: E. Sturm) have revealed far-infrared (FIR) OH features with P Cygni profiles indicative of massive molecular outflows in a number of ULIRGs (Fischer et al. 2010; Sturm et al. 2011, hereafter F10 and S11, respectively; also Gonzalez-Alfonso et al. 2012, 2013). In a few of these objects, particularly those with dominant AGNs such as Mrk 231, the nearest quasar known, the terminal outflow velocities exceed 1000 km s⁻¹ and the outflow rates (up to ∼1000 M_☉ yr⁻¹) are several times larger than the infrared (IR) based SFRs. (2) Independent, spatially resolved CO-emission observations of Mrk 231 with the IRAM/PdB millimeter-wave interferometer (Feruglio et al. 2010; Cicone et al. 2012; it has also been mapped in HCN, HCO⁺, and HNC by Aalto et al. 2012) have confirmed this outflow and deduced mass outflow rates of ∼700 M_☉ yr⁻¹, far larger than the ongoing IR-based SFR (∼160 M_☉ yr⁻¹) in the host galaxy. Further millimeter CO-observations have now been done by our group in a number of the objects discussed here, finding similar CO outflows (Cicone et al. 2013). Remarkably, the CO outflow in Mrk 231 partially overlaps spatially and kinematically with blueshifted optical Na i D λλ5890, 5896 absorption features detected out to ∼2–3 kpc from the nucleus (Rupke et al. 2005c). (3) Our high-resolution Gemini/Integral Field Unit (IFU) observations have revealed that the neutral-atom Na i D outflow in Mrk 231 is wide-angle, thus driven by a QSO wind rather than a jet (Rupke & Veilleux 2011, hereafter RV11; for short, neutral-atomic Na i D outflows are simply called neutral outflows in the rest of this paper, unless noted otherwise). Similar wide-angle high-velocity outflows have now been mapped at high resolution in a number of other local ULIRGs (Rupke & Veilleux 2013a, hereafter RV13, 2013b).

These powerful outflows may be the long-sought "smoking gun" of quasar mechanical feedback that clears out the molecular disk formed from dissipative collapse during the merger. The tentative trend of increasing OH terminal outflow velocity with increasing AGN luminosity identified by S11 adds support to this idea, but it is based on only six ULIRGs (and NGC 253, a regular starburst galaxy). To properly test this idea, one needs to examine a statistically representative sample of ULIRGs and study the properties of their molecular outflows uniformly and systematically. In this paper, we report the first results from such a study, focusing on the directly measured outflow velocities (in a later paper, we will report on our analysis of the model-dependent dynamical quantities (masses, momenta, and energies) in the subset of objects for which we have multiple OH transitions). The sample used for this analysis is described in Section 2. The observations and methods used to reduce and analyze the data are discussed in Section 3. The results from this analysis are presented in Section 4 and their implications are discussed in Section 5. The main conclusions are summarized in Section 6. Throughout this paper, we adopt H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and $\Omega _\Lambda$ = 0.7. We also adopt the standard convention that approaching material has a negative velocity with respect to the systemic velocity of the host galaxy due to Doppler shift.

2. SAMPLE

2.1. Sample Selection

As mentioned in Section 1, the first results from the SHINING GTO survey of ULIRGs (S11) revealed a tentative trend between OH outflow velocities and AGN luminosities, but this was largely driven by the velocities measured in Mrk 231 and F08572+3915, the only two warm (large f₂₅/f₆₀ flux ratio), quasar-dominated, late-stage mergers in the S11 sample. It quickly became clear that the rest of the SHINING GTO data would not be able to address this question adequately because the GTO sample largely probes the early stages of the merger after the first peri-passage: 20 of the 23 (87%) GTO ULIRGs are cold (small f₂₅/f₆₀ ratio) pre-mergers or recent mergers, where the bolometric and mechanical luminosities are still dominated by the nuclear starburst (Veilleux et al. 2009b, hereafter V09). A cycle 1 open-time program (OT1 $\_$ sveilleu $\_$ 1, PI: S. Veilleux) has allowed us to expand the original GTO sample to include a set of warm quasar-dominated ULIRGs. These objects were selected using four criteria. (1) The targets had to be part of the QUEST sample of local (z < 0.3) luminous mergers: either ULIRGs or QSOs (Veilleux 2012, and references therein). (2) Their bolometric luminosity had to be dominated by the quasar based on the Spitzer data or, equivalently, IRAS-band 25 to 60 μm flux ratio f₂₅/f₆₀ ≳ 0.15 (this criterion also automatically selects late-stage, fully coalesced mergers; see V09). (3) High signal-to-noise ratio (S/N ∼ 30–40) Na i D absorption spectra had to be planned or available for each object. This criterion did not bias the sample in any way since these objects were selected randomly based on visibility at the time of the ground-based observations. (4) They had to have f₁₀₀ > 1 Jy so that high S/N in the continuum could be reached in a reasonable amount of time with Herschel-PACS. These criteria resulted in a sample of 18 objects, 3 of which were already part of the SHINING GTO survey of ULIRGs. This sample size was a good match to the existing 20 cold pre-merger ULIRGs in SHINING.

Since all ULIRGs in the SHINING survey and its OT1 extension were selected to have IRAS-band 100 μm fluxes f₁₀₀ > 1 Jy to more easily probe OH 119 μm, this sample did not contain any "classic" IR-faint QSOs in the critical late merger phases when the quasar has finally gotten rid of its natal "cocoon" and the effects of feedback are predicted to subside (Narayanan et al. 2008; Hopkins et al. 2009). Therefore we put in a successful open-time request in Cycle 2 (OT2 $\_$ sveilleu $\_$ 4, PI: S. Veilleux) to further expand the sample by including objects meeting the following three additional criteria. (1) IR-to-bolometric luminosity ratios smaller than 0.6 i.e., IR-fainter than all ULIRGs in the SHINING and OT1 programs (IR here and throughout this paper refers to 8–1000 μm). (2) High-S/N (∼30–40) Na i D and FUV absorption spectra had to be planned or available for each object. The new FUV criterion favored the "classic" UV-bright QSOs. (3) The f₁₀₀ threshold was lowered by a factor of two down to >0.5 Jy to accommodate these IR-fainter systems. This factor of two was chosen to provide enough objects in the OT2 sample without requiring unrealistically long exposure times per target. These criteria resulted in a sample of five objects. While this sample size is small, recall that the main objective of this OT2 extension is to anchor the merger sequence at the latest (IR-faintest) stages that could be probed with Herschel. This is a merger phase that has not been explored with Herschel until now so each of these objects provides important new information.

In summary, the sample consists of 43 objects, 23 from the GTO sample, 15 from the OT1 program, and 5 from the OT2 program. All of these objects show signs of ongoing or recent interactions (NGC 4418 = F12243−0036, a GTO target, is not a major merger but it is IR luminous so we include it in the present sample). We include here the ULIRGs from S11 again, applying the same fitting method as for the rest of the sample, in order to be complete and consistent. The main objective of the present study is to search for trends with AGN and host properties as well as age across the merger sequence (after the first encounter, post-LIRG phase). In Section 3, we will subdivide these 43 mergers into various subcategories (e.g., morphology, starburst power, AGN luminosity) to try to find the driving parameter of the wind properties. Our Spitzer results (V09) have shown that BH accretion, in addition to depending on the merger phase, also has a strong chaotic/random component, so a sample size substantially smaller than 43 objects would make it hard to draw statistically significant conclusions.

2.2. Properties of the Sample Galaxies

The properties of our sample galaxies are listed in Table 1. The notes to Table 1 briefly explain the meaning of each of these quantities. The luminosities listed in Table 1 rely on a number of assumptions that deserve further explanation. Following V09, the bolometric luminosities for ULIRGs were estimated to be L_BOL = 1.15 L_IR, where L_IR is the IR luminosity over 8–1000 μm (Sanders & Mirabel 1996), and L_BOL = 7 L(5100 Å) + L_IR for the Palomar–Green (PG) QSOs (Netzer et al. 2007). The starburst and AGN luminosities were next calculated from

$\begin{eqnarray} L_{\rm BOL} &=& L_{\rm AGN} + L_{\rm SB}\\ \ms &=& \alpha _{\rm AGN} L_{\rm BOL} + L_{\rm SB}, \end{eqnarray} \tag{ 1 }$

where α_AGN is the fractional contribution of the AGN to the bolometric luminosity, hereafter called the "AGN fraction" for short. In V09, we measured α_AGN in the QUEST ULIRGs and QSOs using six independent methods that span a range in wavelength and give consistent results within ∼±10%–15% on average. These six methods relied on (1) the [O iv]/[Ne ii] ratio, (2) the [Ne v]/[Ne ii] ratio, (3) the polycyclic aromatic hydrocarbon (PAH) 7.7 μm equivalent width, (4) the modified version of the Laurent et al. (2000) diagram which plots the 14–16 μm to 5.3–5.8 μm continuum ratio versus the PAH 6.2 μm equivalent width, (5) the PAH-free, silicate-free 5–25 μm to FIR continuum ratio, and (6) the 15 to 30 μm continuum ratio, f₁₅/f₃₀.

Table 1. Galaxy Properties

Name	Other Name	z	$\frac{f_{15}}{f_{30}}$	α_AGN	log L_bol	log L_SB	log L_AGN	σ_*	M_H	IC	Type	Program (t_exp)
Name	Other Name	z	$\frac{f_{15}}{f_{30}}$	(%)	(L_☉)	(L_☉)	(L_☉)	(km s⁻¹)	(mag)	IC	Type	(hr)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
F00509+1225	I Zw 1	0.0610	0.434	90.1	12.01	11.00	11.96	188	⋅⋅⋅	S	AGN 1	OT1 (1.3)
F01572+0009	Mrk 1014	0.1631	0.161	64.6	12.68	12.23	12.49	200	−25.39	S	AGN 1	OT1 (1.0)
F05024−1941		0.1920	0.050	7.3	12.43	12.40	11.30	⋅⋅⋅	−24.49	B	AGN 2	OT1 (4.6)
F05189−2524		0.0426	0.198	71.7	12.22	11.67	12.07	137	−23.96	S	AGN 2	GTO (0.5)
07251−0248		0.0876	0.075	30.0	12.45	12.29	11.92	⋅⋅⋅	⋅⋅⋅	S	⋅⋅⋅	GTO (1.1)
F07599+6508		0.1483	0.377	87.6	12.59	11.68	12.53	⋅⋅⋅	−24.55	S	AGN 1	OT1 (3.7)
F08572+3915		0.0584	0.191	70.4	12.20	11.67	12.05	⋅⋅⋅	−23.58	B	L	GTO (1.0)
09022−3615		0.0596	0.125	54.9	12.35	12.00	12.09	⋅⋅⋅	⋅⋅⋅	S	⋅⋅⋅	GTO (0.5)
F09320+6134	UGC 5101	0.0394	0.130	56.4	12.05	11.69	11.80	⋅⋅⋅	⋅⋅⋅	S	L	GTO (0.3)
F10565+2448		0.0431	0.105	47.1	12.10	11.82	11.77	⋅⋅⋅	⋅⋅⋅	B	H	GTO (0.5)
F11119+3257		0.1890	0.266	79.9	12.71	12.02	12.61	⋅⋅⋅	−23.44	S	AGN 1	OT1 (4.6)
Z11598−0114		0.1507	0.210	73.5	12.60	12.03	12.47	⋅⋅⋅	−25.06	S	AGN 1	OT1 (0.8)
F12072−0444		0.1286^a	0.220	74.8	12.46	11.86	12.33	⋅⋅⋅	−24.22	B	AGN 2	OT1 (0.8)
F12112+0305		0.0733	0.060	17.8	12.38	12.30	11.63	⋅⋅⋅	−24.36	B	L	GTO (0.5)
F12243−0036	NGC 4418	0.0071	0.128	55.7	11.06	10.71	10.81	⋅⋅⋅	⋅⋅⋅	S	AGN 2	GTO (0.2)
F12265+0219	3C 273	0.1583	0.840	98.5	12.86	11.04	12.86	⋅⋅⋅	−28.06	S	AGN 1	OT1 (1.0)
F12540+5708	Mrk 231	0.0422	0.272	80.5	12.60	11.89	12.51	120	−24.52	S	AGN 1	GTO (0.3)
13120−5453		0.0308	0.080	33.4	12.30	12.12	11.83	⋅⋅⋅	⋅⋅⋅	S	AGN 2	GTO (0.3)
F13305−1739		0.1484	0.377	87.5	12.32	11.42	12.26	⋅⋅⋅	−25.47	S	AGN 2	OT1 (4.9)
F13428+5608	Mrk 273	0.0378	0.081	34.2	12.21	12.03	11.74	285	−24.32	S	AGN 2	GTO (0.3)
F13451+1232		0.1217	0.273	80.6	12.38	11.67	12.29	146,167^b	−25.52	B	AGN 2	OT1 (1.0)
F14348−1447		0.0830	0.060	17.4	12.40	12.31	11.64	160	−24.99	B	L	GTO (1.1)
F14378−3651		0.0676	0.064	21.1	12.17	12.07	11.50	153	⋅⋅⋅	S	AGN 2	GTO (1.1)
F14394+5332		0.1045	0.152	62.5	12.17	11.74	11.96	⋅⋅⋅	−24.88, −23.56^c	B^b	AGN 2	OT1 (1.0)
F15206+3342		0.1244	0.158	63.9	12.30	11.86	12.10	⋅⋅⋅	−24.47	S	H	OT1 (1.5)
F15250+3608		0.0552	0.095	42.3	12.10	11.86	11.73	150	⋅⋅⋅	S	L	GTO (0.5)
F15327+2340	Arp 220	0.0181	0.049	5.8	12.22	12.19	10.98	164	⋅⋅⋅	B	L	GTO (0.9)
F15462−0450		0.0998	0.145	60.6	12.27	11.86	12.05	169	−23.4	S	AGN 1	OT1 (1.0)
F16504+0228	NGC 6240	0.0245	0.130	56.2	11.90	11.54	11.65	229	⋅⋅⋅	B	L	GTO (0.3)
F17207−0014		0.0428	0.038	⩽5%	12.45	12.43	⩽11.15	229	⋅⋅⋅	S	L	GTO (0.3)
F19297−0406		0.0857	0.066	23.4	12.44	12.33	11.81	⋅⋅⋅	⋅⋅⋅	S	H	GTO (0.3)
19542+1110		0.0624^d	0.069	25.5	12.12	11.99	11.52	⋅⋅⋅	⋅⋅⋅	S	L	GTO (1.1)
F20551−4250		0.0430	0.132	56.9	12.11	11.74	11.87	140	⋅⋅⋅	S	L	GTO (0.5)
F22491−1808		0.0778	0.057	14.4	12.90	12.83	12.05	⋅⋅⋅	−24.21	B	H	GTO (0.8)
F23128−5919	ESO 148-IG 002	0.0446	0.154	63.0	12.09	11.66	11.89	⋅⋅⋅	⋅⋅⋅	B	H	GTO (0.3)
F23233+2817		0.1140	0.204	72.5	12.11	11.55	11.97	⋅⋅⋅	−24.54	S	AGN 2	OT1 (1.0)
F23365+3604		0.0645	0.100	44.6	12.22	11.97	11.87	145	⋅⋅⋅	S	L	GTO (0.5)
F23389+0300		0.1450	0.066	22.7	12.19	12.07	11.54	⋅⋅⋅	−24.15	B	AGN 2	OT1 (4.4)
PG 1126−041	Mrk 1298	0.0620	0.329	88.7	11.52	10.57	11.47	194	−24.15	S	AGN 1	OT2 (3.1)
PG 1351+640		0.0882	0.389	92.2	12.04	10.93	12.00	⋅⋅⋅	⋅⋅⋅	S	AGN 1	OT2 (10.6)
PG 1440+356	Mrk 0478	0.0780	0.361	90.7	11.80	10.77	11.76	⋅⋅⋅	−24.65	S	AGN 1	OT2 (5.4)
PG 1613+658	Mrk 0876	0.1290	0.398	92.6	12.29	11.16	12.26	⋅⋅⋅	−25.68	S	AGN 1	OT2 (4.1)
PG 2130+099	UGC 11763	0.0630	0.563	98.5	11.77	9.95	11.76	156	−24.14	S	AGN 1	OT2 (14.0)

Notes. Column 1: galaxy name. Coordinate-based names beginning with F are sources in the IRAS Faint Source Catalog. Column 2: another name. Column 3: redshift. Column 4: f₁₅/f₃₀ values from V09 and Stierwalt et al. (2013). When no 15 to 30 μm data is available the value is estimated by using the f₂₅/f₆₀ ratio and the good correlation between the 25 to 60 versus 30 to 15 μm flux ratios (see V09). Column 5: α_AGN, fractional contribution of the AGN to the bolometric luminosity based on the f₁₅/f₃₀. See Section 2.2. Column 6: bolometric luminosity. See Section 2.2 for more detail. Column 7: starburst bolometric luminosity. See Section 2.2 for more detail. Column 7: AGN luminosity. See Section 2.2 for more detail. Column 8: σ_*, stellar velocity dispersion from Dasyra et al. (2006a, 2006b, 2007). Column 10: M_H, absolute H-band magnitudes from Veilleux et al. (2006, 2009b). For the sake of completeness we also included M_H absolute magnitudes derived from $M_{k^{\prime }}$ values assuming H–K' = 0.35 mag (Veilleux et al. 2002). Column 11: interaction class, "S" for singles and "B" for binaries, based on classifications of Veilleux et al. (2002, 2006, 2009b) and, when not available, our own visual classification from HST/ACS F814W images (Iwasawa et al. 2011). Column 12: optical spectral types of our sample, H ii (for H ii galaxies), L (for LINER-like), AGN 2 and AGN 1 (for type 2 AGN and type 1 AGN, respectively). For the spectral type we adopted the classification from Veilleux et al. (1999) and V09, when not available, we used the values from NED/SIMBAD. Column 13: program, GTO SHINING program (PI: E. Sturm), QUEST OT1 and OT2 programs (PI: S. Veilleux). In parenthesis is the total PACS exposure time. ^aRedshift from Spoon & Holt (2009). This value is based on [Ne ii] 12.8 μm and is slightly larger than that derived from optical data (e.g., Strauss et al. 1992). ^bThe values presented correspond to the eastern and western nucleus of F13451+1232, respectively. ^cF14394+5332 is a triple system with an eastern pair and a western component. ^dAverage redshift derived from our fit to the [Ne ii], [Ne iii] and [O iv] emission lines in the Spitzer/IRS spectra.

Download table as: ASCII Typeset image

Unfortunately, several objects in the SHINING GTO sample are not part of the QUEST sample and therefore do not have such accurately measured AGN contribution to the bolometric luminosity. In the following discussion, we rely solely on the restframe f₁₅/f₃₀ continuum ratio, which is available for all the objects in the Herschel sample, to quantify the AGN contributions to the bolometric luminosities of our systems. This ratio was found by V09 to be more tightly correlated with the PAH-free, silicate-free mid-IR (MIR)/FIR ratio and the AGN contribution to the bolometric luminosity than any other Spitzer-derived continuum ratio.

Following V09, we use f₃₀/f₁₅ as a surrogate of the PAH-free, silicate-free MIR/FIR ratio and adopt log(f₃₀/f₁₅) = 0.20 and 1.35 as the zero points for the "pure" AGN and the "pure" starburst ULIRG, respectively (note that V09 make a distinction between "pure" starburst ULIRGs and "normal" star-forming galaxies). The AGN contribution is calculated from a linear interpolation between these two extremes. The zero point for the pure AGN corresponds to the average f₃₀/f₁₅ ratio of FIR-undetected PG QSOs, while the zero point for the pure starbursts is calculated from the 10 ULIRGs with the largest f₃₀/f₁₅ ratios. We used the correction factors listed in Table 10 of V09 to transform 15 μm luminosities into bolometric luminosities: log[L_ν(15 μm)/L_BOL] = −14.33 for a pure AGN and −14.56 for a pure starburst ULIRG. The uncertainty on α_AGN is estimated to be ±20% on average, but is likely higher for strongly buried sources with significant FIR extinctions (e.g., Arp 220 = F15327+2340). We return to the issue of FIR extinction in Sections 4.6 and 5.

Figure 1 shows the distributions of redshifts, starburst luminosities, AGN fractions, and AGN luminosities for all 43 objects in our sample. This figure emphasizes the broad range of properties of our sample galaxies. The starburst and AGN luminosities are independent of each other i.e., we find no correlation or anti-correlation between the starburst and AGN luminosities in our sample of galaxies, but perhaps not surprisingly objects with the largest starburst (AGN) luminosities (≳10¹² L_☉) tend to have lower (higher) AGN fractions.

3. OBSERVATIONS, DATA REDUCTION, AND SPECTRAL ANALYSIS

3.1. Observations

All of the data in this paper were obtained with the PACS FIR spectrometer (Poglitsch et al. 2010) on board Herschel (Pilbratt et al. 2010). Several systems in the OT1 and OT2 sub-samples are relatively faint in the FIR so an exhaustive multi-line PACS survey was not possible. Contrary to the GTO data, which generally cover the ground-state OH 119 μm ²Π_3/2 J = 5/2–3/2 rotational Λ-doublet transitions, the high-lying (E_lower = 290 K) 65 μm ²Π_3/2 J = 7/2–9/2 and (E_lower = 120 K) 84 μm ²Π_3/2 J = 7/2–5/2 rotational Λ-doublet transitions, and the ground-state cross-ladder 79 μm ²Π_1/2–²Π_3/2 J = 1/2–3/2 rotational Λ-doublet transitions, the OT1 and OT2 data focus on a single OH feature, the OH 119 μm doublet (and ¹⁸OH 120 μm counterparts). The value of this feature as a wind diagnostic was beautifully demonstrated in F10 and S11 and our early GTO data. It is the strongest OH transition in Mrk 231 and most of the GTO targets. It is positioned in wavelength near the peak spectroscopic sensitivity of PACS. Finally, it is seen in absorption more often than the 79 μm feature and thus is more likely to provide an unambiguous signature of gas motion (rotation, inflow, or outflow).

PACS was used in range scan spectroscopy mode in high sampling centered on the redshifted OH 119 μm + ¹⁸OH 120 μm complex with a velocity range of ∼8000 km s⁻¹ (restframe 118–121 μm) to provide enough coverage on both sides of the OH complex for reliable continuum placement. The resulting PACS spectral resolution is ∼270 km s⁻¹. The total amount of time spent to carry out these observations (including overheads) was 13.0, 32.5, and 37.2 hr during GT, OT1, and OT2, respectively. The on-target exposure times were adjusted according to the continuum level under OH, estimated from an extrapolation of the Spitzer MIR + IRAS FIR energy distribution. A small chopper throw of 1 farcm 5 was used in most cases. The observing time for each target, including all overheads, is listed in the last column of Table 1.

3.2. Data Reduction

All of the PACS data were reduced in the same way as in S11. As described there, the data reduction was done using the standard PACS reduction and calibration pipeline (ipipe) included in HIPE 6.0. The spectra were normalized to the telescope flux (which dominates the total signal) and recalibrated it with a reference telescope spectrum obtained from dedicated Neptune observations during the Herschel performance verification phase. All of our objects are point sources for PACS. In the following, we use the spectrum of the central 9'' × 9'' spatial pixel (spaxel) only, applying the point-source correction factors (PSF losses) as given in the PACS documentation. We have verified this approach by comparing the resulting continuum flux density level to the continuum level of all 25 spaxels combined (which is free of PSF losses and pointing uncertainties). In all cases the agreement is excellent, however the central spaxel alone provides better S/N.

The reduced spectra were next smoothed using a Gaussian kernel of width 0.05 μm (i.e., about half a resolution element) to reduce the noise in the data before the spectral analysis. A spline was fit to the continuum and subtracted from the spectra; these continuum-subtracted spectra were used for the subsequent spectral fitting.

3.3. Spectral Analysis

The profiles of the OH 119.233, 119.441 μm doublet were modeled using four Gaussian components. The fits were carried out using PySpecKit, an open-source, extensible spectroscopic analysis toolkit for astronomy (Ginsburg & Mirocha 2011). This tool uses the Levenberg–Marquardt technique to solve the least-squares problem in order to find the best fit for the observations.¹⁴ Each line of the OH doublet was fitted with two Gaussian components characterized by their amplitude (either negative or positive), peak position, and standard deviation (or, equivalently, FWHM). The separation between the two lines of the doublet was set to 0.208 μm in the restframe (∼520 km s⁻¹) and the amplitude and standard deviation were fixed to be the same for each component in the doublet. Overlap effects between the various components were ignored for simplicity (e.g., the emission component in P Cygni profiles is not affected by the absorption component in the foreground. This is discussed further in the next paragraph). The ¹⁸OH 120 μm doublet is detected in a number of objects but is too far to the red to significantly affect the profiles of the OH 119 μm feature. Likewise, contamination from ¹⁷OH, H₂O⁺, CH, CH⁺, and CO (22–21) is considered unlikely since these lines are either too far to the red or blue of the OH 119 μm feature or undetected even in the high-S/N spectrum of Mrk 231 (see Figure 2(a) of Gonzalez-Alfonso et al. 2013).

We find that three distinct scenarios apply to our data: (1) pure OH absorption, (2) pure OH emission, and (3) P Cygni profiles (our data reveal no new evidence for inverted P Cygni OH profiles; this interesting result is discussed further in Section 4.3). In scenario 1, there is no evidence for any OH emission and each line of the OH doublet is fitted with two absorption components, tracing both the stronger low-velocity component of the outflow and the fainter high-velocity component (e.g., Mrk 231; Gonzalez-Alfonso et al. 2013). Scenario 2 is treated similarly. In this scenario, there is no evidence for any OH absorption feature and two Gaussian components are used to model each of the line in the OH doublet. Scenario 3 is more subtle. We first tried to fit the P Cygni profiles with three or four components per line (1–2 blueshifted absorption components and 1–2 redshifted emission components), but we quickly realized that these fits were underconstrained. We had to settle on fitting these profiles with a single blueshifted absorption component and a single redshifted emission component i.e., no attempt is made to capture the fainter high-velocity component of the outflow. This has the potential to underestimate the outflow velocities derived from this method (making the velocities derived from the absorption feature more positive than they should be). There is also the risk of a degeneracy between the emission and absorption components since in principle they each could be unrealistically strong as long as they cancel each other to reproduce the spectrum. In the radiative transfer models of Gonzalez-Alfonso et al. (2012, 2013), overlap effects between the background emission component and the foreground absorption component prevent this degeneracy, virtually zeroing out the emission component(s) over the range of velocities of the absorption component(s). To approximate these overlap effects, the fit in some cases is further constrained by requiring the absorption feature to not be deeper than the observed absorption. This assumption was made for six objects (F09022−3615, F10565+2448, F12072−0444, F12112+0305, F14348−1447, and F23389+0300), but in the end there are only two cases where the emission component affects the absorption from the blue (119.233 um) line of the doublet: F12112+0305 and F23389+0300. However, even in these objects, the emission contribution is weak at the position of this line and does not affect the measured velocities by more than the typical measurement uncertainties (±50 km s⁻¹; see below), based on comparisons with unconstrained fits. In Section 4.3, we compare the velocities measured in scenario 2 with those measured in scenario 1 and find no significant shift, suggesting that the assumptions made in scenario 2 have no systematic effect on the velocity measurements. Nevertheless, we use caution and distinguish the absorption velocities derived under scenario 2 from those measured in scenarios 1 and 3 in many of the figures presented in Section 4.

These fits were first used to quantify the strength and nature (absorption versus emission) of this feature: (1) the total flux and equivalent width of the OH 119.441 μm line, adding up all of the absorption and emission components, (2) the flux and equivalent width of the absorption component(s) used to fit this line, and (3) the flux and equivalent width of the emission component(s) used to fit this line.

We next characterized the OH profile by measuring a few characteristic velocities from our fits: (1) v₅₀(abs), the median velocity of the fitted absorption profile, i.e., 50% of the absorption takes place at velocities above (more positive than) v₅₀(abs), (2) v₈₄(abs), the velocity above which 84% of the absorption takes place, (3) v₅₀(emi), the median velocity of the fitted emission profile, i.e., 50% of the emission takes place at velocities below (less positive than) v₅₀(emi), and (4) v₈₄(emi), the velocity below which 84% of the emission takes place. We did not extend our fit-based analysis beyond 84% of the absorption/emission profiles because the wings of the OH profiles may not be well captured by our simple two-component fits (especially in the case of P Cygni profiles where a single Gaussian was used to model the absorption (or emission) feature of each line in the doublet). For comparison with S11, we also estimated the terminal outflow velocity, v_max, from the maximum extent of the blueshifted wing of the OH 119.233 μm absorption profile (it was not measured in the case of pure OH emission profiles). This quantity was measured independently by three members of our team to try to reduce possible biases, but it is admittedly more subjective and uncertain than the other velocity measurements, as noted below.

There are a number of potential sources of uncertainties when measuring the equivalent widths, fluxes, and velocities from our data. The wavelength uncertainty depends on mispointing. Any mispointing may also have an effect on the line profiles. We checked for this and found the effect not to apply in our cases. We find instead that the uncertainties on these measurements are dominated by the placement of the continuum. To estimate this source of uncertainty, we re-measure these quantities assuming different continuum shapes (spline versus polynomial) and positions. We find that the typical uncertainties on the equivalent widths and fluxes are ±20% and ±50 km s⁻¹ on the velocities, except for v_max which is much more uncertain (±200 km s⁻¹). For this reason, we give considerably more weight to results based on v₅₀ and v₈₄ in the following discussion.

4. RESULTS

Figure 2 shows the fits to the OH 119 μm profiles. The OH 119 μm parameters derived from these fits are tabulated in Table 2 along with their uncertainties. The average and median values of these parameters are listed at the end of this table. The meaning of each parameter is discussed in Section 3.3 and the notes to Table 2. Note that the fluxes and equivalent widths in that table need to be multiplied by a factor of two when considering both lines of the doublet. In this section we compare these results with the other galaxy properties tabulated in Table 1.

Table 2. Properties of the OH 119 μm Profiles

(Name)	v₅₀ (abs)	v₈₄ (abs)	v_max (abs)	Flux_abs	EQW_abs	v₅₀ (emi)	v₈₄ (emi)	Flux_emi	EQW_emi	EQW_Total
(Name)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(Jy km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(Jy km s⁻¹)	(km s⁻¹)	(km s⁻¹)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
F00509+1225	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	171	363	150.0	−65	−65
F01572+0009	−724::	−892::	−1100::	−32.8::	23::	195::	345::	32.7::	−24::	0::
F05024−1941	−183	−508::	−850::	−83.6	83	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	83
F05189−2524	−327	−574	−850	−510.0	59	189	351	424.1	−51	10
07251−0248	−63	−255	−550	−453.8	56	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	56
F07599+6508	−459	−652	−1000	−47.7	35	315	808	74.1	−53	−19
F08572+3915	−489	−832	−1100	−347.5	124	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	124
09022−3615	−153	−297	−650	−265.5	33	195	339	405.7	−51	−17
F09320+6134	−9	−225	−1200	−5.0	83	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	83
F10565+2448	−267	−489	−950	−2140.7	156	99	273	1099.2	−82	75
F11119+3257	−423	−814:	−1200:	−141.6	124	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	124
Z11598−0114	−153	−477	−800	−295.2	121	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	121
F12072−0444	−69	−321	−1200	−221.9	135	243	393	137.6	−87	51
F12112+0305	−117	−237	−400	−0.8	52	243	447	0.8	−48	2
F12243−0036	111	−21	⋅⋅⋅	−1406.6	63	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	63
F12265+0219	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
F12540+5708	−237	−610	−1500	−5510.7	243	147	357	2956.9	−134	113
13120−5453	−195	−520	−1200	−5037.9	113	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	113
F13305−1739	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
F13428+5608	−201	−495	−750	−679.2	61	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	61
F13451+1232	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	−51	190	203.2	−136	−136
F14348−1447	−291:	−508	−900	−696.1	112	93	273	542.1	−88	25:
F14378−3651	−219	−556	−1200	−1009.4	206	261	423	423.6	−91	119
F14394+5332	−291	−495	−750	−106.9	62	387	520	43.0	−26	37
F15206+3342	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
F15250+3608	189	−21	⋅⋅⋅	−491.6	125	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	125
F15327+2340	21	−153	−700	−22450.3	211	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	−211
F15462−0450	−225:	−459:	−600:	−161.9	80:	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	80:
F16504+0228	−207	−544	−1200	−2126.8	84	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	84
F17207−0014	51	−165	⋅⋅⋅	−4619.4	148	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	148
F19297−0406	−231	−532	−1000	−750.5	119	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	119
19542+1110	−93:	−489:	−700	−331.9	69	243	441	191.2	−41	29:
F20551−4250	−381	−748	−1200	−2.1	70	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	70
F22491−1808	99	3::	⋅⋅⋅	−0.3	25:	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	25:
F23128−5919	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
F23233+2817	−267	−423	−500	−41.3	30	201	339	79.2	−58	−27
F23365+3604	−243	−604	−1300	−1348.8	197	153	321	619.1	−94	107
F23389+0300	−171	−285	−600	−43.4	48	351	634	86.5	−91	−46
PG 1126−041	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
PG 1351+640	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	33	171	12.8	−40	−40
PG 1440+356	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	−105	69	31.1	−43	−43
PG 1613+658	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	135	333	55.8	−80	−80
PG 2130+099	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅

Average	−194	−444	−927	−1605	103	202	401	439	−72	67
Median	−204	−492	−925	−340	84	198	354	197	−70	73

Notes. The uncertainties on the equivalent widths and fluxes (velocities, except v_max) are typically 20% (50 km s⁻¹), unless the value is followed by a colon (when the uncertainties are 20%–50% and 50–150 km s⁻¹, respectively) or a double colon (when the uncertainties are larger than 50% and 150 km s⁻¹, respectively). For v_max, the typical uncertainty is ±200 km s⁻¹. Column 1: galaxy name. Coordinate-based names beginning with F are sources in the IRAS Faint Source Catalog. Column 2: v₅₀(abs) is the median velocity of the fitted absorption profile, i.e., 50% of the absorption takes place at velocities above—more positive than—v₅₀. Column 3: v₈₄(abs) is the velocity above which 84% of the absorption takes place. Column 4: v_max(abs) is the maximum extent of the blueshifted wing of the OH 119.233 μm absorption profile. Column 5: the total integrated flux for the absorption component(s). Column 6: the total equivalent width for the absorption component(s). Column 7: v₅₀ (emi) is the median velocity of the fitted emission profile. Column 8: v₈₄ (emi) is the velocity below which 84% of the emission takes place. Column 9: the total integrated flux for the emission component(s). Column 10: the total equivalent width for the emission component(s). Column 11: the total equivalent width for the sum of the two components for one line of the OH doublet.

Download table as: ASCII Typeset image

4.1. Detection Rate of the OH 119 μm Feature

The 119 μm OH doublet was detected in 37 of the 43 objects in our sample (86%). The OH detection rate shows no obvious dependence on the properties of the objects (Figure 3). A closer inspection of the data indicates that the six objects without OH detection have systematically lower S/N ratios (≲20–30) in the FIR continuum, and therefore less stringent limits on the OH 119 μm equivalent widths, than the others. Two of them are IR-faint QSOs from the OT2 sample (PG 1126−041 and PG 2130+099). The other four are part of our OT1 sample of warm ULIRGs (F13305−1739, F15206+3342, F23128−5919, and F12265+0219 = 3C 273). These six objects are excluded from the rest of the analysis.

**Figure 3.** Histograms showing the overall detection rates of OH 119 μm as a function of (a) the starburst luminosities, (b) the AGN fraction, and (c) the AGN luminosities. The uncertainties on the detection rates assume a binomial distribution. No significant trend is found with any of these parameters (the statistics are poor at low starburst luminosities).
Download figure:
Standard image High-resolution image

4.2. Nature of the OH 119 μm and Mid-infrared Silicate Features

Of the 37 OH 119 μm detections, 17 are seen purely in absorption,¹⁵ 15 show absorption + emission composite profiles, and 5 are purely in emission. Figure 4 shows no obvious trend between the total equivalent width of the OH feature and AGN fraction, starburst luminosity or AGN luminosity, except for the fact that all four objects with AGN fractions α_AGN ≳ 90% show pure OH emission. This result is consistent with the detection of OH 119, 79, and 163 μm emission in the prototypical Seyfert 2 galaxy NGC 1068 (Spinoglio et al. 2005; Hailey-Dunsheath et al. 2012). Preliminary results on the Swift Burst Alert Telescope AGN sample (OT2 $\_$ sveilleu $\_$ 6, PI: S. Veilleux), where OH 119 μm is often seen in emission but log (L_AGN/L_☉) < 12 (M. Meléndez et al. 2013, in preparation), reinforces the idea that the AGN fraction is more important in setting the character (i.e., strength of emission relative to absorption) of the OH feature than the AGN luminosity. These results suggest that dominant AGN in ULIRGs and QSOs provide a favorable environment to produce the OH molecule and excite it (via radiative pumping or collisional excitation) to the ground-state ²Π_3/2 J = 5/2 level. Interestingly, a similar trend with AGN fraction was seen by Teng et al. (2013), when examining the H i 21 cm feature in a smaller set of QUEST ULIRGs and QSOs with the Green Bank Telescope (GBT).

**Figure 4.** Total (absorption + emission) equivalent widths of OH 119 μm (positive values indicate absorption, negative values indicate emission) as a function of (a) the AGN fractions and (b) the AGN luminosities. The meanings of the symbols are as follows: filled black squares, filled red triangles, and open red triangles represent ULIRGs with OH 119 μm seen purely in absorption, composite absorption/emission, and purely in emission, respectively. The OT2 QSOs, where OH 119 μm is purely in emission, are shown separately as open black circles. The vertical dotted lines indicate the positions of the six sources with undetected OH 119 μm (two of these sources have nearly exactly the same AGN fraction). The horizontal gray line marks the null OH equivalent width.
Download figure:
Standard image High-resolution image

Figure 5 compares the strength of the OH 119 μm absorption feature with that of the 9.7 μm silicate absorption feature, a measure of obscuration in ULIRGs. The strengths of the 9.7 μm feature in our sample galaxies are taken from two separate papers that use two different methods to derive this quantity. In V09, $\tau _{\rm 9.7\, \mu{\rm m}}$ is the effective 9.7 μm silicate optical depth measured relative to the sum of the fitted blackbody components (a larger value of $\tau _{\rm 9.7\, \mu{\rm m}}$ implies a deeper silicate absorption feature). In Stierwalt et al. (2013), $S_{\rm 9.7\, \mu{\rm m}}$ is the logarithm of the ratio of the measured flux at the central wavelength of the silicate absorption feature to the local continuum flux (a more negative value of $S_{\rm 9.7\, \mu{\rm m}}$ implies a deeper silicate absorption feature). In both cases, we note that the more obscured systems ( $\tau _{\rm 9.7\, \mu{\rm m}}$ > 4 in V09 or $S_{\rm 9.7\, \mu{\rm m}}$ < −2 in Stierwalt et al. 2013) more often show OH purely in absorption than the less obscured systems. Conversely, systems with strong OH emission either show very weak silicate absorption or silicate in emission. This last category includes three of the four QSOs seen in emission in both OH and silicate (F00509+1225 = I Zw 1, PG 1440+356, and PG 1613+658; PG 1351+640 is the only QSO that is not part of the Spitzer sample of Schweitzer et al. 2008; these objects are not plotted in Figure 5, but would reinforce the trend between FIR OH and MIR silicate).

**Figure 5.** Total (absorption + emission) equivalent widths of OH 119 μm as a function of (a) the optical depths of the silicate feature at 9.7 μm measured by V09 and (b) the depths, on a logarithmic scale, of the 9.7 μm features relative to the local continuum measured by Stierwalt et al. (2013; the absolute value of $S_{\rm 9.7\, \mu{\rm m}}$ is actually plotted here). In both panels, the strength of the silicate absorption feature increases to the right (see Section 4.2 for definitions of the silicate-related quantities). The sign convention for the equivalent widths and the meanings of the symbols are the same as in Figure 4. The vertical dotted lines indicate the positions of the two sources with silicate measurements but undetected OH 119 μm. The horizontal gray line marks the null OH equivalent width.
Download figure:
Standard image High-resolution image

Figures 4 and 5 paint a consistent picture where both the OH 119 μm and 9.7 μm silicate features change character from absorption to emission as the merger progresses and the AGN becomes dominant.

4.3. Distributions of Velocities and Wind Detection Rates

Figure 6 shows the distributions of velocities derived from both the OH absorption and emission line features (v₅₀(abs), v₈₄(abs), v₅₀(emi), v₈₄(emi), and v_max(abs), as defined in Section 3.3 and listed in Table 2). The velocity distributions of the absorption features show a distinct excess at negative values while the opposite is true of the emission features. As mentioned in Section 3.3, all 15 objects with clear composite absorption + emission features show a distinct P Cygni profile i.e., a blueshifted absorption feature accompanied with a redshifted emission feature, an unambiguous signature of outflow. As listed at the bottom of Table 2, the average (median) v₅₀(abs), v₈₄(abs), and v_max(abs) are −194 (−204), −444 (−492), and −927 (−925) km s⁻¹, respectively. Given possible projection effects, these measurements represent lower limits on the true velocities of the OH outflows in these objects.

Figures 6(c) and (d) show v₅₀(abs) and v₈₄(abs) in objects with pure absorption features and P Cygni profiles. The average (median) v₅₀(abs) and v₈₄(abs) are −135 (−183) and −398 (−477) km s⁻¹ for the objects with pure OH absorption profiles versus −262 (−243) and −495 (−495) km s⁻¹ for the objects with P Cygni profiles. Given the measurement uncertainties (±50 km s⁻¹ or larger in some cases; Table 2), we do not consider these differences significant. A Kolmogorov–Smirnov (K-S) between these distributions supports this conclusion. As a case in point, we generally find excellent agreement between our velocities and those of S11 and Gonzalez-Alfonso et al. (2012; the largest discrepancy is seen for the P Cygni profile of F14387−3651, where v₈₄ in S11 is 240 km s⁻¹ more blueshifted than the value reported here). Nevertheless, to remain aware of possible biases, the velocities derived in objects with OH P Cygni profiles will be labeled differently from those derived in objects with OH purely in absorption or emission in all of the following figures.

We follow Rupke et al. (2005b) and conservatively define a wind as having an OH absorption feature with a median velocity (v₅₀) more negative than −50 km s⁻¹. Objects with OH purely in emission are thus excluded by this definition, even though one of them (PG 1613+658¹⁶) shows broad OH profiles extending in excess of ∼500 km s⁻¹, suggestive of non-gravitational motion.¹⁷ As described in Rupke et al., the −50 km s⁻¹ cutoff is used to avoid contamination due to systematic errors and measurement errors in wavelength calibration, line fitting (see Section 3.3), and redshift determination. The few redshifted absorption components detected in our sample have v₅₀ ≲ 50–100 km s⁻¹ (Figure 6), suggesting that this cutoff is reasonable.

Using this conservative definition, we detect a molecular wind in 26 (70%) of the 37 objects in our sample with measurable OH 119 μm feature. This fraction is ∼65% and 62% if the wind velocity threshold is changed to v₅₀ = −75 and −100 km s⁻¹, respectively. The requirements to detect the wind in absorption against the continuum source and for the absorption velocities to be ⩽−50 km s⁻¹ necessarily bias the sample against winds seen edge-on. This wind detection rate is therefore a lower limit and our results may be consistent with all ULIRGs having such molecular winds. Rupke et al. (2005b) came to a similar conclusion regarding neutral winds traced by Na i, although the sample of Rupke et al. (2005b) is different from the present one in that it is made of starburst dominated (U)LIRGs and therefore primarily traces starburst-driven winds. The opening angle of the molecular outflows inferred for the wind detection rate of 70% is ∼145°, assuming all objects in our sample have an outflow.

Figure 7 shows no clear dependence of the wind detection rate on AGN fraction, AGN luminosity, or star formation luminosity, to within the uncertainties of the data (recall that objects with OH feature in pure emission are excluded from this analysis; these objects are discussed in Section 4.7). In contrast, Rupke et al. (2005b) detected Na i winds slightly more frequently in ULIRGs (∼75%) than in LIRGs (∼45%), and attributed this difference to higher collimation of the outflows in LIRGs (see also Chen et al. 2010). The results of Rupke et al. (2005b) refer to starburst dominated LIRGs, which are not adequately sampled by our set of objects, so we cannot verify if the dependence of the wind detection rate on SFR seen in the neutral gas is also present in the molecular gas. More relevant is the sample of Rupke et al. (2005c), which contains several objects in common with the present sample. The lower neutral wind detection rate found among Seyfert 2 ULIRGs in Rupke et al. (2005c; 45% ± 11%) is not seen in OH, although the comparison is necessarily limited to systems in the present sample with α_AGN ≲ 90%. We return to the multi-phase nature of these winds in Section 5.

Only four objects in our sample (F12243−0036 = NGC 4418, 15250+3609, F17207−0014, and F22491−1808) show clear evidence for inflow based on the detection of a redshifted OH absorption feature with median velocity v₅₀(abs) ⩾ 50 km s⁻¹. This paucity suggests that the infalling gas is either rare or subtends a relatively small fraction of 4 π steradians (e.g., planar or filamentary geometry). The fact that clear inverted P Cygni profiles are not seen in our sample of galaxies adds support to this argument. While Gonzalez-Alfonso et al. (2012) have shown that the OH 119 μm feature in NGC 4418 is best fit as a weak inverted P Cygni profile, this conclusion is driven largely by the clear detection of such a profile in [O i] 63 μm. This result serves as a cautionary tale that weak blueshifted OH emission may be missed in our single-line fits. Nevertheless, this does not change the fact that inverted P Cygni profiles with strong blueshifted emission are absent from our sample of objects (although there are a few known exceptions among other ULIRGs, e.g., J. Gracia Carpio et al. 2013, in preparation). This means that most of the infalling gas in ULIRGs lies in front of the FIR continuum source and therefore does not subtend a large solid angle as seen from the source.

4.4. Kinematics versus Host Galaxy Masses and Merger Phases

The velocities of local neutral gas outflows traced by Na i show a tendency to become more negative with increasing host galaxy stellar velocity dispersions and stellar masses (e.g., Rupke et al. 2002, 2005b, 2005c; Martin 2005; Chen et al. 2010). A similar trend is present when using low-ionization tracers in z ∼ 1 star-forming galaxies (Erb et al. 2012; cf. Steidel et al. 2010; Law et al. 2012 at higher redshifts). Figure 8 compares all of the velocities of the OH absorption features derived from our data (regardless of whether they indicate outflow or inflow) with these same host properties. We detect no obvious trend with near-infrared derived stellar velocity dispersions or stellar masses. Perhaps this is not too surprising since our sample covers a considerably narrower range of host galaxy properties than the Na i samples. Indeed, the low-z trends between Na i outflow velocities and host galaxy properties largely disappear when only ULIRGs are considered. Moreover, Rothberg & Fischer (2010) and Rothberg et al. (2013) have recently argued that the near-infrared derived velocity dispersions of luminous IR merger remnants may more closely reflect the properties of the young stellar disk than the global properties of the host galaxy system as a whole. The noticeably smaller stellar velocity dispersions among systems with OH P-Cygni or pure-emission profiles is unexpected, and is not seen in the right panels of Figure 8 where the host stellar masses are considered. This result may be related to the young stellar population bias in near-infrared derived velocity dispersions noted by Rothberg & Fischer.

**Figure 8.** 50%, 84%, and terminal OH outflow velocities as a function of (a)–(c)–(e) the near-infrared derived stellar velocity dispersions measured by Dasyra et al. (2006a, 2006b, 2007; the typical uncertainties on the velocity dispersions are ∼50 km s⁻¹) and (b)–(d)–(f) the H-band absolute magnitudes of the host galaxies from Veilleux et al. (2006) or from Veilleux et al. (2002) assuming H–K' = 0.35 mag (the typical uncertainties on these magnitudes are 0.1 and 0.3 mag, respectively). The meanings of the symbols are the same as in Figure 4. In the lower panels, data points with red circles around them indicate values based on Veilleux et al. (2002). The two data points joined by a segment in the left (right) panels correspond to F13451+1232 E and W (F14394+5332 E and W). Typical uncertainties on v₅₀ and v₈₄ are ±50 km s⁻¹ and ±200 km s⁻¹ on v_max. The smaller symbols have larger uncertainties (values followed by double colons in Table 2).
Download figure:
Standard image High-resolution image

Since nearly all galaxies in our sample are undergoing a galaxy interaction, it is natural to ask whether there is a dependence of the molecular velocities on merger phase, or equivalently, "interaction class" as defined in Veilleux et al. (2002). Given the relatively small number of objects with detected OH absorption features, we simplified the classification into "binaries," where two distinct nuclei have been detected via arcsecond resolution imaging, and "singles," where the two nuclei have coalesced. The results are shown in Figure 9. The OH absorption velocities in binaries are not statistically different from those in singles. This conclusion is the same whether we use v₅₀(abs), v₈₄(abs), and v_max(abs). Recall, however, that IR-faint QSOs are not included in this analysis since their OH feature is in emission. One of the four IR-faint ULIRGs/QSOs with detected OH (PG 1613+658) shows broad OH emission profiles and may therefore qualify as having an outflow. We return to this point below (Section 4.7).

**Figure 9.** 50%, 84%, and terminal OH outflow velocities as a function of interaction class: binaries or single sources. The meanings of the symbols are the same as in Figure 4. The error bars indicate the 1σ scatter in each distribution. There are only two ULIRGs with OH purely in emission: one of them (F13451+1232) is a binary while the other has a single nucleus (F00509+1225 = I Zw 1). Typical uncertainties on v₅₀ and v₈₄ are ±50 km s⁻¹ and ±200 km s⁻¹ on v_max (see Table 2).
Download figure:
Standard image High-resolution image

4.5. Kinematics versus Starburst Luminosities

Plots of the velocities of local neutral gas outflows traced by Na i (on the vertical axis) versus SFRs (on the horizontal axis) show an "envelope" with the most negative measured outflow velocities increasing with SFRs (e.g., Rupke et al. 2005b; Martin 2005; Chen et al. 2010). This trend seems driven largely by the modest outflow velocities measured in dwarf galaxies (Schwartz & Martin 2004) and the large scatter in the outflow velocities measured in the ULIRGs (Rupke et al. 2005b; Martin 2005). A similar trend may be present at z ≳ 1 when considering wind velocities derived from low-ionization species (e.g., Weiner et al. 2009), although the situation likely is more complicated at these higher redshifts (Steidel et al. 2010; Kornei et al. 2012; Martin et al. 2012). Figure 10 shows a similar plot substituting the OH velocities for those based on Na i, covering ∼2 dex in SFRs. Shown on the horizontal axis are the starburst luminosities (∝ SFR) derived using L_SB ≡ (1 − α_AGN)L_BOL (see Section 2.2). No correlation, upper envelope, or increase in scatter within the more luminous starbursts is seen in this figure, regardless of which velocity (v₅₀(abs), v₈₄(abs), or v_max(abs)) is plotted.¹⁸ If anything, there may even be a tendency for the most luminous starbursts to have more positive velocities (slower outflows), but this trend is not found to be statistically significant: the probability, P[null], that the two populations are drawn from the same parent population or are uncorrelated is larger than ∼10%, regardless of the statistical test used for the analysis (see values listed in Table 3). Recall from Section 2.2 that the ULIRGs with the most luminous starbursts in our sample also tend to have smaller AGN fractions (≲40%). As discussed in the next section, these smaller AGN fractions might help explain the slower outflows in these objects. Finally, we warn the readers that this analysis is necessarily limited by the range in SFRs covered by our sample (∼2 dex), which is narrower than that in low-z Na i D studies.

**Figure 10.** 50%, 84%, and terminal OH outflow velocities as a function of the starburst luminosities. The meanings of the symbols are the same as in Figure 4. The four objects with significant inflows (50% OH velocities above 50 km s⁻¹) are not shown here. None of these quantities is significantly correlated with the starburst luminosities (P[null] > 0.05; see Table 3). F12072−0444, indicated by an open red triangle, was not included in the evaluations of these correlations because it is not clear which of the two nuclei is responsible for the OH absorption feature (both nuclei are included in the PACS entrance aperture). Typical uncertainties on v₅₀ and v₈₄ are ±50 km s⁻¹ and ±200 km s⁻¹ on v_max. The smaller symbols have larger uncertainties (values followed by double colons in Table 2).
Download figure:
Standard image High-resolution image

Table 3. Results from the Statistical Analysis for Sources with v₅₀(abs) < 50 km s⁻¹

Parameter	Number of Objects	Break	P_K-S	ρ_s	P_ρ	τ	P_τ	r	P_r
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
log L_SB–v₅₀	28	12	0.11	0.33	0.087	0.23	0.086	0.18	0.37
log L_SB–v₈₄	28	12	0.54	0.19	0.33	0.14	0.30	0.19	0.32
log L_SB–v_max	28	12	0.54	0.21	0.27	0.15	0.25	0.23	0.25
α_AGN–v₅₀	28	50	0.16	−0.47	0.012	−0.34	0.012	−0.47	0.011
α_AGN–v₈₄	28	50	0.16	−0.39	0.038	−0.26	0.049	−0.45	0.015
α_AGN–v_max	28	50	0.47	−0.32	0.10	−0.22	0.097	−0.35	0.069
log L_AGN–v₅₀	28	11.8	0.032	−0.46	0.013	−0.34	0.012	−0.53	0.0036
log L_AGN–v₈₄	28	11.8	0.071	−0.39	0.038	−0.27	0.041	−0.50	0.0070
log L_AGN–v_max	28	11.8	0.44	−0.24	0.23	−0.19	0.16	−0.33	0.087

Notes. Column 1: quantities considered for the statistical test. Column 2: number of objects with no significant inflows (v₅₀(abs) < 50 km s⁻¹). Column 3: value of the break used to define the two distributions for the Kolmogorov–Smirnov (K-S) test. Column 4: null probability of the K-S test. Column 5: Spearman rank order correlation coefficient. Column 6: null probability of the Spearman rank order correlation. Column 7: Kendall's correlation coefficient. Column 8: null probability of Kendall's correlation. Column 9: Pearson's linear correlation coefficient. Column 10: Two-tail area probability of Pearson's linear correlation. Null probabilities less than ∼0.05 (shown in bold-faced characters) indicate statistically significant trends; those with null probabilities between 0.05 and 0.10 are tentative.

Download table as: ASCII Typeset image

4.6. Kinematics versus AGN Fractions

Figure 11 shows the OH velocities plotted against the AGN fractions α_AGN derived from the f₁₅/f₃₀ ratios. A visual inspection of this figure suggests that objects with dominant AGNs (α_AGN ⩾ 50%) have more negative velocities (faster outflows) than objects with dominant starbursts (α_AGN ⩽ 50%). However, a K-S test between the velocity distributions of dominant AGNs and dominant starburst systems indicates that this difference is not statistically significant.¹⁹ The probability P[null] that the two populations are drawn from the same parent population ≃16% for both v₅₀ and v₈₄ and ∼47% for v_max (Tables 3 and 4). A similar negative result is found when α_AGN is varied from 40% to 60% (beyond this range, fluctuations due to the small number of objects in one of the two distributions dominate the statistics making the results unreliable). However, a search for a correlation between the OH kinematics and the AGN fraction produces more positive results with P[null] ≃ 1%–5% (not considering the more uncertain v_max; see Table 3). The fact that the Pearson's test for a linear correlation produces small P[null] indicates that the data favor a smooth transition with AGN fraction from low-velocity to high-velocity outflows rather than a sudden threshold effect.

Table 4. Parameters from Linear Regression Analysis for Sources with v₅₀(abs) < 50 km s⁻¹

X–Y	Y = aX + b
X–Y	a	b
α_AGN–v₅₀	−3.0 ± 1.0	−92.0 ± 40.3
α_AGN–v₈₄	−3.5 ± 1.2	−332.7 ± 61.5
α_AGN–v_max^a	−4.0 ± 1.9	−732.0 ± 98.4
log L_AGN–v₅₀	−212 ± 72	2286 ± 838
log L_AGN–v₈₄	−236 ± 73	2310 ± 867
log L_AGN–v_max^a	−237 ± 100	1891 ± 1185

Notes. X and Y represent the independent and dependent variables, respectively. a and b represent the regression coefficient (slope) and regression constant (intercept), respectively. For the relationship presented in this table we used the ordinary least-square regression of the dependent variable, Y, against the independent variable X, OLS(Y—X). ^aNote from Table 3 that these correlations are tentative so the parameters listed here should be treated with caution.

Download table as: ASCII Typeset image

The preceding analysis used an equal weighting scheme for all data points. The correlations with AGN fractions get slightly stronger if we remove from the analysis the velocities with the largest uncertainties (those indicated by small symbols in Figure 11 and double colons in Table 2).

Note that a positive trend between wind velocities and AGN fractions does not necessarily imply that the AGN is the driver of these outflows. Once a significant fraction of the obscuring material has been swept away, either by an AGN- or starburst-driven wind, we expect to be able to see down the core of the ULIRG more easily. It therefore becomes easier to observe the central high-velocity outflowing material and detect the AGN via optical–MIR diagnostics. Orientation effects may also be at play in some cases: both the AGN and central high-velocity wind are more easily detectable in systems where the dusty torus/disk are seen face-on. To test whether the AGN plays a role in driving these outflows, we need to compare the outflow kinematics with the AGN luminosities; this is done in the next section.

4.7. Kinematics versus AGN Luminosities

Figure 12 displays the OH outflow velocities versus the AGN luminosities, L_AGN ≡ α_AGNL_BOL (Section 2.2). This figure shows distinctly more negative velocities in systems with larger AGN luminosities. This trend is more evident with v₅₀ and v₈₄ than with the more uncertain v_max. The null probability for a correlation between the OH velocities (v₅₀ or v₈₄) and log L_AGN range from ∼0.4% to ∼4%, depending on the statistical test (Tables 3 and 4). A visual inspection of panels (a) and (b) of Figure 12 also suggests a steepening of the relation between outflow velocities and log L_AGN at high AGN luminosities: systems with v₅₀ ⩽ −400 km s⁻¹ all have AGN luminosities larger than ${\rm log} (L_{\rm AGN}^{\rm break}/L_\odot) = 11.8 \pm 0.3.$ The factor of two uncertainty on this luminosity break takes into account the measurement uncertainty on this break visually estimated from the panels of Figure 12 (estimated to be ±0.2) and the uncertainties (±0.3) on the values of α_AGN and L_AGN for our objects. To verify the value of this break more quantitatively, K-S tests were carried out between the velocity distributions of objects with log (L_AGN/L_☉) > N and those with log (L_AGN/L_☉) < N, where N was varied from 11.7 to 12.0 (beyond this range, fluctuations due to the small number of objects in one of the two distributions dominate the statistics, making the results unreliable). The results from these tests confirm that the most significant difference between the two distributions occurs when N = 11.8 (P[null] ≃ 3% for v₅₀ and 7% for v₈₄; Tables 3 and 4). This break is not statistically significant when using v_max. These results emphasize the need to use more robust values of the outflow velocities, such as v₅₀ and v₈₄, when carrying out these analyses.

Once again, the preceding analysis used an equal weighting scheme for all data points. The conclusions remain the same if we remove from the analysis the velocities with the largest uncertainties (those indicated by small symbols in Figure 12 and double colons in Table 2), although in that case the significance of the correlations with AGN luminosities become virtually the same as with AGN fractions.

Recall that none of the four ULIRGs/QSOs with AGN fraction α_AGN ≳ 90% and detected OH 119 μm features, including all of the OH-detected OT2 IR-faint quasars, shows OH absorption components. These objects thus do not appear in Figures 10–12 and are excluded from our statistical analysis (Table 3). The generally modest widths of the OH emission profiles in these objects (PG 1613+658 is the only exception) do, however, suggest that the winds, if present, have modest velocities. This result seems inconsistent with the trend of increasing outflow velocities with increasing AGN luminosities since all four objects with AGN fraction α_AGN ≳ 90% have AGN luminosities near or above 10^11.8 L_☉. Perhaps we are seeing the end of the fast wind phase predicted by some models (Narayanan et al. 2008; Hopkins et al. 2009).

5. DISCUSSION

Our kinetic analysis of the OH 119 μm velocities in nearby ULIRGs and QSOs described in Section 4 shows trends of increasing molecular wind velocities (becoming more negative) with increasing AGN fractions and luminosities. While the former may be explained as an obscuration effect where both the AGN and central high-velocity outflowing material are more easily detectable when the dusty material has been swept away or is seen more nearly face-on, the trend with AGN luminosity suggests that the AGN is actually playing a role in driving the fast winds.²⁰ Our data favor the existence of a nonlinear (steepening) relation between outflow velocities and the logarithms of the AGN luminosities. Above

$\begin{eqnarray} &&{\rm log} \left(L_{\rm AGN}^{\rm break}/L_\odot\right) = 11.8 \pm 0.3, \end{eqnarray} \tag{ 3 }$

the AGN seems to play a dominant role in driving the outflowing molecular gas in the objects of our sample. This luminosity break coincides approximately with the minimum AGN luminosity traditionally used to define a quasar (M_B ⩽ −23).

This conclusion only applies to objects with absorption signatures of molecular outflows. Fully coalesced ULIRG/QSO mergers with AGN fractions α_AGN ≳ 90%, including all of the OH-detected OT2 IR-faint quasars where OH is seen purely in emission, generally show modest line widths and thus modest outflows despite their large AGN luminosities.

The results from S11 on 6 ULIRGs in common with our sample indicate that the OH outflows are not only fast but also powerful, capable of displacing a significant fraction of the entire molecular interstellar medium (ISM) of the host galaxies. Given our methods of selection, there is no reason to believe that the objects studied here are any different from those in S11. Since the molecular ISM is the fuel for future star formation and galaxy build-up, the quasar-driven outflows detected here may therefore have a significant effect on the evolution of their hosts, potentially quenching star formation in their centers. The modest OH emission line widths in the IR-faint quasars perhaps indicate that AGN feedback subsides once the quasar has poked a substantial hole through its natal "cocoon," although this is based on only four objects.

There has been considerable theoretical support in recent years for the idea of quenching star formation in galaxies via AGN-driven winds (e.g., Narayanan et al. 2008; Hopkins et 2009, and references therein), but relatively little direct observational evidence for it, until recently (e.g., Feruglio et al. 2010; S11; RV11; RV13). This last paper also points out the possibility of increasing outflow velocities at higher AGN luminosities. Using detailed IFU data on a sample of six ULIRGs, all from the QUEST sample, RV13 find that fast AGN-driven neutral (Na i D) and/or ionized (Hα) winds are only present in systems with log $(L_{\rm AGN}^{\rm break}/L_\odot) \ge 11.7$ . This AGN luminosity, while only tentatively identified in that paper since it is based on a very small sample, is consistent to within the uncertainties to that inferred from the larger set of Herschel data. This suggests that the conditions to drive fast neutral/ionized outflows are approximately the same as for the molecular outflows.

In this context, it is instructive to compare the OH velocities measured here with the Na i velocities measured by Rupke et al. (2005a, 2005b, and 2005c) and H. B. Krug et al. (2013, in preparation) from integrated spectra and RV13 from IFU data (Figure 13). Also shown in this figure for comparison are the MIR fine-structure line outflow velocities deduced by Spoon & Holt (2009) from Spitzer integrated spectra, the Hα velocities measured by RV13 from the IFU data, and the H i 21 cm velocities measured by Teng et al. (2013) from GBT integrated spectra. While the various velocities are measured over different scale (≲200 pc for the OH velocities according to S11 and ≲1−2 kpc for the neutral and ionized gas velocities), general trends are observed: the outflow velocities of the molecular gas are often comparable to the velocities of the neutral gas, but they are generally smaller (more positive) than the velocities of the ionized gas (panels (b) and (c) in Figure 13). The similarity between the OH and Na i velocities is particularly evident in panel (b) of Figure 13 where the sometimes very large nuclear Na i and Hα outflow velocities, produced on sub-parsec scales (e.g., Mrk 231; Veilleux et al. 2013; Figure 13(a)), are excluded from the means. This panel also reiterates the result found by RV13 that the ionized gas is often, but not always, the fastest component of the outflow. While dust obscuration may explain some of these variations (especially for the highly ionized gas found near the center of AGN-dominated ULIRGs; e.g., blue wings in [Ne v] emission profiles of Spoon & Holt 2009), it likely does not explain all of the spatial variations found in the IFU data.

**Figure 13.** Multi-phase comparison of the outflow velocities. The 50% OH outflow velocities are compared with (a) the 50% Na i D outflow velocities measured by Rupke et al. (2005a, 2005b, 2005c) and H. B. Krug et al. (2013, in preparation) from integrated spectra (in F05024−1941, F07599+6508, and 12540+5708, there are multiple Na i components seen at different velocities), (b) the spatially averaged 50% Na i D (encircled) and Hα outflow velocities measured by RV13 outside of the nucleus using the IFU data, (c) the outflow velocities measured by Spoon & Holt (2009) from the emission profiles of [Ne iii] 15.5 μm (encircled) and [Ne v] 14.3 μm in integrated *Spitzer* spectra, and (d) the H i 21 cm outflow velocities measured by Teng et al. (2013) from GBT data. The meanings of the symbols are the same as in Figure 4. The red dotted line is the 1:1 line.
Download figure:
Standard image High-resolution image

These multi-phase comparisons may help constrain the processes involved in accelerating the material entrained in the wind. In most cases, the acceleration a = F/m is expected to depend strongly on the ISM phase. For instance, in the case of ram pressure acceleration of clouds in a fast wind or a jet (e.g., Strel'nitskii & Sunyaev 1973; Chevalier & Clegg 1985), the force on a spherical cloud of mass M_cloud and density ρ_cloud is expected to scale as the area of the cloud or (M_cloud/ρ_cloud)^2/3 and thus we expect the acceleration to scale as $(M_{\rm cloud} \rho _{\rm cloud}^2)^{-1/3}$ (we neglect gravity). Naively, one would expect the molecular clouds to be considerably denser and more massive, on average, than the ionized gas clouds so one would expect the molecular material to always move at lower velocities. This simple-minded argument may be relaxed if each accelerated cloud has a multi-phase structure (e.g., ionized skin or halo surrounding a neutral/molecular core) and thus contributes all at once to the molecular, neutral, and ionized phases of the outflow. But the significant difference in velocities often seen between the ionized and molecular gas phases (Figure 13) indicates that this situation cannot apply to the bulk of the outflowing material since the difference in velocities between the various phases would imply severe shearing/shredding and very short cloud lifetimes (∼R_cloud/Δv ∼ 10⁴(R_cloud/pc)(100 km s⁻¹/Δv) years, where R_cloud is the cloud radius and Δv is the difference in velocities between the various gas phases).

Another possible source of acceleration for these winds is UV–IR radiation pressure. It has the distinct characteristic of being less efficient at accelerating ionized material than dusty neutral/molecular material due to the different sources of opacity. Within the dust sublimation radius, $R_{\rm subl} \sim L_{46}^{1/2}$ pc (where a dust sublimation temperature T_subl ≃ 1200 K was assumed), the dominant source of opacity is electron scattering and the corresponding luminosity condition to drive a (spherical) wind is the well-known Eddington luminosity. Estimates of the BH masses M_BH in our sample galaxies range from ∼0.5–5 × 10⁸ M_☉ (e.g., Vestergaard & Peterson 2006; Veilleux et al. 2009a), corresponding to Eddington luminosities L_Edd = 1.3 × 10⁴⁶ (M_BH/10⁸) = 0.65–6.5 × 10⁴⁶ erg s⁻¹. Systems near the AGN luminosity break identified in our data are therefore accreting at sub-Eddington values $\Gamma \equiv L_{\rm AGN}^{\rm break} / L_{\rm Edd}$ ∼ 0.04–0.4.

Outside the dust sublimation radius, the much larger dust opacity prevails and UV–IR radiation pressure from the AGN may drive significant winds. In this situation where dust opacity dominates, Murray et al. (2005) have derived a handy expression for the limiting, Eddington-like luminosity above which momentum deposition from the quasar (and/or starburst) is presumed to be enough to clear most of the gas out of galactic nuclei (i.e., this is a criterion to "blow away" the gas rather than just a "blow out" in the nomenclature of De Young & Heckman 1994):

$\begin{eqnarray} L_M &\simeq & (4 f_g c)/G) \sigma ^4\nonumber\\ &\simeq& 3 \times 10^{46} (f_g / 0.1) (\sigma / 200\ {\rm km\ s}^{-1})^4\, {\rm erg}\ {\rm s}^{-1}, \end{eqnarray} \tag{ 4 }$

where f_g is the gas fraction of the host and σ is the stellar velocity dispersion, a measure of the depth of the gravitational potential well. In our sample galaxies, f_g ≃ 0.1 on average, although there is considerable scatter among ULIRGs (e.g., Downes & Solomon 1998; Combes et al. 2013, and references therein), while σ ranges from ∼120 to ∼280 km s⁻¹ with an average value of ∼200 km s⁻¹ (see Figure 8; as mentioned in Section 4.4, this value can be considered a lower limit; Rothberg & Fischer 2010; Rothberg et al. 2013). The steep dependence of L_M on σ implies that $L_{\rm AGN}^{\rm break}$ ≃ 2%–100% of L_M, with an average value of ∼10%, i.e., the AGN luminosity break above which the AGN plays a dominant role in driving the outflowing molecular gas in local ULIRGs is generally lower than the Eddington-like luminosity derived by Murray et al. (2005).

The fact that $L_{\rm AGN}^{\rm break}$ is lower than L_M is not surprising since the latter quantity is the maximum luminosity an object may have before blowing away all of the gas from a galaxy via radiation pressure on dust. In addition, the theoretically derived L_M is subject to a number of issues that may overestimate its value in a real system: (1) First, it is important to note that the expression for L_M is an order-of-magnitude estimate of the actual value. A number of simplifying assumptions were made to derive this expression, e.g., the ISM is approximated as an isothermal sphere of gas with constant gas fraction. The ISM in ULIRGs is undoubtedly much more complex. Disky ISM distributions, such as those presumed to exist in a number of ULIRGs (e.g., Downes & Solomon 1998), offer less resistance to winds along the disk rotation axis for a given L_AGN. A good example is the wide-angle galactic wind in Mrk 231 (RV11, RV13), which is inferred to be oriented perpendicular to the central molecular disk, yet it entrains a substantial fraction of the gas in the nuclear region (i.e., it is closer to a "blow away" event than just a "blow out," so L_M is relevant here). In situ formation of the accelerated clouds (Faucher-Giguère et al. 2012) and/or energy-conserving acceleration due to inefficient cooling in the shocked winds of the AGN (Faucher-Giguère & Quataert 2012) may also effectively reduce the value of L_M. (2) The expression for the maximum luminosity only takes into account the radiation pressure from the AGN. Momentum injection from loosely collimated jets (e.g., Teng et al. 2013) would reduce the AGN luminosity needed for efficient molecular outflow. (3) While not dominant in systems with AGN luminosities above $L_{\rm AGN}^{\rm break}$ , the powerful starbursts known to exist in these systems also contribute energy and momentum to driving the winds (including possibly significant contributions from cosmic rays). We have neglected the starburst contributions in the present treatment. This would relax the requirements on the AGN, effectively reducing L_M from the AGN by a factor of up to ∼2.

Observationally, the gas mass fraction f_g is difficult to determine in ULIRGs and depends on the uncertain CO-to-H₂ conversion factor. The value of f_g likely varies considerably from one galaxy to the other in our sample. This is particularly true of the IR-faint QSOs, although they are still undoubtedly gas-rich (e.g., Scoville et al. 2003). Given the uncertainties on f_g and σ, the large galaxy-to-galaxy variations in these quantities, and the steep dependence of L_M on σ, it is in fact not surprising to find that our data (Figure 12) favor a trend, albeit nonlinear, between OH velocities and AGN luminosities rather than a clear-cut AGN luminosity threshold for powerful molecular outflows. Moreover, the exact shape of this outflow velocity–AGN luminosity relation likely applies only to local gas-rich ULIRGs. Comparisons with other types of objects with different f_g, σ, and redshifts have to be done with great care. The gas mass fraction of local giant spirals is ∼7%–10% (Leroy et al. 2008; Saintonge et al. 2011), similar to that of local ULIRGs, but there is mounting evidence that it increases with redshifts, reaching ∼34% ± 5% and 44% ± 6% among "main-sequence" galaxies at z ∼ 1.2 and ∼2.3, respectively (Tacconi et al. 2010; Daddi et al. 2010). New observations also indicate that the gas mass fraction increases with redshift among ULIRGs, by a factor of 3 ± 1 from z = 0 to 1 (Combes et al. 2013). These results fit well with the large gas mass fractions (∼40%) previously derived in z ∼ 2–3.4 submillimeter galaxies (Tacconi et al. 2006). For objects with the same stellar velocity dispersion, L_M at high redshifts is thus expected to be three to four times larger than in local objects. Future ALMA observations should be able to verify this assertion.

6. CONCLUSIONS

Herschel/PACS data from the SHINING guaranteed-time key project were combined with data from cycles 1 and 2 open-time programs to carry out a systematic search for molecular outflows in a sample of 43 nearby ULIRGs and QSOs using OH 119 μm. The main results from this surveys are:

1.
The OH 119 μm feature was detected in 37 (86%) of the 43 sample galaxies.
2.
The nature of the OH 119 μm feature (emission, absorption, or both) depends on the properties of the galaxies: OH emission is stronger relative to OH absorption in quasar-dominated systems, becoming completely dominant (i.e., pure OH emission) in objects where the quasar contributes more than ∼90% to the bolometric luminosity (called the AGN fraction for short). This behavior is similar to that of the silicate 9.7 μm feature studied with Spitzer and the H i 21 cm feature measured with GBT in a subset of ULIRGs and QSOs.
3.
Unambiguous evidence for molecular outflows, based on the presence of absorption features with median velocities more blueshifted than −50 km s⁻¹, is seen in 26 (70%) of the 37 targets with detected OH 119 μm. Given possible projection effects, this wind detection rate is a lower limit and is consistent with all ULIRGs having molecular winds with an average opening angle ∼145°. In contrast, absorption features with median velocities above +50 km s⁻¹, indicative of inflow, are detected in only four objects, none of which show clear inverted P Cygni profiles. This result suggests that molecular inflows are either rare or subtend a relatively small fraction of 4 π steradians (e.g., planar or filamentary geometry).
4.
Typical median outflow velocities are ∼−200 km s⁻¹, but terminal outflow velocities of ∼−1000 km s⁻¹ are detected in several objects. Deprojected outflow velocities will be even larger. The kinematics of these molecular outflows do not show any obvious dependence on the properties of the host galaxies (e.g., near-infrared derived stellar velocity dispersions, stellar masses) and SFRs, although our sample covers a range of properties that is narrower than that of low-z non-molecular wind studies, where trends have been observed with these quantities.
5.
There is a tendency for objects with dominant and luminous AGN to show more blueshifted OH velocities (faster outflows). Our data favor a steepening of the relation between outflow velocities and the logarithms of the AGN luminosities above log $(L_{\rm AGN}^{\rm break} / L_\odot) = 11.8 \pm 0.3$ , although better statistics are needed, particularly at lower starburst and AGN luminosities, to confirm this nonlinearity. While the trend between outflow velocities and AGN fractions may be explained as an obscuration effect, where both the AGN and central high-velocity outflowing material are more easily detectable when the dusty material has been swept away or is seen more nearly face-on, the dependence of the velocities on the AGN luminosities suggests that the AGN is playing a role in driving the fast winds. Our results emphasize that the use of the terminal velocities as the sole measures of outflow velocities should be avoided due to the much larger uncertainties on these measurements.
6.
The OH emission profiles in three of the four late-stage ULIRG/QSO mergers with AGN fractions above ∼90% are narrow, despite their large AGN luminosities, suggesting that the winds in these objects, if present, often have modest velocities. Quasar feedback thus seems to subside among mergers that have cleared a path through the dusty circumnuclear material.

We thank the referee for suggesting changes that helped improve this paper. Support for this work was provided by NASA through Herschel contracts 1427277 and 1454738 (S.V. and M.M.) and contracts 1364043, 1435724, and 1456609 (J.F.). S.V. also acknowledges support from the Alexander von Humboldt Foundation for a "renewed visit" to Germany following up the original 2009 award, and thanks the host institution, MPE Garching, where a portion of this paper was written. E.G.-A. is a Research Associate at the Harvard-Smithsonian Center for Astrophysics, and thanks the support by the Spanish Ministerio de Economía y Competitividad under projects AYA2010-21697-C05-0 and FIS2012-39162-C06-01. Basic research in IR astronomy at the Naval Research Laboratory is funded by the U.S. Office of Naval Research. This research made use of PySpecKit, an open-source spectroscopic toolkit hosted at http://pyspeckit.bitbucket.org. This work has made use of NASA's Astrophysics Data System Abstract Service and the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

FAST MOLECULAR OUTFLOWS IN LUMINOUS GALAXY MERGERS: EVIDENCE FOR QUASAR FEEDBACK FROM HERSCHEL

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION