Super Star Clusters in the Central Starburst of NGC 4945

NGC 4945 is a nearby (3.8 Mpc) galaxy hosting a nuclear starburst and Seyfert Type 2 AGN. We use the Atacama Large Millimeter/submillimeter Array (ALMA) to image the 93 GHz (3.2 mm) free-free continuum and hydrogen recombination line emission (H40$\alpha$ and H42$\alpha$) at 2.2 pc (0.12'') resolution. Our observations reveal 27 bright, compact sources with FWHM sizes of 1.4 - 4.0 pc, which we identify as candidate super star clusters. Recombination line emission, tracing the ionizing photon rate of the candidate clusters, is detected in 15 sources, 6 of which have a significant synchrotron component to the 93 GHz continuum. Adopting an age of ~5 Myr, the stellar masses implied by the ionizing photon luminosities are $\log_{10}$($M_{\star}$/M$_{\odot}$) $\approx$ 4.7 - 6.1. We fit a slope to the cluster mass distribution and find $\beta = -1.8 \pm 0.4$. The gas masses associated with these clusters, derived from the dust continuum at 350 GHz, are typically an order of magnitude lower than the stellar mass. These candidate clusters appear to have already converted a large fraction of their dense natal material into stars and, given their small free-fall times of ~0.05 Myr, are surviving an early volatile phase. We identify a point-like source in 93 GHz continuum emission which is presumed to be the AGN. We do not detect recombination line emission from the AGN and place an upper limit on the ionizing photons which leak into the starburst region of $Q_0<10^{52}$ s$^{-1}$.


INTRODUCTION
Many stars form in clustered environments (Lada & Lada 2003;Kruijssen 2012). Bursts of star formation with high gas surface density produce massive Corresponding author: Kimberly L. Emig emig@strw.leidenuniv.nl (> 10 5 M ), compact (FWHM size of 2-3 pc; Ryon et al. 2017) clusters, referred to as super star clusters. Super star clusters likely have high star-formation efficiencies (Goddard et al. 2010;Ryon et al. 2014;Adamo et al. 2011Adamo et al. , 2015Chandar et al. 2017;Johnson et al. 2016;Ginsburg & Kruijssen 2018). They may represent a dominant output of star-formation during the peak epoch of star-formation (z ∼ 1 − 3; Madau & Dickinson 2014). The process by which these massive clusters form now may also relate to the origin of globular clusters.
The earliest stages of cluster formation are the most volatile and currently, unconstrained (Dale et al. 2015;Ginsburg et al. 2016;Krause et al. 2016;Li et al. 2019;Krause et al. 2020). Characterizing properties of young (<10 Myr) clusters is a key step towards understanding their formation, identifying the dominant feedback processes at each stage of cluster evolution, determining which clusters survive as gravitationally bound objects, and linking all of these processes to the galactic environment.
While young clusters of mass ∼10 4 M are found within our Galaxy Longmore et al. 2014;, the most massive, young clusters in the local universe are often found in starbursting regions and merging galaxies (e.g., Zhang & Fall 1999;Whitmore et al. 2010;Linden et al. 2017). Direct optical and even near-infrared observations of forming clusters are complicated by large amounts of extinction. Analyses of optically thin free-free emission and long wavelength hydrogen recombination lines of star clusters offer an alternative, extinction-free probe of the ionizing gas surrounding young star clusters (Condon 1992; Roelfsema & Goss 1992;Murphy et al. 2018). However, achieving a spatial resolution matched to the size of young clusters O(1 pc) (Ryon et al. 2017) in galaxies at the necessary frequencies and sensitivities has only recently become possible thanks to the Atacama Large Millimeter/submillimeter Array (ALMA).
We have recently analyzed forming super star clusters in the central starburst of the nearby (3.5 Mpc) galaxy NGC 253 at ∼2 pc resolution (Leroy et al. 2018;Mills et al. 2020). NGC 4945 is the second object we target in a campaign to characterize massive star clusters in local starbursts with ALMA.
NGC 4945 is unique in that it is one of the closest galaxies (3.8 ± 0.3 Mpc; Karachentsev et al. 2007) where a detected AGN and central starburst coexist. In the central ∼200 pc, the starburst dominates the infrared luminosity and ionizing radiation (Spoon et al. 2000;Marconi et al. 2000), and an outflow of warm ionized gas has been observed (Heckman et al. 1990;Moorwood et al. 1996;Mingozzi et al. 2019). Individual star clusters have not previously been observed in NGC 4945, due to the high extinction at visible and short IR wavelengths (e.g., A V 36 mag; Spoon et al. 2000). Evidence for a Seyfert AGN comes from strong, variable X-ray emission, as NGC 4945 is one of the brightest sources in the X-ray sky and has a Compton thick column density of 3.8×10 24 cm −2 (Marchesi et al. 2018). A kinematic analysis of H 2 O maser emission yields a black hole mass of 1.4 × 10 6 M (Greenhill et al. 1997).
In this article, we use ALMA to image the 93 GHz freefree continuum and hydrogen recombination line emission (H40α and H42α) at 2.2 pc (0.12 ) resolution. This emission allows us to probe photo-ionized gas on star cluster scales and thereby trace ionizing photon luminosities. We identify candidate star clusters and estimate properties relating to their size, ionizing photon luminosity, stellar mass, and gas mass.
Throughout this article, we plot spectra in velocity units with respect to a systemic velocity of v sys = 580 km s −1 in the local standard of rest frame; estimates of the systemic velocity vary by ±25 km s −1 (e.g., Henkel et al. 2018;Chou et al. 2007;Roy et al. 2010). At the distance of 3.8 Mpc, 0.1 corresponds to 1.84 pc.

OBSERVATIONS
We used the ALMA Band 3 receivers to observe NGC 4945 as part of the project 2018.1.01236.S (PI: A. Leroy). We observed NGC 4945 with the main 12 m array telescopes in intermediate and extended configurations. Four spectral windows in Band 3 -centered at 86.2, 88.4, 98.4, and 100.1 GHz -capture the millimeter continuum primarily from free-free emission and cover the hydrogen recombination lines of principal quantum number (to the lower state) n = 40 and n = 42 from the α (∆n = 1) transitions. The rest frequency of H40α is 99.0230 GHz and of H42α is 85.6884 GHz. In this article, we focus on the 93 GHz (λ ∼ 3.2 mm) continuum emission and the recombination line emission arising from compact sources in the starbursting region. We image the data from an 8 km extended configuration, which are sensitive to spatial scales of 0.07 -6 (2-100 pc), in order to focus on the compact structures associated with candidate clusters. We analyze the observatoryprovided calibrated visibilities using version 5.4.0 of the Common Astronomy Software Application (CASA; Mc-Mullin et al. 2007).
When imaging the continuum, we flag channels with strong spectral lines. Then we create a continuum image using the full bandwidth of the line-free channels. We also make continuum images for each spectral window. For all images, we use Briggs weighting with a robust parameter of r = 0.5.
When imaging the two spectral lines of interest, we first subtract the continuum in uv space through a first order polynomial fit. Then, we image by applying a CLEAN mask (to all channels) derived from the fullbandwidth continuum image. Again, we use Briggs weighting with a robust parameter of r = 0.5, which represents a good compromise between resolution and surface brightness sensitivity.
After imaging, we convolved the continuum and line images to convert from an elliptical to a round beam shape. For the full-bandwidth continuum image presented in this article, the fiducial frequency is ν = 93.2 GHz and the final full-width half maximum (FWHM) beam size is θ = 0.12 . The rms noise away from the source is ≈0.017 mJy beam −1 , equivalent to 0.2 K in Rayleigh-Jeans brightness temperature units. Before convolution to a round beam, the beam had a major and minor FWHM of 0.097 × 0.071 .
For the H40α and H42α spectral cubes, the final FWHM beam size is 0.20 , convolved from 0.097 × 0.072 and from 0.11 × 0.083 , respectively. The slightly lower resolution resulted in more sources with significantly detected line emission. We boxcar smoothed the spectral cubes from the native 0.488 MHz channel width to 2.93 MHz. The typical rms in the H40α cube is 0.50 mJy beam −1 per 8.9 km s −1 channel. The typical rms in the H42α cube is 0.48 mJy beam −1 per 10.3 km s −1 channel.
As part of the analysis, we compare the high resolution data with observations taken in a 1 km intermediate configuration as part of the same observing project. We use the intermediate configuration data to trace the total recombination line emission of the starburst. We use the continuum image provided by the observatory pipeline, which we convolve to have a circular beam FWHM of 0.7 ; the rms noise in the full-bandwidth image is 0.15 mJy beam −1 . The spectral cubes have typical rms per channel of 0.24 mJy beam −1 with the same channel widths as the extended configuration cubes. We do not jointly image the configurations because our main science goals are focused on compact, point-like objects. The extended configuration data on their own are well suited to study these objects and any spatial filtering of extended emission will not affect the analysis.
We compare the continuum emission at 3 mm with archival ALMA imaging of the ν = 350 GHz (λ ∼ 850 µm) continuum (project 2016.1.01135.S, PI: N. Nagar). At this frequency, dust emission dominates the continuum. We imaged the calibrated visibilities with a Briggs robust parameter of r = −2 (towards uniform weighting), up-weighting the extended baselines to produce a higher resolution image, suitable for comparison to our new Band 3 data. We then convolve the images to produce a circularized beam, resulting in a FWHM resolution of 0.12 (from an initial beam size of 0.10 × 0.064 ), exactly matched to our 93 GHz continuum image. These data have an rms noise of 0.7 mJy beam −1 (0.2 K). We compare the ALMA data with Australian Long Baseline Array (LBA) imaging of ν = 2.3 GHz continuum emission (Lenc & Tingay 2009). At this frequency and resolution, the radio continuum is predominantly synchrotron emission. We use the Epoch 2 images (courtesy of E. Lenc) that have a native angular resolution slightly higher than the 3 mm ALMA data, with a beam FWHM of 0.080 × 0.032 and an rms noise of 0.082 mJy beam −1 .

CONTINUUM EMISSION
The whole disk of NGC 4945, as traced by Spitzer IRAC 8 µm emission (Program 40410, PI: G. Rieke), is shown in Figure 1. The 8 µm emission predominantly arises from UV-heated polycyclic aromatic hydrocarbons (PAHs), thus tracing the interstellar medium and areas of active star formation. The square box indicates the 8 × 8 (150 pc × 150 pc) starburst region that is of interest in this article. Figure 2 shows the 93 GHz (λ ∼ 3 mm) continuum emission in the central starburst of NGC 4945. Our image reveals ∼30 peaks of compact, localized emission with peak flux densities 0.6-8 mJy (see Section 3.1). On average, the continuum emission from NGC 4945 at this frequency is dominated by thermal, free-free (bremsstrahlung) radiation (Bendo et al. 2016). Freefree emission from bright, compact regions may trace photo-ionized gas in the immediate surroundings of massive stars. We take into consideration the point-like sources detected at 93 GHz as candidate massive star clusters, though some contamination by synchrotrondominated supernova remnants or dusty protoclusters may still be possible. The morphology of the 93 GHz emission and clustering of the peaks indicate possible ridges of star formation and shells. The extended, faint negative bowls flanking the main disk likely reflect the short spacing data missing from this image. We do not expect that they affect our analysis of the point sourcelike cluster candidates.
The large amount of extinction present in this high inclination central region ( i ∼ 72 • ; Henkel et al. 2018) has previously impeded the direct observation of its star clusters. Paschen-α (Pa-α) emission (Marconi et al. 2000) of the n = 3 hydrogen recombination line at 1.87 µm, shown in Figure 3, reveals faint emission above and below the star-forming plane. Corrected for extinction, the clumps of ionized emission traced by Pa-α would give rise to free-free emission below our ALMA detection limit. Pa-α along with mid-infrared spectral lines give support for dust extinction of A V > 160 mag surrounding the AGN core and more generally A V 36 mag in the star-forming region (Spoon et al. 2000).
A large fraction (18/29) of the 93 GHz sources coincide with peaks in dust emission at 350 GHz, as shown in Figure 3. Overall there is a good correspondence between the two tracers. This indicates that candidate clusters are relatively young and may still harbor reservoirs of gas, though in Section 5.7 we find that the fraction of the mass still in gas tends to be relatively small.
In Figure 3, the 93 GHz peaks without dust counterparts tend to be strong sources of emission at 2.3 GHz (Lenc & Tingay 2009), a frequency where synchrotron emission typically dominates. As discussed in Lenc & Tingay (2009), the sources at this frequency are predominantly supernova remnants. The presence of 13 possible supernova remnants -four of which are resolved into shell-like structures with 1.1 to 2.1 pc in diameterindicates that a burst of star-formation activity started at least a few Myr ago. Lenc & Tingay (2009) modeled SEDs of the sources spanning 2.3-23 GHz and found significant opacity at 2.3 GHz (τ = 5 − 22), implying the presence of dense, free-free plasma in the vicinity of the supernova remnants.
At 93 GHz, the very center of the starburst shows an elongated region of enhanced emission (about 20 pc in projected length, or ∼1 ) that is also bright in 350 GHz emission. This region is connected to the areas of highest extinction. Higher column densities of ionized plasma are also present in the region; Lenc & Tingay (2009) observations reveal large free-free opacities at least up to 23 GHz. The brightest peak at 93 GHz, centered at (α, δ) 93 = (13h 05m 27.4798s±0.004s, −49 • 28 05.404 ± 0.06 ), is co-located with the kinematic center as determined from H 2 O maser observations (α, δ) H2O = (13h 05m 27.279s±0.02s, −49 • 28 04.44 ±0.1 ) (Greenhill et al. 1997) and presumably harbors the AGN core. We refer to the elongated region of enhanced emission surrounding the AGN core as the circumnuclear disk. The morphological similarities between 93 GHz and 350 GHz, together with the detection of a synchrotron point source (likely supernova remnant; see Section 3.2 and Source 17) in the circumnuclear disk, indicate starformation is likely present there.

Point Source Identification
We identify candidate star clusters via point-like sources of emission in the 93 GHz continuum image. Sources are found using PyBDSF (Mohan & Rafferty 2015) in the following way. Islands are defined as contiguous pixels (of nine pixels or more) above a threshold of seven times the global rms value of σ ≈ 0.017 mJy beam −1 . Within each island, multiple Gaussians may be fit, each with a peak amplitude greater than the peak threshold, a threshold of ten times the global rms. We chose this peak threshold to ensure that significant emission can also be identified in the continuum images made from individual spectral windows. The number of Gaussians is determined from the number of distinct peaks of emission higher than the peak threshold and which have a negative gradient in all eight evaluated directions. Starting with the brightest peak, Gaussians are fit and cleaned (i.e., subtracted). A source is identified with a Gaussian as long as subtracting its fit does not increase the island rms.
Applying this algorithm to our 93 GHz image yielded 50 Gaussian sources. We remove 5 sources that fall outside of the star-forming region. We also removed 3 sources that appeared blended, with an offset < 0.12 from another source. Finally we remove 13 sources that do not have a flux density above ten times the global rms after extracting the 93 GHz continuum flux density through aperture photometry (see Section 3.2). As a result, we analyze 29 sources as candidate star clusters. In Figure 3, we show the location of each source with the apertures used for flux extraction. The sources match well with what we would identify by eye.

Point Source Flux Extraction
For each source, we extract the continuum flux density at 2.3 GHz, 93 GHz, and 350 GHz through aperture photometry. Before extracting the continuum flux at 2.3 GHz, we convolve the image to the common resolution of 0.12 . We extract the flux density at the location  Figure 2. ALMA 93 GHz (λ ∼ 3.2 mm) continuum emission in the central starburst of NGC 4945. The continuum at this frequency is dominated by ionized, free-free emitting plasma. In this paper, we show that the point-like sources are primarily candidate, massive star clusters. The brightest point source of emission at the center is presumably the Seyfert AGN. The rms noise away from the source is σ ≈ 0.017 mJy beam −1 and the circularized beam FWHM is 0.12 (or 2.2 pc at the distance of NGC 4945). Contours of the continuum image show 3σ emission (gray) and [4, 8, 16, ...256]σ emission (black). of the peak source within an aperture diameter of 0.24 . Then we subtract the extended background continuum that is local to the source by taking the median flux density within an annulus of inner diameter 0.24 and outer diameter 0.30 ; using the median suppresses the influence of nearby peaks and the bright surrounding filamentary features. The flux density of each source at each frequency is listed in Table 1. When the extracted flux density within an aperture is less than three times the global rms noise (in the 2.3 GHz and 350 GHz images), we assign a three sigma upper limit to that flux measurement.
In Figure 4, we plot the ratio of the flux densities extracted at 350 GHz and 93 GHz (S 350 /S 93 ) against the ratio of the flux densities extracted at 2.   Table 1. Circles show apertures (diameter of 0.24 ) used for continuum extraction. Their colors indicate the measured in-band spectral index, as in Figure 4 where dark purple indicates synchrotron dominated emission and yellow indicates dust dominated emission. Top right: HST Paschen-α emission -hydrogen recombination line, n = 3, at 1.87 µm -(courtesy P. van der Werf) tracing ionized gas at ≈ 0.2 resolution (Marconi et al. 2000). Dust extinction, of AV > 36 mag, obscures the Pa-α recombination emission at shorter wavelengths from the starburst region. Contours trace 93 GHz continuum, as described in Figure 2. Bottom left: ALMA 350 GHz continuum emission tracing dust. Bottom right: Australian LBA 2.3 GHz continuum imaging of synchrotron emission primarily from supernova remnants (Lenc & Tingay 2009). in the middle of the plot. One exception is the AGN (Source 18) which, due to self-absorption at frequencies greater than 23 GHz, is bright at 93 GHz but not at 2.3 GHz and therefore has the lowest S 2.3 /S 93 ratio. From the continuum measurements we construct simple SEDs for each source. These SEDs are used for illustrative purposes and do not affect the analysis in this paper. Examples of the SEDs of three sources are included in Figure 5. We show an example of a free-free dominated source (Source 22), which represents the majority of sources, as well as a dust (Source 12) and a synchrotron (Source 14) dominated source. Of the sources with extracted emission of >3σ at 2.3 GHz, nine also have the free-free absorption of their synchrotron spectrum modeled. We plot this information whenever possible. When a source identified by Lenc & Tingay (2009) lies within 0.06 (half the beam FWHM) of the 93 GHz source, we associate the low-frequency modeling with the 93 GHz source. We take the model fit by Lenc & Tingay (2009) and normalize it to the 2.3 GHz flux that we extract -as an example, see the solid purple curve in the middle panel of Figure 5. The SEDs of all sources are shown in Figure 13 in Appendix C.

Free-free Fraction at 93 GHz
The flux density of optically thin, free-free emission at millimeter wavelengths (see Appendix A.2; Draine 2011) arises as S ff =(2.08 mJy) n e n + V 5 × 10 8 cm −6 pc 3 × (1) where EM C = n e n + V is the volumetric emission measure of the ionized gas, D is the distance to the source, and T e is the electron temperature of the medium. Properties of massive star clusters (i.e., ionizing photon rate) can thus be derived through an accurate measurement of the free-free flux density and an inference of the volumetric emission measure. We will determine a free-free fraction, f ff , and let S ff = f ff S 93 .
In this section we focus on determining the portion of free-free emission that is present in the candidate stars clusters at 93 GHz. To do this, we need to estimate and remove contributions from synchrotron emission and dust continuum. We determine an in-band spectral index across the 15 GHz bandwidth of the ALMA Band 3 observations. Using the spectral index 1 , we constrain 1 similar methods have been used by e.g., Linden et al. (2020)   the free-free fraction as well as the fractional contributions of synchrotron and dust (see Figure 4 and Table 1). We found the in-band index to give stronger constraints and, for 2.3 GHz, to be more reliable than extrapolating (assuming indices of −0.8,−1.5) because of the two decades difference in frequency coupled with the large optical depths already present at 2.3 GHz.

Band 3 Spectral Index
We estimate an in-band spectral index at 93 GHz, α 93 , from a fit to the flux densities in the spectral window continuum images of the ALMA Band 3 data. The spectral windows span 15 GHz, which we set up to have a large fractional bandwidth. We extract the continuum flux density of each source in each spectral window using aperture photometry with the same aperture sizes Note-RA and Dec refer to the center location of a Gaussian source identified in the 93 GHz continuum image, in units of hour angle and degrees, respectively. S93 is the flux density in the 93 GHz full bandwidth continuum image. α93 is the spectral index at 93 GHz (S ∝ ν α ), as determined from the best fit slope to the 85-101 GHz continuum emission. S2.3 is the flux density extracted in the 2.3 GHz continuum image. S350 is the flux density extracted in the 350 GHz continuum image. f ff , fsyn, f d are the free-free, synchrotron, and dust fractional contribution to the 93 GHz continuum emission, respectively, as determined from the spectral index. See Section 3.3 for a description of how the errors on these fractional estimates are determined. a A 3σ upper limit to the sources undetected in 2.3 GHz continuum emission is 1.2 mJy. b The error on the 350 GHz flux density measurement is 3.2 mJy. A 3σ upper limit to the undetected sources is 9.6 mJy. c The error on the fractional contributions are the same as for f ff unless otherwise noted. as described in Section 3.2. We fit a first order polynomial to the five continuum measurements -four from each spectral window, one from the full bandwidth image. The fit to the in-band spectral index is listed in Table 1 with the one sigma uncertainty to the fit. The median uncertainty of the spectral indices is 0.13. The spectral indices measured in this way were consistent for the brightest sources with the index determined by CASA tclean; however, we found our method to be more reliable for the fainter sources. Nonetheless, the errors are large for faint sources.

Decomposing the fractional contributions of emission type
From the in-band spectral index fit, we estimate the fractional contribution of each emission mechanismfree-free, synchrotron, and (thermal) dust -to the 93 GHz continuum. To do this, we simulate how mixtures of synchrotron, free-free, and dust emission could combine to create the observed in-band index. We assume fixed spectral indices for each component and adjust their fractions to reproduce the observations (see below).
We consider that a source is dominated by two types of emission (a caveat which we discuss in detail in Section 6): free-free and dust, or free-free and synchrotron. With synthetic data points, we first set the flux density at 93 GHz, S 93 , to a fixed valueand vary the contributions of dust and free-free continua, such that S 93 = S d + S f f . For each model, we determine the continuum flux at each frequency across the Band 3 frequency coverage as the sum of the two components. We assume that the frequency dependence of the free-free component is α ff = −0.12 and the frequency dependence of the dust component is α d = 4.0. We fit the (noise-less) continuum of the synthetic data across the Band 3 frequency coverage with a power-law, determining the slope as α 93 , the in-band index. We express the results in terms of a free-free fraction and dust fraction, rather than absolute flux. In this respect, we explore free-free fractions to the 93 GHz continuum ranging from f ff = (0.001, 0.999) in steps of 0.001.We let the results of this process constrain our free-free fraction when the in-band index measured in the actual (observed) data is α 93 ≥ −0.12.
Next, we repeat the exercise, but we let synchrotron and free-free dominate the contribution to the continuum at 93 GHz. We take the frequency dependence of the free-free component as α ff = −0.12 across the Band the dust emission is optically thin with a wavelengthdependent emissivity so that τ ∝ λ −2 (e.g., see Draine 2011). The low optical dust optical depths estimated in Section 5.6 imply a dust spectral index steeper than 2, though the exact value might not be 4, e.g., if our assumed emissivity power law index is not applicable. The generally faint dust emission indicates that our assumptions about dust do not have a large effect on our results. For the synchrotron frequency dependence, we assume α syn = −1.5. This value is consistent with the best fit slope of α syn = −1.4 found by Bendo et al. (2016) averaged over the central 30 of NGC 4945. Furthermore, the median slope of synchrotron-dominated sources modeled at 2.3-23 GHz is -1.11 (Lenc & Tingay 2009), indicating that even at 23 GHz the synchrotron spectra already show losses due to aging, i.e., are steeper than a canonical initial injection of α syn ≈ −0.8. While our assumed value of α syn = −1.5 is well-motivated on average, variations from source to source are likely present.
Through this method of decomposition, the approximate relation between the in-band index and the freefree fraction is for which α 93 ≤ −0.12, the synchrotron fraction is found to be f syn = 1 − f ff and we set f d = 0, and for which α 93 ≥ −0.12 the dust fraction is found to be f d = 1 − f ff and f syn = 0. This relation is depicted in the bottom panel of Figure 4 using the values that have been determined for each source.
3.3.3. Estimated free-free, synchrotron, and dust fractions to the 93 GHz continuum Table 1 and Figure 4 summarize our estimated fractional contribution of each emission mechanism to the 93 GHz continuum of each source. We use the same relation above to translate the range of uncertainty on the spectral index to an uncertainty in the emission fraction estimates. We find, at this 0.12 resolution, the median free-free fraction of sources is f ff = 0.62 with median absolute deviation of 0.29. Most of the spectral indices are negative, and as a result, we find synchrotron emission can have a non-trivial contribution with a median fraction of f syn = 0.36 and median absolute deviation of 0.32. On the other hand, three of the measured spectral indices are positive. The median dust fraction of sources is f d = 0 and median absolute deviation of 0.10. A 1σ limit on the fractional contribution of dust does not exceed f d = 0.47 for any single source.
Additional continuum observations of comparable resolution at frequencies between 2 GHz and 350 GHz would improve the estimates of the fractional contribution of free-free, dust and synchrotron to the 93 GHz emission.

RECOMBINATION LINE EMISSION
Hydrogen recombination lines at these frequencies trace ionizing radiation (E > 13.6 eV); this recombination line emission is unaffected by dust extinction. The integrated emission from a radio recombination line transition to quantum number n, which we derive for millimeter wavelength transitions in Appendix A.1, is described by where b n+1 is the LTE departure coefficient, EM L = n e n p V is the volumetric emission measure of ionized hydrogen, D is the distance to the source, T e is the electron temperature of the ionized gas, and ν is the rest frequency of the spectral line. Figure 6 shows spectra of the 15 sources with detected radio recombination line emission. We extract H40α and H42α spectra at the location of each source with an aperture diameter of 0.4 , or twice the beam FWHM. We average the spectra of the two transitions together to enhance the signal-to-noise ratio of the recombination line emission. To synthesize an effective H41α profile, we interpolate the two spectra of each source to a fixed velocity grid with a channel width of 10.3 km s −1 , weight each spectrum by σ −2 rms where σ rms is the spectrum standard deviation, and average the spectra together lowering the final noise. The averaged spectrum has an effective transition of H41α at ν eff = 92.034 GHz.
We fit spectral features with a Gaussian profile. We calculate an integrated signal-to-noise ratio for each line by integrating the spectrum across the Gaussian width of the fit (i.e., ±σ Gaus ) and then dividing by the noise over the same region, √ N σ rms , where N is the number of channels covered by the region. We report on detections with an integrated signal of > 5σ rms . Table 2 summarizes the properties of the line profiles derived from the best-fit Gaussian. The median rms of the spectra is σ rms = 0.34 mJy.
In Figure 12 of Appendix B, we show that the central velocities of our detected recombination lines are in good agreement with the kinematic velocity expected of the disk rotation. To do this, we overlay our spectra on  Figure 6. Radio recombination line spectra for sources with significantly detected emission. The thin blue line is the H40α spectrum. The thin green line is the H42α. These spectra have been regridded from their native velocity resolution to the common resolution of 10.3 km s −1 . The thick black line is the weighted average spectrum of H40α and H42α, effectively H41α. In red is the best fit to the effective H41α radio recombination line feature.
H40α spectra extracted from the intermediate configuration observations (0.7 resolution).
In the top panel of Figure 7, we show the total recombination line emission extracted from the starburst region in the 0.2 resolution, "extended" configuration observations. The aperture we use, designated as re-gion T1, is shown in Figure 8. Details of the aperture selection are described in Section 4.1.1. The spectrum consists of two peaks reminiscent of a double horn profile representing a rotating ring. We fit the spectrum using the sum of two Gaussian components. In Table 3, we include the properties of the best fit line profiles. The sum total area of the fits is (2.1 ± 0.6) Jy km s −1 .

Line Emission from 0.7 resolution, Intermediate-configuration Observations
Figure 7 also shows integrated spectra derived from intermediate-resolution (0.7 ) data. We use these data as a tracer of the total ionizing photons of the starburst region. We expect that the intermediate-resolution data includes emission from both discrete, point-like sources and diffuse emission from any smooth component. Note-Vcen is the central velocity of the best fit Gaussian. Peak is the peak amplitude of the Gaussian fit. FWHM is the full-width half maximum of the Gaussian fit. σrms is the standard deviation of the fit-subtracted spectrum.
The total integrated emission in the intermediateconfiguration data is about three times larger than the integrated emission in the extended-configuration data. Spectra representing the total integrated line flux are shown in Figure 7. In Figure 8, we show the integrated intensity map of H40α emission from the intermediateconfiguration (0.7 ) data. In Table 3, we include the best fit line profiles.
We also compare the line profiles of H40α and H42α in the intermediate-configuration (0.7 ) data, see Ta Figure 8). In blue is the extended-configuration ("extended") spectrum extracted from the high-resolution 0.2 data; this spectrum shows the maximum total integrated line flux extracted. In purple, the intermediate-configuration ("intermed") spectra extracted from low-resolution, native 0.7 data; the spectrum from the T2 region is the total maximum integrated line flux from this data. The solid black line represents the sum total of two Gaussian fits. The dotted, black line represents the single Gaussian fits. throughout the starburst region and up to 30 pc in apparent size beyond the region where we detect the bright point sources at high resolution.
Also shown in Figure 8 are the apertures used to extract spectra in Figure 7.
We fit a two dimensional Gaussian to the continuum emission in the 0.7 resolution observations (see Figure 9). This results in a best fit centered at (α, δ) = (13h 05m 27.4896s, −49 • 28 05.159 ), with major and minor Gaussian widths of σ maj = 2.4 and σ min = 0.58 , and an angle of θ = 49.5 • ; we use this fit as a template for the aperture location, position angle, width and height. We independently vary the major and minor axes (in multiples of 0.5σ maj and 0.5σ min , respectively) in order to determine the aperture which maximizes the total integrated signal in channels within ±170 km s −1 . With the extended-configuration cube, we find the largest integrated line emission with an aperture of 8.4 × 1.5 , which we refer to as T1. With the intermediate-configuration cube, the largest integrated  line emission arises with an aperture of 12.0 × 2.9 , which we refer to as T2.
We extract the total H40α line flux from the intermediate-configuration (0.7 resolution) cube using the T2 aperture (see bottom panel of Figure 7). The spectrum shows a double horn profile, indicating ordered disk-like rotation. We fit the features with two Gaussians. The sum total of their integrated line flux is (6.4 ± 1.0) Jy km s −1 .
We also extract H40α line flux from the intermediateconfiguration (0.7 resolution) cube within the T1 aperture in order to directly compare the integrated line flux in the two different data sets using the same aperture regions. We find more emission in the intermediateconfiguration data, a factor of ∼2.3 greater than the extended-configuration data. This indicates that some recombination line emission originates on large scales (>100 pc) to which the high-resolution, long baselines are not sensitive.

H42α Contamination
In this section we compare the line profiles from H40α and H42α extracted from the 0.7 intermediate configuration data. In principle, we expect the spectra to be virtually identical, which is why we average them to improve the signal-to-noise at high-resolution. Here we test that assumption at low-resolution. To summarize, we find evidence that a spectral line may contaminate the H42α measured line flux in broad (typically spatially unresolved) line profiles. Yet we see good agreement between the two lines at the scale of individual cluster candidates.
We extracted spectra in three apertures to demonstrate the constant velocity offset of the contaminants. We approximately matched the locations of these apertures to those defined in Bendo et al. (2016), in which H42α was analyzed at 2.3 resolution; in this way we are able to confirm the flux and line profiles we extract at 0.7 resolution with those at 2.3 . The non-overlapping circular apertures with diameters of 4 designated as North (N), Center (C), and South (S) are shown in Figure 9.
We used our intermediate-configuration data to extract an H40α and an H42α spectrum in each of the regions. We overplot the spectra of each region in Figure 10. We fit a single Gaussian profile to the line emission, except for H42α emission in region N where two Gaussian components better minimized the fit. The total area of the fits are presented in Table 4 as the integrated line flux.
Our line profiles of H42α are similar in shape and velocity structure as those analyzed in Bendo et al. (2016) and the integrated line emission is also consistent (within 2σ). This indicates that we are recovering the H42α total line flux and properties with our data.
On the other hand, the H40α flux we extract is about a factor of ∼1.6 lower than the H42α fluxes in these apertures (see Figure 10 and Table 4). The discrepancy  grows to a factor of 2 in the profile extracted from the total region. The additional emission seen in the H42α spectrum at velocities +100 km s −1 to +250 km s −1 with respect to the bright, presumably hydrogen recombination line peak, is absent in the H40α profile. It is not likely to be a maser-like component of hydrogen recombination emission since the relative flux does not greatly vary in different extraction regions, and densities outside of the circumnuclear disk would not approach the emission measures necessary (e.g., EM v 10 10 cm −6 pc 3 ) for stimulated line emission.
We searched for spectral lines in the frequency range ν rest ∼ 85.617 -85.660 GHz, corresponding to these  Figure 9 as N (top), C (middle), and S (bottom). We find H42α to be contaminated by spectral lines which may include c-C3H2 432 − 423 -shown as a dashed line in the panels at expected velocities with respect to H42α. When contaminant lines are included, the integrated line flux of H42α is over estimated by a factor of 1.5 in these apertures; this grows to a factor of 2 when integrating over the total starburst emission.
velocities and find several plausible candidates, though we were not able to confirm any species with additional transitions in the frequency coverage of these observations. A likely candidate may be the 4 32 − 4 23 transition of c-C 3 H 2 . c-C 3 H 2 has a widespread presence in the diffuse ISM of the Galaxy (e.g., Lucas & Liszt 2000) and the 2 20 − 2 11 transition has been detected in NGC 4945 (Eisner et al. 2019). As an example we plot the velocity of c-C 3 H 2 4 32 − 4 23 relative to H42α in Figure 10.

PHYSICAL PROPERTIES OF THE CANDIDATE STAR CLUSTERS
In this section, we estimate properties of the candidate star clusters, summarized in Table 6. We discuss their size and approximate age. Properties of the ionized gas content, such as temperature (see Table 5), metallicity, density and mass are derived from the continuum and recombination line emission. We estimate the ionizing photon rate of the candidate stars clusters and use it to infer the stellar mass (see Figure 11). From the dust emission at 350 GHz, we estimate gas masses of the candidate star clusters. With a combined total mass from gas and stars, we estimate current mass surface densities and free-fall times.
We exclude Source 5 from the analysis since the freefree fraction is f ff < 0.01. We also remove the presumed AGN core (Source 18) from the analysis.

Size
The sources identified through PyBDSF in the 93 GHz continuum image are fit with two dimensional Gaussians. The average of the major and minor (convolved) FWHM is listed in Table 6 as the FWHM size of the source in units of pc. The Gaussian fits are all consistent with circular profiles within error. FWHM sizes of (1.4-4.0) pc are observed, consistent with typical sizes of young, massive star clusters (Ryon et al. 2017;Leroy et al. 2018). However, the lower end may reflect the resolution limit of our beam, with a FWHM size of 2.2 pc. The uncertainties we report reflect the errors of the Gaussian fit.
Based on high-resolution imaging of embedded clusters in the nucleus of the Milky Way and NGC 253, some of these clusters might break apart at higher resolution (Ginsburg et al. 2018, Levy et al., in prep). If they follow the same pattern seen in these other galaxies, each source would have one or two main components potentially with several associated fainter components.

Age
Throughout our analysis, we assume that the candidate star clusters formed in an instantaneous burst of star-formation roughly 5 Myr ago (with a likely uncertainty of ∼1 Myr). This (approximately uniform) age is supported through the coincident detection of RRLs and supernovae remnants, previous analyses of the global population of the burst (e.g., Marconi et al. 2000;Spoon et al. 2000) and an orbital timescale of ≈3 Myr for the starburst region. We elaborate on this supporting evidence below.
As we discuss in Section 3.3, dust does not significantly contribute to 93 GHz emission (with a median fraction of f d = 0), but synchrotron emission does through supernova remnants. Supernova explosions begin from ∼3 Myr in the lifetime of a cluster and cease around ∼40 Myr when the most massive stars have died out; this puts the loosest bounds on the age of the candidate clusters we observe. The coincident detection of supernovae remnants in a third (6/15) of our recombination line detected sources implies that the burst is likely not at the earliest stage of the supernovae phase. However the ionizing photon rate changes dramatically over 3 Myr to 10 Myr, dropping by about two orders of magnitude (Leitherer et al. 1999). As a result clusters are significantly harder to detect in radio recombination lines or free-free continuum emission after ∼5 Myr.
Properties of star-forming activity in the central starburst have been estimated by combining far-infrared (FIR) and optical/IR tracers. Marconi et al. (2000) discerned an age of 6 Myr and mass of 4 × 10 7 M by using Paα and Brγ to trace the energy distribution of the photon output of the population. However, the dust extinction was underestimated, complicated by the uncertainty in the AGN contribution. Mid-infrared (MIR) observations with the Infrared Space Observatory (ISO; Kessler et al. 1996) of line ratios further constrained this scenario. Spoon et al. (2000) estimated an extinction of A V = 36 +18 −11 , determined that the AGN is not dominating the ionizing radiation field, and found that the star-forming population is consistent with a burst of age ≥5 Myr.
As a sanity check on whether a synchronized burst might be expected, we calculate the orbital timescale associated with the the burst region. Taking the rotation velocity ∼ 170 km s −1 from the integrated spectrum and the radius ∼ 80 pc associated with region T1, we estimate an orbital timescale of ∼ 3 Myr. If we take this as roughly the timescale for the nuclear disk to react to changing conditions, a burst shutting off or turning on in a ∼ 5 Myr timescale is reasonable.

Temperature and Metallicity
The ratio of the integrated recombination line flux (Equation 3) to the free-free continuum flux density (Equation 1) allows the electron temperature to be determined. Dependencies on the distance, emission measure, and (possible) beam-filling effects cancel out under the assumption that the two tracers arise in the same volume of gas. We show in Appendix A.3.1, that when taking the ratio of the integrated line to continuum, R LC , and solving for the temperature, T e , we arrive at T e = 10 4 K b n+1 (1 + y)  where b n+1 is the non-LTE departure coefficient, and y is the abundance ratio of ionized helium to hydrogen number density, y = n He+ /n p , which we fix as y = 0.10 (de Pree et al. 1996;Mills et al. 2020). Table 5 lists the temperatures we derive in the region. We focus on the 5 sources with bright (peak S/N > 4.7σ) and well-fit recombination line emission. Most of these sources have higher free-free fractions than the median. To derive the temperatures, we re-evaluate the continuum (fraction of) free-free emission at the resolution of 0.2 , since the free-free fraction may change with resolution. Therefore, we convolve the Band 3 continuum images to 0.2 resolution. We extract the continuum from the full-bandwidth image through aperture photometry, using an aperture diameter of 0.4 . In order to exactly match the processing of the spectral line data, we do not subtract background continuum emission within an outer annulus. Then, by extracting the continuum in each spectral window (using the same aperture diameters just described), we fit for the in-band spectral index. We use the procedure described in Section 3.3 to constrain the free-free fraction from the spectral index fit.
With the free-free fraction and measured fluxes, we plug in the line to continuum ratio into Equation 4 and take b n = 0.73 (Storey & Hummer 1995) to arrive at the temperature. The departure coefficient at n = 41 is loosely (< 15% variation) dependent on the temperature. We iterate (once) on the input b n and output temperature. b n = 0.73 is the modeled value for this temperature and for typical densities of n e = (10 3 − 10 4 ) cm −3 of ionized gas surrounding young, massive stars and consistent with the ionized gas densities we derive in Section 5.4.
The uncertainties in the electron temperatures we derive in Table 5 are dominated by the uncertainties in the free-free fraction. We take the mean and standard deviation values of T e = (6000 ± 400) K as a representative electron temperature of the ionized plasma in the candidate star clusters. This temperature is consistent with the temperature derived from a lower-resolution analysis of NGC 4945 at 2.3 × 2.6 resolution, which finds T e = (5400 ± 600) K (Bendo et al. 2016).
Our estimated temperature implies a thermal line width of (16 ± 4) km s −1 (Brocklehurst & Seaton 1972). Given that this is smaller than our observed line widths, non-thermal motions from bulk velocities (such as turbulence, inflow or outflow) must contribute to broadening the spectral line profiles.
The electron temperature of free-free plasma surrounding massive stars is related to the metallicity of the plasma, as the metals contribute to gas cooling. Shaver et al. (1983) with the temperatures and metallicities derived with (auroral) collisionally excited lines at optical wavelengths. Furthermore, they showed that these temperatures are consistent with electron temperatures derived from radio recombination lines. We find a representative O/H metallicity of 12 + log 10 (O/H) = 8.9 ± 0.1. This value is in approximate agreement (within 2σ) with the average metallicity and standard deviation of 12+log 10 (O/H) = 8.5±0.1 (Stanghellini et al. 2015) determined in 15 star-forming regions in the galactic plane of NGC 4945 (and which is consistent with no radial gradient) using strong-line abundance ratios of oxygen, sulfur, and nitrogen spectral lines.

Ionized Gas: Emission Measure, Density and Mass
We determine the volumetric emission measure of gas ionized in candidate stars clusters using Equations 1 and 3 together with the mean temperature derived in Section 5.3. In Table 6, we list the results for each candidate star cluster. Emission measures that we determine from the free-free continuum range from log 10 (EM C /cm −6 pc 3 ) ∼ 7.3 -8.7, with a median value of 8.4. We also calculate the volumetric emission measure of ionized hydrogen as determined by the effective H41α recombination line when applicable, noting that EM C = (1 + y) EM L . The line emission measures range from log 10 (EM L /cm −6 pc 3 ) ∼ 8.4 -8.9, with a median value of 8.5. The uncertainty in the emission measures is ∼0.4 dex and is dominated by the errors of the free-free fraction.
Next, we solve for the electron density. We use the emission measure determined from the free-free continuum, assume n e = n + , and consider a spherical volume with r = FWHM size /2. We arrive at densities between log 10 (n e /cm −3 ) = 3.1-3.9 with a median value of 3.5.
We matched (see Section 3.2) five of the candidate star clusters that have recombination line emission detected -Sources 1, 6, 14, 21, 27 -with the 2.3 GHz objects of Lenc & Tingay (2009) which have the freefree optical depth modeled through their low-frequency turnovers. Although the 2.3 GHz objects have nonthermal indices, it is their radio emission which is opaque to free-free plasma. Using the optical depths derived in Lenc & Tingay (2009) and our fiducial electron temperature, we solve for the density through the relation, τ ≈ 3.28 × 10 −7 Te (Condon & Ransom 2016), where EM = n e n + and for which a spherical region the pathlength translates as = 3 4 r. We find densities in the range log 10 (n e /cm −3 ) = 3.3 -3.6. This agrees well with the values we separately derive.
We convert the ionized gas density and source sizes to an ionized gas mass through, where we have assumed a 1.36 contribution of helium by mass and we let r = FWHM size /2. The ionized gas masses of the candidate star clusters range from log 10 (M + / M ) = 2.7 -3.5 with a median value of 3.1. The ionized gas mass is a small fraction ( 1%) of the stellar mass (see Section 5.5).

Ionizing Photon Production and Stellar Mass
We estimate the number of the ionizing photons needed per second to maintain the total free-free emitting content (see Table 6). From the emission measure of ionized gas and the temperature-dependent recombination coefficient for case B recombination, the rate of ionizing photons (see Appendix A.3.2) with E > 13.6 eV is Q 0 = 3.8 × 10 51 s −1 n e n + V 5 × 10 8 cm −6 pc 3 × where EM C = n e n + V is the volumetric emission measure of the total ionized gas which we take from the continuum derived emission measure, and T e is the electron temperature of the ionized gas. Our candidate star clusters have ionizing photon rates in the range log 10 (Q 0 /s −1 ) ∼ 50.4 -51.8. The sum of the ionizing photon rate over all candidate, massive star clusters is 5.3 × 10 52 s −1 . In the top panel of Figure 11, the ionizing photon rates of the candidate clusters are plotted as complementary cumulative fractions.
We use Starburst99 calculations (Leitherer et al. 1999) to infer the stellar mass from the ionizing photon output of a 5 Myr old stellar population, via We arrive at this value by simulating a single 10 6 M stellar population, with the initial mass function (IMF) of Kroupa (2001), a maximum stellar mass of 100 M , and the default stellar evolution tracks and tuning parameters. Then we divide the ionizing photon output at 5 Myr by the initial mass of the stellar population. We note that this is a rough approximation which has not accounted for the amount of ionizing photons absorbed by dust, mass ejected from the system, and/or enhanced emission from stellar binaries. Our candidate star clusters have stellar masses in the range log 10 (M /M ) ∼ 4.7-6.1 (see Table 6) with a median of 5.5. The error on the mass estimate is ∼0.4 dex. The sum of the stellar masses of the candidate stars clusters is ≈ 1.1×10 7 M . In the bottom panel of Figure 11, the estimated stellar masses of the candidate clusters are plotted as cumulative fractions.

Gas Mass from Dust
We estimate the mass of gas associated with each candidate star cluster (see Table 6) from dust emission at 350 GHz. We determine the dust optical depth by comparing the measured intensity with that expected from an estimate of the true dust temperature. Assuming a mass absorption coefficient, we convert the optical depth to a dust column density. We arrive at a gas mass by multiplying the dust column density with the measured source size and an assumed dust-to-gas mass ratio.
We assume a dust temperature of T dust = 130 K, as has been determined for the gas kinetic temperature in the forming super star clusters in NGC 253 (Gorski et al. 2017). This is an approximation, though the uncertainty is linear. Then we convert the 350 GHz flux density into an intensity (I 350 ), and solve for the optical depth through I 350 ≈ τ 350 B ν (T dust ) where B ν (T dust ) is the Planck function evaluated at 350 GHz. We measure optical depths in the range τ 350 ∼ 0.02 -0.10, with a median value of τ 350 ∼ 0.04, justifying our optically thin assumption. We note that the 3σ upper limit of the sources which have not been detected at 350 GHz corresponds to τ 350 < 0.02. Note-Source 5 is not included since its free-free fraction is f ff < 0.01; Source 18 is not included since it is the AGN core. FWHM is the source size as best fit from a Gaussian (to flux that has not been deconvolved); the errors reflect the fit of the Gaussian. EMC is the free-free emission measure derived from the continuum as in Equation 1. EML is the hydrogen free-free emission measure derived from effective H41α as in Equation 3; we note EMC = (1 + y) EML. Q0 is the ionizing photon rate derived from EMC as in Equation 7. M is the stellar mass derived from Q0 as in Equation 8. a The error on these quantities is ∼0.4 dex.
Next, we convert the optical depth to a dust column density using an assumed mass absorption coefficient (κ). We adopt κ = 1.9 cm 2 g −1 which should be appropriate for ν ∼ 350 GHz and dust mixed with gas at a density of ∼ 10 5 − 10 6 cm −3 (Ossenkopf & Henning 1994), but we do note the large (factor of 2) uncertainties on this value. Finally, we combine the dust surface density with an adopted dust-to-gas mass ratio (DGR) of 1-to-100, approximately the Milky Way value and similar to the value found for starburst galaxies by Wilson et al. (2008). Our estimate for the gas surface density is determined with: We determine the gas mass by multiplying the gas surface density by the two dimensional area of the source size, M gas = A Σ gas . The gas masses we estimate are included in Table 6. We find values in the range of log 10 (M gas / M ) = 4.4 -5.1 with a median value of 4.7. Upper limits for the sources which have not been detected in 350 GHz emission are included in the Table. Table 6) and in NGC 253 (Mills et al. 2020). The turn-off at lower values likely reflects completeness limits. Bottom: The stellar mass, M , of our candidate star clusters inferred from the ionizing photon rates of a cluster with an age of 5 Myr, plotted as a complementary cumulative distribution. We also include the stellar masses of clusters in the starburst of NGC 253 (Mills et al. 2020) and in the galaxies LMC, M 51, and the Antennae (Mok et al. 2020). Since the clusters in NGC 253 are likely close to a zero age main sequence, they would produce more ionizing photons per unit mass as compared with the slightly older stellar population in the clusters of NGC 4945.

Total Mass from Gas and Stars
We estimate the current total mass of each candidate star cluster as the sum of the gas and star masses, where M Tot = M gas + M . The total mass is dominated by the stellar mass, as we find low gas mass fractions of M gas /M Tot = 0.04 -0.22 and a median of 0.13. When calculating the total mass of the 10 sources which are not detected at 350 GHz, we do not include the lower limit of the dust mass; we only consider the stellar mass.
We express the total mass of each source in terms of a surface density (in Table 6). This is calculated within the FWHM of the region; we thereby divide the total mass by 2 and divide by the 2D area of the FWHM. The values range from log (Σ Tot ) = 3.1 -4.8 M pc −2 , or alternatively, Σ Tot = 0.3 -13 g cm −2 .
Using the total mass, we estimate the gravitational free fall time of the clusters, t ff = πr 3 8GMTot 1/2 . The values we derive are included in Table 6, and range from log 10 (t ff / yr ) = 4.3 -5.2 with a median of 4.6. This is the gravitational free fall time that would be experienced by gas with no support if all of the cluster mass were gas. The fact that the gas mass fractions are low and that the age appears longer than the free fall time adds support to the idea that these clusters have mostly already formed.
Summing all sources, we find a total mass within the candidate clusters of M Tot ≈ 1.5 × 10 7 M .

DISCUSSION
Having estimated properties of the candidate star clusters, we explore implications of the results and the role of the candidate clusters with respect to the starburst. Before elaborating on that, we discuss the two main sources of uncertainty in the properties of the candidate star clusters: the free-free fractions at 93 GHz and the age of the burst.

Discussion of uncertainties
Both the free-free fractions at 93 GHz and the age of the burst could be better constrained using future observations, and this would result in more accurate estimates of almost all properties of the cluster candidates.
The first source of uncertainty that we discuss is the estimate of the free-free flux, by which we estimate a free-free fraction to the 93 GHz continuum.
In estimating the free-free fractions, we enforce that only two contributions of emission type, free-free and synchrotron, or, free-free and dust, compose the flux at 93 GHz. If we consider three emission types, we would obtain non-unique solutions in estimating the fractions of emission from the spectral index. Assuming two types is the best that we have assessed with the currently available data. Small dust opacities and a majority of negative slopes measured at 93 GHz independently indicate that dust does not substantially contribute in the majority of sources. Therefore assuming two dominant components to the emission at 93 GHz appears to be valid in the majority of sources. In order to decompose the emission at 93 GHz through SED fitting, additional observations at intermediate frequencies (within 2.3 GHz and 350 GHz) of comparable resolution are needed.
To estimate the free-free fraction, we also assume fixed energy distribution slopes for emission from free-free, synchrotron and dust. However, variations from source to source may be expected. If we re-derive the free-free fractions, letting α d = 2 without changing the free-free and synchrotron indices, we find the median dust fraction of sources remains at f d = 0. This is another indication that the spectral index assumed for dust does not have a considerable impact on the majority of sources. The range of reasonable values to consider for the frequency dependence of non-thermal emission is less constrained, where single-injection indices of α syn = −0.5 to −0.8 can be expected all the way up to the dramatic exponential cutoff 2 in a single-injection scenario where the highest energy electrons (at the high frequency end) completely depopulate after a characteristic energy-loss time (e.g., Klein et al. 2018). As we note in Section 3.3, our assumed index of α syn − 1.5 agrees well with the measured value of −1.4 at lower resolution (Bendo et al. 2016). If we re-derive the free-free fractions assuming the canonical value of α syn = −0.8 determined at 10 GHz (Niklas et al. 1997), the median free-free fraction of sources would be f ff = 0.28 and the synchrotron fraction would be f syn = 0.72.
On average the fractional contributions to the 93 GHz emission that we obtain assuming indices of α ff = −0.12, α syn = −1.5, and α d = 4 appear to be reasonable and consistent. Our median free-free fraction of f ff = 0.62 with median absolute deviation of 0.29 agrees within error to the free-free fraction derived by (Bendo et al. 2016) of f ff = 0.84 ± 0.10 at 86 GHz. Our median free-free fraction is also consistent with that derived for NGC 253, f ff = 0.70 ± 0.10, at our same frequency but averaged over 30 (Bendo et al. 2015). Similarly, pilot survey results of the MUSTANG Galactic Plane Survey at 90 GHz and with parsec-scale resolution indicate that > 80% of the emission in (candidate) clusters is composed of synchotron or free-free emission (Ginsburg et al. 2020). We can also perform another self-consistent check on the free-free emission at 93 GHz, by using the recombination line emission and assuming the candidate star clusters have uniform temperatures. A constant temperature implies a constant integrated recombination line to continuum ratio (see Equation 4). Letting R LC ≈ 35 km s −1 (for a temperature of T e = 6000 K), we plug in our measured line and continuum fluxes and solve for f ff . A median free-free fraction of f ff = 0.85 with median absolute deviation of 0.33 is found when considering all sources with detected recombination line emission. In comparison, the relation assumed for the in-band spectral index, also at 0.2 resolution, results in a median value of f ff = 0.66 and a median absolute deviation of 0.24. These values are in reasonable agreement given the uncertainties in the two methods.
Overall, the assumptions made in decomposing the emission at 93 GHz, especially regarding synchrotron 2 where S ∝ e −ν/ν b and ν b is the break frequency emissions, likely affect individual clusters at the ∼ 30% level. This is not enough to bias our overall results, but follow up observations at other frequencies would be extremely helpful.
The second potentially major source of uncertainty is the assumed age of the burst. We discuss in Section 5.2 how we arrive at an adopted age of the clusters of ∼5 Myr. The assumed age has a large impact on the stellar mass inferred from the ionizing photon rate. The ionizing photon output changes substantially (by a factor of 40 from a zero age main sequence to an age of 5 Myr) as the most massive stars explode as supernova. An uncertainty of ∼1 Myr about an age of 5 Myr of a star cluster results in an uncertainty in the inferred stellar mass by a factor of four. This is roughly included in the 0.4 dex uncertainty, though it would represent a systematic offset.

Super Star Clusters
The estimated properties (mass, size, age; Table 6) that we derive for these candidate star clusters meet the criteria for young, massive clusters (e.g., Portegies Zwart et al. 2010), and these sources can be considered super star clusters (SSCs). Stars clusters forming in high gas surface density environments may be able to acquire significant amounts of mass before feedback effects (likely radiation pressure) can disrupt and/or disperse the cluster (e.g., Adamo & Bastian 2016). Super star clusters with stellar masses of M 10 5 M typically have high star formation efficiencies (ε > 0.5) and therefore remain bound.
With the mass surface densities that we estimate and the age of the burst equaling many multiples of the free-fall times, these star clusters will likely remain bound, at least initially. They are being born into a violent environment and clusters often still experience significant mortality after forming. Estimates of the virial mass, escape velocity and momentum driving determined through molecular line observations will help quantity their mass and initial dynamical state.

Cluster Mass Function
In Figure 11, we plot the cluster stellar-mass function of candidate SSCs in NGC 4945 and compare it with the cluster mass distributions in additional galaxies with SSC populations.
A power-law fit to the cluster mass distribution, of 22 sources down to M = 2.4×10 5 M , in NGC 4945 results in a slope of β = −1.8±0.4. Note that the fit to the data results in β = −1.76 ± 0.07, but given the uncertainty in our mass estimates (0.4 dex), we adopt the more conservative figure. We include the SSCs of NGC 253, with stellar mass properties determined on similar spatial scales (∼2 pc resolution) and through H40α recombination lines for a zero-age main sequence population (Leroy et al. 2018;Mills et al. 2020). A power-law fit to the cluster mass distribution of NGC 253, including 12 star clusters down to M = 7.9 × 10 4 M , results in a slope of β = −1.6 ± 0.3. We estimated the completeness limits by eye. The turn-off from power-law distributions likely includes non-physical effects, resulting from the depth/sensitivity of the observations as well as source confusion due to the high inclination in which we view the starburst regions (and which appears to be higher in NGC 253).
We also include recent results from a homogeneous analysis by Mok et al. (2020) of the cluster mass distributions of young (τ < 10 Myr), massive clusters in six galaxies. As representative examples, we include the best fit power-laws from three systems: the Large Magellanic Clouds (LMC), M 51, and the Antennae System. A main result from Mok et al. (2020) is that the cluster mass functions across the six galaxies are consistent with slopes of β = −2.0 ± 0.3. Our measurements for NGC 4945 and NGC 253 are consistent with these findings.
As we discuss in Section 6.6, we estimate the total stellar mass in the burst to be 8.5 × 10 7 M . Extending the best fit of our cluster mass function down to a cluster mass of ∼ 2×10 4 M would account for the additional mass and correspondingly the additional ionizing photon luminosity. As we discuss in Section 6.4, if these lower mass clusters are present they would need to be located in a more-extended region than the super star clusters. Therefore the cluster mass distribution may not reach down to 2 × 10 4 M in the region where the SSCs are located, but perhaps ∼ 4×10 4 M if a third of the additional recombination line emission results from diffuse ionized gas.

Ionizing Photons and Diffuse Ionized Gas
Total ionizing photon rate of the starburst. We measure an H40α recombination line flux of (6.4±1.0) Jy km s −1 integrated within the T2 aperture (see Figure 8) in the intermediate-configuration (0.7 ) data. Assuming the temperature of T e = (6000 ± 400) K from Section 5.3, the total ionized content, which accounts for the full star-formation rate of the starburst, yields an ionizing photon rate of Q T2 = (3.9 ± 0.3) × 10 53 s −1 .
Ionizing photon rate from candidate super star clusters. The sum of the ionizing photon rate over all massive star clusters is ≈ 5.3 × 10 52 s −1 . The total integrated line flux of H40α in the extended-configuration (0.2 ) data is S dv = (2.1 ± 0.6) Jy km s −1 within an aperture (T1) which covers the clusters in the starburst region. With the temperature of T e = (6000 ± 400) K, the integrated line flux corresponds to an ionizing photon rate of Q T1 = (1.2 ± 0.4) × 10 53 s −1 . Given the large uncertainties derived for individual clusters, we consider this value consistent with the sum over individual clusters. Therefore, we conclude that 20% to 44% of the total ionizing photons in the starburst can be attributed to the candidate super star clusters identified in NGC 4945.
Low mass clusters. The intermediate-configuration (0.7 ) data are sensitive to emission on larger physical scales, but they also reach deeper sensitivities per unit area. Consequently, emission from a distribution of many compact, low mass clusters could be traced in these data, but not in the extended-configuration data.
We do not find it likely that the deficit of line emission in the extended-configuration data can be attributed primarily to low mass clusters. Focusing only on the low resolution (0.7 ) observations, the integrated line flux is greater in the larger aperture (T2) than in the aperture covering only the star clusters (T1), indicating that additional line emission originates outside of the area where massive clusters are forming. In order for low mass clusters to account for the difference, these would need to form in a more extended region than the bright SSCs that we see.
AGN and the circumnuclear disk. The circumnuclear disk, including the AGN, could be an origin for ionizing photons observed only on larger scales. The strongest recombination line emission is seen at lowresolution in this region yet no significant line detections are obtained towards Sources 17, 18, and 20. We test this scenario by extracting spectra from the 0.7 intermediate-configuration data and the 0.2 extendedconfiguration data in 7 non-overlapping apertures of 3 diameter consecutively distributed along the major axis of the starburst region. The ratios of the line flux in the intermediate-to extended-configuration data are all consistent within error; the ratios in the seven regions have a mean and standard deviation of 2.3 ± 0.5. This indicates that a deficit of line emission on the physical scales probed in the high resolution data is ubiquitous (and not unique to the region surrounding the AGN), and at most ∼ 20% may originate from the circumnuclear disk and AGN.
Diffuse ionized gas on large-scales. Summarizing the information above, we find that up to ∼70% of the total radio recombination line flux may originate on scales larger than ∼100 pc; this emission cannot be directly attributed to the AGN. Ionizing radiation may be escaping the immediate surroundings of massive stars and reaching larger scales. The approximate age of the burst indi-cates clusters have had time to shed their natal material, dust is not a significant contribution of their emission, low gas mass fractions have been determined, and the current free-fall times of the clusters indicate abundant time to expel gas. Despite these indicators, extinctions may still be considerable and patchy regions may help to leak ionizing photons. Note that if radiation is escaping from star clusters, this would also imply that the stellar masses (derived from the ionizing photon rates) are underestimated.
However, a missing 70% of line flux in the 0.2 resolution data should be considered as an upper limit. Given the low signal-to-noise ratio of these lines, artifacts in both the image and spectral domains can impact the properties of the line profiles. For example, incomplete uv coverage on the scales of the emission can result in incorrect deconvolution and thus systematic underestimates. Deeper observations at high resolution (∼0.1 ) would allow radio recombination lines to be mapped out with adequate signal-to-noise ratio and obtain sensitivities approaching that of the low-resolution (0.7 ) observations.

Role of the AGN
Although the Seyfert 2 AGN in NGC 4945 is one of the brightest in our sky at X-ray energies, obscuration by e.g., A V ≥ 160 (Spoon et al. 2000) has prevented its observation at virtually all other wavelengths (also given the spatial resolutions observed). Our observations at 93 GHz provide another piece of evidence for its existence.
We detect a point source in the 93 GHz image with a deconvolved size determined by PyBDSF of FWHM = 0.84 ± 0.01 pc. The flux of 9.7 ± 1.0 mJy that we extract from the source is ∼ 10% of the total continuum (∼120 mJy) flux in the 0.12 resolution data. The Band 3 spectral index we derive is α 93 = −0.85 ± 0.05, which is consistent with freshly accelerated electrons emitting synchrotron radiation. However, we follow the procedure outlined in Section 3.3, assuming free-free emission may also contribute, and place a limit on the ionizing photon rate. The spectral index at 93 GHz corresponds to f ff ≥ 0.47 ± 0.03. If we assume a temperature of T e = 6000 (30 000) K this would correspond to Q 0 = 1.1 (0.5) × 10 52 s −1 . The low level of escaping ionizing radiation is consistent with previous work from MIR line ratios (Spoon et al. 2000).
We place a limit on the ionizing photons that could be leaking into medium. The bolometric luminosity estimated from the X-ray luminosity is L bol ≈ 6 × 10 43 erg s −1 , which we determined using the relation for Seyfert AGN of L bol ∼ 20L X (Hopkins et al. 2007).
Given the bolometric luminosity, we can now estimate the total expected ionizing photon luminosity. We use the standard relation of Elvis et al. (1994), which posits L bol /L ion ∼ 3 with an average ionizing photon energy of 113 eV. Therefore the expected ionizing photon luminosity of the AGN is L ion ≈ 2 × 10 43 erg s −1 ≈ 5 × 10 9 L , while our limit corresponds to ≈ 2 × 10 42 erg s −1 . This would imply that < 10% of the ionizing photon luminosity of the AGN is escaping into the surrounding nuclear starburst and creating ionized gas.

Total Burst of Star-formation
We can convert the total ionizing photon rate of the burst Q T2 (see Section 6.4), as measured in the T2 region, to a luminosity. Assuming an average ionizing photon energy of hν ≈ 17 eV, the ionizing luminosity is L T2 ≈ 1.1 × 10 43 erg s −1 = 2.8 × 10 9 L . A ratio of bolometric luminosity to ionizing photon luminosity L bol /L 0 ∼ 14 is predicted for a 5 Myr old population, using Starburst99 (Leitherer et al. 1999). Thus, we can expect a bolometric luminosity of 4 × 10 10 L which agrees within error to the bolometric luminosity derived from FIR observations, L bol = (2.0 ± 0.2) × 10 10 L (Bendo et al. 2016).
We note that, if we had adopted a lifetime for the burst of 4 Myr rather than 5 Myr, Starburst99 calculations expect a bolometric luminosity of L bol ≈ 2 × 10 10 L . Also, this calculation assumes the source of the additional recombination line emission in the intermediateconfiguration data is also characterized by a 5 Myr old stellar population.
The ratio of bolometric luminosity to mass of a 5 Myr old burst is estimated at Ψ ∼ 470 L /M . With an L bol = 4 × 10 10 L , the inferred total stellar mass of the population is ∼ 8.5 × 10 7 M . The total stellar mass in our candidate clusters is ≈ 1.1 × 10 7 M , or ∼13% of the expected total stellar mass.

Star Clusters and the Central Wind
An outflow of warm ionized gas in NGC 4945 has been modeled as a biconical outflow with a deprojected velocity of v ≈ 525 km s −1 with emission out to 1.8 kpc (Heckman et al. 1990). The outflow has been traced in [NII], [SII] and Hα in multiple analyses (Heckman et al. 1990;Moorwood et al. 1996;Mingozzi et al. 2019). Heckman et al. (1990) estimated a total energy of 2 × 10 55 erg and momentum flux of 9 × 10 33 dyne = 1400 M km s −1 yr −1 (after correcting for the updated distance). These values do not include possible contributions by colder phases, which are likely associated with the NGC 4945 outflow (Bolatto et al., in prep.). While these values are uncertain, they are use-ful for an order of magnitude comparison with expected properties of the star clusters.
Gas reaching 1.8 kpc and moving at 525 km s −1 would have been launched ∼3 Myr ago assuming a constant velocity, which is within the time-frame of the burst. A ∼3 Myr time-scale and a velocity of 525 km s −1 would put 8 × 10 6 M of warm ionized gas mass into the outflow. Using Starburst99 (Leitherer et al. 1999), an estimated mechanical luminosity of ∼ 3 × 10 41 erg s −1 from the clusters, fairly constant over the 5 Myr age, equates to an injected energy of ≈ 6 × 10 55 erg. It would require 30% of the expected mechanical energy output of the clusters to drive the ionized outflow. From simulations, the total momentum supplied to the ISM per supernova is expected to be 2.8 × 10 5 M km s −1 (Kim & Ostriker 2015, and references therein), with stellar winds contributing an additional 50% to the momentum. Star-burst99 predicts the supernova rate to be fairly constant at ∼0.009 yr −1 for these clusters. From these values, we estimate that the momentum supplied to the ISM by the candidate star clusters would be ∼ 3800 M km s −1 yr −1 . It would require about 40% of the expected momentum output of the clusters to drive the ionized outflow. Therefore, neither the energetics nor the inferred momentum of the outflow prevent it from being driven solely by the star clusters in NGC 4945, although there are very considerable uncertainties associated with this calculation.

Comparison with NGC 253
NGC 253 is a nearly edge-on galaxy located nearby (3.5 ± 0.2 Mpc ; Rekola et al. 2005) with similar properties as NGC 4945. It hosts a central starburst spanning ∼200 pc with young, massive clusters (Leroy et al. 2018;Mills et al. 2020). A major difference between these two galaxies is that NGC 4945 shows unambiguous signatures of harboring an active super-massive black hole. Since NGC 4945 is just the second analysis we have undertaken with millimeter wavelength ALMA observations on scales which resolve the clusters, it is relevant to compare the properties of their super stars clusters.
Analyses of the cluster population of NGC 253 at ∼2 pc resolution with ALMA (Leroy et al. 2018;Mills et al. 2020) revealed 18 sources as (proto-)super star clusters which are in the process of forming or close to a zero age main sequence (∼1 Myr) (see also Rico-Villas et al. 2020). Overall the properties of the clusters are strikingly similar as those in NGC 4945; they have a median stellar mass of 1.4 × 10 5 M and FWHM sizes of 2.5-4 pc. The star clusters of NGC 253 boast slightly higher ionizing photon luminosities and larger gas fractions M gas /M tot ∼ 0.5, reflecting a slightly younger age.
The slopes of the cluster mass function we derive are also consistent (see Section 6.3), just offset by a factor of ∼3.5 in mass. As in NGC 4945, at least 30% of the of the nuclear starburst of NGC 253 originates in a clustered mode of star-formation.
While we have not estimated some of the dynamical properties of NGC 4945, the presence of broad recombination line emission in four sources in NGC 253 indicates the star clusters are operating under similar processes. In NGC 253, the sources are young enough that feedback has not managed to unbind a large fraction of the gas from the clusters. The slightly higher total mass surface densities and smaller free-fall times in NGC 4945 indicate its clusters might be surviving a young, disruptive stage.
The concerted feedback of the young (< 10 Myr), forming SSCs identified in NGC 253 (Leroy et al. 2018) are likely not responsible for the starburst driven outflow, as they are expected to impart momentum that is a factor of 10-100 lower than the outflow momentum estimated in CO (Bolatto et al. 2013;Krieger et al. 2019). Some (global) event may be initiating the active formation of star clusters. On the other hand, the star clusters in NGC 4945 could potentially influence the outflow of warm ionized gas; some event appears to be inhibiting the formation of new star clusters.

SUMMARY
Massive, clustered star-formation is an efficient and possibly common mode of star-formation in high gas density environments. Nearby galaxies with bursts of star-formation in the central O(100 pc) are local laboratories to study this mode, and for NGC 4945, in the presence of a Seyfert AGN. High levels of dust extinction (A V ∼ 40) in the nearly edge-on (i ∼ 72 • ) central region of NGC 4945 have previously prevented the direct observation and characterization of its massive star clusters.
We identify 27 super star cluster candidates in the central starburst of NGC 4945. We derive properties of the candidate clusters through ALMA observations at 2.2 pc resolution of the 93 GHz (3 mm) free-free emission and hydrogen recombination line emission (H40α and H42α) arising in photo-ionized gas. We also use and present high-resolution 350 GHz imaging of the dust continuum observed with ALMA, and supplement our analysis with 2.3 GHz continuum imaging which primarily traces synchrotron emission (Lenc & Tingay 2009).
Our results are as follows: • The 27 point sources identified (≥ 10σ) in 93 GHz continuum emission as candidate super star clusters have FWHM sizes of 1.4-4.0 pc. The 93 GHz emission in these bright, compact regions is dominated by free-free emission, with a median freefree fraction of f ff = 0.62. Synchrotron emission from recent supernova remnants contributes to the 93 GHz emission with a median fraction of f syn = 0.36. Substantial dust emission is found in three sources.
• We average the spectra of the H40α and H42α recombination lines to synthesize an effective H41α profile. Recombination line emission is detected in 15 candidate clusters, generally with narrow (FWHM ∼ 36 km s −1 ) line widths. Six of the detected sources have significant (f syn 0.5) synchrotron emission; three of those have broad line widths with FWHM > 105 km s −1 .
• We estimate an electron temperature of T e = (6000 ± 400) K of the ionized gas using the flux ratio of the integrated line to free-free continuum. This electron temperature implies an average metallicity of 12 + log 10 (O/H) = 8.9 ± 0.1 surrounding these young massive stars. The ionized gas densities of log 10 (n e /cm −3 ) = 3.1-3.9, that we derive are typical of classic H II regions. The ionized gas masses of the clusters are typically < 1% of the estimated stellar mass.
• We determine ionizing photon rates of the candidate SSCs in the range log 10 (Q 0 /s −1 ) ∼ 50.4 -51.8. Adopting an age of ∼5 Myr, the stellar masses implied by the ionizing photon rates are log 10 (M /M ) ∼ 4.7-6.1. The sum of the stellar masses of the candidate SSCs is ≈ 1×10 7 M . The uncertainties on these measurements are 0.4 dex. We discuss the age estimate in Section 5.2.
• We fit the cluster stellar-mass distribution and find a slope of β = −1.8 ± 0.4. The slope of the fit is consistent with our fit to the candidate SSCs in the central starburst of NGC 253 (Mills et al. 2020) and to recent findings in the LMC, M 51, and the Antennae system (Mok et al. 2020).
• We estimate the gas mass of the candidate clusters from dust emission. The total mass, M Tot , is evaluated by the combined stellar and gas masses. Gas mass fractions range from M gas /M Tot = 0.04-0.22, with a median value of 0.13. We calculate the total-mass surface density of the clusters and find a median value of Σ Tot = 3 × 10 4 M pc −2 . The median free-fall timescale is 0.04 Myr.
• With low-resolution (0.7 ) observations of the H40α recombination line, we measure a total ion-izing photon rate of the burst of Q 0 = (3.9±0.3)× 10 53 s −1 . The candidate star clusters that we analyze contribute 20-44% of the total ionizing photon rate. Additional recombination line emission present in the low-resolution data appears to be ubiquitous throughout the starburst region, cannot be directly attributed to the AGN, and is also found outside of the area where the SSCs are located. Diffuse ionized gas may be responsible for some of the additional emission, although low mass clusters and distributed star formation are expected to also contribute. This indicates that ionizing radiation may be escaping the immediate surroundings of massive stars and reaching the larger (> 100 pc) scales traced at lower resolution.
• We compare the candidate SSCs in NGC 4945 with the (proto-)SSCs recently identified in the central starburst of NGC 253 (Leroy et al. 2018;Mills et al. 2020). The age of the burst (1 Myr) is a bit younger in NGC 253, the gas mass fractions (∼0.5) are a bit higher, and the stellar mass of the clusters are slightly smaller (factor of 3.5). While the actively forming clusters are not a major contributor to the starburst driven outflow in NGC 253, the slightly older population of the star clusters in NGC 4945 may contribute to driving a nuclear outflow of warm ionized gas.
• Strong, variable X-ray emission, which is Compton thick, provides evidence for a Seyfert AGN in NGC 4945. The bright, 93 GHz source which we presume to be the AGN is point-like in our 2.2 pc beam. We measure its in-band spectral index at 93 GHz to be α 93 = −0.85 ± 0.05, likely dominated by synchrotron emission. We do not detect recombination line emission from this point source. Our observations support previous findings in which UV ionizing radiation from the AGN is heavily obscured in all directions. We estimate an upper limit for its escaping ionizing photon rate of Q 0 < 1 × 10 52 s −1 , which is 10% of the expected luminosity of ionizing photons for a typical Seyfert AGN.
• Lastly, we report on a shortcoming of H42α as a "low-resolution" tracer of ionizing radiation. When observed with low spatial resolution and/or from broad line components, we find this spectral line is likely contaminated by line emission from other species, resulting in the recombination line flux (and star-formation rates derived from it) being overestimated by a factor of ∼2, when compared to a similar analysis of H40α. where y = n He + /n p is the ratio of singly ionized helium to hydrogen by number and n He + is the singly ionized helium number density. If we assume the emission region is composed of only hydrogen and singly ionized helium, then n p n + = n p n p + n He + = 1 + n He + n p −1 = (1 + y) −1 .
We can rearrange the integrated RRL line to continuum ratio and solve for the electron temperature T e of the emission region, T e = 10 4 K b n+1 (1 + y) The rate of ionizing photons (E > 13.6eV ) is given by, where α B is the case B recombination coefficient (Draine 2011), α B = 2.59 × 10 −13 cm 3 s −1 T e 10 4 K −0.833−0.034 ln (Te / 10 4 K) (A13) which is valid for 3000 K < T e < 30,000 K. Thus we have, Q 0 = 3.805 × 10 51 s −1 n e n + V 5 × 10 8 cm −6 pc 3 T e 10 4 K −0.833−0.034 ln (Te / 10 4 K) . (A14) B. RECOMBINATION LINE SPECTRA OF ALL SOURCES In Figure 12, we show that the central velocities of our detected recombination lines are in good agreement with the kinematic velocity expected of the disk rotation, though we include the recombination line spectra for all 29 sources (including those that are not significantly detected). The recombination line spectra from our sources are overlaid on H40α spectra extracted from the intermediate configuration observations (0.7 resolution). High-resolution recombination line emission is coincident with emission from the intermediate configuration data. Velocity [km s 1 ] Figure 12. Effective H41α recombination line spectra (thick, black line) for sources with significantly detected emission and the best fit line profile (red) -same as in Figure 6. H40α line emission tracing the kinematics is overlaid in thin purple; it represents the H40α spectrum centered at the same source locations and extracted from low resolution (∼0.7 ) observations.  Figure 13. SEDs constructed for each source, same as in Figure 5. The dashed green line represents a dust spectral index of α = 4.0, normalized to the flux density we extract at 350 GHz (orange data point). The dashed black line represents a free-free spectral index of α = −0.12, normalized to the flux density we extract at 93 GHz (black data point). The white data points show the flux densities extracted from the band 3 spectral windows. The gray shaded region is the 1σ error range of the band 3 spectral index fit, except we have extended the fit in frequency for displaying purposes. The purple dashed line represents a synchrotron spectral index of α = −1.5, normalized to the flux density we extract at 2.3 GHz (purple data point); except for Source 14 where the solid purple line represents the matched 2.3-23 GHz fit. Error bars on the flux density data points are 3σ.