© S. K. Stuber et al. 2021

1. Introduction

Understanding the mechanisms driving galactic outflows and their role in galaxy evolution is one of the challenges in modern astronomy and astrophysics. Recent theoretical works have suggested that galactic outflows are essential for regulating galaxy evolution (e.g., the shape of the galaxy stellar mass function; Scannapieco & Oh 2004; Silk & Mamon 2012) by quenching and morphologically transforming “star-forming” late-type galaxies to “red and dead” early-type galaxies (e.g., Di Matteo et al. 2005; Sijacki et al. 2007). Thus, outflows are thought to be capable of shutting down star formation by heating and dragging away molecular clouds that are ready to form stars (negative feedback; e.g., Sturm et al. 2011; Fabian 2012). On the other hand, outflows are also assumed to be able to enhance star formation by shock-driven compression and fragmentation of the surrounding interstellar medium (positive feedback; e.g., van Breugel & Dey 1993; Silk 2013; Maiolino et al. 2017).

Driven by active star formation, active galactic nuclei (AGN) or both, galactic outflows are considered a common feature of these galaxies. Galactic outflows driven by AGN activity are considered to play a crucial role in the coevolution of supermassive black holes and their host galaxies (e.g., Magorrian et al. 1998; Kormendy & Ho 2013; Costa et al. 2014). Also, they are the primary mechanism that pollute the circum- and intergalactic medium with metals. Therefore, they are strongly related to the chemical enrichment across cosmic history (e.g., Veilleux et al. 2005).

Galactic outflows are known to consist of multiple gas phases: hot X-ray gas (e.g., Komossa et al. 2003), warm ionized gas ([O III], Hα; Venturi et al. 2018; López-Cobá et al. 2020; Hogarth et al. 2021), neutral atomic gas (e.g., the Na I doublet; Sato et al. 2009), warm and cold molecular gas (OH, CO; Feruglio et al. 2010; Veilleux et al. 2013; Zschaechner et al. 2018; Lutz et al. 2020), and dense molecular gas (HCN, HCO⁺; Michiyama et al. 2018; Krieger et al. 2019). The molecular gas phase may dominate the outflow mass in certain active galaxies as revealed by multiphase observations of nearby galaxies, such as NGC 0253 (Krieger et al. 2019), quasar hosts at redshift z ∼ 2.4 (Carniani et al. 2015) or AGN hosts when comparing the molecular to their ionized gas phase (Fluetsch et al. 2019). This motivates a detailed characterization of the molecular phase of outflows.

Progress is being made with both numerical simulations and observational surveys to further quantify galactic outflows. As an example, recent cosmological simulations have predicted the statistical properties of galactic outflows in star-forming “main sequence” (MS) galaxies at different redshifts and their relation to observable properties of the host galaxies such as star formation rate (SFR) and stellar mass in the context of galaxy evolution (e.g., Nelson et al. 2019).

Observationally, thanks to advances in recent instrumentation, galactic outflow studies are now reaching statistically significant sample sizes of several tens or even several hundreds of galaxies (e.g., Veilleux et al. 2013; Concas et al. 2019; Förster Schreiber et al. 2019; Roberts-Borsani et al. 2020). These statistical works now provide typical outflow properties at given galaxy properties and redshift, although many of these lack spatial resolution sufficient to separate the non-circular motions that signpost outflowing gas from the rotation-dominated kinematics of the gas settled in the disk, except for very massive outflows with large velocities of several 100 km s⁻¹. Therefore, resolved and multi-phase studies are still rare and have been limited to individual galaxies: a few local MS galaxies (e.g., Combes et al. 2013; López-Cobá et al. 2019; Hogarth et al. 2021), local starburst galaxies (e.g., Bolatto et al. 2013a; Leroy et al. 2015), nearby luminous and ultraluminous infrared galaxies (LIRGs and ULIRGs; e.g., Feruglio et al. 2010; Cicone et al. 2014), as well as the Fermi Bubbles of the Milky Way (e.g., Su et al. 2010; Bordoloi et al. 2017). A more detailed overview on galactic outflows and recent studies is provided by Veilleux et al. (2020).

This implies that our understanding of the kinematic and morphological structures of galactic outflows is largely biased by nearby luminous galaxies. A major next advance will come in understanding the frequency and effect of nuclear outflows in “normal” star-forming galaxies. Even if these molecular outflows are not reaching the circumgalactic medium (CGM), they may play a crucial role determining the fate of the starbursts or activity in galaxy centers. Thus, in order to fully understand the relation between central outflows and their host galaxy’s SFR, high angular resolution observations are needed to see the interaction of the outflowing gas and its surrounding, as well as a large sample to understand statistical properties. Therefore, determining these properties for a representative sample of nearby “normal” galaxies, like MS galaxies, is a key objective. In this paper, we make a first step by estimating the frequency of molecular outflows in nearby MS galaxies.

For the purpose of a systematic search for central molecular outflows in nearby main sequence galaxies, ¹²CO (J = 2−1) (hereafter CO (2–1) line) is suitable, which is a fair tracer of bulk molecular gas. We use sensitive, high (∼100 pc) resolution CO (2–1) ALMA data for a large fraction of nearby, face-on, relatively massive MS galaxies obtained by the Physics at High Angular resolution in Nearby GalaxieS (PHANGS) survey¹ (Leroy et al. 2021b).

The PHANGS-ALMA CO (2–1) data products include position-position-velocity 3D data cubes as well as line-integrated 2D moment maps (Leroy et al. 2021a). Fully utilizing this 2D and 3D information allows us to identify molecular outflows. With these data we can search for outflows in galaxies spanning two orders of magnitude in stellar mass and derive the frequency of outflows based on a representative galaxy sample. Most important for this study, the sample selection does not contain any prior information about the presence of molecular outflows and nuclear activity, making the PHANGS sample unbiased for studying molecular outflows (Sect. 2). The PHANGS survey also includes an IFU survey with MUSE for 19 galaxies, which allows possible future studies to compare molecular and ionized gas outflow properties.

In this paper, we employ three popular methods to identify outflows that have been frequently used in the literature: 1D spectra (e.g., Feruglio et al. 2010), 2D position-velocity diagrams (e.g., Oosterloo et al. 2017), and channel-selected 2D integrated intensity maps (hereafter line-wing maps; e.g., Sakamoto et al. 2014). We describe both the strength and weakness of each method based on our evaluation experiences.

This paper is structured as follows: In Sect. 2, we briefly describe the PHANGS sample properties and ALMA data. Section 3 summarizes the methods used to identify outflows and our evaluation process. We present the results of the evaluation process in Sect. 4, including the final outflow candidates. In Sect. 5, we compare these candidates to the full sample, discuss possible biases and derive outflow rates which are compared to other samples. Finally, Sect. 6 provides a brief summary and conclusion of the main results.

2. Data and sample

The PHANGS-ALMA survey allows for the first time a systematic search for central molecular outflows in a representative sample of nearby galaxies on the star-forming main sequence. Here we briefly describe the PHANGS-ALMA survey and the characteristics of the data we use, and highlight information related to our molecular outflow study. Detailed descriptions are provided in other papers (Leroy et al. 2021a, and in prep.; see also Herrera et al. 2020 and Sect. 2 of Sun et al. 2020 for a comprehensive summary).

2.1. PHANGS-ALMA sample

PHANGS-ALMA is an ALMA Large Program (PI E. Schinnerer) to map the CO (2–1) line in the disks of 90 nearby galaxies when including its pilot and extension projects (Leroy et al. 2021b). This sample comprises close to all nearby (d < 24 Mpc), massive (9.3 ≤ log₁₀(M_⋆/M_⊙)≤11.1), relatively face-on (i < 75°), star-forming (log₁₀ (sSFR/yr⁻¹) > − 11) galaxies that are observable by ALMA (−75 ° < δ < +25°). Although PHANGS includes galaxies with high specific SFR, most starburst galaxies are excluded as they are rare in the nearby (z = 0) universe. For more details on completeness and data selection we refer to Leroy et al. (2021b).

We use 80 (out of a final sample of 90) PHANGS galaxies that were processed and available by September 2020 (i.e., internal data release version 3.4). An updated data release with improved continuum subtraction is publicly available on the homepage². Global properties of the sample galaxies³ are taken from Leroy et al. (2021b) with distances from Anand et al. (2021) and galaxy geometries from Lang et al. (2020) and are listed in Table C.1. Galaxies are classified as AGN based on Véron-Cetty & Véron (2010) and/or barred using the compilation by Querejeta et al. (2021), largely based on the work of Herrera-Endoqui et al. (2015) and Menéndez-Delmestre et al. (2007).

2.2. PHANGS-ALMA data

The PHANGS-ALMA observations used the extended 12m array, compact 7m array, and total power antennas in order to recover emission arising from all spatial scales within nearby galaxies. The data calibration and imaging is generally following standard procedures, we refer the reader for details on data reduction and product generation to Leroy et al. (2021a). The delivered data cubes and images have a median native angular resolution of 1.3″ translating into a spatial resolution of ∼100 pc, a spectral resolution of ∼2.5 km s⁻¹, and a typical rms sensitivity of ∼85 mK at their native resolution (Leroy et al. 2021b). This translates into typical molecular gas mass surface densities of 3σ ≈ 6 M_⊙ pc⁻² per channel (assuming a MW-like CO (1–0)-to-H₂ conversion factor of α_CO = 4.35 M_⊙ pc⁻² (K km s⁻¹)⁻¹ and a CO (2–1)-to-CO (2–1) line ratio of R₂₁ = 0.65 (Bolatto et al. 2013a; Leroy et al. 2013; den Brok et al. 2021).

We degraded the spectral resolution from 2.5 km s⁻¹ to 5 km s⁻¹ by smoothing the data to increase the signal to noise. We do so by applying Hanning smoothing using CASA (version 5.5) task specsmooth⁴ using dmethod=copy to avoid channel-to-channel correlation.

3. Methods

In this section, we present our procedure to systematically identify central molecular outflows in a large sample of galaxies with similar properties. Rather than case studies tailored to individual galaxies (e.g., NGC 0253; Krieger et al. 2019), we present a generic approach to identify galaxies which are likely to possess molecular outflows. This allows us to statistically compare properties of outflows (e.g., frequency, outflow masses, and outflow rates) as well as the non-detection rate of the sample in a robust, quantitative way. Such comparisons will provide a better understanding of which properties are necessary for galaxies to launch outflows.

We begin with a short review of previous identification methods in the literature (Sect. 3.1), before providing a detailed description of our applied methodology (Sect. 3.2).

3.1. Outflow identification in the literature

Identifications of outflows in the literature have used both spatial and kinematic information. As the use of spatial information requires favorable geometries of the outflowing gas, such as the presence of radio jets with a large and therefore more easily detectable extent of the outflowing gas, or a favorable disk inclination, (for instance, as in the case of M82 studied by e.g., Walter et al. 2002; Leroy et al. 2015), we focus on the kinematic approach. This approach leverages the fact that despite being less massive and thus less luminous than the brighter spiral arms, disk or central region, outflowing gas possesses a relatively high outwards velocity component, such that at a given projected position in a galaxy it is kinematically decoupled from the local rotational velocity. In order to leave the central disk of the galaxy, outflows must have velocities that are higher than the circular velocities of gas within the galaxy disk, but may not necessarily reach escape velocity of the galaxy halo. The outflow candidates we identify in this work are therefore consistent with (central) fountain flows. Tracing such low-mass gas components with anomalous velocities can be done by using one dimensional information (spectra for different apertures) or also considering spatial information in the form of position-velocity (pv) diagrams or line-wing maps (lwms; e.g., Cicone et al. 2014; Fluetsch et al. 2019; Lutz et al. 2020).

Outflows in several galaxies have been found via evidence for broad components in spectra from central regions observed via molecular or ionized gas tracers (e.g., Veilleux et al. 2013; Davies et al. 2020; Wylezalek et al. 2020). Often, Gaussian fitting is used to distinguish between the single- or double-horned spectral line expected for gas in rotating disks and an underlying broader extended component that indicates an outflow. However, bar-induced streaming motions, as well as fast-rotating central disks can result in a similar broadening of the spectrum (e.g., high velocity dispersions due to bar streaming motions are described in Kormendy 1982). Also, other non-circular motions (e.g., from spiral arms) may have a major contribution to the observed spectra, which emphasizes the importance of combining this method with the below listed methods. Furthermore, the movement of gas perpendicular to the line of sight or absolute projected outflow velocities that are consistently smaller than the velocities of more luminous components in the disk will not be visible in the spectrum taken over a central region, nor distinguishable from a spectrum taken over the full galaxy disk.

Finally, when using only the integrated spectrum, the extent and exact location of the outflow cannot be determined. These quantities are crucial to estimating timescales and mass outflow rates.

Another method to identify the high-velocity components arising from an outflow is via pv–diagrams along the kinematic major and minor axes. Although observations are limited to projected velocity components only (resulting in possible underestimation of the true three dimensional velocities), this method provides a rough estimate of the position and velocity of the outflowing gas.

The third commonly used method to detect outflows is via line-wing maps (lwm), which is widely used in identifying outflows from protostars (e.g. Dunham et al. 2014; Maud et al. 2015). The data cube is integrated over a previously determined velocity range – the velocities where the outflows appear distinguishable from the emission of the galaxy – and then presented as a map of integrated intensity (moment 0). This technique provides two advantages: First, the integration of intensity over several selected channels allows for locating the usually fainter gas, and secondly, this integrated intensity is then displayed in a two dimensional map, making it easier to see the outflow as a distinct component.

In this paper we applied the three methods described above to identify outflow signatures. The combination then provides the best possible inspection of data cubes for identifying outflow candidates.

3.2. Methodology applied in this paper

As we expect that most of our identifications would require follow-up observations to truly confirm the presence of an outflow, we will use the term outflow candidates instead of referring to them as outflows.

The finding charts of the three methods for each PHANGS galaxy (examples are shown in Figs. 1–3) were independently evaluated by three authors (S.K.S., T.S., and E.S., hereafter inspector A, B, C) and a ranking was generated, as defined in Table 1, for the possible presence of molecular outflows. This ranking number represents the outflow confidence label (OCL) and ranges from 0 (no detected velocity signature) to 3 (velocity signature found, interpreted as very likely outflow signature). An OCL of 3 is assigned when most of the criteria for each method (listed below) are fulfilled. On the other hand, an OCL of 0 is assigned when none of the criteria are fulfilled. OCLs of 1 or 2 correspond to intermediate cases, where only one to two criteria are fulfilled or the velocity signature is weak (for example, if the wings in the spectra are slightly broader than the usual double-horn profile, but are not very pronounced).

Fig. 1. Spectra of NGC 1566 and NGC 1546 used for visual classification of outflow candidates. Integrated and normalized spectra (black solid line) of the central regions, i.e., the inner 300 pc in diameter (top panels), and using annuli covering diameter with 300 pc < dgal < 2 kpc (bottom panels). For comparison, the integrated spectrum from the whole galaxy normalized to a peak of unity is shown by a gray line. The red solid lines represent the interpolated 3σ noise level for the central region/annulus. The mid-points of the 20% intersections of the peak (green triangle) are shown for comparison. Velocity has the systemic velocity subtracted. The transparent blue and red areas represent the velocity ranges over which the blue- and red-shifted line-wing maps (Sect. 3.2.3) are integrated over. These areas are selected to correspond to the velocity at 20% of the peak from the central spectrum out to the velocity at 5% of the peak plus an empirically determined additional 25 km s−1 to take into account faint emission at low significance. The green solid lines mark the 20% and 5% peak levels, respectively. Intersections with the central spectrum are given as red triangles. The spectra of NGC 1566 received OCLs of 3, 3, 2 for the 300 pc spectrum and 2, 2, 3 for the 2 kpc-annulus spectrum from all inspectors due to the pronounced line wings. Both images of NGC 1546 received an OCL of 0 from all inspectors for both apertures as no broad wings are evident.

Fig. 2. Example pv–diagrams along the major and minor kinematic axes of NGC 3627 (top) and NGC 4654 (bottom) used to search for outflows. The color scale represents the intensity at a given position along the offset averaged over a width of 3 pixels. The contours are multiples of the noise at 2.5σ, 3.5σ, followed by niσ with ni = 2ni − 1, for n0 = 2.5σ and i = 1, 2, 3, …. The color range is saturated to emphasize the faint emission we are interested in. The pv–diagrams of NGC 3627 (NGC 4654) correspond to OCLs of 3, 2, 2 (0, 0, 0) from the three different inspectors. The systemic velocity (yellow horizontal line) is added for comparison and is subtracted from the velocity axis.

Fig. 3. Line-wing maps for NGC 1566. Left panels: full moment-0 map available with the integrated lwms (red- and blue-shifted contours, respectively) overplotted. Right panels: zoom-in of the central region with 2 kpc diameter. The beam size is shown as a gray circle in the lower left corner of the zoom images. Center positions are indicated with a white plus. The lwms shown in the upper panels are derived by integrating over the velocity range selected from the 300 pc spectrum as described in Sect. 3.2.3, the lower panels use the 2 kpc annulus spectrum. A white circle indicates the size over which the corresponding spectrum is taken. These 4 panels together received an OCLlwm of 1, 1, 0 from the three inspectors. The contour levels are selected to be the same as for the pv–diagrams, with the noise being calculated for the zoom images, not for the full maps. The underestimation of the noise at the edges of the full map (left panels) is evident, however, it does not affect our search for central outflows.

We assign an uncertainty of ±1 to each inspector’s OCL assignment, as in some cases the inspectors found it difficult to decide between two neighboring OCL, as a result of the evaluation criteria listed below, which allow for some deviations. A full list of inspector-averaged OCLs per method and per galaxy are provided in Table D.1.

In the following sections, we describe the evaluation criteria for each of the methods used (spectra in Sect. 3.2.1, pv–diagrams in Sect. 3.2.2, and lwm in Sect. 3.2.3).

3.2.1. Spectra

Generation of the spectra. For the 80 galaxies that we consider, we extracted integrated spectra from circular regions (within the plane of the sky) centered on the galaxy’s center following previous works (e.g., Fluetsch et al. 2019; Roberts-Borsani et al. 2020; Roberts-Borsani 2020). The adopted center positions and systemic velocities are taken from Lang et al. (2020) as described in Leroy et al. (2021b) and listed in Table C.1.

Given that outflows are reported to extend between a few 100 pc to kpc-scales, we choose two apertures for the integrated spectra. One is a circular aperture on the sky with a diameter of 300 pc, the other one is an annulus with an inner diameter of 300 pc and an outer diameter of 2000 pc (in short hereafter “300 pc spectrum” and “2 kpc spectrum”). This should allow us to detect outflowing gas at the center of a galaxy, as well as to trace fainter outflows at kpc-scales without central emission. For comparison, we also extracted spectra integrated over the full disk. These spectra tend to show relatively higher noise levels than the spectra integrated over smaller areas, as we did not apply any masking, that could include additional biases.

We use circular apertures on the sky in this analysis because (1) galaxy inclinations are mostly low by survey selection criteria, thus the actual projected aperture in the galaxy plane will have a small ellipticity, and (2) gas emission is often asymmetrically structured, which will have a larger impact on the analysis than selecting the perfect ellipticity.

Evaluation of the spectra. The criteria for visual evaluation of the integrated spectra are as follows:

– Broad wings: We visually identify broad wings in the spectra as clear deviations from the expected shape of either a double-horned spectrum, a single-Gaussian spectrum, or shapes in-between. We do not fit the profiles, but even without fitting, all three investigators find OCLs with good agreement (see Sect. 4.1).

– Comparison to the spectrum of the whole galaxy: We compare the spectrum of the central region to the spectrum obtained for the whole galaxy as defined above. If the line wings of the central spectrum exceed the velocity extent of the whole galaxy spectrum at 5% normalized peak intensity, we take this as evidence for the existence of gas that is moving at enhanced velocities. We chose a threshold of 5% of the normalized peak intensity as it traces the spectral wings well, does not directly depend on the noise level, and exceeds the 3σ noise level for most galaxies.

– Absolute velocities: The absolute velocities reached in the line wings are used to identified gas with large velocities as could be generated by outflows. We visually estimate the velocity offset between the midpoint at 20% normalized peak intensity and the 5% normalized peak intensity of the central spectra. If this offset reaches ≳100 km s⁻¹, this is considered unusual and taken as a hint for the existence of an outflow. For most lower mass galaxies in our sample we do not expect to find gas moving at 100 km s⁻¹ or more, as the peak velocities of their (observed) rotation curves typically do not exceed this value (Lang et al. 2020).

We assigned an individual OCL for both apertures (see Sect. 3.2.1), as the spatial extent of the central outflow can vary and, thus, outflow signatures might only be visible in the central 300 pc spectrum or the 2 kpc spectrum, but not necessarily in both spectra.

We experimented with using Gaussian fits for an automated classification (e.g., Saito et al. 2018), but we found that the analysis of the results varied significantly due to large variations in spectral shapes seen in our sample of galaxies.

Two example spectra are shown in Fig. 1, displaying the effect of different apertures. NGC 1566 (left panels) shows broad wings that exceed the spectrum of the whole galaxy at 5% normalized peak intensity up to high offset-velocities in the 300 pc and 2 kpc spectra resulting in high OCLs. However, this galaxy is not part of our outflow candidates as it lacks outflow signatures in the other methods (see Sect. 4.1) demonstrating the importance of combining all three methods. The right panels of the figure display spectra of NGC 1546, where no broad component or other features typical for outflows could be visually identified resulting in OCLs of 0 from all inspectors.

3.2.2. Pv–diagrams

Generation of pv–diagrams.pv–diagrams are created using the CASA tool impv along the major and minor kinematic axes of all PHANGS galaxies using the position angles from Lang et al. (2020, see also Table C.1). Two examples are given in Fig. 2. The width of the slit should be approximately a beam width, and so we used a width of 3× the pixel size. We set the length of the slit up to 100 arcsec (∼1.5−12.5 kpc). We note that the noise in the data cubes varies mildly as function of frequency and towards the edges of the surveyed area (which can be seen in some pv–diagrams). As comparison of the central part to signal further out is most critical for pv–diagrams, the exact physical length of the slit is less important than for aperture-integrated spectra.

Evaluation of the pv–diagrams. Outflowing gas can be identified in a pv–diagram as a high velocity component at certain spatial locations. For a simple rotating disk, we expect a constant velocity (centered at systemic velocity) along the minor kinematic axis and an S-shaped rotation curve along the major kinematic axis. Deviations could be due to outflowing gas. However, galactic structures such as bars, spiral arms, or central disks can also modify the motions of the gas from simple circular rotation, making it difficult to distinguish between actual outflows and deviations due to other kinematic features (e.g., position angle dependence on bar observations; Koda & Wada 2002). We assigned a single OCL for major and minor axes pv–diagrams, as for most viewing orientations outflow signatures should be seen along both axes, due to its expected opening angle which would allow for projected velocity components.

We used the following criteria to assign the OCLs listed in Table 1 to all PHANGS galaxies based on the evaluation of their pv–diagrams:

– Faint emission with high velocity dispersion: We expect the outflow signature to be a faint component with a large extent in velocity. We consider emission at ≳100 km s⁻¹ from the mean value of the rotation curve at a given location as a strong sign for outflowing gas. We consider signals between 40 and 100 km s⁻¹ to be still significant, while components with even lower absolute velocities are not considered⁵.

– Comparing to galactic rotation and other rotation features: Presence of (faint) emission at anomalous velocities compared to the S-shaped rotation curve (along major axis) or a constant velocity (along minor axis) are telling signs for an outflow. Rigid rotation (e.g., from a stellar bar) can be visible as straight lines of strong emission up to high velocities in the center of the galaxy and must be considered.

– Relevance: As for this work we are only searching for central outflows, the potential outflow signatures should lie within the central 1 kpc (in radius).

Two examples are presented in Fig. 2: The S-shaped rotation curve along the major kinematic axis (left panel) as well as the constant velocity along the minor kinematic axis (right panel) of NGC 3627 (upper) and NGC 4654 (lower) can be seen. In NGC 3627, multiple features suggest an outflow, such as, faint red-shifted emission at roughly +16 arcsec along the major axis extending almost to ∼150 km s⁻¹ and faint red- and blue-shifted emission at +5 and +14 arcsec along the minor axis with a velocity extent of ∼100 and ∼150 km s⁻¹, respectively, possibly moving against the galaxy’s primary sense of rotation. These velocities are impossible to reach for gas in a purely rotating disk. The features found in NGC 3627 fulfill most of the criteria listed above making it one of our outflow candidates (see Sect. 4.1). NGC 4654’s pv–diagrams fulfill none of the criteria resulting in OCLs of 0. We show pv–diagrams of all outflow candidates in Appendix E.

3.2.3. Line-wing maps

Generation of lwms. The lwms allow us to localize emission at a certain velocity (range) in the galaxy plane. For this, we overplotted the emission integrated over a selected velocity range on the moment 0 map (which is integrated over the full velocity axis and has a threshold of 5σ) in order to compare to the overall gas distribution. The selection of the velocity range is based on the wings of the central spectra. Assuming that outflowing gas will be present at the tails of the line wings, we set a lower limit to the offset velocities at the intersection at 20% of the peak of the central spectrum. Lower velocities are considered to originate from galactic rotation. We set the upper limit to the offset velocities at the 5% intersections of the central spectrum, which for most of our data still lies above the noise level. To this velocity range we add an additional 25 km s⁻¹ to include faint emission. An example of the velocity ranges used for lwms can be seen in Fig. 1 as blue and red shaded areas, respectively.

This method is insensitive to outflows with velocities close to the systemic one (e.g. the outflow in the nearby galaxy NGC 253, see Bolatto et al. 2013a; Walter et al. 2017; Krieger et al. 2019), as the spectra are insensitive to such velocities as well, and thus the selected velocity ranges miss these velocities, but works well for outflows with a different velocity. This also provides additional information about the spatial location of the outflowing gas contributing to the high-velocity wings.

This procedure is applied to both the 300 pc and 2 kpc spectra. An example of the resulting two images and their close-ups for one galaxy is presented in Fig. 3. We add the integrated emission for the whole map, used to identify whether the high velocity emission is, for instance, something unique in the galaxy or rather corresponds to some other larger feature in the galaxy. For example, if the emission does not arise from outflowing gas, the high velocity emission simply traces the high-velocity part of the rotation curve and we might see emission in the form of a spider diagram in our map. Thus emission that is located on both sides of the center on two positions on a central ring and elongated in oval shapes along the central ring might not trace an outflow. Also, molecular gas co-located with a stellar bar may be experiencing bar streaming motions, rather than outflowing. Alternatively, there could be a narrow outflowing structure with a different inclination, but the same position angle as the bar. This second possibility is statistically rare, and so we do not consider this a possibility in our analysis.

Evaluation. We assigned a single OCL based on the lwms for each galaxy. As the spatial extent and location of the outflowing gas can vary between a couple and several hundreds of pc, we do not necessarily expect that outflowing signatures are present in both spectra. Therefore, we also do not necessarily expect that outflowing signatures are present in both lwms, as we use the velocity ranges from the two spectra. This means that not all of the automatically generated lwms will cover all of the outflowing gas and due to overlap in velocity range the resulting maps are also not independent.

The criteria for an outflow signature are as follows:

– No overlap with galaxy structures (location of emission): Comparing the location of the blue- and red-shifted emission to the overall galaxy structure such as inner ring, gas lanes in bars, and spiral arms allows us to exclude high velocity emission arising from non-circular motions within these galactic structures. For example, bars can induce strong shocks in their gas lanes which might be mistaken as an outflow.

– Comparing with overall rotation (velocity of emission): If the high velocity emission lies at “unexpected” velocities, such as finding red-shifted emission where blue-shifted emission is expected from the rotation curve, then this signature provides good evidence for the presence of an outflow.

– Spatial extent of the emission: Given that our sample galaxies only host low-luminosity AGN and have only moderate global SFRs, we expect low-mass outflows (visible as faint emission relative to our bright galaxy disks) with small to moderate extents.

– Relevance: If the emission of the potential outflow is located outside the central 1 kpc (in radius), it is not considered in this study.

4. Results

4.1. Resulting outflow confidence labels

We compare the OCLs assigned by different inspectors (A, B, C) for each method individually, and find an overall good agreement with Spearman rank correlation coefficients of ρ_S > 0.67. All coefficients are about equally correlated and coincide within 3σ. This allows us to derive reasonable average OCLs for each galaxy and each method by averaging across the inspectors.

We find the uncertainties of the Spearman rank correlation coefficients to be ≤0.07. We estimate these uncertainties as follows: First we perturb the OCLs per inspector per method randomly within their adopted uncertainty of +1, 0, or −1. Next, we calculate Spearman rank correlation coefficients for the perturbed set of OCLs for all combinations of two inspectors for each method (jackknife). Then we repeat this a thousand times and calculate the standard deviation of the resulting Spearman values. The Spearman rank correlation coefficients for the different combinations of inspectors per method and their uncertainties are listed in Table 2.

The agreement of OCLs between different methods is slightly less good, but within 3σ. This is expected, as the different methods are sensitive to different outflow characteristics (see Sect. 3) which can result in deviating OCLs. Next, we define as the method-average of the inspector-averaged OCLs of the pv–diagram, the lwm, and the better ranked one of the spectra (see motivation in Sect. 3.2.1). To ensure consistent treatment of the spectra-OCLs with the other two methods, we require agreement of the inspectors on the same aperture and use the higher ranked inspector-averaged spectrum-OCL.

The resulting inspector-averaged OCLs for each method as well as the for each galaxy are reported in Table D.1. Figure 4 displays the (black cross) as well as the inspector-only-averaged OCLs per method per galaxy (colored filled symbols). Galaxies with are considered outflow candidates according to the definition of the OCLs. Uncertainties for the OCL and are derived by perturbing the OCL per method per inspector randomly within their adopted error of +1, 0, or −1, and performing jackknifing for each method. When averaging, this results in different inspector-averaged OCLs per method as well as different inspector- and method-averaged for each perturbation. We repeat this 1000 times to get the standard deviation of both the OCLs per each galaxy per each method as well as for .

Fig. 4. Average outflow candidate OCLs for PHANGS galaxies. Black crosses represent the per galaxy, which is the average of the inspector-averaged OCL of the pv–diagram, lwm, and the higher ranked spectrum. Galaxies with OCL ≳ 2 are highly likely to exhibit strong central outflow signatures, whereas OCL ≲ 1 correspond to an absence of any outflowing gas signatures. For each galaxy the inspector-averaged OCLs of the individual methods are shown as red circle (central spectrum with 300 pc aperture), purple square (spectrum of annulus with 2 kpc outer diameter), blue triangle (pv–diagram), and light-blue diamond (lwm). The green-shaded area marks all 2 ≤ OCL ≤ 3. A galaxy with (and lwm-OCL ≥ 2) is considered an (secure) outflow candidate and emphasized with green galaxy names (see Sect. 4.2). Uncertainties are derived by perturbing and jackknifing the OCL per method per inspector, as described in the text.

4.2. Frequency of central molecular outflows

In order to determine the frequency of central outflows, we need to establish a set of outflow candidates based on their OCL values. Given our initial definition for assigning an OCL (Table 1), we require that the average OCL ≥ 2, to ensure that possible outflow signatures are visible in all methods (see below). A total of 20 galaxies, or outflow candidates hereafter, meet this requirement resulting in an frequency F of . The statistical uncertainty has been obtained similarly to the uncertainty of the inspector-averaged OCL and the . During each of the 1000 perturbations of the OCLs, we additionally calculate the frequency according to our criteria. From the resulting 1000 frequency estimates, we can take the standard deviation as uncertainty.

For all outflow candidates, their as well as host galaxy properties (stellar mass, SFR, central SFR (cSFR), presence of bar and/or AGN) are listed in Table 3.

Secure outflow candidates. In some cases, like NGC 4321 or NGC 1097, the putative outflow emission seen in the spectra and pv–diagrams could also be explained via shocked regions or bar flows evident in the lwm. Thus by adding a second criterion, we can increase the weight of the lwm-OCLs, which might be more sensitive to the contribution of streaming motions. Thus for secure outflow candidates we require in addition to average that the lwm-OCL is ≥2.

These two criteria are fulfilled by 16 out of 80 galaxies and result in a secure frequency estimate of & OCL_lwm ≥ 2) = 20 ± 2 %. The four outflow candidates with lwm − OCL < 2 are indicated with a star (^⋆) in Table 3. We determine the statistical uncertainty using the same approach as before, but adding the second criteria to the calculations.

Further, we additionally inspect the data cubes of the outflow candidates for the presence of outflow emission without integrating or further reducing the available information per axis confirming the importance of the second criterion for the selection: In all cases where lwm − OCL < 2, emission in the data cube likely represents bar emission or shock emission rather than outflow-like emission.

To unambiguously identify actual outflows among our outflow candidates, an individual analysis of the putative outflow geometry and the galaxy disk will be carried out in a future work. This involves analysing whether the putative outflowing gas can be found in counter-rotating regions, that could only be due to a bulk motion out of the plane. Also, we plan on comparing the outflow position angles to those of galactic stellar bars, which could reveal more insights on the driving mechanisms, as AGN-driven outflows are generally randomly orientated compared to star-formation driven outflows which are usually orientated along the minor axes. This is, however, beyond the scope of this work that focuses on statistical results rather than individual galaxy properties.

5. Discussion

Using standard outflow detection methods based on visual classification, we find the frequency of outflows to be F = 20−25%±2%_stat depending on the number of criteria applied. We divide the sample into outflow candidates and non-outflow galaxies using an average -threshold (and additional an lwm − OCL) of 2 and higher. The outflow frequency as function of -threshold in steps of is given in Appendix A for interested readers.

In Sect. 5.1, we discuss possible differences between outflow candidates and the remaining PHANGS sample galaxies. We test possible relations between the OCLs and host galaxy properties and the reliability of the determined frequency in Sect. 5.1.2. In Sect. 5.2, we estimate outflow masses and mass outflow rates for our candidates and compare them to literature results (see Sect. 5.3).

5.1. Outflow candidates versus non-outflow galaxies

Among our 80 targets, we find 20 outflow candidates with 16 classified as secure candidates. We compare the outflow candidates to the full sample with respect to their host galaxy properties.

Figure 5 displays the fraction of AGN and bars among all and the secure outflow candidates, and the full sample. Out of 20 (16 secure) outflow candidates, 50% (50%) contain an AGN, and 89% (87%) host a stellar bar, which are both significantly higher than the global fraction (24% AGN out of 79 galaxies with AGN/no-AGN classification, 61% barred out of 72 galaxies with good morphological bar/no-bar classification). Thus, we find an enhanced fraction of (secure) outflow candidates among our 19 AGN host galaxies of 53% (42%) compared to 17% (13%) for the 60 galaxies without an AGN. The same applies for barred galaxies, where a higher fraction of (secure) outflow candidates among our 44 barred galaxies of 39% (30%) can be to a fraction of 7% (7%) for the 28 unbarred galaxies.

Fig. 5. Fraction of AGN (violet, left) and bars (blue, right) among PHANGS outflow candidates and the secure outflow candidates, as well as the entire PHANGS sample.

This large fraction of AGN, with about half of all AGN being an outflow candidate, agrees well with the results of Rakshit & Woo (2018), who found that 45% of Type 2 AGN as well as 89% of Type 1 AGN possess a broad wing component in their ionized gas spectrum. Our AGN sub-sample consists of mainly Type 2 AGN.

The preference for barred galaxies among our molecular outflow candidates has not been reported before. This could be due to covariance with galaxy mass, given that both the bar fraction as well as the AGN fraction increase with stellar mass (e.g., Rosas-Guevara et al. 2020; Consolandi 2016; Gavazzi et al. 2015). This higher fraction is likely to be related to the findings that barred galaxies show a higher concentration of molecular gas in their central ∼1 kpc region compared to unbarred galaxies (e.g., Sakamoto et al. 1999; Sun et al. 2020), resulting in a reservoir for the outflow that otherwise might not be present.

To analyse this further, we search for a possible influence of host galaxy properties such as stellar mass, SFR, and central SFR.

5.1.1. Role of host galaxy on outflow candidate selection

As the PHANGS-ALMA sample is selected to represent normal main sequence galaxies at z = 0 (see Sect. 2), we check for trends with SFR, central SFR (within 2 kpc in diameter), and stellar mass (global properties from Leroy et al. (2021b), see same Sect. 2) among the outflow candidates and in the ranking. Figure 6 displays the distribution of s and (secure) outflow candidates on the star-forming main sequence. It also shows the ratio of central-to-global SFR for (secure) outflow candidates as a function of stellar mass and highlights galaxies hosting an AGN.

Fig. 6. (a) Distribution of PHANGS outflow candidates (secure candidate: red hexagon, otherwise: orange star) and non-outflow galaxies (gray hexagon) along the star-forming main sequence. The histograms show the distribution over stellar mass (top) and SFR (left) for outflow candidates (red, all candidates combined) and galaxies without outflows (gray). The main sequence according to Leroy et al. (2019) is added for comparison (gray solid line) with a width of 0.5 dex (gray shaded area) and extrapolated to larger and smaller stellar masses (dotted gray line). (b) Distribution of PHANGS galaxies along the main sequence of star-forming galaxies. The is shown by color. A dark color represents a large (e.g., class 3), thus showing host galaxies with high probabilities to possess molecular outflows. Light colors represent low (e.g., class 0), thus no or only tentative outflow signatures are present in these galaxies. The histograms show the distribution over stellar mass (top) and SFR (right) for three bins. (c) Ratio of central-to-global SFR as a function of stellar mass for PHANGS galaxies likely to host a central molecular outflow (secure candidate: red hexagon, otherwise: orange star) or not (gray hexagon). Histograms show the distribution of all outflow candidates (red) compared to galaxies without outflow signatures (gray) in stellar mass (top) and SFR ratio (left). AGN are marked by an open black diamond. For all plots the median uncertainty for individual data points is shown in the lower right corner.

Our findings are as follows:

a) All molecular outflow candidate galaxies have stellar masses of log₁₀(M_⋆/M_⊙)≳10 and SFRs of log₁₀SFR/, covering roughly the upper halves of the respective ranges spanned by the full sample. This is also clearly evident from the histograms (outflow candidates (red) versus remaining galaxies (gray)). We discuss the possible existence of a minimum threshold in stellar mass and/or SFR, before a molecular outflow can be launched in a main sequence galaxy, in Sect. 5.1.2. Figure 7 shows the frequency of outflows in bins of stellar mass which we find to increase monotonically as a function of stellar mass over the stellar mass range probed by our sample. We use stellar mass bins of width Δlog₁₀(M_⋆) = 0.5 here, but the same trend can be seen for smaller mass bins (e.g., 0.33 or 0.25). Unfortunately, we want to emphasize that the sample size at each bin is small so that possible influences from AGN and bars cannot be investigated at this point. As galactic molecular outflows are thought to play a critical role in regulating the galaxy stellar mass function, the possible correlation between molecular outflow frequency and stellar mass found in local main sequence galaxies represents a future direction for studying galactic molecular outflows.

Fig. 7. Fraction of (secure) outflow candidates as a function of stellar mass. The fraction is determined per stellar mass bin with widths of Δlog10(M⋆/M⊙) = 0.5 that contain a minimum of 3 galaxies. Each level represents the fraction of outflow candidates relative to all galaxies in the mass bin. The total amount of galaxies per mass bin is given at the bottom of each bin.

b) We observe a weak trend in the SFR versus stellar mass plane where galaxies with high stellar mass and high SFR appear more likely to possess a higher value than galaxies with low stellar mass and low SFR indicated by Spearman rank correlation coefficients for SFR (stellar mass) and s of ρ_SFR = 0.40 ± 0.03 (ρ_M⋆ = 0.56 ± 0.03). Uncertainties of the Spearman rank correlation coefficient are determined via jackknifing. We perturb the OCLs as well as the corresponding stellar mass and SFR values by their uncertainties, and calculate the Spearman rank correlation coefficient. We repeat this 1000 times and adopt the scatter among the Spearman rank correlation coefficients as uncertainty. In Sect. 5.1.2, we will come back to this topic and search for secondary dependencies, namely the presence or absence of a stellar bar or AGN.

c) There is a strong preference for outflow candidates (∼60%) which have a significant (> 10%) fraction of their total SFR confined to the central 2 kpc in diameter (Fig. 6). Such a trend that the hosts of our outflow candidates exhibit higher fractional central SFR compared to the galaxies without outflows is not surprising as central SFR tends to correlate with global SFR (see point above). Another possible explanation might be that the central SFR is the more fundamental quantity that enables outflows to be launched compared to the global SFR. Similarly to the global SFR, we find a threshold in the central-to-global SFR ratio fulfilled by our outflow candidates (see discussion in Sect. 5.1.2).

5.1.2. Correlation of OCL and host galaxy properties

We explore here how much the trends seen pertain to outflow host galaxies or present biases in our classification methodology. Overall, we see a trend along the main sequence with significant (ρ > 3σ) Spearman rank correlation coefficients between and global SFR as well as stellar mass.

Since we utilize CO to identify outflows in our sample galaxies, host galaxy properties affecting the amount of molecular gas might lead to a potential bias. A higher stellar mass (implying a higher central stellar mass surface density), likely leads to higher molecular gas mass surface density (ρ = 0.67 for stellar mass and molecular gas mass surface density) and therefore results in higher signal-to-noise observations which increase the chance to identify outflowing gas. In fact, we find a strong correlation between and molecular gas mass surface density (ρ = 0.81).

Similarly, if the conversion factor of CO luminosity to total molecular mass, α_CO, is different in the outflowing gas relative to the galaxy, it may possibly lead to a correlation between and stellar mass. However, α_CO in outflowing gas is currently unconstrained. On the other hand, a higher molecular gas mass surface density can increase the incidence of molecular outflows as they require a molecular gas reservoir to be present. As an example, Hogarth et al. (2021) find that edge-on galaxies that host molecular outflows seem to have their molecular gas more centrally distributed than in their control sample. The correlation between and stellar mass is furthermore supported by recent studies that reported higher stellar mass galaxies to have more and faster outflows (if an AGN is present, e.g., at z = 0.6−2.7; Förster Schreiber et al. 2019), thus a higher stellar mass can actually increase the chance of (detectable) outflows to be launched.

Similarly, Cicone et al. (2016) suggest that for galaxies with stellar masses comparable to those in our sample, a higher SFR increases the incidence of ionized gas outflows. This may suggest an increase of ionized gas outflow incidence with offset from the main sequence, for example, as it was found for galaxies at z = 0.6−2.7 (Förster Schreiber et al. 2019). Although ionized outflows might be more common when molecular gas is rare or not present, the underlying launching mechanisms out of the galaxy will have some similarities, making such ionized gas studies complementary. We do not find an increase in outflow velocity with larger offset from the main sequence (ρ_ΔMS = 0.24, p-value > 0.2) as they do for their candidates. This might be a consequence of the PHANGS sample not extending far above or below the main sequence.

The lack of correlation between galaxy distance and (ρ_dist = 0.164, p-value ∼0.15) argues against any resolution dependency for the OCLs or the relations found.

The observed trend between and both stellar mass and SFR highlights the importance to extend outflow studies to a larger range of main sequence galaxies. Up to now, outflows have been preferentially identified in starburst galaxies with stellar masses and SFRs in excess to those of the PHANGS galaxies (Fig. 10; e.g. Cicone et al. 2014; Fluetsch et al. 2019; López-Cobá et al. 2019; Roberts-Borsani 2020; Roberts-Borsani et al. 2020). Although we see strong evidence for a minimum threshold in stellar mass and/or SFR for a main sequence galaxy to possess a (detectable) molecular outflow (e.g., Fig. 6), our data do not exclude the possibility of a continuous relation between frequency of outflows and stellar mass and/or SFR. We infer a stellar mass threshold of M_⋆ ≳ 10^10.0 M_⊙ and SFR threshold of SFR ≳ 10^−0.5 M_⊙ yr⁻¹, or fractional central SFR of cSFR/SFR ≳ 10^−1.9. Such thresholds have been found for ionized gas outflows as well, with a limiting stellar mass of M_⋆ ∼ 10^10.9 M_⊙ (e.g., Genzel et al. 2014) below which the occurrence of outflow signatures drops sharply. If a stellar mass threshold exists, it would be more sensible to infer the actual outflow frequency for the stellar mass range where outflows can occur. For our sample with log₁₀(M_⋆/M_⊙) ≥ 10.0, this corresponds to a frequency of outflow candidates of 36%. Further, this would imply that mechanisms capable of launching molecular outflows are less effective in galaxies of intermediate stellar mass.

Regarding the correlations found with the ratio of central to global SFR, there is a longstanding idea in the literature that outflow velocity correlates strongly with a higher SFR per area, Σ_SFR. This idea, that instead of global SFR the SFR per area has a larger influence on the presence of outflows was introduced by Heckman et al. (1990, 2015). In addition, Heckman (2003) suggested a threshold of SFR per unit area above which outflows appear ubiquitous. This might be similar to the threshold for central to global SFR found above.

There is still the possibility that our threshold is only a coincidence and outflows appear in galaxies spanning a wider range of stellar mass and SFR, (e.g., over smaller area and/or with smaller outflow velocities) and our outflow detection methods are only sensitive to the high end of SFR and stellar mass. Studies such as López-Cobá et al. (2019) analysing ionized gas emission from 273 highly inclined galaxies find no stellar mass threshold, but outflows over their full stellar mass range 9.5 < log₁₀(M_⋆/M_⊙) < 11.

To investigate whether other host galaxy properties affect the correlations of with stellar mass and SFR, we analyse the presence/absence of a bar or AGN. The full analysis is provided in Appendix B, here we briefly summarize our findings: Overall, barred galaxies show a stronger correlation between and SFR compared to non-barred galaxies, and a similar trend is seen for AGN compared to inactive galaxies. It could be plausible that these trends arise due an enhanced bar and AGN fraction at higher SFR in the sample, and the relevant correlation is between our and SFR. In contrast, we find no significant difference between stellar mass and for galaxies with or without AGN and/or bar. Therefore, the role of bars and AGN in enhancing the outflow incidence above the dependence on stellar mass remains ambiguous for our galaxy sample.

To test if all methods are equally impacted by stellar mass, SFR, and the presence of either bar or AGN, we calculate Spearman rank correlation coefficients for OCLs of the individual methods (300 pc spectra, 2 kpc spectra, pv–diagrams, and lwms, each averaged over the different inspectors) and stellar mass (SFR) for all galaxies (for details and ρ values see Appendix B.1). Overall, we do not find a significant difference to the previously found correlations. As all methods trace gas kinematics in a different way, this is reassuring and lowers the chance of a possible inherent bias. This increases our confidence that the trends between stellar mass and outflow incidence or SFR and outflow incidence are real.

5.2. Outflow masses, rates and other properties

In this paper we have identified outflow candidates in PHANGS using three different diagnostic tools. We treat them as candidates only, as a number of dynamical processes may mimic outflow signatures; follow-up with multi-wavelength data is required to unambiguously confirm the presence of outflowing gas. Nevertheless, here we derive outflow rates (Ṁ_outflow in Table 3) based on some simple assumptions for an initial characterization of these potential outflows.

First, we define an outer boundary to the outflow radial extent () in the plane of the sky from visual inspection of the different diagnostic plots and list them in Table 3. Next, we construct moment maps considering emission in the blue and red line wings of the ALMA cube separately, following the automated criteria described in Sect. 3.2.3 for a 300 pc aperture. For the range of channels in each wing (blue or red), we integrate the total molecular mass inside (i.e., we sum the moment 0 map constructed without any thresholding over these channels). We also derive the CO-weighted average radius () in the plane of the galaxy, and the CO-weighted moment 1 velocity (), again inside and over each line wing. We discuss the connection between 3D outflow geometry and the here defined projected radius and line-of-sight velocity below. Both for the weighting and the computation of the moment 1 map, we use a thresholded version of the ALMA cube (4σ clipping).

With these quantities, we estimate a characteristic mass outflow rate as:

(1)

for each of the wings, and add up both contributions (). We note that this is representative of an outflow made up of thin expelled shells or clumps averaged over time (Veilleux et al. 2005; Fluetsch et al. 2019). For a spherical or multi-conical geometry, the outflow rates would be three times higher (e.g., Maiolino et al. 2012; Cicone et al. 2015). In the conversion from CO luminosity to total molecular mass, we adopt a value of α_CO = 0.8 M_⊙ (K km s⁻¹ pc²)⁻¹, typical of starburst galaxies and often assumed for molecular outflows (Bolatto et al. 2013b; Lutz et al. 2020; Veilleux et al. 2020).

The mass outflow rates in Table 3 are solely based on the velocity projected along the line of sight, whereas their extent is measured in the plane of the sky. Generally, constraining the angle of the outflow is not straightforward, and we do not attempt to make that correction here (thereby matching most literature), as it requires a dedicated case-by-case study. Since we measure the outflow radius vector and velocity vector as projected onto the plane of the sky and the line of sight, respectively, correction factors of sin α and cos α need to be applied to recover deprojected values, where α denotes the angle between the line of sight and radius vector or velocity vector. The correction factor on the outflow rate is then sin α /cos α = tanα (see e.g., Krieger et al. 2019). For example, deviations of up to 30° from α = 45° result in correction factors of up to ∼4. Therefore, in the most optimistic case, the individual outflow rates could be considered accurate within a factor of a few. In pathological cases where the outflow is extremely close to the line of sight or nearly perpendicular to it, the corrections can be much larger. Fortunately, when averaging mass outflow rates over an ensemble of galaxies, the projection corrections should cancel out (because the orientation of the outflows are random relative to the line of sight; Cicone et al. 2015). When averaging over the 20 outflow candidates in Table 3, we find a modest mean mass outflow rate of ⟨Ṁ_outflow⟩ ∼ 3 M_⊙ yr⁻¹.

Another concern when studying such modest outflows is the difficulty to separate all of the outflowing gas from the gas rotating in the galaxy disk. While our definition of line wings attempts to do this, there could be a range of velocities in the outflow (e.g., due to projection), where some CO emission associated with the outflow might be indistinguishable from the rotating disk. In such a case, our estimates from the line wings should be regarded as lower limits.

Finally, we also checked the sensitivity of our results to the outer radial outflow boundary, which we estimate by eye. Perturbing by ±2″, roughly twice the typical resolution of the ALMA data, results in outflow rates that differ on average by 30%, which is a modest variation. For outflow rates, we therefore assign a statistical uncertainty of 0.2 dex. However, the actual uncertainty may be larger given the uncertainty contributions discussed above and might be closer to ∼1 dex.

Figure 8 shows the mass outflow rates calculated according to Equation (1) as function of SFR for PHANGS outflow candidates as well as a compilation of known molecular outflows by Fluetsch et al. (2019) that use the same α_CO conversion factor as well as the same mass outflow rate calculation (Eq. (1)). We indicate where the mass loading factor η = Ṁ_outflow/SFR is unity (gray line). Overall, we find good agreement between both outflow samples, with the PHANGS outflow candidates further expanding the trend to lower mass outflow rates.

Fig. 8. Mass outflow rates Ṁout as function of SFR for the secure outflow candidates (red circle) and possible outflow candidates (orange star) from PHANGS as well as the outflows compiled by Fluetsch et al. (2019) (blue diamond, blue errorbar). The gray line corresponds to a mass loading factor η of unity. PHANGS SFRs have an uncertainty of 0.11 dex and the mass outflow rates have statistical uncertainties of 0.2 dex when perturbing the radius the outflow mass is determined over, with a larger total uncertainty estimated to be ∼1 dex given the unknown actual outflow geometries and properties.

In Fig. 9, we show the mass loading factor η as function of stellar mass for the PHANGS outflow candidates and the literature sample from Fluetsch et al. (2019). The PHANGS outflow candidates have lower mass loading factors and lower stellar masses, but overall no significant difference between both samples is found. The average mass loading factor of the PHANGS sample is log₁₀⟨η⟩ = 0.0, with a statistical uncertainty of 0.23 dex and an estimated total uncertainty of ∼1 dex, indicating that SFR and outflow rate are comparable. Put in another way, the outflows remove as much gas from the galaxy center as is consumed by central SFR when present, and so they can represent an important contribution to the overall future evolution of the gas. We do not find a significant difference between PHANGS outflow candidate galaxies that have an AGN or a bar compared to those that do not. Using the central SFRs derived over the central 2 kpc instead of global SFRs will result in a significant increase of the mass loading factor. However, more meaningful estimates of this factor might be derived from resolution matched SFR maps that allow for tailored measurements in a future work.

Fig. 9. Mass loading factor η as function of stellar mass for PHANGS outflow candidates (secure candidate: red circle, otherwise: orange star) and outflows compiled by Fluetsch et al. (2019) (blue diamond, blue errorbar). The PHANGS outflow candidates populate the lower stellar mass and lower mass loading part of the probed parameter space, but overall no significant difference between both samples is found. PHANGS stellar masses have an uncertainty of 0.11 dex and the statistical uncertainty of the mass loading factor due to error propagation is estimated to be 0.23 dex (red errorbar) with a larger total uncertainty of up to ∼1 dex.

In conclusion, while the outflow rates are subject to large uncertainties, our measurements suggest that these are modest molecular outflows (if the kinematic signatures that we observe indeed trace only outflowing molecular gas). Such outflows might not be able to quench a whole galaxy, but can determine the fate of starbursts or the evolution of galaxy centers. A detailed follow-up study with multi-wavelength data would be necessary to confirm which of these outflow candidates host true outflows and to refine the calculation of individual mass outflow rates. Also, searching for differences in outflow properties such as extent, mass, energetics, and geometry between SF- and AGN-driven outflow will require a significant and well-characterized sample. The molecular outflows of the PHANGS outflow candidates have an average extent of pc, velocity of km s⁻¹ (see also Table 3), and integrated CO luminosity of ⟨L_outflow⟩∼8.6 K km s⁻¹. These numbers emphasize the importance of high resolution, high sensitivity studies such as PHANGS, that enable us to visually detect these outflows on sub-kpc scales. This must be considered when comparing this study to previous results from the literature (see Sect. 5.3).

5.3. Literature comparison

5.3.1. Individual galaxies

To gauge the robustness of our outflow classification methods, we searched the literature for detected outflow signatures in molecular gas for our PHANGS galaxy sample. We separately discuss our outflow candidates and the rest of the sample.

For 7 out of 16 of our secure outflow candidates either molecular and/or ionized central gas outflows have been reported in the literature. Four studies on our candidates explicitly exclude outflows in cold molecular gas. For the remaining secure candidates, no previous studies on molecular outflows were found.

Out of the four candidate outflow galaxies that are marked with a star (^⋆) in Table 3, thus not fulfilling the stricter criterion, none had a cold molecular outflow study in the literature and only some evidence for other gas phases (e.g., warm molecular gas or ionized gas).

For four of the remaining non-outflow candidate galaxies studies report evidence for molecular/ionized outflows. However, for three of them and for three additional non-outflow candidate galaxies also studies are available that present evidence against the presence of outflows. For the remaining galaxies in the non-outflow candidate sample, no studies on molecular outflows were found in the literature.

In summary, we find reasonable agreement for most of our outflow candidates, given that for nearly half of the galaxies no molecular outflow study has been made, and the existing differences in resolution, gas phase, and methodology applied to identify the outflows. Our study may have missed one potential outflow candidate among the PHANGS galaxies, namely NGC 1068, and for two of our candidates (NGC 1365, NGC 1672) high resolution CO (3–2) observations of the galaxy centers provided no clear evidence for molecular outflows.

PHANGS outflow candidates and literature matches.

– NGC 0253: Its central molecular outflow is long known and has been studied in detail (e.g., Sakamoto et al. 2006; Bolatto et al. 2013a; Walter et al. 2017; Zschaechner et al. 2018; Krieger et al. 2019).

– NGC 1365: Combes et al. (2019, CO (3–2) emission) do not find clear evidence for a molecular outflow but rather for an inflow. However, NGC 1365 has long been known to exhibit a biconical outflow in, for example, ionized gas or X-ray (as reported by e.g., Lena et al. 2016; Venturi et al. 2017; Davies et al. 2020). Recent work by Gao et al. (2021) reports evidence for an outflow in the central 5.4 kpc in molecular (CO (2–1)) and ionized gas.

– NGC 1433: Evidence for a central molecular outflow is reported by Combes et al. (2013, in CO (3–2) emission), and re-analyzed with additional observations by Smajić et al. (2014, in molecular hydrogen and CO (3–2) emission)

– NGC 1672: Combes et al. (2019) study the cold molecular gas torus in the central region and do not find evidence for an outflow. Fazeli et al. (2020) rules out the presence of an outflow in hot molecular and ionized gas.

– NGC 3351: Leaman et al. (2019) find that a molecular bubble exists in the center, and Swartz et al. (2006) interpret their results from X-ray as evidence for an outflow confined by cold ambient gas.

– NGC 4579: This galaxy possesses a jet in the center and radio emission indicates the presence of shocked material at R ≳ 100 pc (Contini 2004). García-Burillo et al. (2009) conclude the presence of an outflow to be unlikely based on their multiphase observations (CO (1–0), CO (2–1), HI and several optical, UV and IR images).

– NGC 4941: Tentative evidence for outflowing photoionized gas is presented by Coccato et al. (2005). Furthermore, high-resolution molecular CO (3–2) gas observations by Garcia-Burillo et al. (2021) find tentative evidence for an outflow.

– NGC 5643: The presence of an outflow is recently confirmed by high-resolution molecular CO (3–2) gas observations by Garcia-Burillo et al. (2021). Also, evidence for a cold molecular outflow is reported by Davies et al. (2014, in S(1)H₂ (1–0)) and in Alonso-Herrero et al. (2018, in CO (2–1)) and the outflow region is analysed in more detail by García-Bernete et al. (2021) using CO and ionized gas emission. Several other studies of different gas phases find evidence for an outflow (e.g., Cresci et al. 2015, for ionized gas).

For the four galaxies that did not fulfill the stricter lwm criterion, and are not classified as secure candidates, we find the following:

– NGC 1097 is known to possess an inflow in molecular gas as well as other gas phase studies (e.g., van de Ven & Fathi 2010; Piñol-Ferrer et al. 2011). There is evidence for inflows as well as outflows between the central spiral arms seen in warm molecular hydrogen gas H₂ (Davies et al. 2009).

– For NGC 1637, NGC 4321, and NGC 4569 no cold molecular outflow studies are available, and only indications for hot gas outflow (Immler et al. 2003) for the former and detection of an ionized outflow (Boselli et al. 2016) in the latter.

– For the remaining galaxies NGC 1317, NGC 2566, NGC 3507, NGC 3627, NGC 4293, NGC 4457, NGC 5134, NGC 7496, no studies on molecular outflows were found in the literature.

PHANGS non-outflow candidate literature matches. PHANGS galaxies that are known to possess molecular outflows in the central 1 kpc region that we did not recognize as outflow candidates:

– For NGC 1566 several studies have analyzed the kinematics of molecular and other gas phases. Slater et al. (2019) interpret their findings in CO (2–1) and ionized gas at a resolution of ∼24 pc as outflowing gas with little evidence for an additional inflow. On the other hand, Combes et al. (2014, 2019) argue for the presence of an inflow in their molecular gas observations rather than an outflow in this galaxy. In this work we obtained an of 1.89 for this galaxy.

– Davies et al. (2014) find molecular outflow signatures (S(1)H₂ (1–0)) in NGC 6300 and NGC 7743. Our s are rather low with 1.56 for NGC 6300 and 0.33 for NGC 7743. For NGC 6300, Gaspar et al. (2019) interpret a broad component in ionized gas as a central outflow, but they find no outflow signatures in molecular hydrogen observations.

– Stone et al. (2016) analysed molecular OH emission of several PHANGS galaxies (NGC 1068, NGC 1365, NGC 1566, NGC 4579, NGC 4941, NGC 6300, NGC 7465). They find no evidence for inflows in these galaxies. For the galaxies where OH absorption was available (NGC 6300), a search for galactic-scale outflows provided no evidence for an outflow.

– García-Burillo et al. (2014) find strong evidence for an molecular outflow in NGC 1068 in CO (2–1). The outflow is orientated along the inner stellar bar of the galaxy, as revealed by a careful modeling of the velocity distribution in their work.

– Yukita et al. (2012) find signs that suggest a potential outflow in X-ray emission for NGC 2903.

5.3.2. Sample comparison

The most comprehensive compilation of central molecular outflows has been assembled by Fluetsch et al. (2019). In Fig. 10, we compare our outflow candidates to their compiled sample in the SFR versus stellar mass plane. The galaxy NGC 1433 is in both the sample of Fluetsch et al. as well as our outflow candidate sample, although SFR and stellar mass differ slightly. On average our candidates overlap with the low-mass end of their sample. In this mass range, their galaxies are also mostly consistent with being on the main sequence of star-forming galaxies while most of their galaxies with stellar mass of log₁₀(M_⋆/M_⊙) ≥ 11.0 are in the starburst regime. This is in line with our study where ∼36% of the galaxies in the stellar mass range 10.0 ≤ log₁₀(M_⋆/M_⊙)≤11.0 are outflow candidates, in fact 19 of our outflow candidates are in this mass regime.

Fig. 10. Outflow candidates (red circles, orange stars) compared to galaxies without outflows (gray circles) in the PHANGS sample and the compilation of outflows by Fluetsch et al. (2019) (blue diamonds). The main sequence of star-forming galaxies from Leroy et al. (2019, grey line) with a width of 0.5 dex (gray shaded area) and extrapolated to higher and lower stellar masses (dotted line) is shown for reference.

5.3.3. Literature outflow frequencies

Assuming that all of our (secure) outflow candidates correspond to actual outflows, the number of (secure) outflow candidates translates into an outflow frequency of 25 ± 2% (20 ± 2%) for the full sample of 80 galaxies used out of the 90 PHANGS galaxies. If we restrict it to the stellar mass range of 9.95 ≤ log₁₀(M_⋆/M_⊙)≤11.15, where we have identified outflows, this increases to 36%. This is the first time the frequency of central molecular outflows in high resolution and high sensitivity data has been robustly estimated as frequency determination in molecular gas are rare, and sample sizes are often small. For example, Domínguez-Fernández et al. (2020) analyzed the CO (2–1) emission in five nearby Seyfert galaxies and find two of their targets being likely to host an outflow.

Wylezalek et al. (2020) used the MaNGA sample of galaxies at mean redshift of z ∼ 0.03 to determine the outflow frequency in galaxies. Over matched stellar mass range of 9 ≥ log₁₀(M_⋆/M_⊙) ≥ 11, they find that 25% of low- to intermediate luminosity AGN and 7% of non-AGN MaNGA galaxies exhibit outflow-typical emission line widths. This is a lower fraction than what we expect based on our sample: as mentioned in Sect. 5.1 we find 53% (42% secure) outflow candidates in the 19 AGN-active PHANGS galaxies and 17% (13% secure) outflow candidates in PHANGS galaxies without AGN. However, gas phase, resolution as well as detection methods are different.

Veilleux et al. (2013) studied molecular gas outflows traced by OH in 37 local ULIRGs and QSOs, and find evidence for outflows for 70% of their targets, which unsurprisingly is higher than for our less active sample.

López-Cobá et al. (2019) analyzed ionized gas emission from 203 highly inclined disk galaxies in the CALIFA sample to identify outflow candidates over their full stellar mass range 9.5 < log₁₀(M_⋆/M_⊙) < 11 and in the SFR range of about −0.7 < log₁₀(SFR/M_⊙yr⁻¹) < 0.7, resulting in an ionized outflow frequency of 8% for the highly inclined galaxies or 2% for the extended CALIFA sample. This is in contrast to our galaxy sample, that selects galaxies of low to moderate inclination (i < 75 deg, see Leroy et al. 2021b). This or the difference in gas phase and tracers may explain the discrepancy in frequency estimates.

Roberts-Borsani et al. (2020) studied neutral gas (traced by Na I D) in 405 local MaNGA galaxies with stellar masses (M_⋆ ≥ 10¹⁰ M_⊙) similar to our sample and find outflow evidence in the central regions of 78 galaxies (19%).

Stone et al. (2016) analyzed emission from OH as tracer for warm molecular gas for galactic-scale outflows in local AGN galaxies (BAT-AGN). They find that 24% of their targets, in which an outflow analysis was possible, show evidence for outflows. One of the PHANGS galaxies (NGC 6300) was included in this sub-sample and no outflow was found, in agreement with our analysis.

A study of ionized gas in main sequence galaxies and galaxies near the main sequence with 9 < log₁₀(M_⋆/M_⊙) < 11.7 at redshift 0.6 < z < 2.7 by Förster Schreiber et al. (2019) reports a global ionized gas outflow frequency of ∼25% (with ∼10% being SF-driven and ∼15% AGN-driven) with a strongly increasing fraction of AGN-driven outflows toward higher-mass galaxies (log₁₀(M_⋆/M_⊙) ≳ 10.8) to up to 3/5. Although this frequency value agrees with our results, a direct comparison is not straight-forward, as the outflow tracers, the galaxies, and their environment are very different. For example, our candidate molecular outflows have much smaller spatial extents and also lower outflow velocities.

In summary, our frequency of central molecular outflows of (20−36)% ± 2% depending on the mass range considered agrees reasonably well with previous studies, though a direct comparison is significantly hindered by the use of different gas tracers, methods, galaxy populations, as well as the data quality. This emphasizes the need for studies of more nearby, well resolved galaxies for a direct comparison to our results. Thus, in the context of galaxy evolution our numbers need to be interpreted carefully, as our frequency allows for two interpretations:

• a)

  One possible scenario would be to assume that all galaxies possess outflows that are detectable with our methods. This immediately implies a short outflow lifetime, so that we end up with our 20−36% frequency.

• b)

  Another scenario is that outflows are long-lived, and thus our frequency tells us that outflows are more rare events.

Likely, the answer will be in-between. Total SFR and stellar mass play an important, but yet not fully understood role in outflow launching and future research might uncover unknown processes to answer these questions.

5.4. Caveats and improvements

Although the methods used in this work are commonly applied in the literature, we note some caveats regarding the construction and evaluation of spectra, pv–diagrams, and lwms, which add to the uncertainty of the outflow frequency estimate.

Construction and analysis of the method figures. In general, projection effects, that make it impossible to detect outflow signatures with our methods, can lead to an underestimation of the outflow frequency.

– Spectra and pv–diagrams: Instead of using fixed values for aperture and slit size, observationally motivated values for each galaxy (e.g., depending on galaxy size) may avoid evaluating misleading signals as outflow signs or enhance fainter outflow features.

– lwms: The automated creation of lwms depends on the spectrum and its signal to noise (S/N). If the spectrum has low S/N, the intersections with the 5% and 20% level of the peak will not contain the intended information, but only random velocity ranges between noise peaks. In case of a noise-dominated spectrum, no significant amounts of cold molecular gas are present in the central galaxy region, thus a detection of a molecular outflows is also unlikely. For sparsely distributed gas in the center, summed up noise can hide faint broad wings in the spectrum. Whereas with an optimized velocity range, these features may become visible in the lwms since the spatial information is retained.

Visual evaluation. OCLs were allocated by visual inspection of the figures (Sects. 3.2.1–3.2.3).

– Visual identification of potential outflow components is generally a less rigorously defined approach than fitting the spectra and investigating pv–diagrams. However, fitting proved to be an unreliable method in tests due to the large variety of spectral shapes and S/N in our sample.

– Some of the pv–diagrams and lwms showed potential low-level continuum emission or other artefacts, that may either overlap and hide actual outflow emission, or could be mistaken for outflow emission. As all inspectors were aware of this effect, the impact should be negligible.

– In addition to the lwms, velocity fields and residual velocity fields could potentially improve the outflow identification process, However, this would result in a much more detailed and time-consuming analysis of each galaxy, which is beyond the scope of this paper.

Overall we find a good agreement between the three authors’ assigned OCLs, which indicates that – despite the presence of some caveats – reproducibility is achieved.

6. Summary and conclusion

We use ∼1″ ≈ 100 pc resolution CO (2–1) imaging for 80 nearby, massive, star-forming galaxies from the PHANGS-ALMA survey (Leroy et al. 2021a, see also Sun et al. 2020) to search for molecular outflows in the central regions of 300 pc to 2 kpc in diameter. Based on our applied methods, we identify 20 outflow candidates out of which 16 are classified as secure outflow candidates which translates into an estimated frequency of central molecular outflows in massive (9.17 ≤ log₁₀(M_⋆/M_⊙)≤11.15) star-forming main sequence galaxies of 25 ± 2% (or 20 ± 2% for secure candidates; see Sect. 4.2).

Further results are as follows:

1. We examine central spectra of two different apertures, pv–diagrams along major and minor axes, and automatically created line-wing maps (lwm) using velocity ranges determined from the spectra to identify kinematic signatures possibly tracing outflowing gas (Sect. 3). Combining these methods ensures high robustness against other kinematic features, such as gas flows along bars. The resulting outflow confidence labels (OCLs) allocated via visual inspection agree well among different inspectors for each galaxy and can be used to determine the outflow candidates as well as secure outflow candidates (Sect. 4, Fig. 4). Global host galaxy properties (stellar mass, SFR) show a significant correlation with these confidence labels, which could be pointing to a correlation of these labels with actual outflow incidence or a correlation of outflow properties (e.g., outflow mass and extent) with stellar mass and/or SFR, which enhance the detectability of outflow signatures (Sect. 5.1.1).

2. No outflow candidates are identified in galaxies with stellar masses below log₁₀(M_⋆/M_⊙) = 9.95 (Sect. 5.1.1 and Fig. 6). The estimated frequency of central molecular outflows in the mass range 9.95 ≤ log₁₀(M_⋆/M_⊙) ≤ 11.15 is 36%. This may indicate a possible stellar mass and/or SFR threshold needed to launch outflows or a threshold for our detection methods.

3. The fraction of galaxies hosting an AGN is ∼2× higher among outflow candidates (50%, secure candidates: 50%) than the fraction in the entire sample (24%). The same applies to barred galaxies which are ∼1.5× more prominent among outflow candidates (89%, secure candidates: 87%) than the entire sample (61%; Sect. 5.1, Fig. 5). This may imply that bars transport the gas necessary to build up an outflow towards the central region. Similarly, AGN might provide the necessary energy to launch outflows. However, this may also be related to the fact that both AGN and stellar bars correlate with global host galaxy properties (Sect. 5.1.2).

4. The molecular outflows of our outflow candidates have a typical extent of ∼260 pc with velocities of ∼130  km s⁻¹ as well as integrated CO luminosities of ∼8.6 K km s⁻¹, indicating that high resolution, sub-kpc, high sensitivity studies such as PHANGS are required to detect molecular outflows in main sequence galaxies (Sect. 5.2).

5. Our outflow candidates extend the trend between SFR and mass outflow rates down to lower outflow rates (log₁₀(Ṁ_out/M_⊙ yr⁻¹) < 0) (Fig. 8). Mass loading factors are on average near unity, thus indicating that these small and less-massive outflows remove as much gas from the centers as is consumed by SFR, and can represent an important correction to the overall future evolution of the gas (Fig. 9).

Our outflow candidates significantly expand the previous compilation of molecular outflow galaxies by Fluetsch et al. (2021), in particular for main sequence galaxies with intermediate stellar mass. Our study suggests that central molecular outflows might be a common feature (∼1 out of 3 galaxies may possess an outflow) in star-forming main sequence galaxies with stellar mass log₁₀(M_⋆/M_⊙)≳10.0. Given that the outflow signatures are not very prominent in many of our galaxies, spatial resolution and sensitivity are important to identify faint outflows and to distinguish kinematic signatures from outflows and other non-circular motions in galaxies.

Acknowledgments

We greatly appreciate the very helpful and constructive comments by the anonymous referee. This work was carried out as part of the PHANGS collaboration. A. T. B. and F. B. would like to acknowledge funding from the European Research Council (ERC) under the European Unio’s Horizon 2020 research and innovation programme (grant agreement No.726384/Empire). E. S., H. A. P., T. S., and T. G. W. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). R. S. K. acknowledges financial support from the German Research Foundation (DFG) via the collaborative research center (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2, and B8). He also thanks for funding from the Heidelberg Cluster of Excellence “STRUCTURES” in the framework of Germany’s Excellence Strategy (grant EXC-2181/1, Project-ID 390900948) and for funding from the European Research Council via the ERC Synergy Grant “ECOGAL” (grant 855130). J. M. D. K. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through an Emmy Noether Research Group (grant number KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1), as well as from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). H. A. P acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). M. Q. acknowledges support from the research project PID2019-106027GA-C44 from the Spanish Ministerio de Ciencia e Innovación. E. R. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2017-03987. The work of J. S. is partially supported by the National Science Foundation under Grants No. 1615105, 1615109, and 1653300. A. U. acknowledges support from the Spanish funding grants AYA2016-79006-P (MINECO/FEDER), PGC2018-094671-B-I00 (MCIU/AEI/FEDER), and PID2019-108765GB-I00 (MICINN). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00650.S, ADS/JAO.ALMA#2013.1.00803.S, ADS/JAO.ALMA#2013.1.01161.S, ADS/JAO.ALMA#2015.1.00121.S, ADS/JAO.ALMA#2015.1.00782.S, ADS/JAO.ALMA#2015.1.00925.S, ADS/JAO.ALMA#2015.1.00956.S, ADS/JAO.ALMA#2016.1.00386.S, ADS/JAO.ALMA#2017.1.00392.S, ADS/JAO.ALMA#2017.1.00766.S, ADS/JAO.ALMA#2017.1.00886.L, ADS/JAO.ALMA#2018.1.01321.S, ADS/JAO.ALMA#2018.1.01651.S. ADS/JAO.ALMA#2018.A.00062.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

References

Appendix A: Outflow frequency threshold dependence

We provide the outflow frequency as a function of its threshold in Figure A.1. A steep decline at small thresholds (∼0.0−0.5), is followed by a linear trend for medium thresholds (∼0.5−2.0) and another steeper decline for large thresholds (≳2.0). In this work, we set the threshold for galaxies to be considered an outflow candidate to as well as for the secure outflow candidacy a threshold of . We can see from this plot, that most of our galaxies (> 50%) have , thus exhibit no significant outflow signatures at all.

Fig. A.1. Frequency of galaxies with threshold value, as a function of threshold in steps of 0.1 between 0.0 and 2.9 for the full sample of 80 galaxies.

Appendix B: OCL method analysis - the presence of AGN and bars

In Section 5.1.2 we found a significant (> 3σ) correlation between and stellar mass as well as SFR. In this Section we search for secondary dependencies due to the presence of stellar bars or AGN.

Bars The upper panels in Figure B.1 show the s as a function of stellar mass (top left) and SFR (top right), divided into sub-samples of barred (blue circles) and unbarred galaxies (cyan diamonds). The Spearman rank correlation coefficients and their uncertainties are provided at the bottom of each panel (see also Table B.1), correlating s with stellar mass and SFR of either barred or unbarred galaxies. Uncertainties are derived by perturbing the and stellar mass (SFR) of all sub-samples by their errors, calculating the Spearman correlation coefficient for the perturbed values, and then repeating this process 1000 times. We infer the scatter as uncertainty.

We find equally high Spearman rank correlation coefficients between both stellar mass and as well as SFR and for barred galaxies compared to those of the full sample (ρ_M⋆ = 0.56 ± 0.03 and ρ_SFR = 0.40 ± 0.03). The unbarred systems have an insignificantly (< 3σ) smaller correlation coefficient between and stellar mass. Also, unbarred systems exhibit a significantly different rank correlation coefficient between and SFR compared to both those of the full sample and the barred sub-sample. This indicates that SFR increases the likelihood of outflow signatures in barred galaxies, but not in unbarred galaxies, where the correlation coefficient is close to zero. However, most of the unbarred galaxies are clustered at lower SFR and low s, thus, if the likelihood of outflows to occur increases with SFR, this correlation is not due to the presence of a bar, but rather due to enhanced bar fraction at higher SFRs. We interpret this in the following way:

• a)

  Kinematic patterns of bars can be mistaken for outflows, since both outflows and streaming motions in bars exhibit substantial radial velocity components. However, by combining different methods we minimize the impact on the s. In particular, the OCLs of lwms tend to show the least trends as described above.

• b)

  The higher values are only the result of higher SFR and the presence of a bar is a coincidence, as bar fraction increases with SFR in our sample as seen in the histograms of Figure B.1. This could imply that SFR itself is the driver for the enhanced outflow probability. With a higher fraction of barred galaxies to unbarred galaxies at higher SFR, the chances are higher to find more outflow candidates in barred galaxies than in unbarred galaxies, resulting in a putative correlation.

• c)

  Bar instabilities can be an efficient means to transport gas into the center of galaxies (e.g., bar-fed nuclear starburst in NGC 253, e.g. Chown et al. 2019). The larger gas reservoirs in the centers of barred galaxies can then lead to a higher fraction of detectable outflows relative to unbarred galaxies as well as a higher SFR. This is similar to the previously mentioned enhanced correlation between central molecular gas density and (compare Section 5.1.2).

We find that this correlation, that barred galaxies show stronger trends between and SFR than galaxies without a bar can be explained by the properties of host galaxy parameters, but we cannot exclude that other mechanisms are also at play. The trends for with stellar mass are not significantly different for barred and unbarred galaxies, indicating that stellar mass might directly correlate with the likelihood of outflows to occur, rather than with the presence of a bar.

We also tested for correlations between stellar mass (SFR) and outflows for the individual identification methods for the sub-samples and found no significant differences (see Table B.1).

AGN Here we test if the depends on the presence or absence of an AGN. The bottom panels of Figure B.1 show of host galaxies with (purple circle) and without (orange diamond) AGN as a function of stellar mass (bottom left panel) and SFR (bottom right panel). Again the Spearman rank correlation coefficients between and stellar mass (SFR) are provided in each panel for the sub-samples. Uncertainties are derived in the same way as done for the barred and unbarred galaxies.

Both active and inactive sub-sample reveal significant (> 3σ) correlations between stellar mass and , which are in agreement within 3σ with each other and with the correlation of stellar mass and found in the full sample. For galaxies with an AGN the correlation coefficient between SFR and is significantly (> 3σ) higher than the one for the inactive galaxies and also significantly (> 3σ) higher than the values of the full sample. This implies that a higher SFR leads to a higher outflow confidence when an AGN is present. Such a trend is not seen for the inactive sub-sample, which shows a slightly weaker trend between and SFR, in agreement with the one of the full sample. The inactive and active galaxy rank correlation coefficients differ by 5σ.

Fig. B.1. dependence on the presence of a stellar bar or AGN. Top row: s of PHANGS galaxies as a function of the host galaxy’s stellar mass (left) or SFR (right) for barred (blue circles) and unbarred galaxies (cyan diamonds) (galaxies without a bar classification: grey cross). Histograms of the barred (blue) and unbarred (cyan) sub-sample are added for stellar mass (top left), SFR (top right panel) and (far right panel). Bottom row: Analogous plot for galaxies hosting an AGN (purple circles) or not (orange diamonds). Histograms of the AGN (purple) and non-AGN (orange) sub-samples are added for stellar mass (top left), SFR (top right) and s (far right panel). Spearman rank correlation coefficient between and galaxy property are listed at the bottom of each panel, for the respective sub-samples (ρbar, ρno bar, ρAGN, ρno AGN). For comparison, the Spearman rank correlation coefficient of the full sample between and stellar mass is ρM⋆ = 0.56 ± 0.03 and between and SFR is ρSFR = 0.40 ± 0.03. Typical uncertainties for and the galaxy property are added in the upper left corner of each panel.

However, the AGN fraction of 24% (19 galaxies) in the PHANGS sample of 80 galaxies is too small to perform proper statistical tests. Figure B.1 shows the spare distribution of galaxies with AGN. The uncertainty of the Spearman correlation coefficient utilizes the statistical uncertainties of stellar mass and , which is likely much smaller than the actual systemic uncertainties. The observed trend is thus tentative and could also be the result of two clusters of galaxies with AGN in the SFR– plane.

On the other hand, AGN activity might actually help boost the correlation between and SFR. This is not unexpected, as there is a relation between SFR and AGN activity due to stellar feedback-driven feeding (e.g., Davies et al. 2007). Also, while stellar mass is distributed over large radii and varies slowly in time, SFR and AGN activity can vary significantly in time and can be more centrally concentrated, which is key for our study of central outflows.

We can interpret this trend in the following ways:

• a)

  Bias in WISE SFR estimates: global SFRs are determined using WISE data. In cases where a luminous AGN is present, this can lead to a significant overestimation of the SFR, which might falsify the observed trend. Among the PHANGS galaxies, such a luminous AGN is only present in NGC 1365, which should not strongly contribute to the observed trend.

• b)

  SFR correlates with AGN strength (compare e.g. Davies et al. 2007): The trend between SFR and OCL is actually the trend that higher SFR implies higher AGN strength, and our OCLs are biased towards higher AGN strength, which could again be a trend with stellar mass, because SFR and stellar mass are clearly correlated in this MS-sample.

• c)

  A recent ionized gas study by Woo et al. (2020) find that AGNs with stronger outflow strengths are hosted by galaxies with higher SFR. If some of the AGN galaxies possess outflows, they might show stronger outflow signatures at higher SFR, and thus we expect higher OCLs at higher SFR. As our final outflow candidate sample is small and only half are hosted by an AGN galaxy, this can not explain the correlation seen.

In short, the presence of an AGN or bar seems to affect the , either because of outflow kinematics, revealing an actual correlation, or due to underlying correlations with other galaxy properties instead of outflows, which we can not exclude.

Similarly as done for the correlation between bars and , we perform tests whether the correlation between AGN and OCL depends on the particular method used to identify and classify the outflows. The resulting Spearman rank correlation coefficients can be found in Table B.1. Except for one value, we do not find a significant difference.

B.1. OCL method analysis – Impact on individual methods

To test if all methods are equally impacted by the presence of bars, we calculate Spearman correlation coefficients for OCLs of the individual methods (Spectra with 300 pc diameter, spectra with annulus with outer 2 kpc diameter, pv–diagrams and lwm, each averaged over the different inspectors) and stellar mass (SFR) for all galaxies (see Table B.1 Columns (2), (3)) as well as barred and unbarred galaxies (Columns (4)–(7)). We do the same for AGN and inactive galaxies (Columns (8)–(11)).

Except for one value we do not find any significant (> 3σ) difference between the Spearman rank correlation coefficients of individual methods compared to the previously presented coefficients calculated with the . Still, we observe that lwm-OCL tend to have the least correlation coefficient among all values. The 2 kpc spectra-OCL on the other hand tend to have the highest correlation coefficient.

The only significant difference can be found for the Spearman rank correlation coefficient between the 2 kpc spectra-OCL and stellar mass at 3σ. One possible explanation might be, that some CO velocity curves of our galaxies peak at ∼1 kpc radius from the center (Lang et al. 2020), thus roughly at the outer edge of the 2 kpc (diameter) aperture. This peak can correspond to a small amount of emission which broadens the spectrum and makes it harder to identify potential outflows. This effect then results in a larger OCL. As more massive galaxies tend to have central disks or bulges, they also have higher rotation velocities in their central 2 kpc. With stellar mass correlating with SFR, this could explain the found deviation. This is, however, a tentative assumption, as many galaxies do not show a clear peak in their rotation curve and the peaks can stretch over a range of galactic radii.

Appendix C: PHANGS galaxy properties

We list the global properties of the 80 PHANGS galaxies used in this analysis in Table C.1.

Appendix D: PHANGS galaxy OCL results

We list all inspector-averaged OCLs per method and per galaxy in Table D.1 a well as all total averages (. Outflow candidates are marked in column (2).

Appendix E: Additional pv-diagrams

We present pv-diagrams as in Figure 2 for all remaining outflow candidates.

Fig. E.1. same as Figure 2 for the remaining outflow candidates. Galaxies marked with * have continuum emission at offset = 0 due to imperfect subtraction in the data pipeline. This issue is fixed in the publicly available data release.

Fig. E.2. continued.

All Tables

All Figures

Fig. 1. Spectra of NGC 1566 and NGC 1546 used for visual classification of outflow candidates. Integrated and normalized spectra (black solid line) of the central regions, i.e., the inner 300 pc in diameter (top panels), and using annuli covering diameter with 300 pc < dgal < 2 kpc (bottom panels). For comparison, the integrated spectrum from the whole galaxy normalized to a peak of unity is shown by a gray line. The red solid lines represent the interpolated 3σ noise level for the central region/annulus. The mid-points of the 20% intersections of the peak (green triangle) are shown for comparison. Velocity has the systemic velocity subtracted. The transparent blue and red areas represent the velocity ranges over which the blue- and red-shifted line-wing maps (Sect. 3.2.3) are integrated over. These areas are selected to correspond to the velocity at 20% of the peak from the central spectrum out to the velocity at 5% of the peak plus an empirically determined additional 25 km s−1 to take into account faint emission at low significance. The green solid lines mark the 20% and 5% peak levels, respectively. Intersections with the central spectrum are given as red triangles. The spectra of NGC 1566 received OCLs of 3, 3, 2 for the 300 pc spectrum and 2, 2, 3 for the 2 kpc-annulus spectrum from all inspectors due to the pronounced line wings. Both images of NGC 1546 received an OCL of 0 from all inspectors for both apertures as no broad wings are evident.

In the text

Fig. 2. Example pv–diagrams along the major and minor kinematic axes of NGC 3627 (top) and NGC 4654 (bottom) used to search for outflows. The color scale represents the intensity at a given position along the offset averaged over a width of 3 pixels. The contours are multiples of the noise at 2.5σ, 3.5σ, followed by niσ with ni = 2ni − 1, for n0 = 2.5σ and i = 1, 2, 3, …. The color range is saturated to emphasize the faint emission we are interested in. The pv–diagrams of NGC 3627 (NGC 4654) correspond to OCLs of 3, 2, 2 (0, 0, 0) from the three different inspectors. The systemic velocity (yellow horizontal line) is added for comparison and is subtracted from the velocity axis.

In the text

Fig. 3. Line-wing maps for NGC 1566. Left panels: full moment-0 map available with the integrated lwms (red- and blue-shifted contours, respectively) overplotted. Right panels: zoom-in of the central region with 2 kpc diameter. The beam size is shown as a gray circle in the lower left corner of the zoom images. Center positions are indicated with a white plus. The lwms shown in the upper panels are derived by integrating over the velocity range selected from the 300 pc spectrum as described in Sect. 3.2.3, the lower panels use the 2 kpc annulus spectrum. A white circle indicates the size over which the corresponding spectrum is taken. These 4 panels together received an OCLlwm of 1, 1, 0 from the three inspectors. The contour levels are selected to be the same as for the pv–diagrams, with the noise being calculated for the zoom images, not for the full maps. The underestimation of the noise at the edges of the full map (left panels) is evident, however, it does not affect our search for central outflows.

In the text

Fig. 4. Average outflow candidate OCLs for PHANGS galaxies. Black crosses represent the per galaxy, which is the average of the inspector-averaged OCL of the pv–diagram, lwm, and the higher ranked spectrum. Galaxies with OCL ≳ 2 are highly likely to exhibit strong central outflow signatures, whereas OCL ≲ 1 correspond to an absence of any outflowing gas signatures. For each galaxy the inspector-averaged OCLs of the individual methods are shown as red circle (central spectrum with 300 pc aperture), purple square (spectrum of annulus with 2 kpc outer diameter), blue triangle (pv–diagram), and light-blue diamond (lwm). The green-shaded area marks all 2 ≤ OCL ≤ 3. A galaxy with (and lwm-OCL ≥ 2) is considered an (secure) outflow candidate and emphasized with green galaxy names (see Sect. 4.2). Uncertainties are derived by perturbing and jackknifing the OCL per method per inspector, as described in the text.

In the text

	Fig. 5. Fraction of AGN (violet, left) and bars (blue, right) among PHANGS outflow candidates and the secure outflow candidates, as well as the entire PHANGS sample.
In the text

Fig. 6. (a) Distribution of PHANGS outflow candidates (secure candidate: red hexagon, otherwise: orange star) and non-outflow galaxies (gray hexagon) along the star-forming main sequence. The histograms show the distribution over stellar mass (top) and SFR (left) for outflow candidates (red, all candidates combined) and galaxies without outflows (gray). The main sequence according to Leroy et al. (2019) is added for comparison (gray solid line) with a width of 0.5 dex (gray shaded area) and extrapolated to larger and smaller stellar masses (dotted gray line). (b) Distribution of PHANGS galaxies along the main sequence of star-forming galaxies. The is shown by color. A dark color represents a large (e.g., class 3), thus showing host galaxies with high probabilities to possess molecular outflows. Light colors represent low (e.g., class 0), thus no or only tentative outflow signatures are present in these galaxies. The histograms show the distribution over stellar mass (top) and SFR (right) for three bins. (c) Ratio of central-to-global SFR as a function of stellar mass for PHANGS galaxies likely to host a central molecular outflow (secure candidate: red hexagon, otherwise: orange star) or not (gray hexagon). Histograms show the distribution of all outflow candidates (red) compared to galaxies without outflow signatures (gray) in stellar mass (top) and SFR ratio (left). AGN are marked by an open black diamond. For all plots the median uncertainty for individual data points is shown in the lower right corner.

In the text

	Fig. 7. Fraction of (secure) outflow candidates as a function of stellar mass. The fraction is determined per stellar mass bin with widths of Δlog10(M⋆/M⊙) = 0.5 that contain a minimum of 3 galaxies. Each level represents the fraction of outflow candidates relative to all galaxies in the mass bin. The total amount of galaxies per mass bin is given at the bottom of each bin.
In the text

Fig. 8. Mass outflow rates Ṁout as function of SFR for the secure outflow candidates (red circle) and possible outflow candidates (orange star) from PHANGS as well as the outflows compiled by Fluetsch et al. (2019) (blue diamond, blue errorbar). The gray line corresponds to a mass loading factor η of unity. PHANGS SFRs have an uncertainty of 0.11 dex and the mass outflow rates have statistical uncertainties of 0.2 dex when perturbing the radius the outflow mass is determined over, with a larger total uncertainty estimated to be ∼1 dex given the unknown actual outflow geometries and properties.

In the text

Fig. 9. Mass loading factor η as function of stellar mass for PHANGS outflow candidates (secure candidate: red circle, otherwise: orange star) and outflows compiled by Fluetsch et al. (2019) (blue diamond, blue errorbar). The PHANGS outflow candidates populate the lower stellar mass and lower mass loading part of the probed parameter space, but overall no significant difference between both samples is found. PHANGS stellar masses have an uncertainty of 0.11 dex and the statistical uncertainty of the mass loading factor due to error propagation is estimated to be 0.23 dex (red errorbar) with a larger total uncertainty of up to ∼1 dex.

In the text

Fig. 10. Outflow candidates (red circles, orange stars) compared to galaxies without outflows (gray circles) in the PHANGS sample and the compilation of outflows by Fluetsch et al. (2019) (blue diamonds). The main sequence of star-forming galaxies from Leroy et al. (2019, grey line) with a width of 0.5 dex (gray shaded area) and extrapolated to higher and lower stellar masses (dotted line) is shown for reference.

In the text

	Fig. A.1. Frequency of galaxies with threshold value, as a function of threshold in steps of 0.1 between 0.0 and 2.9 for the full sample of 80 galaxies.
In the text

Fig. B.1. dependence on the presence of a stellar bar or AGN. Top row: s of PHANGS galaxies as a function of the host galaxy’s stellar mass (left) or SFR (right) for barred (blue circles) and unbarred galaxies (cyan diamonds) (galaxies without a bar classification: grey cross). Histograms of the barred (blue) and unbarred (cyan) sub-sample are added for stellar mass (top left), SFR (top right panel) and (far right panel). Bottom row: Analogous plot for galaxies hosting an AGN (purple circles) or not (orange diamonds). Histograms of the AGN (purple) and non-AGN (orange) sub-samples are added for stellar mass (top left), SFR (top right) and s (far right panel). Spearman rank correlation coefficient between and galaxy property are listed at the bottom of each panel, for the respective sub-samples (ρbar, ρno bar, ρAGN, ρno AGN). For comparison, the Spearman rank correlation coefficient of the full sample between and stellar mass is ρM⋆ = 0.56 ± 0.03 and between and SFR is ρSFR = 0.40 ± 0.03. Typical uncertainties for and the galaxy property are added in the upper left corner of each panel.

In the text

	Fig. E.1. same as Figure 2 for the remaining outflow candidates. Galaxies marked with * have continuum emission at offset = 0 due to imperfect subtraction in the data pipeline. This issue is fixed in the publicly available data release.
In the text

	Fig. E.2. continued.
In the text