# Clustered Formation of Massive Stars within an Ionized Rotating Disk

Roberto Galván-Madrid Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia, Michoacán 58089, México.  Qizhou Zhang Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA.  Andrés Izquierdo European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany.  Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands.  Charles J. Law Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA.  Thomas Peters  Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany.  Eric Keto  Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA.  Hauyu Baobab Liu Physics Department, National Sun Yat-Sen University, No. 70, Lien-Hai Road, Kaosiung City 80424, Taiwan, R.O.C.  Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, ROC.  Paul T. P. Ho Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, ROC.  East Asian Observatory (EAO), 660 N. Aohoku Place, University Park, Hilo, Hawaii 96720, USA.  Adam Ginsburg Department of Astronomy, University of Florida, P.O. Box 112055, Gainesville, FL, USA.  Carlos Carrasco-González Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia, Michoacán 58089, México.

(Accepted to ApJ Letters, December 7, 2022)

###### Abstract

We present ALMA observations with a 800 au resolution and radiative-transfer
modelling of the inner part ($r\approx 6000$ au) of the ionized accretion flow
around a compact star cluster in formation at the center of the luminous
ultra-compact (UC) Hii region G10.6-0.4. We modeled the flow with an ionized
Keplerian disk with and without radial motions in its outer part, or with an
external Ulrich envelope. The MCMC fits to the data give total stellar masses
$M_{\star}$ from 120 to $200~{}M_{\odot}$, with much smaller ionized-gas
masses $M_{\mathrm{ion-gas}}=0.2$ to $0.25~{}M_{\odot}$. The stellar mass is
distributed within the gravitational radius $R_{g}\approx 1000$ to 1500 au,
where the ionized gas is bound. The viewing inclination angle from the face-on
orientation is $i=49$ to $56~{}\deg$. Radial motions at radii $r>R_{g}$
converge to $v_{r,0}\approx 8.7$ km s-1, or about the speed of sound of
ionized gas, indicating that this gas is marginally unbound at most. From
additional constraints on the ionizing-photon rate and far-IR luminosity of
the region, we conclude that the stellar cluster consists of a few massive
stars with $M_{\mathrm{star}}=32$ to $60~{}M_{\odot}$, or one star in this
range of masses accompanied by a population of lower-mass stars. Any active
accretion of ionized gas onto the massive (proto)stars is residual. The
inferred cluster density is very large, comparable to that reported at similar
scales in the Galactic Center. Stellar interactions are likely to occur within
the next Myr.

H II regions (694) – Massive stars (732) – Star formation (1569) – Young star
clusters (1833) – Radiative transfer simulations (1967)

††facilities: ALMA, SMA††software:  CASA [36], Astropy [1], sf3dmodels [16],
RADMC-3D [3], emcee [5], corner [4], spectral-cube (https://spectral-
cube.readthedocs.io), pvextractor (https://pvextractor.readthedocs.io), IMF
(https://github.com/keflavich/imf)

##  1 Introduction

Stars with masses greater than 30 $M_{\odot}$ are rare among the stellar
population in galaxies, but they play a main role in shaping the dynamical
evolution of the interstellar medium through feedback from their radiation,
ionization, and supernovae [for a review of feedback effects in star
formation, see 27]. Since these very massive stars are mostly formed in dense
stellar clusters, their violent death may lead to the formation of black holes
or neutron-star binary systems that give rise to gravitational-wave events
[e.g., 8]. Massive protostars must accrete more than 30 $M_{\odot}$ to become
early O-type stars, but they begin core hydrogen burning well before reaching
their final masses [40, 14]. This poses a serious challenge to the current
paradigm that massive stars form through accretion of molecular gas. The
detection of resolved molecular disks around protostars more massive than
$20~{}M_{\odot}$ shows that this scenario is valid for the formation of
individual stars up to about a few tens of solar masses [e.g., 18, 15, 47, 39,
33]; however, it likely breaks down for the case of the formation of more
massive stars, especially in clustered environments [9]. In a clustered
scenario, and after the protostars have reached several tens of solar masses,
the amount of EUV photons is large enough to ionize the otherwise molecular
accretion flow(s) and the clump in which they reside [20, 21, 41, 43, 25].
This raises several outstanding questions: How does star formation proceed in
such environment? Is it through an ionized accretion disk? Does this disk
surround a single or a multiple stellar system? When does active
(proto)stellar accretion stop due to the increasing photoionizing feedback?

We present the first spatially resolved observations of an accretion flow in
ionized gas around a dense star cluster in formation. Our observations were
performed with the Atacama Large Millimeter/submillimeter Array (ALMA) in the
1.3 mm continuuum and hydrogen $30\alpha$ recombination line emission. The
target is G10.6-0.4 (hereafter G10.6), a luminous ($L_{\mathrm{FIR}}\approx
3\times 10^{6}~{}L_{\odot}$) and spatially concentrated [see 29] massive star
formation region located at 4.95 kpc from the Sun [46]. The molecular gas in
the central pc-scale clump of G10.6 presents infall and rotation motions
toward an inner, flattened ultracompact Hii region [e.g., 11, 24, 49, 32, 2].
A characteristic “bullseye” velocity pattern in the most redshifted molecular
absorption against the central UC Hii region reveals that the molecular infall
is coupled with rotation [see 49]. By equating the observed molecular infall
speed to free fall, Sollins et al. [49] estimated a central stellar mass of
$M_{\star}\sim 150$. Keto & Wood [22] reported a velocity gradient in the
H$66\alpha$ recombination line within the central 0.05 pc. Those authors
proposed that the central ionized gas is an accretion flow feeding a nascent
cluster of massive stars.

##  2 Observations

Observations of G10.6 with ALMA (project ID: 2015.1.00106.S, PI Q. Zhang) were
carried out in four execution blocks (EBs) between 2016 and 2017. The total
on-source integration time was 2.2 hours. The two sessions executed on 2016
September 9 and 10 had 36 12-m antennas in the array with physical baselines
from 15 m to 3144 m. The remaining sessions were observed in 2017 July 19 had
42 12-m antennas with baselines from 19 m to 3697 m. The precipitable water
vapor (PWV) in the atmosphere ranged between 0.46 to 0.50 mm during the
observations. System temperatures varied from 49 to 75 K. Quasars J1832-2039,
J1924-2914 and J1733-1304 were used as gain, bandpass, and flux calibrators,
respectively.

The digital correlator was configured to the FDM mode for four spectral
windows, centered at 217.9 GHz, 220.0 GHz, 232.0 GHz and 233.9 GHz,
respectively. Each spectral window has an effective bandwidth of 1.875 GHz,
divided into 1910 channels, providing a uniform spacing of 0.976 MHz per
spectral channel. The pointing center of G10.6 was $\alpha$(J2000) = 18h 10m
28.7s, $\delta$(J2000) = $-19^{\circ}55^{\prime}49.1^{\prime\prime}$.

The calibration of the visibilities was performed in Common Astronomy Software
Applications package [CASA 36] using the pipeline script supplied by the ALMA
observatory. The calibrated visibilities were Fourier transformed and
‘cleaned’ using tclean and Briggs weighting of visibilities with robust
parameter of 0.5. Spectral channels free of significant line emission are used
to construct continuum data. Continuum emission was subtracted from the
calibrated visibilities to produce the spectral line visibilities. The rms
noise in the continuum image is about 35 $\mu$Jy beam-1. The rms noise in the
line free channels of the H30$\alpha$ image cube is 0.5 mJy beam-1 per 1.5 km
s-1 channel, with a beam of $0^{\prime\prime}.15\times 0^{\prime\prime}.13$,
$\mathrm{P.A.}=-65.8\deg$.

##  3 Results

###  3.1 Overview

Figure 1: Distribution of molecular and ionized gas in the dense cluster
G10.6-0.4. (a): CH3OH line emission imaged with the SMA. Red, green, and blue
show CH3OH with upper energies 35, 59, and 97 K [see 32]. The white square
marks the region of the ALMA maps shown in the other panels. (b): Emission of
the CH3CN J=12-11 line [see 28] [Jy beam-1 km s-1]. (c): 1.3-mm continuum [Jy
beam-1]. (d): Integrated H30$\alpha$ line emission [Jy beam-1 km s-1]. The
spatial cut for the P-V diagrams is shown with a red line.

The ALMA observations with a physical resolution of 792 au ($0.16\arcsec$ at a
distance of 4950 pc) resolve the 1.3 mm continuum and H${30}\alpha$ emission
within the inner 6000 au radius of the star-forming clump. Figure 1 presents
an overview. The top-left panel shows CH3OH emission imaged with the
Submillimeter Array. The molecular envelope of a $\sim 0.3$ pc radius flattens
toward its center, most notably in the warmer (blue) CH3OH transition [see
32]. An evacuated bipolar cavity is also seen. The emission from complex
organic molecules (COMs) in the inner 0.05 pc was recently studied by Law et
al. [28] using ALMA. Those authors found highly structured emission in the
form of localized hot cores accompanied by bright, extended emission. The
overall COM emission presents a flattened X-shaped morphology with a
brightness dip toward the H${30}\alpha$ peak (see Fig. 1, CH3CN panel). In
contrast, the 1.3-mm continuum peaks at the central dip of the COM emission
and matches well the H${30}\alpha$ emission (Fig. 1, bottom panels). The
continuum emission is almost entirely due to the free-free radiation in the
region that we analyze. The average continuum intensity within a stripe of
$1\arcsec$ length and 1 beamwidth height is 0.045 Jy beam-1 (40 K), whereas
the free-free continuum expected from the H30$\alpha$ intensity in the same
area, using equation (14.29) in Rohlfs & Wilson [45], is 0.044 Jy beam-1.
Contributions from warm dust to the continuum become important at the location
of hot cores in the periphery of the central ionized structure [30, 28]. This
flattened ionized emission within 0.05 pc radius has been interpreted as the
innermost part of a cluster accretion flow which transitions from being
molecular to ionized at its center [e.g., 24, 22].

The LSR systemic velocity of the H$30\alpha$ line is determined to be
$V_{\mathrm{sys}}\approx 0.5$ km s-1. The H$30\alpha$ panel of Fig. 1 shows
the directional cut at $\mathrm{PA}=307.8^{\circ}$ (East of North) used for
the modelling. Figure 2 shows the resulting position-velocity (P-V) diagram.
High-velocity emission characteristic of Keplerian rotation is seen on either
side from the position center, reaching velocities up to $\pm\sim 40$ km s-1
from $V_{\mathrm{sys}}$ in the inner few $10^{3}$ au. This kinematic structure
is in direct contrast to that of classical H ii regions, in which the
photoionized gas at $10^{4}$ K exerts an outward pressure that drives an
unimpeded expansion. We rather interpret our observations in the context of H
ii regions still dominated by the gravity of the embedded stars [19, 20, 21].
In these models the ionized gas within the gravitational radius
$R_{g}=GM_{*}/c_{s}$ (where $c_{s}$ is the sound speed of the ionized gas)
remains bound, and can continue to accrete onto the stars.

###  3.2 Radiative Transfer and Model Fitting

We use RADMC-3D [3] to model the 1.3-mm free-free continuum and the
H$30\alpha$ emission in non-LTE. We use the version that allows to calculate
recombination lines as presented in Peters et al. [42]. The physical 3D model
grids are created with sf3dmodels [16] and manipulated with the tools within
that package to be compatible with RADMC-3D.

The observational data sets that we model are a 1D intensity cut of the 1.3-mm
continuum and the 2D P-V diagram of the H${30}\alpha$ emission. Both cuts are
centered at ICRS RA $=$ 18h 10m 28.652s, Dec $=$ $-19$d 55m 49.66s,
$\mathrm{PA}=307.8^{\circ}$. The length of the P-V and continuum cuts are
2.32″(11480 au) and 0.72″(3560 au), respectively, and their width is 0.16″.
The continuum cut was truncated because the emission deviates from a single
power law beyond the defined length. This is probably due to a real change in
the density profile of the ionized disk, or due to contributions from dust
emission at larger radii (see Section 3.1). The continuum intensity cut is
only used in the model fitting to constrain the distribution of the electron
density $n_{e}$. We consider fully ionized hydrogen gas ($n_{e}=n_{ion}$). The
kinematic structure is as follows.

A Keplerian rotating disk. The simplest model is a flared Keplerian disk as in
Pringle [44]. The disk electron density $n_{e}$ goes as:

$n_{e}(R,z)=n_{0}\biggl{[}\frac{R}{R_{0}}\biggr{]}^{p}\exp[-z^{2}/2H^{2}],$
(1)

where $R$ is the polar radius normalized at $R_{0}=10$ au, $z$ is the distance
from the midplane, and $H(R)\propto R$ is the scale height. For simplicity, we
fix $H$ to have a linear dependency with radius, but we have verified that
steeper radial variations $H(R)\propto R^{1.25}$ or $\propto R^{1.5}$ only
affect the results minimally. The normalization of $H$ is selected such that
at the gravitational radius $H(R_{g})=R_{g}$, in accordance with the physical
expectation that within $R_{g}$ the ionized gas is bound [13, 21, 50]. These
assumptions leave only two free parameters for the density description in
Equation 1: its normalization $n_{0}=n_{e}(R=10~{}\mathrm{au})$, and its
power-law index $p$.

The circular velocity in the Keplerian model is defined as:

$v(R)=\biggl{[}\frac{GM_{\star}}{R}\biggr{]}^{0.5},$ (2)

where $G$ is the gravitational constant and $M_{\star}$ is the central stellar
mass. This kinematical model has two free parameters: $M_{\star}$ and the
inclination angle $i$, where $i=90\deg$ means an edge-on view.

A Keplerian disk with outer radial motions. The next level of complexity is to
include radial spherical motions to the disk, which are vector-added to the
Keplerian rotating disk. This adds a third kinematical free parameter
$v_{r,0}$.

A Keplerian disk with an Ulrich envelope. We also considered a model
consisting of a Keplerian disk surrounded by an Ulrich-type envelope [52, 37].
This model has been previously used to interpret the kinematics of G10.6 at
$\sim 0.1$ pc scales [22], and is a widely-used option to interpret rotating-
infalling envelopes which settle into a rotationally supported disk [e.g., 23,
16]. By construction, the radial motions in the outer part of this model are
inward. We run models with four free parameters: the stellar mass $M_{\star}$,
the viewing inclination angle $i$, the accretion rate of the Ulrich envelope
$\dot{M}_{\mathrm{env}}$, and an adimensional scaling factor $A$ which
controls the density contrast between the inner disk and the outer envelope.

Table 1: Model Fitting

Density profile | Keplerian disk |  Keplerian $+$ radial |  Keplerian $+$ Ulrich
---|---|---|---
$n_{0}$ | $p$ | $M_{\star}$ | $i$ | $M_{\star}$ | $i$ | $v_{\mathrm{r,0}}$ | $M_{\star}$ | $i$ | $\dot{M}_{\mathrm{env}}$ | $A$
$[10^{12}~{}\mathrm{cm}^{-3}]$ |  | $[\mathrm{M_{\odot}}]$ | $[\deg]$ | $[\mathrm{M_{\odot}}]$ | $[\deg]$ |  [km s-1] | $[\mathrm{M_{\odot}}]$ | $[\deg]$ | $[10^{-5}~{}\mathrm{M_{\odot}~{}yr}^{-1}]$ |
$5.5\pm 0.1$ | $-0.67\pm 0.01$ | $194.0\pm 0.8$ | $51.2\pm 0.1$ | $127.9\pm 0.8$ | $55.5\pm 0.1$ | $8.72\pm 0.04$ | $187.8\pm 0.9$ | $48.9\pm 0.2$ | $3.33\pm 0.02$ | $5.29\pm 0.03$

Note. — Continuum fit for density profile parameters: density normalization
$n_{0}$ and index $p$. H$30\alpha$ fit for kinematical parameters: Keplerian
model – stellar mass $M_{\star}$, viewing inclination angle $i$; Keplerian
model with external radial motions – $M_{\star}$, $i$, and radial velocity
$v_{\mathrm{r,0}}$; Keplerian model with external Ulrich envelope –
$M_{\star}$, $i$, accretion rate of the envelope $\dot{M}_{\mathrm{env}}$,
disk density scaling factor $A$. The statistical errors ($\pm 1\sigma$) from
the MCMC fitting are defined to contain $68.2~{}\%$ of the values.

Table 1 summarizes the free parameters of the different modelling scenarios.
Figures 4, 5, and 6 in the Appendix show diagnostic plots for each of the
best-fit models. The modelling and fitting procedure is as follows:

  * •

Homogenize the continuum and line data to a common circular beam of
$\mathrm{FWHM}=0.16\arcsec$ (792 au), with a pixel size of $0.08\arcsec$.
Nyquist sampling was chosen to avoid overfitting the data.

  * •

Extract the observational continuum intensity and P-V cuts as previously
described.

  * •

Create a 3D model with sf3dmodels [16] and export it to the grid format of
RADMC-3D [3, 42].

  * •

Use RADMC-3D to calculate the radiative transfer of the free-free continuum
and non-LTE H$30\alpha$ emission of the previously computed model.

  * •

Convolve the model images to the same resolution and pixel size of the
observations. Then extract the model continuum and P-V cuts.

  * •

Compare observations and model, iterating over the model parameter space using
the MCMC sampler implemented in emcee [5].

The log-likelihood function that is maximized is defined as:

$\chi^{2}=-0.5\sum\biggl{(}\frac{I_{\mathrm{data}}-I_{\mathrm{model}}}{\sigma_{\mathrm{noise}}}\biggr{)}^{2},$
(3)

where the summation is over the 9 pixels of the continuum cut or the $29\times
60=1740$ pixels of the H${30}\alpha$ P-V image, depending on the respective
maximization111Note that the number of independent measurements is half of the
number of pixels. We have verified that changing the number of pixels per beam
does not affect the results.. To get a conservative estimate of the noise in
the fitted data, for the line modelling we measure the rms in the P-V diagram
$\sigma_{\mathrm{noise}}=1.5$ mJy beam-1. Note that it is larger than the rms
noise measured in a single channel in a region free of emission (see Section
2), since it contains contributions from channel-to-channel variations and
correlated noise due to interferometric sidelobes. For the continuum modelling
we set $\sigma_{\mathrm{noise}}=0.5$ mJy beam-1.

The continuum modelling is used to find the two free parameters of the disk
density profile: $n_{0}=5.5\times 10^{12}$ cm-3, $p=-0.67\approx-2/3$ (see
Table 1). These two density parameters are then fixed in the three different
kinematical models of the H$30\alpha$ line. Figure 2 shows the results for the
continuum and line modelling. The best-fit models give similar residuals
because their overall properties are similar (see figures in the Appendix).
Figure 3 shows the posterior probability “corner” plots for the different
types of models. The statistical errors are very small ($\sim 1~{}\%$) and
should be considered as strict lower limits. Our consideration of different
scenarios shows that the main source of uncertainty is in the model
assumptions.

Figure 2: Results from the radiative-transfer model fitting. Top row, left
panel: continuum modelling. The top cut shows the ALMA data and the bottom cut
shows the best-fit power law model for the electron density. Top row, center
and right panels: Keplerian disk best-fit and residuals. Middle row: ALMA
H$30\alpha$ position-velocity diagram, along with the best fit Keplerian model
with radial motions and residuals. Bottom row: same as middle row for the
Keplerian disk with Ulrich envelope. Residuals are defined as observational
data minus model.

For the case of a purely Keplerian model, the MCMC runs converge to
$[M_{\star},i]=[194.0~{}M_{\odot},51.2~{}\deg]$ . For the scenario of a
Keplerian disk with radial motions, we initially set the radius at which
radial motions start to the gravitational radius of the purely Keplerian
model: $r_{0}=R_{g}(M_{\star}=194~{}M_{\odot})=2320$ au. The resulting central
mass is consistently smaller ($\approx 130~{}M_{\odot}$) when radial motions
are included compared to the purely Keplerian case, therefore we updated the
fixed value of $r_{0}$ to 1550 au, or about $R_{g}$ for a central mass of
$130~{}M_{\odot}$. The resulting values for the free parameters of this fit
are $[M_{\star},i,v_{r,0}]=[127.9~{}M_{\odot},55.5\deg,8.72$ km s-1] (see
Table 1). The residuals from this model are slightly improved because radial
motions help to produce a larger velocity spread at radii beyond $\sim 1000$
au (see Fig. 2). We have verified that this result is insensitive to our
selection of $r_{0}$. Including $r_{0}$ as a fourth free parameter gives
consistent results:
$[M_{\star},i,v_{r,0},r_{0}]=[123.3~{}M_{\odot},56.4~{}\deg,8.0$ km s-1$,972$
au]. The assumption that radial motions start at the gravitational radius
given by our simple prescription might be only approximately valid, since
radiation pressure onto ions could decrease the radius at which outward radial
motions start [51]. Interestingly, the results presented in Table 1 where
$v_{r,0}$ was restricted to be positive (i.e., outward radial motions) are
identical if $v_{r,0}$ is restricted to be negative (inward motions).
Therefore, the modelling cannot distinguish between these two scenarios
because the H30$\alpha$ emission is optically thin. The arguments for the
selection of $r_{0}\sim R_{g}$ suggest an outward interpretation, but the
infalling molecular gas in the exterior suggests otherwise (see Section 4). It
is also remarkable that the fitted radial motions $v_{r,0}$ are identical to
the sound speed of ionized gas ($c_{s}=8.6$ km s-1) at the assumed electron
temperature $T_{e}=9000$ K. This means that, if unbound, the ionized gas in
the disk has not yet accelerated to its terminal supersonic expansion [6].
Finally, the results for the kinematical models of a Keplerian disk with an
external Ulrich envelope are
$[M_{\star},i,\dot{M}_{\mathrm{env}},A]=[187.8~{}M_{\odot},48.9~{}\deg,3.3\times
10^{-5}~{}M_{\odot}$ yr${}^{-1},5.29]$ (see Table 1). The obtained stellar
mass and inclination angle are consistent with the values of the disk-only
models.

Figure 3: Corner plots of posterior probability distributions of the model
parameters fitted by maximizing the log-likelihood $\chi^{2}$ function defined
in Equation 3. The top-left panels show the parameters labeled as “Density
profile” in Table 1. The top-right panels correspond to the “Keplerian disk”
modelling. The bottom-left panels are for the “Keplerian $+$ radial” models.
The bottom-right panels denote the “Keplerian $+$ Ulrich” models.

##  4 Discussion

Our results make the strongest case to date for the central stellar mass of
the G10.6 massive star forming clump to be very large, $M_{\star}=120$ to
$200~{}M_{\odot}$. The modelling results are consistent with the evidence
presented by Sollins et al. [49], who used NH3 inversion lines with a
resolution similar to ours, but originating in the molecular infall zone
exterior to the central ionized disk, to infer a central mass $M_{\star}\sim
150$ M⊙. From our kinematical modelling, this mass is distributed within a
radius $<1000$ to 1500 au, since this is the part of the ionized structure
that is dominated by Keplerian rotation. Moreover, the mass has to be
(proto)stellar, because the amount of ionized gas needed to reproduce the
continuum and H$30\alpha$ brightness is relatively small. The best-fit models
give a total amount of ionized gas of only 0.2 to $0.25~{}M_{\odot}$ within
the entire domain of length 11480 au. Besides constraining the central mass of
the cluster, three main questions remain: $(i)$ what is the distribution of
(proto)stellar masses?; $(ii)$ what is the accretion stage of these
(proto)stars?; and $(iii)$ what is the fate of this compact cluster?

For question $(i)$, the ionizing-photon rate $N_{\mathrm{Ly}}$ inferred from
the free-free continuum gives independent constraints on the most massive
stars. Our models give a range $N_{\mathrm{Ly}}=6.1$ to $6.6\times 10^{48}$
s-1, whereas the spherical estimation by Sollins et al. [49] gives a
significantly larger $N_{\mathrm{Ly}}=2\times 10^{50}$ s-1. The main reason
for this discrepancy is the rapid density decrease in the vertical direction
in our models, which are tailored to fit the averaged midplane cuts of the
flattened ionized structure. The UC Hii region is far from spherical222The
spherical power law with a density index $n\propto r^{-1.5}$ assumed by
Sollins et al. [49] also implies an ionized-gas mass $0.5~{}M_{\odot}$, about
$\times 2$ larger than in our models., therefore the true value for
$N_{\mathrm{Ly}}$ should be somewhere in between, although probably closer to
our estimate. Using the stellar calibrations of Martins et al. [34], the range
$N_{\mathrm{Ly}}=10^{49}$ to $10^{50}$ s-1 is equivalent to the output of one
Zero Age Main Sequence (ZAMS) star with mass $M_{\star}\approx 32$ to $\gtrsim
60~{}M_{\odot}$. The corresponding stellar luminosity is $L_{\star}\approx
2\times 10^{5}$ to $1\times 10^{6}~{}L_{\odot}$, which is well within the
bolometric luminosity of the G10.6 region reported by Lin et al. [29],
$L_{\mathrm{FIR}}\approx 3\times 10^{6}~{}L_{\odot}$.

The kinematical and luminosity constraints are satisfied by a number of
stellar-mass arrangements, e.g., three stars of $50~{}M_{\odot}$ each provide
enough gravity with a total $N_{\mathrm{Ly}}\approx 1\times 10^{50}$ s-1. In
the scenario of having $N_{\mathrm{Ly}}$ closer to our lower estimate $\sim
1\times 10^{49}$ s-1, one star with $M_{\star}\approx 35~{}M_{\odot}$ would be
enough to provide the ionizing photons, but the kinematical constraints would
require the presence of an unknown number or lower-mass stars. Assuming that
this unseen stellar cluster follows a Kroupa [26] IMF, we calculate that the
median cluster mass corresponding to a median maximum stellar mass
$M_{\mathrm{\star,max}}=35~{}M_{\odot}$ is $M_{\mathrm{cl}}\approx
500~{}M_{\odot}$, subject to significant stochasticity. We conclude that it is
unlikely that a cluster sampling the full IMF is forming in such a reduced
volume (see below), but an unseen group of lower-mass stars might coexist with
the ionizing sources. Further observations at higher angular resolution from
the radio to the mid-IR with the JWST are needed to clarify the distribution
of (proto)stellar masses.

For question $(ii)$, the current accretion stage depends on the interpretation
of the ionized-gas kinematics. While the evidence for infall and rotation in
the molecular gas beyond the central ionized structure is strong [49, 31], our
modelling shows that the radial motions in the ionized gas – regardless of
their direction – become important starting somewhere in between $r\sim 1000$
to 1500 au, or about $r\sim R_{g}$. This suggests outward radial motions in
the ionized gas. This scenario would be similar to that of an ionized disk
wind surrounded by collapsing molecular gas. A few cases of ionized disks have
been reported around young massive stars [e.g., 12, 35, 10, 17], but for
objects that are less massive and luminous. Since the aforementioned objects
are less embedded, their most likely interpretation is that of a remnant
ionized disk without active accretion to the central star. The situation in
G10.6 could be different. Even if the mass reservoir in the ionized disk is
smaller than a solar mass, exterior to it exists a reservoir of infalling
molecular gas that surpasses the central stellar mass beyond a radius of $\sim
0.1$ pc, and reaches up to $M_{\mathrm{gas}}\sim 2500$ M⊙ at a radius of 0.5
pc [see the mass profiles in 31]. The molecular gas could replenish the
ionized reservoir and the observed boundary of the ionized and molecular
emission will depend on the asymmetric interactions of the ionizing photons
and the surrounding molecular gas [41, 7]. Under the disk plus Ulrich envelope
scenario, the timescale for the replenishment of available $0.2~{}M_{\odot}$
of ionized gas is $t_{\mathrm{rep}}=M_{\mathrm{ion-
gas}}/\dot{M}_{\mathrm{env}}\approx 6000$ yr, or $\times 10$ shorter than the
expected star-formation timescale. Therefore, we conclude that even if
replenishment of ionized gas occurs, it can provide at most a few $M_{\odot}$
of fresh material. Any active accretion ought to be in a residual stage.

For question $(iii)$, the mass density of this young and compact cluster is
large. Taking the upper and lower limits in mass ($M_{\star}=120$ to
$200~{}M_{\odot}$) and radius ($r_{\star}=1000$ to 1500 au), we estimate
$\rho_{\star}\sim 7.4\times 10^{7}$ to $4.2\times 10^{8}$ $M_{\odot}$ pc-3.
Constraining the stellar number density is challenging because of the
aforementioned lack of knowledge of the lower-mass stellar population, but a
number of stars in the range 5 to 20 translates into $n_{\star}\sim 1\times
10^{7}$ pc-3. Our lower limit to $\rho_{\star}$ is slightly larger than the
value derived by Schödel et al. [48] for the innermost 0.01 pc of the nuclear
star cluster around Sgr A∗, $\rho_{\mathrm{GC}}=2.6\times 10^{7}$ $M_{\odot}$
pc-3. However, the Galactic Center cluster is composed of $\sim 200$ solar-
type stars. Our loose constraints on the number densities indicate that
stellar interactions are likely to occur within the next Myr [38, 53].

##  5 Conclusions

Using ALMA, we report the first kinematically resolved observations of an
ionized rotating disk around a forming star cluster. The target is the
luminous star formation region G10.6-0.4. We use radiative transfer models of
the 1.3-mm free-free continuum and H$30\alpha$ line emission to constrain the
density and velocity structure within a radius of 6000 au. Our preferred best-
fit model is that of an ionized Keplerian disk with radial motions beyond a
radius $R_{g}\sim 1000$ to 1500 au, which corresponds to the radius within
which the ionized gas is expected to be bound. The central stellar mass is
robustly constrained to be in the range $M_{\star}=120$ to $200~{}M_{\odot}$.
The ionized-gas mass is only $M_{\mathrm{ion-gas}}=0.2$ to $0.25~{}M_{\odot}$.
The viewing inclination angle from face-on is in the range $i=49$ to
$56~{}\deg$. The fitted radial motions $v_{r,0}=8.7$ km s-1 correspond exactly
to the sound speed ($c_{s}=8.6$ km s-1) of ionized gas at $T_{e}=9000$ K,
indicating that the outer ionized gas is barely unbound at most. From
constraints on the amount of ionizing photons and FIR luminosity, we conclude
that there are either a few massive stars with $M_{\star}=32$ to
$60~{}M_{\odot}$, or one such massive star accompanied by an unknown number of
lower-mass stars. Any active accretion of ionized gas onto the (proto)stars is
mostly residual. The inferred cluster density is large, which suggests that
stellar interactions are likely to occur within the next Myr .

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00106.S.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. R.G.-M.
and C.C.-G. acknowledge support from UNAM-PAPIIT projects IN108822 and
IG101321, and from CONACyT Ciencia de Frontera project ID: 86372. R.G-M also
acknowledges support from the AAS Chrétien International Research Grant.
H.B.L. is supported by the National Science and Technology Council (NSTC) of
Taiwan (Grant Nos. 108-2112-M-001-002-MY3, 111-2112-M-001-089-MY3, and
110-2112-M-001-069). A.G. acknowledges support from NSF AAG 2008101 and NSF
CAREER 2142300.

## References

  * Astropy Collaboration et al. [2018] Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Beltrán et al. [2011] Beltrán, M. T., Cesaroni, R., Neri, R., & Codella, C. 2011, A&A, 525, A151, doi: 10.1051/0004-6361/201015049
  * Dullemond et al. [2012] Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015.  http://ascl.net/1202.015
  * Foreman-Mackey [2016] Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. [2013] Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Franco et al. [1989] Franco, J., Tenorio-Tagle, G., & Bodenheimer, P. 1989, Rev. Mexicana Astron. Astrofis., 18, 65
  * Galván-Madrid et al. [2011] Galván-Madrid, R., Peters, T., Keto, E. R., et al. 2011, MNRAS, 416, 1033, doi: 10.1111/j.1365-2966.2011.19101.x
  * Giacobbo & Mapelli [2018] Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011, doi: 10.1093/mnras/sty1999
  * Goddi et al. [2020] Goddi, C., Ginsburg, A., Maud, L. T., Zhang, Q., & Zapata, L. A. 2020, ApJ, 905, 25, doi: 10.3847/1538-4357/abc88e
  * Guzmán et al. [2020] Guzmán, A. E., Sanhueza, P., Zapata, L., Garay, G., & Rodríguez, L. F. 2020, ApJ, 904, 77, doi: 10.3847/1538-4357/abbe09
  * Ho & Haschick [1986] Ho, P. T. P., & Haschick, A. D. 1986, ApJ, 304, 501, doi: 10.1086/164184
  * Hoare [2006] Hoare, M. G. 2006, ApJ, 649, 856, doi: 10.1086/506961
  * Hollenbach et al. [1994] Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654
  * Hosokawa et al. [2010] Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478, doi: 10.1088/0004-637X/721/1/478
  * Ilee et al. [2018] Ilee, J. D., Cyganowski, C. J., Brogan, C. L., et al. 2018, ApJ, 869, L24, doi: 10.3847/2041-8213/aaeffc
  * Izquierdo et al. [2018] Izquierdo, A. F., Galván-Madrid, R., Maud, L. T., et al. 2018, MNRAS, 478, 2505, doi: 10.1093/mnras/sty1096
  * Jiménez-Serra et al. [2020] Jiménez-Serra, I., Báez-Rubio, A., Martín-Pintado, J., Zhang, Q., & Rivilla, V. M. 2020, ApJ, 897, L33, doi: 10.3847/2041-8213/aba050
  * Johnston et al. [2015] Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19, doi: 10.1088/2041-8205/813/1/L19
  * Keto [2002] Keto, E. 2002, ApJ, 580, 980
  * Keto [2003] —. 2003, ApJ, 599, 1196
  * Keto [2007] —. 2007, ApJ, 666, 976, doi: 10.1086/520320
  * Keto & Wood [2006] Keto, E., & Wood, K. 2006, ApJ, 637, 850, doi: 10.1086/498611
  * Keto & Zhang [2010] Keto, E., & Zhang, Q. 2010, MNRAS, 406, 102, doi: 10.1111/j.1365-2966.2010.16672.x
  * Keto et al. [1987] Keto, E. R., Ho, P. T. P., & Haschick, A. D. 1987, ApJ, 318, 712, doi: 10.1086/165405
  * Klaassen et al. [2018] Klaassen, P. D., Johnston, K. G., Urquhart, J. S., et al. 2018, A&A, 611, A99, doi: 10.1051/0004-6361/201731727
  * Kroupa [2001] Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Krumholz et al. [2014] Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, Protostars and Planets VI, 243, doi: 10.2458/azu_uapress_9780816531240-ch011
  * Law et al. [2021] Law, C. J., Zhang, Q., Öberg, K. I., et al. 2021, ApJ, 909, 214, doi: 10.3847/1538-4357/abdeb8
  * Lin et al. [2016] Lin, Y., Liu, H. B., Li, D., et al. 2016, ApJ, 828, 32, doi: 10.3847/0004-637X/828/1/32
  * Liu et al. [2010] Liu, H., Ho, P. T. P., Zhang, Q., et al. 2010, ApJ, 722, 262, doi: 10.1088/0004-637X/722/1/262
  * Liu [2017] Liu, H. B. 2017, A&A, 597, A70, doi: 10.1051/0004-6361/201629582
  * Liu et al. [2011] Liu, H. B., Zhang, Q., & Ho, P. T. P. 2011, ApJ, 729, 100, doi: 10.1088/0004-637X/729/2/100
  * Lu et al. [2022] Lu, X., Li, G.-X., Zhang, Q., & Lin, Y. 2022, Nature Astronomy, 6, 837, doi: 10.1038/s41550-022-01681-4
  * Martins et al. [2005] Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049, doi: 10.1051/0004-6361:20042386
  * Maud et al. [2018] Maud, L. T., Cesaroni, R., Kumar, M. S. N., et al. 2018, A&A, 620, A31, doi: 10.1051/0004-6361/201833908
  * McMullin et al. [2007] McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Mendoza et al. [2004] Mendoza, S., Cantó, J., & Raga, A. C. 2004, Rev. Mexicana Astron. Astrofis., 40, 147.  https://arxiv.org/abs/astro-ph/0401426
  * Moeckel & Bally [2007] Moeckel, N., & Bally, J. 2007, ApJ, 661, L183, doi: 10.1086/518738
  * Olguin et al. [2020] Olguin, F. A., Hoare, M. G., Johnston, K. G., et al. 2020, MNRAS, 498, 4721, doi: 10.1093/mnras/staa2406
  * Palla & Stahler [1993] Palla, F., & Stahler, S. W. 1993, ApJ, 418, 414
  * Peters et al. [2010a] Peters, T., Banerjee, R., Klessen, R. S., et al. 2010a, ApJ, 711, 1017, doi: 10.1088/0004-637X/711/2/1017
  * Peters et al. [2012] Peters, T., Longmore, S. N., & Dullemond, C. P. 2012, MNRAS, 425, 2352, doi: 10.1111/j.1365-2966.2012.21676.x
  * Peters et al. [2010b] Peters, T., Mac Low, M.-M., Banerjee, R., Klessen, R. S., & Dullemond, C. P. 2010b, ApJ, 719, 831, doi: 10.1088/0004-637X/719/1/831
  * Pringle [1981] Pringle, J. E. 1981, ARA&A, 19, 137, doi: 10.1146/annurev.aa.19.090181.001033
  * Rohlfs & Wilson [2000] Rohlfs, K., & Wilson, T. L. 2000, Tools of radio astronomy (Tools of radio astronomy / K. Rohlfs, T.L. Wilson. New York : Springer, 2000\. (Astronomy and astrophysics library,ISSN0941-7834))
  * Sanna et al. [2014] Sanna, A., Reid, M. J., Menten, K. M., et al. 2014, ApJ, 781, 108, doi: 10.1088/0004-637X/781/2/108
  * Sanna et al. [2019] Sanna, A., Kölligan, A., Moscadelli, L., et al. 2019, A&A, 623, A77, doi: 10.1051/0004-6361/201833411
  * Schödel et al. [2018] Schödel, R., Gallego-Cano, E., Dong, H., et al. 2018, A&A, 609, A27, doi: 10.1051/0004-6361/201730452
  * Sollins et al. [2005] Sollins, P. K., Zhang, Q., Keto, E., & Ho, P. T. P. 2005, ApJ, 631, 399, doi: 10.1086/432503
  * Tanaka et al. [2013] Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155, doi: 10.1088/0004-637X/773/2/155
  * Tanaka et al. [2017] Tanaka, K. E. I., Tan, J. C., & Zhang, Y. 2017, ApJ, 835, 32, doi: 10.3847/1538-4357/835/1/32
  * Ulrich [1976] Ulrich, R. K. 1976, ApJ, 210, 377
  * Zinnecker & Yorke [2007] Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481, doi: 10.1146/annurev.astro.44.051905.092549

##  Appendix A Diagnostic plots

In this Appendix we show diagnostic plots for the best-fit models of each of
the considered scenarios: purely Keplerian disk (Fig. 4), Keplerian disk with
radial motions starting at $r_{0}$ (Fig. 5), and Keplerian disk with an
external Ulrich-type envelope (Fig. 6).

Figure 4: Diagnostic plots for best fit of purely Keplerian model. The top-
left panel shows the midplane ion density $n_{ion}=n_{e}$, along with a cut in
which the velocity field is measured. The top-right panel shows cartesian
velocity components in the midplane $z=0$ cut. The bottom-left panel shows the
midplane density profile. The bottom-right panel shows the scale-height $H$
profile.  Figure 5: Same as Fig. 4, but for the case of the case of a
Keplerian disk with external radial motions. Figure 6: Same as Fig. 4, but for
the case of the case of a Keplerian disk with an exterior Ulrich envelope.
Note that the Ulrich envelope has a divergence at the radius in which it
settles into a disk [52, 37], which we smooth by averaging a few neighboring
grid points.