ABSTRACT
Mass outflows driven by stars and active galactic nuclei (AGNs) are a key element in many current models of galaxy evolution. They may produce the observed black-hole–galaxy mass relation and regulate and quench both star formation in the host galaxy and black hole accretion. However, observational evidence of such feedback processes through outflows of the bulk of the star-forming molecular gas is still scarce. Here we report the detection of massive molecular outflows, traced by the hydroxyl molecule (OH), in far-infrared spectra of ULIRGs obtained with Herschel-PACS as part of the SHINING key project. In some of these objects the (terminal) outflow velocities exceed 1000 km s−1, and their outflow rates (up to ∼1200 M☉ yr−1) are several times larger than their star formation rates. We compare the outflow signatures in different types of ULIRGs and in starburst galaxies to address the issue of the energy source (AGN or starburst) of these outflows. We report preliminary evidence that ULIRGs with a higher AGN luminosity (and higher AGN contribution to LIR) have higher terminal velocities and shorter gas depletion timescales. The outflows in the observed ULIRGs are able to expel the cold gas reservoirs from the centers of these objects within ∼106–108 years.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
Gas-rich galaxy merging may trigger major starbursts, lead to the formation of elliptical galaxies, and account for the formation and growth of supermassive black holes (BHs; 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 are 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 quench the growth of both the BH and spheroidal component and explain the tight BH–spheroid mass relation (e.g., Di Matteo et al. 2005; Murray et al. 2005). These winds are purported to inhibit star formation in the merger remnants ("negative mechanical feedback"; e.g., Veilleux et al. 2005, for a review), and to create a population of red gas-poor ellipticals, thereby explaining the bimodal color distribution observed in large galaxy surveys (e.g., Kauffmann et al. 2003).
Finding observational evidence of such feedback processes in action is one of the main challenges of current extragalactic astronomy. While outflows have been observed frequently in many starbursts and QSOs, so far they have been detected mostly in the ionized and neutral atomic gas component. To inhibit star formation in the host galaxy, outflows have to affect the molecular gas out of which stars form. Few detections of molecular outflows have been reported so far (e.g., Baan et al. 1989; Walter et al. 2002; Sakamoto et al. 2009). In this Letter, we demonstrate that far-infrared molecular spectroscopy with Herschel-PACS of (ultra-)luminous infrared galaxies is providing a breakthrough in identifying and analyzing massive molecular outflows. Our recent OH-absorption observations have revealed a >1000 km s−1 molecular outflow in the closest quasar known, Mrk 231 (Fischer et al. 2010). Independent, spatially resolved CO-emission observations of Mrk 231 with the IRAM/Plateau de Bure (PdB) mm-wave interferometer (Feruglio et al. 2010) have confirmed this outflow with inferred mass outflow rates of ∼300–2200 M☉ yr−1, significantly larger than the current star formation rate (SFR ∼ 100 M☉ yr−1) in the host galaxy. We now show that molecular outflows are indeed a common phenomenon in many of the luminous major mergers in our sample, reaching outflow velocities of 1000 km s−1 and outflow rates up to ∼1000 M☉ yr−1 in some of them.
2. OBSERVATIONS AND DATA REDUCTION
The data presented here are part of the Herschel guaranteed time key program SHINING, a study of the far-infrared properties of the interstellar medium in starbursts, Seyfert galaxies, and infrared luminous galaxies. Here we present velocity-resolved line profiles of the OH 79 μm cross ladder, ground-state doublet for a small subset of our ULIRG sample, and a comparison starburst galaxy. The observations were obtained with the PACS far-infrared spectrometer (Poglitsch et al. 2010) on board Herschel (Pilbratt et al. 2010). Our data set also contains other OH transitions, but the 79 μm line was observed first and is common to all objects, as it can be obtained simultaneously (in second order) with the [C ii] 157 μm line (in first order). At 79 μm the PACS resolution is ∼140 km s−1. The data reduction was done using the standard PACS reduction and calibration pipeline (ipipe) included in HIPE 5.0. However, for the final calibration we normalized the spectra to the telescope flux (which dominates the total signal, except for NGC 253) and re-calibrated it with a reference telescope spectrum obtained from dedicated Neptune observations during the Herschel performance verification phase. All of our objects (except NGC 253) 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 signal-to-noise ratio (S/N). We note for completeness that for NGC 253 the total OH 79 μm line profile summed over all 5 × 5 spaxels yields emission, consistent with the Infrared Space Observatory Long-Wavelength Spectrometer observations by Bradford et al. (1999).
In a next step we have performed a continuum (spline) fit. Due to the limited wavelength coverage these fits are somewhat subjective. To help define continuum points and potential additional spectral features (such as the H2O absorption line at 78.74 μm, indicated with an arrow in Figure 1), we have used our full range spectra of Arp 220 and NGC 4418. These two sources will be analyzed in detail in forthcoming papers, but preliminary data points for Arp 220 are included in Figures 2 and 3. We note here that NGC 4418 shows signatures of an inflow.
Figure 1. Observed PACS spectra (continuum-normalized) of the OH transition at 79 μm (gray). Overplotted are the low-velocity (dotted) and high-velocity (dashed) fit components and the total fit (solid). The arrow indicates the rest position of H2O 423–312. The dash-dotted line for IRAS 14378 shows the observed spectrum of the OH transition at 119 μm for this object.
Download figure:
Standard image High-resolution image
Figure 2. Maximum outflow velocities (terminal velocities) as a function of star formation rate (upper panel) and AGN luminosity (lower panel). The asterisk denotes NGC 253 and the triangle denotes Arp 220.
Download figure:
Standard image High-resolution image
Figure 3. Upper panel: the ratio of the mass outflow rate to the SFR vs. SFR; lower panel: depletion timescale vs. AGN luminosity. Symbols are as in Figure 2.
Download figure:
Standard image High-resolution image
3. TARGETS
For this first study of outflow signatures in our data we use a sub-sample that is mainly constrained by the observing schedule of Herschel, but that covers a broad range of AGN and starburst activity, including a starburst template (NGC 253), a cold, starburst-dominated ULIRG (IRAS 17208−0014), warm ULIRGs (S25/S60 > 0.1) and/or ULIRGs with strong AGN contributions (Mrk 231, IRAS 13120−5453, IRAS 14378−3651), and a heavily obscured ULIRG (IRAS 08572+3915), which hosts a powerful AGN (e.g., Veilleux et al. 2009, hereafter V09).
4. RESULTS AND DISCUSSION
Figure 1 shows the (continuum-normalized) OH 79 μm line spectra for all objects. For NGC 253, we show the central spaxel only. The Mrk 231 spectrum is taken from Fischer et al. (2010) and is repeated here for completeness. In all cases we detect P-Cygni profiles typical of outflows, with blueshifted absorption and redshifted emission features, of more than 1000 km s−1 in some cases. For comparison with the literature, we list in Table 1 various measures of outflow velocities (relative to the system velocity of the blue component of the OH doublet): the peak velocity (vpeak), the maximum (terminal) velocity (vmax), and 
for which 85% of the outflowing gas has lower (absolute) velocities. The uncertainties of these velocities are dominated by the uncertainties in the continuum definition, the S/N in the spectra, and the spectral resolution. We estimate an overall error of ±150 km s−1. The S/N in the IRAS 14378−3651 spectrum is relatively low, but the high terminal velocity is confirmed by the OH 119 μm transition (see Figure 1).
Table 1. Target Properties, Outflow Rates, and Outflow Velocities (1σ Uncertainties in Parentheses)
| Source | SFR | αa | LAGN | Mgasb |  c |
vpeakd |  e |
vmaxf |
|---|---|---|---|---|---|---|---|---|
| (M☉ yr−1) | (%) | (1011 L☉) | (109 M☉) | (M☉ yr−1) | (km s−1) | (km s−1) | (km s−1) | |
| Mrk 231 | 101 (15) | 71 (11) | 28 (4) | 4.2 (1.3) | 1190+4700 − 890 | −600 | −660 | −1170 |
| IRAS 08572+3915 | 42 (6) | 72 (11) | 12 (2) | 1.3 (0.4) | 970+2900 − 730 | −700 | −740 | −1260 |
| IRAS 13120−5453 | 168 (25) | 9 (1.4) | 1.8 (0.3) | 5.8 (1.7) | 130+390 − 95 | −520 | −600 | −860 |
| IRAS 14378−3651 | >79 | <45 | <7.2 | 4.2 (1.3) | 740+2200 − 550 | −800 | −860 | −1170 |
| IRAS 17208−0014 | 274 (41) | 11 (1.7) | 3.4 (0.5) | 12.2 (3.7) | 90+270 − 65 | −100 | −170 | −370 |
| NGC 253 | 1.7 (0.3) | 0 | 0 | 0.7 (0.2) | 1.6+4.8 − 1.2 | −75 | −130 | −280 |
Notes. Estimated uncertainty for all velocities: ±150 km s−1. aFraction of the AGN contribution to Lbol, where Lbol = 1.15 × LIR. bGas mass (taken from Graciá-Carpio et al. 2011). cMass outflow rate (see the footnote of Table 2). dPeak velocity of the blueshifted high velocity component (relative to systemic velocities). eVelocity for which 85% of the outflowing gas has lower (absolute) velocities. fTerminal velocity.
Download table as: ASCIITypeset image
These high OH outflow velocities may be the long-sought conclusive evidence of powerful mechanical feedback from vigorous star formation and/or accreting central BHs. The possible feedback and outflow mechanisms (e.g., winds from supernovae, radiation pressure) are debated in the literature (see, e.g., the review by Veilleux et al. 2005). It is not clear if such mechanisms could indeed be sufficient to power outflows that are strong enough to significantly affect the host galaxy and to actually quench the star formation in these objects. It is also unclear from the models whether it is possible to distinguish AGN-driven outflows from stellar-driven outflows observationally (see, e.g., Hopkins & Elvis 2010). In the following, we adopt an empirical approach with our new data.
4.1. Are the Strong Outflows We Observe Driven by the AGN Rather than by the Star Formation in These Objects?
Rupke et al. (2005a, 2005b, 2005c) and Krug et al. (2010) have studied large samples of AGN and star-forming galaxies in neutral gas (blueshifted optical Na i D 5890, 5896 Å absorption features). They found that, for fixed SFR, ULIRGs with higher AGN fractions have higher neutral gas outflow velocities, reaching velocities well above 1000 km s−1 in some broad-line AGN (see also, e.g., Heckman et al. 2000; Martin 2005, 2006; Thacker et al. 2006). Theoretical models predict that supernovae-driven outflows cannot reach velocities higher than 500–600 km s−1 (e.g., Martin 2005; Thacker et al. 2006). Predictions of vmax from models of outflows driven by radiation pressure from a starburst (or AGN) are, however, less certain. The terminal velocity we measure in the OH outflow of our starburst template NGC 253 is ∼300 km s−1. IRAS 17208−0014, a starburst-dominated ULIRG with little AGN contribution, has only a slightly higher vmax (370 km s−1). In significant contrast to this, the terminal velocities of the OH outflow in the two AGN-dominated ULIRGs (Mrk 231 and IRAS 08572+3915) are well above 1000 km s−1. Thus, the OH outflow velocity could be a very promising tool to distinguish AGN-driven outflows from starburst-driven outflows, with AGN-dominated outflows reaching much higher velocities.
In Figure 2, we compare the terminal outflow velocities to the SFRs and AGN luminosities of our objects. SFRs are calculated from the IR luminosities (based on the calibration of Kennicutt 1998, but using a Chabrier initial mass function (IMF) rather than a Salpeter IMF), i.e., SFR = (1 − α) × 10−10LIR, applying AGN correction factors α (the fraction of the contribution from the AGN to Lbol, where Lbol = 1.15 × LIR as in V09). These factors are individually derived from Spitzer mid-IR diagnostics provided by the QUEST programme (V09). For IRAS 13120−5453 and IRAS 14378−3651 (which are not part of V09), we calculated AGN fractions using method 1 from V09, based on the observed [O iv]/[Ne ii] ratio (and the observed upper limit for IRAS 14378−3651) from Farrah et al. (2007). AGN luminosities are then calculated as LAGN = α × Lbol. AGN luminosities and SFRs are listed in Table 1.
For star-formation-driven outflows, one might expect the outflow velocities to scale with the SFRs (e.g., Tremonti et al. 2007). We do not see such a correlation (Figure 2, upper panel). Instead, we see a rough correlation of vmax with LAGN (Figure 2, lower panel), consistent with the idea of the high velocity outflows being powered mostly by (radiation pressure from) the AGN.
Very energetic outflows have been found in radio galaxies at high z (e.g., McCarthy et al. 1996; Best et al. 2000; Nesvadba 2009). Radio jets could in principle be driving the outflows we see. However, nearly all of the objects presented here are radio quiet, so this energy source can safely be assumed to be negligible. Only in Mrk 231 does a radio jet contribute to the Na i D outflow (Rupke & Veilleux 2011).
4.2. Does the Outflow Carry Sufficient Molecular Gas to Remove the Star Formation Fuel and Actually Quench the Star Formation?
To compute mass outflow rates (
) for comparison with SFRs, we have modeled the observed spectra using the radiative transfer code described in González-Alfonso & Cernicharo (1999). The outflow is modeled as concentric expanding shells around a nuclear continuum source, allowing for each source one to three components with different velocity gradients and distances to the central source. Such geometry reproduces quite naturally the redshifted emission (which is produced in the receding cocoons) and the blueshifted absorption (from the approaching parts).
The density profile for each velocity component is determined through mass conservation assuming a stationary outflow. With the exception of NGC 253, where we fit the nuclear continuum emission from the observed PACS spectrum, the nuclear FIR in all other sources is modeled as the warm component fitted in Mrk 231 (González-Alfonso et al. 2010). Free parameters are the inner and outer radii, the velocity field of each velocity component, the OH density at the inner radius, the covering factor of the continuum FIR source, and the solid angle of the outflow. Because the OH energy levels are radiatively pumped in the outflows, transitions at different wavelengths and energy levels, in combination with continuum component fits, yield crucial information about the radial location where the lines are formed. Besides the 79 μm doublet, the 119 (also ground state) and 65 μm (Elower = 290 K) doublets were also observed in all our sources, except for IRAS 17208. The excitation of the high-lying 65 μm doublet requires high far-IR radiation densities, indicating that the line is tracing the inner, highly excited, molecular region of the outflowing material, close to the nuclear source (González-Alfonso et al. 2008). Models that account for the 65 μm line absorption simultaneously reproduce the two ground-state lines as well, indicating that the three observed OH lines—in contrast to the atomic Na i D lines, which are extended on galactic scales in some cases—are mostly sensitive to the inner ∼0.5 kpc of the outflow. For IRAS 17208, only the 79 μm line is currently available, but we have assumed a similar compact outflow in OH. The OH columns are typically several 1016 cm−2. In Table 2, we list the fit results. Figure 1 shows the various components of the spectral fits on top of the observed spectra.
Table 2. (Preliminary) Model Fit Results (Errors are Discussed in the Text)
| Source | Component | Rcin | n(OH)din | fe | θf | 〈N(OH)〉g |
|---|---|---|---|---|---|---|
| (pc) | (10−4 cm−3) | (deg) | (1016 cm−2) | |||
| Mrk 231 | HVCa | 105 | 7 | 0.6 | 90 | 3 |
| LVCb | 115 | 5 | 1.0 | 61 | 5 | |
| IRAS 08572+3915 | HVC | 110 | 20 | 0.2 | 77 | 4 |
| LVC | 110 | 40 | 0.3 | 69 | 4 | |
| IRAS 13120−5453 | HVC | 210 | 9 | 0.3 | 43 | 3 |
| LVC | 210 | 1 | 0.5 | 90 | 2 | |
| IRAS 14378−3651 | HVC | 100 | 3 | 0.5 | 90 | 1 |
| LVC | 120 | 2 | 1.0 | 90 | 1 | |
| IRAS 17208−0014 | LVC | 110 | 6 | 1.0 | 57 | 2 |
| NGC 253 | LVC | 90 | 0.5 | 1.0 | 51 | 1 |
Notes. These values are used to compute the mass outflow rate in Table 1: 
, where tdyn ∼ R/v, g is a function of the opening angle θ, and χ(OH) = 5 × 10−6 is the OH abundance relative to H2. Outflow rates in Table 1 are the sum of the components (HVC+LVC). A more detailed description of a refined model will be given in a future paper.
aHigh velocity component.
bLow velocity component.
cInner shell radius.
dOH density at Rin.
eFilling factor.
fHalf opening angle (i.e., θ = 90° for a full 4π coverage).
gOH column density N(OH)× f.
Download table as: ASCIITypeset image
We derive mass-loss rates 
(see Table 1 and footnote of Table 2) from the OH density profile of the outflowing gas, the gas velocity in each component, and the adopted OH abundance. The resulting mass-loss rates are higher than 700 M☉ yr−1 in those sources with high vmax. The uncertainty in 
is dominated by the adopted OH abundance, the properties of the underlying FIR continuum source, and the outflow geometry. We adopt an overall uncertainty of four, derived from a study of the dependency of the fit results on the various model assumptions. The OH abundance relative to H2 is assumed to be 5 × 10−6, based on modeling of multi-transition OH observations of the Galactic giant molecular cloud Sgr B2 (Goicoechea & Cernicharo 2002). While OH abundances could reach values significantly lower in different environments (e.g., in some of the models of Sternberg & Dalgarno 1995 and Meijerink & Spaans 2005, or the observations by Watson et al. 1985), we expect a relatively high OH abundance close to the nuclear region (González-Alfonso et al. 2008) and note that 
increases with decreasing OH abundance. On the other hand, even a high OH abundance of 10−5 leads to mass-loss rates of several hundred M☉ yr−1. Further details on the modeling will be given in a forthcoming paper.
The starburst-dominated objects of our sample (NGC 253, IRAS 13120−5453, and IRAS 17208−0014) have values of 
similar to their SFRs. In contrast to this, the AGN-dominated and/or warm ULIRGs, as well as the heavily obscured ULIRG, have outflow rates that are factors ∼4–20 larger than their SFRs (see Table 1 and Figure 3, upper panel). The outflows we observe in the ULIRGs—if continued at the current rate—are able to completely expel the cold gas reservoirs from the centers of these objects within ∼106–108 years (Figure 3, lower panel). The gas depletion time is roughly inversely correlated with the AGN luminosity (and AGN fraction). This is consistent with the merger scenario where the highest outflow rates are a short-lived, late, AGN-dominated stage in the merger evolution. We note, however, that the lowest timescales we find (few ×106 years) appear to be inconsistent with star formation/galaxy evolution models. We may have underestimated the timescales: a higher gas mass by a factor two and a higher OH abundance by a factor two would result in four times longer timescales (which would still be within the error bars of Figure 3). A fundamental change in the basic model assumptions, like the geometry of the outflow, could also change the outflow rates. Our current models are consistent with wide opening angles for the outflowing gas, but this is hard to prove with PACS spectra alone (see below). For Mrk231, however, the high outflow rate is consistent (within a factor two) with the independent CO observations (PdBI mm-interferometry) by Feruglio et al. (2010). Very recently the high rates and the wide angle of the outflow have also been confirmed by optical spectroscopy through the analysis of broad wings of Hα and Na i D absorption (Rupke & Veilleux 2011). Finally, the outflow may not continue at such a high rate, since—as gas and dust are removed—the outflows may be quenched as the radiation pressure decreases. Therefore, our depletion timescales are time scales rather than absolute durations.
We do not have the spatial resolution with Herschel-PACS to study the spatial distribution of the outflowing material directly (although observations of additional OH transitions with higher energy level, e.g., at 84 and 163 μm, could help to further constrain the geometry of the models). True spatial information on the molecular component can only come from interferometric millimeter observations (e.g., Feruglio et al. 2010). While the observed spectra shown here clearly exhibit high velocity molecular outflows, independent of our modeling, the conclusions based on outflow rates and AGN fractions need further confirmation from increased statistics. Our ongoing observations will cover a much larger sample of both starburst- and AGN-dominated (U)LIRGs than shown here. Combined with millimeter interferometric follow-up observations we will be able to better constrain our OH outflow models and further investigate the surprisingly low depletion timescales for some of the objects. The final sample should then allow us to study potential trends in the outflow characteristics. Examples of such trends could be a different or tighter correlation of the outflow velocity with AGN luminosity than with SFR, or the extent to which AGN fractions and outflow velocities reflect different merger stages with evolving outflows creating lower and lower covering factors to the AGN.
We thank Dave Rupke for helpful discussions. Basic research in IR astronomy at NRL is funded by the US ONR; J.F. also acknowledges support from the NHSC. E.G.-A. is a Research Associate at the Harvard-Smithsonian Center for Astrophysics. A.S. thanks the DFG for support via German-Israeli Project Cooperation grant STE1869/1-1.GE625/15-1. 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); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); 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).
Footnotes
- *
Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.