Outflows from Super Star Clusters in the Central Starburst of NGC 253

In this paper, we focus on the line profiles, specifically the CS 7−6 and H¹³CN 4−3 lines toward three of the SSCs. We will present the continuum data more fully in a forthcoming paper (R. C. Levy et al. 2021, in preparation). Some of the data are presented and used here, and we provide a brief summary.

To extract the 350 GHz continuum data, we flagged channels that may contain strong lines in the band, assuming a systemic velocity of 243 km s⁻¹ (Koribalski et al. 2004). Lines included in the flagging are ¹²CO 3−2, HCN 4−3, H¹³CN 4−3, CS 7−6, HCO⁺ 4−3, ²⁹SiO 8−7, and ¹²CO 3−2 v = 1 (though this line is not detected), and channels within ±200 km s⁻¹ of the rest frequencies of these lines were flagged. We imaged the line-flagged visibilities using the CASA version 5.4.1 tclean task with specmode = 'mfs', deconvolver='multiscale', scales=[0,4,16,64], Briggs weighting with robust = 0.5, and no uv-taper. The baseline was fit with a linear function (nterms = 2) to account for any change in slope over the wide band. The continuum was cleaned down to a threshold of 27 μJy beam⁻¹, after which the residuals resembled noise. This map has a beam size of 0.024'' × 0.016'' and a cell (pixel) size of 0.0046'' (0.078 pc), so that we place around 4 pixels across the minor axis of the beam. Finally, the map was convolved to a circular beam size of 0.028'' (0.48 pc). The rms noise of the continuum image (away from emission) is 26 μJy beam⁻¹ (0.3 K).

From this high-resolution 350 GHz continuum image (Figure 1), we identify approximately three dozen clumps of dust emission. These are co-spatial with the 14 proto-SSCs identified by Leroy et al. (2018), with many sources breaking apart into smaller clumps that likely represent individual clusters. We note that there is appreciable overlap in the spatial scales probed by these observations and those presented by Leroy et al. (2018). This means that the larger cluster sizes measured by Leroy et al. (2018) are due to the lower resolution, with the true cluster sizes better traced by these observations. We further verify this by computing the fluxes in the high-resolution continuum image in the same apertures used by Leroy et al. (2018, see their Table 1). For SSCs 4, 5, and 14 (the focus of this analysis), we recover >80% of the reported flux.

Zoom In Zoom Out Reset image size

We follow the SSC naming convention of Leroy et al. (2018), adding letters to denote subclusters in decreasing order of brightness. For clusters that break apart, there is always a main, brightest cluster. For example, SSC 4 from Leroy et al. (2018) breaks apart into six dust clumps, where the peak intensity of SSC 4a is >4× that of any of the smaller clumps. For this paper, we will focus only on the primary clusters since those are the most massive, and smaller dust clumps may not be true clusters. We also remove SSC 6 from this analysis, as it is extended and more tenuous in the high-resolution continuum map and so may not be an SSC. Leroy et al. (2018) also noted that this clump is weak and may instead be a supernova remnant, which is further supported by a non-negligible synchrotron component to this cluster's spectral energy distribution (SED; E. A. C. Mills et al. 2021, in preparation).

The position and orientation of each primary cluster is measured by fitting a 2D Gaussian to each cluster in the continuum image using emcee (Foreman–Mackey et al. 2013). The reported values are the medians of the marginalized posterior likelihood distributions (Table 1). As a nonparametric measurement of the cluster sizes, we compute the half-flux radius (r_half−flux) for each cluster. Given the cluster positions and orientations from the 2D Gaussian fit, we construct elliptical annuli in steps of half the beam FWHM. The half-flux radius is the median of the cumulative flux distribution. We use this radius measurement as opposed to those from the 2D Gaussian fit because the cluster light profiles are not necessarily well described by a Gaussian, which will be discussed further in a forthcoming paper (R. C. Levy et al. 2021, in preparation). We deconvolve the beam from the radius measurements by removing half the beam FWHM from the half-flux radius in quadrature. The beam deconvolution will be treated more fully in a forthcoming paper focusing on the continuum properties of the SSCs, which will also include the lower resolution data (R. C. Levy et al. 2021, in preparation).

We image the CS 7−6 and H¹³CN 4−3 lines by selecting a 800 km s⁻¹ window around the line center (rest frequencies of 342.883 and 349.339 GHz, respectively), assuming a systemic velocity of 243 km s⁻¹ (Koribalski et al. 2004). We clean the lines of interest using tclean in CASA version 5.4.1 with specmode='cube', deconvolver='multiscale', scales=[0], Briggs weighting with robust=0.5, and no uv-taper. No clean mask was used. The baseline was fit with a linear function (nterms=2) to account for any change in slope over the wide band. The continuum is not removed from the visibilities or the cleaned cubes. The cell (pixel) size is 0.0046'' (0.078 pc), so that we place around 4 pixels across the minor axis of the beam. The final elliptical beam is 0.027'' × 0.019''. Finally, we convolve to a circular 0.028'' (0.48 pc) beam. The rms noise is 0.5 mJy beam⁻¹ (6.5 K) per 5 km s⁻¹ channel. Figure 2 shows the peak intensity maps for CS 7−6 and H¹³CN 4−3 within ±200 km s⁻¹ about the galaxy's systemic velocity to avoid possible contamination by other strong lines (especially by HC₃N in the case of H¹³CN 4−3).

Zoom In Zoom Out Reset image size

Though the analysis is focused on the CS 7−6 and H¹³CN 4−3 lines, we also make a "dirty" cube covering the entire band and imaged area by running tclean with niter=0. This dirty cube has a cell (pixel) size of 0.02'' and an rms noise away from the emitting regions of 0.57 mJy beam⁻¹ per 5 km s⁻¹ channel. The synthesized beam is 0.05'' × 0.025''. The continuum is not removed from the visibilities or the dirty cube.

In Figure 3, we show the full band spectrum of SSCs 4a, 5a, and 14, extracted from the dirty cube and averaged over the FWHM continuum source size. There are several immediately striking features in these spectra. The brightest lines in SSCs 4a, 5a, and 14 show deep, blueshifted absorption features, extending down to ∼0 mJy beam⁻¹, and redshifted emission. This line shape is commonly referred to as a P-Cygni profile—for the star in which it was first detected—and is indicative of outflows. As a comparison, we also show the spectrum of SSC 8a in Figure 3, which does not have these P-Cygni line shapes. We note that the CO 3−2 line has a complicated structure, likely because it traces lower density gas than other detected lines (n_crit ∼ 2 × 10⁴ cm⁻³; e.g., Turner et al. 2017) and so it arises in many places along the line of sight to the SSCs.

Zoom In Zoom Out Reset image size

There are many detected lines in addition to the usual dense gas tracers. These include shock tracers such as SO 8−7 and SO₂ lines and vibrationally excited HC₃N and HCN 4−3, which are primarily excited through IR pumping (e.g., Ziurys & Turner 1986; Krieger et al. 2020). Many of these have been previously identified by Krieger et al. (2020), some of which are marked in Figure 3.

Because the SSCs are embedded in the nucleus of NGC 253, there is dense gas along the line of sight that may contribute to the observed spectra but is not associated with the SSCs themselves. This can be seen in the line channel maps around each SSC, shown in Figure 4 for CS 7−6. This gas is most evident in SSCs 5a and 14, which both show changes in the gas morphology with velocity (relative to the cluster systemic velocity; see Section 2.5). The first panel of each set of channel maps shows the CO 3−2 velocity field in the same region as a tracer of the bulk gas motions not associated with the clusters themselves. Because we are interested in the large-scale gas motions, we smooth the CO 3−2 velocity maps presented by Krieger et al. (2019) to 5× the beam major axis (5 × 0.17'' = 0.85'' = 14.4 pc). The velocities shown in Figure 4 are relative to the cluster systemic velocities. All three clusters are located at different velocities compared to the bulk gas motions (i.e., CO 3−2 velocity maps are not centered on zero). In the case of SSC 14, the distinct morphological evolution with velocity matches qualitatively with the CO 3−2, though there is a velocity offset. It is difficult to know what fraction of this is associated with the clusters themselves or with dense gas in the environment of the clusters but not necessarily bound to them. Here, we will assume that all of the emission in the green annuli in Figure 4 is not bound to the clusters and should be corrected to obtain the intrinsic spectrum of the clusters. The spectrum of the environmental gas (T_env) is averaged in an annulus with an inner radius of 3× the half-flux radius and an outer radius of 5× the half-flux radius, as shown in Figure 4. This effectively assumes the environmental dense gas forms a uniform screen in between the cluster and our line of sight. This is unlikely to be the case and is one of the uncertainties of this correction.

Zoom In Zoom Out Reset image size

We calculate the optical depth (τ_ν ) of the environmental gas following Mangum & Shirley (2015). Combining their Equations (24) and (27), we have

$$\begin{eqnarray}&&{\tau }_{\nu }=-\mathrm{ln}\left[1-\displaystyle \frac{{{kT}}_{\mathrm{env}}}{h\nu }{\left(\displaystyle \frac{1}{{e}^{h\nu /{{kT}}_{\mathrm{ex}}}-1}-\displaystyle \frac{1}{{e}^{h\nu /{{kT}}_{\mathrm{bg}}}-1}\right)}^{-1}\right],\end{eqnarray} \tag{ 1 }$$

where the filling fraction is assumed to be unity, the excitation temperature (T_ex) is 102 ± 28 K, and the background temperature (T_bg) is measured from the continuum level of the environmental spectrum (but is ≈0 K; Figure 5). The assumed value of T_ex is the average of the excitation temperatures found by Meier et al. (2015) over larger scales in the nucleus (74 K) and by Krieger et al. (2020) for regions more localized to the SSCs (130 K), since the true value of the excitation temperature of this dense gas near the clusters likely falls somewhere between these two measurements. The uncertainty reflects half of the difference of these two values.

Zoom In Zoom Out Reset image size

We assume that half of the environmental emission comes from the foreground, so that the optical depth of the material along the line of sight in front of the cluster is τ_ν,fg = τ_ν /2. We also assume that the environmental gas is colder than the gas in the cluster, so it will absorb emission from the cluster. We then derive the corrected (intrinsic) SSC spectra

$$\begin{eqnarray}\,\begin{array}{rcl}{T}_{\mathrm{obs}} & = & {T}_{\mathrm{intrinsic}}{e}^{-{\tau }_{\nu ,{fg}}}+{T}_{\mathrm{env}}\left(1-{e}^{-{\tau }_{\nu ,{fg}}}\right)\\ \Rightarrow {T}_{\mathrm{intrinsic}} & = & {e}^{{\tau }_{\nu ,{fg}}}\left[{T}_{\mathrm{obs}}-{T}_{\mathrm{env}}\left(1-{e}^{-{\tau }_{\nu ,{fg}}}\right)\right],\end{array}\,\end{eqnarray} \tag{ 2 }$$

where the second equation is a rearrangement of the first. The observed, environmental, and intrinsic spectra of CS 7−6 for SSCs 4a, 5a, and 14 are shown in Figure 5. This process is repeated for H¹³CN 4−3 and the full band spectra in Figure 3.

As discussed earlier, it is difficult to determine what fraction of the extended emission seen in the channel maps (Figure 4) is environmental or associated with the cluster. Moreover, our method for the foreground correction is quite simplistic, and there are uncertainties related to the fraction of gas that may be associated with the clusters, the geometry of the environmental material, and the assumed excitation temperature. As can be seen in Figure 5, however, the correction mostly affects the peak of the emission components. After the foreground correction, the emission peak for SSC 14 increases by 7% (17 K), whereas the absorption trough only increases by 4.5% (1.5 K). There is no appreciable change in the spectra for SSCs 4a and 5a. The absorption features are, therefore, largely unaffected by this correction, and the derived properties of the outflow (presented in Sections 3.1 and 3.2) are robust. The intrinsic (corrected) spectra are used throughout the analysis of this paper.

We constrain the systemic velocity of all of the detected clusters using the full band spectra and the detected line list presented by Krieger et al. (2020). Using the full band spectrum for each cluster, we simultaneously fit each line in the line list using a series of Gaussians of the form

$$\begin{eqnarray}&&I(\nu )={I}_{\mathrm{cont}}+\displaystyle \sum _{\mathrm{line}}{I}_{\mathrm{line}}\exp \left[\displaystyle \frac{-{\left(\nu -{\nu }_{\mathrm{line} \mbox{-} \mathrm{rest}}\right)}^{2}}{2{\sigma }_{\mathrm{line}}^{2}}\right],\end{eqnarray} \tag{ 3 }$$

where I_cont is the fitted continuum level over the whole spectrum, I_line is the fitted intensity of each line, ν_line‐rest is the fixed rest frequency of each line, σ_line is the fitted width of each line, and ν is the rest frequency of the observed spectrum, which depends on the assumed systemic velocity. The primary lines used to determine the systemic velocity are marked in orange in Figure 3; these lines tend to have strong emission and simple line shapes. Other strong lines in the band (e.g., CO 3−2, HCN 4−3, HCO⁺ 4−3) tend to have complicated line shapes that make accurately determining the systemic velocity difficult. Due to the presence of lines with complicated shapes and line blending, the cross correlation of each full band spectrum and the multi-Gaussian fit is done by eye. The systemic velocities of for all the clusters are listed in Table 1. The uncertainty is estimated to be ±5 km s⁻¹ . For SSCs 4a, 5a, and 14, the systemic velocities are further refined and confirmed using other lines in the cleaned CS 7−6 and H¹³CN 4−3 cubes (Figure 6). Several of these lines—the vibrational transitions, in particular—have sharp peaks and are spatially localized to the clusters themselves. They, therefore, provide tighter constraints on the cluster systemic velocities. For SSCs 4a, 5a, and 14, we estimate that our cluster systemic velocities are accurate to ±1 km s⁻¹ .

Zoom In Zoom Out Reset image size

In addition to cluster properties derived from these new 0.5 pc resolution data described in the previous subsections, we use stellar masses, total H₂ gas masses, and escape velocities of the clusters calculated by Leroy et al. (2018). Here, we briefly describe how these quantities were derived. These quantities are reproduced in Table 2 for the three clusters that are the focus of this work.

Table 2. Cluster and Outflow Properties from Absorption Line Fits

		SSC 4a		SSC 5a		SSC 14
Quantity	Unit	CS 7−6	H13CN 4−3	CS 7−6	H13CN 4−3	CS 7−6	H13CN 4−3
${\mathrm{log}}_{10}({M}_{{{\rm{H}}}_{2},\mathrm{tot}})$ a	log10(M⊙)	5.1		5.3		5.7
${\mathrm{log}}_{10}({M}_{* })$ b	log10(M⊙)	5.0		5.4		5.5
Vescape c	km s−1	22.0		33.2		45.5
τ d		2.3	2.3	0.2	0.0	1.8	0.6
${V}_{\max \mbox{-} \mathrm{abs}}$ e	km s−1	6	6	22	20	24	24
${V}_{\mathrm{out},\max }$ f	km s−1	40	46	42	34	51	56
ΔVout,FWHM g	km s−1	40	46	24	17	32	38
log10(tcross) h	log10(yr)	5.0	5.0	4.5	4.6	4.4	4.4
${\mathrm{log}}_{10}($ ${N}_{{{\rm{H}}}_{2},\mathrm{out}})$ i	log10(cm−2)	23.9	25.3	22.7	23.2	23.7	24.6
${\mathrm{log}}_{10}({M}_{{{\rm{H}}}_{2},\mathrm{out}})$ j	log10(M⊙)	4.7	6.1	3.8	4.3	4.6	5.5
${\mathrm{log}}_{10}($ ${M}_{{{\rm{H}}}_{2},\mathrm{out}}$/${M}_{{{\rm{H}}}_{2},\mathrm{tot}}$) k		−0.4	1.0	−1.5	−1.0	−1.1	−0.2
${\mathrm{log}}_{10}($ ${M}_{{{\rm{H}}}_{2},\mathrm{out}}$/M*) k		−0.3	1.1	−1.6	−1.1	−0.9	−0.0
${\mathrm{log}}_{10}({\dot{M}}_{{{\rm{H}}}_{2},\mathrm{out}})$ l	log10(M⊙ yr−1)	−0.3	1.1	−0.7	−0.3	0.2	1.1
log10(tremove‐gas) m	log10(yr)	5.4	4.0	6.0	5.6	5.5	4.6
${\mathrm{log}}_{10}($ pr/M*) n	log10(km s−1)	0.7	2.1	0.0	0.5	0.7	1.6
log10(Ekin) o	log10(erg)	49.2	50.6	49.6	49.9	50.4	51.3

Note. Values (aside from the top three) are derived from fits to the absorption component of the line profiles. See Section 3.1 and Appendix A for details on how these values were calculated.

^a The cluster gas mass measured from the dust continuum and assuming a gas-to-dust ratio of 100 (Leroy et al. 2018). Uncertainties are ∼0.4−0.5 dex, with these measurements likely biased low due to assumptions about the dust temperature, dust-to-gas ratio, and the dust opacity. See Section 2.6 for a detailed discussion of how these values were calculated. ^b The cluster zero age main-sequence (ZAMS) stellar mass measured from the 36 GHz free–free emission (Leroy et al. 2018). See Section 2.6 for a detailed discussion of the uncertainties on these stellar masses. ^c The escape velocity from the cluster (Leroy et al. 2018). Uncertainties are dominated by those from the gas and stellar masses. See Section 2.6 for a detailed discussion of how these values were calculated. ^d The optical depth in the outflow calculated from the absorption to continuum ratio (Equation (A4)). ^e The velocity where the absorption is maximized, corresponding to the velocity at which the bulk of the material traced by CS 7−6 and H¹³CN 4−3 is outflowing. ^f The maximum outflow velocity, defined as 2σ from the mean outflow velocity (Equation (A3)). For a Gaussian line, this means that 95% of the dense material has an outflow velocity slower than ${V}_{\max ,\mathrm{out}}.$ ^g The FWHM line width from the Gaussian fit to the absorption feature. ^h The gas crossing time, or the time it would take a parcel of gas to travel from the center of the cluster to r_SSC at ${V}_{\max \mbox{-} \mathrm{abs}}$ (Equation (A2)). ⁱ The H₂ column density in the outflow derived from τ (Equations (A5)–(A9)). This calculation assumes local thermal equilibrium (LTE) with an excitation temperature of 130 ± 56 K and abundances of CS 7−6 and H¹³CN 4−3 relative to H₂ from Martín et al. (2006). The abundances may vary by an order of magnitude on these small scales, so all quantities that depend on the column density are limited to order-of-magnitude precision. The blueshifted outflow component only probes the portion of the outflow on the approaching side, so the values are multiplied by two to account for the receding side of the outflow, assuming it is identical to the approaching side. Uncertainties are ±1 dex. ^j The H₂ mass in the outflow derived from ${N}_{{{\rm{H}}}_{2}},\mathrm{out}$ and the projected cluster size (Equation (A10)). This calculation and others that depend on it assume the outflow is spherical, which is supported by the results of our modeling (Section 3.2). Uncertainties are ±1 dex. ^k The H₂ mass in the outflow compared to the total gas or stellar mass in the cluster. Uncertainties are ±1 dex. ^l The mass outflow rate derived over one crossing time (Equation (A11)). Uncertainties are ±1 dex. ^m The gas-removal time, or the time it would take for the current mass outflow rate to expel all of the cluster gas mass (${M}_{{{\rm{H}}}_{2}},\mathrm{tot}$) from the clusters (Equation (A12)). Uncertainties are ±1dex. ⁿ The radial momentum carried by the outflow per unit stellar mass, which also assumes the outflow is spherical (Equation (A13)). Uncertainties are ±1 dex. ^o The kinetic energy of the outflow calculated from ${M}_{{{\rm{H}}}_{2}},\mathrm{out}$ and ${V}_{\max \mbox{-} \mathrm{abs}}$ (Equation (A14)). Uncertainties are ±1 dex.

Download table as:  ASCIITypeset image

Stellar Masses (M_* )—The stellar masses are calculated based on the 36 GHz image from the Very Large Array (VLA) convolved to 0.11'' (1.9 pc) resolution (Gorski et al. 2017, 2019; Leroy et al. 2018). At 36 GHz, the emission is assumed to be entirely due to free–free (Bremsstrahlung) emission. From the 36 GHz luminosity, Leroy et al. (2018) derived M_* for a zero-age main sequence (ZAMS) population with a standard initial mass function (IMF; see their Section 4.3.1 and Table 2). As described below, we consider additional sources of uncertainty beyond those described by Leroy et al. (2018). First, we assume that the primary clusters retain all of the stellar mass measured by Leroy et al. (2018), which may lead to an overestimate of M_* for those clusters that break apart. We estimate that this is at worst 50% in those cases (although the higher resolution data presented here shows a number of satellite structures, the original one remains dominant in most). Second, Leroy et al. (2018) estimated that there is a ±20% systematic uncertainty due to assumptions about the Gaunt factor. Finally, E. A. C. Mills et al. (2021, in preparation) measured the M_* of these clusters at 5 pc resolution using hydrogen RRLs. These agree with the M_* estimated by Leroy et al. (2018) to within ±0.3 dex on average. For SSCs 4a, 5a, and 14, the agreement is even better, and the RRL analysis finds M_*, which differ by +0.1, −0.1, and −0.2 dex from those derived from the 36 GHz emission, where + (−) means the RRL measurements produce a larger (smaller) M_*. Included in the calculations by E. A. C. Mills et al. (2021, in preparation) is the effect of a synchrotron component. If synchrotron contaminates the 36 GHz flux (which is assumed to be entirely free–free in the M_* calculation of Leroy et al. 2018), this would lead to an overestimate of M_*. There is negligible synchrotron emission in SSCs 4a, 5a, and 14, but this could affect the M_* of other clusters (especially SSCs 1a, 10a, 11a, and 12a; E. A. C. Mills et al. 2021, in preparation). To determine the lower error bar on M_*, we take the largest of the above three sources of uncertainty for each cluster, where the uncertainty from the RRL measurements are included only if they yield a smaller M_*. To determine the upper error bars, we take the larger of the 20% systematic uncertainty and RRL-measured M_* where they yield a larger M_*. In addition to these quantifiable uncertainties on M_*, there are unquantified uncertainties relating to ionizing photons absorbed by dust and evolution beyond the ZAMS, both of which would result in the reported stellar masses being underestimates (Leroy et al. 2018). While we do not attempt to calculate these unquantifiable uncertainties, we caution that they could be important. In this work, we use these stellar masses primarily to place the clusters in a mass–radius diagram and to evaluate the possible mechanisms powering the observed outflows.

Total H₂ Gas Mass ( ${M}_{{{\rm{H}}}_{2}},\mathrm{tot}$ )—The cluster gas masses are based on the 350 GHz image from ALMA at 0.11'' (1.9 pc) resolution (Leroy et al. 2018). At 350 GHz, the majority of the emission is due to thermal dust emission, though there may be a small level of free–free contamination. Leroy et al. (2018) quantified this free–free contamination and corrected the measured 350 GHz fluxes to determine a more accurate flux due to thermal dust emission (see their Section 4.1 and Table 1). Assuming a dust temperature of 130 K, dust-to-gas ratio of 1:100, and a dust mass absorption coefficient of 1.9 cm² g⁻¹, Leroy et al. (2018) calculated the total H₂ mass of the clusters. Uncertainties are ∼0.4−0.5 dex, with these measurements likely biased low due to assumptions about the dust temperature, dust-to-gas ratio, and the dust opacity. See their Section 4.3.3 and Table 2 for more details. In this work, we consider these to be the total gas masses of the clusters, including gas still bound within the cluster and outflowing from it.

Escape Velocity ( V_escape )—Leroy et al. (2018) calculated the escape velocity from the clusters using their measured stellar and gas masses and their measured sizes from Gaussian fits to the 0.11'' (1.9 pc) resolution 350 GHz image. Uncertainties are dominated by those from M_* and ${M}_{{{\rm{H}}}_{2}},\mathrm{tot}$. See Section 4.3.6 and Table 2 of Leroy et al. (2018) for more details. In this work, we use these escape velocities to compare against the outflow velocities we measure to place constraints on whether the outflowing material will escape the cluster.

We analyze the spectra toward all of the clusters (primary or otherwise) identified in the high-resolution continuum data (Figure 1) in search of spectral outflow signatures. We determine our final list of outflow candidates based on robust detections of blueshifted absorption and redshifted emission in the CS 7−6 and H¹³CN 4−3 lines. We present SSCs 4a, 5a, and 14 as our conservative list of detected P-Cygni profiles indicative of outflows. The CS 7−6 and H¹³CN 4−3 spectra toward SSCs 4a, 5a, and 14 are shown in Figure 6, which are from the cleaned line cubes (Section 2.2) and have been corrected for foreground gas (Section 2.4).

While other SSCs show absorption features, they are not due to outflowing material as the absorption occurs at the line center. Absorption at the line center is likely due to self-absorption of the cold material surrounding the cluster against the embedded continuum source(s). Redshifted absorption would be indicative of material inflowing toward the cluster; a hint of inflow is perhaps seen toward SSC 8a (Figure 3, especially the CS 7–6, HCN 4−3, and HCO⁺ 4−3 lines). Blueshifted absorption, on the other hand, means that the cold foreground material is outflowing from the cluster. For this analysis, we are focused on the blueshifted absorption indicative of outflows.

We also verify that these absorption features are not induced by the spatial filtering by the interferometer. Due to the incomplete sampling of the Fourier plane and the lack of short spacings, emission on larger scales is filtered out and cannot be properly deconvolved, causing an increase in the rms of the data at the velocities where bright extended emission would be present. This problem manifests itself as negative "bowls" around bright emission, but also as negative features that have velocity structure. While this mainly affects the continuum (which is not removed from these data), it can also induce false absorption features in the affected velocity range. This can be seen, for example, in the CO 3−2 lines in Figure 3, where the absorption profiles are very complex. Our choice to focus on CS 7−6 and H¹³CN 4−3 mitigates this problem, as these transitions require high densities to be excited and are hence fairly localized to the clusters themselves, with not much of an extended component. As shown in Figure 5, the environmental CS 7−6 spectra around the SSCs are either flat or show the line in emission; the H¹³CN 4−3 spectra show similar profiles. Therefore, we conclude that the CS 7−6 and H¹³CN 4−3 absorption features are not a result of the spatial filtering by the interferometer, although it may contribute to uncertainties at the 10% level.

A careful measurement of the cluster systemic velocities (Section 2.5) is critical to determine whether the absorption is blueshifted with respect to the cluster or not. Critically, for a few clusters (e.g., SSCs 3a and 13a) there are pathological combinations of absorption at the line center coupled with strong vibrationally excited lines in emission at the redshifted edge of the absorption feature, which together could masquerade as P-Cygni profiles. We identify other line candidates present in the spectra around CS 7−6 and H¹³CN 4−3 using the line list presented by Krieger et al. (2020) as well as Splatalogue, ¹⁴ which are listed in the yellow shaded regions in Figure 6. It is outside the scope of this paper to verify precisely which other lines are present in the spectra, so we present them simply as possible candidates (or combinations of candidates), which confuse the CS 7−6 and H¹³CN 4−3 spectra. This line blending can lead to false detections of P-Cygni profiles. For example, H¹³CN 4−3 v₂ = 1 is located 100 MHz (≈88 km s⁻¹) redward of H¹³CN 4−3, as can be seen in the right column of Figure 6. Similarly, HCN 4−3 v₂ = 1 is located only 45 MHz (≈38 km s⁻¹) redward of HCN 4−3. Given the broad line profiles (in emission and absorption), these lines can easily blend together to produce a facsimile of a P-Cygni profile. A careful measurement of the cluster systemic velocities (using other lines with simpler profiles), cross referencing with the line list presented by Krieger et al. (2020), and further verification using Splatalogue (see footnote 14) reveals that these are not true P-Cygni profiles. It is also worth noting that there is an emission feature (likely OS¹⁸O, an SO₂ isotopologue) seen in the CS 7−6 absorption trough of SSC 4a (upper left panel of Figure 6).

That we detect outflows from SSCs 4a, 5a, and 14 is in large part due to our ability to spatially resolve the individual clusters, since the spectral signatures of the outflows are localized roughly to the central half-flux radius of the clusters. In retrospect, there are indications of these outflows in the CS 7−6 and H¹³CN 4−3 line profiles presented by Leroy et al. (2018, see their Figure 3) in SSCs 4 and 14. The much weaker outflow in SSC 5a is not apparent in these lower resolution data where the clusters are only marginally resolved. There are also kinematic offsets between HCN 4−3 and H40α measurements toward these clusters which are consistent with the effects of these outflows at lower resolution (E. A. C. Mills et al. 2021, in preparation).

Here, we briefly investigate whether SSCs 4a, 5a, and 14 exhibit signs of internal rotation. Observational detections of rotation in star and globular clusters in the Milky Way are mixed (e.g., Kamann et al. 2019; Cordoni et al. 2020), though simulations predict that massive star clusters (>1000 M_⊙) should have appreciable rotation (Mapelli 2017).

We investigate the internal kinematics by constructing position–velocity (PV) diagrams over a range of angles as well as peak-velocity and intensity-weighted velocity (moment 1) maps around each of the clusters. While we see no clear signs of internal rotation, the absorption hinders this analysis significantly as it affects the (effective) intensity weighting of the peak- and intensity-weighted velocity (moment 1) maps. It similarly confuses the interpretation of the PV diagrams. Therefore, we cannot claim that the SSCs are not rotating, only that, given the presence of absorption, we see no clear evidence for rotation. We also check clusters without absorption or outflow signatures, and also find no evidence of rotation.

For a first estimate, we measure the outflow velocities, column densities, gas masses, and momentum of each outflow by fitting the profiles of H¹³CN 4−3 and CS 7−6 with a two-component Gaussian (blue dashed curves in Figure 6). We exclude the portions of the spectra that may be contaminated by other lines, as marked by the yellow shaded regions in Figure 6. The fit to the absorption feature yields the velocity at the maximum absorption (${V}_{\max \mbox{-} \mathrm{abs}}$), or the velocity at which the bulk of the material is outflowing. We define the maximum outflow velocity as 2σ from the average velocity (${V}_{\max ,\mathrm{out}}\equiv {V}_{\max \mbox{-} \mathrm{abs}}+2\tfrac{{\rm{\Delta }}{V}_{\mathrm{out},\mathrm{FWHM}}}{2.355}$), which means that 95% of the material traced by CS 7−6 and H¹³CN 4−3 has an outflow velocity slower than ${V}_{\max ,\mathrm{out}}$ (for a truly Gaussian line). Given the mean outflow velocity and the cluster size, we calculate the crossing time, or the time for a parcel of gas to travel from the center of the cluster to the edge (${t}_{\mathrm{cross}}\equiv {r}_{\mathrm{SSC}}/{V}_{\max -\mathrm{abs}}$), where r_SSC is equivalent to the half-flux radius (Table 1). From the fitted absorption to continuum intensity ratio, we derive the optical depth and hence the column density in absorption. We convert from the column density of the molecule to H₂ (${N}_{{{\rm{H}}}_{2}},\mathrm{out}$) using abundances from Martín et al. (2006). From the column density and projected size (Table 1), we estimate the H₂ mass in the outflow (${M}_{{{\rm{H}}}_{2}},\mathrm{out}$) assuming the outflow is spherical. Constraints on the opening angle of the outflows derived from modeling of the line profiles are discussed in Section 3.2 and show that the outflows must be nearly spherical. Given the mass, crossing time, and mean outflow velocity, we calculate the mass outflow rate (${\dot{M}}_{{{\rm{H}}}_{2},\mathrm{out}}$), the gas-removal timescale (t_remove‐gas), the radial momentum injected per unit stellar mass in the cluster (p_r /M_*), and the kinetic energy of the outflow (E_kin). Further details and equations used to derive these quantities are presented in Appendix A. The outflow parameters derived from the absorption profile fits are reported in Table 2.

The clusters have outflows with maximum velocities (${V}_{\max ,\mathrm{out}}$) of ∼40–50 km s⁻¹. For all three SSCs with outflows, ${V}_{\max ,\mathrm{out}}$ is larger than the escape velocity (V_escape; reproduced in Table 2 from Leroy et al. 2018 and see also Section 2.6). The bulk of the material traced by CS 7−6 and H¹³CN 4−3, however, has velocities less than V_escape (as traced by ${V}_{\max \mbox{-} \mathrm{abs}}$). For SSCs 4a, 5a, and 14, 20%, 7%, and 7% of the outflowing material has velocities larger than the escape velocity and will be able to escape from the cluster. Gas that is outflowing with velocities below V_escape may be reaccreted by the cluster. The gas crossing time (t_cross) is ∼few × 10⁴ yr, which is a lower limit on the age of the outflow. These timescales and the possibility of reaccretion will be discussed further in Section 4.1.

The masses in the outflows are large. The molecular hydrogen column densities and masses derived from H¹³CN 4−3 are larger than those derived from CS 7−6 in all cases. One possibility is that CS 7−6 is more optically thick than H¹³CN 4−3, which may be supported by the depth of the absorption in SSC 4a, but this cannot explain the discrepancy in SSC 5a. A more likely possibility is that the relative molecular abundances we assume are not correct for these small, extreme regions. We adopt molecular abundances for H¹³CN and CS with respect to H₂ of [H¹³CN]/[H₂] = (1.2 ± 0.2) × 10⁻¹⁰ and [CS]/[H₂] = (5.0 ± 0.6) × 10⁻⁹ from a study of the the center of NGC 253 at 14''–19'' (∼240–320 pc) resolution by Martín et al. (2006), where the brackets refer to the abundance of that species. Adopting these abundances on the parsec scales of our clusters, however, is highly uncertain. Observations of envelopes around high-mass young stellar objects in the Milky Way, for example, reveal order-of-magnitude variations in the abundances of molecules, especially nitrogen and sulfur bearing species (van Dishoeck & Blake 1998, and references therein). The impact of the uncertainty on the fractional molecular abundances is that they translate into order-of-magnitude accuracy for the H₂ column density measurements, as well as the subsequent calculations, which depend on the column density, as reported in Table 2. In the following subsection we discuss in more detail what we know about chemical abundances and their variation: our conclusion is that the masses derived from H¹³CN 4−3 are likely more accurate.

3.1.1. Fractional Abundance Variations

The short gas crossing times, large outflow masses, and outflow mass rates (Table 2) suggest the outflow activity in these objects cannot be sustained for a long time. From this standpoint, it is reasonable that we detect outflows in ≲10% of the SSCs, so it is possible that we are indeed catching a small fraction of the SSCs in this short-lived phase. The analysis in Section 3.1, however, implicitly assumes that the outflows are spherical. Another possibility, however, is that the outflows are biconical with a more-or-less narrow opening angle. In that case the outflows from the three clusters we identify are serendipitously pointed close enough to the line of sight to make them detectable. If the observed outflows are not spherical, geometric correction factors will need to be applied to the measured quantities in Table 2, and more outflows could exist that we do not detect because their geometry is unfavorable.

In order to determine whether the observed outflows are spherical or biconical and how they are oriented with respect to the line of sight, we build a simple radiative transfer model to model the spectrum through the outflow. ¹⁵ We consider spherical and biconical outflows, where the opening angle (θ) and orientation to the line of sight (Ψ) of the biconical outflows can be varied. Opening angles are defined as the full angle for one hemisphere; the maximum opening angle is θ = 180°, corresponding to a sphere. Technical details and equations are presented in Appendix B, but we summarize the basic scheme here. We construct a four-dimensional box (x, y, z, ν), where radii are measured in spherical coordinates from the center of the box. We define the measured cluster radius such that r_SSC = r_half−flux (Table 1). The simulated box is scaled to 4 × r_SSC in each spatial dimension (x, y, z) to fully encompass the line emission (e.g., Figure 4) especially since the geometry along the line of sight (z-axis) is unknown. The velocity axis is scaled based on the input outflow velocity and velocity dispersion, so that the velocity resolution is optimized over the velocities relevant for the cluster and outflow; this is described more fully in Appendix B. Three physical components are required to adequately model the CS 7−6 and H¹³CN 4−3 spectra, as shown in Figures 7 and 8. These components and the input parameters for each are described below. The input parameters are denoted with table footnote "a" in Table 3, which lists input parameter values that yield the best model fits. Radial profiles of these components and the input parameters are summarized in Figure B2 in Appendix B.

• 1.  
  Dust continuum component: Shown in green in Figures 7 and 8, this component is a sphere with r = r_SSC and a constant (in space and frequency) temperature (T_cont). The optical depth in a single cell is a maximum (${\tau }_{\mathrm{cont},\max }$) at the center, then decreases like a Gaussian with FWHM = 2 × r_SSC. In other words, the FWHM of the dust continuum optical depth profile matches the diameter of the 350 GHz continuum source. The temperature and optical depth are set to zero for r > r_SSC.
• 2.  
  Hot gas: Shown in yellow in Figure 7 (and encompassed within the green sphere in Figure 8), this spherical component is required to reproduce the strong emission component of the P-Cygni profiles. This component is defined by an input hot gas temperature (T_hot), a H₂ volume density (n_hot) for the central (r = 0) pixel, and a velocity dispersion (ΔV_hot,FWHM). T_hot is constant (spatially) for r ≤ r_SSC, and is set to zero outside. The density falls off ∝r⁻² from the center and is set to zero for r > r_SSC. The line is centered on zero velocity along the frequency axis, and the Gaussian line width is given by ΔV_hot,FWHM. Using the equations in Appendix B, this produces a spectrum at every pixel in the box.
• 3.  
  Cold, outflowing gas: Shown in blue in Figures 7 and 8, this is the outflow component that produces the absorption features. This component is defined by an input gas temperature (T_out), a H₂ volume density (n_out) at the cluster boundary (r = r_SSC), a constant outflow velocity (V_out), a velocity dispersion (ΔV_out,FWHM), an opening angle (θ), and an orientation to the line of sight (Ψ). The gas temperature is constant (spatially) within this component. The density is a maximum at r = r_SSC and deceases ∝r⁻² until the edges of the box; the density is set to zero inside the cluster (r < r_SSC). In the spectral dimension, the line has a centroid velocity given by V_out and an FWHM line width of ΔV_out,FWHM. Using the equations in Appendix B, this produces a spectrum at every pixel in the box. Since the outflow velocity is constant and the density ∝r⁻², the outflow conserves mass, energy, and momentum. To create a biconical outflow with the input opening angle (θ), the velocity of the pixels outside the outflow cones is set to zero. This creates an ambient gas component (shown in light blue in Figure 8), which has the same temperature, density, and velocity dispersion properties as the outflowing gas (but with V_out = 0). The box is then rotated to the input orientation from the line of sight (Ψ).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Line Profile Modeling Input Parameters and Calculated Values

		SSC 4a		SSC 5a		SSC 14
Quantity	Unit	CS 7 − 6	H13CN 4 − 3	CS 7 − 6	H13CN 4 − 3	CS 7 − 6	H13CN 4 − 3
V${}_{\max \mbox{-} \mathrm{abs}}$ a , b	km s−1	7	7	25	22	25	28
ΔVout,FWHM a , c	km s−1	25	25	20	15	20	25
${V}_{\mathrm{out},\max }$ d	km s−1	28	28	42	35	42	49
log10(tcross) e	log10(yr)	4.9	4.9	4.5	4.6	4.4	4.3
ΔVhot,FWHM a , c	km s−1	10	10	50	50	80	85
Tout a , f	K	7	7	25	20	35	40
Thot a , g	K	900	900	800	800	800	800
Tcont a , h	K	105	105	200	185	340	340
${\mathrm{log}}_{10}($ nout) a , i	log10(cm−3)	7.8	7.5	5.5	5.9	6.2	6.8
${\mathrm{log}}_{10}($ nhot) a , i	log10(cm−3)	8.6	9.9	8.8	9.8	8.8	10.0
${\mathrm{log}}_{10}({\tau }_{\mathrm{cont},\max })$ a , j		−1.3	−1.3	−1.3	−1.3	−1.3	−1.3
${\mathrm{log}}_{10}($ ${N}_{{{\rm{H}}}_{2},\mathrm{out}})$ k	log10(cm−2)	25.4	25.1	23.2	23.7	24.0	24.5
${\mathrm{log}}_{10}({M}_{{{\rm{H}}}_{2},\mathrm{out}})$ k	log10(M⊙)	6.3	6.0	4.4	4.9	5.3	5.4
log10(${\dot{M}}_{{{\rm{H}}}_{2},\mathrm{out}}$) k	log10(M⊙ yr−1)	1.4	1.1	−0.1	0.3	0.9	1.1
log10(tremove‐gas) k	log10(yr)	3.7	4.0	5.4	5.0	4.8	4.6

Notes. The input and output or derived parameters for the best-fitting spherical outflow models. See Section 3.2 and Appendix B for details.

^a Best-fit input parameters to the model. The uncertainties listed below are reflective of the grid of tested parameters. ^b The outflow velocity, which is assumed constant. Uncertainties are ±1 km s⁻¹. ^c The FWHM velocity dispersion of the given component. Uncertainties are ±2 km s⁻¹. ^d The maximum outflow velocity defined as ${V}_{\max \mbox{-} \mathrm{abs}}+2\tfrac{{\rm{\Delta }}{V}_{\mathrm{out},\mathrm{FWHM}}}{2.355}$. Propagated uncertainties are ±2 km s⁻¹. ^e The gas crossing time defined as $\tfrac{{r}_{\mathrm{SSC}}}{{V}_{\max \mbox{-} \mathrm{abs}}}$. Propagated uncertainties are 4% for SSCs 14 and 5 and 14% for SSC 4a. ^f The gas temperature in the outflow at the line peak, which is constant spatially. Uncertainties are ±2 K. ^g The gas temperature in the hot component at the line peak, which is constant spatially. Uncertainties are ±25 K. ^h The continuum temperature, which is constant spatially and spectrally. Uncertainties are ±5 K. ⁱ The log of the peak H₂ number density for the given component. Uncertainties are ±0.25 dex. ^j The peak optical depth of the continuum component. Uncertainties are ±0.025. ^k Uncertainties are ±1 dex due to the molecular abundance ratios relative to H₂.

Download table as:  ASCIITypeset image

Together, these three components are integrated from the back of the box forward (e.g., along the −z-axis in Figures 7 and 8). To obtain the final spectrum, only pixels within a cylinder along the line of sight with r = r_SSC are integrated (shown as the black cylinder in Figure 8) to best compare with the measured spectra, which are extracted only over an area corresponding to the continuum source half-flux radius. We adjust the input parameters component-by-component to find the model spectrum that best matches the observed CS 7−6 and H¹³CN 4−3 spectra for SSCs 4a, 5a, and 14. There are degeneracies among input parameters, which are described below and in Appendix B. These best-fit models are shown in red in Figure 6, and the best-fit parameters for CS 7−6 and H¹³CN 4−3 are listed in Table 3. For all spectra and sources, the spherical model provides the best fit, implying that the opening angles of the outflows need to be broad to explain the observed line profiles. A wide opening angle is in agreement with recent magnetohydrodynamic simulations, which show that cluster outflows are asymmetric and chaotic, but still wide angle in general and regardless of the precise feedback mechanism (e.g., Skinner & Ostriker 2015; Kim et al. 2018; He et al. 2019; Lancaster et al. 2021a, 2021b; Geen et al. 2021). From these best-fit models, we also calculate the H₂ column density and mass in the outflows, which are also listed in Table 3.

We tested models with a fourth physical component representing a fast outflowing component. This was mainly motivated by SSC 14 and the mismatch between the spectrum and model at the blueward edge of the absorption feature for both CS 7–6 and H¹³CN 4−3 (Figure 6). In the model, this component is otherwise identical to the "slow" outflow component described above but with a larger outflow velocity and velocity dispersion and a different maximum H₂ volume. While including this component did marginally improve the fits—especially for SSC 14—the improvement was not enough to justify the additional three parameters introduced into the model.

Our model assumes a constant outflow velocity and r⁻² density profile to conserve mass, energy, and momentum in the outflow. In reality, a constant outflow velocity is unlikely to be precisely the case, due to turbulence within the outflow itself (e.g., Raskutti et al. 2017) and because the outflow may decelerate as it encounters the surrounding medium. From the data, we investigate the location of the absorption trough around and across the sources. We see no strong evidence for systematic velocity shifts around or across SSCs 4a, 5a, or 14.

There are ranges of opening angles (θ) and orientations (Ψ) that are degenerate and will produce similar output spectra. To investigate how well we can constrain θ and Ψ, we run a grid of models with the same input parameters as in Table 2 with varying θ in steps of 20° and Ψ in steps of 5°. The results for CS 7−6 in SSC 14 are displayed in Figure 9, which shows that wide outflows and/or those pointed close to the line of sight are strongly favored. Wide-angle outflows from clusters are in agreement with the results of numerical simulations as mentioned previously. Though the outflows in simulations are clumpy and highly nonuniform, they cover nearly 4π steradians and hence approach the spherical limit of our simple modeling. Narrow and/or off-axis opening angles substantially increase the required input density in our models and result in unphysical solutions for the outflowing mass. In any case, we cannot pin down the precise opening angle and/or line-of-sight orientation from this modeling, but we place limits on them that suggest the dearth of clusters with observed outflows is not a selection bias due to geometrical effects. It is possible that we miss outflows if the dense gas is very optically thick and therefore obscures underlying outflow signatures (e.g., Aalto et al. 2019). We could also miss weak outflows below our detection limit of sensitivity and cluster mass.

Zoom In Zoom Out Reset image size

The two methods of measuring the outflow properties described have calculated quantities in common. In Figure 10, we compare the velocity where the absorption is maximum (${V}_{\max -\mathrm{abs}}$), the FWHM of the absorption feature (ΔV_out,FWHM), the maximum outflow velocity $({V}_{\mathrm{out},\max })$, the gas crossing time (t_cross), the H₂ column density in the outflow (${N}_{{{\rm{H}}}_{2}},\mathrm{out}$), the H₂ mass in the outflow (${M}_{{{\rm{H}}}_{2}},\mathrm{out}$), the mass outflow rate (${\dot{M}}_{{{\rm{H}}}_{2},\mathrm{out}}$), and the time to remove all of the gas mass of the cluster at the current ${\dot{M}}_{{{\rm{H}}}_{2},\mathrm{out}}$ (t_remove‐gas). The equations used to calculate these parameters are given in Appendices A and B for the absorption line fits and line profile modeling, respectively. The vertical axis of Figure 10 shows the logarithm of the ratio of the quantities derived from each method; each major tick mark and horizontal grid line corresponds to an order-of-magnitude difference. In general, there is good agreement between the two methods for the quantities that depend on the velocity of the outflow (${V}_{\max -\mathrm{abs}}$, ΔV_out,FWHM, ${V}_{\mathrm{out},\max }$, and t_cross). The uncertainties in Figure 10 reflect the adopted order-of-magnitude uncertainties around the abundance ratios, as discussed in Section 3.1.

Zoom In Zoom Out Reset image size

Given the uncertainties the agreement between the absorption line fits and the line profile modeling is good, especially for H¹³CN 4–3. A notable exception is for CS 7–6 in SSC 4a, where the line profile modeling yields 1.5 dex larger H₂ columns and masses than derived from the absorption line fits. This is likely because the CS 7–6 absorption in SSC 4a is the most saturated (i.e., the closest to zero). For the absorption line fits, this renders a more uncertain optical depth (from the absorption to continuum ratio) as small changes in the absorption depth leads to large changes in the optical depth and hence the column density. For the line profile modeling, the depth of the absorption depends on the assumed gas temperature and H₂ volume density in the outflow. The bottom of the absorption trough sets an upper limit on the assumed gas temperature in the outflow (T_out); in the modeling, the temperature of the absorption trough cannot be lower than the assumed T_out. In the case of SSC 4a, this temperature is low (∼7 K), well below the temperature needed to excite the J = 7 level of CS of 66 K (Schöier et al. 2005). ¹⁶ As a result, the H₂ density in the outflow must be increased, leading to a large outflowing mass. We assume a constant temperature in the outflowing gas, whereas a temperature gradient is likely needed to produce this absorption depth. Because the other sources and lines have shallower absorption depths, they do not suffer from this effect. Other strong lines—HCN 4–3 and HCO⁺ 4–3—shown in Figure 3 also suffer from this saturation, which is why they are not the focus of this analysis. We, therefore, suggest that the column density, mass, and mass outflow rate from the modeling are overestimated in the case of CS 7–6 in SSC 4a. This effect is not seen in the absorption line fitting because the absorption to continuum ratio is used to infer the optical depth, but an excitation temperature of 130 K is assumed to derive the level populations (Appendix A). In the line profile modeling, the input gas temperature (7 K in the case of CS 7–6 in SSC 4a) is used to derive the level populations instead (Appendix B). Discrepancies between the two methods are also seen in both lines toward SSC 5a, though they are not as extreme as for SSC 4a. The discrepancies in SSC 5a cannot, however, be explained by saturation. The same abundance ratios are used for both methods and these enter into the calculations in the same way, so this discrepancy cannot be remedied by changing [CS]/[H₂].

Given the number of assumptions of the line profile modeling and the parameter covariances, we suggest that the outflow properties from the absorption line fits presented in Table 2 are more reliable and should be adopted. The assumption of spherical outflows in computing those numbers is substantiated by the modeling, which strongly favors wide outflows (Figure 9). It is encouraging that the line profile modeling generally finds similar values for the outflowing mass, but there are cases where the model likely overestimates the outflowing mass (e.g., SSC 4a). Both sets of measurements are limited in the same way by the uncertainty on the molecular abundances with respect the H₂. As discussed in Section 3.1.1, discrepancies between quantities derived from CS 7–6 and H¹³CN 4–3 may be driven primarily by changes in [CS]/[H₂] relative to the assumed abundance, so that quantities derived from H¹³CN 4–3 may be more reliable. Studies, such as the ALMA Comprehensive High-Resolution Extragalactic Molecular Inventory (ALCHEMI) ¹⁷ that measure abundances at higher spatial resolution than currently available in these extreme environments may improve the accuracy of our column density, mass, mass outflow rate, and gas-removal timescale estimates.

There are several timescales and ages calculated from this and previous analyses for these clusters, which we summarize here and list in Table 4:

Table 4. A Comparison of Relevant Cluster Timescales and Ages

SSC Number	${\mathrm{log}}_{10}({t}_{\mathrm{ZAMS} \mbox{-} \mathrm{age}})$ a	${\mathrm{log}}_{10}({t}_{\mathrm{ff}})$ b	tZAMS‐age/tff	${\mathrm{log}}_{10}({t}_{\mathrm{cross}})$ c	${\mathrm{log}}_{10}({t}_{\mathrm{remove} \mbox{-} \mathrm{gas}})$ d	${\mathrm{log}}_{10}({t}_{\mathrm{dep}})$ e
4a	4.95	4.9	1.1	5.0	4.3 ± 0.7	6.6
5a	≳5	4.7	≳2.0	4.6	5.5 ± 0.4	6.4
14	4.88	4.5	2.4	4.4	4.9 ± 0.4	6.8

Notes. For details, see the discussion in Section 4.1.

^a The age of the cluster since the ZAMS (t_ZAMS‐age) calculated by the ratio of the luminosity in protostars to that in ZAMS stars from Rico-Villas et al. (2020). ^b The freefall time (t_ff) calculated by Leroy et al. (2018). For SSC 4a, the value for SSC 4 is used because SSC 4a is the dominant component of the SSC 4 complex. The uncertainty introduced by this assumption is likely minor compared to other uncertainties, which are ∼0.4–0.5 dex (see Table 2 of Leroy et al. 2018). ^c The time for gas to travel from the center of the cluster to the radius of the continuum source (r_SSC) at the typical outflow velocity (${V}_{\max \mbox{-} \mathrm{abs}}$). This timescale (t_cross) is a proxy for the age of the outflow. For each SSC, this is the average of the values in Tables 2 and 3. ^d The time for the entire gas mass of the cluster (from Leroy et al. 2018) to be depleted at the current ${\dot{M}}_{{{\rm{H}}}_{2},\mathrm{outflow}}$. For each SSC, this is the average of the values in Tables 2 and 3. Uncertainties reported in the table are the standard deviations of the mean values from Table 2 and the values in Table 3. Propagated systematic uncertainties (due to the uncertainty in the molecular abundances) are 0.4 dex. ^e The gas depletion timescale, defined as ${t}_{\mathrm{dep}}\equiv \tfrac{{M}_{{H}_{2},\mathrm{tot}}}{{\mathrm{SFR}}_{36\mathrm{GHz}}}$. For SSC 4a, the gas mass for SSC 4 is used because SSC 4a is the dominant component of the SSC 4 complex.

Download table as:  ASCIITypeset image

ZAMS Age ( t_ZAMS‐age )—Rico-Villas et al. (2020) used the ratio of the luminosity in protostars to that of ionizing ZAMS stars to estimate the ages of the clusters (t_ZAMS‐age), finding that SSC 4a is ∼10^4.95 yr old, SSC 5a is ≳10⁵ yr old, and SSC 14 is ∼10^4.88 yr old (Table 4). These measurements are based on the same 36 GHz emission used by Leroy et al. (2018) to calculate the stellar masses. The stellar masses and ages associated with the ionizing photon rates derived from the ZAMS assumption may be underestimated if the cluster stellar population has evolved beyond the ZAMS stage, or if some fraction of the ionizing photons are absorbed by dust. These three SSCs have negligible synchrotron components of their SEDs, so synchrotron contamination of the 36 GHz emission is a small effect. These ages are, therefore, likely lower limits on the true "age" of the cluster.

Cluster Formation Timescales—Given the t_ZAMS‐age and freefall times (t_ff) calculated by Leroy et al. (2018), we can try to place the clusters in a relative evolutionary sequence. An important caveat is that Leroy et al. (2018) estimated t_ff based on the marginally resolved data. With these spatially resolved data, the radii of the clusters decreased (e.g., Figure 11), meaning that these values of t_ff may be overestimated. SSCs 4a, 5a, and 14 have t_ZAMS‐age/t_ff = 1.1, ≳2.0, and 2.4, respectively (Table 4). Modeling by Skinner & Ostriker (2015) suggests that gas is actively collapsing to form stars on timescales ∼1–2t_ff, the typical timescale for cluster formation is ∼5t_ff, and the gas is completely dispersed by ∼8t_ff. These clusters should be nearing the end of the period of active gas collapse. SSC 5a may possibly be transitioning to the initial stages of gas dispersal, especially because it is the only cluster visible in the NIR (Section 4.3), though the gas dispersal has not yet finished because we still see evidence for an outflow. It is important to note, however, that these evolutionary stages are not clear-cut divisions, as gas accretion can continue while the cluster is forming and while outflows are present, and that the t_ZAMS‐age are likely lower limits on the cluster ages. Moreover, the possibility that expelled gas is reaccreted onto the clusters may mean that this cluster formation sequence is more cyclic.

Zoom In Zoom Out Reset image size

Crossing Time (t_cross )—The crossing times (t_cross) we report in Tables 2 and 3 are short: 10^4.5–5.2 yr. Similarly short crossing times are also seen in a SSC candidate in the Large Magellanic Cloud (∼10^4.8 yr; Nayak et al. 2019) and in simulations (e.g., Lancaster et al. 2021a, 2021b). At least one crossing time has passed since the outflows turned on, as P-Cygni line profiles are detected out to ≳ r_SSC. If the outflows are present beyond r_SSC, they are increasingly difficult to detect in absorption away from the continuum source. Therefore, this timescale places a lower limit on the age of the outflow.

Gas-removal Time ( t_remove‐gas )—The gas-removal times (t_remove‐gas) are longer in general than the crossing or freefall times, though the uncertainties on t_remove‐gas are large. The average t_remove‐gas ∼ 10^5.0±0.6 yr, where the uncertainty is the standard deviation of the mean t_remove‐gas for each SSC and line. The gas-removal times are ≳t_ZAMS‐age, except for SSC 4a. Assuming a constant mass outflow rate, this would imply that there is still gas in the clusters to be removed, though a constant mass outflow rate is unlikely (e.g., Kim et al. 2018). This timescale also assumes that none of the expelled gas is reaccreted later on. Given that the bulk of the gas has outflow velocities below the escape velocity (Section 3.1), reaccretion of material is a likely scenario.

Gas Depletion Timescale ( t_dep )—This timescale is the duration of future star formation, assuming a constant star formation rate (SFR) for each cluster and no mass loss: t_dep $\equiv {M}_{{H}_{2},\mathrm{tot}}/\mathrm{SFR}$, where ${M}_{{{\rm{H}}}_{2}},\mathrm{tot}$ is from Leroy et al. (2018). The SFRs we use are also from Leroy et al. (2018) and are based on the measured 36 GHz fluxes that trace the free–free emission from each cluster (Gorski et al. 2017, 2019). This estimate of t_dep based on the SFR assumes continuous star formation (over ∼10 Myr; Murphy et al. 2011), whereas we would expect the actual star formation in these clusters to be bursty. These are the longest timescales for each cluster listed in Table 4. Compared to the clusters' t_ZAMS‐age, this may suggest that the clusters are early in their star formation process and that there is plenty of fuel to form new stars and for the clusters to continue grow. This assumes, however, that all of the molecular gas remains in the cluster. The current gas-removal times of the outflows (t_remove‐gas) are much shorter than t_dep, indicating that these outflows will substantially affect the cluster's SFE. The possibility that gas is reaccreted by the cluster will affect the available gas reservoir for future star formation.

That we detect outflows only in three sources, or ∼8% of the three dozen SSCs in the center of NGC 253 (Figure 1), gives credence to the idea that this outflowing phase must be short-lived. It is unlikely that we miss many sources with outflows due to their orientation and geometry because the modeling presented in Section 3.2 as well as simulations (e.g., Geen et al. 2021) suggest that the outflows are wide. We could, however, be missing outflows if the outer layers of dense gas are very optically thick, which could obscure the outflows (e.g., Aalto et al. 2019) or if there are weak outflows below our sensitivity or cluster mass detection limits. Given that it is expected that the SSCs begin disrupting their natal clouds after ∼10^5–6 yr (Johnson et al. 2015) and taking t_cross = 10^4.5 yr as the lower limit on the age of the outflow, we would expect to find outflows in at least 3%–30% of SSCs, which agrees well with our detection rate of 8%. This percentile range is a lower limit because t_cross is the minimum possible age of the outflow and there could be additional outflows below our detection limit, though they would be weak. This also assumes that the clusters formed at the same time, which is also unlikely.

In general the chemistry-based age sequences presented by Krieger et al. (2020) lead to different relative cluster ages than the dynamical progression presented here, which are also different from the ZAMS age sequence of Rico-Villas et al. (2020). It is important to keep in mind, however, that the oldest clusters are not necessarily the most evolved, and vice versa. Using HCN/HC₃N as a relative age tracer, Krieger et al. (2020) suggested that SSCs 4 and 14 are in the younger half of the SSCs studied, while SSC 5 is among the oldest. An age sequence using the chemistry of sulfur bearing molecules suggests instead that SSCs 5 and 14 are younger, whereas SSC 4 is older (Krieger et al. 2020). This is in disagreement with the age progression suggested by Rico-Villas et al. (2020), who suggested an inside out formation with SSCs 4−12 being the oldest and SSCs 1−3, 13, and 14 being the youngest. The detections of outflows toward SSCs 4a, 5a, and 14 would suggest that they are the most evolved clusters in the young burst, in the simplest model where the clusters completely and finally clear their gas at the end of their formation periods. As described in the following section, SSC 5a may be among the most evolved clusters, as it is the only one of these clusters visible in the NIR. Krieger et al. (2020) also found the lowest dense gas ratios in SSC 5a, suggesting that it has expelled and/or heated and dissociated much of its natal molecular gas. Given the gas-rich environment surrounding these clusters, however, it is possible that other clusters are older and more evolved, but have reaccreted gas from the surrounding medium or that was not completely expelled. Given that the mean velocities of the outflows in SSCs 4a, 5a, and 14 are less than the escape velocities, this is perhaps a likely scenario.

There are a handful of feedback mechanisms relevant for setting the SFE of star clusters. These include proto-stellar outflows, supernovae, photoionization, UV (direct) radiation pressure, dust-reprocessed (indirect) radiation pressure, and stellar winds. Each of these processes is efficient in driving outflows for different cluster masses, radii, and ages. One way to visualize this is through a mass–radius diagram (e.g., Fall et al. 2010; Krumholz et al. 2019), as shown in Figure 11. There is a locus where none of the feedback mechanisms considered by Krumholz et al. (2019) are efficient, and so clusters with those masses and radii should grow with high SFEs. An important caveat of this figure is that there are other parameters relevant to whether a feedback mechanism is efficient, such as the momentum carried by each mechanism and the timescales over which it operates, which are not shown in this representation.

We show in Figure 11 the primary SSCs in NGC 253, where SSCs without outflows are shown as circles and those with outflows are shown as diamonds. The radii are our measured half-flux radii listed in Table 1. Because of the factor of 4 improvement in the linear spatial resolution with respect to the observations reported by Leroy et al. (2018), the radii of these clusters are smaller than reported by Leroy et al. (2018). As a result, the clusters are systematically shifted down in Figure 11 as shown by the vertical dotted gray lines. The stellar masses are taken from Leroy et al. (2018), as described in Section 2.6. We note, in particular, that the reported stellar masses may be underestimated (beyond the error bars) due to uncertainties about dust absorption and evolution beyond the ZAMS.

Many of the SSCs in NGC 253 fall inside the white area of inefficient feedback in Figure 11, which suggests that none of those mechanisms may be efficient for these clusters. The detection of outflows in SSCs 4a, 5a, and 14, however, means there is direct evidence of strong feedback. The H₂ masses in the outflows are significant compared to the total mass in the cluster itself (Table 2). From the previous section (and quantified in Table 4), the gas-removal times (t_remove‐gas) based on the current mass outflow rates are much shorter than the timescale for future star formation (t_dep). Comparison of these timescales implies that the outflows will remove the molecular gas faster than it would otherwise be used up by star formation. Note that these timescales assume either a constant mass outflow rate or a constant SFR, respectively, neither of which is likely to be the case in reality. Moreover, t_remove‐gas assumes there is no new infall replenishing the reservoir, and t_dep assumes no molecular gas is removed. Nevertheless, the outflows will have a non-negligible effect on the reservoir of gas available to the cluster to form stars and hence on the cluster's SFE.

What mechanism, then, is powering these outflows, given that SSCs 4a and 5a lie in the locus where no mechanism is expected to be efficient and SSC 14 is near a boundary? We explore four plausible mechanisms shown in Figure 11 below (excluding proto-stellar outflows that are only important for much lower-mass clusters; e.g., Guszejnov et al. 2021). We also consider winds from high-mass stars that may be important for clusters of these masses and ages (e.g., Gilbert & Graham 2007; Agertz et al. 2013; Geen et al. 2015; Lancaster et al. 2021a, 2021b). In addition to their location in the mass–radius diagram (a reframing of their surface densities), we also compare the momentum expected to be carried by each of these processes and the timescales over which they operate to the values estimated for these clusters. It is also possible that a synergistic combination of mechanisms is at work (e.g., Rahner et al. 2017, 2019). These potential scenarios are discussed in the reminder of this section.

4.2.1. Supernovae

4.2.2. Photoionization

4.2.3. UV (Direct) Radiation Pressure

4.2.4. Dust-reprocessed (Indirect) Radiation Pressure

Given the large dust columns in these SSCs, dust-reprocessed (indirect) radiation pressure is a promising mechanism to power the outflows. In a study of massive star clusters in the Milky Way's Central Molecular Zone, Barnes et al. (2020) found that indirect radiation pressure is important at early times (<1 Myr). Similarly, Olivier et al. (2021) found that dust-reprocessed radiation pressure is the dominant feedback mechanism in ultra-compact H ii regions in the Milky Way. Moreover, given the possible underestimate of the stellar masses due to absorption of ionizing photons by dust or evolution beyond the ZAMS (which are not reflected in the error bars in Figure 11), indirect radiation pressure is a plausible mechanism to drive the observed outflows.

Whether the dust-reprocessed radiation pressure can drive an outflow fundamentally depends on the balance between the outward force of the radiation pressure and the inward force of gravity (i.e., the Eddington ratio, f_Edd). This depends on the dust opacity and the input luminosity from the cluster, which can be parameterized by the light-to-mass ratio of the assumed IMF. The dust opacity (κ_d ) quantifies a dust grain's ability to absorb IR radiation. As explored by Semenov et al. (2003), this quantity varies with temperature, dust composition, grain shape, grain size distribution, and the opacity model. For example, grains in dense molecular clouds exhibit larger κ_d than those in the diffuse ISM, which is thought to be due to the growth of mantles (e.g., Ossenkopf & Henning 1994). Given the high-density and intense radiation environment in these clusters, the dust opacity may be quite different than found in Galactic regions, though we have no precise observational constraints on how much κ_d can vary in these environments. Oftentimes, a dust-to-gas ratio (DGR) is assumed to convert the dust opacity to the opacity in the gas, and the standard Solar Neighborhood value is DGR = 0.01. Finally, the light-to-mass ratio (Ψ) depends on the assumed IMF, where a top-heavy IMF will have a larger value of Ψ. Top-heavy IMFs have been claimed in massive clusters with ${\rm{\Psi }}\sim 2000\,{L}_{\odot }\,{M}_{\odot }^{-1}$ (e.g., Turner et al. 2017; Schneider et al. 2018) compared to either a Kroupa (2001) or Chabrier (2003) IMF, which have ${\rm{\Psi }}=883\,{L}_{\odot }\,{M}_{\odot }^{-1}$ and ${\rm{\Psi }}=1100\,{L}_{\odot }\,{M}_{\odot }^{-1},$ respectively. Models have also found that the IMFs of globular cluster progenitors—thought to be SSCs—are more top heavy with increasing density and decreasing metallicity (Marks et al. 2012). For clusters with the stellar masses and densities of those in NGC 253 (Leroy et al. 2018), we would expect a typical high-mass IMF slope. Toward constraining the nature of the IMF in these clusters, E. A. C. Mills et al. (2021, in preparation) measured the fraction of ionized helium toward these SSCs using H and He RRLs at 5 pc resolution. This ratio is expected to be somewhat dependent on the IMF, as a top-heavy IMF will result in more massive stars capable of ionizing He. E. A. C. Mills et al. (2021, in preparation) measured mass-weighted He⁺ fractions that are consistent with He⁺ fractions found in H ii regions in the center of the Milky Way (e.g., Mezger & Smith 1976; de Pree et al. 1996; Lang et al. 1997), which is not thought to have a globally top-heavy IMF (e.g., Löckmann et al. 2010), though there are individual clusters that do seem to favor a top-heavy IMF (e.g., Lu et al. 2013). While the measured He⁺ fractions do not completely rule out the possibility of a top-heavy IMF (it is unclear precisely how much the He⁺ fraction changes with the IMF), it is unlikely that the IMF in these clusters is extremely top heavy.

Below we describe two approaches taken by simulations in determining the efficiency of dust-reprocessed radiation pressure in driving outflows, and discuss how adjustments to the fiducial assumptions on κ_d , DGR, and Ψ may help explain the outflows observed in SSCs 4a, 5a, and 14.

Numerical simulations of dust-reprocessed radiation pressure by Skinner & Ostriker (2015) assumed that the dust opacity (κ_d ) is constant with temperature (and hence distance from the UV source). Skinner & Ostriker (2015) found

$$\begin{eqnarray}&&{f}_{\mathrm{Edd}}=0.68\left(\displaystyle \frac{\mathrm{DGR}}{0.01}\right)\left(\displaystyle \frac{{\kappa }_{d}}{0.1\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}}\right)\left(\displaystyle \frac{{\rm{\Psi }}}{883\,{L}_{\odot }\,{M}_{\odot }^{-1}}\right),\end{eqnarray} \tag{ 4 }$$

where we have explicitly included the dependence on the DGR. For the fiducial values of the DGR and Ψ, Skinner & Ostriker (2015) found f_Edd > 1 only for κ_d > 0.15 cm² g⁻¹, which is unphysically large for Solar Neighborhood–like dust properties (typically 0.03 cm² g⁻¹; Semenov et al. 2003). If, however, the dust is different in these environments than in the Solar Neighborhood, κ_d could be larger, though we have no observational constraints to evaluate this. Assuming Solar Neighborhood–like dust and a Kroupa (2001) IMF, we can achieve f_Edd > 1 if the DGR is >0.05. Although the center of NGC 253 is known to have a super-solar metallicity (Z = 2.2 Z_⊙; Davis et al. 2013) and the clusters may be even more dust rich (Turner et al. 2015; Consiglio et al. 2016), a DGR > 0.05 seems quite high. Turning instead to the possibility of a top-heavy IMF, Equation (4) yields f_Edd > 1 for ${\rm{\Psi }}\gt 4300\,{L}_{\odot }\,{M}_{\odot }^{-1}$ (for Solar Neighborhood–like values of κ_d and DGR,), which is much more top heavy than has been found in other SSCs (e.g., Turner et al. 2017). Given that E. A. C. Mills et al. (2021, in preparation) did not find strong evidence for an increased He⁺ fraction in these clusters, a very top-heavy IMF is unlikely. Finally, Skinner & Ostriker (2015) noted that for clustered sources of UV photons—almost certainly the case in the SSCs in NGC 253—there can be appreciable cancellation of the radiation pressure terms from each source, so that the net momentum to drive a cluster-scale outflow will be lower, independent of increases in the DGR, κ_d , or Ψ.Therefore, while there may be some combination of increased κ_d , DGR, and Ψ that enables dust-reprocessed radiation pressure to efficiently drive outflows in these clusters, it is unclear how much internal cancellation will affect the net momentum.

Crocker et al. (2018a) took a slightly different approach to study the efficiency of dust-reprocessed radiation pressure. Modeling of the Rosseland mean opacity by Semenov et al. (2003) found κ_d ∝ T² for T < 100 K. Crocker et al. (2018a) converted this temperature dependence into a radial dependence from the central source of UV photons, introducing a gradient in κ_d . This has the effect of boosting the effective dust opacity for clusters that are very dense, allowing for outflows to be driven more easily compared to Skinner & Ostriker (2015). Second, this allows Crocker et al. (2018a) to express their f_Edd in terms of a critical surface density:

$$\begin{eqnarray}&&\,\begin{array}{l}{{\rm{\Sigma }}}_{* ,\mathrm{IR}}=\left(1.3\times {10}^{5}{M}_{\odot }\,{\mathrm{pc}}^{-2}\right)\times \\ {\left(\displaystyle \frac{{\kappa }_{d}}{0.03{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{{\rm{\Psi }}}{1100{L}_{\odot }{{M}_{\odot }}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{\mathrm{DGR}}{0.01}\right)}^{-1},\end{array}\,\end{eqnarray} \tag{ 5 }$$

where we have again explicitly shown the dependence on the DGR. The pink dashed boundary in Figure 11 assumes the fiducial values of κ_d , Ψ, and DGR. While SSC 14 lies on this boundary, SSCs 4a and 5a fall above it, suggesting that dust-reprocessed radiation pressure is not sufficient for driving outflows in these clusters.

Note, however, that there are considerable uncertainties in this picture. The cluster stellar masses may be underestimated if appreciable ionizing photons are absorbed by dust or if the stellar population is evolved beyond the ZAMS. Therefore, they may be closer to the region where dust-reprocessed radiation pressure is more efficient than shown in Figure 11. Moreover, we do not fully understand the properties of dust in these conditions. As discussed above, the DGR in their surrounding gas may be ≳0.022 given the observed super-solar metallicity in the center of NGC 253 (Davis et al. 2013) and the high-density conditions in the starburst molecular gas, which may favor dust formation. Assuming DGR = 0.022 in Equation (5), moves the Σ_*,IR boundary up to encompass SSCs 4a and 5a. This difference compared to the results of Skinner & Ostriker (2015) comes from the assumed temperature dependence in κ_d , which provides a boosted dust opacity for very compact sources. This itself is very uncertain, since the dust models are designed for proto-planetary disks, and the growth of κ_d with T saturates at T ∼ 100 K (see Semenov et al. 2003). Considering changes to the IMF using the fiducial value of κ_d and the DGR would require ${\rm{\Psi }}\gtrsim 3000{L}_{\odot }{M}_{\odot }^{-1}$ to explain SSCs 4a and 5a, or 1.5× more top heavy than the IMF in NGC 5253 (Turner et al. 2017). Even with the boost in κ_d , there would still be cancellations in the radiation pressure due to clustered UV sources, although these cancellations may not be as severe as in the constant κ_d case.

Therefore, it is possible that dust-reprocessed (indirect) radiation pressure could drive the outflows observed in SSCs 4a, 5a, and 14, though this hinges critically on the behavior of the dust opacity (κ_d ) for which there are no observational constraints in these extreme environments. A likely elevated DGR in these sources helps, but cancellations from clustered UV sources hinders the efficiency with which dust-reprocessed radiation pressure can drive outflows. Therefore, whether dust-reprocessed radiation plays a dominant role in powering the outflows observed from SSCs 4a, 5a, and 14 is an open question.

4.2.5. Winds from High-mass Stars

Given the stellar masses of these clusters (M_* = 10^5.0–5.5 M_⊙), we would expect ≳1000–3000 O stars in each cluster, assuming a Kroupa (2001) IMF. It has been suggested that outflows from young, massive SSCs in the Antennae are driven by a combination of O and Wolf–Rayet (WR) stellar winds (Gilbert & Graham 2007), although other possible mechanisms are not evaluated. The combined power of the winds from these massive stars could, therefore, play an important role in powering the outflows from the clusters in NGC 253.

It has been suggested that winds from WR stars significantly impact how a cluster clears its natal gas, as they impart ∼10× the energy of O-star winds over a shorter period of time (e.g., Sokal et al. 2016). The mass-loss rates of WR stars are metallicity dependent, with higher mass-loss rates at higher metallicity for both carbon- and nitrogen-rich WR stars (Vink & de Koter 2005). Because WR stars are evolved O stars, it is thought that they should not contribute much toward the cluster feedback until after ∼3–4 Myr, when other processes such as supernovae are becoming important and when much of the natal gas has already been dispersed. In an observational study of massive embedded clusters, however, Sokal et al. (2016) found that WR winds are important even at earlier stages, though the clusters they studied all have ages >1 Myr. Moreover, they found that clusters without WR stars tend to stay embedded longer than those with WR stars, indicating that WR winds may accelerate the gas-clearing stage.

It is unclear, however, whether the clusters in NGC 253 harbor WR stars yet, given their very young ages (t_ZAMS‐age ≈ 0.1 Myr). There is evidence of a WR population toward SSC 5—perhaps the most evolved cluster—but higher spatial and spectral resolution observations are needed to confirm this (Kornei & McCrady 2009; Davidge 2016). There is a known WR X-ray binary in NGC 253, but it is outside of the nuclear region studied here by ∼250 pc in projection (Maccarone et al. 2014). Given the young ages of these clusters, it is unlikely that there are many, if any, WR stars present in these clusters, especially SSCs 4a and 14. Once a portion of the O star population evolves into WR stars, however, their winds could potentially strongly contribute to driving outflows.

Simulations of stellar wind feedback on the clearing of a cluster's natal gas have mixed conclusions. Some simulations show that stellar winds are important at early times in a cluster's life, especially for clearing the natal dense gas before supernovae start occurring at around 3 Myr (e.g., Agertz et al. 2013; Geen et al. 2015, 2020, 2021). Given that the clusters in NGC 253 are substantially younger than this (t_ZAMS‐age ≈ 0.1 Myr), stellar winds may play a prominent role in driving the outflows we observe from these SSCs. Other simulations, however, have found that stellar winds are most effective after 3 Myr (Calura et al. 2015). There are also simulations that have found that stellar winds (and photoionization) alone cannot expel the natal gas for massive (>5 × 10³ M_⊙) clusters at any time point unless they form with an SFE ≡ M_*/(M_gas + M_*) > 1/3 (Dinnbier & Walch 2020).

In the absence of gas cooling, stellar winds can impart momentum into the surrounding material typically p_r /M_* ≈ 50–65 km s⁻¹ (Weaver et al. 1977), assuming a wind luminosity of 10³⁴ erg s⁻¹ (Starburst99; Leitherer et al. 1999), the ages of these clusters to be ≈ 10⁵ yr (Rico-Villas et al. 2020) and n_H = 10⁵ cm⁻³. These estimates are on the high side of our observed range of p_r /M_* ≈ 5–126, 1–3, and 5–40 km s⁻¹ for SSCs 4a, 5a, and 14, respectively (Table 2), and are fairly insensitive to the assumed average gas density (an order of magnitude in n_H results in a factor of ≈1.5 change in the expected momentum).

If the gas can cool, however, the momentum imparted will be substantially lower, and this is especially relevant in the large ambient densities found near these SSCs (e.g., Silich et al. 2004; Palouš et al. 2014; Lochhaas & Thompson 2017; Wünsch et al. 2017; El-Badry et al. 2019; Gray et al. 2019; Lancaster et al. 2021a, 2021b). An outflow stalled by cooling has been claimed in a SSC in NGC 5253 (Cohen et al. 2018). Nonetheless, recent simulations by Lancaster et al. (2021a, 2021b) found that although cooling is important for the gas densities in the central starburst of NGC 253, the momentum imparted with significant cooling still provides a modest enhancement over a momentum-conserving wind. The computed enhancement factor (α_p ) is sufficient to power outflows even with efficient cooling (with α_p ≈ 1–4 assuming a normal IMF and solar metallicity). For the typical ages of these SSCs (∼0.1 Myr), Lancaster et al. (2021a, 2021b) predicted a momentum injection of p_r /M_* ≈0.8 km s⁻¹ , lower than momenta measured for SSCs 4a, 5a, and 14 of p_r /M_* ≈ 5–126, 1–3, and 5–40 km s⁻¹ (Table 2). Given the predicted shell velocity from Lancaster et al. (2021a, 2021b), we can estimate the value of α_p implied by the observed outflows:

$$\begin{eqnarray}&&{\alpha }_{p}={\left(\displaystyle \frac{{V}_{\max \mbox{-} \mathrm{abs}}}{2.0\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}\left(\displaystyle \frac{{M}_{{{\rm{H}}}_{2}},\mathrm{out}}{{M}_{* }}\right){\left(\displaystyle \frac{{{\rm{r}}}_{\mathrm{SSC}}}{1\mathrm{pc}}\right)}^{-1}.\end{eqnarray} \tag{ 6 }$$

For the measured outflow properties of SSCs 4a, 5a, and 14, this suggests α_p ≈ 9–227, 4–10, and 34–272, respectively (using the parameters in Tables 1 and 2) compared to α_p ≈ 1–4 from the simulations. For all SSCs, our values of M_* may be underestimated (Section 2.6), reducing the values of α_p inferred from Equation (6). A top-heavy IMF could provide up to a factor of 4 enhancement in the strength of the wind (Lancaster et al. 2021a, 2021b). As described in Section 4.2.4, we do not expect the IMF in these clusters to be particularly top heavy (e.g., Marks et al. 2012, E. A. C. Mills et al. 2021, in preparation), although this possibility cannot be ruled out entirely. The wind strength can also be enhanced if the metallicity is super solar, as is likely the case for these clusters (Davis et al. 2013; Turner et al. 2015; Consiglio et al. 2016).

Given the uncertainties in our measured masses, the IMF, and the metallicity, it is possible that the outflow in SSC 5a is powered by stellar winds. Given the concerns about saturation in SSC 4a, the outflowing mass may be overestimated and so our inferred α_p from Equation (6) may also be overestimated, meaning that the outflow from SSC 4a may also be driven (at least in part) by stellar winds. The outflow from SSC 14, on the other hand, is unlikely to be solely by stellar winds.

Therefore, O star winds are unlikely to be driving the outflow observed toward SSC 14 in NGC 253, although it is possible they play a role in driving the current outflows in SSCs 4a and 5a, especially if the current stellar masses of these clusters are underestimated and if the metallicities are super solar. The O star populations in these clusters are likely too young to host many WR stars, except perhaps in the case of SSC 5a, so the effect of WR star winds is likely negligible at this stage in the SSCs evolution.

4.2.6. Summary: What Mechanisms Power the Outflows?

Several clusters and one SSC have been previously identified in the center of NGC 253 based on HST NIR images (Watson et al. 1996; Kornei & McCrady 2009), with the NIR-detected SSC corresponding to SSC 5 (Leroy et al. 2018). To be able to see the NIR emission, the cluster must have dispersed much of its natal gas and dust or it aligns with a serendipitous hole in the extinction screen. If the presence of NIR emission at this cluster location is from gas clearing, this implies that the cluster must be older than the other, highly extincted clusters. As shown in Figure 6, we see evidence for a weak outflow in the dense gas tracers from this cluster, which may be the tail end of this gas dispersal process from the timescale arguments in Section 4.1. The measured mass outflow rate and momentum injection are also lower than for SSCs 4a and 14. It is important to note, however, that SSC 5a still has a high overall gas fraction (${M}_{{{\rm{H}}}_{2},\mathrm{tot}}/{M}_{* }$) of 0.8, though it is not as high as SSCs 4a and 14 (1.3 and 1.6, respectively).

Further support for this picture comes from the discovery of a shell near SSC 5a in HCN 4–3 as shown in Figure 12. This feature is too faint to be seen in CS 7–6 or H¹³CN 4–3 at this resolution. The shell is visible from ∼35–70 km s⁻¹ relative to the cluster systemic velocity (Table 1) and has a projected radius of ∼6 pc, though it is not perfectly centered on SSC 5a. Using these projected velocities, this implies an age of ∼8–16 × 10⁴ years. At its largest extent, the shell reaches the location of SSC 4a (in projection), which also has a dense gas outflow. This shell-like structure is not seen around any other SSCs.

Zoom In Zoom Out Reset image size

Kornei & McCrady (2009) noted that the Pa-α emission around SSC 5a is asymmetric (Figure 12) and could be indicative of an outflow. We investigate how this asymmetric geometry corresponds to the observed shell in HCN 4−3. There are systematic offsets between the location of SSC 5a in the ALMA data, the NIR clusters identified in the HST Near Infrared Camera and Multi-Object Spectrometer (NICMOS) data by Watson et al. (1996) and the F187N and F190N HST NICMOS data used by Kornei & McCrady (2009). Leroy et al. (2018) corrected the positions of the Watson et al. (1996) NIR clusters by Δα, Δδ = +0.32'', −0.5'', so that the NIR cluster corresponding to SSC 5a has α, δ = 00^h47^m32.985ˢ, −25°17^m19.7ˢ (J2000). The positions of the F187N and F190N mosaics ¹⁸ were corrected by matching the locations of the NIR clusters identified by Watson et al. (1996). A linear shift of Δα, Δδ = +1.48'', − 0.85'' brings the HST images into agreement with the NIR clusters and the ALMA data sets. Figure 12 shows the Pa-α (F187N-F190N) image around SSC 5a with the HCN 4−3 contours overlaid. The expanding HCN 4−3 shell is not particularly aligned with the asymmetry seen in the Pa-α emission, though it does seem to align with the north-western edge.

It is therefore possible that this shell is a previous (stronger) version of the outflow detected spectrally in this work. The shell age of ∼10⁵ yr is in good agreement with the minimum ZAMS age and the end of the period of active gas accretion (Section 4.1). It is unlikely that the outflow is due to a supernova explosion both because synchrotron emission is a negligible component of SSC 5's SED (E. A. C. Mills et al. 2021, in preparation) and because the current momentum injection is well below what is expected from supernovae (e.g., Kim & Ostriker 2015; Kim et al. 2017a).