PHIBSS: Searching for Molecular Gas Outflows in Star-forming Galaxies at z = 0.5–2.6

The Plateau de Bure High-z Blue Sequence Survey (PHIBSS) is a molecular gas survey of 174 typical MS galaxies spanning 0.5 ≤ z ≤ 2.6. Observations were carried out in two large programs, PHIBSS 1 (2009–2011; PIs: L. Tacconi & F. Combes) and PHIBSS 2 (2013–2017; PIs: F. Combes, S. García-Burillo, R. Neri & L. Tacconi), using the Plateau de Bure Interferometer (PdBI)/NOEMA. Observations cover a range of beam sizes from 065–94 × 065–525, with spectral resolutions from 7 km s⁻¹ to 88 km s⁻¹. For all sources, the galaxy properties were derived using ancillary data, mainly Hubble Space Telescope (HST)/WFC3 and Advanced Camera for Surveys (ACS), using spectral energy distribution modeling to get the stellar masses, and Sérsic profile fitting to estimate the effective radii. The SFRs were estimated using UV and IR luminosities, or extinction-corrected Hα luminosities for the highest redshift sources. In all cases, the photometry comes from apertures encompassing the flux from the entirety of the galaxy. Similarly, we extract the CO flux from apertures designed to cover the whole galaxy for the majority of sources. Thus, all our measurements should reflect the global galaxy properties. Detailed descriptions of the observations, sample, and galaxy properties’ estimations can be found in L. J. Tacconi et al. (2013, 2018) and J. Freundlich et al. (2019).

The observations cover the CO(3–2) line for 99 galaxies at redshift 1.00 ≤ z ≤ 2.55, the CO(2–1) line for 70 lower redshift (0.50 ≤ z ≤ 1.53) galaxies, and the CO(6–5) line for the remaining five galaxies (2.01 ≤ z ≤ 2.33). For each galaxy, we compute the integrated signal-to-noise ratio (SNR) of the narrow emission line and discard detections below SNR = 1.5. Although this threshold is lower than what is typically considered significant, the stacking process aims to uncover the signal hidden within the noise, so we retain these detections in the stack. This process brings the sample down to 154 galaxies, with median 〈z〉 = 1.04, , 〈R_e〉 = 4.4 kpc, and〈ΔMS〉 = 0.14 dex (see Table 1, and the corresponding property distributions are plotted in Appendix A). The individual galaxy properties are reported in Table A1 in Appendix A. The distribution is plotted in Figure 1 as a function of their distance from the MS (computed following the parameterization from K. E. Whitaker et al. 2014), stellar mass, and redshift. In addition, we verified that the SNR cut does not change our results or improve the SNR of the stacked spectra by repeating the cut with higher threshold values (SNR > 3 and >5) and performing the same analysis as described in Section 3, finding that using higher SNR thresholds yields the same results.

Table 1. Median Properties of the Subsamples Stacked in Section 4.4, and Their Respective Upper Limits on Mass Outflow Rates and Mass Loading Factor

Sample	〈z〉	〈log(M⋆/M⊙)〉	〈Re〉	〈νobs〉	〈SFR〉	vout		ηUL
			(kpc)	(GHz)	(M⊙ yr−1)	(km s−1)	(M⊙ yr−1)
All	1.05	10.8	4.35	136.9	50.1	450–2300†	135–1567	2.7–31.3
AGN	1.09	11.0	3.25	140.6	112.9	1360†	1619	14.3
M* < 1010.7M⊙	1.01	10.5	4.30	134.2	31.6	380†	78	2.5
M* > 1010.7M⊙	1.01	11.0	4.45	140.7	63.1	420–1360†	237–1383	3.8–21.9
ΔMS > 0.2	1.13	10.7	3.60	137.7	79.4	390–1500†	373–2814	4.7–35.4
ΔMS < 0.2	1.02	10.9	5.00	137.1	31.6	450–1520†	85–525	2.7–16.6
z > 1.7	2.21	10.8	3.50	158.5	108.1	460–1000†	345–1105	2.2–6.9
z < 1.7	0.76	10.8	5.25	144.0	31.6	460–1990†	106–956	3.4–30.2
Subset Sample	1.01	11.0	5.30	140.7	50.1	1857 ± 905	1528 ± 793	30 ± 16

Note. Outflow velocities from the ionized gas population trends reported in N. M. Förster Schreiber et al. (2019) are marked with a † symbol. In cases where outflow velocities are reported for both SF- and AGN-driven outflows, we use two velocity values corresponding to the average outflow velocity for these subsamples.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Since AGNs can be the main driver of strong outflows in the stellar mass range covered by PHIBSS (e.g., N. M. Förster Schreiber et al. 2019; R. L. Davies et al. 2024), we identify the AGNs among the sample using publicly available multiwavelength catalogs. We classify the galaxies as AGN based on X-ray, mid-IR, optical, and radio flux diagnostics.

First, we cross-correlate the PHIBSS catalogs with the Chandra X-ray source catalog.¹⁷ We then apply a luminosity cut based on the galaxies’ expected 2–10 keV X-ray luminosity from their SFR (M. Symeonidis et al. 2014), such that galaxies with an associated X-ray detection 1 dex higher than the expected luminosity purely from SF are flagged as AGNs. From the full sample, 136 galaxies are in Chandra coverage, 24 of which have an associated X-ray detection. Of these 24, 17 fulfill the aforementioned criteria. Next, we crossmatch the PHIBSS survey with the MOSDEF public line-emission catalogs (M. Kriek et al. 2015; N. A. Reddy et al. 2015; A. E. Shapley et al. 2015), finding 15 galaxies in both samples. We investigate the [O iii]λ5007/Hβ versus [N ii]λ6584/Hα BPT diagram (J. A. Baldwin et al. 1981), using the upper limit star-forming abundance sequence reported by L. J. Kewley et al. (2013) at the average redshift of the 15 galaxies. We find six galaxies lying above the limit, which are thus classified as AGN. Following this, we use the Spitzer/IRAC photometry available for 116 of the galaxies to identify AGN based on the color criteria presented by J. L. Donley et al. (2012). Eight of the galaxies are classified as AGN based on their mid-IR colors. Lastly, we crossmatch the sample with public source catalogs from JVLA radio surveys. Following I. Delvecchio et al. (2017), we classify as AGN all those sources with an excess in L_1.4GHz compared to the value expected for SFGs according to their far-IR-radio scaling relation. Seven of the PHIBSS galaxies satisfy this criterion. In total, 24 galaxies fulfill one or more of the above criteria. The galaxies are reported as AGN in Table A1. We note that due to the sample selection aiming for typical SFGs, none of the sources in our sample are luminous AGNs or have emission dominated by AGN contribution.

We use the following method to extract the spectra of individual galaxies and maximize the SNR. First, we define a beam-sized aperture at the center of the data cube, sum the spectra of each pixel within the aperture, and perform a Gaussian fit to the line profile. Using the fit results, we collapse the initial data cubes along the velocity channels, summing all the velocity bins within the FWHM of the line and centering on the best-fit velocity. We then run SEXTRACTOR (E. Bertin & S. Arnouts 1996) on the collapsed image to determine the area of detection of the galaxy and create a spatial mask of this detection. Since, in most cases, the detected area by SEXTRACTOR is smaller than the beam size, we use the beam size as the aperture for the mask, which we center on the SEXTRACTOR coordinates. When the detected area is larger than the beam, we use this as the aperture for the mask. In most cases, the source lies at the center of the field of view and is not resolved above the beam size. However, for some sources, the iterative process helps improve the detection substantially by either identifying the correct source position or assessing the detection area properly. An example of such a case is shown in Appendix B, where the source is slightly offset from the center of the frame, and Figure B1 shows the initial extracted spectrum compared to the spectrum extracted once correcting for the source position, clearly showcasing the difference in detection.

Finally, we sum all nonmasked pixels in the 3D cube and fit the resulting spectrum with a single Gaussian. We iterate over this process three times to obtain the best possible mask (i.e., the mask for which the integrated spectrum returns the line with the highest SNR). Iterating over this process optimizes the detection area and the estimates of the central velocity and FWHM, which are both crucial to our stacking analysis for precise alignment and to correctly rebin the spectra (see Section 3.2). The final apertures have median major axis kpc and minor axis kpc. The final FWHMs are reported in Table A1, and their distribution is shown in Appendix A, Figure A1.

After the final iteration, we extract the spectra from the masked cube. We also scale each spectrum’s flux by normalizing it to the best-fit amplitude of the Gaussian fit to the emission line, such that the brightest galaxies do not dominate the stacked spectrum, as we are interested in the population average.

The spectral stacking of emission lines with different widths, due to the different ranges of projected velocities for different galaxies, can lead to the creation of a broadlike component in the stacked line profile, which can be misinterpreted as an outflow signature in the stacked profile (see, e.g., Section 5.2; F. Stanley et al. 2019 and A. Concas et al. 2022 for detailed discussions).

To mitigate this shortcoming of the stacking method, we rebin the spectral channels of each galaxy spectrum according to their estimated FWHM such that each line would have an FWHM covering the same number of spectral channels. Upon inspection of the distribution of FWHMs in our sample (see Appendix A, Figure A1), we choose a width of seven channels as the common line width for rebinning, to favor downsampling but to limit the bin size change to a moderate amount (with a maximum factor of 6).¹⁸ After rebinning, the average channel width is 37 ± 25 km s⁻¹. This method impacts the accuracy of the velocity information one can retrieve from the spectra but preserves the line shapes such that any underlying broad component revealed by the stacking could not be attributed to numerical effects. We investigate in Section 5 whether rebinning impacts our ability to retrieve real underlying outflows.

After rebinning, we stack the normalized spectra of the individual galaxies using the LINESTACKER package presented in J.-B. Jolly et al. (2020). Given a central velocity for each spectrum (here we use the best-fit velocity retrieved in the spectrum extraction process; see Section 3.1), the code aligns the spectra and computes the mean of the flux in each spectral channel, where each spectrum is given the same weight in the mean. We stack both the full sample and the following subsamples (median properties are listed in Table 1): above and below , above and below ΔMS = 0.2, and above and below z = 1.7. These choices are motivated by the trends in ionized gas outflow incidence of N. M. Förster Schreiber et al. (2019), corresponding to the stellar mass (ΔMS) above which AGN (SF)-driven outflows become more frequent. We also stack a subsample of galaxies with specific SFR, sSFR > 0.1 Gyr⁻¹, using the value of local sSFR for which an increase in outflows is observed (B. Reichardt Chu et al. 2025). Additionally, we separately stack the subsample of 24 AGN-hosting galaxies selected in Section 2.2. Finally, since outflow detection is also dependent on galaxy orientation, we create and stack subsamples divided by inclination for the 117 galaxies for which we have the inclination: below 30°, between 30° and 60°, and above 60°. The stacking results for the full sample, the AGN sample, and the subsamples are presented and discussed in Section 4.2.

To estimate the noise level in the stacked spectra, we extract a “blank” spectrum from an empty sky region using the same aperture as used to extract the corresponding galaxy spectrum. After stacking them following the same procedure described above, we take the resulting standard deviation as the noise estimate for each stack. Using this method instead of the noise directly in the spectra allows us to ensure we do not include any “hidden” signal in the noise estimation.

Finally, we repeat the stacking procedure on the same sample and subsamples but weigh each spectrum by its noise, such that noisier spectra weigh less in the stack. There are no differences in the results from both stacking prescriptions, and for the rest of the analysis, we use the results from the nonweighted stack.

As a first step of the analysis, we examine the individual spectra with integrated SNR > 7 to search for the presence of outflow signatures in these deeper data sets, as this SNR value indicates robust emission line detection. Out of the full sample of 154 galaxies, 20 have emission lines with SNR above this value (SNR =7.3–20.8, black circles in Figure 1), four of which are identified as AGNs (all galaxies have coverage by either Chandra, MOSDEF, Spitzer, or Very Large Array or VLA; see Section 2.2). This subsample of high SNR spectra has median 〈ΔMS〉 = 0.28, 〈z〉 = 1.16, and (see Table A1 and Figure 1), thus preferentially probing higher stellar masses (85% have ), and slightly higher z and ΔMS.

Three of these galaxies exhibit notable asymmetry in their CO line profile. For two of them, the line profile is explained by the fact that the galaxy is well detected over more than one beam and that the mask applied in Section 3.1 does not cover the full galaxy. This can happen during the spectrum extraction process, either if the fitting fails to encompass the full line, or if the aperture misses some flux. In such cases, when creating a collapsed cube to estimate the spatial extent of the detection, we do not sum all of the channels of the line, which results in a truncated detection area, and thus the final extracted spectra will also be missing part of the flux. The standardized iterative process described in Section 3.1 is designed to mitigate this effect as best as possible, but failed in those two cases of galaxies detected well above the beam size. In this case, the line from the initial mask has an asymmetric shape and becomes the double-peaked profile expected from the galactic rotation after adjusting the mask (see Appendix C, Figure C1).

For the remaining spectrum, EGS13003805, a z = 1.23 galaxy associated with an AGN, it is difficult to ascertain whether the spectrum shows signs of molecular gas outflows, as the line profile displays signs of asymmetry similar to what is described above (see Appendix C, Figure C2). The issue does not appear to come from the masking. Available HST/ACS I- and V-band imaging reveals EGS13003805 as a spiral galaxy (L. J. Tacconi et al. 2013). In addition, the galaxy is detected over an area larger than the beam size. With this information, and comparing the Atacama Large Millimeter/submillimeter Array detection size and the HST size (which are comparable), we favor the interpretation that the double peak in the integrated spectrum is due to galaxy rotation, rather than being an outflow signature, although a more in-depth analysis of the kinematics is required to fully rule out the outflow scenario.

Thus, we observe that prominent outflow signatures are not prevalent at high SNR, even for a subsample with galaxy properties for which we would expect higher outflow incidences and strength (high mass, high redshift, high SFR, AGN hosting). We also investigate later the presence of extended emission beyond the initial masking process (see Section 4.3). Moreover, we do not expect the masking issue presented above to affect our results more broadly as, in most cases, the detected area is smaller than the beam size, so the majority of our masks encompass more than the detection.

To fully exploit the PHIBSS data set, we now use the stacked spectra to search for potential outflow signatures that cannot be identified on an individual galaxy basis.

The results from the stacking of the spectra are shown in Figure 2 for the full sample and in Figures 3 and 4 for subsamples defined in Section 3.3. The stacked lines are fit both with a single and double-Gaussian component, where the fit has the following constraints: (1) one component cannot have both the larger amplitude and the larger width, (2) the amplitude of both components should be positive, (3) the center of the second component should be within 10 channels of the center of the first component, to avoid fitting noise peaks on the edges of the spectrum. In addition, for each channel, we define as the error on the flux the error derived in Section 3.3 from the stacked “blank” spectra, scaled by the number of sources stacked in that particular channel ( where N_obj is the number of sources stacked). This number varies from channel to channel, as all spectra do not have the same number of spectral channels (even prior to rebinning). When aligning them on the emission line and averaging channel-per-channel, the edges, in particular, will be the average of fewer channels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all cases, the single-Gaussian fit residuals, normalized by the noise in the spectra, are displayed below the stacked spectrum in Figures 2, 3, and 4. Since some of the individual integrated galaxy spectra display the double-peaked profile expected from inclined rotating disks, this feature sometimes propagates in the stacked spectra and can be seen in the residuals. However, no clear outflow signature is visible.

For completeness, we perform a two-Gaussian fit to the spectra and observe that it does not improve the goodness-of-fit. Indeed, the reduced χ² of the single-Gaussian fit is systematically lower than that of the double-Gaussian fit. We also perform an analysis of the Bayesian information criterion (BIC, G. Schwarz 1978)¹⁹ of our single and double-Gaussian fits. To do so, we define ΔBIC = BIC_sing − BIC_doub as the difference between the BIC of the single-Gaussian fit and the double-Gaussian fit, respectively, where ΔBIC > 0 favors the double-Gaussian model to fit the data, and ΔBIC > 10 is strong evidence that a double-Gaussian fit is needed to fit the data (A. R. Liddle 2007). For all stacked spectra, ΔBIC < 0, in line with the result that there are no outflow signatures in the stacks. In addition, the flux of the very low amplitude broad component is systematically below the noise level (see Section 4.5). Therefore, there is no statistically significant evidence for a broad component associated with outflows in the stacked spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since outflows in various phases have been observed in extended regions in the outskirts of their host galaxies (e.g., R. Genzel et al. 2011, 2014; A. C. Fabian 2012; S. F. Newman et al. 2012a, 2012b; N. M. Förster Schreiber et al. 2014, 2019; A. K. Leroy et al. 2015; R. Herrera-Camus et al. 2019; R. L. Davies et al. 2020), we also looked into the extended emission around our sample. To do so, we investigate the spectra extracted from both larger apertures and annuli apertures centered on the galaxies. The larger apertures are designed by “dilating” the original mask by 10 kpc in every direction, while the “annuli” apertures are the difference between the bigger and original apertures. The 10 kpc increase in size is chosen based on the extent of the molecular gas outflow in R. Herrera-Camus et al. (2019). Finally, we perform the same stacking analysis and fitting as previously on these newly extracted spectra. While both methods do reveal spatially extended emission that is undetected in individual cubes and thus missed by the nominal masks employed, fitting the double-Gaussian profile to the lines returns the same result as for the initial aperture: we observe no faint broad component in the stacked profile, and the fit residuals show no improvement from performing a two-component Gaussian fit (see Appendix D for the stacked profiles from the extended and annuli apertures respectively).

We created the subsamples described in Section 3.3 based on observed correlations between ionized gas outflow incidence and galaxy properties. However, if molecular gas outflows have different correlations or are less prevalent in galaxies, a nonmatched choice of parameter space leaves too much “dilution” from sources with weak or absent outflow signatures. We apply the subset sampling analysis described in F. Stanley et al. (2019) to search for potential outflows in subsamples that may not match the subsets defined by M_⋆, ΔMS, or AGN as defined in Section 3.3 based on ionized gas outflow studies.

Following the methodology, we design a subset of spectra by randomly selecting between three and 154 galaxies from the full sample. These are then stacked and fitted with our two-component Gaussian model. Each spectrum in the subset is assigned a grade based on the strength of the broad Gaussian fit component (i.e., how much flux is contained in the secondary broad component, where the flux is defined as , with A_out the best fit broad component amplitude and σ_out is the best fit broad component width). The process is repeated 100,000 times, each time randomly selecting spectra from the full sample, such that every spectrum is graded. This allows us to select a final subsample of galaxies with the highest grade, i.e., the best candidates to have an outflow. We then repeat the stacking process on the galaxies with the highest grades and evaluate if we retrieve any underlying broad component.

Using this method, we find a subset of 41 galaxies (26% of the full sample, marked with a † symbol in Table A1) for which the stacked spectrum displays a tentative outflow signature, evaluated at 4σ, shown in Figure 5. The residuals of the best-fit show improvement when using a double-Gaussian model, and its reduced χ² is closer to one than that of the single-Gaussian fit. Repeating the analysis of the BIC presented in Section 4.2, we find ΔBIC = 6, which favors the double-Gaussian fit as better for the spectrum. All but one galaxy have coverage by either Chandra, MOSDEF, Spitzer, or VLA, and 13 are flagged as AGN (Section 2.2). We investigate the subsample properties in Appendix E and find no striking correlation between the tentative outflow detection and the galaxies’ physical properties, except for a slight bias toward higher masses (〈log(M_⋆/M_⊙)〉 = 10.95). This might be consistent with the trend in M_⋆ reported by, e.g., N. M. Förster Schreiber et al. (2019) in the case of predominant AGN-driven winds. We also investigate the outflow properties (outflow mass, mass outflow rate, mass loading factor) from the detection in Section 4.5 and report the values in Table 1. Still, the low significance of the result and lack of outflow signature in the stacks of high-mass and AGN-hosting galaxies prevent any firm conclusion with the data in hand. In addition, a caveat of the subsampling analysis is that it can result in selecting the galaxies whose noise properties will stack positively, mimicking an outflow detection. We try varying the rebinning of the spectra and stacking the spectra extracted from the extended apertures of the galaxies, and observe that the feature remains.

From our nondetections, we measure an upper limit on the flux an outflow could have based on the noise level measured in the empty stacked spectrum. Defining a detection above the noise level as 3σ_noise, where σ_noise is the standard deviation per channel of our empty spectrum (see Section 3.3), we compute the upper flux limit of a putative broad component using:

where N_ch = 2 × 3 × σ_out,ch is the number of channels spanned by this hypothetical Gaussian, defined as 3 standard deviations σ_out,ch (in channels) on both sides of the Gaussian center to encompass >99% of the flux.²⁰ As we do not have a value of σ_out, we compute this upper limit for values of σ_out between 0.07 and 10 times the narrow component width σ_ν,systemic (see Figure 6).

To investigate what this limit represents for the outflow and what parameter space this hypothetical outflow might span, we create a grid of 500 flux values for a set of broad-to-narrow amplitude ratio (0.01–0.8) and broad-to-narrow line width (0.07–10; which is also translated to physical values using the average channel width and average line amplitude of the stacked spectra). We plot this grid as the background of Figure 6, as a function of the input parameters. We note that for outflow widths below 2 times the galaxy line width, it will be almost impossible to separate the galaxy emission from the outflow emission, if the lines have the same central velocity (see also Section 5, Figure 7). On the other hand, if the outflow central velocity is shifted with respect to the galaxy emission lines, it should be possible to detect an outflow with a small width.

Zoom In Zoom Out Reset image size

This upper flux limit is computed by assuming the maximum possible flux value a stacked outflow component might have under the noise level. If all the galaxies in our sample have an outflow participating in the average—i.e., if the outflow incidence is 100%—the flux limit computed above corresponds to the gray shaded area in Figure 6. However, previous studies (e.g., N. M. Förster Schreiber et al. 2019) report lower incidence values, which also depend on stellar mass, redshift, and SFR. In this case, the stacking of spectra with and without outflow components will “average out” the strength of the potential outflows, decreasing the stacked flux. Consequently, the corresponding maximum flux of outflows in the sample increases as the incidence of outflows decreases. We show this effect by overlaying on the plot upper flux limits computed for incidences of 50%, 25%, and 12.5% (dashed white lines). This increases the flux limit as it would require the stacked outflow to have 2, 4, and 8 times, respectively, the flux of the outflow computed for 100% incidence. Finally, using the ionized gas outflow incidence values from N. M. Förster Schreiber et al. (2019), we compute the overall incidence of the sample by taking into account the stellar mass and distance from the main sequence of the galaxies in the sample. This corresponds to an overall 39.5% outflow incidence in the sample of 154 galaxies, plotted in Figure 6 as the solid red line. In addition, we add to the plot the results from the subset sampling analysis (purple circular marker), where the stack of 41 galaxies reveals a tentative outflow signature. Compared to the full sample, this corresponds to an outflow incidence of 26%, which is in marginal agreement with the upper limits derived from the noise.

For comparison with other molecular gas outflows studies, we also overlay on Figure 6 the positions of cosmic noon detections of molecular gas outflows in individual galaxies by C. Feruglio et al. (2017), M. Brusa et al. (2018), R. Herrera-Camus et al. (2019), and A. Vayner et al. (2021), in addition to two local Universe individual detections reported in C. Cicone et al. (2012) and S. Veilleux et al. (2017). All but the lowest redshift detection in Mrk 231 by C. Cicone et al. (2012) lie above the flux limit for 100% outflow incidence. This could be a result of the proximity of the target, which allows the detection of weaker molecular gas outflows. Unfortunately, very few CO detections of outflows exist at cosmic noon, and all originate from luminous AGNs or QSOs (except for the one reported in R. Herrera-Camus et al. 2019), which are not representative of the main galaxy population (see S. Veilleux et al. 2020 for a review).

5.1.1. Creating Mock Outflows

5.1.2. Output Flux Recovery

5.1.3. Comparison with the Real Data

A possible caveat of the method is rebinning before stacking the spectra. As mentioned in Section 3.2, stacking line profiles with a wide range of FWHM can result in the creation of artificial broad wings, which can be misinterpreted as an outflow component. This effect is illustrated in Figure 8 for the stack of 62 ΔMS > 0.2 galaxies, where we stack the spectra before and after rebinning to compare the effect of the varying line widths of the emission lines on the stacked profile. On the other hand, rebinning may reduce the ability to detect a broad component through dilution.

Zoom In Zoom Out Reset image size

To investigate how rebinning affects the retrieval of outflows in our sample, we try varying the values we rebin the spectra to, which all return similar results to those displayed in Figures 2 and 4, and where ΔBIC < 0 favors a single-Gaussian component fit in all cases. In addition, we also try stacking together subsamples of galaxies with similar line widths (in which case there is no need for rebinning) and observe no outflow component once again. Since the integrated line width can also serve as a proxy for galaxy inclination (for galaxies of the same mass, i.e., rotating at the same velocity, the more inclined galaxies’ line profiles will be wider), the stacking of galaxies with similar line widths without rebinning also ensures that we do not miss outflows based on inclination.

Finally, the procedure described in Section 5 in which we produce mock outflow emission, shows that rebinning does not prevent the fitting from detecting a broad secondary component when a clear signature is present. As a result, if the PHIBSS sample had broad outflow signatures present, which could be revealed by stacking, the rebinning process would not hinder our ability to do so.

We derive upper limits on outflow properties such as the mass outflow rate and the mass loading factor for all stacked samples. For the samples for which we do not have an outflow detection, we use as outflow widths and velocities the average values from ionized gas outflows reported in N. M. Förster Schreiber et al. (2019), extracted from their stacked spectra. We note that molecular gas outflows traced by CO are generally reported to be slower than their ionized gas counterparts, and thus, these computed values should be interpreted as very conservative upper limits. For the stacked spectrum from the subset sampling analysis in which we have a detection above the noise level (see Section 4.4), we define the outflow velocity as v_out = ∣Δv − 2σ_out∣, where Δv is the central velocity differential between the narrow and broad component, and σ_out is the width of the best-fit broad component (S. Veilleux et al. 2005; R. Genzel et al. 2011, 2014; A. Concas et al. 2022). All values are reported in Table 1.

To translate the relative limits discussed in Section 5 into physical limits, we compute the CO line luminosity in K km s⁻¹ pc² using the following relation (P. M. Solomon et al. 1997; S. Veilleux et al. 2017; L. J. Tacconi et al. 2020):

where the luminosity distance D_L is in Mpc, the observed frequency ν_obs is in GHz, and the integrated line flux S_COΔV is in Jy km s⁻¹. We use the median values of the galaxies that go into the stacking for the redshift, observed frequency, and luminosity distance.

This luminosity can then be converted into an outflow molecular gas mass estimate:

where r₃₁ = 0.77 is the ratio of temperature brightness to correct for the fact that we are not observing the CO(1–0) transition line but the CO(3–2) (L. A. Boogaard et al. 2020), and α_CO = 0.8 is the ULIRG-like H2–CO conversion factor in units of M_⊙ (K km s⁻¹ pc²)⁻¹ (S. Veilleux et al. 2017; A. Fluetsch et al. 2019). We use the ULIRG α_CO in this study, which is commonly used in the literature to obtain the outflow molecular gas mass, but this value relies on the assumption that outflows are well mixed and have high metallicity, which might not be the case. Other commonly invoked values of α_CO range from 0.3 for the optically thin case (e.g., A. D. Bolatto et al. 2013b; A. J. Richings & C.-A. Faucher-Giguère 2018) to 4.3 for the Milky Way value (A. D. Bolatto et al. 2013a). Studies in the local Universe indicate that α_CO in outflows is variable from case to case; however, the lack of consensus and the resolution of our data do not allow us to take this into account (see S. Veilleux et al. 2020 for a recent review). Since the molecular gas mass estimate is directly proportional to α_CO, this represents one of our estimates’ major sources of uncertainty.

Using the estimated molecular gas mass, we can compute an upper limit for the outflowing mass rate, using a simplistic yet plausible assumption commonly used in the literature (e.g., E. Sturm et al. 2011; F. Fiore et al. 2017; R. Herrera-Camus et al. 2019; with the details of the derivation in D. S. Rupke et al. 2005):

Here, R_out is the outermost radius reached by the outflow (we adopt the median effective radius of the galaxies stacked derived from HST imaging of the data (L. J. Tacconi et al. 2013), by lack of further spatial information and following what was done in S. F. Newman et al. 2012b and N. M. Förster Schreiber et al. 2019) and v_out is the outflow velocity. We compute for each stacked subsample.

Using this mass outflow rate, we compute the mass loading factor SFR, where SFR is the median SFR of the sample. The resulting values for each subsample are in the ranges and η_UL = 2.2–35.4. All values and average galaxy properties used for these estimates are reported in Table 1. In all cases, the derived upper limits on the mass loading factor allow for molecular gas outflows to carry a substantial amount of mass, enough to dominate the gas depletion in the galaxy. However, these values rely on numerous assumptions made to yield very conservative upper limits.

These limits are in line with some results from state-of-the-art simulations, which find that the cool (<10⁴ K) gas dominates the mass budget of AGN- and SF-driven outflows, with mass loading factors η ∼ 1–10 (T.-E. Rathjen et al. 2023; S. R. Ward et al. 2024). In more extreme cases, C.-G. Kim et al. (2020) and L. E. Porter et al. (2024) find η ∼ 100, which appears to be ruled out by our upper limits. However, the mass loading factors in C.-G. Kim et al. (2020) are measured at one scale height above the midplane, and the majority of the winds have very low velocities (∼10–100 km s⁻¹). As a consequence, the mass loading factors at larger distances are much lower and in better agreement with our upper limits.

However, in most cases, the simulations cannot resolve gas colder than 10⁴ K, limiting the conclusions that can be drawn on the dominant gas phase. Comparing with theoretical predictions, “bathtub models” usually suggest an η_out of unity (A. Dekel & M. R. Krumholz 2013). In addition, A. Dekel & N. Mandelker (2014) estimate that to match a simplistic bathtub toy model to observations, the net mass loading factor (the mass loading factor of the outflow minus the mass loading factor of recycled gas falling back in the galaxy) at z ∼ 2 should be 0, i.e., none of the outflowing gas escapes the galaxy, and all is recycled. This result can be in line with our upper limits on η as we only take into account the outflowing gas mass loading factor η_out.

The vast majority of cold molecular gas outflow detections are in the local Universe, found in starbursts (A. D. Bolatto et al. 2013b; A. Fluetsch et al. 2019), (U)LIRGs (R. Herrera-Camus et al. 2020b, 2020c; D. Lutz et al. 2020; A. Fluetsch et al. 2021) and AGNs (C. Cicone et al. 2012; F. Fiore et al. 2017; A. Fluetsch et al. 2019; see C. M. Harrison & C. Ramos Almeida 2024 for a recent review). The reported outflow velocities range from a few tens of km s⁻¹ (A. Fluetsch et al. 2019) to >10³ km s⁻¹ (F. Fiore et al. 2017; A. Fluetsch et al. 2019; D. Lutz et al. 2020), although most have velocities of the order 10² km s⁻¹ (F. Walter et al. 2002; A. D. Bolatto et al. 2013a; F. Fiore et al. 2017; A. Fluetsch et al. 2019; N. Krieger et al. 2019; D. Lutz et al. 2020). In line with the nearby Universe, the few CO detections at cosmic noon have similar velocities: M. Brusa et al. (2018) report an outflow detected with velocity ∼700 km s⁻¹, R. Herrera-Camus et al. (2019) find velocities between ∼300 and ​500 km s⁻¹, and the quasars in A. Vayner et al. (2021) host outflows at ∼400–1100 km s⁻¹. In general, all studies find that molecular gas outflow strength scales with AGN luminosity when one is present. The bulk of the mass budget in the outflows is dominated by the molecular gas phase, although there are no CO detections of SF-driven outflows at cosmic noon to compare with (see also I. Langan et al. 2025 in preparation). In contrast, outflows in warmer gas phases, such as ionized gas, display systematically higher velocities: for AGN-driven outflows, velocities are ∼1000–2000 km s⁻¹, against a few ∼100 km s⁻¹ for SF-driven outflows (A. Fluetsch et al. 2019; N. M. Förster Schreiber et al. 2019; D. Lutz et al. 2020). With these differences in outflow properties depending on gas phase and power source mechanism, it is not surprising that detecting molecular gas outflows in typical SFGs and low luminosity AGNs at cosmic noon remains a challenge.

However, the challenging nature of the detection may not be the only explanation for the lack of outflow detections in this study. Indeed, the cold molecular phase of outflows might not be the dominant one: most of the outflow could be in a warmer molecular phase as a consequence of energy injection by stellar/AGN feedback in the galactic interstellar medium (ISM). In the nearby Universe, studies have already investigated the hot molecular gas phase of outflows (≳1000 K; e.g., R. I. Davies et al. 2014; C. Ramos Almeida et al. 2019; R. A. Riffel et al. 2023) through observations of H₂ rovibrational lines, finding that the AGN could drive tens of solar masses per year of warm molecular gas, which is unlikely to affect large-scale properties of the host galaxy (R. I. Davies et al. 2014). More recently, JWST/MIRI has opened the door to the study of warm molecular gas at ∼100–1000 K through purely rotational H₂ lines in the mid-IR, finding low (∼150 km s⁻¹) outflow velocities (R. Davies et al. 2024; K. Y. Dan et al. 2025). In some cases, the molecular gas outflows probed in these local AGNs come from the intersection of the ionized gas outflow with the galactic disk, which perturbs the gas within the galaxy but will not expel it from the galaxy (C. Ramos Almeida et al. 2022; R. Davies et al. 2024).

The H₂ molecules in the gas could also be dissociated by the galaxies’ radiation fields into atomic hydrogen, particularly in MS galaxies. Using a simple theoretical model, A. Vijayan & M. R. Krumholz (2024) showed that the dissociation of molecular to atomic gas was mostly the result of radiative processes. Consequently, molecular gas in outflows can survive in starburst galaxies, thanks to the combination of denser environments and shorter dynamical times, which allows the molecular gas to escape the galaxy’s radiation field before getting dissociated. In more typical SFGs, however, the molecules get dissociated, and the cool phase of outflows is mostly composed of atomic gas. In this scenario, the absence of cold molecular gas outflow signatures in our sample might be explained by the fact that the majority of the outflow is in the atomic phase.

The neutral atomic phase of outflows has also been observed, using low-ionization interstellar absorption lines in the rest UV and optical (e.g., Mg ii; Na i D), both in the local Universe (e.g., K. H. R. Rubin et al. 2014; G. B. Zhu et al. 2015; S. Cazzoli et al. 2016; D. S. N. Rupke 2018; G. W. Roberts-Borsani & A. Saintonge 2019; C. R. Avery et al. 2022), and at higher redshift (e.g., A. E. Shapley et al. 2003; C. C. Steidel et al. 2010; D. K. Erb et al. 2012; Y. Sugahara et al. 2017). Specifically, using facilities such as JWST or Keck/LRIS, studies find prevalent neutral gas outflow signatures in cosmic noon galaxies (A. Fluetsch et al. 2019; A. Fluetsch et al. 2021; S. Belli et al. 2024; R. L. Davies et al. 2024; E. Kehoe et al. 2024; E. Taylor et al. 2024). In particular, S. Belli et al. (2024) find a significant neutral outflow component with velocity ∼200 km s⁻¹ in a poststarburst galaxy at z = 2.45, which could explain the rapid quenching of SF in the galaxy. In line with this finding, R. L. Davies et al. (2024) observe neutral outflows through Na i D absorption in z = 1.7–3.5 galaxies, with velocities 200–1100 km s⁻¹, which are found in star-forming and quenching systems alike. In a similar (1.7 < z < 2.7) redshift range, E. Kehoe et al. (2024) report neutral outflows detected through low-ionization interstellar absorption lines, albeit with lower mean velocities (〈Δv_LIS〉 = −56 km s⁻¹). These results indicate that, while the exact properties and impact on galaxy evolution are still uncertain, the neutral gas phase of outflows might play a significant role in the evolution of their host galaxy and represent a major fraction of the outflowing material.

Going along in the direction of these recent observations and theoretical models, T.-E. Rathjen et al. (2023) investigate the gas phase of stellar feedback using state-of-the-art magnetohydrodynamic simulations, where a nonequilibrium chemical network allows one to look into gas at temperatures down to ∼5 K. They find that 50%–90% of the outflowing material is in the so-called warm phase, defined as gas between 300 and 10⁴ K, which encompasses gas from the warm molecular gas phase to the ionized gas. As observations find that ionized gas does not dominate the mass budget of outflows (W. R. Freeman et al. 2019; N. M. Förster Schreiber et al. 2019; R. L. Davies et al. 2020), and in light of the recent studies discussed previously, this result seems in line with the picture in which gas phases other than the cold molecular and hot ionized gas are major contributors of the total mass budget of outflows.