The Molecular Outflow in NGC 253 at a Resolution of Two Parsecs

The data presented in this paper are based on observations in ALMA cycles 2, 3, and 4 in bands 3, 6, and 7 that cover the redshifted emission in NGC 253 of CO(1–0), CO(2–1), and CO(3–2), as well as other molecular lines. For the data reduction and imaging of band 3 and 6 data, see Bolatto et al. (2013), Leroy et al. (2015a), Meier et al. (2015), and Zschaechner et al. (2018). Table 1 gives an overview of the data sets used in this analysis.

Table 1.  Details of the Data Sets Used in This Analysis

	CO(1–0)	CO(2–1)	CO(3–2)
ALMA ID	2011.1.00172.S	2012.1.00108.S	2015.1.00274.S
Spatial Resolution	185 × 132	170 × 102	017 × 013
	31.4 pc × 22.4 pc	28.8 pc × 17.3 pc	2.9 pc × 2.2 pc
Spectral Resolution	5.0 km s−1	5.0 km s−1	2.5 km s−1
rms Noise Per Channel	1.99 mJy beam−1	2.19 mJy beam−1	0.81 mJy beam−1
	75 mK	29 mK	0.37 K

Download table as:  ASCIITypeset image

For the band 7 observations, we tuned the lower sideband to 342.0–345.8 GHz and the upper sideband to 353.9–357.7 GHz (total bandwidth 7.6 GHz) with a 976.6 kHz channel width (corresponding to 0.8 km s⁻¹). We targeted the central ∼750 pc of NGC 253 in a linear four-pointing mosaic with two configurations of the 12 m array (12 m compact and 12 m extended, half-power beam width ∼30'') and a five-pointing mosaic of the 7 m array (ACA, half-power beam width ∼50''). Additional single-dish observations with the total power array (TP) recovers emission on large spatial scales. The baseline ranges covered by this setup are 8.9–49.0 m, 15.1–783.5 m, and 15.1–1813.1 m for the ACA and the two 12 m setups, respectively.

The observations were carried out primarily in the first half of 2016 (TP: 2015 December 7 to 2016 August 2; ACA: 2015 December 7 to 2016 November 23; 12 m compact configuration: 2016 April 16, 2016 April 23, 2016 June 17, 2016 June 27; 12 m extended configuration: 2016 August 30, 2016 September 3). The total on-source observation time is 48 hr and 45 minutes split across 27 hr and 23 minutes (TP), 14 hr and 57minutes (ACA), 2 hr and 37 minutes (12 m compact), and 3hr and 59 minutes (12 m extended). The calibrators were J0006–0623 (bandpass); J0038–2459 (complex gain); the asteroid Pallas (absolute flux density); and J0104–2416 and J0106–2718 (both WVR). Visibilities of the 12 m data are calibrated using the ALMA cycle 3 pipeline in casa 4.6.0 and the delivered calibration script. The other data sets are calibrated in casa 4.7.2 and the cycle 4 pipeline.

In order to image the spectral lines, we subtract the continuum in the (U, V) plane using a first-order polynomial fitted to the channels that do not contain strong spectral lines. We reliably detect >25 lines in the range 342.0–345.0 GHz and 353.9–357.7 GHz aside from the four strong lines CO(3–2), HCN4–3), HCO⁺(4–3), and CS(7–6). Most of these lines are weak and only detected in small spatial regions so they do not affect the overall continuum fit and subtraction.

Combined imaging of the interferometric data is done with the tclean task in casa 5.4.0, which includes crucial bug fixes for ALMA mosaics.¹² We regrid the visibilities during deconvolution to a spectral resolution of 2.5 km s⁻¹. Applying a Briggs weighting scheme with robust parameter 0.5 results in a synthesized beam of 017 × 013 (pixel scale 005). The images are cleaned to a level of 2.5× the rms noise in line-free channels of 2.5 × 0.81 mJy beam⁻¹ (2.5 × 0.37 K) using a clean mask derived from a low-resolution image of the compact 12 m array CO(3–2) data only.

We correct the cleaned images for the mosaic sensitivity pattern (mosaic primary beam response pattern), combine them with the TP images using feather, and finally convert the units to brightness temperature.

For the final images, we do not consider the ACA data as they introduce large-scale noise fluctuations toward the edge of the mosaic, which we attribute to decreasing sensitivity of the 12 m data relative to the ACA data. These fluctuations obscure the regions where outflows were previously found. This work requires accurate integrated flux measurements and correct representation of the small-scale structure, which are defined by single-dish observations (TP) and long baselines (extended 12 m), respectively. By checking the images without ACA data against the images including ACA data, we can confirm that neither the overall flux scale nor the small-scale structure is significantly altered.

Data products for CO(1–0) are shown in Bolatto et al. (2013), Meier et al. (2015), and Leroy et al. (2015a), and Zschaechner et al. (2018) present the CO(2–1) data. Imaging results for CO(3–2) are presented in the following section.

In order to keep the amount of detail and contrast in the high-resolution data, we do not match the spatial resolution to that of the data with the lowest resolution, but perform our analysis at the native resolution of each data set. All further steps work on the data cubes masked at 5.0σ (see Table 1) and further masks where necessary. To generate the masks, we do not consider the nonuniform noise level caused by the mosaic sensitivity pattern but use the per-channel rms noise in the center of the field of view.

In this section, we present the CO(3–2) data in different representations. Channel maps (Figure 1), moment maps (Figure 2), and a pV diagram (Figure 3) show the spatial and kinematic structures to be discussed and highlight the data quality.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows channel maps of the image cube. To retain the intrinsic resolution, only every 16th channel (40 km s⁻¹ spacing) is shown here. Aside from the rotating disk of molecular gas, we clearly detect the prominent southwest (SW) streamer (Walter et al. 2017) in the range 180–250 km s⁻¹ (Figure 1, panels 220 and 260 km s⁻¹). Additional gas streamers are apparent between ∼60 and ∼350 km s⁻¹ toward the north and south of the disk, as can be seen, for example, in the panels at 260 or 340 km s⁻¹. Several notable molecular shells are present between 180 and 340 km s⁻¹. Aside from the (super)shells at the eastern (left) and western (right) edges of the map that were previously identified by Sakamoto et al. (2006) and Bolatto et al. (2013), further smaller shell-like structures are located along the molecular disk.

We calculate image moments (Figure 2) with immoments in casa for emission above 5σ for the moment 0 (integrated intensity), moment 1 (intensity-weighted line-of-sight velocity), and moment 2 (intensity-weighted velocity dispersion) maps. Note that due to the complex line shapes, the moment 2 map does not directly correspond to velocity dispersion, which is only the case for Gaussian line profiles. The maps are further constrained to the region defined by the collapsed clean mask to limit them to emission that has been processed by the clean algorithm.

Figure 3 shows the kinematic structure of NGC 253 in a pV diagram along the major axis of the disk (PA = 55°) averaged over the full width of the field of view (∼30'') centered on the kinematic center. The pV cut shows several high-velocity dispersion structures extending from a rotating disk, indicative of outflows.

Our goal is to account for all molecular wind, separating outflowing molecular gas from foreground or background disk emission. A clean separation in 2D position–position space cannot be easily accomplished, due to the inclination of 78° of NGC 253. At this high inclination, outflows and disk emission are cospatial in projection. Kinematic information from line-of-sight velocities, however, makes it possible to disentangle the outflow. Note that this becomes increasingly difficult as the velocity vector aligns with the plane of the sky, resulting in line-of-sight velocities that are systemic. From Hα kinematic modeling, the NGC 253 outflow is approximately biconical with an axis normal to the disk and an opening angle of ∼60° (Westmoquette et al. 2011), and thus the range of possible projection angles is large (see Meier et al. 2015 for a sketch). Note that because the cone-opening angle is larger than the angle between the axis of the cone and the plane of the sky, gas in the approaching and receding cones can have both blue- and redshifted velocities with respect to systemic.

The launching of molecular gas occurs within the disk through star formation feedback, thus the outflows originate from the same location in ppV space as disk molecular clouds. Outflows will therefore blend into the disk near their launching sites, which makes disentanglement increasingly difficult closer to the starburst region.

The complexity of systematically separating the emission corresponding to the disk and the outflow in ppV space is challenging. Algorithmically, this separation is simpler in lower dimensional space, obtained by slicing the data cube into a collection of 2D pV diagrams. In what follows, we identify kinematic components in these diagrams, which we then project back to 3D ppV space. In order to avoid introducing biases, we model the large-scale disk velocity field and use this model as the basis of the kinematic separation.

Images of the center of NGC 253 on large scales show an elongated gas structure (Figure 2 top) with a regular velocity field (Figure 2 middle) that roughly matches a rotating disk disturbed by streaming motions from a bar (e.g., Paglione et al. 2004). The elongated gas structure is consistent with a highly inclined disk of molecular gas, or possibly a ring-like structure as observed in other galaxies (for example, NGC 1512, NGC 1808; Salak et al. 2016; Ma et al. 2018). Similar structures break up at higher spatial resolution into two embracing spiral arms or complex nonclosed orbits in the Milky Way center (Krumholz & Kruijssen 2015; Henshaw et al. 2016; Sormani et al. 2018).

Superimposed on this large-scale structure are smaller features that are not part of the large-scale pattern of rotation and streaming motions. Some of them have high aspect ratios in channel maps and line-of-sight velocity gradients, both typical for outflows. Local deviations from the large-scale velocity field can also be due to infalling gas, or clumps of gas that do not follow the global pattern, due perhaps to cloud–cloud collision.

Henceforth, we will refer to the bulk of the molecular gas that moves according to the large-scale velocity field in the central regions of NGC 253 as the disk. We assume that this large-scale velocity field consists of rotation and streaming motions. The term nondisk refers to any gas that is not following the ppV structure of the disk. By this definition, nondisk gas encompasses material from features that may be attributed to a variety of physical processes, including outflow and infall.

Structures of outflowing gas are frequently referred to by names that describe their kinematic or spatial appearance, such as "streamer." We will use the term outflow to denote localized structures with morphology and kinematics consistent with gas moving away from the disk, as inferred from their location in ppV space. Typical signatures are a velocity that is inconsistent with rotation in the plane of the disk and a high aspect ratio oriented roughly perpendicular to the disk major axis. Note that similar kinematic and structural properties can arise in infalling gas clouds. We will assume that all molecular gas with these characteristics is outflowing, which is likely the case for the majority of the material in NGC 253.

Kinematic analyses typically depend on high signal-to-noise ratios (S/N) because faint features can easily drown in noisy spectra. As a trade-off between necessarily high S/N and also trying to include as much faint emission as possible, we conduct the following analysis on data cubes masked at the 5σ level (see Table 1).

We split the ppV cubes into pV slices along the major axis of NGC 253 as shown in Figure 4. The slices assume the kinematic center is α, δ = 00^h47^m33134, −25°17^m1968 (Müller-Sánchez et al. 2010), and are oriented along the major axis of the projected CO emission with PA = 55°. The area sliced is chosen to cover the region for which we have overlapping CO(1–0), (2–1), and (3–2), and also cover the full length of the SW streamer outflow feature (175; Walter et al. 2017).

Zoom In Zoom Out Reset image size

These requirements are fulfilled by slices of 50'' (850 pc) length (major axis) and covering 50'' (850 pc) along the minor axis (Figure 4). To reduce the problems introduced by splitting features across slices, each slice is 50 (85 pc) wide, and we overlap slices by half their width (25, 42 pc). A sample pV diagram is shown in Figure 5 for the central slice, which runs along the major axis (offset 00). A complete set of pV diagrams is given in Appendix C. The resolution differences between our three transitions, a factor of ∼100 in beam solid angle, are apparent in Figure 5. In the high angular resolution CO(3–2), small features with large line widths are common. These features are blurred out in the lower resolution CO(1–0) and (2–1).

Zoom In Zoom Out Reset image size

We derive a model for the velocity of the disk component from the CO(1–0) observations using the kinematic fitting tool diskfit (Spekkens & Sellwood 2007; Sellwood & Sánchez 2010; Sellwood & Spekkens 2015). Because the CO(1–0) observations cover the largest area among our observations, we use them to derive the model; the additional information provided by the CO(2–1) and/or CO(3–2) data is negligible in terms of the bulk motions of the gas. We obtain a CO(1–0) velocity field by computing the first moment of the cube after masking it at 20σ (1.26 K), in order to represent the velocity of the bright emission. We show the details of the fit parameters and a comparison to the CO(1–0) velocity field in Appendix A.

In each pV slice, we use the velocity profile of the diskfit model to define the local disk velocity. We consider the CO emission to be consistent with the disk component of the emission when the velocity difference is within the local observed velocity range, Δv. This velocity range varies spatially and depends on the distance x from the major axis, increasing toward the center due to the combined effects of higher intrinsic velocity dispersion and projection. For the success of this analysis, it is crucial that Δv is broad enough to cover the observed velocity range of the disk but also narrow enough in order to not classify potential outflows as a disk. The definition of Δv is thus a crucial source of uncertainty for the derived quantities. We parameterize the velocity range of the disk as

$$\begin{eqnarray}&&{\rm{\Delta }}v(x)=120\,\exp \left(-{\left(\displaystyle \frac{x}{2.5}\right)}^{2}\right)+100,\end{eqnarray} \tag{ 1 }$$

with Δv in km s⁻¹ and x in arcseconds. We find this empirical relation to fit the pV data best, and small variations of order 10%–20% already show a noticeable mismatch as discussed in Appendix B. The quality of this definition can be assessed from Figure 5 and Appendix C: the velocity ranges are wide enough to include obvious emission of the disk but do not extend into the kinematically distinct features (potential outflows) that appear as spikes. This is most apparent for CO(3–2) as this line offers the highest spatial resolution. The effect of a 10% change in the velocity range Δv corresponds to up to 0.1 dex variations in the derived quantities (see Appendix B).

We use the modeled velocity field and the Δv relation together to define a "disk mask" over ppV space, corresponding to emission that is consistent with disk rotation. We show in Figure 5 the central pV slices for CO(1–0), (2–1), and (3–2). We show in Appendix C the complete set of pV slices.

Note that the CO emission extends beyond the disk mask. These extensions are not symmetric and due to nondisk gas and projection effects. At ∼78° inclination, gas flowing perpendicular to the galaxy disk toward the south (negative slice offsets) is primarily approaching us and seen at lower velocities relative to the disk emission. Similarly, outflow emission toward the north is primarily at velocities higher than the disk emission. Consequently, emission in pV slices shifts from lower to higher velocities relative to the disk model when the offset from the major axis increases (see Figure 12 in Appendix C). We designed the disk mask to be wide enough to capture the disk emission but exclude the asymmetric component caused by outflows.

We define a "nondisk" mask that is the mathematical complement of the disk mask, with the addition of removing emission from known sources (portions of the spiral arms) that were not included in our model of the central disk and are not of interest for this analysis.

We identify three different types of structures in the nondisk component (Figure 6; contours in Figure 7 highlight these structures):

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

(1) Emission that is colocated with the central disk and bar in projection. This is visible as a ridge in CO(3–2) (inner contour in Figure 7) and also present but less apparent in CO(1–0) and CO(2–1). The structure is unlikely to be an outflow. It appears more likely to be an additional kinematic component of the disk/bar that is not included in the model we used for the separation. We therefore do not consider this gas to contribute to the total mass outflow rate.

(2) Emission associated with the so-called western superbubble, located to the west and north of the central starburst region (Sakamoto et al. 2006; Bolatto et al. 2013; shown by the western contour in Figure 7). This feature is already known to be kinematically distinct from the surrounding gas. Part of it is likely the base of the northern outflow cone (and giving rise to the NW streamers identified by Bolatto et al., for example), but it is difficult to know what portion of the emission should be associated with a net outflow. In our calculations below, we exclude this feature from the total outflow rate of NGC 253, although it likely has some contribution to the outflow.

(3) The remaining gas associated with the nondisk component is organized in small clumps along the edge of the disk region or beyond it. Some of this gas is not discernible as individual structures, particularly in the CO(1–0) and CO(2–1) cubes, perhaps due to the resolution but maybe also due to the excitation conditions, constituting extended regions with diffuse emission. Some of the emission is located in well-defined structures known to be part of the outflow, such as the SW streamer, which is apparent in all CO transitions.

In summary, the nondisk component consists of these three subcomponents: structures that we associate with a net "outflow," structures that are part of the "western superbubble," and structures that are colocated with the "disk." The latter is not associated with the outflow, while parts of the western superbubble may contribute to it. Below, we calculate properties for the two components, disk and nondisk, and its subcomponents individually where it is feasible to do so.

We quantify in Table 2 the CO luminosities of the disk and nondisk components. We measure luminosities of 2.8 × 10⁸ K km s⁻¹ pc², 2.3 × 10⁸ K km s⁻¹ pc², and 1.8 ×10⁸ K km s⁻¹ pc² for CO(1–0), (2–1), and (3–2), respectively, in the central disk of NGC 253. The nondisk component is, naturally, much fainter with luminosities of ∼4.2 × 10⁷ K km s⁻¹ pc² for CO(1–0), ∼4.2 × 10⁷ K km s⁻¹ pc² for (2–1), and ∼4.2 ×10⁷ K km s⁻¹ pc² (3–2). These correspond to approximately 12.9%, 16.4%, and 6.5% of the total luminosity. These luminosities are measured over the sliced area (see Figure 4) for which the coverage is not the same among the data sets. We therefore also measure luminosities integrated over the same spatial region, here defined as the overlap between the data sets. This overlap amounts to 885 arcsec² (2.55 × 10⁵ pc²). The luminosities in the overlap area are disk: 2.6 × 10⁸ K km s⁻¹ pc², 2.1 ×10⁸ K km s⁻¹ pc², and 1.8 × 10⁸ K km s⁻¹ pc²; nondisk: 2.6 ×10⁷ K km s⁻¹ pc², 2.4 × 10⁷ K km s⁻¹ pc², and 1.2 × 10⁷ K km s⁻¹ pc² for CO(1–0), (2–1), and (3–2), respectively.

Table 2.  Results of Separating Disk from Nondisk Emission in NGC 253

Quantity	Units	CO(1–0)		CO(2–1)		CO(3–2)
		Disk	Nondisk	Disk	Nondisk	Disk	Nondisk
Luminosity
LCO	K km s−1 pc2	2.8 × 108	4.2 × 107	2.3 × 108	4.5 × 107	1.8 × 108	1.2 × 107
Fraction	%	87	13	84	16	94	6.5
Molecular Gas Mass Mmola
Totalb	M⊙	3.1 × 108	4.5 × 107	3.1 × 108	6.1 × 107	2.9 × 108	2.0 × 107
Outflowc	M⊙		2.7 × 107		4.1 × 107		8.3 × 106
Superbubbled	M⊙		8.9 × 106		8.9 × 106		5.8 × 106
Other-diske	M⊙		7.6 × 106		7.4 × 106		5.8 × 106
Molecular Mass Outflow Rate $\dot{M}$ f
Outflow (continuous)g	M⊙ yr−1		14		20		2.7
Outflow (constant)h	M⊙ yr−1		29		39		4.8
Kinetic Energy Ekini
Outflow (continuous)g	erg s−1		3.9 × 1054		4.5 × 1054		6.5 × 1053
Outflow (constant)h	erg s−1		2.5 × 1054		3.1 × 1054		4.3 × 1053
Momentum Pj
Outflow (continuous)g	M⊙ km s−1		6.9 × 108		8.7 × 108		1.2 × 108
Outflow (constant)h	M⊙ km s−1		4.8 × 108		6.4 × 108		8.0 × 107

Notes. Uncertainties for these quantities are discussed in the corresponding subsections of Section 4. Sources of error are discussed and quantified in the respective subsections of Section 4.

^aMolecular gas mass derived using a conversion factor for CO(1–0) emission of ${X}_{\mathrm{CO}}=0.5\times {10}^{20}{\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}$, including the contribution of helium, and assuming CO brightness temperature line ratios of r₂₁ = 0.80 and r₃₁ = 0.67 for CO(2–1) and CO(3–2) relative to CO(1–0). ^bCO line luminosity of all emission considered consistent with disk rotation (disk) and not consistent with disk rotation (nondisk), respectively. ^cNondisk excluding the western superbubble and the gas that is cospatial with the projected disk. ^dNondisk emission belonging to the western superbubble as defined by Sakamoto et al. (2006). ^eNondisk gas that is cospatial with the disk in projection. See Section 4.3 for the definition. ^fDeprojected molecular mass outflow rate. Fiftieth percentile best estimate assuming a flat distribution of outflow inclinations for the unknown geometry. ^gOutflowing gas as defined by note (c) under the assumption of continuous mass ejection without accelerations to the gas after ejection. ^hOutflowing gas as defined by note (c) under the assumption of approximately constant starting mass outflow rate over the lifetime of the starburst. ⁱDeprojected kinetic energy of the molecular gas. Fiftieth percentile best estimate assuming a flat distribution of outflow inclinations for the unknown geometry. ^jDeprojected momentum of the molecular gas. 50th percentile best estimate assuming a flat distribution of outflow inclinations for the unknown geometry.

Download table as:  ASCIITypeset image

An interesting result coming out of our decomposition is that not all of the material we identify as "outflow" is in well-defined structures such as the streamers identified by Bolatto et al. (2013). Correctly estimating the outflow rate requires also accounting for a diffuse extended component.

It is important to compare our fluxes to measurements in the literature. Mauersberger et al. (1996) found a CO(2–1) luminosity of 1.2 × 10⁶ K km s⁻¹ arcsec², which translates¹³ to 3.5 × 10⁸ K km s⁻¹ pc² or 1.3 times our measurement. Their observations cover 80'' × 60'', an area similar to our CO(1–0) observations (but ∼4 times larger than the area of our CO(3–2) observations). For the outflow CO(1–0) luminosity, Bolatto et al. (2013) derived an estimate of 2.0 × 10⁷ K km s⁻¹ pc² by summing over individual identified molecular outflow features. This includes flux from the "superbubble" component, so it is probably better compared to the sum of our "outflow" and "superbubble" components of ∼3.6 × 10⁷ K km s⁻¹ pc². Given the large methodological differences and the importance of the diffuse emission, these numbers are in reasonable agreement.

The total gas mass M is estimated from the CO line luminosity, using the conversion factor ${X}_{\mathrm{CO}}=0.5\times {10}^{20}\ {\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}$ cm⁻² corresponding to ${\alpha }_{\mathrm{CO}}=1.1\,{M}_{\odot }{\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{-1}\right)}^{-1}$ discussed by Leroy et al. (2015a) for the central starburst region. This value accounts for the effects of moderate optical depth, high-velocity dispersion, and warm gas temperatures that are likely to dominate the central regions of NGC 253. The masses we report include the contribution of helium to the total mass. To compute masses using the CO(2–1) and CO(3–2) transitions, we assume typical line ratios of r₂₁ = 0.80 and r₃₁ = 0.67 relative to CO(1–0) as implied by Zschaechner et al. (2018). Note that we do not measure line ratios from the images but adopt a uniform factor to keep the mass measurements from the three observed CO lines independent.

Table 2 lists the masses corresponding to the disk and the nondisk components. Uncertainty in the mass estimates arises primarily from the assumed conversion factor and the apportioning of emission among the different components. The calibration uncertainty for the flux measurements is ∼10%–15% for the ALMA observations. Overall, we adopt a systematic error of factor ∼2 for the the derived masses.

The molecular masses derived from the three CO transitions are very similar. They match within 10% for the disk component, and within 50% for the nondisk component. We estimate the total gas mass in the center of NGC 253 to be ∼3.5 × 10⁸ M_⊙ (adding the disk and nondisk components), with estimates in the range of 3.1–3.6 × 10⁸ M_⊙ for the different transitions. About 85% of the total mass is in the disk component. The masses estimated in the nondisk components using CO(1–0) and CO(2–1) are fairly similar at 4.5 × 10⁷ M_⊙ and 6.1 × 10⁷ M_⊙, whereas in CO(3–2) we detect a lower mass, 2.0 × 10⁷ M_⊙, a consequence of the lower luminosity. The nondisk masses are primarily contributed by the outflow component (∼50%). About 20%–30% of mass is in the western superbubble and 12%–30% is cospatial with the disk but kinematically distinct.

It is important to compare these mass estimates to previous results for the total molecular gas mass in NGC 253, noting that our analysis covers the central 45'' × 25'' (750 pc × 400 pc). Toward the east, ∼10% of the known molecular gas close to the center is not covered by our CO(3–2) observations and thus not considered in this analysis. The agreement with previous measurements is very good. Mauersberger et al. (1996) reported a mass of 1.3 × 10⁸ M_⊙ over a similar area (80'' × 60'' in the center of NGC 253), but this was based on a different distance and X_CO. After correcting for those differences, their luminosity corresponds to 4.2 × 10⁸ M_⊙, consistent with our measurement. Using the same distance and the same 1–0 observations, Leroy et al. (2015a) measured a molecular mass of 3.5 × 10⁸ M_⊙. Pérez-Beaupuits et al. (2018) reported a total gas mass of 4.5 × 10⁸ M_⊙ derived from the submillimeter dust spectral energy distribution, which is very consistent with our result given the very different methodologies.

No estimates in the literature separate the "disk" and "nondisk" components as we do above. Previous estimates of the outflowing mass range from a lower limit of 6.6 × 10⁶ M_⊙ calculated for the optically thin limit (Bolatto et al. 2013), to 2–4 × 10⁷ M_⊙, when accounting for optical depth (Zschaechner et al. 2018). Because we identify an outflowing mass of ∼5 × 10⁷ M_⊙, the agreement with the latter estimate is fairly good. Note, however, that these studies derive the outflowing mass from individual features rather than using the pV information as we do here in a systematic way.

Mass flow rate is defined as the flux of mass per unit time through a surface. In our case, we are interested in the flow of molecular gas mass through a virtual closed surface around the center of NGC 253 at a given distance. Note that an individual outflow feature observed over a certain length, such as the SW streamer, can develop in at least two ways: as the distance of an outflow from its origin corresponds to time since ejection times velocity, continuously outflowing gas results in extended streaming structures. Gas ejected at a single ejection event in the past with a distribution of ejection velocities, on the other hand, will also result in an extended streamer. In reality, gas will be ejected with a distribution of velocities at a varying rate over a period of time, and in order to interpret the measurements, we need to make some simplifying assumptions. We chose two edge cases to span the range of different interpretations: (1) the gas does not experience accelerations after being launched (however, see Walter et al. 2017), and (2) that the gas outflow rate is approximately constant with time. For all calculations, however, we assume that the projected direction of flow is perpendicular to the central plane of the bar and that CO emission traces the mass with a constant conversion factor.

We compute the outflow rate both as a function of distance and as a function of time. If we assume that the mass outflow rate has been approximately constant over the lifetime of the starburst, for example, a diminishing outflow rate as a function of distance would suggest that gas is either launched with or somehow develops a distribution of velocities. Conversely, if we assume that the present-day velocity has been constant since the gas was ejected, we can derive a history of the mass outflow rate as a function of time and account for a variable mass outflow rate. Both interpretations are equally, but not simultaneously, valid.

For the detailed calculation of the mass outflow rate, we proceed as follows. The mass outflow rate is $\dot{M}=M\,{t}^{-1}$ with mass M and relevant timescale t. For each image element i (3D pixel or sometimes also called voxel), we calculate the outflow rate ${\dot{m}}_{i}$ of the gas that was ejected at time ${t}_{\mathrm{eject},i}$ over the time interval ${\rm{\Delta }}{t}_{\mathrm{cross},i}$, as the ratio of mass of the pixel m_i to the pixel-crossing interval t_cross,i. The ejection time and pixel-crossing interval are functions of the outflow velocity v_i and the distance s_i between the current pixel position and the launching site, and the pixel size in the direction of the flow Δs, respectively. We therefore compute

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{cross},i}=\displaystyle \frac{{\rm{\Delta }}s}{{v}_{i}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\dot{m}}_{i}\ ({t}_{\mathrm{eject},i})=\displaystyle \frac{{m}_{i}}{{\rm{\Delta }}{t}_{\mathrm{cross},i}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{eject},i}=\displaystyle \frac{{s}_{i}}{{v}_{i}},\end{eqnarray} \tag{ 4 }$$

obtaining a mass outflow rate ${\dot{m}}_{i}$, a distance s_i, and an ejection time t_eject,i for each pixel in the "outflow" component. Note that this approach takes the 3D phase space information into account by treating pixels independently. Typically, a sight line shows multiple pixels with emission at different velocities that all contribute an outflow rate with their respective mass, distance, and velocity. We then bin the outflow rates ${\dot{m}}_{i}$ by ejection time ${t}_{\mathrm{eject},i}$ and integrate over the time range $\left[{T}_{1},{T}_{2}\right]$ to obtain the average outflow rate in this time interval,

$$\begin{eqnarray}&&\dot{M}\left({T}_{1},{T}_{2}\right)=\displaystyle \frac{{\displaystyle \sum }_{i}{\dot{m}}_{i}\left({T}_{1}\lt {t}_{\mathrm{eject},i}\lt {T}_{2}\right){\rm{\Delta }}{t}_{\mathrm{cross},i}}{{T}_{2}-{T}_{1}}.\end{eqnarray} \tag{ 5 }$$

Similarly, binning by distance results in the average outflow rate at a given distance,

$$\begin{eqnarray}&&\dot{M}\left({D}_{1},{D}_{2}\right)=\displaystyle \frac{{\displaystyle \sum }_{i}{\dot{m}}_{i}\left({D}_{1}\lt {s}_{i}\lt {D}_{2}\right){\rm{\Delta }}s}{{D}_{2}-{D}_{1}}.\end{eqnarray} \tag{ 6 }$$

Performing binning on a sequence of time intervals yields the outflow rate history, while binning in distance tells us how far from the launching site a given fraction of the mass is able to escape.

Calculating velocity v and distance s requires knowledge about the geometry and origin of each outflowing gas parcel. The simplest assumption, used here, is that on average, outflows are launched in the plane of the central region of the galaxy, which corresponds to launching on the major axis. The distance s is thus the projected distance to the major axis on the edge of an outflow cone with a given opening angle. Note that the outflow originates from an extended region in the disk and the term cone thus refers to a cutoff cone (called a frustum in geometry). The velocity v is the velocity difference between the launching site and current velocity of the outflow parcel, i.e., the velocity difference over the distance s. The velocities of both the launching site and the projected distance are uncertain. The velocity changes by ±25 km s⁻¹ when an outflow originates from the northern/southern edge of the observed CO disk (above/below the plane), while the projected distance traveled by the gas changes by ±125 (±20 pc).

Distance s, ejection timescale t_eject, and pixel-crossing timescale t_cross are measured as projected quantities that need to be deprojected to account for the outflow geometry. The bright molecular streamers (the SW and SE streamers) seem to lie at the edge of the ionized outflow cone with a ∼60° opening angle (Bolatto et al. 2013). Assuming that all molecular outflows are along this cone and that the axis of the cone is oriented perpendicular to the disk (i = 78°), the range of the effective inclination of outflowing gas can be anywhere between θ = 48° and θ = 108°. Deprojected velocity, ${v}_{\mathrm{depro}}\,={v}_{\mathrm{obs}}/\sin \theta$, and distance, ${s}_{\mathrm{depro}}={s}_{\mathrm{obs}}/\cos \theta$, have a direct effect on the deprojected outflow rate, ${\dot{m}}_{\mathrm{depro}}={\dot{m}}_{\mathrm{obs}}\tan \theta$, and also on the inferred time and distance evolution of the outflow rate. We use a Monte Carlo approach to derive the errors introduced by deprojection, assuming that the outflow direction has an equal probability of being in any direction along the surface of the outflow cone.

Figure 8 shows the molecular mass outflow rate as a function of time or distance, corresponding to the two alternative interpretations we discuss above: a flow where the distribution of material is interpreted as resulting from the history of the mass outflow rate (top panel) and one where we show the mass outflow rate as a function of distance, which under the assumption of a constant outflow rate over the last several megayears can be interpreted as the efficiency of ejection to a given distance (bottom panel). Indeed, for an outflow with a distribution of velocities, the slower material will not travel as far in a given time, and neither will it escape the galaxy if it does not have a high-enough velocity. Close to the starburst region (or at small times since ejection), the mass outflow rate drops to zero, because it becomes increasingly difficult to separate the "outflow" component from the "disk" component. At large values of distance or time, it also drops to zero, due to a decreasing amount of outflowing molecular material detected far from the starbursts (and the fact that the observations have a limited field of view). The constant outflow rate out to ∼300 pc is in tension with Kim & Ostriker (2018), who find a steeply dropping "cold" (T < 5050 K) component in their TIGRESS simulation. Within 400 pc, they find the average mass-loading factor to drop by two orders of magnitude. It is unlikely that the SFR in NGC 253 has increased by two orders of magnitude within the past 1–2 Gyr, which would alter the observed constant outflow rate profile to be consistent with the Kim & Ostriker (2018) simulation. Note, however, that their simulation recreates solar neighborhood-like conditions instead of a starburst. A direct comparison may thus not be possible.

Zoom In Zoom Out Reset image size

Our data show that the average outflow rates within 20'' (340 pc) from the major axis are 29 M_⊙ yr⁻¹ (${}_{-0.35}^{+0.48}$ dex), 39 M_⊙ yr⁻¹ (${}_{-0.34}^{+0.49}$ dex), and 4.8 M_⊙ yr⁻¹ (${}_{-0.39}^{+0.50}$ dex) for CO(1–0), (2–1), and (3–2), respectively. Similarly, within the past 1.0 Myr, the average outflow rates are 14 M_⊙ yr⁻¹ (${}_{-0.29}^{+0.25}$ dex), 20 M_⊙ yr⁻¹ (${}_{-0.37}^{+0.27}$ dex), and 2.7 M_⊙ yr⁻¹ (${}_{-0.56}^{+0.22}$ dex) for CO(1–0), (2–1), and (3–2), respectively. The uncertainties, indicated by the 16th–84th percentiles in the Monte Carlo described above, are substantial at a factor of 2–3. Real systematic uncertainties are even larger, because there can be conversion of molecular into atomic material (see Leroy et al. 2015b), or in general variations in the CO-to-H₂ conversion.

Note that the average outflow rates quoted above differ between the two representations, with the median outflow rate as a function of distance being about twice as high as a function of time. The outflowing mass is identical in both cases, and the difference arises solely from binning. Comparing between lines, it is apparent that as measured in CO(3–2), the outflow rate is roughly one order of magnitude lower than for the lower two transitions. This is a direct consequence of the lower mass detected in CO(3–2) and the smaller field of view of those observations. The CO(3–2) observations cover only ∼125 (∼210 pc projected) above/below the disk and thus miss significant amounts of nondisk gas. Their lower surface brightness sensitivity means we also fail to detect a diffuse nondisk component, as we see in the two lower lines. The measurements in CO(3–2) should thus be interpreted as a lower limit, and in that sense they are consistent with those for the lower two transitions.

Overall, the deprojected total mass outflow rate in the starburst of NGC 253 is most likely in the range ∼14–39 M_⊙ yr⁻¹ as derived from CO(1–0) and CO(2–1) with ∼0.4 dex uncertainty. The large spread arises due to different interpretations of the kinematics of the observed gas while the errors are due to unknown geometry. The majority of this outflow rate is contributed by massive outflows alongside the disk like the SW/SE streamers, with a significant contribution by diffuse molecular gas.

The present-day SFR in the central region of NGC 253 is 1.7–2.8 M_⊙ yr⁻¹, derived from radio continuum and far-infrared measurements (Ott et al. 2005; Bendo et al. 2015; Leroy et al. 2015a). This results in a mass-loading factor $\eta ={\dot{M}}_{\mathrm{out}}/{\dot{M}}_{\mathrm{SFR}}$ in the range η ∼ 5.4–23.5. Note that this is for gas ejected as far as 340 pc. We do not currently know what fraction of the gas makes it to the far regions of the halo or reaches escape velocity from the system. Theoretical works suggest that most of the molecular outflow will not escape but rain back down on the galaxy (e.g., Shapiro & Field 1976 up to recent work by Kim & Ostriker 2018 or Tollet et al. 2019). In our data, no gas reaches the escape velocity of v_esc = 500 km s⁻¹ (Walter et al. 2017). The uncertainty on v_esc is substantial, so allowing a factor of 2 is still plausible. At v_esc = 250 km s⁻¹, the fraction of gas above v_esc by mass is 0.5%, 0.5%, and 6.0% for CO(1–0), CO(2–1), and CO(3–2), respectively. The mismatch between the lower transitions and CO(3–2) implies that some high-velocity gas can be found on small scales that is blurred out in the low-resolution observations.

This estimate of the molecular mass outflow rate is higher than the lower limit found by Bolatto et al. (2013) for optically thin emission. Zschaechner et al.'s (2018) analysis of the CO line ratios in the SW streamer shows that the emission there is optically thick, which the authors used to rescale the Bolatto et al. (2013) measurements, finding a NGC 253 galactic outflow rate of 25–50 M_⊙ yr⁻¹. The result presented here, a mass outflow rate of ∼14–39 M_⊙ yr⁻¹, is consistent with this number using an independent and a more complete methodology than the original work.

From Hα observations by Westmoquette et al. (2011), we can estimate the ionized outflow rate to be ∼4 M_⊙ yr⁻¹ using their ionized mass (M = 10⁷ M_⊙) and typical velocity (200 km s⁻¹) at a mean deprojected distance (510 pc). X-ray observations yield comparable values. Strickland et al. (2000) found an upper limit of 2.2 M_⊙ yr⁻¹ assuming a standard outflow velocity of 3000 km s⁻¹. The upper limit reported in Strickland et al. (2002) translates to 2.3 M_⊙ yr⁻¹ when assuming 3000 km s⁻¹ outflow velocity and a reasonable 10% filling factor. These estimates scale linearly with the unknown velocity and also depend on the unknown metallicity and filling factor in the outflow. Estimates of the outflow rate in neutral gas are not known in the literature but are arguably at a similar level. The molecular phase thus clearly dominates the mass budget in the outflow close to the disk as found in other galaxies (e.g., M82, Leroy et al. 2015b; and simulations, e.g., Kim & Ostriker 2018).

Similar to the mass outflow rate, energy and momentum can be calculated as a function of time and distance, which is shown in Figure 9. In Equations (5) and (6), the molecular outflow rate ${\dot{m}}_{i}$ is replaced by kinetic energy ${E}_{\mathrm{kin},i}=\tfrac{1}{2}{m}_{i}\,{v}_{i}^{2}$ or momentum ${P}_{i}={m}_{i}\,{v}_{i}$. As with the outflow rate, the dominant sources of error are the uncertainty in the launching site and the geometry for which we Monte Carlo the errors as described before. Our median estimate with 16th and 84th percentile uncertainties are given below and in Table 2.

Zoom In Zoom Out Reset image size

The kinetic energy in the outflow integrated over the past 1.0 Myr is 3.9 × 10⁵⁴ erg (${}_{-0.75}^{+0.91}$ dex) in CO(1–0), 4.5 × 10⁵⁴ erg (${}_{-0.80}^{+0.94}$ dex) in CO(2–1), and 6.5 × 10⁵³ erg (${}_{-0.83}^{+0.58}$ dex) for CO(3–2). Within 20'' (340 pc), the kinetic energies amount to 2.5 × 10⁵⁴ erg (${}_{-0.65}^{+0.96}$ dex) in CO(1–0), 3.1 × 10⁵⁴ erg (${}_{-0.65}^{+0.96}$ dex) in CO(2–1), and 4.3 × 10⁵³ erg (${}_{-0.64}^{+0.98}$ dex) for CO(3–2). For the reasons described above, the CO(3–2) measurement is a lower limit, thus the results for the lower transitions are consistent.

NGC 253 does not appear to host an energetically important AGN, and the outflow is driven by the starburst. It is interesting then to compare our results for the kinetic energy to the energy released by the starburst. We assume the current SFR of ∼2.8 M_⊙ yr⁻¹ in the central region (Ott et al. 2005; Bendo et al. 2015), which has been approximately constant over the last few megayears.

The total energy E_bol produced by the starburst is simply the time-integrated bolometric luminosity L_bol, which depends¹⁴ on the bolometric magnitude M_bol:

$$\begin{eqnarray}&&{E}_{\mathrm{bol}}={L}_{\mathrm{bol}}\times {\rm{\Delta }}t\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\,3\times {10}^{35-0.4{M}_{\mathrm{bol}}}{\rm{\Delta }}t\left(\displaystyle \frac{\mathrm{SFR}}{1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 8 }$$

According to Starburst99 (Figure 46 in Leitherer et al. 1999), the bolometric magnitude of a starburst at an age of 10⁷–10⁸ yr is M_bol ∼ −20.5 for SFR = 1 M_⊙ yr⁻¹. The total energy output of the starburst over the past 1 Myr is thus 4.2 × 10⁵⁷ erg. The observed kinetic energy ∼3.9–4.5 × 10⁵⁴ erg in the outflow is a factor of ∼10³ lower, which places the coupling efficiency of the outflow kinetic energy to starburst energy at ∼0.1%.

In terms of only kinetic energy, the fraction is higher. It is primarily supernovae and winds that supply kinetic energy to the ISM, which can be estimated from the energy deposition rate according to Leitherer et al. (1999) as given in Chisholm et al. (2017) and Murray et al. (2005):

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{SN}}=3\times {10}^{41}\left(\displaystyle \frac{\mathrm{SFR}}{1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 9 }$$

Each SN releases approximately 10⁵¹ erg in kinetic energy, with the progenitor releasing a similar amount of kinetic energy by winds during its lifetime (e.g., Leitherer et al. 1999). The approximate total kinetic energy released by SNe in the past 1 Myr is then ∼5.3 × 10⁵⁵ erg, compared to the ∼3.9–4.5 ×10⁵⁴ erg we observe in the outflow. Hence, the observed starburst is sufficient to kinetically power the measured molecular outflows with ∼8% efficiency.

The commonly adopted 50% relative contribution of wind feedback is a first-order estimate that is subject to environmental dependence and requires careful modeling to determine precisely (e.g., Leitherer et al. 1999). Furthermore, it should be noted that the observed outflow energy and its error are based on a fixed mass conversion factor that may vary. The uncertainty on the energy coupling efficiency is thus substantial, and it should be understood as an order-of-magnitude comparison.

The above calculation ignores the contribution of other energies, such as the turbulent energy within the molecular outflow and the kinetic energy of the neutral and ionized gases. Matsubayashi et al. (2009) derived a kinetic energy of the ionized wind in NGC 253 of 1.3 × 10⁵³ erg or more than one order of magnitude lower than the molecular outflow kinetic energy. The molecular outflow is slower (50–100 km s⁻¹ on the scales we observed here) than the ionized outflow (up to ∼400 km s⁻¹; Matsubayashi et al. 2009) but also more massive. The ionized outflow thus has only a very small effect on the total kinetic energy and the coupling efficiency.

Deprojected outflow momenta integrated over the past 1.0 Myr are 6.9 × 10⁸ M_⊙ km s⁻¹ (${}_{-0.49}^{+0.50}$ dex), 8.7 ×10⁸ M_⊙ km s⁻¹ (${}_{-0.57}^{+0.57}$ dex), and 1.2 × 10⁸ M_⊙ km s⁻¹ (${}_{-0.59}^{+0.33}$ dex) for CO(1–0), (2–1), and (3–2), respectively. Within 20'' (340 pc) deprojected distance from the launching site, the outflow momenta integrate to 4.8 × 10⁸ M_⊙ km s⁻¹ (${}_{-0.35}^{+0.48}$ dex) in CO(1–0), 6.4 × 10⁸ M_⊙ km s⁻¹ (${}_{-0.34}^{+0.49}$ dex) in CO(2–1), and 8.0 × 10⁷ M_⊙ km s⁻¹ (${}_{-0.39}^{+0.50}$ dex) in CO(3–2).

The momentum released initially by SNe is given in Murray et al. (2005):

$$\begin{eqnarray}&&{\dot{P}}_{\mathrm{SN}}=2\times {10}^{33}\left(\displaystyle \frac{\mathrm{SFR}}{1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\,{\rm{g}}\,\mathrm{cm}\,{{\rm{s}}}^{-2}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&=\,317\,\left(\displaystyle \frac{\mathrm{SFR}}{1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,\mathrm{km}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$$

In 1 Myr, a constant SFR of 2.8 M_⊙ yr⁻¹ yields 8.9 ×10⁸ M_⊙ km s⁻¹. Assuming a contribution by stellar winds of the same order (Leitherer et al. 1999), the total momentum is 1.8 × 10⁹ M_⊙ km s⁻¹, or roughly twice the observed outflow momentum. SNe, however, gain significant amounts¹⁵ of momentum by sweeping up surrounding material. From simulations, the total momentum supplied to the ISM is expected to be 2.8 × 10⁵ M_⊙ km s⁻¹ per SNe (Kim & Ostriker 2015 and references therein). For a constant SFR of 2.8 M_⊙ yr⁻¹ over 1 Myr, this amounts to 1.0 × 10¹⁰ M_⊙ km s⁻¹ or 2.0 × 10¹⁰ M_⊙ km s⁻¹ when adopting 50% contribution by stellar winds, which is about four times the observed momentum. The efficiency of transferring feedback momentum to outflow momentum is thus in the range of 27%–49%, considering the initially available momentum, or a 2.5%–4% efficiency for total to outflow momentum transfer.

These outflow momenta are much higher than the momentum currently produced by young (<10 Myr) super star clusters in the starburst. Leroy et al. (2018) list 14 candidate clusters that together produce 1.5 × 10⁷ M_⊙ km s⁻¹ measured from gas kinematics, a factor 10–100 lower than the observed outflow momentum. The currently forming (super) star cluster thus could not have launched the outflow but the feedback of another population of stars is needed to explain the observed outflows. This is indicative of the time delay of SF feedback.

Energy and momentum curves in Figure 9 differ only by a factor v but follow a similar evolution. This implies that the median velocity at a given distance must be roughly constant along the outflow. As the curves as a function of distance are roughly constant within 50 pc < s < 300 pc, especially for kinetic energy, the outflow mass at a given distance must also be approximately constant along the outflow. The decline in energy and momentum below 50 pc is caused by a decrease in outflow mass, again because both curves follow a similar trend. This is at least partially related to the difficulty of separating the outflow from the disk, where the former emerges from the latter. The decrease could also be interpreted physically as the outflow sweeping up mass while emerging from the disk. An estimation of the relative importance of these effects requires high-resolution modeling of the outflow that are not possible yet because we do not know the detailed outflow geometry. The drop beyond ∼300 pc (∼200 pc in CO(3–2)) is partially related to reaching the edge of the field of view. Discerning this effect from an actual decrease is not possible with our data as we do not know the inclination at every location in the outflow. The edge of the field of view thus corresponds to a range of deprojected distances from the disk, which gradually depresses the curve rather than showing a sudden drop. A physical reason for the decrease could be the destruction of the molecular gas, e.g., photodissociation by the intense starburst radiation or ionization.

The kinetic energy and momentum evolution in Figure 9 thus suggest both energy and momentum conservation along the outflow from ∼50 pc to ∼300 pc, as well as approximately constant molecular gas mass.

When additionally assuming no acceleration of the outflow after launch, it becomes possible to study the time evolution. The corresponding plots (Figure 8 top and 9 top) all show a peak within the past 0.5 Myr and steady decrease toward earlier gas ejection times. Corresponding to the decline toward zero distance, the decrease toward zero ejection time is most likely a methodological complication. From the peak at t_eject = 0. 2–0.3 Myr, kinetic energy and momentum in the outflow drop by a factor of 10 within ∼2 Myr. This decline would be physically plausible if the starburst in NGC 253 is very young and takes into account a time delay between the start of star formation, feedback, and efficient outflow driving (superbubble breakout). For the observed age of the starburst of 20–30 Myr (Rieke et al. 1980; Engelbracht et al. 1998) this scenario is implausible. Time delays of >20 Myr are longer than the lifetime of the most massive stars. A younger generation of massive stars at an age of ∼6 Myr (Kornei & McCrady 2009) may, however, drive the currently visible molecular outflows. If this were to be true, a time delay between star formation and outflow launching of ∼4 Myr is implied. Outflow launching in this context means the time after which the outflow reaches a mass loading η > 1. The time delay is 2 Myr until the outflow carries more energy (momentum) than the feedback kinetic energy (momentum) of a single high-mass star. Note that these rough estimates depend on the assumption of no acceleration (neither positive nor negative) of the outflow after being launched from the disk, which might be a close enough approximation on these scales of a few hundred parsecs. The very young (≲1 Myr) and still deeply embedded super star cluster discussed recently by Ando et al. (2017) and Leroy et al. (2018) are most likely too young to have affected the observed molecular outflow.