# Observational signatures of forming young massive clusters:
continuum emission from dense H ii regions

Mutsuko Inoguchi,1 Takashi Hosokawa,2 Hajime Fukushima,3, Kei E. I. Tanaka,4,
Hidenobu Yajima,3, Shin Mineshige,1
1Department of Astronomy, Graduate School of Science, Kyoto University, Sakyo,
Kyoto 606-8502, Japan
2Department of Physics, Graduate School of Science, Kyoto University, Sakyo,
Kyoto 606-8502, Japan
3Center for Computational Sciences, University of Tsukuba, Ten-nodai, 1-1-1
Tsukuba, Ibaraki 305-8577, Japan
4Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
Meguro, Tokyo, 152-8551, Japan E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Young massive clusters (YMCs) are the most massive star clusters forming in
nearby galaxies and are thought to be a young analogue to the globular
clusters. Understanding the formation process of YMCs leads to looking into
very efficient star formation in high-redshift galaxies suggested by recent
JWST observations. We investigate possible observational signatures of their
formation stage, particularly when the mass of a cluster is increasing via
accretion from a natal molecular cloud. To this end, we study the broad-band
continuum emission from ionized gas and dust enshrouding YMCs, whose formation
is followed by recent radiation-hydrodynamics simulations. We perform post-
process radiative transfer calculations using simulation snapshots and find
characteristic spectral features at radio and far-infrared frequencies. We
show that a striking feature is long-lasting, strong free-free emission from a
$\sim$ 10 pc-scale H ii region with a large emission measure of $\gtrsim
10^{7}\mathrm{cm}^{-6}\ \mathrm{pc}$, corresponding to the mean electron
density of $\gtrsim 10^{3}~{}\mathrm{cm}^{-3}$. There is a turnover feature
below $\sim$ 10 GHz, a signature of the optically-thick free-free emission,
often found in Galactic ultra-compact H ii regions. These features come from
the peculiar YMC formation process, where the cluster’s gravity effectively
traps photoionized gas for a long duration and enables continuous star
formation within the cluster. Such large and dense H ii regions show distinct
distribution on the density-size diagram, apart from the standard sequence of
Galactic H ii regions. This is consistent with the observational trend
inferred for extragalactic H ii regions associated with YMCs.

###### keywords:

keyword1 – keyword2 – keyword3

††pubyear: 2015††pagerange: Observational signatures of forming young massive
clusters:  continuum emission from dense H ii regions–Observational signatures
of forming young massive clusters:  continuum emission from dense H ii regions

##  1 Introduction

Young massive clusters (YMCs), also known as superstar clusters, are the most
massive star clusters forming in the present-day universe (Portegies Zwart et
al., 2010; Longmore et al., 2014). Their typical mass ($\gtrsim
10^{4}~{}\mathrm{M}_{\odot}$) and density ($\gtrsim
10^{3}~{}\mathrm{M}_{\odot}$ pc-3) are much larger than those of open clusters
(OCs), typical star clusters in our Galaxy (Beck, 2015). YMCs are often found
in nearby starburst and interacting galaxies (e.g. Whitmore et al., 1993,
1999; Whitmore, 2000), and they are considered a young analogue to the
globular clusters. Therefore, revealing the YMC formation processes in the
nearby galaxies leads to glimpsing the star formation in the early universe.
Some latest JWST observations resolve individual YMCs in distant galaxies even
at redshifts $z\simeq 4-6$ (e.g. Vanzella et al., 2022b, a), suggesting that
the YMCs should contribute larger fractions of the star formation activities
in such galaxies.

Star clusters in their early formation stage are called “embedded clusters"
(Lada & Lada, 2003). Previous optical and infrared (IR) observations report
candidates of YMCs in the embedded stage in nearby galaxies (e.g. Gorjian et
al., 2001; Beck et al., 2002; Galliano et al., 2005). The cold gas and dust
associated with these candidates are confirmed by recent high-resolution
observations using Atacama Large Millimeter/sub-millimeter Array (ALMA)
(Johnson et al., 2015; Turner et al., 2017; Leroy et al., 2018; Finn et al.,
2019). More recently, the latest JWST observations report discoveries of new
populations of YMCs, many of which seem to be in an early embedded stage
(Whitmore et al., 2023). Recent ALMA observations by He et al. (2022) report
possible candidates of individual $\sim 10$ pc-scale H ii regions associated
with embedded YMCs in the Antennae galaxy. Theoretical models suggest that
YMCs and progenitors of globular clusters often form in highly pressurized
environments (Elmegreen & Efremov, 1997; Kruijssen, 2015), and some
observations point to cloud-cloud collisions as a key physical process (Tsuge
et al., 2021a, b). However, the exact mechanism of the YMC formation,
including the evolution of a cluster under radiative feedback from high-mass
cluster-member stars (Krumholz et al., 2019), remains uncertain.

Recently, Fukushima & Yajima (2021) (FY21 hereafter) studied conditions
required for the YMC formation, systematically performing a suite of
radiation-hydrodynamic (RHD) simulations of the cluster formation and cloud
destruction by the stellar radiative feedback (see also, e.g. Dale et al.,
2012; Kim et al., 2018; He et al., 2019; Dobbs et al., 2020; Grudić et al.,
2021; Menon et al., 2022). As a result of examining cases with different
parameters, such as the initial cloud mass, size, metallicity, etc., they
found the critical gas surface density above which the YMC formation occurs,
$\Sigma\simeq 200~{}\mathrm{M}_{\odot}{\rm pc}^{-2}$. Recent observations
report molecular clouds with such high surface densities are ubiquitously
found in distant galaxies at redshifts $z\sim 1$ (Dessauges-Zavadsky et al.,
2023). In the cases of the YMC formation, a star cluster rapidly becomes
massive enough to gravitationally trap the hot photoionized gas, before the
entire cloud is blown away by an expanding H ii bubble. The cluster mass
continues to increase even after the emergence of the H ii bubble for such a
case (Bressert et al., 2012), resulting in a high star-formation efficiency
(SFE), or the mass ratio of the cluster to the natal molecular cloud.

FY21 suggest characteristic features of the formation stage of the YMC, i.e.,
the stage when the cluster mass increases through accretion from the natal
cloud. The gravity from the growing massive cluster keeps the H ii bubble
compact and dense around it for a long duration (Keto, 2003). We aim to derive
observational signatures of such a characteristic evolutionary stage. To this
end, we conduct various post-process radiation transfer calculations using
simulation snapshots. In this paper, we particularly investigate the continuum
spectra expected for the YMC formation stage. We study the free-free radio
spectrum from dense H ii regions associated with forming YMCs, with which the
mean density of the photoionized gas has been observationally inferred (e.g.
Garay & Lizano, 1999; Kim & Koo, 2001). Intriguingly, Hunt & Hirashita (2009)
(hereafter HH09) demonstrate that H ii regions found in Blue Compact Dwarf
galaxies obey the density-size relation apart from that for Galactic H ii
regions (see also Gilbert & Graham, 2007). They show that such extragalactic H
ii regions are distinctively denser than their Galactic counterparts with
similar sizes, suggesting a different population of H ii regions in extreme
starburst environments. We show that our post-process calculations for OC- and
YMC-forming cases well explain such observational trends regarding the
density-size relations.

We organize the rest of the paper as follows. In Section 2, we briefly review
our numerical methods on the radiation hydrodynamic simulations by FY21. We
also describe the method for post-process radiation transfer calculations to
give the continuum spectra based on the simulation data. In Section 3.1, we
compare the two representative cases of the OC and YMC formation, which we
mainly consider throughout the paper. In Section 3.2, we consider the
evolution of the continuum spectra from radio to IR wavelengths, with which we
consider possible observational signatures of the YMC formation. Finally,
Sections 4 and 5 provide discussions and concluding remarks.

##  2 Methods

###  2.1 Radiation-hydrodynamics simulations

We here briefly describe the method of three-dimensional (3D) RHD simulations
of the cluster formation and cloud destruction (see also FY21, ). They use a
modified version of the adaptive mesh refinement (AMR) code sfumato
(Matsumoto, 2007; Matsumoto et al., 2015), for which the M1 method (e.g.
Levermore, 1984) is implemented to handle the radiative transfer (dubbed
sfumato-m1). They adopt the reduced speed of light approximation with
$\tilde{c}=3\times 10^{-4}c$, where $c$ is the light speed (Rosdahl et al.,
2013). Photoionization of atoms, photodissociation of molecules, and radiative
heating of gas and dust around a light source are solved by considering the
transfer of extreme ultraviolet (EUV; $13.6~{}{\rm eV}<h\nu$), Lyman-Werner
(LW; $11.2~{}{\rm eV}<h\nu<13.6~{}{\rm eV}$), far-ultraviolet (FUV; $6~{}{\rm
eV}<h\nu<13.6~{}{\rm eV}$), and infrared photon (IR) photons (e.g., Hosokawa &
Inutsuka, 2006). The chemistry solver is based on the scheme of Sugimura et
al. (2020), and it is extended to include the network of Nelson & Langer
(1997) for CO formation. The chemical network also includes O0, O+, and O2+,
whose abundances are solved with the same procedure as in Fukushima et al.
(2020).

We insert sink particles representing a small star cluster when the density
exceeds a threshold value and other conditions are satisfied (Federrath et
al., 2010). We assign photon emissivity to each sink particle in the following
two ways. One is the same as in FY21, where the luminosity and spectrum are
given by taking the averages of the stellar isochrone of Chen et al. (2015)
and the Chabrier initial mass function (IMF, Chabrier, 2003). The other is
stochastic sampling, where the total cluster mass is distributed into stellar
mass bins in a probability-weighted manner based on the IMF (Fukushima &
Yajima, 2022).

###  2.2 Cases examined

Table 1: Cases considered Model | $M_{\rm cl}$ | $\Sigma_{\rm cl}$ | $R_{\rm cl}$ | $n_{0}$ | $t_{\mathrm{ff}}$
---|---|---|---|---|---
| (M⊙) | (M${}_{\odot}/$pc2) | (pc) | (cm-3) | (Myr)
M6R28st (YMC) | $10^{6}$ | 400 | 28.2 | 309 | 2.5
M6R56st (OC) | $10^{6}$ | 100 | 56.4 | 38.6 | 7.0
M6R10 (YMC) | $10^{6}$ | 3200 | 10.0 | 7000 | 0.52
M5R20 (OC) | $10^{5}$ | 80 | 20.0 | 87 | 4.7
M6R60 (OC) | $10^{6}$ | 88 | 60.0 | 32 | 7.7
Figure 1: Birth of a star cluster and destruction of a cloud simulated in the
representative case of the OC formation, M6R56st. The columns of panels
display snapshots at the different epochs, (1) $t=7.2$ Myr, (2) 9.2 Myr, (3)
10.8 Myr, and (4) 14.6 Myr from left to right. The top, middle, and bottom
rows display the distributions of the total gas surface density $\Sigma$,
column density of electrons $N_{\mathrm{el}}$, and that of hydrogen molecules
$N_{\mathrm{H}_{2}}$ measured along the line of sight, $z$-axis. The white
dots represent the distributions of the star (cluster) particles.  Figure 2:
The same as Fig. 1 but for the representative case of the YMC formation,
M6R28st. The columns of panels correspond to the different epochs of (1)
$t=2.5$ Myr, (2) 3.2 Myr, (3) 3.7 Myr, and (4) 5.0 Myr from left to right, the
same times as in Fig. 1 if normalized by the initial free-fall timescale
$t_{\mathrm{ff}}$. Figure 3: Star formation histories against the elapsed time
normalized by the free-fall timescale for the representative cases of the OC
and YMC formation, M6R56st and M6R28st (panels a and b). In each panel, the
solid and dashed lines represent the cluster mass normalized by the initial
cloud mass $M_{*}/M_{\mathrm{cl}}$ and mass fraction of gravitationally bound
particles $f_{\mathrm{bd}}$, respectively.

Table 1 summarizes the cases considered. Each simulation run starts from a
homogeneous cloud with the mass $M_{\mathrm{cl}}$, surface density
$\Sigma_{\mathrm{cl}}$, and radius $R_{\mathrm{cl}}$. We also assume that a
turbulent velocity field fills the cloud. Throughout the paper, we
particularly focus on two representative cases, M6R28st and M6R56st, for which
the cloud mass is the same as $M_{\rm cl}=10^{6}~{}\mathrm{M}_{\odot}$. Case
M6R28st starts from a relatively compact cloud with $\Sigma_{\rm
cl}=400~{}\mathrm{M}_{\odot}\mathrm{pc}^{-2}$, or with the radius
$R_{\mathrm{cl}}\simeq 28$ pc. The free-fall timescale for the initial cloud
is $t_{\mathrm{ff}}\simeq$ 2.5 Myr. The other case M6R56st starts from a cloud
with $\Sigma_{\rm cl}=100~{}\mathrm{M}_{\odot}\mathrm{pc}^{-2}$, or with the
radius $R_{\mathrm{cl}}\simeq 56$ pc. The corresponding free-fall timescale is
$t_{\mathrm{ff}}\simeq$ 7 Myr. As shown in Section 3.1 below, M6R28st and
M6R56st represent the typical cases of yielding YMC-like and OC-like clusters,
respectively. We newly perform these simulation runs using the stochastic
stellar sampling developed in Fukushima & Yajima (2022) (Section 2.1) as
indicated by the label "st". Kim et al. (2016) show that the stochasticity of
the stellar population impacts the emissivity only when the total stellar mass
is smaller than $10^{4}~{}\mathrm{M}_{\odot}$ (see also, Fukushima & Yajima,
2022). In our work, the cluster mass is much higher than this value for all
the cases examined, and these cases provide the essentially same results as in
FY21.

Apart from the above cases, we consider three additional cases M6R10, M5R20,
and M6R60, identical to those studied in FY21, where the model names have
extensions of "Z0A1". For instance, case M6R10 corresponds to M6R10Z0A1 in
FY21.

###  2.3 Post-process continuum radiative transfer calculations

For post-process radiative transfer calculations, we convert the original AMR
data of a simulation run into a data cube composed of $(128)^{3}$ cells. These
cells homogeneously distribute in the whole computational domain. We map all
the physical quantities on the original AMR grids to the cartesian grids by
linear interpolation. We solve the transfer equation at a given frequency,
which ranges from $\nu_{\rm min}=0.1$ GHz to $\nu_{\rm max}=10^{5}$ GHz, along
a line of sight chosen as an axis of the data cube,

$\frac{dI_{\nu}}{ds}=-(\kappa_{\nu,{\rm ff}}+\kappa_{\nu,{\rm
d}})I_{\nu}+\kappa_{\nu,{\rm ff}}B_{\nu}(T)+\kappa_{\nu,{\rm d}}B_{\nu}(T_{\rm
d}),$ (1)

where $I_{\nu}$ is the intensity, $\kappa_{\nu,{\rm ff}}$ the free-free
opacity, $\kappa_{\nu,{\rm d}}$ the dust opacity, $B_{\nu}$ the Planck
function, and $T$ and $T_{\rm d}$ are the gas and dust temperatures. The free-
free opacity is written as

$\kappa_{\nu,{\rm ff}}=\frac{j_{\nu,{\rm ff}}}{B_{\nu}(T)},$ (2)

where $j_{\nu,{\rm ff}}$ is the free-free emission coefficient, which is given
by

$\displaystyle j_{\nu,{\rm ff}}$ $\displaystyle\simeq$ $\displaystyle
5.4\times 10^{-41}{\rm erg}{\rm s}^{-1}{\rm cm}^{-3}{\rm Hz}^{-1}{\rm
rad}^{-1}$ $\displaystyle\times\ \
T_{4}^{-1/2}n_{e}n_{p}\exp\left(-\frac{h\nu}{k_{\rm B}T}\right)g_{\rm ff},$

where $T_{4}\equiv(T/10^{4}~{}{\rm K})$, $n_{e}$ and $n_{p}$ the electron and
proton number densities, $h$ the Planck constant, $k_{\rm B}$ the Boltzmann
constant, $g_{\rm ff}$ the Gaunt factor (e.g. Rybicki & Lightman, 1986;
Draine, 2011). Regarding the dust opacity $\kappa_{\nu,{\rm d}}$, we use

$\kappa_{\nu,\rm d}=0.1\left(\frac{\lambda}{300\mu{\rm m}}\right)^{-2}{\rm
cm}^{2}{\rm g}^{-1},$ (4)

where $\lambda$ is the wavelength. The functional form of Eq. (4) is valid for
the far-IR range, and it has been typically assumed for the Herschel Infrared
Galactic plane survey (Hi-GAL, e.g. Molinari et al., 2010; Elia et al., 2013),
which covers $70~{}\micron\lesssim\lambda\lesssim 500~{}\micron$. Note that in
our RHD simulations the dust temperature $T_{\rm d}$ is computed by solving
the monochromatic IR transfer using the Planck mean opacity given by Laor &
Draine (1993).

##  3 Results

###  3.1 Bimodal evolution: OC and YMC formation

Here we briefly review the evolution for cases M6R56st and M6R28st,
representatives of the OC and YMC formations observed in our RHD simulations.
FY21 also provide similar descriptions in more detail.

Fig. 1 shows the time evolution of case M6R56st, where the surface density of
the cloud is relatively low, $\Sigma_{\rm
cl}=100~{}\mathrm{M}_{\odot}\mathrm{pc}^{-2}$. In the early phase, the initial
turbulent motion controls the gas dynamics and induces filamentary structures
along which the star formation occurs. An H ii region appears around the star
cluster by the epoch of $t\simeq t_{\mathrm{ff}}$. The natal cloud gradually
disperses owing to the expanding H ii region during $t_{\mathrm{ff}}\lesssim
t\lesssim 2t_{\mathrm{ff}}$. We see that the distribution of star cluster
particles extends during that period. Most hydrogen molecules are destroyed by
stellar FUV radiation (Inoguchi et al., 2020). The star cluster is only
surrounded by photoionized gas at the final snapshot of $t\simeq
2t_{\mathrm{ff}}$.

Fig. 2 shows the evolution of case M6R28st, where the initial cloud surface
density is higher than the above as $\Sigma_{\rm
cl}=400~{}\mathrm{M}_{\odot}\mathrm{pc}^{-2}$. While the basic evolution looks
similar to that presented in Fig. 1, a denser and more compact star cluster
finally appears in this case. As studied in FY21, the gravitational force from
the star cluster overwhelms the thermal pressure gradient caused by the H ii
regions. The electron density near the cluster centre remains much higher than
in case M6R56st. We derive observational signatures of such large and dense H
ii regions in Section 3.2 below.

Fig. 3 presents the time evolution of the ratio of the cluster mass $M_{*}$ to
the initial cloud mass $M_{*}/M_{\mathrm{cl}}$ and the bound fraction of the
star cluster $f_{\mathrm{bd}}$, the mass fraction of star-cluster particles
that are gravitationally bound. Both cases show a common evolution that
$M_{*}/M_{\mathrm{cl}}$ reaches 0.1 at $t\simeq 1.3~{}t_{\mathrm{ff}}$.
However, the subsequent evolution is very different among these models.
Whereas $M_{*}/M_{\mathrm{cl}}$ saturates at $\simeq 0.13$ at $t\simeq
1.5~{}t_{\mathrm{ff}}$ for case M6R56st, it continuously increases to reach
0.6 until $t\simeq 4~{}t_{\mathrm{ff}}$ for case M6R28st. Early theoretical
studies generally predict that the bound fraction of the star cluster is a
sensitive function of the SFE; high $f_{\mathrm{fd}}$ is realized with the SFE
exceeding a few $\times$ 10 % after the dispersal of a natal cloud (e.g.,
Adams, 2000; Baumgardt & Kroupa, 2007; Shukirgaliyev et al., 2017). Similarly,
FY21 show that the evolution when $M_{*}/M_{\mathrm{cl}}\sim 0.1$ determines
the subsequent fate, whether it leads to the YMC formation. In the former case
of M6R56st (see also Fig. 1), $M_{*}/M_{\mathrm{cl}}$ is relatively small at
$t\simeq 1.3~{}t_{\mathrm{ff}}$, so that the bound fraction remains less than
0.1 throughout the simulation. The UV radiative feedback effectively
suppresses the star formation no later than $t\simeq 1.5~{}t_{\mathrm{ff}}$.
In the latter case of M6R28st (see also Fig. 2), $M_{*}/M_{\mathrm{cl}}$
becomes relatively large before the stellar UV feedback becomes effective,
i.e., before the H ii bubble sweeps a large part of the original cloud. The
bound fraction exceeds 0.9 at $t\simeq 1.5~{}t_{\mathrm{ff}}$, and the
gravitational potential around the star cluster becomes deep enough to capture
the ambient gas.

As mentioned in Section 1, FY21 show that different initial cloud properties
lead to the bimodal evolution overviewed for the above two cases M6R56st and
M6R28st. We thus particularly study observational signatures of these cases,
supposing that M6R56st and M6R28st demonstrate the typical YMC and OC
formation processes, respectively.

###  3.2 Continuum emission from ionized gas and dust around forming clusters

####  3.2.1 Continuum spectra

Figure 4:  Time evolution of the continuum spectrum emitted from a central
part of cluster-forming regions. Panels (a) and (b) represent the cases of
M6R56st and M6R28st, where the YMC and OC formation eventually occur. The
vertical axis represents the intensity averaged over a part where the
intensity is higher than 0.1 % of the peak value at a given frequency. In both
panels, the darker line colors represent the snapshots in the later stages,
the same snapshots as in Figs. 1 and 2. The direction of the vertical arrows
denotes how the evolution proceeds. The thin dashed curves represent the
reference spectrum of the Planck function $B_{\nu}(T_{\mathrm{eff}})$ with the
effective temperatures $T_{\mathrm{eff}}=1000$ K, 100 K, and 10 K.  Figure 5:
EM map for the representative cases of the OC and YMC formation, M6R56st (top
row) and M6R28st (bottom row). In each row, the left and right panels
illustrate snapshots at the same epochs as in the second and fourth columns of
Figs. 1 and 2. The lines of sight are also the same as these figures. The
dashed circle in each panel denotes the half-mass radius of the star cluster
measured from its mass center.

Fig. 4 shows the temporal evolution of the continuum spectra for the
representative cases of the normal OC formation (top, case M6R56st) and YMC
formation (bottom, case M6R28st), corresponding to the same snapshots as in
Figs. 1 and 2. These panels show the contributions from regions where the
local intensity is more than 0.1 % of the peak values for given frequencies.
As shown later in Section 3.2.2, the size of the emitting regions is $\sim 10$
pc for the YMC-forming case and $\sim$ a few $\times$ 10 - 100 pc for the OC-
forming case. The overall shape of the spectrum is common for these cases; one
component of the dust thermal emission at $\nu\gtrsim 10^{2}-10^{3}$ GHz, and
the other component of the free-free emission at the lower frequencies.

Despite the similarity in the spectrum shape, there are striking differences
among these cases, for instance, in the time evolution. In case M6R56st of the
OC formation, the emission gradually declines at all frequencies for $1\leq
t/t_{\mathrm{ff}}\leq 2.1$, during which the SFE saturates to $\simeq 0.2$
(Fig. 3). The brightness temperature at $\simeq 30$ GHz is always less than 10
K throughout the evolution. This is expected with the standard picture of the
H ii bubble expansion; the mean density (and also the column density) of the
photoionized gas decreases as the bubble expands.

In case M6R28st of the YMC formation, in contrast, the emission oppositely
continues to get stronger for almost the same duration of $1\leq
t/t_{\mathrm{ff}}\leq 2$, particularly at $\nu\gtrsim 10^{3}$ GHz
corresponding to the dust thermal emission. The free-free component at the
lower frequencies never decreases but rather increases slightly. The
brightness temperature at $\simeq 30$ GHz is much higher than for case
M6R56st, staying at $\sim 100$ K. These features come from the characteristic
evolution of the YMC formation described in Section 3.1. The strong gravity of
the forming massive cluster prevents the density of the photoionized gas from
decreasing. Recall that the cluster mass continues to increase even after
$t=2~{}t_{\mathrm{ff}}$ (Fig. 3).

Another difference is found at the low-frequency end of the free-free
continuum spectrum for $\nu\lesssim 10$ GHz. Whereas the spectrum is almost
flat for case M6R56st of the OC formation, it decreases with decreasing the
frequency for case M6R28st of the YMC formation. This corresponds to the
optically thick regime of the free-free emission. Such a feature is known to
appear below the turnover frequency

$\nu_{\rm to}\simeq 16.0~{}{\rm GHz}\left(\frac{{\rm
EM}}{10^{9}\mathrm{cm}^{-6}\
\mathrm{pc}}\right)^{0.48}\left(\frac{T_{\mathrm{i}}}{10^{4}~{}\mathrm{K}}\right)^{-0.64}$
(5)

(e.g. Mezger & Henderson, 1967; Kurtz, 2005), where $T_{\mathrm{i}}$ is the
temperature of the ionized gas and EM represents the emission measure defined
as

$\mathrm{EM}\equiv\int n_{e}^{2}\ ds,$ (6)

for which the integration is performed along the lines of sight. Fig. 5 shows
the evolution of the EM distribution for the cases considered above. In case
M6R56st of the OC formation, the peak value of the EM is $\simeq
10^{7}~{}\mathrm{cm}^{-6}\ \mathrm{pc}$ in an early stage at $t\simeq
1.3~{}t_{\mathrm{ff}}$, and it decreases by orders of magnitude by the epoch
of $\simeq 2.1~{}t_{\mathrm{ff}}$. This is consistent with Fig. 4 (a), where
the absorption feature only appears at $\nu\lesssim 1$ GHz in the earliest
stage and gradually disappears. In contrast, Fig. 5 shows that for case
M6R28st of the YMC formation the EM takes the much higher peak values at $\sim
10^{9}~{}\mathrm{cm}^{-6}\ \mathrm{pc}$, which hardly decreases until the
latest stage of $t=2~{}t_{\mathrm{ff}}$. The high EM is due to the trapping
effect of the photoionized gas by the gravity of the forming massive cluster.
The persistent turnover feature at $\nu\lesssim 10$ GHz in Fig. 4 (b) agrees
with the EM evolution presented in Fig. 5.

####  3.2.2 Density-size diagram of H ii regions

Figure 6:  Mean densities and diameters of H ii regions estimated from
emission-measure maps at the frequency of 36 GHz. Panels (a) and (b) represent
different threshold intensities for taking a central part, 0.1% and 1% of the
peak value. In both panels, the large blueish and reddish circles represent
cases of M6R28st and M6R56st, where the YMC and OC formation eventually occur.
The darker colors represent the later evolutionary stages as in Fig. 4. We
also present additional cases of M6R10, M5R20, and M6R60 with different
symbols of the small squares, diamonds, and triangles. M6R10 corresponds to
the YMC-formation case, and M5R20 and M6R60 to the OC-formation cases. The
darker colors indicate the later stages of $t=1.0$, 1.3, 1.5, and 2.1
$t_{\mathrm{ff}}$ also for these cases. In each panel, the dashed and dotted
lines represent the density-size relations of Galactic H ii regions compiled
by Garay & Lizano (1999) and Kim & Koo (2001), respectively. The solid line
corresponds to the relation for the extragalactic ultra-dense H ii regions
given by HH09. The open red star symbols represent actual samples obtained by
radio interferometric observations.

One of the important properties provided by radio observations is the density-
size relation of H ii regions, which we consider with Fig. 6. It has been well
established that the Galactic H ii regions obey the density-size relation
$n_{\mathrm{e}}\propto D^{-1}$, where $n_{\mathrm{e}}$ is the mean electron
number density of an H ii region and $D$ is its diameter (Garay & Lizano,
1999; Kim & Koo, 2001). The same correlation continues over H ii regions with
sizes that differ by more than four orders of magnitude, from $\sim 0.1$ pc of
ultra-compact H ii regions to over 100 pc of giant regions. The black dashed
and dotted straight lines in Fig. 6 represent the best fits to the Galactic
samples given by Garay & Lizano (1999) and Kim & Koo (2001), respectively. Its
inclination is somehow shallower than that expected by the simple Strömgren-
sphere argument assuming a constant emissivity of ionizing photons,
$n_{\mathrm{e}}\propto D^{-1.5}$.

Whether the similar correlation exists for extragalactic H ii regions has been
also studied for decades (Kennicutt, 1984). Interestingly, some H ii regions
associated with YMCs are known to locate in a distinct region in the density-
size diagram. Early studies have already suggested the presence of pc-scale H
ii regions with high densities above
$n_{\mathrm{e}}>10^{3}~{}\mathrm{cm}^{-3}$ (Kobulnicky & Johnson, 1999; Beck
et al., 2002), which are only found for ultra-compact ($D\lesssim
0.1~{}\mathrm{pc}$) regions in the Galaxy (Churchwell, 2002). HH09 compile a
large sample of such H ii regions in the literature and investigate their
statistical properties. The open star symbols in Fig. 6 represent their radio
samples, many of which are H ii regions found in blue compact dwarf galaxies.
HH09 also point out that their density-size relation also approximately follow
$n_{\mathrm{e}}\propto D^{-1}$, which corresponds to the solid straight line
in Fig. 6.

We superpose our simulation data points in Fig. 6. To do that, we extract an
area where the radio intensity at 36 GHz is higher than some threshold values,
0.1% and 1% of the peak value for each snapshot (upper and lower panels). We
have confirmed that our results do not change if we shift the frequency to 150
GHz, which is expected by the flat spectrum in this range (Fig. 4). We compute
the effective radius of the area as

$r=\sqrt{\frac{S}{\pi}},$ (7)

where $S$ is the area surface, and the mean emission measure
$\langle\mathrm{EM}\rangle$ by simply taking the arithmetic average of
$\mathrm{EM}$ over cells contained within the area. The mean density $\langle
n_{\mathrm{e}}\rangle$ is accordingly given by

$\langle n_{\mathrm{e}}\rangle=\sqrt{\frac{\langle\mathrm{EM}\rangle}{D}},$
(8)

where $D$ is the effective diameter $D=2r$.

As expected with the bimodal evolution for the OC- and YMC-forming cases
(Section 3.1), the simulation data points for these cases separately
distribute in Fig. 6. For case M6R56st of the OC formation (large reddish
circles), for instance, the points reside in the relatively lower right
portion of each panel. The points go to the lower-right side as time goes on,
representing the standard evolution of an H ii bubble expansion where the
electron number density gradually decreases with increasing time (Spitzer,
1978). These points locate near the density-size relation for the Galactic H
ii regions. Although the inclination of these points is steeper than
$n_{\mathrm{e}}\propto D^{-1}$, reducing the threshold intensity makes them
spread widely around the Galactic correlation (dotted and dashed) lines.

In stark contrast, data points for case M6R28st representing the YMC formation
(large bluish circles) are in the relatively upper left portion of each panel,
near the extragalactic radio samples by HH09. Moreover, these data points
hardly move for $t\lesssim 2~{}t_{\mathrm{ff}}$, unlike the OC-forming case
described above. These remarkable features come from the peculiar evolution of
the YMC formation suggested by RHD simulations; an H ii bubble bound by the
strong gravity of a newborn cluster postpones the dynamic expansion and
retains the high electron density. Fig. 6 suggests that our numerical
simulations successfully explain the YMC formation occurring in nearby
starburst galaxies.

In addition to the above two cases representing the OC and YMC formation, we
also perform the same post-process calculations for other cases studied in
FY21 (Table 1). As indicated by small filled symbols representing these cases,
the distribution of the simulation data is bimodal, reflecting the bimodal
evolution of the OC- and YMC-forming cases. We conclude that the distinct
distribution of extragalactic H ii regions on the density-size diagram is the
signature of the YMC formation, i.e., that of the gravity-bounded long-lasting
dense H ii regions.

####  3.2.3 Dust thermal emission

Figure 7: Dust continuum emission maps for the representative case of the OC
formation, M6R56st. The top, middle, and bottom rows show the time sequence of
the distributions of the column-density weighted dust temperature, and
emission intensities at the wavelengths $\lambda=70$ and 300µm. The columns of
panels show the snapshots at the same epochs as in Fig. 1. In each panel, the
dashed circle denotes the half-mass radius of the star cluster measured from
its mass center.  Figure 8: Mean radial distributions of the dust continuum
emission intensities for the representative case of the OC formation, M6R56st.
In each panel, the points represent the intensities averaged over rings
centered on the cluster mass center. The blue and red colors correspond to the
different wavelengths of $\lambda=70$ and 300µm. The panels show the snapshots
at the same epochs as in Figs. 1 and 7. The vertical dashed line in each panel
represents the half-mass radius of the star cluster.  Figure 9: The same as in
Fig. 7 but for the representative case of the YMC formation, M6R28st. The
columns of panels show the snapshots at the same epochs as in Fig. 2.  Figure
10: The same as Fig. 8 but for the representative case of the YMC formation,
M6R28st. The panels show the snapshots at the same epochs as in Figs. 2 and 9.

As already suggested by Fig. 4, the dust thermal emission at far-IR
wavelengths also shows the characteristic features of the OC and YMC
formation. Fig. 7 shows the evolution of the column-density-weighted dust
temperature $\bar{T}_{\mathrm{d}}$ and continuum emission maps at
$\lambda=70~{}\micron$ and $300~{}\micron$ for the representative case of the
OC formation, M6R56st. The dust temperature is an indicator of the strength of
the local stellar radiation field, the heating source for grains. In fact, the
dust temperature is highest near the cluster centre at each snapshot. The peak
value gradually rises for $t_{\mathrm{ff}}\lesssim t\lesssim
1.5~{}t_{\mathrm{ff}}$, during which the cluster mass increases (Fig. 3).
However, it decreases for $1.5t_{\mathrm{ff}}\lesssim t\lesssim
2~{}t_{\mathrm{ff}}$, because the cluster spatially expands only with a slight
increase in mass (Fig. 1). The dust temperature at $\sim 100$ pc away from the
cluster centre monotonically increases as the cluster expands and the cloud
disperses instead. For a given snapshot, the emission at
$\lambda=70~{}\micron$ is more centrally concentrated than at $300~{}\micron$
because it traces the distribution of warmer grains. The radial extent of the
emission at these wavelengths gradually expands with time as the cluster
expands. The emission is always strongest near the cluster centre, but its
peak values substantially decrease for $1.5t_{\mathrm{ff}}\lesssim t\lesssim
2~{}t_{\mathrm{ff}}$. This is because the H ii bubble expansion reduces the
column density of the gas and dust, in addition to the decrease of
$\bar{T}_{\mathrm{d}}$. Note that the dust-to-gas mass ratio is homogeneous
within the computational domain in our simulations as in FY21. Fig. 8 also
clearly shows such evolution, presenting mean radial distributions of the dust
continuum emission intensities at each snapshot.

Figs. 9 and 10 are the same as Figs. 7 and 8 but for the representative case
of the YMC formation, M6R28st. Reflecting the different evolution outlined in
Section 3.1 among these cases, these figures also show the different features
from those described above. For instance, Figs. 9 shows that the dust
temperature $\bar{T}_{\mathrm{d}}$ near the cluster centre is much higher than
in case M6R56st for $t\gtrsim 1.3~{}t_{\mathrm{ff}}$, because the cluster mass
continues to increase until the final snapshot. Since the cluster does not
expand, unlike case M6R56st, the stellar radiation field at the cluster center
gets stronger and stronger while the cluster grows in mass. The evolution of
the emission maps at $\lambda=70$ and $300~{}\micron$ also differs from the
above. Most remarkably, the maps in the later stages of $t\gtrsim
1.5~{}t_{\mathrm{ff}}$ show strong and centrally concentrated emission at both
wavelengths. In addition to the high $\bar{T}_{\mathrm{d}}$, this is also
because of the difference in the gas density structure. Since the cluster’s
gravity is strong enough to even trap the photoionized gas, there remains the
dense gas in the vicinity of the cluster. The column densities of the gas and
dust hardly decrease for a long time. Fig. 10 suggests that for $t\gtrsim
1.5~{}t_{\mathrm{ff}}$ the peak intensities in case M6R28st are higher than
those in case M6R56st by a few orders of magnitudes. At the final snapshot at
$t\simeq 2~{}t_{\mathrm{ff}}$, there are sharp cusps within the central $\sim$
10 pc, which is in stark contrast to the flat core-like distributions seen for
case M6R56st (c.f. Fig. 8).

##  4 Discussions

###  4.1 Alternative signatures: line emission and absorption

Our focus has been on the continuum spectra to derive characteristic
signatures of the YMC formation, but this should not be the only available
observational clue. Alternative ways include the line emission against the
continuum spectra, which is to be further studied. Various lines spanning from
radio to infrared wavelengths are expected to emerge from different gas phases
surrounding YMCs during their formation stage. For instance, infrared emission
lines from a photodissociated region (PDR) are promising candidates (e.g.
Berné et al., 2022, and references therein). Similar to the dust thermal
emission studied in Section 3.2.3, these line emission maps will exhibit
distinct features that differentiate between OC- and YMC-forming cases. This
is one of our next projects in progress.

An advantage of the line features is allowing us to infer the gas kinematics
around the forming YMC. Our preparatory analyses of the simulation data show
that, for the YMC-forming cases, the cluster’s strong gravity causes the
infall motion of the surrounding gas toward the centre of an H ii region. Such
a signature of the infall motion has been really reported on Galactic ultra-
compact H ii regions (Sollins et al., 2005; Beltrán et al., 2006). These
studies make use of the absorption feature of foreground molecular (NH3)
inversion lines against the free-free continuum emission (Ho & Townes, 1983).
According to our simulations, the H ii region that emerges in the YMC
formation is much larger than the typical ultra-compact H ii regions by two
orders of magnitude. Nonetheless, they might show similar observational
signatures suggesting the infall motion in future observations.

Throughout this paper, we have mainly considered the YMC formation in the
nearby universe. As mentioned in Section 1, however, YMC formation may play a
more significant role in star formation in the early universe. Recent
observations with the JWST have found more galaxies at $z>10$ than expected
(e.g. Naidu et al., 2022; Harikane et al., 2023), which prompts intensive
studies and discussions (e.g. Boylan-Kolchin, 2022; Yajima et al., 2022;
Lovell et al., 2023; Dekel et al., 2023). One possible explanation is that YMC
formation, which achieves a high SFE, is the dominant mode of star formation
(Inayoshi et al., 2022). In observations of distant galaxies, restframe far-IR
line emissions such as [O i] 63$\micron$ and [C ii] 158$\micron$ coming from
PDRs are often used because they redshift into the ALMA bands. In the same
vein, [O iii] 88$\micron$ from H ii regions is a well-known tracer (Inoue et
al., 2016; Hashimoto et al., 2018; Witstok et al., 2022). All of these line
emissions can be predicted for both OC- and YMC-forming cases using our
simulation data. Although individual H ii regions cannot be resolved in
distant galaxies, if YMC formation is dominant for the overall star formation
in galaxies at $z>10$, it will be essential to comprehensively examine these
line emissions. In such extreme environments, however, other effects that we
have not considered in this paper, such as very low metallicity or a different
IMF (Zackrisson et al., 2011; Chon et al., 2021; Chon et al., 2022; Fukushima
& Yajima, 2023), may also come into play. It is also interesting to study
emission line diagnostics for YMC formation, taking into account these
effects.

###  4.2 Possible signatures in an earlier evolutionary stage

While our simulation data predict observational signatures consistent with the
radio observations (e.g. Fig. 6), the simulation runs start from the idealized
initial conditions. Such situations of isolated molecular clouds are not
necessarily realized in the YMC formation. Although the actual triggers of the
YMC formation are still uncertain, cloud-cloud collisions or large-scale
converging flows have been proposed as a possible process (e.g. Maeda et al.,
2021; Sameie et al., 2022; Dobbs et al., 2022). Our lack of knowledge about
the realistic initial conditions of YMC formation prevents us from predicting
observational signatures during the early stages of YMC formation. Rico-Villas
et al. (2020) and Rico-Villas et al. (2022) conduct diagnostic analyses with
many molecular lines detected toward “proto-superstar clusters” candidates in
NGC 253, and consider possible evolutionary sequences in the YMC formation.
They suggest that there is the “super hot core” stage in the YMC formation
before the appearance of an H ii region, similar to the hot core stage in the
individual high-mass star formation (Garay & Lizano, 1999; Hoare et al.,
2007). To improve the accuracy of our predictions and better compare them with
such observations, we need to update our simulations to reflect more realistic
initial conditions. Nevertheless, we expect that our results on the
observational signatures of the gravitationally-bound H ii regions should not
depend on details of the initial conditions, as long as the key process
enabling continuous star formation against the stellar UV feedback is
accurately captured. This also remains to be examined in future studies.

##  5 Summary and Conclusions

We have studied possible observational signatures of YMCs in their formation
stage, i.e., the evolutionary stage where the cluster mass grows through
accretion from the surrounding medium. In particular, we have considered the
continuum spectra from $\sim$ 10 pc-scale dense H ii regions created around
the accreting YMC, as suggested by recent RHD simulations (FY21). We have
performed post-process radiative transfer calculations using the simulation
snapshots to derive the spectra from radio to far-IR frequencies. For
comparisons, we have also performed the same post-process calculations for the
cases where normal OCs eventually appear instead of YMCs. Our findings are
summarized as follows.

For both simulation runs that represent OC and YMC formation, the continuum
spectrum commonly consists of two components: one dominated by the thermal
dust emission at $\nu\gtrsim 10^{2}-10^{3}$ GHz, and the other dominated by
the free-free emission at the lower frequencies. However, there are remarkable
differences between these cases, as illustrated in Fig. 4. The spectrum for
the normal OC-forming case (M6R56st) represents the standard picture of an H
ii bubble expansion, where the electron density gradually decreases as the
bubble expands. The intensities also decline from the radio to infrared for
$t_{\mathrm{ff}}\lesssim t\lesssim 2t_{\mathrm{ff}}$, during which the
radiative feedback by the expanding bubble destroys the natal cloud and
quenches the star formation. For the YMC-forming case (M6R28st), in contrast,
the intensities gradually rise for the same time interval normalized by
$t_{\mathrm{ff}}$ and are always much stronger than those for the OC-forming
case. Moreover, there is a turnover feature below $\sim 10$ GHz resulting from
the large emission measure. These all come from the peculiar evolution in the
YMC-forming case, i.e., the rapid star formation under weak radiative
feedback, which leads to the formation of a dense H ii region that remains
trapped by the cluster’s gravitational field.

Previous radio observations provide density-size relations of Galactic and
extragalactic H ii regions, which have been investigated to infer the variety
of high-mass star- and cluster-forming environments (Garay & Lizano, 1999; Kim
& Koo, 2001). Remarkably, it has been pointed out that some extragalactic H ii
regions distribute separately from Galactic H ii regions (HH09). They are
large ($\sim 10$ pc), dense ($\sim 10^{4}~{}\mathrm{cm}^{-3}$), and associated
with YMCs. We have superposed our simulation data points on the density-size
diagram based on the synthetic continuum emission maps at 36 GHz (Fig. 6).
Reflecting the qualitatively different evolution between the OC- and YMC-
forming cases, the corresponding data points distinctly distribute on the
diagram. For the normal OC-forming cases, the simulation points are situated
close to the Galactic density-size relationships. As the evolution progresses
to later stages, the points shift to the less-dense and larger-size portion,
in accordance with the standard expansion law of H ii bubbles. The points for
the YMC-forming cases, in contrast, scatter separately from the OC-forming
cases and near the observational data of extragalactic radio samples. These
points hardly move on the diagram regardless of $t/t_{\mathrm{ff}}$,
reflecting that the cluster’s gravity traps an H ii bubble and prevents it
from dynamically expanding for a while. We propose that previous radio
observations have already captured the signatures of the YMC formation
suggested by recent RHD simulations.

We suggest that similar trends to the above should be also found in the far-IR
continuum dust thermal emission maps (Figs. 7–10). We have investigated the
evolution of the maps at $\lambda=70~{}\micron$ and $300~{}\micron$. For the
normal OC-forming case, the radial emission distribution gradually broadens
with time. As a result, the emission peak associated with the cluster center
shows a decrease in intensity. A flat, core-like distribution persists within
the half-mass radius of the cluster. This can be attributed to the dynamic
expansion of the H ii bubble, as well as to the outward expansion of the star
cluster induced by the dispersal of the surrounding molecular cloud. The YMC-
forming case, in contrast, shows that the central emission peak does not
weaken but rather becomes stronger over time. The radial emission distribution
evolves to show a remarkably high degree of central concentration,
characterized by a sharp peak at the center of the cluster. This is again due
to the emergence of a dense H ii bubble gravitationally bound by the nascent
YMC, which continues to accrete mass.

Our study suggests that the peculiar YMC formation process found in recent RHD
simulations, i.e., the formation of gravitationally-bound dense H ii regions,
should indicate observational signatures. We conclude that such signatures
have been captured by extragalactic radio observations. As discussed in
Section 4, our work can be easily stretched to further studies that link
simulations and upcoming observations. Particularly, considering observational
signatures of emission and absorption lines associated with the YMC formation
is of great importance, and it is also a target of our next work.

## Acknowledgements

The authors sincerely appreciate Jeong-Gyu Kim, Hiroyuki Hirashita, Eric Keto,
Rolf Kuiper, Yurina Nakazato, Akino Inoue, and Shu-ichiro Inutsuka for the
comments and discussions. The numerical simulations were carried out on XC50
Aterui II at the Center for Computational Astrophysics (CfCA) of the National
Astronomical Observatory of Japan. This research could never be accomplished
without the support by Grants-in-Aid for Scientific Research (TH:19H01934,
21H00041, HF:23K13139) from the Japan Society for the Promotion of Science.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Adams (2000) Adams F. C., 2000, ApJ, 542, 964
  * Baumgardt & Kroupa (2007) Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589
  * Beck (2015) Beck S., 2015, International Journal of Modern Physics D, 24, 1530002
  * Beck et al. (2002) Beck S. C., Turner J. L., Langland-Shula L. E., Meier D. S., Crosthwaite L. P., Gorjian V., 2002, AJ, 124, 2516
  * Beltrán et al. (2006) Beltrán M. T., Cesaroni R., Codella C., Testi L., Furuya R. S., Olmi L., 2006, Nature, 443, 427
  * Berné et al. (2022) Berné O., et al., 2022, PASP, 134, 054301
  * Boylan-Kolchin (2022) Boylan-Kolchin M., 2022, arXiv e-prints, p. arXiv:2208.01611
  * Bressert et al. (2012) Bressert E., Ginsburg A., Bally J., Battersby C., Longmore S., Testi L., 2012, ApJ, 758, L28
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chen et al. (2015) Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068
  * Chon et al. (2021) Chon S., Omukai K., Schneider R., 2021, MNRAS, 508, 4175
  * Chon et al. (2022) Chon S., Ono H., Omukai K., Schneider R., 2022, MNRAS, 514, 4639
  * Churchwell (2002) Churchwell E., 2002, ARA&A, 40, 27
  * Dale et al. (2012) Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 424, 377
  * Dekel et al. (2023) Dekel A., Sarkar K. S., Birnboim Y., Mandelker N., Li Z., 2023, arXiv e-prints, p. arXiv:2303.04827
  * Dessauges-Zavadsky et al. (2023) Dessauges-Zavadsky M., et al., 2023, MNRAS, 519, 6222
  * Dobbs et al. (2020) Dobbs C. L., Liow K. Y., Rieder S., 2020, MNRAS, 496, L1
  * Dobbs et al. (2022) Dobbs C. L., Bending T. J. R., Pettitt A. R., Bate M. R., 2022, MNRAS, 509, 954
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Elia et al. (2013) Elia D., et al., 2013, ApJ, 772, 45
  * Elmegreen & Efremov (1997) Elmegreen B. G., Efremov Y. N., 1997, ApJ, 480, 235
  * Federrath et al. (2010) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, ApJ, 713, 269
  * Finn et al. (2019) Finn M. K., Johnson K. E., Brogan C. L., Wilson C. D., Indebetouw R., Harris W. E., Kamenetzky J., Bemis A., 2019, ApJ, 874, 120
  * Fukushima & Yajima (2021) Fukushima H., Yajima H., 2021, MNRAS, 506, 5512
  * Fukushima & Yajima (2022) Fukushima H., Yajima H., 2022, MNRAS, 511, 3346
  * Fukushima & Yajima (2023) Fukushima H., Yajima H., 2023, arXiv e-prints, p. arXiv:2303.12405
  * Fukushima et al. (2020) Fukushima H., Hosokawa T., Chiaki G., Omukai K., Yoshida N., Kuiper R., 2020, MNRAS, 497, 829
  * Galliano et al. (2005) Galliano E., Alloin D., Pantin E., Lagage P. O., Marco O., 2005, A&A, 438, 803
  * Garay & Lizano (1999) Garay G., Lizano S., 1999, PASP, 111, 1049
  * Gilbert & Graham (2007) Gilbert A. M., Graham J. R., 2007, ApJ, 668, 168
  * Gorjian et al. (2001) Gorjian V., Turner J. L., Beck S. C., 2001, ApJ, 554, L29
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Harikane et al. (2023) Harikane Y., et al., 2023, ApJS, 265, 5
  * Hashimoto et al. (2018) Hashimoto T., et al., 2018, Nature, 557, 392
  * He et al. (2019) He C.-C., Ricotti M., Geen S., 2019, MNRAS, 489, 1880
  * He et al. (2022) He H., Wilson C., Brunetti N., Finn M., Bemis A., Johnson K., 2022, ApJ, 928, 57
  * Ho & Townes (1983) Ho P. T. P., Townes C. H., 1983, ARA&A, 21, 239
  * Hoare et al. (2007) Hoare M. G., Kurtz S. E., Lizano S., Keto E., Hofner P., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. p. 181 (arXiv:astro-ph/0603560), doi:10.48550/arXiv.astro-ph/0603560
  * Hosokawa & Inutsuka (2006) Hosokawa T., Inutsuka S.-i., 2006, ApJ, 646, 240
  * Hunt & Hirashita (2009) Hunt L. K., Hirashita H., 2009, A&A, 507, 1327
  * Inayoshi et al. (2022) Inayoshi K., Harikane Y., Inoue A. K., Li W., Ho L. C., 2022, ApJ, 938, L10
  * Inoguchi et al. (2020) Inoguchi M., Hosokawa T., Mineshige S., Kim J.-g., 2020, MNRAS, 497, 5061
  * Inoue et al. (2016) Inoue A. K., et al., 2016, Science, 352, 1559
  * Johnson et al. (2015) Johnson K. E., Leroy A. K., Indebetouw R., Brogan C. L., Whitmore B. C., Hibbard J., Sheth K., Evans A. S., 2015, ApJ, 806, 35
  * Kennicutt (1984) Kennicutt R. C. J., 1984, ApJ, 287, 116
  * Keto (2003) Keto E., 2003, ApJ, 599, 1196
  * Kim & Koo (2001) Kim K.-T., Koo B.-C., 2001, ApJ, 549, 979
  * Kim et al. (2016) Kim J.-G., Kim W.-T., Ostriker E. C., 2016, ApJ, 819, 137
  * Kim et al. (2018) Kim J.-G., Kim W.-T., Ostriker E. C., 2018, ApJ, 859, 68
  * Kobulnicky & Johnson (1999) Kobulnicky H. A., Johnson K. E., 1999, ApJ, 527, 154
  * Kruijssen (2015) Kruijssen J. M. D., 2015, MNRAS, 454, 1658
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland -Hawthorn J., 2019, ARA&A, 57, 227
  * Kurtz (2005) Kurtz S., 2005, in Cesaroni R., Felli M., Churchwell E., Walmsley M., eds, IAU Symposium Vol. 227, Massive Star Birth: A Crossroads of Astrophysics. pp 111–119, doi:10.1017/S1743921305004424
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Laor & Draine (1993) Laor A., Draine B. T., 1993, ApJ, 402, 441
  * Leroy et al. (2018) Leroy A. K., et al., 2018, ApJ, 869, 126
  * Levermore (1984) Levermore C. D., 1984, J. Quant. Spectrosc. Radiative Transfer, 31, 149
  * Longmore et al. (2014) Longmore S. N., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 291–314 (arXiv:1401.4175), doi:10.2458/azu_uapress_9780816531240-ch013
  * Lovell et al. (2023) Lovell C. C., Harrison I., Harikane Y., Tacchella S., Wilkins S. M., 2023, MNRAS, 518, 2511
  * Maeda et al. (2021) Maeda R., Inoue T., Fukui Y., 2021, ApJ, 908, 2
  * Matsumoto (2007) Matsumoto T., 2007, PASJ, 59, 905
  * Matsumoto et al. (2015) Matsumoto T., Dobashi K., Shimoikura T., 2015, ApJ, 801, 77
  * Menon et al. (2022) Menon S. H., Federrath C., Krumholz M. R., 2022, MNRAS, 517, 1313
  * Mezger & Henderson (1967) Mezger P. G., Henderson A. P., 1967, ApJ, 147, 471
  * Molinari et al. (2010) Molinari S., et al., 2010, PASP, 122, 314
  * Naidu et al. (2022) Naidu R. P., et al., 2022, ApJ, 940, L14
  * Nelson & Langer (1997) Nelson R. P., Langer W. D., 1997, ApJ, 482, 796
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Rico-Villas et al. (2020) Rico-Villas F., Martín-Pintado J., González-Alfonso E., Martín S., Rivilla V. M., 2020, MNRAS, 491, 4573
  * Rico-Villas et al. (2022) Rico-Villas F., González-Alfonso E., Martín-Pintado J., Rivilla V. M., Martín S., 2022, MNRAS, 516, 1094
  * Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188
  * Rybicki & Lightman (1986) Rybicki G. B., Lightman A. P., 1986, Radiative Processes in Astrophysics
  * Sameie et al. (2022) Sameie O., et al., 2022, arXiv e-prints, p. arXiv:2204.00638
  * Shukirgaliyev et al. (2017) Shukirgaliyev B., Parmentier G., Berczik P., Just A., 2017, A&A, 605, A119
  * Sollins et al. (2005) Sollins P. K., Zhang Q., Keto E., Ho P. T. P., 2005, ApJ, 624, L49
  * Spitzer (1978) Spitzer L., 1978, Physical processes in the interstellar medium, doi:10.1002/9783527617722.
  * Sugimura et al. (2020) Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ, 892, L14
  * Tsuge et al. (2021a) Tsuge K., Fukui Y., Tachihara K., Sano H., Tokuda K., Ueda J., Iono D., Finn M. K., 2021a, PASJ, 73, S35
  * Tsuge et al. (2021b) Tsuge K., Tachihara K., Fukui Y., Sano H., Tokuda K., Ueda J., Iono D., 2021b, PASJ, 73, 417
  * Turner et al. (2017) Turner J. L., Consiglio S. M., Beck S. C., Goss W. M., Ho P. T. P., Meier D. S., Silich S., Zhao J.-H., 2017, ApJ, 846, 73
  * Vanzella et al. (2022a) Vanzella E., et al., 2022a, arXiv e-prints, p. arXiv:2211.09839
  * Vanzella et al. (2022b) Vanzella E., et al., 2022b, ApJ, 940, L53
  * Whitmore (2000) Whitmore B. C., 2000, arXiv e-prints, pp astro–ph/0012546
  * Whitmore et al. (1993) Whitmore B. C., Schweizer F., Leitherer C., Borne K., Robert C., 1993, AJ, 106, 1354
  * Whitmore et al. (1999) Whitmore B. C., Zhang Q., Leitherer C., Fall S. M., Schweizer F., Miller B. W., 1999, AJ, 118, 1551
  * Whitmore et al. (2023) Whitmore B. C., et al., 2023, ApJ, 944, L14
  * Witstok et al. (2022) Witstok J., et al., 2022, MNRAS, 515, 1751
  * Yajima et al. (2022) Yajima H., Abe M., Fukushima H., Ono Y., Harikane Y., Ouchi M., Hashimoto T., Khochfar S., 2022, arXiv e-prints, p. arXiv:2211.12970
  * Yang et al. (2021) Yang A. Y., et al., 2021, A&A, 645, A110
  * Zackrisson et al. (2011) Zackrisson E., Rydberg C.-E., Schaerer D., Östlin G., Tuli M., 2011, ApJ, 740, 13