Forming Super Star Clusters in the Central Starburst of NGC 253
Adam K. Leroy, Alberto D. Bolatto, Eve C. Ostriker, Fabian Walter, Mark Gorski, Adam Ginsburg, Nico Krieger, David S. Meier, Elisabeth Mills, Juergen Ott, Erik Rosolowsky, Todd A. Thompson, Sylvain Veilleux, Laura K. Zschaechner
aa r X i v : . [ a s t r o - ph . GA ] O c t Draft version October 26, 2018
Typeset using L A TEX twocolumn style in AASTeX62
Forming Super Star Clusters in the Central Starburst of NGC 253
Adam K. Leroy, Alberto D. Bolatto, Eve C. Ostriker, Fabian Walter, Mark Gorski, Adam Ginsburg,
Nico Krieger, Rebecca C. Levy, David S. Meier,
Elisabeth Mills, J¨urgen Ott, Erik Rosolowsky, Todd A. Thompson,
1, 11
Sylvain Veilleux, and Laura K. Zschaechner Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, Ohio 43210, USA Department of Astronomy, University of Maryland, College Park, Maryland 20742, USA Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544, USA Max-Planck-Institut f¨ur Astronomie, K¨onigstuhl 17, D-69117, Heidelberg, Germany Department of Physics and Astronomy, University of Western Ontario, London, Ontario N6A 3K7, Canada National Radio Astronomy Observatory, PO Box O, 1003 Lopezville Road, Socorro, New Mexico 87801, USA Jansky Fellow New Mexico Institute of Mining & Technology, 801 Leroy Place, Socorro, New Mexico 87801, USA Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, Massachusetts 02215, USA Department of Physics, University of Alberta, Edmonton, AB T6G 2E1, Canada Center for Cosmology & Astro-Particle Physics, The Ohio State University, Columbus, Ohio 43210 University of Helsinki, Physicum, Helsingin Yliopisto, Gustaf H¨allstr¨omin katu 2, 00560 Helsinki, Finland Finnish Center for Astronomy with ESO
Submitted to ApJABSTRACTNGC 253 hosts the nearest nuclear starburst. Previous observations show a region rich in moleculargas, with dense clouds associated with recent star formation. We used ALMA to image the 350 GHzdust continuum and molecular line emission from this region at 2 pc resolution. Our observations reveal ∼
14 bright, compact ( ∼ − Hubble near-infrared (IR)imaging. Our observations imply that gas still constitutes a large fraction of the overall mass in thesesources. Their high brightness temperature at 350 GHz also implies a large optical depth near the peakof the IR spectral energy distribution. As a result, these sources may have large IR photospheres andthe IR radiation force likely exceeds
L/c . Still, their moderate observed velocity dispersions suggestthat feedback from radiation, winds, and supernovae are not yet disrupting most sources. This modeof star formation appears to produce a large fraction of stars in the burst. We argue for a scenarioin which this phase lasts ∼ INTRODUCTIONVigorous bursts of star formation in galaxy cen-ters and merging galaxies produce “super” star clus-ters (SSCs, e.g., Holtzman et al. 1992; Whitmore 2003;Portegies Zwart et al. 2010). The SSCs in well-knownlocal starbursts like M82 and the Antennae galaxieshave been studied for decades (e.g., Whitmore 2003;McCrady et al. 2005). These massive (M ⋆ > M ⊙ ),compact ( R ∼ M ⋆ & M ⊙ , age .
100 Myr) have been found in theMilky Way and many nearby galaxies and (see re-view by Portegies Zwart et al. 2010; Longmore et al. 2014). Overall, the fraction of stars formed in clus-ters appears to increase with the surface density ofstar formation (Kruijssen 2012; Johnson et al. 2016;Ginsburg & Kruijssen 2018). Because higher levels ofstar formation activity were prevalent at z ∼ − Leroy, Bolatto et al. the Large Magellanic Cloud (Ochsendorf et al. 2017).But despite decades of optical and near infrared stud-ies, only a pair of candidate forming
SSCs have been re-solved in cold gas and dust emission(Turner et al. 2017;Oey et al. 2017). In both cases, CO (3-2) emission hasbeen seen associated with an SSC in a starburst dwarfgalaxy. This CO (3-2) emission may trace moderatelymore excited and dense gas than the CO (1-0) line.In the Milky Way, ∼ . M ⊙ , somewhat lower than the extragalacticproto-SSC candidates. They appear likely to form M ⋆ ∼ × M ⊙ clusters, assuming ∼
30% efficiency(see Bressert et al. 2012).In this paper, we report the identification of 14 can-didate proto-SSCs in NGC 253. This galaxy hostsone of the nearest nuclear starbursts ( d ∼ . ∼ ⊙ yr − (e.g., Leroy et al. 2015; Bendo et al.2015). This burst is fed by the galaxy’s strong bar(Sorai et al. 2000), making NGC 253 a prototype forthe common phenomenon of bar-fed nuclear starbursts(see Kormendy & Kennicutt 2004).Watson et al. (1996) and Kornei & McCrady (2009)used the Hubble
Space Telescope to discover a young( ∼ A V ∼
17 mag) SSC inthe nuclear region of NGC 253. Any other SSCs in thenuclear region must be too heavily embedded to appearprominent in
Hubble images, including the near-infrared(IR) images presented in Walter et al. (2017). ButUlvestad & Antonucci (1997), following Turner & Ho(1985), showed the presence of ∼
30 flat spectrum, com-pact ( ∼ Hii regions. One of these coincides with the SSC ofWatson et al. (1996) and Kornei & McCrady (2009).Previous mm- and submm-wave observations showthat the NGC 253 nuclear region hosts massive, densemolecular clouds (Sakamoto et al. 2011; Leroy et al.2015). The whole region resembles a heavily scaledup version of the Milky Way’s Central Molecular Zone(Sakamoto et al. 2011). Observations at θ ≈ . ′′ res-olution show that these clouds host compact <
10 pcsized clumps of gas and dust, which have the appropri-ate masses to form massive clusters and are associatedwith signatures of massive star formation (Ando et al.2017). Given these candidate SSC-forming structuresand the presence of at least one known SSC, NGC 253is the ideal target to catch SSC formation in the act.In this paper, we use the Atacama Large submillime-ter/Millimeter Array (ALMA) to map dust emissionfrom the NGC 253 starburst at 0 . ′′ ≈ . OBSERVATIONSWe used ALMA to observe NGC 253 at ν ∼
350 GHz( λ ∼ µ m) as part of project 2015.1.00274.S (P.I. A.Bolatto). We observed with the main array in both anintermediate configuration and a 2-km extended configu-ration. We also used 7-m array in the Atacama CompactArray (ACA) to recover short spacing information. Thebandpass captures the sub-mm continuum from dustemission and covers several molecular rotational tran-sitions, including CO (3-2), HCN (4-3), HCO + (4-3),CS (7-6), and H CN (4-3). The full suite of line im-ages and the extended molecular gas distribution andkinematics traced by CO (3-2) will be presented by N.Krieger et al. (in preparation).We combined the observatory-provided calibrated vis-ibilities for two 12-m configurations and the ALMACompact Array (ACA) 7-m array and imaged them inversion 5 . . tclean task. The in-clusion of the ACA means that scales out to 19 ′′ arerecovered.We are interested in the compact structures at theheart of the starburst. Therefore, when imaging thecontinuum, we adopted a Briggs robust parameter r = − CN (4-3), sensitivity remainsa concern. Therefore in the line images we emphasizedsurface brightness sensitivity and used a standard Briggsweighting with robust parameter r = 0 . ν = 350 GHz and the finalFWHM beam size is θ = 0 . ′′ . Before the convolutionto a round beam, the beam of the continuum image hasmajor and minor FWHM 0 . ′′ × . ′′ .The rms noise away from the source in the cleaned,round beam image is 0 . ≈ . − .For the H CN (4-3) and CS (7-6) line images the beamsize is 0 . ′′ (convolved from ∼ . ′′ × . ′′ ) and thetypical rms in the cube is 0 . − channel.The ancillary CO (3-2) and HCN (4-3) observations havesimilar resolution and noise. More details of the lineimaging appear in N. Krieger et al. (in preparation).We compare the ALMA data to Karl G. JanskyVLA imaging of ν = 36 GHz continuum emission(Gorski et al. 2017, and M. Gorski et al. in preparation).At this frequency and resolution, the radio continuumis likely to be predominantly free-free emission (e.g., seeMurphy et al. 2011). These data have native resolution orming Super Star Clusters in NGC 253 Figure 1. ( top left ) LVL IRAC 8 µ m image of NGC 253 (Dale et al. 2009; Lee et al. 2009). The square shows the regionhighlighted in the bottom left panel. ( bottom left ) The field analyzed in this paper (square box) plotted over the centralpart of the 8 µ m image. We focus on the innermost region of the galaxy. This region hosts ∼
10 dense molecular clouds(e.g., Sakamoto et al. 2011; Leroy et al. 2015; Ando et al. 2017), a large amount of high density gas (e.g., Paglione et al. 2001;Knudsen et al. 2007; Leroy et al. 2015; Meier et al. 2015), and forms stars at a rate of ≈ ⊙ yr − (see Bendo et al. 2015;Leroy et al. 2015). ( right ) ALMA ν = 350 GHz continuum emission from the inner region of NGC 253 at 0 . ′′ ∼ . . ν = 350 GHz ( ∼ µ m), kinematics, and associationwith dense, excited gas suggest that many of these peaks may represent forming super star clusters (see Sections 3 and 4). Notethat we show a 10 ′′ × ′′ cutout covering the region of interest, but that the FWHM of the primary beam of ALMA’s 12 mantennas is 18 ′′ at 350 GHz. Contours in the continuum image show 0.6 K (gray), and 1, 2, 4, ... K (black). slightly better than the ALMA continuum image, witha FWHM beam of 0 . ′′ × . ′′ . We convolve them tothe match the resolution of ALMA, θ = 0 . ′′ , for anal-ysis. After convolution the VLA data have rms noise ∼ .
03 mJy beam − .We also compare to Hubble
Space Telescope imagingof the near-IR ( λ = 1 . µ m) continuum. These were ob-tained to serve as an off-line continuum measurement forthe Paschen β observations presented by Walter et al.(2017). CANDIDATE FORMING SUPER STARCLUSTERSThe top left panel in Figure 1 shows whole disk ofNGC 253 seen at 8 µ m by the Local Volume Legacy sur-vey (LVL Dale et al. 2009; Lee et al. 2009). The 8 µ m image shows the location of UV-heated small dust grains(likely polycyclic aromatic hydrocarbons (PAHs)), andso illustrates the overall morphology of the ISM in thegalaxy. The bottom left panel zooms in on the squareregion indicated in the top panel. The square in the bot-tom left panel shows our regions of interest in this paper.This is a square field, 10 ′′ × ′′ across that includes mostof the active star formation and dense clumps in the nu-clear starburst. ALMA’s 12 m antennas have a primarybeam of 18 ′′ at 350 GHz, so the ALMA observationscover a somewhat larger field of view than we show inthe Figure. In total, the nuclear burst in NGC 253 con-tains ∼ × M ⊙ of molecular material and formsstars at ∼ ⊙ yr − (Leroy et al. 2015; Bendo et al.2015). Leroy, Bolatto et al.
Figure 2.
Our region of interest in gas, radio continuum, and near-IR emission. The region is rich in molecular gas andsignatures of recent star formation. Panel (a) shows a sprawling, high column density distribution of CO (3-2) emission (N.Krieger et al. in preparation). However, the CO integrated intensity on its own only roughly suggests the dense peaks seenin the continuum images. On the other hand, panel (b) shows that the continuum sources do coincide with peaks of emissionfrom CS (7-6) and H CN (4-3), tracers of dense, excited molecular gas. Panel (c) shows that the continuum peaks are alsomostly coincident with bright ν = 36 GHz radio continuum emission. Radio continuum at this frequency and resolution is mostlikely to be free-free, tracing ionizing photon production by young stars. Despite the prodigious concentrations of gas and thelikely presence of embedded massive stars, these sources are mostly missing from the HST near-IR continuum image shown inpanel (d) . In that panel, we mark the four clusters identified from HST imaging by Watson et al. (1996). Among our candidateforming clusters, only source . µ m. It coincides with the previously known SSC from Watson et al.(1996) and Kornei & McCrady (2009). Contour levels: (a) CO (3-2) integrated intensity at 500 K km s − (gray), then 1000,2000, ... K km s − (black); (b) H CN integrated intensity 50, 100, 200, ... (black); (c) ALMA 350 GHz continuum image (i.e.,Figure 1) at 1, 2, 4, ... K (gray); (c) ALMA 350 GHz continuum image (i.e., Figure 1) at 1, 2, 4, ... K (black). orming Super Star Clusters in NGC 253
Dust Continuum and Gas
The right panel of Figure 1 shows ν = 350 GHz( λ ∼ µ m) continuum emission from this inner regionat θ = 0 . ′′ ≈ . T b ∼ . − σ above any contourshared with another peak, or 5 σ above 0 K if there is noshared contour.These peaks have brightness of a few K up to a fewtens of K and FWHM sizes of ∼ . − ∼ − M ⊙ suggest that these structures are formingSSCs.The bright peaks are still associated with large sur-rounding reservoirs of gas. We show this in the topleft panel of Figure 2, where we plot the line-integratedCO (3-2) intensity. The peaks sit at the hearts of themassive clouds and clumps studied by Sakamoto et al.(2011), Leroy et al. (2015), and Ando et al. (2017).They are not conspicuous in the integrated CO (3-2)intensity, although N. Krieger et al. (in preparation)show that they can be identified from the CO kinemat-ically.The continuum peaks stand out in lines that trace highdensity molecular gas. The top right panel of Figure 2shows our region of interest in line-integrated CS (7-6)as a color image with contours showing line-integratedH CN (4-3) intensity. The H CN (4-3) line emits mosteffectively at densities n & cm − (Shirley 2015) and T &
40 K. The CS (7-6) emission, which also traceswarm, dense gas, has critical density ∼ × cm − and requires T &
60 K.At coarser resolution, HCN (4-3) and CS (7-6) emis-sion correlate with IR emission in star-forming galaxies,with a linear relationship relating IR and line luminos-ity (Zhang et al. 2014; Tan et al. 2018). Here we seeH CN (4-3) and CS (7-6) emission directly associatedwith the sites of massive star and cluster formation on ∼ Signatures of Massive Star Formation
Are there actually signatures of young, massive starsassociated with these peaks of gas and dust emission?The bottom left panel of Figure 2 shows ν ∼
36 GHzcontinuum emission from our target field (Gorski et al.2017, and M. Gorski et al. in preparation) with thedust continuum contours overlaid. Our dust continuumpeaks are coincident with, or near to, peaks of brightradio continuum emission.At this frequency and resolution, most of the sourcesin the ν = 36 GHz map arise from free-free emission.Modulo loss of ionizing photons to dust and a mild de-pendence on the electron temperature and Gaunt factor,free-free emission directly traces ionizing photon produc-tion in a manner similar to optical recombination lines.But unlike optical and near-IR recombination line emis-sion, 36 GHz emission is almost totally unaffected byextinction. Thus the bottom left panel of Figure 2 sug-gests that young, heavily embedded massive stars lie ator near most of our observed dust clumps.The bottom right panel of Figure 2 shows that thesesignatures of massive star formation are almost totallyobscured by dust even in the near infrared. We plotthe near infrared (1 . µ m) continuum as seen by Hub-ble ’s Wide Field Camera 3 (filter F130N). We also in-dicate the position of the four clusters identified byWatson et al. (1996) from earlier
Hubble near-IR imag-ing. To match the astrometry of our near-IR data, whichaligns well with the ALMA and VLA observations, wefound it necessary to shift the measured positions fromWatson et al. (1996) by ∆ α, ∆ δ ≈ +0 . ′′ , − . ′′ .The image shows bright stellar continuum emis-sion coincident with the brightest SSC known fromWatson et al. (1996) and Kornei & McCrady (2009).Otherwise our dust peaks do not correspond to clearenhancements in the near-IR continuum. Given thepresence of free-free emission, these sources are likely tobe bright, compact, massive collections of young stars.But the extinction in the inner region of the galaxy istoo severe to pick them out even in the near-infrared.This overwhelming extinction is striking, but not sur-prising. From the top left panel in Figure 2, we seethat our peaks all lie at I CO − & ,
000 K km s − .Under the conservative assumptions of thermalized COlines, a low α CO = 0 . ⊙ pc − (K km s − ) − , and aGalactic dust-to-gas ratio, this amount of gas still cor-responds to E ( B − V ) ∼
34 mag, or A J ∼
30 mag. Evenwithout accounting for the dense concentrations withinthe clouds, the central region of NGC 253 is heavilyextinguished and capable of hiding luminous clusters atnear-IR wavelengths.
Leroy, Bolatto et al.
This 36 GHz view of the NGC 253 nucleus resem-bles the ∼
43 GHz, 3 pc resolution view of M82 byTsai et al. (2009). In M82, another starburst at d ∼ . ∼ Hii regionspowered by massive clusters. Tsai et al. 2009 showedthese compact
Hii regions to exist within dense gasstructures observed at ∼
45 pc resolution. Based onour observations of NGC 253, it seems plausible, evenlikely, that some of the individual Tsai et al. (2009)sources will still be in the process of formation andthat high resolution sub-mm observations of M82 wouldshow associated pc-scale concentrations of gas and dust.Schinnerer et al. (2007) observed similar coincidence ofdense gas tracers and continuum signatures of em-bedded massive star formation at ∼
10 pc resolutionin the inner region of NGC 6946, though there weresome detailed differences between their HCN map andthe NGC 6946 continuum emission seen by Tsai et al.(2006). Turner et al. (2003) and Turner & Beck (2004)found similar compact
Hii regions surrounding the SSCin NGC 5253. Turner et al. (2017) showed CO emissionfrom the same source, though that emission appears op-tically thin in CO, perhaps indicating that the NGC5253 cluster is at a later evolutionary stage than theones that we observe. PROPERTIES OF THE CANDIDATE FORMINGSUPER STAR CLUSTERSWe estimate the size, line width, and fluxes associatedwith each peak and use these to gauge the masses of thecandidate proto-SSCs in several ways.4.1.
Size, Line Width, and Flux Measurements
Size, Peak Temperature, and Flux at
GHz:
To measure the sizes associated with each peak, we buildan azimuthally averaged radial profile centered on thepeak. In each half-beam thick ring centered on the peak,we calculate the median intensity. Using the mediansuppresses the influence of nearby peaks and the brightsurrounding filamentary features, and so emphasizes theprofile of the central peak. We further reject 3 σ outliersabout this median profile. Figure 3 shows the resultingprofiles appear as black points, with error bars showingthe scatter about the profile. Blue lines show a Gaussianfit to these profiles; this fit includes a background term,which is small in all cases.To compute deconvolved sizes we subtract the beamsize in quadrature from the Gaussian fit to the profile.We also correct the peak temperature for the effects ofthe beam by scaling the measured peak temperature bythe ratio of measured source area to the deconvolvedsource area. The deconvolved profiles appear as red linesin Figure 3.We report the measured sizes in Table 1. As a check,we also fit two dimensional Gaussians to each source.The deconvolved FWHM size from the Gaussian fits agree with our measured sizes with a scatter of ± . T peak to be the sum in quadrature of the fractional un-certainty due to statistical noise and the fractional un-certainty in the deconvolved beam area.From the fits to the profiles, we also calculate the fluxof each source at ν = 350 GHz, which we give in Table1. The statistical uncertainties on this flux are low, be-cause the peaks are all detected at high signal to noise.In this case, we quote a 10% uncertainty on the over-all flux for each source, reflecting a mixture of calibra-tion uncertainty (which should be covariant among allsources), uncertainty in the image reconstruction, anduncertainty due to the adopted methodology.For reference, from a lower resolution, robust-weightedversion of the continuum map, we calculate a total350 GHz flux of ≈ . S/N = 3. The sources in Table 1 have total flux0 .
28 Jy, and so account for ∼
15% of the total 350 GHzemission from the nuclear region.
Line Widths:
We measure the line width associatedwith each peak. To do this, we define a series of aper-tures. The aperture associated with a peak has radiusequal to the FWHM fit (not deconvolved) size of thepeak and sits centered on the peak. The backgroundregion associated with each aperture extends from ra-dius 1 to 3 times the FWHM fit size of the centralsource. The background excludes apertures associatedwith other peaks. We calculate the source spectrumby subtracting the average spectrum in the backgroundfrom the average spectrum in the aperture. Note thatthe central apertures that we use are never less than0 . ′′ across (diameter). Thus the measurement regionis always at least moderately extended compared to the0 . ′′ beam of the H CN (4-3) and CS (7-6) cubes.We fit a Gaussian to the background-subtracted forCS (7-6) and H CN (4-3) spectra, fitting over a velocityrange picked by eye to cover the emission line. Figure 3shows both background-subtracted spectra and the fitsfor each peak. We take the linear average of the two linewidths as the characteristic line width for the source. Weadopt one half the difference in the line width derivedbetween the two lines as our uncertainty.
36 GHz fluxes:
We measure fluxes for each sourcefrom the VLA 36 GHz map. To do this, we subtract theaverage intensity in the local background region fromthe region inside the aperture. Then we sum the fluxinside the aperture. We use the same apertures used toderive the line widths.Similar to the case of the ALMA fluxes, the statisticaluncertainty in the 36 GHz flux is small (compare thefluxes in Table 1 to the 0.02 mJy beam − noise). Weadopt an uncertainty of 10% to reflect calibration andimage reconstruction uncertainties. orming Super Star Clusters in NGC 253 Figure 3.
Spatial ( left ) and spectral ( right ) profiles of our 14 peaks. The left column shows the binned median-based radialprofile of 350 GHz emission about each peak (black bins, with a blue Gaussian fit). The black, shaded profile in each panelindicates the shape of the synthesized beam. The red profile shows the inferred shape of the peak after deconvolving the beam.The spectra show the background subtracted CS (7-6) emission in red and H CN (4-3) emission in blue. Lines indicate Gaussianfits to the profiles. Spectrum
Leroy, Bolatto et al.
Table 1.
Measured Properties Candidate Young Clusters in NGC 253 T pk FWHM σ v F F , appa F ,app a f ff a ( ◦ ) ( ◦ ) (K) (pc) (km s − ) (mJy) (mJy) (mJy)1 11.886669 -25.289232 6.1 ± . ± . ± . ± . ± .
45 0.51 ± .
05 0.03 ± .
012 11.886749 -25.289240 23.5 ± . ± . ± . ± . ± .
21 0.55 ± .
06 0.04 ± .
013 11.886829 -25.289200 11.1 ± . ± . ± . ± . ± .
44 0.32 ± .
03 0.01 ± .
014 11.887265 -25.288950 14.2 ± . ± . ± . ± . ± .
15 2.57 ± .
26 0.07 ± .
025 11.887444 -25.288816 29.8 ± . ± . ± . ± . ± .
28 6.86 ± .
69 0.13 ± .
026 11.887542 -25.288733 4.2 ± . ± . ± . ± . ± .
62 4.91 ± .
49 0.63 ± .
407 11.887579 -25.288628 2.5 ± . ± . ± . ± . ± .
06 0.87 ± .
09 0.07 ± .
058 11.887978 -25.288244 22.7 ± . ± . ± . ± . ± .
85 1.65 ± .
16 0.05 ± .
029 11.887984 -25.288389 8.5 ± . ± . ± . ± . ± .
98 8.02 ± .
80 0.32 ± . ± . ± . ± . ± . ± .
66 4.85 ± .
48 0.11 ± . ± . ± . ± . ± . ± .
89 9.81 ± .
98 0.41 ± . ± . ± . ± . ± . ± .
00 26.72 ± .
67 0.71 ± . ± . ± . ± . ± . ± .
28 1.66 ± .
17 0.04 ± . ± . ± . ± . ± . ± .
50 7.43 ± .
74 0.07 ± . a Values measured in apertures centered on the peaks. The apertures have radius equal to the FWHM size of the sourcebefore deconvolution (i.e., to recover this add 1 . Note —R.A., Dec. refer to peak position in the 350 GHz continuum image. FWHM size assumes a distance of 3 . . ′′ ≈ . T pk reports the peak intensity, after deconvolvingthe beam, at 350 GHz in Rayleigh Jeans brightness temperature units. σ v is the linear average of the CS (7-6) andH CN (4-3) rms line width. F is the flux at 350 GHz estimated from the Gaussian fit to the profile. F is the fluxof 36 GHz emission obtained from aperture photometry, this is scaled into a luminosity using the distance and used toestimate Q and M ⋆ using the equations in the text. f ff refers to the estimated fractional free-free contribution to the350 GHz emission based on comparing fluxes measured in matched apertures. orming Super Star Clusters in NGC 253 Figure 4.
Profile similar to those in Figure 3 but forthe known Milky Way protoclusters in Sgr B2 (e.g., seeGinsburg et al. 2018). We show the profile of the sourcesas seen by ATLASGAL (see Urquhart et al. 2018) matchedto our 1 . N (24-23). In both cases, we have plottedthe data similar to how we show our profiles, removing alocal background from the radial profile and subtracting thecontinuum from the spectrum. The combined Sgr B2 pro-toclusters show peak continuum brightness ∼ ∼ . σ line width ∼ . − . Thus, theywould appear as among the least bright and least compactof our sources, with the narrowest line widths. But they dooverall resemble the sources that we see in NGC 253, so thatour observations pick out sources that resemble scaled upversions of known protoclusters in the Milky Way. Similarto Sgr B2, we might expect some of our sources to break upinto two or more protoclusters within our 1 . Based on the 36 GHz emission, we estimate the frac-tional contribution of free-free flux to the ALMA bandvia: f ff ≈ (cid:18) (cid:19) . F F (1)where the first factor reflects the expected − . F and F refer to the observed total flux at 36 GHz and350 GHz. For this application only, we measure fluxesfrom the ALMA 350 GHz map in exactly the same waythat we measure the 36 GHz fluxes (i.e., using aperturephotometry and the same aperture definitions). We re-port both sets of fluxes in Table 1. f ff is the fraction of the 350 GHz flux in the aperturethat can be attributed to free-free emission, assumingthat all of the 36 GHz flux comes from optically thinfree-free emission. A value ≪ ∼ f ff around four sources: peaks f ff = 0 . . . . ∼ −
200 GHz and at ∼ −
25 GHz will helpresolve the nature of the emission (some observationsat slightly coarser resolution already exist Mohan et al.2005; Ulvestad & Antonucci 1997).Peak ∼
10 pc away mightcontribute to ionization in the region. In any case, weapply f ff as a correction to the gas mass estimates, view-ing this as the most conservative option.4.2. Resemblance to a Known Milky Way Protocluster
As a check, we construct profiles similar to those inFigure 3 for the known Galactic protoclusters Sgr B2.This pair of bright sources near the Galactic center is re-garded as very likely to be forming young massive clus-ters (see Ginsburg et al. 2018; Urquhart et al. 2018, andreferences therein). We degrade the ATLASGAL 500 µ mdata to a resolution of 1 . β = 2). We also extract a spectrum ofHC N (24-23) at 3 pc (FWHM) resolution to serve asa proxy for our H CN and CS measurements, thoughnote that HC N (24-23) has larger excitation require-ments than the H CN (4-3) or CS (7-6).The resulting profile and spectrum, shown in Figure4, show that Sgr B2 would have a slightly larger sizeand narrower line width than our candidate clusters. Itwould also have among the lowest brightness tempera-tures of our sources. But overall, the structure in Figure4 does resemble what we see for our NGC 253 sources(Figure 3). The comparison gives us confidence that wedetect moderately more compact, scaled up versions ofa known Galactic protocluster.Sgr B2 appears as a single extended source in this ex-ercise, but also breaks into two massive protoclusters athigher resolution (e.g., see Figure 1 in Ginsburg et al.2018). Therefore, this comparison also highlights thelikelihood that despite our high (for extragalactic work)1 . Leroy, Bolatto et al. observed at higher resolution. From a first look at ∼ Gas, Stellar, and Dynamical Masses
Based on their size, line width, and fluxes, we estimatethe gas, stellar, and dynamical masses of these clustercandidates. We report these in Table 2. To do this, weassume that (1) the free-free corrected 350 GHz emis-sion arises from dust with some adopted temperatureand emissivity, which we take to be well mixed withthe gas with some characteristic dust-to-gas ratio; (2)the 36 GHz represents free-free emission emitted by ayoung stellar population on the zero age main sequence(ZAMS); and (3) our objects are in virial equilibrium, sothat their sizes and line widths together indicate theirdynamical mass.4.3.1.
Zero Age Main Sequence Stellar Mass
Assuming that all of the 36 GHz emission is producedby free-free interactions, we can estimate the ionizingphoton production rate of each source. From this, wecan calculate the mass of ZAMS stars needed to producethis number of ionizing photons.Following Murphy et al. (2011) and Caplan & Deharveng(1986), a 36 GHz luminosity, L , implies an ionizingphoton production rate of: Q ∼ . × L s − . (2)Where we have assumed an electron temperature T e =7 ,
000 K (slightly higher than the estimate for NGC 253by Bendo et al. 2015). Here Q is the ionizing pho-ton production rate per second and L is measured inerg s − Hz − .Based on Starburst99 calculations (Leitherer et al.1999), for a ZAMS population the mass ( M ⋆ ) relatesto the ionizing photon production rate ( Q ) via M ⋆ ∼ Q × M ⊙ . (3)We arrive at this value by simulating 10 M ⊙ single stel-lar population, with the initial mass function of Kroupa(2001), a maximum stellar mass of 100 M ⊙ , and thedefault stellar evolution tracks and tuning parameters.Then we divide the ionizing photon output at time zeroby mass of the stellar population. Together, Equations3 and 2 yield the mass of the embedded stellar popu-lation needed to produce the observed 36 GHz flux viafree-free emission. Our candidate proto-SSCs have median log M ⋆ [ M ⊙ ] =5 . M ⋆ [ M ⊙ ] ∼ . − .
0. Basedon this calculation, all of our sources already qual-ify as young massive clusters (i.e., M ⋆ & M ⊙ Portegies Zwart et al. 2010; Longmore et al. 2014).Note that the two highest values, for peaks f ff found for these objects (Table 1) and the uncer-tain nature of the ionization in this region (see above).Our M ⋆ is a linear translation of the 36 GHz flux, andour sources are detected at high signal-to-noise. Thisyields small statistical uncertainties, ∼ ± M ⋆ . First,uncertainties in the Gaunt factor imply a systematic un-certainty of ≈ ±
20% (Rybicki & Lightman 1986). Sec-ond, if dust absorbs a significant fraction of the ioniz-ing photons, our M ⋆ will represent an underestimate.Loss of ionizing photons to dust already appears to be asignificant effect in massive star forming regions in theMilky Way (Binder & Povich 2018), so we do expectthis to also be important in the dustier, dense nucleusof NGC 253. But the magnitude of the effect is notclear. More, if there is ongoing accretion and many ofthe stars in the cluster have not yet reached the mainsequence, we would also expect a higher mass per ion-izing photon produced. Most of these effects have thesense that our quoted M ⋆ likely represents a moderateunderestimate. The calculation also has the usual un-certainties related to the upper mass cutoff and shape ofthe IMF. Our estimates also take no account of the in-fluence of binary stars (e.g., Eldridge et al. 2017). Notethat, following Xiao et al. (2018), we expect the inclu-sion of binaries to have the largest effect after a few Myr.Thus we expect binarity to be a minor concern for thispaper, which focuses on young sources. Finally, notethat we estimate ZAMS mass in an aperture centeredon the cluster. Future high resolution comparison of the36 GHz structure, dust, and gas will help us understandhow much of this mass is, in fact, directly associatedwith the gas and dust peaks.4.3.2. Implications for the Ionized Gas Content
From Q and a plausible size, we estimate the ionizedgas mass and density associated with each source. Weposit an Hii region with radius r S at the heart of eachsource. Assuming Case B recombination, complete ion-ization of hydrogen, and 1 .
36 contribution of helium bymass, we expect M ion ≈ (cid:18) Q α B / π r S (cid:19) . . m H (4) M ion ≈ M ⊙ (cid:18) Q s − (cid:19) . (cid:18) r S (cid:19) . . Here α B ≈ . × − cm s − is the adopted recom-bination rate coefficient corresponding to an Hii region orming Super Star Clusters in NGC 253 Table 2.
Estimated Physical Properties Candidate Young Clusters in NGC 253 M VTa log M gasb log M ⋆ c log Σ gas+ ⋆ b,c log ρ gas+ ⋆ b,c log τ ffb,c p r /M ⋆ v esc (M ⊙ ) (M ⊙ ) (M ⊙ ) (M ⊙ pc − ) (M ⊙ pc − ) (yr) (km s − ) (km s − )1 5.6 4.9 4.3 3.9 3.4 5.2 79.4 13.22 5.4 4.7 4.3 4.6 4.4 4.7 64.7 18.23 6.2 5.1 4.1 4.1 3.7 5.1 443.8 16.74 5.1 5.1 5.0 4.4 3.9 4.9 18.6 22.05 5.7 5.3 5.4 4.8 4.5 4.7 20.1 33.26 5.9 3.6 5.3 4.5 4.1 4.9 0.8 21.87 5.5 4.5 4.5 3.7 3.2 5.3 17.9 10.78 5.5 5.2 4.8 4.6 4.2 4.8 51.3 23.39 5.7 4.7 5.5 4.5 4.1 4.9 3.5 26.510 5.7 5.2 5.3 4.3 3.7 5.1 17.1 22.411 6.2 4.5 5.6 4.5 4.0 4.9 3.5 26.812 6.4 4.1 6.0 4.6 3.9 5.0 0.5 35.313 5.7 5.2 4.8 4.8 4.5 4.6 89.4 27.414 5.7 5.7 5.5 5.1 4.8 4.5 53.3 45.5 a Dynamical mass from the virial theorem. Dominant uncertainty is sub-beam structure, including whether the sourcebreaks up into multiple smaller sources, ∼ . ≈ . b Gas-mass based quantities. Uncertainties from the assumed dust temperature, dust-to-gas ratio, and opacity. Likelymagnitude is ∼ . − . c Zero age main sequence stellar mass needed to produce the observed 36 GHz emission as free-free following Equations2 and 3. Mild uncertainty due to the assumed temperature, Gaunt factor, and possible contamination by synchrotron.Larger uncertainties regarding the amount of ionizing photons absorbed by dust and the possibility of pre-main sequencestars in the source.
Note —Masses estimated following Section 4.3. The uncertainty in all of the mass estimates are dominated by systematicuncertainties. We note the dominant uncertainty for each quantity in the associated footnote. M VT refers to thedynamical mass estimated from the virial theorem. M gas refers to gas mass estimated from dust emission at 350 GHz. M ⋆ refers to the zero age main sequence stellar mass needed to match the ionizing photon production rate estimatedfrom the 36 GHz emission. Σ gas+ ⋆ is the estimated gas plus stellar surface density within the 2-d FWHM given our sizeand mass estimates. ρ gas+ ⋆ is mass volume density within the 3-d FWHM given our given our size and mass estimates. τ ff is the gravitational free-fall time implied by that density. p r /M ⋆ is the equivalent radial momentum per unit stellarmass calculated from the gas velocity dispersion, gas mass, and stellar mass (Equation 11). temperature of 7000 K (Draine 2011). A higher ionizingphoton flux requires more ionized gas to be present, anda larger Hii region implies more ionized gas mass. Theorder of magnitude for our M ion agrees with the cal-culations by Ulvestad & Antonucci (1997), though theadopted distances and other details do vary.Measuring the sizes of the Hii regions will help con-strain this measurement, and is a natural next direc-tion. For r S ∼ M ion /M gas ∼
50% for our fiducialassumptions), source ∼ ∼ <
10% and usually . n ion ≈ − (cid:18) Q s − (cid:19) . (cid:18) r S (cid:19) − . (5) yields densities mostly in the range n ion ∼ − cm − for r S = 1 pc. But this depends stronglyon r S . A measured value of r S will allow us to determinif the Hii regions are overpressured and evolving (and,e.g., to compare to Krumholz & Matzner 2009) or inapproximate pressure equilibrium with the surroundinggas (e.g., as in the center of IC 342 Meier et al. 2011).4.3.3.
Gas Mass From Dust
We estimate the mass of gas associated with each pro-tocluster candidate from the 350 GHz dust emission.To do this, we estimate the optical depth at the peak bycontrasting the measured brightness with an estimate ofthe true dust temperature. Then we convert the opticaldepth to a dust column using an assumed mass absorp-tion coefficient. We convert the dust to a gas columnvia an adopted dust-to-gas ratio. Finally, we scale thegas column at the center of the proto-SSC by the areaof the peak to calculate a total gas mass. In the future,we hope to constrain the dust-to-gas ratio by comparingour dust mass estimates to gas mass estimates based onthe gas emission lines. At the moment, we consider the2
Leroy, Bolatto et al.
Figure 5.
Estimated mass budget in our candidate forming clusters. ( left ) Stellar mass ( y -axis), estimated from the 36 GHzcontinuum emission assuming that it all arises from free-free emission and is produced by a population of stars on the zero agemain sequence, as a function of gas mass ( x -axis), estimated from the ALMA-observed dust continuum. ( right) Combined gasplus stellar mass ( y -axis) as a function of dynamical mass ( x -axis) estimated from the measured size and line width of eachsource. In both panels the bold line shows equality with dotted lines offset by successive factors of two. The sources show arange of apparent gas richness, but most show at least half of their mass in gas. Given the signatures of dense, compact gas,these structures seem likely to still be forming. On average, we estimate a dynamical mass a factor of ∼ dust a more reliable estimate of the gas content than themolecular lines that we observe.We assume a fiducial temperature of T dust = 130 K,assuming that the clusters coincide with the warmcomponent seen in ammonia spectroscopy (Gorski et al.2017; Mangum et al. 2013) and that the gas and dust arecollisionally coupled. We convert our measured 350 GHzintensity at the peak, I , to a dust optical depth via I = [1 − exp ( − τ )] B ν ( T dust ) . (6)Here I is our measured 350 GHz intensity, correctedfor free-free contamination using the value in Table 1,and expressed in cgs units. B ν is the Planck functionevaluated at 350 GHz for our adopted dust temperature.This formulation deals better with mild optical deptheffects than assuming the emission to be optically thin.However, if the emission is strongly clumped within ourbeam then these optical depth effects will be underesti-mated.Equation 6 yields optical depths at 350 GHz mostlyin the range 0 .
035 to 0 .
35, with τ ∼ .
09 on aver-age. The dust appears to be moderately optically thinat 350 GHz. Because of the frequency-dependent opac-ity of dust, τ ν ∝ ν β with β ∼ . −
2, these values implythat these sources will be quite optically thick at higherfrequencies (shorter wavelengths) where most of the en-ergy is emitted. After calculating τ , we convert to a dust mass sur-face density using an assumed mass absorption coeffi-cient, κ . We adopt κ = 1 . g − . According toOssenkopf & Henning (1994), this should be appropri-ate for ν ∼
350 GHz and dust mixed with gas at density ∼ − cm − , but this value is uncertain by a factorof ∼ gas = 1DGR κ τ (7)We then scale this Σ gas by the physical area of thepeak, assuming each source to be a two dimensionalGaussian with the size quoted in Table 1. Thus, M gas = A Σ gas .We report the the results in Table 2. We find me-dian log M gas [M ⊙ ] ∼ . M gas [M ⊙ ] ∼ . − . M ⋆ , M gas represents a nearly linear transfor-mation of the measured source flux at 350 GHz. Becausewe detect the sources at high S/N, the statistical uncer-tainties are quite low. Systematic uncertainties in the orming Super Star Clusters in NGC 253 M gas .Based on the excitation requirements of the lines thatwe see, and on the global spectral energy distribution(SED), T dust seems unlikely to be lower than ∼ −
60 Kin these dense, heated regions. Because the clusters arelikely to be optically thick near the peak of the IR SEDswe can ask what temperature, along with our measuredsizes, would place all of the luminosity of the burst inour targets. Assuming L = 4 πr σ SB T , and half of thebolometric IR luminosity from (Sanders et al. 2003) tobe in the burst, we find that T dust must be <
160 K. Weconsider a reasonable plausible range T dust ∼ −
160 K;given that the ammonia temperatures for the “hot”components lie in the intermediate part of this range, T dust ∼
130 K with 50% uncertainty seems like a rea-sonable assessment.As noted, κ appears uncertain by a factor of ∼ ∼
30% uncertainty in the dust-to-gas ra-tio, the overall uncertainty on the gas mass is likely ∼ . − . ≈
30 timeslower gas mass than what one would calculate from the350 GHz light-to-gas-mass conversion of Scoville et al.(2016). That is, we take the dust in these proto-SSCsto be more emissive and much hotter than typical dustin galaxies. Bearing this in mind, we consider that ourgas masses are most likely to be underestimates.4.3.4.
Dynamical Masses
We estimate dynamical mass of each source via M VT = 892 ℓ FWHM σ v . (8)Here σ v the measured velocity dispersion (in km / s), ℓ FWHM is the full width half max deconvolved size ofthe source (both from Table 1), and M VT is the virialtheorem-based dynamical mass in units of solar masses.The prefactor here assumes a density profile ρ ∝ r − .Based on this calculation, we find median dynami-cal mass log M VT [M ⊙ ] ∼ . M VT [M ⊙ ] ∼ . − . M VT to represent an overestimate. We esti-mate the systematic uncertainty due to unresolved sub-structure to be . . Comparison of Mass Estimates
Figure 5 compares our mass estimates. Our sourcesshow gas and stellar masses ∼ up to ∼ M ⊙ . Just as the sizes that we measure are typical of youngcluster sizes (e.g., Ryon et al. 2017), these masses re-semble those seen for massive clusters in nearby star-bursts (e.g., Whitmore 2003; McCrady et al. 2005). Oursources thus do meet the definition of young massiveclusters suggested by, e.g., Portegies Zwart et al. (2010)and Longmore et al. (2014).Our observations suggest a range of gas richness forthe targets, but the left panel of Figure 5 shows thatgas often contributes a large fraction of the mass. Inall but four sources gas contributes &
50% of the mass(and bear in mind that we are suspicious of the M ⋆ forsources M gas / ( M gas + M ⋆ ) across the sample is ∼ ∼ .
5. Given the uncertaintiesin the gas and stellar mass estimation, this still repre-sents reasonable agreement. The dynamical mass es-timate reinforces that these structures mostly contain ∼ − M ⊙ in a region a few pc across, with largecontributions from both stars and gas. The discrepancybetween the two total mass estimates could, in prin-ciple, reflect out-of-equilibrium motions (outflows or in-flows). Or it might indicate that the sources are in a non-virialized dynamical state, for example, the line widthsmight reflect blended motion of unassociated sources.Just as likely, the discrepancy between the virial andgas plus stellar masses between the large uncertaintiesin our mass estimates, especially ionizing photons ab-sorbed by dust and our uncertainties in κ and T dust .4.3.6. Density and Free Fall Time
Table 2 also reports the total (gas plus stellar) surfacedensity, volume density, and implied gravitational freefall time. These are calculated within the FWHM, sothat in two dimensions we divide the mass by 2 anddivide by the area at FWHM. In three dimensions, wedivide the mass by 3 . Σ gas+ ⋆ [M ⊙ pc − ] ∼ . . − .
1. Recasting these val-ues in terms of mass surface density from the edgeto the center of the structure (i.e., converting unitsand dividing by 2), our calculations imply a median ∼ . − from the center to the cluster edge anda range 0 . − . − . These values resemble thosefound in the Milky Way for other regions of high massstar formation (e.g., see McKee & Ostriker 2007). Thehigh end of our range of measured surface densities ap-proaches the ∼
20 g cm − maximum surface density(now through the whole object, not center to edge) forstellar systems found by Hopkins et al. (2010). On av-erage, these proto-clusters are a factor of ∼ Leroy, Bolatto et al.
The median gas plus stellar volume density in our tar-gets is log ρ gas+ ⋆ [M ⊙ pc − ] ∼ . . − .
5) inunits of M ⊙ pc − . This would correspond to a median n H ∼ cm − if all of the material were in moleculargas. The gravitational free fall times implied by thesedensities will be log τ ff ∼ . . − .
3) years.Considering only the gas mass, the implied surfacedensities for our sources would correspond to a median ∼
500 mag (range 20-4200 mag) of V -band extinctionfor a Milky Way dust-to-gas ratio (Bohlin et al. 1978).This helps explain why most of our targets do not appearas distinct sources in the HST imaging.The final column of Table 2 quotes the escape veloc-ity p GM/r calculated within the 3-d FWHM of thesource using the gas plus stellar masses. Using the datareported in the tables: v esc = vuut G (cid:16) M ⋆ + M gas . (cid:17)(cid:0) ℓ FWHM (cid:1) (9)where factors of 3 . ℓ FWHM refersto the FWHM, deconvolved size of the source from Table1. These v esc for all of our sources exceeds the ∼
10 km s − sound speed expected for photoionized gas.As a result, the clusters should match the definition foryoung massive protoclusters from Bressert et al. (2012).4.4. Notes on Individual Sources
As already mentioned, sources f ff is quite high (see Table 1). More, thecontrast with local background or extended features as-sociated with the bright, nearby source < DISCUSSIONWe identify 14 candidate young super star clusters inthe inner region of NGC 253. How much of the starformation in the burst can these sources account for?What can we say about the efficiencies and timescalesassociated with these sources? And given such intenseconcentrations of gas and young stars, what can we inferabout feedback in these sources?5.1.
Timescales and Relation to the Starburst
The central burst in NGC 253 forms ∼ ⊙ yr − (Leroy et al. 2015; Bendo et al. 2015). How much ofthat can be attributed to these sources? Fraction of Ionizing Photons and IR Luminos-ity From These Sources:
We estimate a total ionizingphoton production of Q ∼ . × s − from our tar-gets, with half of this coming from sources Q = 3 . ± . × s − for the whole burst. Thus, our sources may produce be-tween 20 and 40% of the total ionizing photons in theburst. Accounting for absorption of ionizing photonsby dust, which must be more significant in our targetsthan any less embedded population, would increase thisfraction.An analogous case holds for the bolometric luminos-ity. If 50% of the total infrared luminosity measuredby Sanders et al. (2003) arises from the nuclear region,then L TIR ∼ . × L ⊙ for this region. We esti-mate the contribution of our sources to this value bytaking the light-to-mass ratio of a ZAMS population tobe Ψ ≡ L ⋆ /M ⋆ ∼ ⊙ M − ⊙ and adopting the M ⋆ calculated from the 36 GHz emission. In this case ourclusters together contribute ∼
17% of the bolometric lu-minosity of the burst; neglecting sources ∼ Relevant Timescales:
Several distinct timescalesshould combine to produce our observations. First, thetimescale for cluster formation is the time spent in thiscompact, gas-rich phase. Skinner & Ostriker (2015) finda typical timescale of ∼ τ ff for cluster formation. Theforming clusters simulated by Skinner & Ostriker (2015)do evolve over this time. The phase in which gas is ac-tively collapsing to make stars is even shorter, 1 − τ ff ,and the surrounding gas dispersed by ∼ τ ff . Theymeasure their τ ff averaged over their r = 10 pc cloud,with τ ff ∼ . ∼ yr (Table 2). Assuming theSkinner & Ostriker (2015) results to scale with the free orming Super Star Clusters in NGC 253 ∼ × yr.Second, ionizing photon production declines rapidlyafter ∼ − × yr. This should be the timescaleto produce the overall Q in the burst. Third, in-frared or bolometric luminosity is produced over a longertimescale than ionizing photons, with a single stellarpopulation still producing significant light for many tensof Myr.This short cluster formation timescale, ∼ yr, im-plies that a large amount fraction of stars in the burstare born in clusters (see next section). Our observationsdo support the idea that the clusters are young. Below,we show that the total radial momentum in these clus-ters appears low relative to their stellar mass. This im-plies that feedback has not yet unbound the protoclus-ters, consistent with an age young enough that a largeamount of supernovae have not yet gone off. The clus-ters also show signs of ionizing photon production fromcompact regions. Still, these radial momentum limitsand the presence of ionizing photons only place hard lim-its of . −
10 Myr on the age of the clusters. The valueof ∼ τ ff ∼ Are Most Stars in the Burst Born in Clusters?
The estimated timescale for cluster formation, ∼ yr,is ∼ −
30% of the timescale for ionizing photon pro-duction. If all stars are born in clusters, we expect20 −
30% of the ionizing photons coming from the burstat any given time to arise from still-forming clusters. Inthis case, our observations agree with a large fraction( ∼ M gas+ ⋆ ≈ × M ⊙ , split approximately equallybetween gas and stars. Combining this total masswith the 8 τ ff ∼ Myr cluster formation timescalebased on Skinner & Ostriker (2015), then our observedsources can already account for more than the total ∼ ⊙ yr − SFR in the burst. This assumes contin-uous star formation at the time-average rate, leveragesour assumed cluster formation timescale, adopts a 100%gas to star conversion efficiency, and relies on our some-what uncertain mass estimates. All of these assumptionslikely break down in detail. But the calculation showsthat the clusters represent a large fraction of the massthat the burst has likely formed over the last ∼ yr.Observations and theory both predict a larger frac-tion of star forms born in clusters in regions of highstar formation surface density (e.g., Kruijssen 2012;Johnson et al. 2016; Ginsburg & Kruijssen 2018). Thenucleus of NGC 253 has one of highest star formationsurface densities in the local universe. Finding ∼ A Plausible Scenario:
Our observations appearconsistent with a scenario in which most of the starsin the burst form in massive young clusters. The for-mation process last for ∼ Q in the burst inferred from free-free and radiorecombination line emission (Bendo et al. 2015). Mean-while their corresponding Hii regions would grow in size,fading in surface brightness and becoming much moredifficult to pick up in our interferometric radio contin-uum maps. Eventually, they might evolve in analogsto the older, visible clusters seen at larger radii byFern´andez-Ontiveros et al. (2009).We expect the strong feedback that drives the X-ray and molecular gas winds (Strickland et al. 2002;Bolatto et al. 2013; Walter et al. 2017) to occur afterthe embedded young cluster phase that we observe. Af-ter ∼
10 Myr, many massive stars will explode as su-pernovae. These explosions may trigger both the hotand cold outflows. A scenario in which strong feedbackoccurs well after the embedded phase agrees with ourobservations, which show that the protocluster candi-dates are approximately gravitationally bound at 2 pcscales, with no evidence for high velocity line wings intheir spectra. In this case the clusters that we observenow are not the immediate drivers of the outflows ob-served in Bolatto et al. (2013) Walter et al. (2017), andZschaechner et al. (2018). They may, however, drivesimilar outflows in the future.
Lower Mass Clusters:
We observe candidateprotoclusters down to a combined gas plus stel-lar mass of log M gas+ ⋆ ∼ .
8. The cluster massfunction is often taken to have equal power perdecade (e.g., Portegies Zwart et al. 2010). Takingthe cutoff for young massive clusters as ∼ M ⊙ (Portegies Zwart et al. 2010; Longmore et al. 2014;Bressert et al. 2012), there may be as much mass inlow mass, unidentified clusters as in the sources that westudy. Likely many of these lower mass sources will besubstructure still unresolved by our 1 . Leroy, Bolatto et al.
Of course, if twice as much mass — and ionizing pho-tons and bolometric luminosity — are present in clustersoutside our sources then our bookkeeping above breaks.This could indicate a longer timescale for cluster for-mation or it might reflect that we have systematicallyoverestimated the mass of the clusters.
Lower Limit:
As emphasized, the cluster formationtimescale remains uncertain. Our mass estimates arealso uncertain at the factor of two level. Given that allof the sources show 36 GHz flux and that at least 20%of the ionizing photon production occurs in our sources,a reasonable limiting case is that the visibility lifetimefor the clusters equals the ∼ − ∼
20% of the stars in the burstform in these structures. Even in this limiting case, theburst represents a prodigious cluster production factory,far more extreme than what we see around us in theMilky Way.5.2.
Likely High Efficiency Per Free Fall Time
We find M ⋆ ∼ M gas and free fall times t ff ∼ yrbased on the combined gas plus stellar mass in the clus-ters. Assuming that no mass has escaped from the clus-ter since its initial formation, then M ⋆ ∼ M gas impliesan overall efficiency of ∼ ∼ t ff , thisimplies an efficiency per free fall time of > ∼ ∼
5% of the gasmass converted to stars per free fall time.Note that some mass-loss may already have occurred(e.g., as might be expected following Thompson & Krumholz2016), in which case the efficiency per free-fall timewould be lower than we calculate. In order to have a ∼
1% efficiency per free fall time as is observed at largerscales (e.g., Krumholz & Tan 2007; Krumholz et al.2012; Utomo et al. 2018), several times the currentlyobserved gas must have already been expelled; however,this seems unlikely based on observed linewidths.5.3.
High Infrared Opacity
The deconvolved 350 GHz peak brightness tempera-tures associated with our sources are high, often ∼
10 Kand in a few cases 20 −
40 K or more. We do not knowthe true dust temperature, but our arguments abovesuggest that it cannot be much more than ∼
130 K onaverage. In this case, the dust optical depth at 350 GHzis already τ ∼ . Opaque at IR wavelengths:
For dust, τ ∝ ν β with β ∼ . − τ at 350 GHz, this im-plies that dust continuum emission from our compactsources will be optically thick for wavelengths shorterthan λ ∼ − µ m. They will have a factor of ∼ −
100 higher optical depth at λ ∼ µ m com- pared to 350 GHz (855 µ m). This yields optical depths τ ∼ −
10 at 100 µ m, and much larger near the impliedpeak of the dust SED at ∼ − µ m.With such high optical depths, these cluster-formingstructures might be expected to have large IR photo-spheres. The regions could appear much larger at IRwavelengths than at sub-mm wavelengths, so that re-solving them is only possible with ALMA. Clumpy sub-structure might render this a more local effect, so thatthe gross morphology of the sources does not change, butthe opacity effects must be present at some scale. Thisclumpy substructure might be expected from compar-ing the mean particle densities, n ∼ cm − , with thetypical densities needed to excite the bright H CN (4-3) and CS (7-6) emission we measure, which requiresdensities n ∼ cm − . Significant Infrared Radiation Pressure Force:
This high opacity at IR wavelengths also implies a strongradiation pressure force exerted by the cluster stars onthe surrounding gas (Murray et al. 2010). For spheri-cal systems optically thick to stellar radiation, the stel-lar radiation creates an outward force L ⋆ /c , with L ⋆ the bolometric luminosity. For systems that are opti-cally thick in the infrared, the reprocessed infrared lightalso contributes to this outward force. This force dueto infrared radiation force exceeds the force associatedwith the primary stellar radiation by a factor equal tothe Rosseland mean optical depth ∼ τ IR . The high τ IR in our clusters thus implies a strong radiation pressureforce on the surrounding gas. Effect of Radiation Pressure on Cluster For-mation:
How does this high radiation pressure affectcluster formation? Following Murray et al. (2010) andSkinner & Ostriker (2015), for a spherical system cen-tered on a star cluster, the ratio of the IR radiationforce to the gravitational force from the stars is f Edd , IR = κ IR , gas F IR /cGM ∗ /r = κ IR , gas Ψ4 πcG , (10)where F IR = L ⋆ / (4 πr ) is the IR flux, assuming allstarlight to be reprocessed into the IR. Here κ IR , gas isthe mass absorption coefficient per unit gas mass ; notethe difference from above where we discuss the mass ab-sorption coefficient of dust alone. Thus κ IR , gas includesboth the dust properties and the dust-to-gas ratio. HereΨ ≡ L ⋆ /M ⋆ refers to the light-to-mass ratio of the cen-tral stellar population.If f Edd , IR > f Edd , IR > ∼ f Edd , IR < orming Super Star Clusters in NGC 253 κ IR , gas is set bydust properties and the dust abundance relative to gas.Adding gas to the system increases the total dust opac-ity, leading to a higher τ IR and more support from ra-diation pressure. But at the same time, adding gas tothe system increases the weight of gas. Because thesetwo effects balance, a high τ IR does not necessarily im-ply anything about force balance in the cluster (thoughthere can be an indirect dependence of κ IR , gas on τ IR through the dust temperature Krumholz & Thompson2012).For a zero age main sequence with a Kroupa initialmass function, Ψ ∼ − g − ∼ L ⊙ M − ⊙ .For temperature range and dust abundance relevantto our clusters, Semenov et al. (2003) find a Rosse-land mean opacity κ IR , gas . g − . In this case f Edd . . &
2. In this case, ra-diation pressure would help support the cloud againstcollapse, but not supply all of the support nor tearthe cloud apart. Including gas self-gravity would onlystrengthen the effects of gravity relative to radiationpressure.This situation could change if Ψ > L ⊙ M − ⊙ ,e.g., due to a top-heavy IMF. Top-heavy IMFs havebeen claimed in 30 Doradus (Schneider et al. 2018) andthe proto-SSC in NGC 5253 (Turner et al. 2017). Al-ternatively, if the gas associated with the clusters hasa higher than Galactic dust-to-gas ratio, or unusuallyopaque grains, κ IR would be higher than assumed above.We find virial masses within a factor of ∼ M ⋆ + M gas . This supports a scenario in which the clus-ters are gravitationally bound in approximate equiib-rium. It appears that radiation forces, though certainlyenhanced by a high τ IR , at most balance gravity at thepresent time, consistent with the expectations above. Atpresent, we lack independent constraints on the dust togas ratio, nor do we independently measure L ⋆ and M ⋆ .5.4. Limits on Feedback From Observed RadialMomentum
The correspondence between virial masses and M ⋆ + M gas implies that gravity approximately balances theoutward force in our clouds. Given enough time, bothsupernovae (SNe) and stellar winds can inject enoughmomentum to unbind the gas and drive a radial expan-sion. The contrast between the observed radial motionsin our sources and the expected momentum injectionfrom SNe and stellar winds provides additional indirectsupport for the idea that our sources are young. Observed Limits on Radial Momentum:
Takingall motions to be radial and outward, the momentumper unit stellar mass for an expanding spherical systemis p r M ⋆ ≡ √ σ v M gas M ⋆ . (11) This p r will be the maximum radial momentum com-patible with an observed velocity dispersion σ v and gasmass M gas . Normalizing by M ⋆ allows a direct compar-ison with input from SNe and stellar winds, which bothscale with stellar mass.We report p r /M ⋆ limits for our sources in Table 2.We find mostly p r /M ⋆ <
100 km s − , with the largestvalue ∼
400 km s − for source Momentum From Supernova Feedback:
Nu-merical simulations considering clustered SNe explod-ing in an inhomogenous medium find a momentuminjection per SN (after cooling and shell formation)of ∼ M ⊙ km s − (e.g. Kim et al. 2017, and refer-ences therein). For a Kroupa IMF, with roughly onesupernova per 100 M ⊙ formed, we expect p ⋆ /M ⋆ ∼ km s − at late times ( ∼
10 Myr). This is an or-der of magnitude higher than what we observe for mostsources, suggesting that supernovae have not yet had asignificant effect on internal motions.The large values of p ⋆ /M ⋆ are associated with along timescale, t SN &
10 Myr. Spreading p ⋆ /M ⋆ ∼ km s − across this t SN , the mean momentum in-jection rate, ( p ⋆ /M ⋆ ) /t SN , may not exceed the gravita-tional force, GM gas+ ⋆ M gas / ( M ⋆ r ) ∼ ( p r /M ⋆ ) /t ff , espe-cially at early times.This argument does not preclude any supernovae hav-ing gone off. Kornei & McCrady (2009) note the pres-ence of iron lines in their SSC (our source .
10 Myr.
Momentum From Stellar Winds:
For stellarwind feedback, the pressure-driven bubble solution ofWeaver et al. (1977) yields a ratio of shell momentumto central cluster mass of p ⋆ M ⋆ = 65 km s − ˙ E / w, n / M − / ⋆, t / . (12)Here ˙ E w, is the average wind luminosity injected per M ⊙ of stars in units 10 erg s − ; n is the mean hydro-gen density in units 10 cm − ; M ⋆, is the cluster massin units 10 M ⊙ ; and t is the cluster age in units 10 yr.From Starburst99, ˙ E w, = 1 (Leitherer et al. 1999).This calculation also predicts the radius of the wind-driven bubble. If there are no energy losses and gravityis negligible, r b = 3 pc ˙ E / w, n − / M / ⋆, t / . (13)Though we are not yet in a position to measure therelative structure of the ionized gas, molecular gas, anddust, we do not expect r b to exceed our observed source8 Leroy, Bolatto et al. size. In that case, the stellar winds would have clearedout the cold gas.If the sources in NGC 253 are young ( t < ∼ t ff , then the pre-dicted bubble radius and wind momentum input maysignificantly exceed our observed limits.Again, it seems likely that the momentum injectionfrom the wind has been balanced by gravity. Comparingthe predicted momentum input rate from a wind-blownbubble, ( p ⋆ /M ⋆ ) /t ∝ t / to the gravitational force,( p r /M ⋆ ) /t ff , this calculation suggests that the force fromwinds should exceed that of gravity for t ∼ yr.Again, the lack of strong signatures of gas expulsionargue that our sources are young, . yr. In thiscase, the effective wind luminosity must be reduced be-low the expected input value, either by mixing and cool-ing or by other processes. This is reasonable basedon the low X-ray emission observed in somewhat moreevolved systems, where the energy in hot gas seems tobe far below the value nominally expected from winds,consistent with a reduction in ˙ E w, well below unity(e.g. Harper-Clark & Murray 2009; Lopez et al. 2011;Rosen et al. 2014). SUMMARYWe present new, ∼ Massive, Compact Sources:
We measuresource sizes from the ALMA 350 GHz contin-uum (Table 1). We estimate gas masses basedon dust emission, and calculate stellar masses as-suming that the observed 36 GHz continuum isfree-free emission from a zero-age main sequencepopulation (Table 2).We find sizes of a few pc (FWHM) and estimatetotal masses M gas + M ⋆ & M ⊙ . We also es-timate dynamical masses from the measured sizesand line widths and assuming virialization. Thevirial masses are typically ∼ . M gas + M ⋆ estimates, which represents rea- sonable within the substantial uncertainties on themass estimation.2. Likely Young Super Star Clusters:
Thesemasses and sizes resemble those of young massiveclusters seen in the Milky Way and other galax-ies (Portegies Zwart et al. 2010; Longmore et al.2014). More, these masses and radii imply es-cape speeds >
16 km s − . This is larger than thesound speed of photo-ionized gas, ∼
10 km s − , sothat the sources also match the criteria for youngmassive protoclusters laid out by Bressert et al.(2012). Clusters in this mass range are often re-ferred to as super star clusters.3. Still in the Process of Formation:
Our esti-mates of the gas and stellar mass, while uncertain,suggest that gas still contributes a large fraction ofthe total mass in these objects (Table 1). We ob-serve that the dust emission coincides with H CN(4-3) and CS (7-6) emission, both tracers of dense,excited gas. Thus, many of these objects seemlikely to still be in the process of formation.4.
Short Free-Fall Times and High Efficiency:
The free fall times implied by the gas plus stel-lar masses of our sources is short, τ ff ∼ yr.Given theoretical expectations of ∼ τ ff for thecluster formation timescale (Skinner & Ostriker2015), this implies that the sources are young.Stars typically represent half of the mass in oursources, implying both a high net efficiency anda high efficiency per free fall time. Note that thisstatement does not take into account possible massloss. Any mass lost from the system would de-crease the both the net star formation efficiencyand the efficiency per free fall time.5. A Large Fraction of Stars Form in theseSources:
At least 20% of the ionizing photonproduction in the burst appears associated withthese sources. This represents a firm lower limiton the fraction of stars that form in such sources.If the cluster formation timescale is short com-pared to the time for stars to produce ionizing pho-tons, then an even larger fraction of star formationmay proceed through this phase. Accounting for ashort cluster formation timescale and the possibil-ity of lower mass, still-unidentified clusters, order ∼ Opaque in the Infrared:
These sources havehigh brightness temperature. Given plausible dusttemperatures, they also have moderate ( τ ∼ . orming Super Star Clusters in NGC 253 M gas + M ⋆ ,this force may help support the clouds but is notunbinding them. This agrees with theoretical ex-pectations.7. Young Based on Large Gas Fraction and Be-ing Approximately Bound:
Our sources retaina large fraction of their mass in gas (as evidencedby the dust continuum). They also appear tobe approximately gravitationally bound. We alsocalculate limits on the radial momentum in oursources and compare them to expectations fromsupernova and stellar wind feedback. Our sourceshave lower radial momentum and smaller sizesthan expected from either clustered supernova orstellar winds acting over many Myr, though lossesin wind energy may be important. All of thesepieces of evidence suggest that the sources areyoung enough that feedback has not managed tounbind the gas from the cluster.Given the brightness of these sources, ALMA and theJansky VLA both offer the prospect for even more de-tailed detailed follow up. Higher resolution dust ob-servations are already underway, as is the construc-tion of full radio-to-mm SEDs for each source (buildingon Ulvestad & Antonucci 1997; Mohan et al. 2005, andleveraging new ALMA and VLA work).It will also be important to link these structures largercontext of the burst. In the Milky Way’s Central Molec-ular Zone, star formation has been linked to the orbitalpaths of individual clouds (e.g., Kruijssen et al. 2014).The linear distribution of the sources we see suggestsan underlying bar-like structure (see Leroy et al. 2015;Paglione et al. 2004) or loosely wound arms. It may bepossible to link this structure to the triggering of starformation. More generally, we do not see clear analogsfor these structures in the Milky Way. This might be because NGC 253 sits at a different part of some longterm nuclear fueling cycle (e.g., Krumholz et al. 2017).A more detailed comparison of the two systems (buildingon Sakamoto et al. 2011) could help reveal the overalltriggers and likely duty cycle of the burst. This mightalso help reveal the fate of the proto-clusters after theydisappear from our ALMA and VLA imaging, and per-haps link them to clusters seen on larger scales outsidethe area we study (Fern´andez-Ontiveros et al. 2009).We thank the anonymous referee for a constructive re-port that improved the paper. We also thank GerhardtMeurer and Mark Krumholz for useful feedback duringrevision.This paper makes use of the following ALMA data:ADS/JAO.ALMA (Tully et al. 2009), the HyperLedadatabase (Makarov et al. 2014), the NASA/IPAC Ex-tragalactic Database , and the SAO/NASA Astro-physics Data System . Facilities:
ALMA, VLA, HST
Software:
CASA (McMullin et al.2007),
IDL ,CPROPS (Rosolowsky & Leroy 2006; Leroy et al. 2015)REFERENCES Ando, R., Nakanishi, K., Kohno, K., et al. 2017, ArXive-prints, arXiv:1710.01432Bendo, G. J., Beswick, R. J., D’Cruze, M. J., et al. 2015,MNRAS, 450, L80 http://edd.ifa.hawaii.edu/index.html http://leda.univ-lyon1.fr http://ned.ipac.caltech.edu https://github.com/akleroy/cpropstoo Binder, B. A., & Povich, M. S. 2018, ArXiv e-prints,arXiv:1808.00454Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ,224, 132Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013,Nature, 499, 450Bressert, E., Ginsburg, A., Bally, J., et al. 2012, ApJL, 758,L28Caplan, J., & Deharveng, L. 1986, A&A, 155, 297Dale, D. A., Cohen, S. A., Johnson, L. C., et al. 2009, ApJ,703, 517 Leroy, Bolatto et al.
Draine, B. T. 2011, Physics of the Interstellar andIntergalactic MediumEldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017,Publications of the Astronomical Society of Australia, 34,e058Fern´andez-Ontiveros, J. A., Prieto, M. A., &Acosta-Pulido, J. A. 2009, MNRAS, 392, L16Fukui, Y., Torii, K., Ohama, A., et al. 2016, ApJ, 820, 26Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271Ginsburg, A., Bressert, E., Bally, J., & Battersby, C. 2012,ApJL, 758, L29Ginsburg, A., & Kruijssen, J. M. D. 2018, ApJL, 864, L17Ginsburg, A., Bally, J., Barnes, A., et al. 2018, ApJ, 853,171Gorski, M., Ott, J., Rand, R., et al. 2017, ApJ, 842, 124Harper-Clark, E., & Murray, N. 2009, ApJ, 693, 1696Herrera, C. N., Boulanger, F., Nesvadba, N. P. H., &Falgarone, E. 2012, A&A, 538, L9Holtzman, J. A., Faber, S. M., Shaya, E. J., et al. 1992, AJ,103, 691Hopkins, P. F., Murray, N., Quataert, E., & Thompson,T. A. 2010, MNRAS, 401, L19Johnson, K. E., Leroy, A. K., Indebetouw, R., et al. 2015,ApJ, 806, 35Johnson, L. C., Seth, A. C., Dalcanton, J. J., et al. 2016,ApJ, 827, 33Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834,25Knudsen, K. K., Walter, F., Weiss, A., et al. 2007, ApJ,666, 156Kormendy, J., & Kennicutt, Jr., R. C. 2004, ARA&A, 42,603Kornei, K. A., & McCrady, N. 2009, ApJ, 697, 1180Kroupa, P. 2001, MNRAS, 322, 231Kruijssen, J. M. D. 2012, MNRAS, 426, 3008Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G.,et al. 2014, MNRAS, 440, 3370Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ,745, 69Krumholz, M. R., Kruijssen, J. M. D., & Crocker, R. M.2017, MNRAS, 466, 1213Krumholz, M. R., & Matzner, C. D. 2009, ApJ, 703, 1352Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304Krumholz, M. R., & Thompson, T. A. 2012, ApJ, 760, 155Lee, J. C., Gil de Paz, A., Tremonti, C., et al. 2009, ApJ,706, 599Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999,ApJS, 123, 3Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015,ApJ, 801, 25 Longmore, S. N., Kruijssen, J. M. D., Bastian, N., et al.2014, Protostars and Planets VI, 291Longmore, S. N., Walsh, A. J., Purcell, C. R., et al. 2017,MNRAS, 470, 1462Lopez, L. A., Krumholz, M. R., Bolatto, A. D., Prochaska,J. X., & Ramirez-Ruiz, E. 2011, ApJ, 731, 91Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., &Vauglin, I. 2014, A&A, 570, A13Mangum, J. G., Darling, J., Henkel, C., et al. 2013, ApJ,779, 33McCrady, N., Graham, J. R., & Vacca, W. D. 2005, ApJ,621, 278McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565McMullin, J. P., Waters, B., Schiebel, D., Young, W., &Golap, K. 2007, in Astronomical Society of the PacificConference Series, Vol. 376, Astronomical Data AnalysisSoftware and Systems XVI, ed. R. A. Shaw, F. Hill, &D. J. Bell, 127Meier, D. S., Turner, J. L., & Schinnerer, E. 2011, AJ, 142,32Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ,801, 63Mohan, N. R., Goss, W. M., & Anantharamaiah, K. R.2005, A&A, 432, 1Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011,ApJ, 737, 67Murray, N., Quataert, E., & Thompson, T. A. 2010, ApJ,709, 191Ochsendorf, B. B., Zinnecker, H., Nayak, O., et al. 2017,Nature Astronomy, 1, 784Oey, M. S., Herrera, C. N., Silich, S., et al. 2017, ApJL,849, L1Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943Paglione, T. A. D., Yam, O., Tosaki, T., & Jackson, J. M.2004, ApJ, 611, 835Paglione, T. A. D., Wall, W. F., Young, J. S., et al. 2001,ApJS, 135, 183Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M.2010, ARA&A, 48, 431Rekola, R., Richer, M. G., McCall, M. L., et al. 2005,MNRAS, 361, 330Rosen, A. L., Lopez, L. A., Krumholz, M. R., &Ramirez-Ruiz, E. 2014, MNRAS, 442, 2701Rosolowsky, E., & Leroy, A. 2006, PASP, 118, 590Rybicki, G. B., & Lightman, A. P. 1986, RadiativeProcesses in Astrophysics, 400Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ,841, 92Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, ApJ,735, 19 orming Super Star Clusters in NGC 25321