The Morpho-Kinematic Architecture of Super Star Clusters in the Center of NGC253

The center of the nearby galaxy NGC\,253 hosts a population of more than a dozen super star clusters (SSCs) which are still in the process of forming. The majority of the star formation of the burst is concentrated in these SSCs, and the starburst is powering a multiphase outflow from the galaxy. In this work, we measure the 350~GHz dust continuum emission towards the center of NGC\,253 at 47~milliarcsecond (0.8 pc) resolution using data from the Atacama Large Millimeter/submillimeter Array (ALMA). We report the detection of 350~GHz (dust) continuum emission in the outflow for the first time, associated with the prominent South-West streamer. In this feature, the dust emission has a width of $\approx$~8~pc, is located at the outer edge of the CO emission, and corresponds to a molecular gas mass of $\sim~(8-17)\times10^6$~M$_\odot$. In the starburst nucleus, we measure the resolved radial profiles, sizes, and molecular gas masses of the SSCs. Compared to previous work at somewhat lower spatial resolution, the SSCs here break apart into smaller substructures with radii $0.4-0.7$~pc. In projection, the SSCs, dust, and dense molecular gas appear to be arranged as a thin, almost linear, structure roughly 155~pc in length. The morphology and kinematics of this structure can be well explained as gas following $x_2$ orbits at the center of a barred potential. We constrain the morpho-kinematic arrangement of the SSCs themselves, finding that an elliptical, angular momentum-conserving ring is a good description of the both morphology and kinematics of the SSCs.


INTRODUCTION
Nuclear starburst regions in galaxies are thought to be fueled by inflows of cold gas to their centers. These gas inflows can be driven by a strong bar, a merger, or tidal interaction. In the case of a barred system, the bar efficiently funnels gas toward to center, where it set-and evolution of massive star clusters. However relatively few studies have dissected the star forming central molecular zones of galaxies at the high (∼ 1 pc scales) resolution needed to distinguish the locations of cluster formation. In this paper we conduct such an analysis targeting NGC 253-the nearest bar-fed nuclear starburst to the Milky Way-by measuring the sizes and masses of the forming massive star clusters and constraining their 3D distribution and kinematics in relation to the bar. We also measure the properties of dust detected in the large scale outflow, a result of a central starburst.
In the case of the nearby galaxy NGC 253, a strong bar fuels the nuclear starburst (e.g., Sorai et al. 2000;Paglione et al. 2004). As a result of the inflowing gas along the bar, the central few hundred parsecs of the galaxy contains ≈ 3.5 × 10 8 M of H 2 and is forming stars at a rate of ∼ 2 M yr −1 , resulting in a molecular gas depletion time of ≈ 300 Myr (Leroy et al. 2015a;Krieger et al. 2019). The nuclear starburst is responsible for launching a massive, multiphase wind which is detected across the electromagnetic spectrum, in Xrays (e.g., Strickland et al. 2000Strickland et al. , 2002, ionized gas (e.g., Heckman et al. 2000;Sharp & Bland-Hawthorn 2010;Westmoquette et al. 2011), polycyclic aromatic hydrocarbons (e.g., McCormick et al. 2013), molecular gas (e.g., Sugai et al. 2003;Sturm et al. 2011;Bolatto et al. 2013;Walter et al. 2017;Zschaechner et al. 2018;Krieger et al. 2019), and radio continuum emission (e.g., Turner & Ho 1985;Heesen et al. 2011). The molecular gas in the outflow is concentrated along the edges of the biconical wind, and the most prominent feature is the so-called South-West (SW) streamer (Bolatto et al. 2013;Walter et al. 2017;Zschaechner et al. 2018;Krieger et al. 2019). The deprojected molecular mass outflow rate of the wind is ∼ 14 − 39 M yr −1 , with large uncertainties due to the geometry (Krieger et al. 2019). With a mass loading factor (defined as the ratio of the mass outflow rate to the SFR) of ∼ 7 − 20, the outflow plays a critical role in regulating the star formation activity in the nucleus.
The nuclear region of NGC 253 is a chemically-rich environment (Krieger et al. 2020a;Martín et al. 2021;Haasler et al. 2022) and hosts a number of massive, dense molecular clouds (e.g., Sakamoto et al. 2011;Leroy et al. 2015a), radio continuum sources (likely HII regions and supernova remnants; Turner & Ho 1985;Watson et al. 1996;Ulvestad & Antonucci 1997;Kornei & Mc-Crady 2009), and masers (e.g., Gorski et al. 2017Gorski et al. , 2019Humire et al. 2022). The overwhelming majority of the star formation in the nuclear starburst is concentrated in massive forming "super" star clusters (SSCs; Ando et al. 2017;Leroy et al. 2018;Rico-Villas et al. 2020;Mills et al. 2021;Levy et al. 2021). These proto-SSCs are responsible for 3% of the total infrared emission of the galaxy (Martín et al. 2021). The gas content and star formation activity in this region resemble a scaledup version of that found in the Central Molecular Zone (CMZ) of the Milky Way (MW; Sakamoto et al. e.g., 2011;Martín et al. e.g., 2021). The SSCs in NGC 253 have stellar masses ≈ 10 4.0−6.0 M and gas masses ≈ 10 3.6−5.7 M (Leroy et al. 2018;Mills et al. 2021). At the resolution of these studies (∼ 2 − 5 pc), however, multiple SSCs are blended together, as revealed by very high (0.5 pc) resolution data of this region .
From high resolution images of the SSCs taken using ALMA, the SSCs are embedded in an extended background of dust and molecular gas, and this structure measures ≈ 155 pc × 15 pc in projected length and width (e.g., Ando et al. 2017;Leroy et al. 2018;Levy et al. 2021;Mills et al. 2021). The observed nearly linear arrangement of the SSCs is almost certainly a projection effect, as NGC 253 has an inclination of ≈ 78 • (e.g., Pence 1980;Westmoquette et al. 2011;Krieger et al. 2019). In 3D, it is possible that the SSCs trace a circumnuclear ring which may be connected to the bar. A promising hint in this direction is that the location of the inner inner Lindblad resonance (IILR) -inside which gas is expected to concentrate -is located at a radius of ∼ 350 pc from the center, though the uncertainty on this radius is likely substantial (Sorai et al. 2000) 1 . Qualitatively, the IILR is on the same scale as the SSC and dense gas structures. While Paglione et al. (2004) find weaker evidence of an inner Lindblad resonance (ILR) than Sorai et al. (2000), they do find that the dense molecular gas in the center is consistent with the locations of x 2 orbits (see e.g., their Figure 11). The x 2 orbits are expected to lie between the outer ILR (OILR) and IILR and are oriented perpendicular to the bar major axis (e.g., Contopoulos & Grosbol 1989;Buta & Combes 1996;Das et al. 2001).
Given the nearly edge-on inclination of NGC 253, inferring a connection with the bar and constraining the geometry of the SSC structure from the 2D morphology alone would be nearly impossible. In this study, we use new images of the dust continuum and dense gas emission in the center of NGC 253, which combine multiple ALMA configurations to cover a wide range of spatial 1 Sorai et al. (2000) assumed a galaxy distance of 2.5 Mpc, whereas a more recent and accurate determination of the distance is 3.5 Mpc (Rekola et al. 2005). We re-calculate all physical sizes from Sorai et al. (2000) using this updated galaxy distance.
scales. These data allow us to simultaneously resolve the compact SSCs and the more extended molecular gas and dust emission. We combine the dust continuum images with the systemic velocities of the clusters measured from very high resolution spectral line data by Levy et al. (2021). This velocity information adds a third dimension to the data, allowing us to better constrain the morpho-kinematic architecture of the SSCs and their connection to the larger scale gas flows in this galaxy.
This article is organized as follows. We describe the observations, data processing steps, and imaging in Section 2. We report the detection of dust emission associated with the SW streamer of the molecular outflow in Section 3. The methods used to measure the cluster sizes, fluxes, and gas masses are described in Section 4. In Section 5, we quantitatively compare the arrangement and kinematic structure of the clusters to an elliptical, angular momentum-conserving ring. We summarize our findings in Section 6.
The precise center of NGC 253 is not known and the location of its central supermassive black hole has not been definitively identified. Throughout this paper, we refer to the galaxy center at (α, δ) J2000 = (0 h 47 m 33.06 s , −25 • 17 18.3 ) measured using ionized gas kinematics traced by H92α at ≈ 1.8 resolution (Anantharamaiah & Goss 1996). The 1-σ uncertainty on this center position is ∼ 0.3 .

OBSERVATIONS AND DATA PROCESSING
Data for this project were taken with ALMA as part of projects 2015.1.00274.S and 2017.1.00433.S (P.I. A. Bolatto). We observed the central 16.64 (280 pc) of NGC 253 at Band 7 (ν ∼ 350 GHz, λ ∼ 0.85 mm) using the main 12-m array in the C43-9, C43-6, and C43-4 configurations and the 7-m (ACA) array. These configurations cover baselines of 113 m − 13.9 km, 15.1 m − 1.8 km, 15.1 m − 783.5 m, and 8.9 m − 49.0 m, respectively (Krieger et al. 2019;Levy et al. 2021). Including the ACA data, the maximum recoverable scale is 24.8 (421 pc); excluding the ACA data, the the maximum recoverable scale is 3.9 (66 pc). The spectral setup spans frequency ranges of 342.08 − 345.78 GHz in the lower sideband and 353.95 − 357.66 GHz in the upper sideband. The visibilities were pipeline calibrated (L. Davis et al. in prep.) using the Common Astronomy Software Application (casa; McMullin et al. 2007). More information on these observations has been published previously (see Leroy et al. 2018;Krieger et al. 2019;Levy et al. 2021).

Imaging the Multi-Configuration Data Sets
In this work, we make two different combinations of the multi-configuration datasets. First, we combine the ACA data and three 12-m configurations together to make what we will refer to as the "12-m+ACA map" or the "2.55 pc resolution map." The objective of this map is to recover the most extended dust continuum emission in the nuclear region. We also make a second multi-configuration data set using only the three 12-m configurations, which we will refer to as the "12-m map" or the "0.81 pc resolution map." The objective of this data set is to recover the dust emission associated with the clusters. The calibrated and line-flagged visibilities were combined for imaging using the concat task in casa. We spectrally averaged the combined measurement set to have 10 channels per sideband, so that each channel covers ∼40 MHz.
Since the 12-m+ACA and 12-m maps have different objectives, we used different deconvolution strategies to produce the final images, which we describe below. All of the visibilities were imaged using the casa version 5.4.1 tclean task.

Imaging the 12-m+ACA Data
We imaged the central 48 ×48 of the line-flagged, channel averaged, combined 12-m+ACA visibilities interactively using tclean. Since we are interested in the more extended dust continuum emission, we choose a coarser cell (pixel) size of 0.04 than Levy et al. (2021) used to image the high resolution continuum data. In all iterations, we used specmode='mfs', deconvolver='multiscale', Briggs weighting with robust=0.5, and apply the primary beam correction. The clean components were modeled using a linear spectral slope (nterms=2) to account for any continuum slope over the bandpass. The "dirty" image (niter=0) is shown in Figure 1 (top left) for the inner 20 ×20 . The dirty map has a FWHM Gaussian beam size of 0.110 ×0.095 . This image is convolved to a circular 0.15 beam, to match the resolution of the cleaned, tapered image described below.
Before cleaning the extended emission, it was necessary to carefully clean the point source-like clusters, otherwise the algorithm had a tendency to over-subtract these regions leaving deep negative bowls. We cleaned the emission from the clusters using scales= The top left image shows the dirty map and the top right shows the cleaned map both convolved to a circular 0.15 beam as described in Section 2.1.1. The blue ellipses highlight the SW streamer seen in the dust continuum emission (see Section 3). Bottom left: The 350 GHz continuum image in the central 10 made by combining the three 12-m configurations, which has a final resolution of 0.047 . The cleaning has been optimized for the cluster-scales, as described in Section 2.1.2. Bottom right: The 350 GHz continuum image in the central 10 of NGC 253 using only the highest resolution (0.028 ) data from Levy et al. (2021). interactively controlled the threshold and number of iterations to avoid over-cleaning. We cleaned the point sources until they were no longer point-like in the residual map and so that the extended residual emission near the point sources was similar to the larger scale emission in the map.
Due to the range of spatial scales covered by these combined data sets, the algorithm tends to favor small scales, making cleaning the extended emission time consuming. Since, for the 12-m+ACA map, we are interested in the larger scale more diffuse emission, we used a uv-taper of 0.2 , scales= [0,8,16], smallscalebias=0 which gives equal weight to all scales to more efficiently clean the map. We used a circular 0.15 restoring beam. To avoid over-cleaning, we reduced the gain of each major cycle to 0.05 and interactively lowered the thresh-old. We interactively cleaned the map until the residuals stopped changing. The final cleaned 12-m+ACA map is shown in Figure 1 (top right), which has an rms of 1.1 mJy beam −1 (0.5 K) in regions away from emission.

Imaging the 12-m Data
We imaged the central 48 ×48 of the line-flagged, channel averaged, combined 12-m visibilities interactively using tclean. Since we are interested in the dust continuum emission associated with the clusters, we use a different imaging strategy from the one described above. We use a cell (pixel) size of 0.0046 , the same as was used for the high resolution continuum map  which is shown for comparison in Figure 1 (bottom right). In all iterations, we used specmode='mfs', deconvolver='multiscale', Briggs weighting with robust=0.5, and apply the primary beam correction. The baseline was fit with a linear function (nterms=2) to account for any continuum slope over the bandpass.
As with the 12-m+ACA map, before cleaning the extended emission, it was necessary to carefully clean the point source-like clusters, otherwise the algorithm had a tendency to over-subtract these regions leaving deep negative bowls. We cleaned the emission from the clusters using scales=[0] and interactively controlled the threshold and number of iterations to avoid overcleaning. We cleaned the point sources until they were no longer point-like and so that the emission in those regions was similar to the larger scale emission in the map.
Once the point sources were "cleaned", we carefully cleaned the more extended emission, starting from large scales and moving to smaller ones. For each iteration, we used smallscalebias=0 and no uv-taper. We first started with scales=[16,32,64] (corresponding to ≈0.07 , 0.15 , and 0.30 ). We interactively cleaned these scales until the maximum residual and cleaned flux no longer changed significantly. We then added scales=[8] (≈0.04 ) to the existing scales and continued to clean interactively as before. After this scale was cleaned, the overall residuals resembled noise. The map had a FWHM Gaussian beam size of 0.045 ×0.041 . We convolved the image to a circular 0.0475 beam. The final cleaned 12-m map is shown in Figure 1 (bottom left), which has an rms of 0.2 mJy beam −1 (0.8 K) in regions away from emission.

SOUTH-WEST STREAMER DETECTED IN DUST CONTINUUM EMISSION
In Figure 1 (top), we detect the SW streamer in 350 GHz dust continuum emission for the first time. The SW streamer is the brightest feature of the molecular outflow component and corresponds to the SW edge of the outflow cone. It was first detected in CO 1 − 0 by Bolatto et al. (2013) and has subsequently been detected in other CO transitions and dense molecular gas tracers (e.g., Walter et al. 2017;Zschaechner et al. 2018;Krieger et al. 2019). High line ratios of HCN/CO in the SW streamer indicate that this component of the outflow originates from the central starburst (Walter et al. 2017). The SW streamer has an estimated age of ∼ 1 Myr (Walter et al. 2017), in good agreement with the approximate ages of the massive star-forming regions in the starburst nucleus (see Section 5.4).
We compare the location of the dust component of the SW streamer to the CO 3 − 2 from Krieger et al. (2019) in Figure 2. The CO 3 − 2 data have a similar beam size (0.17 = 2.88 pc) as the dust map (0.15 = 2.55 pc; Figure 1, top right). The dust continuum emission has a similar morphology as the CO 3 − 2, but is offset to the southwest. In the context of the larger-scale outflow, this means that the dust is found towards the outer edge of the outflow cone.
That the dust emission is primarily at the edge of the outflow cone may be an optical depth effect. We would expect the dust to be distributed as a hollow cone, like the outflowing molecular gas, which confines the hot outflowing material (e.g., Leroy et al. 2015b;Meier et al. 2015). In projection, the dust will have the highest optical depth along the line of sight in a streamer-like structure at the very edge of the projected outflow. This effect is illustrated in Figure 8 of Bolatto et al. (2021, see especially the purple regions of this figure). It is still unclear, however, precisely why there is such a large offset between the dust and CO 3−2 in the SW streamer in NGC 253. While both CO and dust show temperature effects, the dust is more likely to remain a simple optically thin column density tracer and hence may better trace the true "spine" of the cone.

Inferred H 2 Column Density and Mass in the SW Streamer
We estimate the flux of dust emission in the SW streamer using the 12-m+ACA maps within the blue ellipse shown in Figure 1 (top right). The flux density in this region is ≈ 540 ± 180 mJy. We convert this flux density to an estimated average H 2 column density and mass, assuming a dust-to-gas ratio of 1/100, a dust mass absorption coefficient of 1.9 cm 2 g −1 following Leroy et al. (2018). We assume a dust temperature of 34 K (Gao & Solomon 2004;Weiß et al. 2008;Mangum et al. 2013), but we note that the dust temperature is not well constrained in the outflow itself. We adopt a minimum dust temperature of 11 K (Zschaechner et al. 2018) and a maximum dust temperature of 50 K (Walter et al. 2017), which we propagate into our uncertainty on the column density and mass. This calculation yields N H2 ∼ (4 +6 −2 ) × 10 23 cm −1 and M H2 ∼ (1.7 +2.6 −0.8 ) × 10 7 M . The uncertainties are dominated by uncertainties on the dust temperature in the outflow. Since the outflow is expected to be warm (e.g., Leroy et al. 2015b;Walter et al. 2017;Zschaechner et al. 2018), the true column density and molecular gas mass are more likely to be N H2 ∼ (2 − 4) × 10 23 cm −1 and M H2 ∼ (8 − 17) × 10 6 M Walter et al. (2017) estimated the H 2 column density in the SW streamer in two ways: from the CO 1 − 0 and from the Hα/Paβ line ratio. From the CO, they found N H2 ∼ 4 × 10 21 cm −2 , which is consistent with their extinction-based estimate of N H2 ∼ 5 × 10 21 cm −2 . They note, however, that the detection of bright emission from HCN and other molecules in the SW streamer implies a larger density of N H2 ∼ 5×10 22 cm −2 , in better agreement with our estimate. Our dust-based estimate of the H 2 mass in the SW streamer is consistent with the minimum mass of the SW streamer of ∼ 10 6 M found by Walter et al. (2017) and other CO-based measurements (i.e., Bolatto et al. 2013;Zschaechner et al. 2018;Krieger et al. 2019).

Width of the SW Streamer in Dust and CO 3 − 2
We estimate the width of the dust emission in the SW streamer, summing the emission along the major axis of the streamer (PA ≈ 140 • ). This yields a profile of the summed intensity along the minor axis of the SW streamer (red histogram in Figure 3). We fit a two component Gaussian to this width profile, with one narrow component for the streamer and a broad component for any disk emission. We remove the beam in quadrature from the width of the narrow component and show this beam-deconvolved narrow component corresponding to the outflow in Figure 3. The dust streamer has a beamdeconvolved FWHM of 8 pc.
We compare the width of the dust streamer to that of the CO 3 − 2. First, we kinematically identify the com-ponents of the CO 3 − 2 emission associated with the SW streamer, using a method similar to that of Walter et al. (2017) for the CO 1 − 0 (see their Section 3.2). We take position-velocity slices through the CO data cube over the field-of-view (FOV) shown in Figure 2. Each slice is 10 pc wide along the major axis of the SW streamer. Pixels within each slice are summed to produce the spectra shown in Figure 4 (left). Away from the midplane, the CO 3−2 has two velocity components, where the blueshifted component traces the outflow and the redshifted component primarily traces emission from the central starburst and disk. We fit the outflow component with a Gaussian at each offset, where the mean velocity and FWHM are indicated by the vertical and horizontal lines for each spectrum in Figure 4 (left). Using these fits, we calculate the integrated CO 3 − 2 intensity over the FWHM velocity range of the outflow component for each slice. This integrated intensity map is shown in Figure 4 (right) where the color coding of the image indicates the offset from the midplane as in the left panel. We overplot contours of the dust continuum emission to again highlight the offset between the dust and CO in the SW streamer. We measure the width of the CO 3 − 2 SW streamer using this integrated intensity map following the same procedure as for the dust described in the previous paragraph. We note that the beam sizes of the 350 GHz continuum and CO 3 − 2 are very similar (2.55 pc and 2.88 pc respectively). The width profile and Gaussian fits for the CO 3 − 2 are shown in Figure 3 (gray). The beamdeconvolved FWHM of the CO 3 − 2 associated with the SW streamer is 13 pc. The SW streamer seen in the dust continuum is narrower than in the CO by a factor of ∼ 1.6. The peak of the dust continuum is offset to the southwest of the peak of the CO 3 − 2 by 8 pc (∼ 3× the FWHM beam size).

CLUSTER SIZE, FLUX, AND MASS MEASUREMENTS
Accurate sizes and masses for the SSCs are crucial to understand their physical properties and compare to numerical simulations (e.g., Grudić et al. 2021). Previous studies which have measured these parameters used data which marginally resolved the SSCs (e.g., Leroy et al. 2018;Rico-Villas et al. 2020;Mills et al. 2021). However, using 0.5 pc resolution data, Levy et al. (2021) showed that the SSCs seen in the 350 GHz dust continuum emission break apart into multiple components once they are spatially resolved. Those data, however, used only the most extended ALMA configurations and hence lacked the short spacings which are sensitive to extended emission. This short spacing information is important to accurately measure the flux of the clusters as well as the background of extended emission in which they are embedded. Accurate measurements of the SSC sizes and masses, therefore, require both high spatial resolution and complete sampling of the Fourier plane.
As described in Section 2, in this work we combine three configurations of 12-m data from ALMA. These data have a resolution of 0.8 pc, which allows us to spatially resolve the SSCs, and a maximum recoverable scale of 66 pc, allowing us to measure the flux of the SSCs and background emission. We, therefore, remeasure the cluster positions and sizes using the 0.8 pc resolution map (Figure 1 bottom left) as described below.

Identifying the Continuum Sources
From the 0.48 pc resolution continuum data, many of the candidate SSCs identified by Leroy et al. (2018) at 2 pc resolution break apart into smaller structures (Figure 1 bottom right; Levy et al. 2021). We find more than three dozen dust clumps by-eye in the 0.48 pc resolution dust image. The SSCs identified by Leroy et al. (2018) remain the largest and brightest structures. We, therefore, follow the SSC nomenclature of Leroy et al. (2018), but add letters to sources that break apart in order of decreasing brightness, as described by Levy et al. (2021). From there, we match the locations of the dust clumps from the 0.48 pc resolution image to the 0.81 pc resolution map (Figure 1 bottom left). The 0.81 pc resolution map combines three configurations of ALMA data and, therefore, better recovers the extended emission than the 0.48 pc resolution map. From the 0.81 pc map, we identify 33 clumps of dust emission, which are listed in Table 1 and are shown in Figure 5. Some of the very small sources previously identified in the 0.48 pc resolution image are no longer visible in the 0.81 pc resolution image due to the slightly lower resolution and extended emission (i.e., SSCs 1c, 3c, 5c, 7b, 9b). For SSCs 7a and 9a, we retain the "a" lettering to indicate that these clusters do break apart in the 0.48 pc resolution image, though these smaller clusters are not visible in the 0.i81 pc resolution map. Clusters without letters (SSCs 2, 6, 14) do not break apart even in the 0.48 pc resolution dust continuum map.  Figure 4. Left: CO 3 − 2 spectra in slices along the SW streamer. The CO 3 − 2 data are binned in 10 pc bins along the major axis of the outflow, as indicated by the colorbar. The data are summed in each bin and normalized based on the maximum intensity. The spectra are offset artificially along the y-axis. The velocity component corresponding to the outflow is fit with a Gaussian, where the mean and FWHM velocity are indicated by the vertical and horizontal line segments. Right: The SW streamer seen in CO 3 − 2 (colored background) and the 350 GHz dust emission (grayscale contours), rotated so the major axis is vertical and the top of the image is closest to the midplane. For the CO 3 − 2, data in each 10 pc bin is color coded based on its offset from the midplane according to the colorbar in the left panel. We find the integrated intensity over the channels within one FWHM from the central velocity (i.e., over velocities indicated by the horizontal line segments in the left panel). The dust continuum contours are the same as in Figure   Intensity (K) Figure 5. The continuum sizes of the SSCs over the high resolution combined 350 GHz continuum map. The map has been rotated counterclockwise by 42 • so that the SSC structure is horizontal. The colored circles show the deconvolved radii from the Gaussian fits to the radial profiles (Table 1). The cluster groups are labeled and colored following the nomenclature of Leroy et al. (2018). Clusters that break apart in the higher resolution images are denoted with letters in order of decreasing brightness. To measure the precise centers of the SSCs, we fit the continuum intensity map with a 2D rotated elliptical Gaussian function. We include a constant background component since the clusters are embedded within more extended emission. Before fitting, we mask out other sources in the images. This is especially important for clusters in crowded fields. We automatically mask out primary clusters based on their half-flux radii (r half−flux ) measured from the high resolution data , removing pixels within 2×r half−flux from the cluster centers. We remove contaminating subclusters byeye and remove pixels within 1.5× the beam half-widthhalf-maximum (HWHM) from the cluster centers.

Cluster Positions
For some of the weaker clusters, the 2D elliptical Gaussian fit does not converge. In these cases, we determine the center position based on the brightest pixel in the dust continuum near the center of the SSC. The best-fitting SSC centers are listed in Table 1 and shown in Figure 5. We estimate that the positional accuracy of this image is ≈ 2.5 milliarcseconds (0.04 pc) 2 .

Radial Profiles
We construct radial profiles for each cluster. Before extracting the radial profiles, we mask the images in the same way as for the 2D Gaussian fitting. We use concentric circular 3 annuli centered on the R.A. and Decl. from the 2D Gaussian fitting (Table 1). The width of the annuli is the beam HWHM (0.024 = 0.40 pc) and the last ring has a radius of 3× the beam FWHM (3 × 0.0475 = 2.4 pc). We measure the median intensity in each annulus, which is shown in Figure 6 for SSC 14; the uncertainty is the standard error in each annulus.
We model the cluster radial profiles using a Gaussian of the form SB(r) = ae − r 2 2σ 2 + c (1) to model the radial surface brightness (SB) profile of the clusters. We also include a constant background component (c) since the clusters sit in an extended background of dust emission. An example of this model fit to the cluster radial profile is shown for SSC 14 in Figure 6, 2 See Section 10.5.2 of the the most recent version of the ALMA Technical Handbook: https://almascience.nrao.edu/ documents-and-tools/cycle9/alma-technical-handbook. Since the SSCs in the image have SNR 20, the positional accuracy is ≈ 5% of the synthesized beam. We note that the actual positional accuracy may be a factor of ≈ 2 poorer than this value due to degradation of atmospheric phase stability in the most extended configurations. 3 From the 2D Gaussian fitting, the median axis ratio of the rotated elliptical Gaussian fits is 0.9, so clusters only deviate from circular by 10%. This means that extracting the radial profiles in circular annuli (rather than in ellipses) will not introduce major systematic errors. where the fitted background level (c) has been removed. We determine the uncertainty on the Gaussian fit using a Monte Carlo simulation where we randomly vary the data points within the errorbars. The dark blue shaded region in Figure 6 reflects the standard deviation of 500 trials.
In addition to a Gaussian function, we also modeled the radial profiles using King (1962) and Plummer (1911) profiles. These provide equally good fits to the cluster radial profiles. We proceed using the Gaussian fits to the radial profiles.

Deconvolved Sizes, Fluxes, and Gas Masses
We deconvolve the beam size from the fitted Gaussian profile by removing the (Gaussian) beam HWHM in quadrature. We produce deconvolved Gaussian radial profiles using the deconvolved radii and conserving the flux. An example of the deconvolved Gaussian radial profile is shown for SSC 14 in Figure 6. We report the deconvolved cluster radii (r deconv ), peak intensities (I peak ), and fluxes in Table 1. For I peak , we subtract the background level so this value reflects the peak intensity of the cluster above the background of surrounding material. We show the distribution of intrinsic cluster radii in Figure 7. Our intrinsic radii cover a narrow range of . The light blue histogram shows the distribution of the intrinsic SSC radii, measured from the beam-deconvolved Gaussian fit to the radial profile. The left y-axis shows the histogram normalized to unit area, whereas the right y-axis shows the number of SSCs in each bin. The vertical blue dashed line shows the beam HWHM. The gray KDE shows the distribution of star cluster intrinsic effective radii measured in LEGUS galaxies normalized to unit area over the plotted radius range (Brown & Gnedin 2021). We select clusters with ages < 2 Myr and stellar masses > 10 4 M to be most comparable to our sample of very young, massive clusters. The inset in the upper right shows the full distribution of the LEGUS radii with the same selection criteria, which peaks around radii of 2 − 3 pc. The vertical black line marks 1 pc, the radial extent shown in the main panel.
radii from ≈ 0.25 − 0.70 pc. We also show the deconvolved radii for each of the clusters in Figure 5 over the continuum image.
We compare our measured cluster radii to those measured by Brown & Gnedin (2021) from the Legacy Extragalactic UV Survey (LEGUS) survey (Calzetti et al. 2015). These clusters are identified in 31 nearby galaxies in five bands from the near-UV to near-IR. Brown & Gnedin (2021) measure the intrinsic stellar half-light (effective) radii of the young star clusters in these galaxies from the "white light" (i.e. combined 5-filter) images. From their cluster catalog 4 , we select clusters with reliable radius and mass measurements, ages ≤ 2 Myr, and stellar masses ≥ 10 4 M ; see Brown & Gnedin (2021) for the definitions of these quantities. We show a kernel density estimator (KDE) of the LEGUS cluster radii in gray in Figure 7, where the inset shows the KDE over their full radius range. The LEGUS clusters tend to be larger than the SSCs studied here. The peak in the LEGUS radius distribution for clusters with the above selection criteria is between 2 − 3 pc.
It is perhaps not unexpected that the clusters identified by Brown & Gnedin (2021) are larger. The radii we measure for the clusters in NGC 253 correspond to the size of the dust (and molecular gas) envelopes, whereas the radii measured for the LEGUS clusters come from the stellar light. The clusters in NGC 253 are still in the process of forming (Leroy et al. 2018;Rico-Villas et al. 2020;Mills et al. 2021) and are, therefore, still very compact. Since the LEGUS clusters are typically older than most of the SSCs in this work and are no longer (deeply) embedded in their natal molecular clouds, it is possible that the stellar light would extend to larger radii than the compact dust emission from the SSCs. Simulations of star cluster evolution show an increase in the radius with age due primarily to mass loss from stellar winds of young massive stars (e.g., Portegies Zwart et al. 2010, and references therein).
We estimate the gas masses of the clusters based on their 350 GHz dust continuum emission, following Leroy et al. (2018, see their Section 4.3.3 for more details). Assuming a fiducial dust temperature of T dust = 130 K, we convert the deconvolved peak intensity at 350 GHz to a dust optical depth via where I 350 GHz is the peak intensity in Table 1 and B ν (T dust ) is the Planck function evaluated at 350 GHz for our adopted value of T dust . Though we do not yet have measurements of the dust temperatures towards these SSCs, the peak intensities of 25 K we measure towards some of these clusters (Table 1) supports a high value of T dust . We convert the dust optical depth to a gas surface mass density where where the dust-to-gas ratio (DGR) is assumed to be 1/100. The central 300 pc of NGC 253 is known to have a somewhat super-solar gas-phase metallicity (Z = 2.2Z ; Galliano et al. 2008;Davis et al. 2013). κ 350 GHz is the mass absorption coefficient; we assume a value of 1.9 cm 2 g −1 , but this value is uncertain by a factor of ∼2 (Ossenkopf & Henning 1994;Leroy et al. 2018). Finally, we convert the gas mass surface density to a gas mass by multiplying by the area of the cluster: where is the area of a 2D Gaussian whose HWHM is equal to the beam-deconvolved radius measurement (r deconv ; Table see 1).

Flux Density and Gas Mass Distributions
From the fluxes (Table 1), we construct the cluster flux density function, which is shown in Figure 8 (left). This is essentially a cumulative distribution function (CDF), where the ordinate counts the number of clusters with flux densities larger than the value on the abscissa. The horizontal error bars come from the uncertainties on the measured fluxes. To determine the vertical error bars, we use a Monte Carlo calculation, allowing the measured fluxes to vary uniformly within their uncertainties. This can change the ordering of the flux densities and hence the CDF. We perform 100 trials of the Monte Carlo, and report the standard deviation of the CDF for each point over those trials as the vertical error bars. We repeat this procedure with the gas masses which are shown in Figure 8 (right).
Both distributions appear to follow a broken power law, with a break at the very low flux/mass end; however, this break is very likely due to incompleteness at the low flux end of the distribution (e.g., Emig et al. 2020). The SSC identification is done by-eye based on the 0.5 pc resolution continuum image (Figure 1 bottom right). The radial profiles are measured from the 0.8 pc resolution image (Figure 1 bottom left) at the locations of the SSCs identified from the high resolution image, as long as they are still apparent in the lower resolution image. The major uncertainty matching the SSCs between these images are from the "speckles" seen in the 0.8 pc resolution continuum image ( Figure  1 bottom left). These speckles arise in the imaging by modeling the extended emission as Gaussians matched to the beam size, and they can resemble small, compact clusters. We, therefore, use these speckles as our test particles to evaluate the completeness and confusion of our SSC flux and gas mass functions. We choose a speckle in the image and measure its radial profile, beam-deconvolved size, and peak intensity as described in Sections 4.3 and 4.4. We find that a typical speckle has a flux of ≈ 3.5 mJy and hence an inferred gas mass of ≈ 5×10 4.5 M . We show these values as the gray shaded regions in Figure 8. Clusters with fluxes or masses near or below this limit are likely more uncertain than represented by the error bars, and we may be missing SSCs in this flux and mass regime.
For values above our completeness threshold, we fit the cluster flux density and gas mass functions using a broken power law of the form which is implemented using BrokenPowerLaw1D from Astropy. As a comparison, we also fit the cluster flux density function with a single power law of the form which is implemented using PowerLaw1D from Astropy. We obtain the uncertainties on the fitted parameters using the same Monte Carlo approach described above. As the flux density or mass points are allowed to vary within their uncertainties, we re-fit Equation 6 at each iteration. The uncertainties listed in Figure 8 reflect the 16 th and 84 th percentiles (i.e., the inner 68%) of the parameter distributions after the Monte Carlo. The same strategy is used to obtain the uncertainties on the model curves, shown as the colored shaded regions in Figure 8.
When the completeness is accounted for, both the broken and single power law fits are able to reproduce the flux and gas mass distributions equally well (χ 2 r = 0.3 − 0.5; Figure 8). We find a single power law slope of 1.25 (1.09) for the flux (gas mass) functions. When we fit a broken power law, we measure a slope of 2.01 (1.58) at the high flux (gas mass) end.
Previous literature studies typically investigate the cluster stellar mass function. Recently, Mok et al. (2019Mok et al. ( , 2020 studied the cluster stellar mass function of a sample of star-forming galaxies. The young clusters included in their studies are older on average and span a wider range of stellar masses than the SSCs in NGC 253 (Leroy et al. 2018;Mills et al. 2021). Mok et al. (2019) found no evidence of a high-mass cut-off of the cluster stellar mass functions though other studies have found evidence of a high-mass cut off (e.g., Gieles et al. 2006;Whitmore et al. 2010;Adamo et al. 2015Adamo et al. , 2017Hollyhead et al. 2016;Johnson et al. 2017;Messa et al. 2018a,b).
Although our measurements are of the gas mass (not the stellar mass), we also do not see strong evidence for a high-mass cut-off up to 2 × 10 5 M , which would be indicated by a more apparent break in the gas mass distribution (Figure 8 right). While our two highest mass clusters may hint at a break around 10 4.9 M , statistically the single and broken power law fits are equally good (as indicated by χ 2 r ). This break (or truncation) mass is similar to what has been measured in other nearby galaxies (see e.g., discussion and references in Messa et al. 2018a   From their fits using a minimum mass of 10 4 M (similar to our completeness limit), Messa et al. (2018a) found a single power law slope of 2.67 ± 0.03, steeper than found by Mok et al. (2020). The power law slope is reduced to 2 (within the uncertainties) if the power law is truncated and if the minimum mass is reduced (see their Table 8).
The average slopes for the youngest clusters found by both Mok et al. (2020) and Messa et al. (2018a) are steeper than, but close to, the slope we measure for the gas mass distribution in NGC 253. We caution, however, that we are measuring the gas mass whereas Mok et al. (2019Mok et al. ( , 2020 and Messa et al. (2018a) measured the stellar mass of the clusters. Both Leroy et al. (2018) and Mills et al. (2021) found a median M gas /M * ≈ 1 for the clusters, but with appreciable scatter (average ≈ 2, standard deviation ≈ 2.5) using different tracers of the stellar and gas masses. This scatter from cluster to clus- 5 We note that the galaxy samples used by Mok et al. (2019Mok et al. ( , 2020 overlap but are not the same.
ter means that the slope of the stellar function function in NGC 253 may be different than measured for the gas mass.

THE MORPHO-KINEMATIC ARCHITECTURE OF THE SSCS
The quasi-linear arrangement of the SSCs in the center of NGC 253 is striking (e.g., Figures 1 and 5). In projection, this structure measures ∼ 155 pc × 15 pc in diameter with a major axis position angle (PA) of ≈ 48 • east-of-north. This axis ratio of ∼ 10 may suggest that the structure is intrinsically very thin and/or that we are seeing this structure at a high inclination. The galactic disk of NGC 253 is nearly edge-on, with an inclination of ≈ 78 • and a PA of ≈ 50 • (e.g., Pence 1980;Westmoquette et al. 2011;Krieger et al. 2019). On ∼ kpc scales, the bar also has an inclination of ≈ 78 • but has a major-axis PA of 68 • (Scoville et al. 1985;Sorai et al. 2000). The quasi-linear arrangement of SSCs has approximately the same PA as the galaxy disk and is offset from the PA of the bar.
Using CO observations at 35 pc resolution, Leroy et al. (2015a) measure the geometry of the GMC structures in which these SSCs are embedded. They build 3D models of the GMC geometry as a disk, a linear bar-like arrangement, and a hybrid model. They find that the hybrid model provides the best fit to the data, where the inner ∼ 100 − 150 pc (diameter) is more disk-like and regions beyond this extending out to ∼ 850 − 1400 pc (diameter) have a more linear structure. They find that the maximum vertical thickness of the GMC structure is < 100 pc for the molecular gas traced by CO and < 55 pc for the denser molecular gas. Our measurement of the minor axis width of the SSC structure sets a maximum vertical extent of < 15 pc, similar to the vertical extent of the MW CMZ (e.g., Molinari et al. 2011;Kruijssen et al. 2015;Shin et al. 2017;Henshaw et al. 2022).

Bar Resonances
Zooming out in scale, Sorai et al. (2000) used CO 1−0 data with 16 (∼270 pc) resolution from the Nobeyama 45-m telescope to study the bar resonances and kinematics in the center of NGC 253. By fitting the CO rotation curve, they were able to estimate the locations of several key bar resonances. From their CO rotation curve, they measure the angular velocity, Ω(R), and the epicyclic frequency, κ(R). We note that Sorai et al. (2000) used an older distance for NGC 253 of 2.5 Mpc; here we have converted all of their measurements assuming a galaxy distance of 3.5 Mpc 1 . Sorai et al. (2000) used previous estimates of the pattern speed of the bar, Ω p ≈ 35 km s −1 kpc −1 . The radius at which Ω p = Ω is called the corotation radius (CR) and occurs at a radius of 5.6 kpc in NGC 253 (green ellipse in Figure 9). Other resonances occur at harmonics of κ. For example, the outer and inner inner Lindblad resonances (OILR and IILR respectively) occur where Ω p = Ω − κ/2.
These resonances are important because gas is expected to collect inside these resonance locations, with gas inside the IILR possibly collapsing to trigger the nuclear starburst (e.g., Goldreich & Tremaine 1979;Wada & Habe 1992;Sorai et al. 2000;Paglione et al. 2004). In NGC 253, the OILR is located at a radius of 1.8 kpc, well matched to the extent of the bar where it intersects with the circumnuclear "2 kpc" ring ( Figure 9). The IILR is located at a radius of 336 pc.
We caution, however, that the location of the IILR determined by Sorai et al. (2000) is especially uncertain because it is on the same scale as their spatial resolution. Moreover, Sorai et al. (2000) determined their CO rotation curve by finding the terminal velocities (see e.g., Section 3.3 of Sofue & Rubin 2001). While this method is often used for highly inclined systems, rotation velocities within the central few resolution elements (especially where there are steep velocity gradients) are especially uncertain (e.g., Sofue & Rubin 2001) Moreover, the central bar will produce strong non-circular motions that are not taken into account in the terminal velocities method. This can produce dramatic artifacts in the determination of the rotation curve (see e.g., Section 5.3 of Binney et al. 1991). Uncertainties on the CO rotation curve will propagate into the determination of the resonance locations. Therefore, all of the resonance locations reported by Sorai et al. (2000) are likely uncertain, with the IILR being the most uncertain.
As an example of the uncertainty in the position of the IILR in NGC 253, we compare the IILR location measured by Sorai et al. (2000) to that inferred by Das et al. (2001). Das et al. (2001) model the central regions of NGC 253 using a logarithmic bar potential (see their Equation 1). From this potential and the values in their Table 1, we calculate Ω(R) and κ(R) to determine the CR, OILR, and IILR from this model. From this, we find that the CR is 5.8 kpc, the OILR has a radius of 1.8 kpc, and the IILR has a radius of 27 pc. While the locations of the CR and OILR are similar between Sorai et al. (2000) and Das et al. (2001), the IILR differs by an order of magnitude. Therefore, because the position of the IILR is so uncertain, we should not attach particular meaning to its location relative to the SSCs and dense molecular gas in NGC 253.

Families of Orbits in a Barred Potential
There are families of closed orbits in a bar potential. The x 1 (bar) orbits are extended along the major axis of the bar whereas the x 2 (antibar) orbits are perpendicular to the bar major axis (e.g., Contopoulos & Mertzanides 1977;Athanassoula 1992a,b;Binney & Tremaine 2008). The x 2 orbits are closely related to the ILR (or an OILR and IILR; e.g., van Albada & Sanders 1982;Athanassoula 1992a). At the intersections between these two orbital families, gas can collide and shock, lose angular momentum, and transition from the x 1 to the x 2 orbits, bringing it closer to the galactic center. Das et al. (2001) calculate the x 1 and x 2 orbit families in NGC 253 assuming a logarithmic bar potential. In Figure 9 (bottom) we show the outermost and innermost x 1 orbits (red) and the outermost and innermost x 2 orbits (pink), for a non-axisymmetry parameter of 0.8 (see also Figure 3 of Das et al. 2001). The x 2 orbits intersect the x 1 orbits and extend down to ∼100 pc scales. Therefore, these orbits can facilitate the transfer of gas from large to small scales, fueling the nuclear starburst.

Molecular Gas
From the discussion above, it is interesting to determine how the SSCs are arranged with respect to the bar orbits and resonances. Due to the nearly edge-on inclination, constraining the arrangement purely from the cluster locations is difficult. A similar challenge is faced in studying the MW CMZ, where the massive star forming regions are viewed nearly edge-on and thought to be arranged as spirals, streams with either open or closed orbits, or a ring (e.g., Sofue 1995;Sawada et al.  200 pc x 1 x 2 Figure 9. NGC 253 as seen in the J-band from 2MASS (grayscale). The images on the left are 16 kpc (15.7 ) on a side; the images on the right are 3 kpc (2.9 ) on a side. The gray squares mark the inner 280 pc, the FOV of the ALMA data. In the right column, the gray × shows the kinematic center determined by Anantharamaiah & Goss (1996). Top row: Overplotted are the locations of important resonances and features, determined by Sorai et al. (2000). Working inward, the solid ellipses show the CR (green), the OILR (dark blue), and the IILR (light blue). We note that the location of the IILR is highly uncertain (see the discussion in Section 5.1). The gray circles in the lower right corners show the 16 (∼270 pc) beam of the CO 1 − 0 observations used by Sorai et al. (2000). Bottom row: The outermost and innermost x1 orbits are shown in red, and the outermost and innermost x2 orbits are shown in pink (Das et al. 2001).
2004; Molinari et al. 2011;Krumholz & Kruijssen 2015;Kruijssen et al. 2015;Henshaw et al. 2016;Ridley et al. 2017;Tress et al. 2020;Sormani et al. 2022). Henshaw et al. (2016) and Henshaw et al. (2022) provide excellent reviews and testing of these models in the CMZ (see especially Figures 18 and 19 of Henshaw et al. 2016). Without knowledge of how the SSC structure is dynamically linked to the bar (or not), it is impossible to tell whether the SSC structure is tilted in the same manner as the galaxy disk (i.e., where the near side of the disk is towards the northwest). We can, however, use the kinematic information from the cluster velocity measurements presented by Levy et al. (2021) together with the morphology to constrain the cluster arrangement.

SSC and Dense Molecular Gas Kinematics
In Figure 10, we show the locations and systemic velocities of the primary SSCs. The radius of each circle shows the deconvolved SSC size (r deconv , see Table  1 and Section 4.3). The systemic velocities were measured by Levy et al. (2021) using a multi-Gaussian fit to many spectral lines detected towards these clusters. The uncertainties on the SSC systemic velocities are better than ±5 km s −1 . The velocity color scale is centered on the systemic velocity of the galaxy (250 km s −1 ; Müller-Sánchez et al. 2010;Krieger et al. 2019Krieger et al. , 2020b. We compare the distribution and kinematics of the primary SSCs to CS 7 − 6 observations from Krieger et al. (2019Krieger et al. ( , 2020a in Figure 10. Their cleaned cubes are a combination of 12-m, 7-m, and total power data and have a spectral resolution of 2.2 km s −1 and a spatial resolution of 0.17 × 0.13 (2.9 pc × 2.2 pc). We fit the CS 7−6 line with a Gaussian at each pixel. We mask the velocity map based on the peak intensity, where pixels with a signal-to-noise ratio (SNR) < 5 are removed. The top panel of Figure 10 shows a 3D position-positionvelocity (R.A.-Decl.-Velocity) diagram of the CS 7 − 6 data and the primary SSCs. The bottom row of Figure  10 shows R.A.-Decl. and position-velocity projections. The panels in Figure 10 are cropped to the FOV of the ALMA data of the SSCs described in Section 2, which cover the central 280 pc of the galaxy.
The positions and velocities of the SSCs agree well with the CS 7 − 6 emitting dense molecular gas. Overplotted in the R.A.-Decl. panel of Figure 10 are the outermost and innermost x 2 orbits (pink dashed ellipses; see also Figure 9). The SSC structure and the dense molecular gas in which they are embedded are remarkably well aligned with the innermost x 2 orbits, especially given the spatial resolution of data with which the x 2 orbits were derived (≈300 pc; Das et al. 2001), though there appears to be a slight PA offset. We note that the CS 7 − 6 data include short-and zero-spacing data and so recover emission on large spatial scales. This means that the concentration of CS 7 − 6 emission around the innermost x 2 orbit is not a result of spatial filtering by the interferometer. This is compelling evidence that the SSCs and the dense molecular gas from which they form are located at the nexus of the family of orbits which transfer gas from the bar on large scales to small scales where the gas can become dense enough to form massive young stellar clusters.

3D Structure of the SSCs and Dense Molecular Gas
Nuclear rings following x 2 orbits are formed spontaneously in a barred potential, as revealed by hydrodynamic simulations (e.g., Tress et al. 2020), and are able to power galaxy-scale outflows (e.g., Nguyen & Thompson 2022). x 2 orbits are mildly elongated perpendicular to the bar and have a nearly-elliptical shape. Their orbital velocity is larger at the pericenter and lower at the apocenter, qualitatively similar to what one would get by assuming that the angular momentum is constant along the orbit. Although the angular momentum is not exactly conserved along x 2 orbits (because a bar potential is non-axisymmetric -it oscillates around a mean value), a reasonable approximation is to model x 2 orbits as elliptical orbits on which the angular momentum is conserved (see also Peters 1975). This type of model has been applied to the MW CMZ and can explain the arrangement of dense molecular gas features and star forming regions ). In the MW, several other types of models have been developed to explain the orbits of gas, stars, and massive star-forming regions in the CMZ, including twisted rings, spirals, and crossing streams (see e.g., Henshaw et al. 2022, and references therein). We apply some of these models to the CMZ of NGC 253 in Appendix A.

A Plausible Model
Here, we construct a simple kinematic model of the nuclear ring in NGC 253 by assuming that the gas follows elliptical orbits on which the angular momentum is conserved. This model is built in Cartesian coordinates where the x-axis corresponds to galactic longitude, the y-axis to the line of sight from the observer, and the z-axis to galactic latitude (Figure 11, top). The elliptical, angular momentum-conserving ring is described by four parameters: the semi-major axis length (a), the semi-minor axis length (b), the orbital velocity at the pericenter of the ring (V orb,0 ), and the position angle of the major axis of the ellipse with respect to the x-axis  Figure 9 and the FOV of the ALMA data). Pixels with SNR < 5 in peak intensity are masked out in the CS 7 − 6 velocity map. The positional offsets are calculated based on the center determined by Anantharamaiah & Goss (1996, shown as the gray × assuming this point has a velocity equal to the galaxy systemic velocity). The color scale is the same in all panels and shows the LSRK velocity. The radii of the SSC markers are 2 × r deconv (Table 1) the ellipse is determined by conservation of angular momentum starting from V orb,0 . The elliptical orbit in this frame is shown in the top panel of Figure 11. To compare this model with the SSCs, we project it into the sky plane assuming some position angle (PA; defined counterclockwise of north to the receding side the the ring) and an inclination (i). The projection of this ring into R.A.-Decl. coordinates is shown in the left column of Figure 11. We also construct a positionvelocity diagram, where the position is taken along the major axis of the ring given by the PA (Figure 11, right column).
We adjust the ring parameters by-eye to best fit the arrangement and kinematics of the SSCs. A ring with a ∼ 110 pc, b ∼ 60 pc, V orb,0 ∼ 115 km s −1 , and θ p ∼   Figure 11. The observer shown in the top panel of Figure 11 is located at (x, y) = (0, −∞) in this plot. The model-dependent deprojected positions of the SSCs along this ring are shown (except for SSC 7a which is not well fit by this model). Arrows indicate the direction of the orbit. Other annotations and color-coding are as in Figure 11. 55 •6 , PA ≈ 235 • , and i ≈ 85 • is a good representation of the SSCs. We reiterate that these parameters were fit by-eye and that there are degeneracies among them. This is not meant to be an exact measurement of the size of the SSC structure, but rather to show that this is a possible configuration with reasonable parameters. As shown in Figure 11, the SSCs agree very well with this simple model. In Figure 12, we show a top-down ("face-on") view of this ring. Using this model, we place the SSCs along the ring; the positions along the line of sight (y-axis) are entirely model dependent. Under this model, SSCs 1a, 2, 3a, 4a, 5a, 9a, 11a, and 12 would be on the front (near) side of the ring, shown as the blue outlines in Figures 11 and 12. SSCs 8a, 10a, 13a, and 14 would be on the back (far) side of the ring, shown as the red outlines.
The orbital period along the elliptical ring shown in Figures 11 and 12 is ≈ 6 Myr. Given the limits on the SSC ages (see Section 5.4.1), this means that the SSCs have not completed a full orbit. On average, the SSCs will complete half an orbit before the massive stars within them explode as supernovae.
SSC 7a is not well fit by this model. Although this cluster is relatively weak compared to others, its systemic velocity was well constrained by Levy et al. (2021) and its velocity agrees well with the CS 7 − 6 (Figure 10). Unlike Leroy et al. (2018) and Mills et al. (2021), we find evidence for a single velocity component towards this SSC, which may be due to the increased spatial resolution. There are multiple spatial components near SSC 7a (see e.g., Figure 1 bottom right and Figure 1 of Levy et al. 2021), which can be blended together at lower resolution leading to a second velocity component. Moreover, a shift in velocity cannot fully bring SSC 7a into agreement with the model; a change in the cluster's velocity would shift its position vertically in position-velocity diagram (Figure 11; right column) which does not bring it into agreement with the model for reasonable adjustments to the cluster's velocity. Because SSC 7a is not well described by this model, we do not show it in Figure 12. Hydrodynamical simulations of the MW CMZ find gas and star formation inside the CMZ ring, and it is possible that this star formation may be happening in the vicinity of the supermassive black hole Sagittarius A * or associated with the nuclear star cluster (Sormani et al. 2020, and references therein). In NGC253, we also detect dense molecular gas inside the elliptical ring model (e.g., Figure 11). SSC 7a could be evidence of massive star formation close to the galactic center (whose precise position is unknown) though, with only a single cluster, this is only speculative.
Along the ring, there are more extreme radial velocities than represented by the SSCs or the dense molecular gas traced by CS 7 − 6, particularly on the redshifted side of the structure (Figure 11; bottom row). In the case of the CS 7 − 6 data, this is not an effect of the interferometer filtering out emission on large scales as these data include the zero-spacing (total power) data (Krieger et al. 2019(Krieger et al. , 2020a. The CS 7 − 6 data cube covers LSRK velocities from 88 to 373 km s −1 , or -162 to +123 km s −1 about the systemic velocity of the galaxy. This could mean that we are missing some of the most redshifted CS 7 − 6 emission if it falls outside the bandpass. As a check, we examine the spectra in the CS 7−6 cube, but we do not find evidence that the line profiles are cut off by the edge of the bandpass. As noted above, simulations show that the angular momentum in an x 2 orbit is approximately constant in a time-averaged sense (e.g., Sormani et al. 2018;Tress et al. 2020). As the major-to-minor axis ratio of the orbit increases, the amplitude of the oscillations in angular momentum increase and the assumption of angular momentum conservation breaks down. For the CMZ of the MW, an x 2 orbit with an axis ratio of 110/60 ≈ 1.8 would yield variations of 25% in the angular momentum over time (Sormani et al. 2018). This means that for an axis ratio like we estimate for NGC 253, the variation in angular momentum are not negligible; this is a limitation of this simple model.
As described in Section 1, we assume that the dynamical center of NGC 253 is the center determined from ionized gas kinematics by Anantharamaiah & Goss (1996). This center position has a 1-σ uncertainty of ∼ 0.3 , which corresponds to 5 pc. This is a negligible uncertainty in terms of the model fit shown in Figure 11, so we conclude that uncertainties in the center position determined by Anantharamaiah & Goss (1996) will not affect the results from the model. Other determinations of the dynamical center of the galaxy, however, differ by more than this measurement uncertainty. For example, the center determined from stellar kinematics by Müller-Sánchez et al. (2010) is 1.3 (22 pc) away to the northeast in projection (between SSCs 8 and 10/11). From the standpoint of the modeling, varying the location of the dynamical center will change how the SSCs are distributed along the ring. If, for example, we instead choose the center determined by Müller-Sánchez et al. (2010), a larger ring with a higher orbital velocity is required. More challenging, however, is that a single (symmetric) ring cannot fit clusters on both the nearand far-sides of the ring simultaneously. In other words, compared to Figure 11 (right panels), a model using the center from Müller-Sánchez et al. (2010) can either fit clusters on the near (blue outlined) side of the ring or on the far (red outlined) side, but not both simultaneously. Therefore, if the true dynamical center of NGC 253 is that from Müller-Sánchez et al. (2010), then this model of a symmetric angular momentum-conserving ring is not a good fit to the data.

Previous Estimates of the SSC Ages in NGC 253
In terms of the relative ages of the SSCs, both Rico-Villas et al. (2020) and Mills et al. (2021) found evidence for an inside-out formation scenario, where clusters towards the ends of the structure are younger and those towards the middle are older (i.e., SSCs 1, 2, 3, 13, and 14 are among the youngest and SSCs 4, 5, 7, 8, 10, 11, and 12 are among the oldest). On the other hand, Krieger et al. (2020a) used line ratios of HCN/HC 3 N to construct a relative chemical age gradient, which does not follow a pattern with distance from the center and where the progression from youngest to oldest is SSC 13, 14, 1, 8, 3, 2.
The absolute ages of the SSCs in NGC 253 are highly uncertain, but limits can be placed on them. Rico-Villas et al. (2020) estimated that the SSCs in NGC 253 have ages ≈ 0.01 − 1 Myr based on the ratio of luminosity of in stars (from free-free emission) and protostars. These ages are likely a lower limit, however, since the ionizing photons that produce the free-free emission may be absorbed by dust within the SSCs (see Levy et al. 2021, for a further discussion). Rico-Villas et al. (2020) also found a weak trend between their estimated ages and the HNCO/CS line ratio, where younger clusters tend to have higher line ratios. Hydrodynamical simulations of the Milky Way show that the stars and dense molecular gas are well coupled until the stars are ∼5 Myr old . In NGC 253, the SSCs and CS 7 − 6 are well matched in terms of their locations and kinematics (Figure 10), indicating that the SSCs are younger than ∼5 Myr. In a separate study, Mills et al. (2021) found that the clusters in NGC 253 cannot be older than ∼ 3 Myr due to the presence of He recombination lines and the lack of appreciable synchrotron emission towards most of the clusters (indicating a lack supernovae in the SSCs). Therefore, the ages of the SSCs are likely ∼ 0.01 − 3 Myr.

Gradients in a Circumnuclear Ring
There are three main models that describe where and when star formation occurs in a circumnuclear ring. These models make predictions for age gradients (or lack thereof) along the orbit. First, the so-called "pearlson-a-string" scenario predicts that star formation occurs just downstream of the contact points between the gas inflow from the bar and the circumnuclear ring (i.e., downstream of the apocenters of the orbit; e.g., Böker et al. 2008;Mazzuca et al. 2008;Sormani et al. 2020). Moreover, gas in the ring will have the slowest orbital velocities at the apocenters, more easily allowing it to pile up and become dense. The star clusters age as they orbit along the ring, so the youngest clusters should be found near the apocenters with an increase in age downstream. Alternately, it has been suggested that star formation could be triggered when clouds are compressed due to close pericenter passages (Longmore et al. 2013;Kruijssen et al. 2015). Under this scenario, the youngest stars should be found closest to pericenter and stellar (or cluster) ages should increase downstream. Finally, the "popcorn" model describes a scenario in which the ring forms from gravitational collapse or turbulence and star formation is distributed uniformly along the ring (Böker et al. 2008). In this scenario, no age gradients are expected as either the clusters have approximately the same age or the ages are randomly distributed along the ring.
It is tempting to use the elliptical angular momentum conserving ring model to infer the relative ages of the SSC. However, the observations of the SSCs in NGC 253 are a single snapshot in time. Hydrodynamical simulations show that the instantaneous star formation distribution and ring morphology vary substantially about the time-averaged behavior (e.g., Tress et al. 2020;Sormani et al. 2020). While the time-averaged simulation results favor the "pearls-on-a-string" model, it is much less clear if signatures of the latter can be seen in the instantaneous simulation results due to large timefluctuations (see especially Figure 9 of Sormani et al. 2020). Therefore, with only this single snapshot in time and a small number of SSCs, we cannot simply use the "pearls-on-a-string" (or other) model to infer the expected cluster age gradients from their locations along the elliptical ring.

SUMMARY
The SSCs in the center of NGC 253 are bright continuum sources at 350 GHz, which primarily traces thermal emission from warm dust (e.g., Leroy et al. 2018;Mills et al. 2021). Here, we combine ALMA data from three 12-m configurations and the 7-m array to construct maps of the dust emission covering scales from 0.028 −24.8 (0.48 pc−421 pc). This enables us to measure the compact dust emission associated with the clusters themselves as well as the more extended emission in which they are embedded (Figure 1). We summarize our main results below, indicating the relevant figures and/or tables.
1. For the first time, we detect the galaxy-scale outflow in dust continuum emission. As shown in Figures 1, 2, and 4, we find dust emission along the SW streamer that is located at the edge of the CO emission and hence on the exterior edge of the outflow cone. The dust streamer has a FWHM of 8 pc.
From the dust, we estimate that the SW streamer has a molecular gas mass of ∼ (8 − 17) × 10 6 M , consistent with other measurements (e.g., Walter et al. 2017).
2. We measure the sizes of the SSCs using Gaussian fits to the cluster radial profiles. We deconvolve the Gaussian beam from the size measurements to provide beam-deconvolved cluster sizes, finding radii of 0.4 − 0.7 pc (Table 1; Figures 6 and 7). Compared to star clusters in the LEGUS survey, the SSCs in NGC 253 tend to be smaller, likely because they are younger (Figure 7).
3. We investigate the morpho-kinematic arrangement of the SSCs and their possible connection to the bar. The SSC structure is on the same scale as the x 2 orbits, suggesting that gas is transported down to these scales by the bar where it can then become dense enough to form massive star-forming regions (Figures 9 and 10). We find that the SSCs have a similar distribution and kinematics as the dense molecular gas (Figure 10). We are able to describe the SSC morphology and kinematics with a simple elliptical, angular momentum-conserving model (Figures 11 and 12), which has a semi-major axis of ∼ 110 pc, a semi-minor axis of ∼ 60 pc, and an orbital period of ≈ 6 Myr. From our perspective, this ring would appear nearly edge-on, leading to the observed nearly linear SSC distribution.
As described in Section 5.4.1, estimates of the (relative) SSC ages in NGC 253 do not all agree (Rico-Villas et al. 2020;Krieger et al. 2020a;Mills et al. 2021). Constraints on the absolute ages of the clusters are relatively weak, though the SSCs must be young (∼ 0.01 − 3 Myr; Rico-Villas et al. 2020;Mills et al. 2021). In the future, approved observations with the MIRI integral field unit onboard JWST will measure the ionizing radiation field of these clusters and hence provide independent constraints on the cluster ages.  Figure 13. Similar to Figure 11, but for the CMZ models; see the caption of Figure 11 for more details. The top row shows a twisted ring model (Molinari et al. 2011), the middle row shows two spirals (e.g., Sofue 1995; Sawada et al. 2004), and the bottom row shows a crossing streams model . Molinari et al. 2011). In the CMZ of the MW, this configuration is thought to follow the x 2 orbits and the twisted shape is induced by a vertical oscillation. The twisted ring model is described by a semi-major axis (a), a semi-minor axis (b), a constant orbital velocity (V orb ), and a vertical oscillation with an amplitude (z 0 ), frequency (ν z ), and phase (θ z ). z 0 is the amplitude of the oscillation, measured from the midplane to the peak (so the ring has a thickness of 2×z 0 . We follow Molinari et al. (2011) and assume that ν z is twice the orbital frequency and θ z = 0. Because this twisted ring model has a constant orbital velocity, it does not conserve angular momentum. This twisted ring model and the projections are shown in Figure 13 (top). The best fitting twisted ring model for the CMZ of the MW has a = 100 pc, b = 60 pc, z 0 = 15 pc, V orb = 80 km s −1 , and θ p = −40 • . While such a ring in NGC 253 would have the same major-and minor-axis lengths, the amplitude of the vertical warp is < 10 pc. Moreover, we require a higher V orb = 110 km s −1 in NGC 253.

A.2. Open Orbit Models
As pointed out by Kruijssen et al. (2015) and others for the MW CMZ, non-circular stable closed orbits are only possible if the potential is not axisymmetric at these scales. In NGC 253, structural modeling by Leroy et al. (2015a) suggests that there are asymmetries on larger scales (as is also expected for the CMZ), but their best fitting structural model is axisymmetric on the 100 − 150 pc scales we study here. We note, however, that their data are not particularly constraining on these scales due to the 35 pc resolution. Therefore, since we know the SSCs are embedded in larger-scale gas structures (e.g., Sakamoto et al. 2011;Leroy et al. 2015a;Krieger et al. 2019Krieger et al. , 2020a, it is perhaps more likely that the SSC orbits will be open. One possible open-orbit model describes streams stemming from the bar-ends. Several of these kinds of models have been proposed for the CMZ in the MW (e.g., Sofue 1995;Sawada et al. 2004;Kruijssen et al. 2015;Henshaw et al. 2016;Ridley et al. 2017;Henshaw et al. 2022). Kruijssen et al. (2015) develop a model of the MW CMZ that consists of crossing streams that form a ring-like structure. Unlike the closed ring model, this stream model is open and Sagittarius A* is located at one of the foci of the eccentric orbits. This twisted pattern precesses with time (see also Tress et al. 2020).
First we compare against a model with two spiral arms (Sofue 1995;Sawada et al. 2004). Following Henshaw et al. (2016), we select two of the streams from the Kruijssen et al. (2015) model, and we fix their vertical positions to 0 (i.e., the streams are in the midplane).
This model and the projections are shown in Figure 13 (middle).
Next, we compare against the full four-stream model developed by Kruijssen et al. (2015). To better fit the SSCs in NGC 253, we stretch this model in the xdirection by a factor of 1.5 and in velocity by a factor of 1.25, and we show this model in Figure 13 (bottom). Unlike any of the other models described here or in Section 5.3, this model does reproduce the position and velocity of SSC 7a; however, we cannot claim that this model is superior to others based on one SSC.