Hidden Gems on a Ring: Infant Massive Clusters and Their Formation Timeline Unveiled by ALMA, HST, and JWST in NGC 3351

We use 0.1″ observations from the Atacama Large Millimeter Array (ALMA), Hubble Space Telescope (HST), and JWST to study young massive clusters (YMCs) in their embedded “infant” phase across the central starburst ring in NGC 3351. Our new ALMA data reveal 18 bright and compact (sub-)millimeter continuum sources, of which 8 have counterparts in JWST images and only 6 have counterparts in HST images. Based on the ALMA continuum and molecular line data, as well as ancillary measurements for the HST and JWST counterparts, we identify 14 sources as infant star clusters with high stellar and/or gas masses (∼105 M ⊙), small radii (≲ 5 pc), large escape velocities (6–10 km s−1), and short freefall times (0.5–1 Myr). Their multiwavelength properties motivate us to divide them into four categories, likely corresponding to four evolutionary stages from starless clumps to exposed H ii region–cluster complexes. Leveraging age estimates for HST-identified clusters in the same region, we infer an evolutionary timeline, ranging from ∼1–2 Myr before cluster formation as starless clumps, to ∼4–6 Myr after as exposed H ii region–cluster complexes. Finally, we show that the YMCs make up a substantial fraction of recent star formation across the ring, exhibit a nonuniform azimuthal distribution without a very coherent evolutionary trend along the ring, and are capable of driving large-scale gas outflows.

1. INTRODUCTION Massive star clusters (M ⋆ ≳ 10 4 M ⊙ ) form in dense, turbulent, high-pressure environments (Portegies Zwart et al. 2010;Longmore et al. 2014;Krumholz et al. 2019).This process appears distinct from star formation under lower gas density conditions, with the massive clusters quickly and efficiently convert their gas content into stars before the newly formed stellar population exerts strong feedback (e.g., Krumholz et al. 2014).Understanding the formation of massive clusters can thus provide unique constraints on theoretical models of star formation and stellar feedback.It can also shed light on how the majority of stellar populations were built up earlier in the cosmic history, when high gas density and turbulence conditions were prevalent (e.g., Madau & Dickinson 2014).
To achieve a more complete, observation-grounded understanding of massive cluster formation, it is critical to probe the earliest ("infant") phase, during which star formation and feedback processes are still ongoing.This has been challenging for several reasons.The relevant physical processes happen on a short timescale (≲ a few Myr; see e.g., Krumholz et al. 2014), making it hard to catch forming YMCs in the infant phase.Those YMCs that are in this phase are still deeply embedded in their natal gas and dust "cradles" and thus remain basically invisible at short wavelengths (UV, optical, and sometimes even IR; see Kornei & McCrady 2009;Leroy et al. 2018).Last but not least, the physical conditions required for creating YMCs (high gas density, turbulence, and pressure) are rare in the Milky Way and the nearest extragalactic systems (Portegies Zwart et al. 2010).Building a sizable sample of infant YMCs thus necessitates obtaining sensitive and resolved observations for these compact objects in more distant galaxies, which can be a tall order even with the latest facilities.

Sun et al.
Recent studies based on deep, long-baseline, targeted observations with the Atacama Large Millimeter Array (ALMA) have identified a promising avenue for addressing these challenges.With its exquisite sensitivity and resolving power in (sub-)millimeter bands, ALMA can detect extinction-free tracers of star formation (via freefree continuum and radio recombination lines) as well as the associated gas reservoir (via molecular lines and dust continuum) for YMCs even at extragalactic distances.Thanks to these unique capabilities, we have already identified and characterized individual forming YMCs or entire YMC populations in our Galaxy (Schmiedeke et al. 2016;Ginsburg et al. 2018) and a handful of nearby galaxies (e.g., the Large Magellanic Cloud: Nayak et al. 2019;NGC 5253: Turner et al. 2017;NGC 253: Leroy et al. 2018;Levy et al. 2021;Mills et al. 2021;NGC 4945: Emig et al. 2020; Henize 2-10: Costa et al. 2021; and the Antennae: Finn et al. 2019;He et al. 2022).
Building on these pioneering studies, the logical next steps are to connect the infant YMCs identified in the (sub-)millimeter bands with more evolved clusters visible at shorter wavelengths, and to put them into the larger-scale context of the host galaxy.Doing so requires multiwavelength observations of a sizable YMC population at matched spatial resolution, supported by rich ancillary data for the host galaxy.It is also important that the host galaxy is at a favorable viewing angle for source localization and cross-wavelength source matching, which have proven difficult in edge-on systems including our Galaxy, NGC 253, and NGC 4945 (Stolte et al. 2014;Emig et al. 2020;Levy et al. 2022).
In this paper, we use nearly matched-resolution ALMA, HST, and JWST observations to examine a population of forming YMCs in a nearby Milky Way analog galaxy, NGC 3351 (a.k.a.M95, see Table 1 for its basic properties).This galaxy features a prominent, dusty, inner ring structure, which is likely fed by large-scale gas inflows induced by a strong stellar bar (e.g., Regan et al. 2006, see Figure 1).This "starburst ring" hosts a large number of optically-visible massive clusters, and previous IR observations show signs of embedded YMC formation as well (e.g., Ma et al. 2018;Calzetti et al. 2021;Turner et al. 2021).Its favorable viewing angle (i ≈ 45 • , Table 1) and clean orbital configurations (Figure 1) enable cross-wavelength source matching as well as examinations of systematic trends along the ring.
The structure of this paper is as follows.Section 2 describes the ALMA, HST, and JWST datasets used in this paper.Section 3 details our source identification, characterization, and cross-matching schemes.Section 4 presents a set of key physical properties measured for all ALMA-identified YMC candidates.Section 5 synthesizes the observational results from ALMA, JWST, and HST, constructs an evolutionary timeline for YMC formation, and puts the YMC population in the large-  scale context of the starburst ring in NGC 3351.We summarize all our findings in Section 6.
2. DATA We use a new set of very high resolution ALMA observations in this work.Here we lay out the observational design, reduction procedures, and data characteristics for this new ALMA dataset (Section 2.1).We then briefly describe ancillary HST and JWST images and data products used in this work (Section 2.2).
Our Band 7 data combine C-6 and C-3 configurations of the 12-m array and supplementary 7-m array observations (all from the same Cycle 8 project) to achieve a similar resolution and MRS as our Band 3 observations.The 12-m observations cover a ∼ 30 ′′ FoV using a 7pointing mosaic and reach on-source integration times of 0.4 h (C-6) and 0.2 h (C-3); the 7-m observations use  Lee et al. 2022).Right: JWST/MIRI composite image for the same field (F1000W+F1130W+F770W; Lee et al. 2023).In each panel, a white square marks the central 20"×20" area, which covers the starburst ring and is roughly the field-of-view of our ALMA observations.a 3-pointing mosaic to cover a similar FoV, with 0.5 h of on-source integration.The Band 7 spectral tuning covers the 342-357 GHz continuum and the CO (3-2), HCO + (4-3), and CS (7-6) lines at a native velocity resolution of 0.95 km s −1 .

Calibration and Imaging
We calibrate the raw visibility data with observatorysupplied scripts and the appropriate version of CASA pipeline (6.2.1 for Cycle 8).The calibrated data show no obvious pathologies upon visual inspection.
From the calibrated measurement sets, we extract and image a relevant subset of visibility data for each molecular line and continuum using a modified version of the PHANGS-ALMA imaging pipeline (Leroy et al. 2021b).Here, we outline the workflow and highlight a few deviations from the original pipeline.
To prepare for continuum imaging, we extract all linefree channels from the calibrated measurement set and regrid them into a small number of channels per spectral window (SPW).Here, rather than collapsing each SPW into one channel (default behavior of the PHANGS pipeline), we manually choose a small enough output channel width to prevent bandwidth smearing.We then combine the continuum-only data from all array configurations (now on a common spectral grid) to make a joint dataset for the continuum in each band.
To prepare for line imaging, we model the continuum emission with a 1st-order polynomial based on all linefree channels across all SPWs of a particular band.We then subtract this model off from the calibrated mea-surement set, extract the subset of channels within the velocity range of interest (±300 km s −1 around the systemic velocity of 778 km s −1 ) for each emission line, and regrid the line spectrum to a desired channel width.Here, we choose a 10 km s −1 channel width for all high critical density molecular lines (HCN, HCO + , CS) to ensure a reasonable signal-to-noise (S/N) ratio per channel; we keep the native 0.85 km s −1 channel width for CO (3-2) as S/N ratio is not a concern there.We similarly combine the continuum-subtracted line data from all array configurations (now on a common spectral grid) to make a joint dataset for each emission line.
We image the joint, calibrated visibility dataset for each emission line and continuum following broadly the PHANGS imaging scheme (Leroy et al. 2021b).We first run a shallow, multiscale tclean to identify both compact and extended emission across the entire FoV down to S/N = 4.After this step, we run a deeper, singlescale tclean to pick up remaining emission down to S/N = 1 for the CO (3-2) line and S/N = 2 for all other lines and continua.To avoid divergence, we restrict this second step to within a cleaning mask, which is constructed based on archival, lower-resolution CO (2-1) data (Sun et al. 2020;Leroy et al. 2021a) and only includes ppv locations with significant CO (2-1) emission.In both tclean calls, we weigh the visibility data using the Briggs method with a robustness parameter of 1.0, which offers the best trade-off between sensitivity and resolution for our science goal (i.e., discerning individual YMCs).This robust weighting is applied to all continua and lines except for CO (3-2), which is bright enough to afford a uniform weighting that minimizes side-lobes and facilitates high dynamic range imaging.

Postprocessing and Data Characteristics
After obtaining the cleaned continuum images and line data cubes, we correct them for the primary beam response pattern and then mildly convolve them so that the beam shapes become round.We also convert all line data cubes to brightness temperature units (K) according to the corresponding line frequencies.
The output data cube for the CO (3-2) line reaches a spatial resolution of 0. ′′ 10 thanks to the uniform weighting.The noise level at the native 0.85 km s −1 channel width is 1.9-2.5 K, which already enables detection of CO (3-2) emission across a large number of channels at high significance in most regions.To further increase the image dynamic range, we also generate another version of the CO (3-2) cube by rebinning to a 5 km s −1 channel width, which lowers the noise down to 0.9-1.2K and allows us to pick up more diffuse emission (Table 2).

HST Data
We make use of multi-band Hubble Space Telescope (HST) images of NGC 3351 acquired by the LEGUS (Calzetti et al. 2015) and PHANGS-HST (Lee et al. 2022) surveys.These are broad-band imaging data obtained with the Wide Field Camera 3 (WFC3), together covering the full galaxy from UV to optical wavelengths (Figure 1 left panel; also see Figure 1 in Turner et al. 2021 for the sky footprint covered by each survey).The point spread functions (PSFs) of these data are ≲ 0. ′′ 1 in FWHM, slightly better than the beam sizes achieved with our ALMA observations.
To facilitate source cross-matching with ALMA, we take advantage of an updated PHANGS-HST star cluster catalog presented in Maschmann et al. (2024) and Thilker et al. (2024).This new catalog builds upon the previous PHANGS-HST catalog release (i.e., DR3/CR12 in Feburary 2022; see Lee et al. 2022).Aside from improved cluster classification schemes (Hannon et al. 2023, building on Whitmore et al. 2021;Thilker et al. 2022), a major improvement is in the spectral energy distribution (SED) fitting, for which a physically motivated "decision tree" has been introduced to help address photometric degeneracies between age, extinction, and metallicity (e.g., Turner et al. 2021;Whitmore et al. 2023a).Specifically, this decision tree uses source morphology, broadband color information, and location relative to Hα-bright complexes (based on ground-based narrowband Hα observations; A. Razza et al., in preparation) to pre-select candidates of young reddened clus-Figure 2. Top left: ALMA CO (3-2) moment-0 map of the central 1×1 kpc (∼20"×20") of NGC 3351, tracing the distribution of the bulk molecular gas.This moment-0 map is made from the CO (3-2) cube with 0. ′′ 10 beam size and 5 km s −1 channel width (see Section 2.1.2).The dotted ellipse in the plot corresponds to a galactocentric radius of 500 pc, which encircles the region of interest in this study.Top right: HST image for the same area overlaid with CO (3-2) integrated intensity contours in yellow (dashed -30 K km s −1 ; solid -400 K km s −1 ).Note the clear correspondence between the dust lanes and the CO contours, both of which trace gas distribution.Bottom: ALMA 93 GHz (band 3) and 350 GHz (band 7) continuum images overlaid with only the outer CO (3-2) contour (30 K km s −1 ) in light grey.The black circles mark the continuum sources detected in both bands, whereas the square/diamond symbols mark those detected only in band 3/7.

Sun et al.
ters.Their SEDs are subsequently fit with a separate set of stellar population models with young ages (≤6 Myr) and potentially higher extinction (up to E(B−V) = 3.0).This refined scheme results in more complete identification and better age estimates for young clusters with associated Hα emission throughout the PHANGS-HST sample.We refer interested readers to Thilker et al. (2024) for more thorough descriptions of the rationale, implementation, and results of this new SED fitting scheme.
In this work, we use this updated cluster catalog to cross-match with sources detected in the ALMA data and compare their estimated physical properties.Within our region of interest (i.e., the central region of NGC 3351), the new SED fitting scheme gives refined age estimates for ∼65% of clusters.This fraction is higher than typical due to the dusty and crowded nature of this region, which makes SED fitting with older methods particularly challenging.Further inspection suggests that the new age estimates are either similar to or more reasonable than the old solutions in ∼80−90% of the cases (see Section 6.1 in Thilker et al. 2024).The remainder are mostly caused by incorrect association of Hα emission to slightly older clusters (due to the much coarser ∼1.′′ 2 PSF of the ground-base Hα data used in the decision tree), which affect the accuracy of cluster age estimates primarily between 1−10 Myr.To mitigate this minor issue, we generally avoid using age estimates for individual clusters and only rely on more robust aggregate statistics, such as the total number of clusters younger than 10 Myr.We anticipate future works based on higher-resolution HST narrowband Hα observations (PIs: F. Belfiore, R. Chandar, D. Thilker) to thoroughly address this issue.

JWST Data
We make use of multi-band JWST NIRCam/MIRI images of NGC 3351 acquired as part of the PHANGS-JWST Cycle 1 Treasury survey (Lee et al. 2023, see Figure 1 right panel).These broad-and medium-band images are available from the PHANGS DR1.0.1 public image release (Williams et al. 2024).This work primarily uses data in the 2.0-11.3µm wavelength range, with PSFs ranging from 0. ′′ 066 (F200W with NIRCam) to 0. ′′ 38 (F1130W with MIRI) in FWHM.
For source cross-matching with ALMA, we also use a catalog of compact 3.35 µm sources constructed by J. Rodríguez et al. (in preparation) following the methodology of Rodríguez et al. (2023).In short, the sources are identified in the F335M image as local maxima above three times the estimated background level in the neighboring <100 pc area.This is done with a peak-finding algorithm find peaks in PHOTUTILS (Bradley et al. 2022).In more crowded regions, the software SExtractor (Source-Extractor; Bertin & Arnouts 1996) is also employed to deblend the sources with a Mexican-hat filter.This F335M-based source identifica-tion scheme allows for more systematic detections of embedded star clusters (see Rodríguez et al. 2023), which is especially important for cross-matching with ALMA sources in this work.
3. ANALYSES From the rich ALMA data set (see Table 2), we identify compact (sub-)millimeter continuum sources, characterize their observable properties, and cross-match them with sources in the HST and JWST data.We detail these analyses in the following subsections.

Identifying ALMA Sources
We identify bright, compact sources in the 93 GHz and 350 GHz continuum images as YMC candidates, following Leroy et al. (2018) and Emig et al. (2020).These sources are selected based on a minimum S/N = 5 at peak brightness.Considering the number of independent measurements (beams) across the region of interest (r < 500 pc, Figure 2), this detection threshold corresponds to a false positive rate of < 0.01, i.e., there is <1% chance for any of our selected sources to be spurious.We choose this rather stringent criterion to ensure sample purity, which is critical for following analysis based on the source number counts (Section 5.2).This is however at the expense of a higher false negative rate, i.e., missing many potentially real sources with lower peak brightness (as can be seen in Figure 2).
We identify 15 sources at 93 GHz and 11 at 350 GHz.Among these, 8 sources at 93 GHz coincide spatially with counterparts at 350 GHz.These are likely the same objects seen at both frequencies.Considering these as duplicates, we identify in total 18 ALMA sources, all of which are marked and labeled in Figure 2 (lower panels; also see Table 3).Almost all sources appear to be compact and well-isolated, with a single peak and a relatively round shape.The only exceptions are two closely connected sources in the 93 GHz image (sources #6 and #7 in Figure 2), which we separate by visual inspection.
All the 350 GHz sources and the majority of the 93 GHz sources are located in regions with bright CO (3-2) emission (see contours in Figure 2 bottom panels), consistent with forming YMCs surrounded by substantial molecular gas reservoir.There are however four 93 GHz-only sources (#1, #4, #13, #16) that are located in regions with no detectable CO (3-2) emission.We will discuss the nature of these sources in Section 5.1.

Characterizing ALMA Sources
We measure fluxes and sizes of the 18 sources in the 93 GHz and 350 GHz continua using the CASA task imfit.This task fits a 2D Gaussian function to the flux distribution within a user-defined region.First, we define a circular aperture around each identified source with a diameter of 0. ′′ 34, which is twice the beam 12.0 ± 3.9 18 160.989975,+11.701884 0.153 ± 0.028 5.9 ± 0.9 1.44 ± 0.35 5.3 Note-(1) source ID as appeared in Figure 2; (2) source RA and Dec coordinates; (3) continuum flux density at 93 GHz; (4) half-light radius at 93 GHz, defined as the geometric mean of the beam-deconvolved semi-major and semi-minor axes (for sources with unsuccessful beam deconvolution, we report the original values in parentheses); (5) continuum flux density at 350 GHz; (6) half-light radius at 350 GHz (note that all detected sources have successful beam deconvolution); ( 7) HCN (1-0) line integrated luminosity (sources without CO 3-2 emission appear as "• • • " due to the lack of prior information for determining velocity range; §3.2); (8) HCN (1-0) line velocity dispersion (i.e., effective width; §3.2).
FWHM3 .We then perform the imfit task on the continuum image within the circular aperture for each source.Note that we set the imfit parameter dooff=True in order to fit for and subtract a smooth background that may be present at the location of some sources.We also ask imfit to try deconvolving the beam size from the best-fit 2D Gaussian profile so that we can recover the intrinsic size of each source.We have visually inspected all imfit results to ensure reliability.
From the imfit results, we record the measured fluxes and radii (original or deconvolved) for all sources and list them in Table 3.The source fluxes range over S 93 = 0.03−0.15mJy at 93 GHz and S 350 = 0.6−1.9mJy at 350 GHz.The beam-deconvolved source half-light radii (calculated as the geometric mean of the semi-major and semi-minor axes) ranges over R hl, 93 = 1.4−6.0pc and R hl, 350 = 2.6−7.6 pc.We note that the deconvolution is unsuccessful for several 93 GHz sources because their observed major/minor axis lengths are smaller than the beam size, whereas none of the 350 GHz sources has similar issues.For the unsuccessful cases, we report the original 2D Gaussian sizes in Table 3 and treat these numbers as upper limits on the intrinsic source radii.
In addition to the continuum measurements, we also extract for each source a set of molecular line spectra (see Table 2 for a list of included lines).We define elliptical apertures for spectra extraction based on the best-fit 2D Gaussian for the continuum sources.For sources with both 93 GHz and 350 GHz detections, we use the apertures determined from the 350 GHz detections, as we expect molecular line emission to be more spatially associated with thermal dust emission than free-free emission.
We also perform background subtraction on the extracted line spectra.To do this, we define a "background  4).
annulus" around the source, calculate the mean emission spectrum within that annulus, and subtract this mean spectrum from the extracted spectrum at the location of the source.For isolated sources, the background annuli are set to have inner/outer radii of 2×/3× the aperture radii for source extraction.For clustered sources, we manually draw an annulus safely enclosing all the sources and set the outer radius to be 2× the inner radius.For source #2 that is in between two CO (3-2) blobs (Figure 2), we omit background subtraction since there is no obvious way to reliably determine the local background level in this case, and that level can be zero for some choices of the background aperture.
After extracting the emission line spectra and subtracting the local background, we measure the line integrated luminosities and line widths for all sources and all lines.We use the brightest line in our sample, CO (3-2), as a prior to determine the velocity range of each source and then measure line integrated flux and effective width (with the latter corrected for broadening due to finite channel width; see Heyer et al. 2001;Sun et al. 2018) for all high critical density lines.This work mainly uses measurements for the HCN (1-0) line, which we report in Table 3.The spectra of all molecular lines are presented in Appendix A.

Cross-matching ALMA-JWST-HST Sources
Note-Meaning of various markers used above: • -unique, unambiguous cross-match; ∞ -a close ALMA source pair (#6 and #7) cross-matched to the same HST cluster.
For the 18 sources detected in the (sub-)millimeter ALMA observations, we further check if they have counterparts in the UV-optical HST images and the nearto mid-IR JWST images.The presence or absence of each source at different wavelengths can inform us the presence or absence of various components in a forming cluster (e.g., gas versus stars) as well as the level of dust extinction towards that cluster.Such qualitative information is especially important for inferring the time evolution of YMCs, which we will discuss in depth below in Section 5.1.
For ALMA-HST source cross-matching, we take the star cluster catalog generated from the multi-band HST images (Section 2.2; also see Lee et al. 2022, and references therein) and focus on class 1 and 2 objects (i.e., compact clusters).As the HST images have higher resolution (sharper PSF) than the ALMA continuum images, we consider an HST cluster to be cospatial with an ALMA source if the distance between their central coordinates is smaller than the major axis half width at half maximum (HWHM) of the ALMA source.This crossmatching scheme yields six ALMA sources with HST cluster counterparts (Figure 3 and Table 4).These HST clusters have SED fitting-based mass estimates ranging Bottom right: MIRI F1000W-F1130W color image highlighting the PAH to continuum contrast.These two panels allow for identifying compact hot dust sources associated with H II regions and/or powered by clusters (Hassani et al. 2023).We find 8 ALMA sources with potential JWST counterparts based on either cross-matching to the 3.35 µm source catalog and/or visual inspection of the images (see Section 3.3 and Table 4).over 5×10 4 −4×10 5 M ⊙ and ages ranging over 1-4 Myr, though the ages for some clusters may be underestimated for the reasons described in Section 2.2.
For ALMA-JWST source cross-matching, we first rely on the 3.35 µm source catalog (see Section 2.3 and Rodríguez et al. 2023).Similar to the ALMA-HST source matching criterion, the center-to-center distance between a JWST 3.35 µm source and an ALMA source needs to be within the major axis HWHM of the ALMA source to be considered a match.For ALMA sources with multiple matched 3.35 µm sources, we use the one with a smaller separation.
One caveat with using only the 3.35 µm source catalog for cross-matching is that it misses a small number of sources that are only visible in other IR bands, partly due to different PSF sizes and sensitivity.To address this issue, we also visually inspect the multiband JWST images (Figure 4) near the location of each ALMA source to determine (1) if an IR counterpart is clearly visible in other bands but not included in the 3.35 µm source catalog, and (2) if a cross-matched 3.35 µm source is consistently present across other IR bands.We find the F200W data particularly helpful in this regard, thanks to its smaller PSF (FWHM 0. ′′ 066 versus 0. ′′ 11) and longer integration time (1200s versus 390s).Combining the catalog matching and visual inspection procedures, we find in total eight unique crossmatched sources (Table 4).

PHYSICAL PROPERTIES OF THE YOUNG MASSIVE CLUSTER CANDIDATES
From the measured YMC candidate sizes and flux densities in the ALMA (sub-)millimeter continua, their molecular line fluxes and line widths, as well as ancillary information available for their HST and JWST counterparts (if any), we estimate their gas masses (Section 4.1), stellar masses (Section 4.2), sizes (Section 4.3), and other key properties (Sections 4.4-4.8).We detail the derivations below and summarize key quantitative results in Table 5.

Gas Mass
We use the 350 GHz continuum emission as a dust tracer, from which we infer the masses and sizes of the YMC gas reservoirs.Generally speaking, the 350 GHz continuum includes not only thermal dust emission, but also free-free and potentially synchrotron emission.To estimate the fractional contribution of these components, we assume the 93 GHz and 350 GHz flux densities for each object are a mixture of thermal dust and freefree emission, with intrinsic spectral indices of α = 4.0 (dust) and α = −0.15(free-free), consistent with what previous studies adopted for similar systems (Emig et al. 2020, and references therein).With this simple, twocomponent decomposition, we derive the fractional contribution of dust and free-free emission at each frequency and find that thermal dust contribution indeed domi-Figure 5. Flux densities at 93 GHz versus that at 350 GHz for all identified continuum sources (the symbols are the same as in Figure 2).The dashed and dotted lines mark two series of expected relations assuming spectral indices of α = −0.15for free-free and α = 4.0 for dust.Almost all sources lie in the unshaded middle region, which suggests that their 93 GHz emission is dominated by free-free and their 350 GHz emission dominated by thermal dust contribution.
nates the 350 GHz continuum for all identified objects (see Figure 5).This conclusion agrees with previous YMC studies in other galaxies (e.g., Leroy et al. 2018;Emig et al. 2020), for which the free-free contribution to the 350 GHz continuum was also found to be < 10% in most cases.
The above analysis ignores the contribution of synchrotron emission, which may be non-negligible at 93 GHz (as found for some YMCs in other systems; see Emig et al. 2020;Mills et al. 2021).Qualitatively, we expect the inclusion of a synchrotron component to lower the estimated fraction of free-free at 93 GHz and raise the fraction of thermal dust at 350 GHz.On the one hand, this makes it even more reasonable to assume the 350 GHz continuum is dust-dominated for all our sources.On the other hand, it becomes an issue when interpreting the 93 GHz continuum, as discussed below in Section 4.2.
After verifying that the 350 GHz continuum is dominated by thermal dust emission, we convert the 350 GHz flux density into a gas mass in two steps.First, we calculate the dust optical depth at 350 GHz via . (1) Here Ω 350 is the solid angle extended by the 350 GHz source on the sky, and B ν (T dust ) the blackbody function for a given dust temperature T dust .We use a dust temperature of T dust = 130 K following Leroy et al. (2018), but there can be a ∼50% uncertainty on this value.We note that Equation 1 is a good approximation only in the optically thin limit.Nonetheless, our estimated optical depths of τ 350 ≈ 0.001-0.005for all sources make this approximation appropriate.
We then estimate the gas mass associated with each source via (2) Here κ 350 = 1.9 cm 2 /g is the adopted dust opacity at 350 GHz (Ossenkopf & Henning 1994, appropriate for dust mixed with gas at densities of ∼10 5−6 cm −3 ), and D/G = 1/100 the adopted dust-to-gas ratio (Draine et al. 2007).Both values are chosen for consistency with previous YMC studies (e.g., Leroy et al. 2018;Emig et al. 2020), but they can introduce systematic uncertainties on the order of ∼0.3 dex.A 350 = Ω 350 d 2 is the physical area of the 350 GHz source, which is needed to convert gas surface density into a total gas mass.The Ω 350 term effectively cancels out when combining Equations 1 and 2. We find gas masses of 3 × 10 4 −2 × 10 5 M ⊙ for the 11 YMCs detected at 350 GHz.For the non-detections, we put 5σ upper limits of ≈ 2×10 4 M ⊙ (see Table 5).These gas masses are consistent with expectations for YMC progenitors or the gas reservoir associated with forming YMCs (e.g., Portegies Zwart et al. 2010).Considering uncertainties on the adopted parameters in Equations 1 and 2, we expect an overall systematic uncertainty of ∼0.5 dex on these gas mass estimates.
As an independent sanity check, we also derive alternative gas masses from the HCN (1-0) line luminosities (Section 3.2 and Table 3) with a commonlyadopted HCN-to-H 2 conversion factor of α HCN ≈ 10 M ⊙ pc −2 (K km s −1 ) −1 (Gao & Solomon 2004), although lower values have been advocated for systems with more extreme conditions (García-Burillo et al. 2012;Shimajiri et al. 2017, c.f. Barnes et al. 2020) that in part resemble the local conditions in YMCs.Our HCN-based gas masses show clear correlation and broad agreement with the 350 GHz-based estimates (Figure 6).Although the former tend to yield slightly higher values, the observed discrepancy is attributable to the large systematic uncertainties on the adopted α HCN and several coefficients in Equations 1 and 2.

Stellar Mass
We use the 93 GHz continuum to constrain the freefree emission from ionized gas in H II regions, from which we infer the ionizing photon production rate and subsequently the mass of the stellar population associated with the YMC candidates.As mentioned in Section 4.1, the 93 GHz continuum may also include contribution from thermal dust and synchrotron emission.5), with the black and gray dotted lines marking the identity relation and a factor of three offset to either side.In addition to statistical uncertainties shown by the error bars, these mass estimates are also subject to ∼0.5 dex systematic uncertainties stemming from our assumptions of dust temperature, emissivity, and abundance as well as HCN conversion factor ( §4.1).The two gas mass estimates broadly agree within the range allowed by those uncertainties.
Our simple two-component modelling (dust+free-free; Figure 5) suggests that the dust contribution at 93 GHz is less than 20% for all but one source (#12).However, the issue of ignoring synchrotron emission in the modelling becomes more concerning at 93 GHz, as synchrotron may have a non-trivial contribution to the observed 93 GHz flux density (e.g., Emig et al. 2020  It is possible to better constrain the synchrotron fraction at 93 GHz given sensitive, lower frequency observations at similar angular resolution.For NGC 3351, previous observations at 1.4 GHz with the Multi-Element Radio Linked Interferometer Network (MERLIN) provide the best lever arm for this purpose (Hägele et al. 2010).With no point source detected in the central region at 0. ′′ 29 × 0. ′′ 17 resolution, they put a 6σ upper limit of 0.30 mJy for any compact source at 1.4 GHz.Adopting a characteristic spectral index of −0.75 for synchrotron emission, this translates to an upper limit of ∼0.013 mJy at 93 GHz.In this case, the synchrotron fraction is constrained to be ≲ 10% for the brightest sources at 93 GHz (S 93 ≈ 0.15 mJy) and ≲ 40% for the fainter ones (S 93 ≈ 0.03 mJy, Table 3).
Considering the lack of more sensitive low-frequency data to further constrain the synchrotron spectral index and fraction, we choose to use our observed 93 GHz flux densities directly as free-free flux densities, which may introduce a bias of at most 40% per our estimates.We convert the 93 GHz free-free flux densities into ionizing photon production rates, following Emig et al. (2020): (3) Here we use a characteristic electron temperature of T e = 6000 K that is typical of galaxy centers (e.g., Bendo et al. 2016;Emig et al. 2020;Henshaw et al. 2023).Equation 3 also implicitly assumes that the escape fraction of ionizing photons from the YMCs are negligible, which should be a reasonable assumption for the young, deeply embedded clusters studied in this work.
We further infer a stellar mass from the ionizing photon production rate via This conversion is based on a Starburst99 simulation (Leitherer et al. 1999) of a 3 Myr old stellar population with a Kroupa (2001) initial mass function and solar metallicity (Williams et al. 2022).If we had assumed a zero-age stellar population (see Leroy et al. 2018), the estimated masses would be 2× smaller; an age of 5 Myr would instead yield 5× larger values (see Emig et al. 2020).Our adopted 3 Myr is consistent with the 1-4 Myr age range given by UV-optical photometric SED fitting for the few sources with HST counterparts (Section 3.3, modulo caveats described in Section 2.2).Such young ages seem reasonable since most 93 GHz sources are deeply embedded and have substantial gas reservoirs around them (Section 4.1), which suggests that supernova explosions have not gone off (i.e., age ≲ 4 Myr).The adopted 3 Myr is also consistent with our inferred YMC evolutionary timeline below in Section 5.2.With Equations 3 and 4, we find stellar masses of 2 × 10 4 −1 × 10 5 M ⊙ for the 15 YMCs with detectable 93 GHz emission.For non-detections, the 5σ upper limit on the 93 GHz flux translates to ≈ 2.5 × 10 4 M ⊙ under the same assumptions (see Table 5).These stellar mass estimates are consistent with typical definition of YMCs (or super star clusters, SSCs; Portegies Zwart et al. 2010).The main source of systematic uncertainty is the adopted stellar age for Equation 4, which can introduce 0.3-0.6 dex of variations as discussed above.
For the few 93 GHz sources with HST counterparts, we find that the 93 GHz-based stellar mass estimates are systematically lower than the HST SED-based mass estimates (median offset ∼0.5 dex, Table 5).On the one hand, these sources are less embedded (visible in HST data), so it is possible that (1) they tend to be slightly older than 3 Myr on average, and (2) a nonnegligible fraction of the ionizing photons do escape from the clusters.Both of these would lead to underestimated stellar masses via Equation 3. On the other hand, the UV-optical SED-based mass estimates also have nonnegligible uncertainties as discussed in Section 2.2.We aim to address this issue in a follow-up study by modelling the full UV-optical-IR SED for all clusters identified in the joint HST+JWST dataset, potentially including the ALMA measurements whenever available.

Size
We use the estimated (sub-)millimeter source sizes as indicators of the YMC candidate physical sizes.At the distance of NGC 3351, the (beam-deconvolved) halflight radii correspond to 2-8 pc for the 11 sources detected at 350 GHz (Table 3).This range broadly overlaps the typical range of radii measured for YMCs in the Milky Way (R hl ∼ 0.5−5 pc; Portegies Zwart et al. 2010) and in other extragalactic systems (R hl ∼ 2−5 pc; Ryon et al. 2017).
For the YMCs that are only detected at 93 GHz, we instead use their estimated sizes from the 93 GHz image, which spans a similar range (1-6 pc).To be more precise, the 93 GHz free-free source sizes should reflect the sizes of the H II regions associated with the YMCs rather than the star clusters themselves.Nonetheless, for the eight sources detected in both bands, we find reasonable agreement between their 93 GHz and 350 GHz sizes (Table 3), which justifies this choice.
We note that the YMC size estimates may be affected by systematic effects due to finite data resolution (in addition to the statistical uncertainties quoted in Table 3).Our estimated source sizes are all within a factor of three of the beam size, as seen in similar studies for other galaxies over a range of data resolution (Figure 7; also see Leroy et al. 2018;Emig et al. 2020;Levy et al. 2021Levy et al. , 2022;;He et al. 2022).When examining our 350 GHz image at the 0. ′′ 13 ∼ 6 pc native resolution, we also find tentative evidence for substructures within some sources identified at the 0. ′′ 17 ∼ 8 pc working resolution.These observations suggest that the estimated source sizes may become smaller as one pushes to smaller beam sizes, which was indeed the case with studies of YMCs in NGC 253 (Leroy et al. 2018;Levy et al. 2021).That being said, the differences between our 350 GHz images at 0. ′′ 13 and 0. ′′ 17 are only marginal, so we choose to work with the latter for more straightforward flux and size comparisons with the 93 GHz data.We expect to test the resolution effect in a follow-up study with ALMA Cycle 9 observations (partially executed, PI: J. Sun), which will improve the data resolution by a factor of three.

Total Mass & Gas Fraction
Adding together the stellar and gas mass estimates (or upper limits) for each YMC candidate, we find total masses of M tot ≡ M gas + M ⋆ ≈ 0.3−3 × 10 5 M ⊙ (Table 5).The estimated gas fraction, f gas ≡ M gas /(M gas + M ⋆ ), ranges from ≲25% to ≳70%.This suggests that the YMC candidates identified from our ALMA data span a wide range of evolutionary stages, from the gas mass-dominated early phase to the stellar massdominated late phase of cluster formation.
It is worth emphasizing that the total mass and gas fraction estimates can be affected by the ∼0.5 dex systematic uncertainties on M gas and M ⋆ mentioned in Section 4.1 & 4.2.For example, because M ⋆ as calculated from Equation 4 can vary by a factor of 2-5 depending on the assumed stellar age, the f gas estimates should be similarly sensitive to this assumption.Indeed, Leroy et al. (2018) and Mills et al. (2021) assumed zero-age stellar population for the YMCs in NGC 253 and found comparable stellar and gas masses (modulo synchrotron contamination); in contrast, Emig et al. (2020) assumed a 5 Myr age for the YMCs in NGC 4945 and found them to be stellar mass-dominated.While the different ages assumed in these works are well motivated by known differences of the host galaxies, it is still possible that these different assumptions may partly cause the differing results between studies.This underlines the importance of accurate YMC stellar age estimates, which we expect to improve with the joint HST+JWST SED fitting.

Virial Mass
We use the measured radius (R hl ; see Section 4.3) and gas velocity dispersion (HCN 1-0 line width, σ HCN ; see Section 3.2) to derive the virial mass for each object: Here the geometrical factor f depends on the density profile of the YMCs (Bertoldi & McKee 1992).We adopt a density profile of ρ ∝ r −2 following Leroy et al. (2018), which corresponds to f = 5/3.
We find virial masses of 10 5 −10 6 M ⊙ for the 12 YMC candidates with HCN (1-0) line measurements available (Table 5).Most of them have M vir ≳ 2 M tot , which seems to imply their gas reservoirs are not gravitationally bound or only marginally bound (Figure 8).However, the estimated virial masses have large statistical and systematic uncertainties due to the uncertain size and velocity dispersion measurements, and they also tend to be overestimated because (1) the actual YMC sizes may be smaller than those measured at our data resolution (Section 4.3), and (2) the HCN velocity dispersion may also be biased high due to contamination by diffuse emission from the lower-density ambient gas, despite our best attempt in removing them (Section 3.2).With these caveats in mind, we conclude that the measured sizes and gas velocity dispersion for the YMC candidates are consistent with them having high total masses of ≳ 10 5 M ⊙ .

Escape Velocity
We estimate the escape velocity for all YMC candidates from their total masses and radii: The factor of 3.4 in the numerator accounts for the fraction of mass inside the FWHM of a 3D Gaussian (following Leroy et al. 2018), which is appropriate since we are using R hl (i.e., the deconvolved source HWHM; see Section 3.2) in the denominator.We find escape velocities of 6−10 km s −1 among all YMC candidates (Table 5).These estimates are uncertain by ∼0.3 dex due to systematic errors on the mass estimates; they may also be biased low if R hl is overestimated at the current data resolution.Considering these factors, the true escape velocities are likely comparable to, if not larger than, the characteristic sound speed of photoionized gas (∼10 km s −1 ).In this case, we expect photoionization feedback from massive stars to be not as effective in removing gas and stopping star formation in these dense structures as they would be in less dense structures (consistent with a defining criterion for young massive protoclusters; see Bressert et al. 2012).In addition, the estimated v esc almost certainly exceed the previously estimated threshold of ∼1 km s −1 for effective gas ejection by protostellar outflows (Matzner & Jumper 2015), in line with the notion that this particular feedback channel has limited impact on cluster formation (e.g., Nakamura & Li 2014).These arguments imply relatively high star formation efficiency during the birth of these YMCs and the need for other feedback mechanisms to remove the remaining gas (e.g., radiation pressure and stellar winds; see Levy et al. 2021;Menon et al. 2023;Polak et al. 2023, also see Section 4.8).

Volume Density & Free-fall Time
We also derive the volume density and gravitational free-fall time from the YMC total masses and radii: We similarly include a factor of 3.4 in Equation 7as in Equation 6.That is, we are calculating the mean mass volume density within a sphere defined by the (deconvolved) FWHM of each source, and then the free-fall time corresponding to that volume density.
We find volume densities of 50−200 M ⊙ pc −3 , which is equivalent to hydrogen nuclei number densities of 1.5−6 × 10 3 cm −3 .The corresponding free-fall times are 0.5−1 Myr (Table 5), with ∼0.3 dex of systematic uncertainty coming from the mass estimates and a potential bias towards larger values due to marginally resolved sizes (see Section 4.3).Given that the YMC formation process can last a couple of free-fall times (Skinner & Ostriker 2015), it is possible that these YMCs are still forming and possess a substantial gas reservoir at ages of ∼3 Myr (also see related discussions in Sections 4.6-4.8).Moreover, the short free-fall time relative to the typical supernova explosion delay time (≳3 Myr) means that the initial gravitational collapse may be too fast for supernova feedback to play a role (Fall et al. 2010).This inference again supports fairly high star formation efficiencies when forming these dense YMCs.

Gas Surface Density
We calculate the gas surface density near the center of each YMC candidate via: Here A eff = 2πσ 2 xy = πR 2 hl / ln 2 is the effective area of the 2D Gaussian fit for each detected YMC in 350 GHz (Section 3.2).The second equality follows from the conversion between the Gaussian HWHM (i.e., R hl ) and dispersion: σ xy = R hl / √ 2 ln 2. We find gas surface densities of 500−2000 M ⊙ pc −2 for the YMC candidates with 350 GHz continuum detections (Table 5), with a ∼0.5 dex systematic uncertainty associated with M gas and a potential bias towards lower values due to possibly overestimated sizes at current data resolution.The estimated as surface densities are much higher than typical values of giant molecular clouds (∼10 2 M ⊙ pc −2 ; Heyer & Dame 2015; Rosolowsky et al. 2021;Sun et al. 2022) and also higher than the median value for dense clumps in the Milky Way (∼500 M ⊙ pc −2 ; Urquhart et al. 2018).However, they are close to surface density thresholds above which one expects high star formation efficiency, as shown by analytical and numerical studies (≳ 10 3 M ⊙ pc −2 ; see e.g., Fall et al. 2010;Kim et al. 2018;Menon et al. 2023).
The gas surface density can play a key role in regulating the star formation efficiency because it determines the effectiveness of various forms of feedback in destroying the natal molecular clouds (e.g., see section 5 of Chevance et al. 2023).In particular, radiation pressure is believed to be the most important feedback process at high surface densities.For direct (UV) radiation pressure, one can derive an Eddington ratio, or the ratio of radiation pressure force to gravitational force, from the gas surface density and gas mass fraction  7) gas mass fraction (or upper/lower limit for non-detection at 350/93 GHz; §4.4); (8) total mass (or 5σ range for non-detection in either band; §4.4); (9) virial mass (for those with measured gas velocity dispersion, §4.5); (10) escape velocity (or 5σ range / lower limit according to columns 2 & 8; §4.6); (11) free-fall time (or 5σ range / upper limit according to columns 2 & 8; §4.7); (12) gas surface density (or upper limit based on 350 GHz continuum sensitivity, §4.8).We note that most quantities reported here have large systematic uncertainties.The YMC radii are likely affected by finite data resolution ( §4.3); the gas and stellar mass estimates have ∼0.5 dex systematic error from the assumptions involved in their derivations ( §4.1-4.2);other columns are also affected by the propagation of these uncertainties ( §4.4-4.8).* These four sources are detected only in 93 GHz continuum.They may be supernova remnants with strong synchrotron emission rather than YMCs (see §5.1), in which case the M⋆, 93 values and subsequent calculations would not be reliable.
where Ψ ≈ 10 3 L ⊙ /M ⊙ is the light-to-mass ratio of a central stellar population with a Kroupa IMF and an age of ≲ 3 Myr (based on a Starburst99 simulation).The above is derived for a uniform gas sphere with a central star cluster, but other spatial distributions give similar results (e.g., Raskutti et al. 2017;Reissl et al. 2018;Krumholz et al. 2019).From this equation, we estimate f Edd ∼ 0.2−0.7 for the YMC candidates with associated gas, which implies that their direct radiation pressure may not (yet) be enough to expel the gas (assuming the gas column density variation from sightline to sightline is small; see Thompson & Krumholz 2016;Raskutti et al. 2016Raskutti et al. , 2017)).Furthermore, at Σ gas ≳ 10 3 M ⊙ pc −2 , the star formation efficiency is expected to exceed 50% when limited only by direct FUV radiation pressure on dust in combination with the rocket effect and photoevaporation from EUV (Kim et al. 2018;He et al. 2019).Stellar winds, since they have a similar momentum flux to that from radiation, would not significantly change the expected star formation efficiency (Lancaster et al. 2021;Polak et al. 2023).Unless the mass function is top-heavy or the dust abundance is enhanced, reprocessed infrared radiation has f Edd = κ IR Ψ/(4πGc) < 1 and therefore does not aid in limiting star formation either (Skinner & Ostriker 2015;Menon et al. 2022).To conclude, the Sun et al.
high gas surface densities found for many of our YMC candidates suggest that the remaining gas may continue to feed star formation, leading to an overall high star formation efficiency.
Other than the implications on feedback and star formation efficiency, the measured gas surface densities also imply V -band extinctions of ∼30-110 mag (assuming Galactic dust abundance and extinction curve; Draine 2011), or ∼15-55 mag if half of the gas is in front of the stellar body.Such high extinction helps explain the lack of optical (or even IR) counterparts for many of our ALMA sources, as was found for YMCs in other systems (e.g., Leroy et al. 2018;He et al. 2022).

DISCUSSION
In Section 4, we report the measured physical properties for the 18 YMC candidates identified from our new ALMA data and complemented by JWST and HST data.Our measurements suggest that most of these objects are massive, compact, in the early phase of formation, and likely associated with high star formation efficiency.Here we build on the quantitative results and address a few key remaining questions: (1) What do the rich multiwavelength measurements tell us about the nature and evolutionary stages of the YMC candidates?
(2) How do the ALMA-identified YMCs relate to the exposed, more evolved clusters?(3) How do the YMCs fit in the large-scale context of the central starburst ring in NGC 3351?

Nature and Evolutionary Stages of the ALMA-identified YMC Candidates
Based on the multiwavelength observational properties of the 18 ALMA sources and their cross-matching results with the JWST and HST data, we can classify most of them into four categories.These categories are inspired by the classification schemes for star clusters introduced in Johnson (2005) and Whitmore et al. (2014) as well as a similar scheme for molecular clouds used in Kawamura et al. (2009).
• Type 1: Starless Clumps (#9, #10, #11).These objects have substantial gas reservoirs indicated by thermal dust emission at 350 GHz, yet they show no signs of star formation through UV-optical stellar photospheric emission or 93 GHz free-free emission from H II regions.They are likely dense gas clumps on their way to becoming YMCs, representing the earliest phase of YMC formation.
• Type 2: Clump-H II region complexes (#3, #5, #12): These objects have both gas reservoirs producing dust emission at 350 GHz and associated H II regions producing free-free emission at 93 GHz, but the stellar content remains largely invisible in UV-to-IR bands due to high extinction.They represent the deeply embedded phase of YMC formation, when a substantial stellar body (including massive stars) has formed but the gas reservoir is neither expelled nor exhausted.
• Type 3: Clump-H II region-cluster complexes (#6, #7, #14, #15, #18): These objects are detected simultaneously in 350 GHz dust emission, 93 GHz freefree emission, as well as PAH and stellar photospheric emission in the near-to mid-IR; three out of five are also visible in optical bands.These are exposed but still forming clusters as they still have ample gas left and H II regions glowing around massive stars.
• Type 4: Exposed H II region-cluster complexes (#2, #8, #17): These objects no longer have detectable 350 GHz emission, suggesting that much of the local gas reservoir has been expelled or exhausted.They are still visible in both 93 GHz free-free emission and UV-to-IR stellar photospheric emission.These are emerging young clusters that have stopped forming and likely survived the violent gas expulsion process near the end of their formation process.
There are four ALMA sources (#1, #4, #13, and #16) that do not belong in any of the above categories.They are visible at 93 GHz but almost completely missing in all other wavelengths.Considering the stringent upper limits we can put on the gas surface density based on non-detections of the 350 GHz continuum, the HCN (1-0) line, and even the CO (3-2) line (see Appendix A), these objects cannot be deeply embedded and should be visible at least in the near-IR if they are indeed forming YMCs with H II regions glowing intensely in free-free emission.
One possible explanation is that these 93 GHz-only sources are not forming YMCs, but rather supernovae remnants (SNRs) with their 93 GHz emission dominated by synchrotron radiation.As mentioned in Section 4.2, previous 1.4 GHz MERLIN observations (Hägele et al. 2010) put an upper limit of 0.3 mJy for any source more compact than the ∼10 pc MERLIN beam.Combining this constraint at 1.4 GHz with the measured 93 GHz flux density of 0.03−0.04mJy for the four ALMA sources in qustion, the implied radio spectral index is ≳ −0.5.This is consistent with a synchrotron spectrum that turns over at an intermediate frequency due to synchrotron self-absorption, which is very probable given the small sizes and strong magnetic field of SNRs (e.g., see Lenc & Tingay 2009 for spectra of SNRs in NGC 4945 with turnover frequencies of ∼5 GHz).Furthermore, Hägele et al. (2010) showed that a SNR at their detection threshold of 0.3 mJy is expected to have a diameter of ∼2 pc.SNRs slightly fainter and larger in size would remain undetected in the MERLIN observations, while they can be detected but unresolved in our ALMA observations, just like the four 93 GHz-only sources.We thus suggest that these sources may actually be SNRs with synchrontron emission dominating the 93 GHz continuum.Future observations with the Karl G. Jansky Very Large Array (VLA) at 3-30 GHz with matched resolution can help verify this hypothesis.

YMCs versus Evolved Clusters
In the previous section, we classified YMCs into four categories, which likely correspond to four evolutionary stages in the cluster formation process.However, there is a fifth stage, i.e., evolved clusters without associated gas or H II regions, that has been omitted from the discussion above.This omission is obviously due to our target selection based on ALMA data -older clusters without 93 GHz free-free or 350 GHz dust emission should not be in our list of ALMA-identified sources to begin with.These older clusters can nonetheless be identified with the HST and JWST data; they are also of great interest as they represent the end products of the cluster formation process.Here, we explicitly connect the ALMA-identified YMCs to more evolved clusters included in the HST cluster catalog (Section 2.2) to fill in the last missing piece.
One important consideration for such cross-dataset comparison is the difference in data sensitivity, which maps into different mass sensitivity for cluster detection.The ALMA continuum observations were designed to probe forming YMCs down to ∼ 10 4.5 M ⊙ in both gas and stellar mass.The HST broad-band observations, in contrast, can detect exposed clusters with much lower masses (∼ 10 3.5 M ⊙ ; Lee et al. 2022).In order to compare the cluster populations detectable by these datasets in a sensible way, we choose to focus on HST clusters above some matching mass threshold.Given the systematic uncertainties in the mass estimates from both sides, we use two thresholds to bracket a sensible range: one at 10 4.5 M ⊙ to directly match the nominal sensitivity limit of the ALMA data, the other at 10 5 M ⊙ to account for the median ∼0.5 dex offsets between the 93 GHz-based and UV-optical SED-based stellar mass estimates for the ALMA-HST cross-matched sources (see Table 5).
The HST cluster catalog includes in total 42 star clusters in the central region of NGC 3351 (r gal < 500 pc).The vast majority (37) of them have SED-based ages < 10 Myr and are located between r gal = 200−400 pc (Figure 9 right panel).The same r gal range covers all 14 ALMA-identified YMCs (including the six ALMA-HST cross-matches and excluding the four sources that are possibly SNRs) and almost all gas structures associated with the starburst ring.Since the orbital period along this ring is also ∼10 Myr (Rubin et al. 1975;Hägele et al. 2007;Leaman et al. 2019), we do not expect much radial migration for these <10 Myr clusters.We thus view these young HST clusters along the ring as either counterparts or direct descendants of the ALMA YMCs (modulo caveats with mass range matching).
From the number counts of ALMA YMCs and young HST clusters within the matched mass range, we can further infer the timescales of the various evolutionary stages during YMC formation (similar to the method used for determining molecular cloud evolution timescales by Kawamura et al. 2009).Here we assume that (1) the SFR of the starburst ring does not change drastically over a 10 Myr timescale; (2) the combined ALMA+HST YMC sample is complete over the matched mass range up to an age of ∼10 Myr; and (3) there is no substantial loss of YMCs over this period.Under these assumptions, the number count of YMCs in each evolutionary stage should be proportional to the average duration of that stage.Furthermore, the combined sample of ALMA YMCs with observable signs of active star formation (i.e., Types 2-4 with free-free emission) plus HST clusters with age ≤10 Myr should span a full 10 Myr window from the onset of star formation.This latter inference provides a necessary reference timescale for converting the relative duration of all stages to absolute values in Myr unit.
We show our inferred evolutionary timeline of YMC formation in Figure 10.Depending on the mass range of HST clusters (≥10 4.5 M ⊙ or ≥10 5 M ⊙ ) used as the reference sample, we find that the typical duration of the four stages (mapping to Types 1-4) are 1-1.7 Myr, 1-1.7 Myr, 1.8-2.8Myr, and 1-1.7 Myr, respectively.That is, a starless clump lasts ∼1-2 Myr (or ∼2 t ff ) before forming the first massive stars capable of creating H II regions.From that point on, it takes ∼1-2 Myr for the star cluster to become visible in the near-IR and ∼2-3 Myr to become visible in the optical.The cold gas reservoir disappears over ∼3-4 Myr (or 4-6 t ff ), and the associated H II region fades away in ∼4-6 Myr.Overall, these timescales agree well with previous observations and simulations of various types of star-forming regions giving birth to star clusters (e.g., Whitmore & Zhang 2002;Tan et al. 2006;Reggiani et al. 2011;Whitmore et al. 2014Whitmore et al. , 2023b;;Skinner & Ostriker 2015;Grasha et al. 2018Grasha et al. , 2019;;Kim et al. 2018Kim et al. , 2019;;Hannon et al. 2019Hannon et al. , 2022;;Kruijssen et al. 2019;Li et al. 2019;Chevance et al. 2020;Grudić et al. 2021;Kim et al. 2021Kim et al. , 2023)).
We would like the emphasize that there are substantial uncertainties on the inferred timescales, such that they should be viewed as preliminary, order-of-magnitude estimates.For example, it is possible that some of the assumptions mentioned above are not appropriate for the starburst ring in NGC 3351.The 10 Myr orbital period of this system means that the SFR can fluctuate moderately on a similar timescale (see observational constraints by Calzetti et al. 2021; also see demonstrations of such behavior in simulations by Armillotta et al. 2019;Sormani et al. 2020;Moon et al. 2022).Star clusters may also get destroyed within 10 Myr due to either violent gas removal at the end of their formation process or mass loss due to stellar evolution (see Krumholz et al. 2019, for a thorough review).These concerns could be addressed, in theory, by choosing even younger HST clusters (e.g., <5 Myr) as the reference sample, but various sources of systematic effects on the cluster age estimates especially at ≲10 Myr (Section 2.2) would render such Sun et al. analysis unreliable in practice.On this front, follow-up studies probing more diverse environments (e.g., those with longer dynamical timescales) will be crucial to verify our results in different physical conditions and to improve the currently limited statistics.
The reliability of the inferred timeline is also fundamentally tied to the reliability of the HST cluster measurements themselves.As discussed in Section 2.2, the revised SED fitting scheme presented in Thilker et al. (2024) reflects the latest and by-far the most systematic efforts in dealing with photometric degeneracies, but we expect follow-up studies to further refine the fitting results and address remaining issues (especially with the youngest clusters).For example, these goals can be achieved by employing HST Hα narrowband and JWST IR data to complement the HST broadband photometry (e.g., K. Henny et al., in preparation).Additional narrowband Paα imaging with JWST in Cycle 3 (PI: A Leroy) and/or radio continuum observations with ALMA and VLA in the future would also provide unique constraints to help pin down the ages.

YMCs in Large-scale Context
After examining the plausible evolutionary phases and timescales of the ALMA YMCs, we now put their properties in the context of the large-scale environment, i.e., the central starburst ring in NGC 3351.
We can use the estimated masses and timescales for the ALMA-identified YMCs to estimate the fractional SFR contributed by these sources alone.The total stellar mass of all YMCs (excluding the four sources that are possibly SNRs) is 7 × 10 5 M ⊙ , and the inferred duration of the 93 GHz-bright phase is 4-6 Myr (Section 5.2).Dividing these two numbers gives us an instantaneous SFR of 0.1-0.2M ⊙ yr −1 , which is already ∼50% of the total SFR of the starburst ring.Such a high fraction is consistent with similar estimates for the central starburst regions in NGC 253 and NGC 4945 (Leroy et al. 2018;Emig et al. 2020).Note that this simple calculation ignores the considerable gas mass associated with the YMCs, and a fraction of that gas may be converted into stars in the future.Nor does it include the contribution from less massive clusters below our detection threshold.Given these omissions, the actual fraction of stars formed in clusters may be even higher.• Testing ring star formation models with the YMCs: Another noteworthy feature of NGC3351's starburst ring is its clean orbital configuration at a favorable viewing angle (see Figure 11 left panel), which stands out among all systems whose YMC populations have been studied by ALMA so far.As an iconic bar-fed nuclear ring with unambiguous orbital streamlines, this system offers an opportunity to test existing models of star formation in similar systems.
There are at least two influential theoretical models of star-forming rings in the literature (see Böker et al. 2008, for a nice summary).A first scenario, dubbed the "popcorn model", considers star formation triggered by gravitational instability as gas accumulates on the ring and reaches a critical density (Elmegreen 1994).As star formation happens stochastically across the system, there may not be any preferable locations along the ring or systematic azimuthal trends in the properties of the young star clusters.
A second scenario, namely the "pearls on a string model" (Böker et al. 2008), suggests that star forma-tion is mostly triggered near or slightly downstream from the contact points, where the gas inflow enters the ring.An alternative model proposed for the Galactic Central Molecular Zone (CMZ) argues that star formation should be triggered by gas tidal compression near the pericenters of the ring orbit (e.g., Longmore et al. 2013, also see Callanan et al. 2021 for similar analyses on M83's center).In either case, an evolutionary sequence is expected downstream from the triggering positions, with progressively older clusters further down the orbit.
As the YMCs identified in our ALMA observations probably have ages younger than the orbital period (see Section 5.2), their distribution and azimuthal trends may be particularly helpful for differentiating the competing scenarios described above.We find that the ALMA YMCs appear to concentrate near the two contact points (Figure 11   tative evidence for a progression in YMC evolutionary stages from ϕ = −180 • to 0 • (i.e., the Eastern side of the ring), with Type 1 sources located right at the Northern contact point and subsequent types further down the ring orbit; a similar trend could be argued for the YMC gas to stellar mass ratio.Nonetheless, these trends become less clear on the other (Western) side of the ring.There is no clear trend or preferential distribution in the HST clusters either.
A subset of our observational results, especially the concentration of YMCs near the contact points and the possible azimuthal trends, seems to favor the "pearls on a string" model.Notably, the apparent progression from Type 1 to Type 4 YMCs spans half of the ring orbit from the Northern to the Southern contact point, and our estimated 4-5 Myr timescale for such progression (Section 5.2) does agree with the ∼5 Myr orbital time across half of the ring (Leaman et al. 2019).Nonetheless, these trends are only based on a small number (∼8) of ALMA YMCs and do not show up consistently across the entire ring or in the (more evolved) HST clusters.
We can partly make sense of the somewhat mixed results in light of recent numerical studies of star-forming rings.For example, Seo & Kim (2013) showed that the presence or absence of a cluster age sequence may depend on the total SFR of the ring.Sormani et al. (2020) found that the instantaneous distribution of clusters can vary substantially with time due to stochasticity of star formation along the ring, and that one may only see clear azimuthal age trends either through time-averaging or when focusing on the youngest clusters (≤0.25 Myr in their system with ∼5 Myr orbital period).Moon et al. (2021Moon et al. ( , 2022) ) showed that local concentration of young clusters upstream from contact points can happen following a temporary, asymmetric boost of gas inflow rate.Such asymmetry and non-steadiness can naturally arise from clumpiness in bar-driven inflows as well as varying accretion efficiency onto the ring (Sormani & Barnes 2019;Hatchfield et al. 2021).Together, these studies highlighted key aspects of star-forming rings that are not entirely captured by the simple "popcorn" or the "pearls on a string" models, but can help explain some of the observed behaviors in NGC 3351's central ring and its YMC population.
We anticipate more sensitive observations in the future (including partially executed ALMA Cycle 9 observations) to improve the statistics especially for the youngest clusters and allow for more detailed comparisons with simulations.Studying a large sample of starforming ring systems (with favorable viewing angles) will be critical for differentiating time-varying effects from persistent trends.
• YMCs as potential drivers of large-scale gas outflows: Central starburst rings are believed to be powerful drivers of multi-phase galactic winds and outflows (e.g., Armillotta et al. 2019;Nguyen & Thompson 2022), and there are indications of such outflows from NGC 3351's central ring from existing multiwavelength observations.By analyzing Chandra X-ray imaging data, Swartz et al. (2006) found evidence of hot gas expanding beyond the starburst ring and likely above the plane of the galaxy.This hot gas is estimated to contain thermal energy of ∼f 1/2 10 54 erg (here f ≳ 10 −4 is the volume filling factor of the gas).A more recent study based on VLT/MUSE data (Leaman et al. 2019) identified and modelled a stream of warm ionized gas outflow, with a radial velocity of ∼70 km s −1 and kinetic energy of ∼6 × 10 52 erg.Both studies highlighted a dust lane visible in optical images to the southeast of the ring (see Figure 1 left panel, right outside the white box) and interpreted it as a gas shell confining the outflow.The new JWST/MIRI images (Figure 1 right panel) clearly show that a similarly shaped gas shell exists on the opposite side, even though it is less obvious from the optical image.Therefore, gas outflow may be present on all sides of the ring and contain more kinetic energy than estimated before.
The YMCs studied in this work are likely too young to be the main energy source for driving these hot and warm gas outflows.With our inferred timeline of ≲4 Myr after birth, most of them have not produced many supernovae to heat the hot X-ray emitting gas; their arrival was also too late to accelerate the aforementioned gas shell to its current location (which was estimated to take ≳10 Myr by Leaman et al. 2019).Besides, we do not see clear evidence of localized outflow around the YMCs in their molecular line spectra either (see Appendix A), unlike those found for some of the YMCs in NGC 253 (Levy et al. 2021, modulo the different spatial resolutions of the observations).
Nonetheless, it is still interesting to compare the winddriving capability of the YMCs to the multi-phase gas outflow seen at the moment.Our Starburst99 simulation suggests that the YMCs can together produce a mechanical luminosity of ∼2 × 10 40 erg s −1 , or a total deposited mechanical energy of ∼3 × 10 54 erg over 4 Myr.A mere 2% energy retention factor (easily reachable for stellar wind and supernova-driven outflows; see Kim et al. 2020;Sirressi et al. 2024) would be enough to power the ionized gas outflow reported by Leaman et al. (2019) and to heat the X-ray emitting gas studied by Swartz et al. (2006).These order-of-magnitude calculations suggest that the current population of YMCs will be more than capable of drive gas outflow at a similar level in the future.Therefore, the multi-phase gas outflow may get enhanced over the next ∼10 Myr and become closer to those seen in NGC 253 and NGC 4945 (e.g., Westmoquette et al. 2011;Krieger et al. 2019;Bolatto et al. 2021).

CONCLUSIONS
In this paper, we examine a population of embedded YMCs in the central 1×1 kpc region of the nearby galaxy NGC 3351.This system features a prominent central starburst ring fed by stellar bar-driven gas inflows from the outer disk (Regan et al. 2006;Leaman et al. 2019).The proximity (9.96 Mpc; Anand et al. 2021), favorable viewing angle (i ≈ 45 • ; Lang et al. 2020), and clean orbital configuration of this system (Figure 1 and 2) make it ideal for a multiwavelength YMC study in the full context of the large-scale host galaxy properties.
To this end, we acquire new ALMA data in Band 3 and 7 (project code 2021.1.00059.S; PI: J. Sun), targeting the 93 GHz and 350 GHz continua along with various molecular lines.The long-baseline observations ensure that the ALMA images reach similarly high resolution (0. ′′ 1−0.′′ 2, equivalent to 5−10 pc) as existing HST and JWST images (Lee et al. 2022(Lee et al. , 2023)), while the shorterbaseline observations guarantee robust imaging by recovering emission on larger spatial scales.The jointlyimaged ALMA data are sensitive enough to detect thermal dust continuum and molecular line emission from the gas reservoir of individual forming YMCs, as well as free-free emission from H II regions created by the YMCs.By cross-matching ALMA sources with those in JWST and HST data, we probe the (often embedded) YMC stellar population, thereby achieving a complete, multiwavelength view of YMC formation.
Our joint analyses of the ALMA, HST, and JWST datasets yield the following key results: 1. We find 18 bright, compact sources in the ALMA continuum images, with 15 detected at 93 GHz and 11 detected at 350 GHz (Figure 2).Subsequent source cross-matching shows that only 8 of them have potential counterparts in JWST images (Figure 4) and 6 have counterparts in HST images (Figure 3; also see Table 4).
Sun et al.
3. The estimated size, mass, and velocity dispersion also imply high total mass (≳10 5 M ⊙ ), a wide range of gas fractions (from ≲25% to ≳70%), large escape velocity (6−10 km s −1 ), short free-fall time (0.5−1 Myr), and high gas surface density (500−2000 M ⊙ pc −2 ), as summarized in Table 5.The last three quantities suggest that various forms of feedback (photoionization, supernova, and radiation pressure) may be less effective in regulating star formation for sources with such extreme conditions (see Sections 4.6-4.8).
4. The multiwavelength properties of these ALMAidentified sources motivate a classification scheme in which most of them belong to one of the following categories: starless clumps (N = 3), clump-H II region complexes (3), clump-H II region-cluster complexes (5), and exposed H II region-cluster complexes (3).These four categories likely represent four phases of YMC formation.The remaining 4 sources (detected only by ALMA at 93 GHz) are possibly supernova remnants with strong synchrotron emission rather than forming YMCs, though follow-up observations are necessary to verify this interpretation.
5. Comparing the number counts of ALMA-identified YMCs versus young HST clusters in a matched mass range, we infer an evolutionary timeline for forming YMCs (Figure 10).Modulo various sources of uncertainty, we estimate a duration of ∼1−2 Myr (or ∼2 free-fall times) for the starless clump phase.It then takes ∼1−2 Myr and ∼2−3 Myr for the newly formed cluster to become visible in the IR and optical bands, respectively.The cold gas reservoir disappears over ∼3−4 Myr (or 4−6 free-fall times), and the H II region disappears over ∼4−6 Myr.These numbers represent our best constraints on YMC formation timeline in an extragalactic, starburst-like environment.They also agree quantitatively with previous estimates by observational and numerical studies of various types of cluster-forming environments.
6. Putting the YMCs in the context of the entire central starburst region of NGC 3351, we find that the YMCs alone can account for at least 30-50% of the total SFR of the ring (0.3−0.5 M ⊙ yr −1 ).While the YMCs exhibit an uneven azimuthal distribution and concentrate towards the two "contact points," there is no consistent azimuthal trend in the inferred YMC evolutionary stages or other properties, as one may naively expect if YMC formation is only triggered by colliding flows near the "contact points."Last but not least, the estimated total mechanical luminosity of the YMCs is large enough to power the previously reported multi-phase gas outflow from this system.
The quantitative measurements presented in this study will likely be improved in follow-up studies based on existing and/or future observations.For example, we plan to revisit the cluster selection and SED fit-ting by jointly analyzing the HST and JWST broadand medium-band data shown in this work as well as archival HST narrow-band Hα data (PI: R. Chandar; also see Calzetti et al. 2021).The inclusion of deeper 93 GHz continuum observations with ALMA (partly executed in Cycle 9), lower frequency observations with VLA, or narrow-band imaging targeting IR recombination lines with JWST will certainly provide much better constraints on stellar age and extinction, eliminating a major source of uncertainty.
Beyond the YMCs, we also expect our rich ALMA dataset to support many other science goals.The deep, highly resolved CO (3-2) image makes it possible to characterize the complex, multi-scale gas structures present throughout NGC 3351's central region.Detailed analysis of the gas kinematics may also shed light on the baryonic cycle across this system, from the large-scale gas inflow, to gas depletion and expulsion due to star formation and feedback, to the feeding (or lack thereof) of the central supermassive black hole.Extracted molecular line spectra for two ALMA sources (#1 and #3).The former represents a non-detection across all lines while the latter is detected in all lines at the same velocity (yellow shaded area).In each panel, a grey line shows the original spectra extracted within the aperture defined for each source, whereas a black line shows the background-subtracted spectra (see Section 3.2).The spectra for all 18 sources are available online as a figure set.

Figure 4 .
Figure 4. JWST NIRCam and MIRI images of the central 1×1 kpc, with locations of ALMA sources marked by black symbols.Top left: NIRCam F335M+F300M+F200W composite image showing the near-IR SED.Bottom left: NIRCam F300M-F335M color image highlighting the 3.35µm PAH emission.These two panels are useful for probing young, embedded clusters (Rodríguez et al. 2023).Top right: MIRI F1000W+F1130W+F770W composite image showing the mid-IR PAH and hot dust emission.Bottom right: MIRI F1000W-F1130W color image highlighting the PAH to continuum contrast.These two panels allow for identifying compact hot dust sources associated with H II regions and/or powered by clusters(Hassani et al. 2023).We find 8 ALMA sources with potential JWST counterparts based on either cross-matching to the 3.35 µm source catalog and/or visual inspection of the images (see Section 3.3 and Table4).

Figure 6 .
Figure6.Gas mass estimates from 350 GHz continuum versus those from HCN (1-0) line (symbols are the same as in Figure5), with the black and gray dotted lines marking the identity relation and a factor of three offset to either side.In addition to statistical uncertainties shown by the error bars, these mass estimates are also subject to ∼0.5 dex systematic uncertainties stemming from our assumptions of dust temperature, emissivity, and abundance as well as HCN conversion factor ( §4.1).The two gas mass estimates broadly agree within the range allowed by those uncertainties.

Figure 7 .
Figure 7.Estimated YMC sizes depend strongly on data resolution.The four datasets being plotted are NGC 3351 at 8 pc resolution (quoted in beam FWHM; this work), the Antennae galaxies at 12 pc (He et al. 2022), NGC 4945 at 2.2 pc (Emig et al. 2020), and NGC 253 at 0.5 pc (Levy et al. 2021).The dashed and dotted lines mark the identity relation and a factor of three above it.

Figure 8 .
Figure 8. YMC total (gas+star) mass versus virial mass, with the latter being calculated from the continuum sizes and HCN (1-0) line widths.The gas reservoirs around most YMCs appear gravitationally unbound or only marginally bound, though the large error budgets on both axes (including statistical and systematic errors; the former shown by the error bars) mean that such inference is inevitably uncertain.

Figure 9 .
Figure 9.Left: ALMA-identified YMCs labeled according to their assigned categories in Section 5.1.The background CO (3-2) image is shaded in a way to highlight structures along the starburst ring (r gal = 200−400 pc).Right: HST-identified young clusters (age < 10 Myr) labeled according to their estimated stellar mass, with the background image similarly shaded.These young clusters (especially the more massive ones) are likely the direct descendants of the ALMA-identified YMCs.

Figure 10 .
Figure10.Two plausible ways of inferring the YMC formation timeline by matching ALMA YMCs (Types 1-4) with young HST clusters (age < 10 Myr).Top: Estimated timeline using HST clusters ≥ 10 4.5 M⊙ as the reference sample, which matches the nominal mass range probed by the ALMA continuum observations.The four phases of YMC formation observable with ALMA are inferred to last for 1.0-1.0-1.8-1.0Myr, respectively.Bottom: Estimated timeline using HST clusters ≥ 10 5 M⊙ as the reference sample, which accounts for a ∼0.5 dex median offset between the 93 GHz-based and optical SED-based stellar mass estimates (see §4.2).The inferred durations of the four phases are 1.7-1.7-2.8-1.7  Myr in this case.
left), though with a slight preference towards their upstream side.When examining trends with the deprojected azimuthal angle ϕ (Figure 11 right), we do not see a clear, coherent trend when considering all YMCs along the ring.There may be ten-

Figure 11 .
Figure 11.Left: Key structural features around the starburst ring, including two streams of bar-driven gas inflow from large radii and the two "contact points" where the inflow gas collides with the ring material.Markers and background image are the same as the left panel in Figure 9. Right: Distributions and properties of the ALMA YMCs and HST clusters, as a function of the deprojected azimuthal angle (increasing counter-clockwise, with zero-point defined by the galaxy position angle).The blue shaded regions indicate azimuthal angle ranges for the Northern/Southern contact points highlighted in the left panel.There is tentative evidence for a progression in YMC evolutionary stage (from Type 1 to Type 4) and gradual decrease in YMC gas to stellar mass ratio along the Eastern side of the ring orbit (i.e., from -180 • to 0 • in azimuthal angle), but the trends are less obvious for the other side.No clear correspondence or variation is seen in the azimuthal distribution of the HST clusters either.

Figure
Figure12.Extracted molecular line spectra for two ALMA sources (#1 and #3).The former represents a non-detection across all lines while the latter is detected in all lines at the same velocity (yellow shaded area).In each panel, a grey line shows the original spectra extracted within the aperture defined for each source, whereas a black line shows the background-subtracted spectra (see Section 3.2).The spectra for all 18 sources are available online as a figure set.

Table 1 .
Basic Properties of the Target Galaxy

Table 4 .
Source Cross-matching Results reported a median synchrotron fraction of 0.36 for YMCs in NGC 4945, Mills et al. 2021 found similar values for those in NGC 253).

Table 5 .
Physical Properties of the ALMA YMC Candidates ID R hl Mgas, 350 Mgas, HCN JS acknowledges support by the National Aeronautics and Space Administration (NASA) through the NASA Hubble Fellowship grant HST-HF2-51544 awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under contract NAS 5-26555.JS also acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference number 568580] through a CITA National Fellowship.The research of HH and CDW is supported by grants from NSERC and the Canada Research Chairs program.The research of KB is supported by a Mitacs Globalink Research Award.RCL acknowledges support for this work provided by the National Science Foundation (NSF) Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102625.The work of AKL is partially supported by the NSF under Grants No. 1615105, 1615109, and 1653300.ECO acknowledges support from grant 510940 from the Simons Foundation.MC gratefully acknowledges funding from the DFG through an Emmy Noether Research Group (grant number CH2137/1-1).COOL Research DAO is a Decentralized Autonomous Organization supporting research in astrophysics aimed at uncovering our cosmic origins.KG is supported by the Australian Research Council through the Discovery Early Career Researcher Award Fellowship (project number DE220100766) funded by the Australian Government.KG is supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (AS-TRO 3D), through project number CE170100013.JDH gratefully acknowledges financial support from