Ubiquitous Molecular Outflows in z > 4 Massive, Dusty Galaxies. II. Momentum-driven Winds Powered by Star Formation in the Early Universe

Where necessary throughout this work, we assume a spherical "time-averaged thin shell" geometry widely used in the literature (see Rupke et al. 2005), in which a mass-conserving outflow with constant outflow rate expands following a density profile $n\propto {r}^{-2}$. In this geometry, the outflow rate and mass are related through

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}=4\pi {R}_{\mathrm{out}}^{2}\ \overline{\mu }\,{m}_{{\rm{H}}}\ {N}_{{\rm{H}}}\ {v}_{\mathrm{out}}/{R}_{\mathrm{out}}\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{M}_{\mathrm{out}}={\dot{M}}_{\mathrm{out}}{R}_{\mathrm{out}}/{v}_{\mathrm{out}},\end{eqnarray} \tag{ 1b }$$

with ${v}_{\mathrm{out}}$ the characteristic outflow velocity, ${R}_{\mathrm{out}}$ the outflow inner radius, $\overline{\mu }=1.4$ the mean mass per hydrogen atom (including the cosmological helium abundance) and m_H the mass of a hydrogen atom. These quantities are fundamentally linked to the column density of gas along the line-of-sight N_H responsible for generating the observed absorption profiles. If we drop the assumption of spherical symmetry, these quantities are also linearly proportional to the covering fraction ${f}_{\mathrm{cov}}$, the fraction of the full 4π sr covered by wind material as seen from the source. For our sample, in which the outflows are detected in absorption, the inferred energetics obviously depend strongly on the orientation of the outflow since no absorption can be detected for outflowing material that does not intersect the line of sight to the galaxy. Redshifted receding material is also difficult to constrain for our high-redshift sample. The ALMA bandwidth probes only a limited range of redshifted velocities, and it is possible for the galaxy itself to be optically thick to emission from the receding material even if the spectral coverage were extended to more redshifted velocities.

As has been discussed extensively in the literature, this assumed geometry leads to more conservative outflow energetics than other simple geometries (e.g., Veilleux et al. 2020, and references therein). In particular the outflow rates are a factor of 3 lower than if the outflow volume is filled with uniform density (which implies a decreasing outflow rate over time for constant flow velocity), and a factor ${\rm{\Delta }}r/{R}_{\mathrm{out}}$ lower than the "local" or "instantaneous" rate if the wind arises from a thin shell of width Δr. We take the characteristic velocity to be ${v}_{84}$, the velocity above which 84% of the absorption occurs. This is also fairly conservative: clearly the maximum velocity is not a "characteristic" outflow velocity, but ${v}_{84}$ is more robust to uncertainties in the systemic redshift than the median absorption velocity ${v}_{50}$, given the deep absorption at systemic velocities present in most of our sources (Paper I). Finally, we take ${R}_{\mathrm{out}}$ to be ${r}_{\mathrm{dust}}$, the effective radius of the dust emission at rest-frame ≈100 $\mathrm{\mu m}$. This radius has been directly measured for the low-redshift literature reference sources by Herschel/PACS (Lutz et al. 2016) and from our lensing reconstructions of the ALMA OH continuum data, and is motivated by our observation that the OH outflow absorption is frequently strongest in equivalent width not in the nuclear regions but toward the outskirts (Paper I).

The momentum and kinetic energy outflow rates are then given by

$$\begin{eqnarray}&&{\dot{p}}_{\mathrm{out}}={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}},\ {\dot{E}}_{\mathrm{out}}=\displaystyle \frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^{2},\end{eqnarray} \tag{ 2 }$$

where the expression for kinetic power assumes negligible contribution from turbulent (i.e., non-bulk) sources.

Although maps of the molecular outflows at ∼500 pc resolution are available for our SPT sample due to their gravitationally lensed nature (Paper I), we do not attempt to match, or otherwise account for, the structures seen in these maps when estimating outflow rates.

In order to determine whether and how well the outflow properties for our high-redshift sample can be estimated given the sole available OH 119 μm transitions, we compare extensively to low-redshift ULIRGs and obscured QSOs with rich OH data and radiative transfer models from Herschel/PACS. In particular, GA17 self-consistently model 12 nearby IR-luminous galaxies with detections of the OH transitions at 119, 84, 79, and 65 μm, all of which showed either P Cygni profiles or blueshifted line wings in the 119 μm doublet. The 84 and 65 doublets are highly excited lines with lower levels 120 and 300 K above the ground state that require an intense and warm IR radiation field to be detected, and the cross-ladder 79 μm doublet has an optical depth ≈40× lower than the 119 μm transition. We supplement this sample with one additional ULIRG with OH radiative transfer modeling (Tombesi et al. 2015; Veilleux et al. 2017), and an additional four sources with outflow rates based on the detection of high-velocity CO line wings that were also observed in OH 119 μm. These four sources have lower typical ${L}_{\mathrm{IR}}$ and lower outflow rates than the primary OH-based sample. We consider these final sources because Lutz et al. (2020) find reasonable agreement between CO-based and OH-based outflow properties, but we exclude them from our later empirical outflow rate estimation out of an abundance of caution. These sources can at some level be considered akin to a small cross-validation sample, although the dynamic range in outflow properties they span is small.

For all sources, we remeasured various properties of the OH 119 $\mathrm{\mu m}$ spectra in the same way as for our high-redshift sample (Paper I), including the broadening from the FWHM ≈ 300 $\mathrm{km}\,{{\rm{s}}}^{-1}$ PACS instrumental spectral resolution at these wavelengths. As for our sample, we use the fits to the spectra to measure the velocities above which 50% and 84% of the absorption takes place, ${v}_{50}$ and ${v}_{84}$, and the "maximum" outflow velocity ${v}_{\max }$ that we take to be the velocity above which 98% of the absorption takes place. We also measure the total equivalent widths of the absorption components as well as the equivalent widths integrated over various blueshifted velocity ranges; for example, ${\mathrm{EW}}_{v\lt -200}$ refers to the equivalent width integrated over velocities more blueshifted than −200 $\mathrm{km}\,{{\rm{s}}}^{-1}$. We note that all these quantities are nonparametric and therefore depend little on the exact methods used to fit the PACS spectra.

We first consider a simple analytic calculation of the outflow rates for our sources and the literature reference sample assuming the OH 119 $\mathrm{\mu m}$ absorption is optically thin. As already discussed, we expect this to be a very bad assumption, but this calculation does at least provide a hard lower bound on the true outflow rate and an opportunity to determine if some overall correction factor to the optically thin outflow rates could allow a more realistic estimate. While for many nearby galaxies other OH transitions with far lower line opacities can be observed (e.g., the 79 $\mathrm{\mu m}$ doublet, or lines of the less abundant ¹⁸OH isotopologue), we must instead attempt to find some other quantity that can provide an empirical correction.

Under the assumption that the absorption is optically thin, the minimum column density of OH molecules ${N}_{\mathrm{OH}}$ is given by

$$\begin{eqnarray}&&{N}_{\mathrm{OH}}=\displaystyle \frac{8\pi {Q}_{\mathrm{rot}}({T}_{\mathrm{ex}})}{{\lambda }^{3}{g}_{u}{A}_{{ul}}}\displaystyle \frac{\exp ({E}_{l}/{T}_{\mathrm{ex}})}{1-\exp \left(\tfrac{-h\nu }{{k}_{{\rm{B}}}{T}_{\mathrm{ex}}}\right)}\int \tau (v){dv},\end{eqnarray} \tag{ 3 }$$

where λ and ν are the wavelength and frequency of the transition, h and k_B are the Planck and Boltzmann constants, A_ul is the Einstein "A" coefficient of the transition, g_u the degeneracy of the upper energy level, E_l the lower energy level in temperature units, Q_rot the rotational partition function evaluated at excitation temperature T_ex, and $\int \tau {dv}$ the integrated optical depth of the absorption profile (e.g., Mangum & Shirley 2015). For the OH 119 $\mathrm{\mu m}$ doublet transitions, E_l = 0 K, A_ul = 0.138 s⁻¹ and g_u = 6 (Müller et al. 2001, 2005). Tabulated values of Q_rot are available from the NASA JPL spectroscopic database (Pickett et al. 1998). We assume an excitation temperature ${T}_{\mathrm{ex}}=100$ K as found in literature OH studies (e.g., GA17); the combination of T_ex-dependent terms in Equation (3) varies by about a factor of 3 for $50\lt {T}_{\mathrm{ex}}\lt 150$. It is then straightforward to calculate the total (H) column density assuming an OH abundance, ${N}_{{\rm{H}}}={N}_{\mathrm{OH}}/[\mathrm{OH}/{\rm{H}}]$. We adopt an OH abundance $[\mathrm{OH}/{\rm{H}}]=2.5\times {10}^{-6}$, as commonly assumed in the literature based on OH studies of the Milky Way star-forming region Sgr B2 (Goicoechea & Cernicharo 2002).

In order to isolate only outflowing material, it is common to integrate the optical depth over a limited range of velocities. Here we calculate the integrated optical depth over velocities more blueshifted than −200 $\mathrm{km}\,{{\rm{s}}}^{-1}$, a commonly adopted threshold. Although this is not fast enough to be guaranteed to trace only outflowing material, we expect it to largely trace the outflows even in our sample DSFGs that often show broad CO or [C ii] emission line profiles (Paper I). Under the assumption of optically thin absorption, the integrated optical depth over this velocity range is equal to the equivalent width over the same range. In practice, because the outflow rate itself is proportional to ${v}_{\mathrm{out}}$ (Equation (1)), we instead calculate ${\int }_{-\infty }^{-200}\tau (v)\,v\,{dv}$ from our spectral fitting procedure, consequently incorporating the ${v}_{\mathrm{out}}$ term of Equation (1) into the column density calculation directly. This allows us to include a first-order consideration of the shape of the absorption profiles while also removing the need to adopt a characteristic velocity in the outflow rate calculation. Finally, because we expect that the wind material does not fully cover the source, we adopt a covering fraction ${f}_{\mathrm{cov}}=0.3$, the average value determined for the low-redshift reference sample (GA17). We discuss covering fractions in detail in Paper I, where we estimate covering fractions ranging from a hard lower bound of ∼0.1 to upper limits of ∼0.7. These covering fractions are also not directly comparable: GA17 estimate ${f}_{\mathrm{cov}}$ using their multi-transition OH radiative transfer analysis, while our estimates are based on our lensing reconstructions with hard lower limits based on spectral analysis. While assuming a different value for ${f}_{\mathrm{cov}}$ would linearly rescale our optically thin outflow rate estimates, this has no impact on our subsequent results, as we explain further below.

We use Equation (1) to calculate the outflow rates for our own and the literature reference samples. For the reference sample, we use far-IR continuum sizes from PACS 100 $\mathrm{\mu m}$ imaging (Lutz et al. 2016), or assume the average size ≈1 kpc if no data are available. We derive minimum optically thin outflow rates spanning 8−370 ${M}_{\odot }$ yr⁻¹ for the literature sources and 5−120 ${M}_{\odot }$ yr⁻¹ for the SPT sample, and corresponding outflow masses $\mathrm{log}({M}_{\mathrm{out}}/{M}_{\odot })\approx 7\mbox{--}8.5$. We emphasize that these values are strong lower limits given the expected high OH line opacities.

Figure 1 (left) shows the ratio of the optically thin outflow rates to the published values from the OH radiative transfer models and CO line wings for the literature sample; upper limits in this plot correspond to those sources that were spatially unresolved by PACS and therefore have upper limits on ${R}_{\mathrm{out}}$. As expected, the optically thin assumption likely underestimates the true outflow rate by a large factor, ≈4–30× for most sources. The fact that the outflow rates can be so drastically underestimated is at some level a testament to the sensitivity of OH 119 $\mathrm{\mu m}$ to even minute amounts of outflowing material: with ALMA at z > 4, in principle it is possible to detect molecular outflow rates of just ∼10 ${M}_{\odot }$ yr⁻¹ in less than an hour of observing time. This is a consequence of both the relatively high OH abundance and especially the large value of the Einstein A_ul for the ground-state 119 $\mathrm{\mu m}$ transition, ∼10⁵ times larger than for [C ii] 158 $\mathrm{\mu m}$ or ∼10⁶ times larger than for low-order CO transitions. The drawback to this sensitivity, of course, is that the line opacities are high and it is not easy to determine by what exact factor the true outflow rate has been underestimated, as Figure 1 shows.

Zoom In Zoom Out Reset image size

There is no obvious trend between the ratio of optically thin to published outflow rate and the published outflow rate itself; evidently the OH 119 $\mathrm{\mu m}$ line opacity varies by a large amount in different galactic winds. We do, however, identify correlations between this ratio of outflow rate estimates and various measures of the outflow velocity. Sources with the fastest outflows are also the closest to being consistent with optically thin absorption. Figure 1 (right) shows this in terms of ${v}_{\max }$, but we obtain results consistent within the uncertainties from ${v}_{84}$ and ${v}_{50}$ as well; while ${v}_{\max }$ is somewhat more difficult to measure (Paper I), it shows the largest dynamic range among the outflow velocity metrics. Such a correlation makes intuitive sense: for a given column density of absorbing gas, if the total gas column extends over a larger range in velocity, the line opacity per unit velocity interval must necessarily be lower and therefore the line opacity averaged over the full absorption profile must also be lower. This leads to a less extreme "correction factor" needed for sources with very fast outflows.

We fit a simple log-linear function to the data in Figure 1 (right), determining uncertainties on the fit using the same Monte Carlo resampling method we use for all outflow rate techniques. We find a best-fit relation for the outflow rates "corrected" from the optically thin values of

$$\begin{eqnarray}&&\mathrm{log}({\dot{M}}_{\mathrm{out}}^{\mathrm{thin}\,\mathrm{corr}.}/{\dot{M}}_{\mathrm{out}}^{\mathrm{thin}})=m({v}_{\max }+1000)+b,\end{eqnarray} \tag{ 4 }$$

with $m=-{6.4}_{-1.7}^{+1.8}\times {10}^{-4}$ ($\mathrm{km}\,{{\rm{s}}}^{-1}$)⁻¹ and $b={0.91}_{-0.06}^{+0.07}$, and ${v}_{\max }$ in $\mathrm{km}\,{{\rm{s}}}^{-1}$. The analysis indicates an intrinsic dispersion of ∼0.15 dex around the best-fit relation in addition to the statistical uncertainties. We use this relation and the measured values of ${v}_{\max }$ to estimate the true outflow rates empirically corrected from the optically thin assumption.

We note again that our prior assumption of ${f}_{\mathrm{cov}}=0.3$ in our calculation of ${\dot{M}}_{\mathrm{out}}^{\mathrm{thin}}$ has no impact on our "corrected" outflow rates, as a different assumed value propagates directly into b in the equation above. There is also no evidence for a correlation between ${f}_{\mathrm{cov}}$ and other galaxy properties in the low-redshift training sample that could influence our outflow rates given the differences between the samples (GA17). While we do find a tentative correlation of ${f}_{\mathrm{cov}}$ with ${L}_{\mathrm{IR}}$ (Paper I), we expect those covering fractions to be upper limits on the true values and stress again that the methods used between low- and high redshift are not directly comparable. Our assumed ${f}_{\mathrm{cov}}=0.3$ lies well within the lower and upper limits we expect for the true values, so we do not expect this to add substantial additional uncertainty beyond the present estimates. Both the optically thin and the corrected outflow rates are given in Table 1.

As we expected, the optically thin outflow rates almost certainly severely underestimate the true outflow rates. While we derived a method to correct these values to more realistic outflow rates, the correction factors remain highly uncertain and in any case the general methodology deserves to be cross-checked by other methods. We now consider two simple empirical methods to provide alternative estimates of the outflow rates before moving to a more complex empirical method.

In S18 we made a simple estimate of the true outflow rates using a subset of the present literature reference sources for which outflows had also been detected in CO emission. In that work we took the OH 119 $\mathrm{\mu m}$ equivalent widths for the low-z literature sources integrated over the velocity ranges where high-velocity wings of CO emission had been detected (GA17), under the philosophy that both traced molecular outflows and that the outflows should appear over the same velocity range in both tracers. Here we follow a similar vein, now including an expanded reference sample. Instead of individually choosing velocity ranges over which to measure the OH equivalent widths, here we simply fit a linear relationship between the published literature ${\dot{M}}_{\mathrm{out}}$ values and ${\mathrm{EW}}_{v\lt -100}$ from our remeasured 119 $\mathrm{\mu m}$ spectral fits, which gave equivalent widths similar to those we used in S18.

Figure 2 (left) shows the results of this analysis. We do not force this fit to have a zero intercept. Although this allows the unphysical scenario of positive outflow rates in the absence of any absorption or even negative outflow rates for low equivalent widths, allowing this freedom in the model yields a better characterization of the uncertainty at low ${\dot{M}}_{\mathrm{out}}$. We find a best-fit expression for the outflow rate

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}^{{\rm{S}}18}=m{\mathrm{EW}}_{v\lt -100}+b,\end{eqnarray} \tag{ 5 }$$

with $m={2.7}_{-0.5}^{+0.4}$ ${M}_{\odot }$ yr⁻¹/($\mathrm{km}\,{{\rm{s}}}^{-1}$) and b = 55 ± 50 ${M}_{\odot }$ yr⁻¹, the outflow rate in ${M}_{\odot }$ yr⁻¹ and the equivalent width in $\mathrm{km}\,{{\rm{s}}}^{-1}$. We find an intrinsic dispersion of ±150 ${M}_{\odot }$ yr⁻¹ around this relation in addition to the statistical uncertainties that at least applies in the low-${\mathrm{EW}}_{v\lt -100}$, low-${\dot{M}}_{\mathrm{out}}$ regime. At higher ${\dot{M}}_{\mathrm{out}}$ there are too few sources to quantify any additional scatter beyond the statistical uncertainties; we assume a constant 150 ${M}_{\odot }$ yr⁻¹ scatter for all values of ${\mathrm{EW}}_{v\lt -100}$.

Zoom In Zoom Out Reset image size

The OH equivalent width is not expected to be the sole controlling parameter that predicts outflow rates, of course. Herrera-Camus et al. (2020) (abbreviated HC20) explored an alternative simple parameterization, fitting the outflow rates to the product ${\mathrm{EW}}_{v\lt -200}\times \sqrt{{L}_{\mathrm{FIR}}}$. This was motivated by an expectation that the outflow rate should depend on both the column density of outflowing gas (related to ${\mathrm{EW}}_{v\lt -200}$) and the size of the source (proportional to $\sqrt{{L}_{\mathrm{FIR}}}$ through a Stefan–Boltzmann type relation), as in Equation (1).

We repeat a similar analysis as HC20, with a couple small modifications. First, we use $\sqrt{{L}_{\mathrm{IR}}}$ instead of $\sqrt{{L}_{\mathrm{FIR}}}$, which is more readily available for all literature reference sources. Second, we do not fit a line forcing the y-intercept to be zero as done in HC20. Again, this allows us to better understand the uncertainties at low ${\dot{M}}_{\mathrm{out}}$. Our best-fit relation for the outflow rate in this way is

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}^{\mathrm{HC}20}=m({\mathrm{EW}}_{v\lt -200}\sqrt{{L}_{\mathrm{IR}}/{10}^{12}{L}_{\odot }})+b,\end{eqnarray} \tag{ 6 }$$

again with the outflow rate in ${M}_{\odot }$ yr⁻¹ and ${\mathrm{EW}}_{v\lt -200}$ in $\mathrm{km}\,{{\rm{s}}}^{-1}$, $m={1.40}_{-0.25}^{+0.21}$ ${M}_{\odot }$ yr⁻¹/($\mathrm{km}\,{{\rm{s}}}^{-1}$) and $b={180}_{-30}^{+40}$ ${M}_{\odot }$ yr⁻¹. We find an essentially identical intrinsic scatter around this relation as before, ≈150 ${M}_{\odot }$ yr⁻¹, where this is assumed to be constant due to the lack of sources with very high outflow rates.

Aside from a more physically justified parameterization, this method also has a slightly higher dynamic range in the abscissa than the S18-style fit. Between the two methods, we have some preference for the HC20 parameterization. Outflow rates derived from both methods are given in Table 1.

Finally, we consider a more complex empirical model to derive outflow rates. While the analyses in the previous subsection relied on specific linear correlations between observables and published outflow rates, there is no particular reason to choose those specific observables over others apart from some physical intuition about the likely important parameters. The fact that we find significant additional intrinsic scatter beyond the inferred uncertainties in the previous fits is a clue that a more complex model connecting the observables and the outflow rates is warranted. Indeed, both our reference and high-redshift samples have many more known properties than we have yet utilized, both from the OH spectra themselves as well as ancillary measurements from other data. Here we perform one final analysis that attempts to discern the most predictive relationship between all available measurements and the outflow rates, at the expense of linking the resulting relationship to any particular physical meaning.

To explore the complex relationship between outflow rates and all available measurements, we use a "partial least squares" (PLS) technique (Wold 1966). PLS is both a regression and dimensionality reduction technique, and can be thought of as somewhat of a hybrid between standard multivariate linear regression and principal component analysis (PCA) or singular value decomposition. PLS is well-suited to cases such as ours where the number of objects in the reference sample is relatively few but the number of measured quantities for each sample object is large, with many of the measured quantities correlated with each other. In our case, for example, while ${v}_{84}$ and ${v}_{\max }$ encode slightly different information about the shape of the OH absorption profile, they are still strongly correlated: fast outflows are fast regardless of the metric used. While PCA techniques are capable of describing the variance in the observables, not all principal components need be predictive of some other quantity (${\dot{M}}_{\mathrm{out}}$ in our case). PLS addresses this by maximizing the covariance between the space of observables and the space of desired predicted quantities. PLS also performs better than some other techniques (e.g., random forest estimators) when some measured observables of the target sample (in this case the SPT objects) lie outside the dynamic range of the training sample—that is, when extrapolation is required for one or more observables. In the end this has little influence on our application because the most predictive observables (see below) are well-sampled by the reference objects and extrapolation is not generally required.

For our purposes, we use PLS to predict the outflow rates ${\dot{M}}_{\mathrm{out}}$ from a variety of (sometimes strongly correlated) observed properties: several metrics of the OH velocity profiles and equivalent widths integrated over various velocity ranges, as well as ancillary galaxy properties such as ${L}_{\mathrm{IR}}$, ${r}_{\mathrm{dust}}$, the AGN contribution to the bolometric luminosity ${f}_{\mathrm{AGN}}$, and the effective dust temperature ${T}_{\mathrm{dust}}$. We experimented extensively with various numbers and combinations of observables and found consistent results for the predicted outflow rates of the SPT sources in almost all cases. Generally regardless of the observables used, a maximum of four or five PLS components minimized the mean squared error in the predicted outflow rates of the reference sample (that is, the dimensionality of the problem could be reduced from the number of observables used to four or five, due to covariances between the observables employed). PLS also allows us to understand which observables are most responsible for driving predictions for the outflow rates. Of those we explored, the outflow velocity and equivalent width were the most predictive of the measured outflow rates, while ${f}_{\mathrm{AGN}}$ and r_dust generally had little predictive power, possibly due to the relatively small dynamic range in these quantities in the training and target samples.

Figure 3 shows the comparison between predicted and published outflow rates for the combination of parameters that includes ${v}_{50}$, ${v}_{84}$, ${v}_{\max }$, ${\mathrm{EW}}_{v\lt -100}$, ${\mathrm{EW}}_{v\lt -200}$, ${\mathrm{EW}}_{\mathrm{total}}$, ${L}_{\mathrm{IR}}$, and ${f}_{\mathrm{AGN}}$. Interestingly, unlike the previous methods, there is no longer any detectable intrinsic scatter between the predicted and published outflow rates; the gray shaded region in Figure 3 illustrates the approximate upper limit on the scatter we can set with the data available, ≈100 ${M}_{\odot }$ yr⁻¹. Although it remains rather unsatisfying to necessarily discard all physical interpretation of the resulting predictions, clearly PLS is capable of translating the complex measurement space into the desired output outflow rates. Predicted outflow rates from this method are provided in Table 1.

Zoom In Zoom Out Reset image size

We now have four different estimates for the molecular outflow rates applied to the high-redshift SPT objects—one corrected from the optically thin assumption, two simple empirical estimators, and one more complex empirical estimate. Figure 4 compares these estimates for both the low-redshift reference sample and as applied to our z > 4 objects. Note that for the literature sources this figure only compares the outflow rates predicted from each method to each other; comparisons with the "true" published values can be found in the preceding figures.

Zoom In Zoom Out Reset image size

Essentially by definition this figure shows good agreement between the various methods for the reference sample, since this sample was used to derive the conversions between observables and outflow rates in the first place. We also find generally good agreement between the estimators for the SPT objects, in particular between the multivariate PLS analysis and the simpler approach of HC20. Evidently these methods make use of the most salient predictive measurements from the OH spectra and ancillary galaxy properties.

Figure 5 compares the outflow rates for the SPT sources in more detail. This figure shows the outflow rates derived from each method for each source. This figure again demonstrates the generally good agreement between methods, although the PLS and HC20-like methods tend to yield slightly higher values than the other methods. We also show a joint distribution of the outflow rates created by equally combining the Monte Carlo trials from each method. While this should not be considered a true joint probability distribution of the outflow rates—the methods to derive input distributions are hardly independent, for example—it both highlights the level of agreement or disagreement between methods and summarizes the constraints we place on the outflow rates. These joint outflow rates are listed in Table 1, referred to as ${\dot{M}}_{\mathrm{out}}^{\mathrm{joint}}$. We use these joint estimates and associated uncertainties throughout the remainder of the text as our "best" estimates of the outflow rates, subsequently dropping the "joint" superscript from the notation for simplicity.

Zoom In Zoom Out Reset image size

The joint distributions from each method suggest that we are able to estimate the outflow rates for our sources at about the factor of 2 level. The OH-based reference sources have quoted uncertainties at the ∼50% level; the higher level of uncertainty for our sources reflects the lack of additional OH data for our sample that propagates into the scatter seen in the four individual methods and thus into our final joint estimates.

We emphasize that for both samples these uncertainties are likely underestimated due to systematics in many of the assumptions, from the OH abundance to the assumed geometry and outflow history. For our high-redshift objects, while our estimates are empirically based, the methods we have described presume that low-redshift IR-luminous galaxies are sufficiently similar to our targets as to not render these calculations meaningless. While the observed characteristics of our sample are contained within the parameter space probed by the reference sample (Figures 1–3 and Paper I), it is certainly possible that some other unmeasured quantity has a strong influence on the outflow rates that is not accounted for by our methods. Thus, while we propagate the uncertainties on the joint outflow rates in the remainder of the text, it is important to remember that these are probably more uncertain by some difficult to quantify amount.

In low-redshift samples, correlations between outflow velocities and host galaxy properties such as SFRs or AGN luminosities have been used to shed light on the physical mechanism(s) responsible for launching the outflows. There are good theoretical reasons to believe that the energy and momentum imparted to the gas from star formation and/or AGN activity should play a role in driving galactic winds, and should then manifest in the properties of the outflows launched. Nevertheless, observations of neutral and low-ionizations species show at best weak correlations between outflow velocities and SFR from the local universe to z ∼ 1 (e.g., Weiner et al. 2009; Rubin et al. 2014; Chisholm et al. 2015; Roberts-Borsani & Saintonge 2019), with any trend mostly due to the weak outflows seen in very low-SFR galaxies (e.g., Heckman & Borthakur 2016).

While this could plausibly be because the neutral outflows are less strongly coupled to the driving source, similarly weak correlations have also been seen for the molecular phase traced by OH in nearby ULIRGs and QSOs (e.g., Veilleux et al. 2013). While still relatively weak, the strongest correlations between outflow velocities and galaxy properties are found with ${L}_{\mathrm{AGN}}$ and ${f}_{\mathrm{AGN}}$, suggesting a connection between the AGN and wind launching at least in these extreme nearby systems.

Figure 6 shows three outflow velocity metrics, ${v}_{50}$, ${v}_{84}$, and ${v}_{\max }$, as a function of ${L}_{\mathrm{IR}}$, the IR surface density ${{\rm{\Sigma }}}_{\mathrm{IR}}$, ${f}_{\mathrm{AGN}}$, and ${L}_{\mathrm{AGN}}$, where we now compare our high-redshift objects to the combined sample of low-redshift ULIRGs and QSOs and nearby AGN-dominated systems (described in detail in Paper I). We distinguish between objects with outflows (filled symbols) and those without (empty), as determined by the original authors. We also note that the subset of low-redshift ULIRGs selected for OH radiative transfer modeling by GA17 is skewed toward sources with the fastest outflows, presumably because these were a more viable sample for multi-transition modeling. We return to this point several more times because it propagates into many of the differences we see with the low-z ULIRGs in other outflow properties as well.

Zoom In Zoom Out Reset image size

The left column of Figure 6 first shows outflow velocities as a function of ${L}_{\mathrm{IR}}$. In agreement with Veilleux et al. (2013), we see no evidence of a correlation in the expanded sample of nearby galaxies. Interestingly, however, we do see hints of a trend within the z > 4 SPT DSFGs when considered alone, with the most luminous sources also driving the fastest outflows. Whether this is a genuine difference between the outflows driven in low- and high-redshift galaxies remains to be seen; a larger sample of high-redshift objects that spans a wider range in ${L}_{\mathrm{IR}}$ and other properties will be required to understand these tentative differences further.

Instead of ${L}_{\mathrm{IR}}$ alone one might instead expect the outflow velocity to depend more strongly on the IR surface density ${{\rm{\Sigma }}}_{\mathrm{IR}}$ (or similarly the SFR surface density), for example, in cases in which radiation pressure on dust grains drives the outflows (e.g., Thompson et al. 2015). The second column of Figure 6 shows these quantities for the low- and high-redshift OH molecular outflow samples, where we used far-IR sizes measured from Herschel/PACS imaging (Lutz et al. 2016, 2018) for the low-redshift samples and the sizes from our lensing reconstructions for the SPT sample (Paper I). We find no convincing evidence of correlation between these quantities even for the SPT sample considered alone.

Of the parameters investigated by Veilleux et al. (2013) for low-redshift ULIRGs and QSOs, the strongest correlations with outflow velocities were found with ${f}_{\mathrm{AGN}}$ and ${L}_{\mathrm{AGN}}$,¹⁹ which those authors argued could be due to obscuration effects whereby the fastest-moving material was more easily visible in the AGN that had already cleared the nuclear regions or were oriented face-on. The subsequent addition of far less luminous AGN-dominated systems by Stone et al. (2016) agreed with this picture although the number of sources with definite outflows was small. The third column of Figure 6 shows outflow velocities against ${f}_{\mathrm{AGN}}$. We also fit a simple linear function to the low-redshift sources with molecular outflows, finding a marginally significant correlation with ${v}_{50}$ that becomes weaker with ${v}_{84}$ and ${v}_{\max }$; the scatter is clearly large. The limits on ${f}_{\mathrm{AGN}}$ for the SPT sources based on rest-frame mid-IR photometry do not clearly result in these objects being outliers, and they certainly would not be outliers even if we have underestimated ${f}_{\mathrm{AGN}}$ by a substantial amount (Section 3).

Finally, the right column of Figure 6 shows outflow velocities as a function of ${L}_{\mathrm{AGN}}$, which Stone et al. (2016) found to be strongly correlated in low-redshift sources in agreement with Veilleux et al. (2013). We also see some relationship between these quantities, namely, sources with low AGN luminosities rarely drive fast outflows. However, we note that while the Stone et al. (2016) sample certainly extends the dynamic range in ${L}_{\mathrm{AGN}}$ probed, this now conflates samples selected in very different ways, with many other possible confounding variables (mass, for example). Regardless, we again find that the limits we can place on ${L}_{\mathrm{AGN}}$ for the SPT sample again do not make them obvious outliers.

In summary, among the SPT DSFGs alone, the total ${L}_{\mathrm{IR}}$ appears to be most strongly correlated with outflow velocity, although a larger sample size will be required to investigate whether this is genuine. While we recover correlations previously noted in low-redshift work with our larger combined literature sample, the OH outflow velocities appear to be at best weak indicators of the driving source of molecular outflows, with substantial scatter. While we currently have no evidence for AGN activity in the SPT DSFGs and only weak limits on ${f}_{\mathrm{AGN}}$, our objects are not obvious outliers in plots of outflow velocity and AGN properties given the substantial scatter seen among the low-redshift objects, and we thus cannot rule out that AGN are responsible for driving the molecular outflows we have observed.

A number of recent works have explored scaling relations between molecular outflow properties and host galaxy properties, compiling samples of now dozens of objects (e.g., Cicone et al. 2014; GA17; Fluetsch et al. 2019; Lutz et al. 2020). While these studies focused exclusively on low-redshift galaxies, we now include our measurements for the first sample of molecular outflows in the early universe. Our primary comparison samples are the OH-based outflow measurements in nearby ULIRGs from GA17, as before, supplemented with the CO-based sample of Lutz et al. (2020), which extends to lower-luminosity systems. All samples assume the same outflow geometry as discussed in Section 2.1. We note that Lutz et al. (2020) found that OH-based outflow rates tended to be ≈0.5 dex higher than CO-based rates in their comparison of galaxies observed in both tracers (while the total outflow masses ${M}_{\mathrm{out}}$ were very similar). Because the CO-based sample spans a different range of parameter space than the other samples, we also comment on how our inferences in this section would change if the CO outflow rates were increased by 0.5 dex. We also detail changes to our interpretation that would result from doubling our present upper limits on ${f}_{\mathrm{AGN}}$ to try to account for the effects of any heavily obscured AGN that may not be detectable even in the rest-frame mid-IR. It is important to note that none of these samples at any redshift are complete or unbiased; the galaxies typically targeted for molecular outflow observations are highly biased toward luminous star-forming systems and/or quasars.

Figure 7 shows the molecular outflow rate ${\dot{M}}_{\mathrm{out}}$ and mass loading factor ${\eta }_{\mathrm{out}}\equiv {\dot{M}}_{\mathrm{out}}/$ SFR as a function of SFR. We find uniformly sub-unity mass loading factors for the high-redshift DSFGs, although the uncertainties of course remain significant. We would still find loading factors ≲1 even if we have underestimated our limits on ${f}_{\mathrm{AGN}}$ by a factor of 2 (which would consequently lower the SFR). This is a perhaps surprising result—these galaxies are among the most luminous, highest-SFR objects known, yet drive relatively weaker outflows than many less luminous nearby galaxies (though again, none of these samples is complete or unbiased). In particular, despite SFRs a few times higher than the low-redshift OH sample of GA17, the outflow rates we derive do not increase accordingly and the loading factors are consequently lower. At least some of this difference is likely due to the selection for fast outflows in the low-redshift work, but because the outflow rate depends only linearly on the velocity this is insufficient to explain the full difference.

Zoom In Zoom Out Reset image size

Our sample does appear, however, to follow the slightly sublinear relationship seen in some low-redshift studies (e.g., Fluetsch et al. 2019), extended now to an order of magnitude higher SFR and the high-redshift universe. Fitting a power law to the combined samples, we find a best-fit relationship $\mathrm{log}({\dot{M}}_{\mathrm{out}})=(0.72\pm 0.05)\mathrm{log}(\mathrm{SFR})+(0.5\pm 0.1)$, with ${\dot{M}}_{\mathrm{out}}$ and SFR in ${M}_{\odot }$ yr⁻¹ and an additional intrinsic scatter on the outflow rates of ≈0.25 dex. Thus, we find a transition from ${\eta }_{\mathrm{out}}\gt 1$ to sub-unity values near an SFR of ∼ 100 ${M}_{\odot }$ yr⁻¹. Because of the distribution of the CO-based Lutz et al. (2020) sample in SFR, increasing the CO-based outflow rates by 0.5 dex would further flatten the power-law slope to ≈0.5 but increase the transition SFR at ${\eta }_{\mathrm{out}}=1$ to 500 ${M}_{\odot }$ yr⁻¹. On the other hand, if we lowered the SFRs of our sample by doubling ${f}_{\mathrm{AGN}}$ to estimate the effect of a possible highly obscured AGN, the power-law slope would marginally increase to ≈0.8.

Interestingly, we find no evidence that the dominance of an AGN plays a secondary role in setting the outflow rate, once the overall trend with SFR (or ${L}_{\mathrm{IR}}$) has been accounted for. The points in Figure 7 are color-coded by ${f}_{\mathrm{AGN}}$, and it is clear that the remaining scatter in the ${\dot{M}}_{\mathrm{out}}$–SFR relationship is uncorrelated with ${f}_{\mathrm{AGN}}$. This remains a point of some contention in the recent literature: Cicone et al. (2014) and Fluetsch et al. (2019) claim a correlation between ${\eta }_{\mathrm{out}}$ and ${f}_{\mathrm{AGN}}$ for ${f}_{\mathrm{AGN}}\gtrsim 0.7$, while Lutz et al. (2020) found no such correlation, even though both studies use a largely overlapping set of literature outflow detections (since we use the Lutz et al. 2020 CO-based literature compilation it is no surprise that we also find no correlation given the fairly small increase in dynamic range in SFR afforded by our sample). While a thorough analysis of this low-redshift discrepancy is beyond the scope of this paper, part of the difference may lie in how the outflow rates were calculated from the CO line wings, as Lutz et al. (2020) included only the wings of the broad CO component, while Fluetsch et al. (2019) included the entire broad component (i.e., including emission at systemic velocities that may not actually be part of the outflow).

Figure 8 shows ${\dot{M}}_{\mathrm{out}}$ and ${\eta }_{\mathrm{out}}$ as a function of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ instead. The SPT sample lies well within the scatter but shows typically lower values of ${\eta }_{\mathrm{out}}$ at a given ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ compared to the low-redshift samples. As noted above, this is at least partially explained by the overall sublinear trend between ${\dot{M}}_{\mathrm{out}}$ and SFR. Similar to our investigation of outflow velocities with ${{\rm{\Sigma }}}_{\mathrm{IR}}$ above, we again find no significant correlation between these properties, a conclusion that would not change if we increased the CO-based outflow rates by 0.5 dex or adopted 2× higher limits on ${f}_{\mathrm{AGN}}$ for our sample. As noted by Lutz et al. (2020), however, the combined literature sample (and our own, clearly) is not complete in ${{\rm{\Sigma }}}_{\mathrm{SFR}}$. Severe selection effects stemming from the diverse selection criteria in individual studies comprising the combined literature sample may exist, particularly at low ${{\rm{\Sigma }}}_{\mathrm{SFR}}$.

Zoom In Zoom Out Reset image size

Finally, Figure 9 shows ${\dot{M}}_{\mathrm{out}}$ and ${\eta }_{\mathrm{out}}$ as a function of ${L}_{\mathrm{AGN}}$. As seen in previous works, the low-redshift combined sample shows a clear relationship between ${\dot{M}}_{\mathrm{out}}$ and ${L}_{\mathrm{AGN}}$, with substantial scatter that increases at low ${L}_{\mathrm{AGN}}$. The distribution of low-redshift objects in this parameter space has been discussed extensively in the literature (e.g., Lutz et al. 2020). For our sample, we find that the limits on ${L}_{\mathrm{AGN}}$ from the available rest-frame mid-IR data do not result in the high-redshift objects being clear outliers in Figure 9, and they would not be outliers in the event that our limits on ${f}_{\mathrm{AGN}}$ were underestimated by a factor of 2. As before in Section 4.1 but from a different perspective, the outflows we have detected do not require an AGN based on mass loading factor scaling relations, but AGN activity also cannot be ruled out in our objects.

Zoom In Zoom Out Reset image size

We now turn to estimates of the total molecular gas mass contained within the outflows. As described in Section 2.1, we use Equation (1) to calculate the outflow masses ${M}_{\mathrm{out}}$ for our sample. For the low-redshift samples, we use the original published masses for both the OH-based and CO-based outflows. Recalculating the masses for the low-redshift samples using our assumed geometry results in <10% differences in the median compared to the published values. Additionally, Lutz et al. (2020) found only a 0.06 dex offset and 0.3 dex dispersion between the masses for the low-redshift sources with outflows observed in both OH and CO. We make use of the total molecular gas masses for the low-redshift samples assembled by the original studies, all of which are based on low-J transitions of CO (CO(3–2) or lower). The conversion factor between CO luminosity and ${M}_{{{\rm{H}}}_{2}}$ is known to vary based on various galaxy properties (e.g., Bolatto et al. 2013); we accept the values used by the original studies. For the SPT sample, we assume ${\alpha }_{\mathrm{CO}}=0.8$ ${M}_{\odot }$ (K $\mathrm{km}\,{{\rm{s}}}^{-1}$ pc²)⁻¹, which we have previously found to be appropriate for the IR-luminous galaxies in our sample (e.g., Spilker et al. 2015; Aravena et al. 2016).

Figure 10 shows the molecular outflow masses as a function of ${L}_{\mathrm{IR}}$, as well as the fraction of the total galaxy molecular gas mass contained in the outflows. Not unexpectedly, ${M}_{\mathrm{out}}$ is clearly correlated with ${L}_{\mathrm{IR}}$, as has been previously noted many times in the literature. The z > 4 SPT DSFG sample has molecular outflow masses in the range $\mathrm{log}{M}_{\mathrm{out}}/{M}_{\odot }=8.6\mbox{--}9.1$, unsurprisingly on the high end of the local samples. The masses of the two most intrinsically luminous sources in our sample, SPT2132-58 and SPT2311-54, are perhaps somewhat low in comparison to the extrapolation of the low-redshift samples, but are well within the observed scatter. From the lower panel of Figure 10, meanwhile, we find that the molecular outflows in our sample contain 1%–10% of the total molecular gas masses of the galaxies. These values are well within the range typically seen in low-redshift galaxies. Further, we find no discernible trend between ${M}_{\mathrm{out}}$/${M}_{{{\rm{H}}}_{2}}$ and ${L}_{\mathrm{IR}}$ despite the increase in dynamic range in ${L}_{\mathrm{IR}}$ afforded by our sample.

Zoom In Zoom Out Reset image size

Figure 11 shows these same quantities as a function of ${L}_{\mathrm{AGN}}$. Together, Figures 10 and 11 are effectively the corresponding versions of Figures 7 and 9 for ${M}_{\mathrm{out}}$ instead of ${\dot{M}}_{\mathrm{out}}$. As with the outflow rates previously, the current limits on ${L}_{\mathrm{AGN}}$ for the SPT DSFGs do not make them obvious outliers in Figure 11. There are no indications from either of these figures that the dominance of the AGN plays any role either in determining ${M}_{\mathrm{out}}$ in general or in defining the scatter in ${M}_{\mathrm{out}}$ at a given ${L}_{\mathrm{IR}}$ or ${L}_{\mathrm{AGN}}$, as evidenced by the lack of secondary trends in these figures with ${f}_{\mathrm{AGN}}$.

Zoom In Zoom Out Reset image size

Meanwhile, Figure 12 compares the molecular gas depletion timescale due to outflows with that due to star formation. Here these depletion timescales are defined as the time it would take for the entire molecular gas reservoir of the galaxies to be removed by outflows or consumed by star formation, assuming the outflow rate or SFR remain constant. That is, ${t}_{\mathrm{dep},\mathrm{out}}\equiv {M}_{{{\rm{H}}}_{2}}/{\dot{M}}_{\mathrm{out}}$ and ${t}_{\mathrm{dep},\mathrm{SF}}\equiv {M}_{{{\rm{H}}}_{2}}$/SFR. These depletion times are only approximate estimates of the important timescales in the evolution of galaxies, given that both outflows and star formation operate simultaneously (which would give shorter depletion times), we ignore molecular gas destruction due to, e.g., photo-heating or shocks (which would also shorten the depletion times), gas accretion and/or cooling into the molecular phase are neglected (which would give longer depletion times), and ${\dot{M}}_{\mathrm{out}}$ and SFR are not in fact constant over time (which could push the depletion times either higher or lower depending on the time variability in ${\dot{M}}_{\mathrm{out}}$ and SFR). Note that changing estimates of ${M}_{{{\rm{H}}}_{2}}$ for any object in Figure 12 moves objects diagonally parallel to the one-to-one line, since ${M}_{{{\rm{H}}}_{2}}$ is incorporated in both axes.

Zoom In Zoom Out Reset image size

For the z > 4 sample, we find for all sources that ${t}_{\mathrm{dep},\mathrm{out}}\gtrsim {t}_{\mathrm{dep},\mathrm{SF}}$, a straightforward consequence of the sub-unity wind mass loading factors we determine (Section 4.2). This conclusion would also hold if we artificially decrease the SFRs of our sample by doubling ${f}_{\mathrm{AGN}}$ as a crude approximation of the effects of a heavily obscured AGN, though the depletion times would be about equal in that case. As before with ${\eta }_{\mathrm{out}}$, this places our sample with a distinct minority of the low-redshift samples. We stress again that all of these samples are biased toward galaxies and quasars that do host powerful outflows, and these results may not hold for objects with less extreme outflows that would be difficult to detect. Unsurprisingly given their very high SFRs and outflow rates, both depletion times are very short, ∼10–100 Myr, and the fact that the two timescales are comparable points to the important role that outflows must play in regulating star formation in these galaxies.

Galactic winds are often classified as either "energy driven" if radiative losses in the outflowing gas are negligible or "momentum driven" if they are not (regardless of the ultimate source(s) of the energy driving the wind). In the former case, the outflow is thought to be launched by the adiabatic expansion of a bubble of hot gas (e.g., Chevalier & Clegg 1985; Silk & Rees 1998) that either lofts cold gas entrained in the expanding hot wind or (re-)forms cold gas from the swept-up shocked material at larger radii where the gas can radiatively cool (e.g., Faucher-Giguère & Quataert 2012; Costa et al. 2014; Richings & Faucher-Giguère 2018). Similar to the energy-conserving Sedov–Taylor phase of supernova expansion, the resulting momentum in the outflowing gas can be "boosted" well above the radiative momentum flux driving the wind. In the momentum-driven case, in which radiative cooling is significant, momentum transferred to the gas from ram pressure or radiation pressure on dust grains results in gentler acceleration that may allow cold gas to reach large radii and high velocities before it is destroyed (e.g., Murray et al. 2005, 2011; Thompson et al. 2016; Brennan et al. 2018). Real winds can of course be intermediate between these cases. For both momentum- and energy-driven winds and for winds driven by an AGN or star formation, theoretical models provide estimates of the coupling efficiency between the input momentum and energy and the outflowing gas.

We calculate estimates of the outflow momentum and energy for our SPT sample as described in Section 2.1. For both the low-redshift OH- and CO-based samples we use the original published values of the outflow momentum and energy instead of those derived from our own assumptions in Section 2.1. A comparison between the published values and our estimates indicates that we may be overestimating the outflow momentum and energy by ∼30% and 70%, respectively, while we find no systematic difference in the mass outflow rates. While still within the substantial uncertainties, this probably means that the values of ${v}_{84}$ we use in Equation (2) are higher than the outflow "characteristic" velocity, and some lower velocity between ${v}_{50}$ and ${v}_{84}$ would provide more accurate outflow energetics. We continue with our use of ${v}_{84}$; our conclusions here would be further strengthened if the outflow momentum and energy were systematically lower.

In Figure 13 we show the fraction of the estimated outflow momentum rate compared to the total radiative momentum rate (${\dot{p}}_{\mathrm{rad}}=L/c$) as a function of the estimated luminosities due to AGN and star formation for the sample of molecular outflows assembled from the literature and our high-redshift sample. In momentum-driven winds the AGN may provide a momentum rate up to ∼2 ${L}_{\mathrm{AGN}}$/c, treating both radiation pressure on dust grains and the AGN inner winds as ${L}_{\mathrm{AGN}}$/c. Meanwhile, a continuous starburst can generate a maximum of ∼3.5 ${L}_{\mathrm{SF}}$/c (Veilleux et al. 2005; Heckman et al. 2015) through a combination of radiation pressure and the pressure of hot wind material driven by supernova ejecta. As seen in Figure 13, molecular outflows in low-redshift galaxies frequently show large momentum boosts ∼2–30 above the radiative momentum provided by the AGN and/or star formation, often taken as evidence that an energy-driven wind phase is required to achieve such large boosts (though see also Thompson et al. 2015, who argue that radiation pressure on dust grains can also achieve large momentum boosts in conditions possibly realized in very dusty and gas-rich galaxies).

Zoom In Zoom Out Reset image size

We find much more modest momentum ratios in our sample of high-redshift DSFGs, with maximum momentum boosts of ∼2 compared to the luminosity due to star formation, and all sources consistent with no momentum boost above the radiative momentum injection at all. This momentum boost is well within the range achievable by radiation pressure on dust in cases where the effective IR optical depth is of order unity (Murray et al. 2005; Thompson et al. 2015). Further, the momentum injection due to star formation alone is fully consistent with the observed outflow momentum fluxes; no additional radiative momentum from AGN is required. Indeed, it is not clear if the AGN alone could provide sufficient momentum to explain the observed outflows given the current limits on ${L}_{\mathrm{AGN}}$; at least some substantial contribution from star formation would be required if AGN are relevant to the outflow energetics. This result would not change if we redistributed the total luminosity arising from the AGN and star formation by doubling ${f}_{\mathrm{AGN}}$ compared to our current limits, although in that case the momentum flux from the AGN would also be sufficient to drive the outflows we observe. All sources would still show momentum boosts ≲3.5 ${L}_{\mathrm{SF}}$/c, would still be consistent with momentum-driving due to star formation, and would not show momentum boosts as large as those seen in the local ULIRGs. For our sources to be >1σ inconsistent with the rough maximum ∼3.5 ${L}_{\mathrm{SF}}$/c would require ${f}_{\mathrm{AGN}}\gt 0.8\mbox{--}0.99$ depending on the source, far above our current limits from the rest-frame mid-IR. Sources with such high ${f}_{\mathrm{AGN}}$ typically show OH solely in emission in nearby objects (see discussion in Stone et al. 2016), while none of our sources show OH in emission. This could be taken as evidence that no source in our sample has ${f}_{\mathrm{AGN}}\gtrsim 0.9$.

Figure 14 shows a similar plot for the outflowing kinetic power, following our outflow calculations in Section 2.1. Hot energy-driven winds from the AGN are thought to be capable of supplying up to about ∼5% of the AGN power to the outflows, of which some fraction ∼1/2 can plausibly be converted into bulk kinetic energy in the wind (e.g., Faucher-Giguère & Quataert 2012; King & Pounds 2015). The mechanical luminosity generated by supernovae during a starburst, meanwhile, may reach ∼2% of the total starburst luminosity, with perhaps ∼1/4 of this luminosity converted into kinetic motion in the interstellar matter (ISM; e.g., Veilleux et al. 2005; Harrison et al. 2014). The outflow energetics in many low-redshift molecular winds exceed the expected coupling efficiency to the starburst luminosity, while the AGN energetics are in better agreement (Figure 14). This has been taken as evidence that the AGN must be primarily responsible for driving the low-redshift molecular outflows, and in combination with the momentum rates in Figure 13, that these winds must be at least partially energy driven.

Zoom In Zoom Out Reset image size

In contrast to these low-redshift results, we find that the outflow kinetic energy rates in our z > 4 DSFGs are uniformly below the threshold coupling efficiency for supernova-driven winds, and would still be consistent with this coupling efficiency if we adopt limits on ${f}_{\mathrm{AGN}}$ twice as high (or more) as current data indicate. As with the momentum rates, AGN are not required in order to explain the observed outflow energetics. Moreover, the AGN in our sample could be an order of magnitude less luminous than the current limits without the outflow kinetic power approaching the theoretical maximum of ∼few percent of the AGN luminosity.

Taken together, we conclude from Figures 13 and 14 that (1) the high-redshift molecular outflows we have observed are fully consistent with expectations for momentum-driven winds, with no need for partially or fully energy-conserving phases, and (2) the observed outflow energetics can be fully explained by the momentum and energy provided by star formation alone in these galaxies, with no need for additional driving by AGN. We emphasize that we do not conclude that AGN are not responsible for driving the observed outflows, merely that AGN are not required to explain the energetics. We note again that these conclusions are further strengthened if we adopt a somewhat lower characteristic velocity in Equation (2) as it appears may be appropriate by comparison to the OH-based outflow energetics (GA17). Similarly, our conclusions are also not changed in the event that our present limits on ${f}_{\mathrm{AGN}}$ are underestimated by a factor of 2 (or more) due to AGN so heavily obscured they are not detectable in the mid-IR. In that case either the AGN or star formation could be the ultimate driving source, but the outflow energetics would still not require AGN momentum or energy injection or energy-conserving phases.

Both these results are surprising and counter to the conclusions typically reached in low-redshift studies. Conventional wisdom dictates that AGN are necessary to regulate galaxy growth in massive galaxies, in part due to scaling relations such as that between black hole mass and galaxy or bulge mass. Yet in our high-redshift rapidly star-forming galaxies, we find no need for an AGN in order to explain the molecular outflow energetics we have measured. While our sample objects are more luminous than almost all of the low-redshift sources, we have no reason to expect that the outflow energetics should not also increase concomitantly with luminosity. Additionally, the energy-conserving wind mode is generally thought to have the highest coupling efficiency with the ISM, capable of sweeping up a large fraction of the gas in the ISM (e.g., Zubovas & King 2012). Our results, however, show that such high-efficiency energy-driven winds are not necessary to explain the observed outflow momenta and kinetic energy rates in our sample.

The differences between the z > 4 DSFGs and nearby ULIRGs—both with outflow properties from OH spectroscopy—are particularly striking given their general similarities as highly dust-obscured and IR-luminous galaxies. GA17 found that additional energy injection from an AGN is required to explain the energetics of most of the low-redshift molecular outflows in ULIRGs (Figure 14) and that at least partially energy-conserving wind phases are likely necessary to explain the large momentum boosts (Figure 13). Neither of these appears to be true for the high-redshift DSFG outflows. We also note that a similar conclusion appears to be true for the only other z > 4 object with detected OH absorption, a z = 6.1 quasar where we expect the total luminosity to be dominated by the AGN (HC20) in contrast to our own sample with only upper limits on ${L}_{\mathrm{AGN}}$. Although an estimate of the AGN and starburst luminosities separately is not available for this source and the OH detection was low significance, applying the same outflow property calculations to this source as our sample would also place it in the general vicinity of our sample objects as long as ${f}_{\mathrm{AGN}}\gtrsim 0.1$, a condition easily met for luminous quasars.

It is tempting to ascribe at least some of the differences we see compared to the nearby ULIRGs to the overall difference in luminosities between the low- and high-redshift sources. Due to observational limitations the high-redshift objects are typically several times more luminous than the low-redshift ULIRGs. Increasing the luminosity of the ULIRGs would move them down and right in Figures 13 and 14, in the direction that would be required to unify the low- and high-redshift objects. However, this would imply that the outflow momentum rates and kinetic power have essentially reached their maximum in low-redshift ULIRGs and no longer continue to increase in more luminous systems as observed locally (GA17). It is also possible that the physics of outflows is qualitatively different between the low- and high-redshift samples. Multiple simulation efforts have found that star formation-driven outflows become inefficient in massive galaxies at z ≲ 1, so it could be that the low-redshift samples are predisposed toward AGN-driven winds by virtue of the fact that they have outflows detected at all (e.g., Muratov et al. 2015; Hayward & Hopkins 2017). It is clear that a larger sample at high redshift that spans a wider range in parameter space than current observations will be required to understand the dependencies of outflow energetics on galaxy properties.

There is also a probable selection effect that appears to be at play in the low-redshift sample. As shown in Figure 6, while our sample overlaps with the low-redshift samples by most metrics, the subset of low-z sources selected for detailed OH radiative transfer modeling by GA17 have preferentially higher outflow velocities than the low-redshift sample overall, likely because these sources presented a more tractable sample for their modeling. This may weight the low-redshift sample toward AGN-driven (fast) outflows. Additionally, a bias toward fast winds can sharply skew the outflow energetics because the outflow velocity enters at least linearly in the outflow momentum rates and at least quadratically in the kinetic power (the outflow rates themselves are also proportional to ${v}_{\mathrm{out}}$). We thus expect that the local ULIRGs with slower outflows would show substantially lower momentum and kinetic energy outflow rates that extend down to the values we find for the SPT sources. For the majority of nearby ULIRGs, then, we expect that the outflow energetics would also be consistent with momentum-driven winds that do not require additional energy injection from the AGN.

The molecular outflows we have observed could plausibly affect the host galaxies over cosmological timescales, especially if large fractions of the cold gas in the outflows travel at sufficiently high velocity to escape the galaxy or even the dark matter halo virial radius. In the latter case, now unbound, the gas may never again be available for star formation. In the former, the gas becomes part of the CGM and could recycle back into the galaxy unless continued energy injection or shock heating prevents the gas from cooling and condensing (see Tumlinson et al. 2017, for a recent review).

We make a simple estimate of the fraction of the outflowing molecular gas that will escape the host galaxies by assuming the outflowing mass as a function of velocity is directly proportional to the equivalent width as a function of velocity, excluding the absorption components centered on systemic velocities (Paper I). To estimate the galaxy escape velocity for each source, we assume a spherical isothermal mass distribution truncated at a maximum radius ${r}_{\max }/r=10$, following Arribas et al. (2014). Because the detection of absorption requires the presence of continuum emission, we take r to be the circularized effective size of the dust emission from our lensing reconstructions, ${r}_{\mathrm{dust}}$. We estimate galaxy masses from total molecular gas masses based on CO(2–1) observations, assuming a typical gas fraction for DSFGs at these redshifts (Aravena et al. 2016). These masses are in reasonable agreement with simple dynamical mass estimates using the available [C ii] or CO line widths and the lens model sizes (e.g., Spilker et al. 2015). We find escape velocities for our sources ranging from ∼400–1000 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (median ∼700 $\mathrm{km}\,{{\rm{s}}}^{-1}$), which agree reasonably well with other simple estimates scaling from the CO or [C ii] line widths or assuming pointlike mass distributions within ${r}_{\mathrm{dust}}$. Given the uncertainties in mass and shape of the gravitational potential, we estimate typical uncertainties on the galaxy escape velocities of ≈40%.

In this calculation, we ignore any additional deceleration of the outflow caused by sweeping up additional material. We have also implicitly assumed that the outflowing material is located at a typical distance from the galaxy center equal to the dust continuum emitting size, which seems reasonable based on our lensing reconstructions of the outflow material (Paper I), but we cannot rule out that much of this material is located deeper within the gravitational potential wells of the host galaxies. Both these effects would lower the fraction of outflowing gas that escapes the galaxies. On the other hand, we also assume that the outflowing gas is no longer being accelerated, which may not be the case if the winds are driven by the outward radiation pressure on dust grains, especially in the event of high far-IR optical depths and/or cosmic-ray pressure. This would result in higher outflow escape fractions than our estimates.

Figure 15 shows the cumulative outflow mass for each object in our sample as a function of the outflow velocity, normalized to the estimated escape velocity. We find typical galaxy escape fractions ∼20% with large variation within the sample. The three objects with estimated escape fractions >25% are SPT2132-58 and SPT2311-54, which have the fastest outflows of our sample, and SPT2319-55, which has an atypically low mass given its outflow velocity (or an atypically fast outflow given its mass). Only ≲10% of the outflowing gas is traveling at 1.5 times the escape velocity or faster, and essentially none is traveling at twice the escape velocity.

Zoom In Zoom Out Reset image size

Figure 15 also shows these escape fractions as a function of ${L}_{\mathrm{IR}}$, now including the local galaxy samples with stellar velocity dispersion measurements available to estimate the escape velocities. The uncertainties on the escape fractions include only those due to the uncertain escape velocities, but do not include the (unknown) contribution from any variations in the equivalent width to outflow mass proportionality (or the CO–H₂ conversion factor for the CO-based masses). We find a very wide range in estimated galaxy escape fractions from ≈0% up to 60%, reflective of the large range in outflow velocities and a rather limited dynamic range in galaxy mass. The escape fraction shows no obvious correlation with ${L}_{\mathrm{IR}}$ or other observables, in agreement with our conclusions in Section 4.1.

The galaxy outflow escape fractions in Figure 15 are substantially higher than those found by Fluetsch et al. (2019) even for the same objects. As discussed in Section 4.2, this is because those authors count the full broad CO component as belonging to the molecular outflows, thereby including a substantial amount of CO flux at systemic velocities that need not actually be outflowing. This additional flux (and therefore mass) artificially lowers the galaxy escape fractions well below the values we obtain following the more conservative definition of Lutz et al. (2020), who only considered the flux in the broad line wings in the outflow definition (excluding the core emission at systemic velocities). This more conservative definition results in total outflow masses a factor of ≈5 lower on average for the CO-based objects in Figure 15 and consequently higher escape fractions compared to those found by Fluetsch et al. (2019).

While the uncertainties are large, the nearby ULIRGs in Figure 15 tend to show somewhat larger escape fractions than our own sample of high-redshift objects. These sources have a mean and median escape fraction ≈40%, about double that for our own sample. As previously discussed, this is most likely due to the fact that the sources with available OH-based radiative transfer models are preferentially also those with the fastest outflows (Figure 6). Given the lack of correlation between outflow velocity and stellar velocity dispersion or stellar mass over the limited dynamic range probed by these samples (e.g., Veilleux et al. 2013), this results in outflow escape fractions skewed toward larger values. As in Section 4.4, we expect that a more complete sample of local ULIRGs would show substantially more overlap with the lower escape fractions we find for the high-redshift DSFGs.

The bulk of the molecular gas in the outflows is destined to remain within the galaxies, where it can become available for future star formation through a galactic fountain flow. At least in the cold molecular phase, most of the gas will not be permanently expelled and therefore these outflows cannot really be responsible for the very low gas fractions that are one of the hallmarks of quenched galaxies at lower redshifts (e.g., Young et al. 2011; Davis et al. 2016; Spilker et al. 2018a; Bezanson et al. 2019). Moreover, without continuous injection of thermal energy or turbulence over the long term, the CGM gas will develop a cooling flow resulting in significant gas accretion (e.g., Su et al. 2020).

Finally, we consider the impact of the outflowing molecular gas that probably will escape in the context of the CGM surrounding these high-redshift DSFGs.²⁰ The top panel of Figure 16 shows the mass of the molecular outflows traveling at speeds greater than the galaxy escape velocity in each source. We assume all of this material enters the CGM, ignoring the loss of any material that escapes the larger dark matter halos (we expect this to be an exceedingly small fraction given the outflow velocity distributions in Figure 15). For most of the SPT DSFGs, we expect ∼few 10⁸ ${M}_{\odot }$ of the outflowing molecular gas to become incorporated into the CGM of the host halos.

Zoom In Zoom Out Reset image size

The typical CGM properties of DSFGs are virtually unknown. Based on a sample of three $z\sim 2$ DSFGs with background quasar sightline absorption spectra, Fu et al. (2016) speculated that the CGM of DSFGs may be less massive and/or that DSFGs inhabit somewhat less massive dark matter halos than co-eval quasars. However, given the much better statistics available for quasars at these redshifts, Figure 16 shows the typical range of total cool ($T\lesssim {10}^{4}$ K) CGM gas mass within the virial radius of $2\lt z\lt 3$ quasar host galaxies thought to reside in $\mathrm{log}{M}_{h}/{M}_{\odot }\sim 12\mbox{--}13$ mass halos (Prochaska et al. 2014; Lau et al. 2016).²¹ Given the possible differences between the CGM of DSFGs and quasars and an expectation that the CGM grows in mass from z ≥ 4 to z = 2.5, we expect this range to be an approximate upper bound on the total cool CGM mass surrounding the higher-redshift DSFGs in our sample.

The bottom panel of Figure 16 shows the total mass in metals being ejected into the CGM, under the simplifying assumption that the molecular outflows have approximately solar metallicity. If, as the outflow energetics suggest (Section 4.4), processes related to star formation are responsible for driving the molecular outflows, we may expect the outflowing gas to be enriched significantly beyond solar, moving the points upwards in the lower panel of Figure 16. In comparison, the metallicity of the cool CGM gas surrounding $2\lt z\lt 3$ quasars is subsolar, Z ∼ 0.1–0.3 Z_⊙, likely because it is a mixture of metal-enriched outflow gas, less metal-rich material stripped or ejected from infalling satellites, and metal-poor material accreting from the cosmic web (e.g., Muratov et al. 2017; Hafen et al. 2019). The range of total CGM metal mass for the same z ∼ 2–3 quasar samples is also shown in Figure 16.

Taken together, the two panels of Figure 16 give an intriguing (if admittedly speculative) picture of the relationship between molecular outflows and the CGM surrounding these galaxies. If the DSFGs in our sample will evolve to become like the quasars observed at slightly lower redshift, the current molecular outflow episodes will contribute only a small fraction ≲1%–10% of the total cool CGM mass. Evidently the total CGM mass must be assembled from some combination of outflowing gas in warmer phases than we have observed, many repeated outflow events through the lifetime of the galaxies, and accretion of additional gas into the CGM from the cosmic web or infalling satellites. While observations of the multiphase components of outflows are rare even in the nearby universe, it appears that in general the molecular phase contains a significant if not dominant portion of the total outflow mass (Fluetsch et al. 2019), so additional mechanisms beyond accounting for the unobserved warmer phases are likely required. On the other hand, the current outflow episodes can contribute some substantial fraction ∼10% or more of the total metals present in the CGM at lower redshift. This fraction would rise further if the outflowing molecular gas is enriched beyond solar metallicity. X-ray observations of the hot plasma in nearby winds typically find α/Fe elemental abundance patterns (i.e., including oxygen, of relevance to our OH observations) enhanced to several times the solar value (e.g., Nardini et al. 2013; Veilleux et al. 2014; Liu et al. 2019), though the composition of the molecular gas in outflows is unknown, even at low redshift.

The outflow metallicities of these highly obscured galaxies are conceivably observable with future observations of far-IR fine structure lines (e.g., Nagao et al. 2011; Pereira-Santaella et al. 2017). Indeed, the [C ii] 158 $\mathrm{\mu m}$ line has recently been detected on 10–30 kpc spatial scales surrounding co-eval lower mass galaxies through stacking and, in a few cases, direct individual detections (Fujimoto et al. 2019, 2020; Ginolfi et al. 2020). These studies conclude that metal-enriched outflows are the most likely source of the extended [C ii] emission, as generally expected from simulations (e.g., Muratov et al. 2015; Hayward & Hopkins 2017; Pizzati et al. 2020).