#  Radiation hydrodynamics simulations of massive star cluster formation in
giant molecular clouds

Hajime Fukushima1, Hidenobu Yajima1
1Center for Computational Sciences, University of Tsukuba, Ten-nodai, 1-1-1
Tsukuba, Ibaraki 305-8577, Japan
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

By performing three-dimensional radiation hydrodynamics simulations, we study
the formation of young massive star clusters (YMCs,
$M_{*}>10^{4}~{}M_{\odot}$) in clouds with the surface density ranging from
$\Sigma_{\rm cl}=80$ to $3200~{}M_{\odot}\;{\rm pc^{-2}}$. We find that
photoionization feedback suppresses star formation significantly in clouds
with low surface density. Once the initial surface density exceeds $\sim
100~{}M_{\odot}\;{\rm pc^{-2}}$ for clouds with $M_{\rm
cl}=10^{6}~{}M_{\odot}$ and $Z=Z_{\odot}$, most of the gas is converted into
stars because the photoionization feedback is inefficient in deep
gravitational potential. In this case, the star clusters are massive and
gravitationally bounded as YMCs. The transition surface density increases as
metallicity decreases, and it is $\sim 350~{}M_{\odot}\;{\rm pc^{-2}}$ for
$Z=10^{-2}~{}Z_{\odot}$. We show that more than 10 percent of star-formation
efficiency (SFE) is needed to keep a star cluster gravitationally bounded even
after the disruption of a cloud. Also, we develop a semi-analytical model
reproducing the SFEs obtained in our simulations. We find that the SFEs are
fit with a power-law function with the dependency $\propto\Sigma_{\rm
cl}^{1/2}$ for low-surface density and rapidly increases at the transition
surface densities. The conditions of the surface density and the metallicity
match recent observations of giant molecular clouds forming YMCs in nearby
galaxies.

###### keywords:

stars: formation - stars: massive - stars: Population II - H II regions -
galaxies: star clusters: general - galaxies: star formation

††pubyear: 2020††pagerange:  Radiation hydrodynamics simulations of massive
star cluster formation in giant molecular clouds –E

##  1 Introduction

Young massive clusters (YMCs) form in giant molecular clouds (GMCs). Its mass
$(\gtrsim 10^{4}~{}M_{\odot})$ and density $(\gtrsim 10^{3}~{}M_{\odot}{\rm
pc^{-3}})$ are larger than typical open clusters in the Milky Way (Portegies
Zwart et al., 2010). YMCs have been frequently observed in local starburst or
merger galaxies (Longmore et al., 2014). Also, YMCs are likely to form in the
early Universe, and some of them remain as globular clusters (GCs) in the
present-day galaxies (e.g., Kruijssen, 2014). These massive dense clusters can
be main formation sites of massive stars which regulate the star formation in
galaxies via stellar feedback such as photoionization, stellar wind, and
supernovae (e.g., Schaye et al., 2010; Wise et al., 2012; Hasegawa & Semelin,
2013; Yajima et al., 2017; Yajima et al., 2020). Besides, YMCs are
gravitationally bounded and can produce black hole binaries via multi-body
gravitational interaction (e.g., Portegies Zwart & McMillan, 2000; Fujii et
al., 2017; Kumamoto et al., 2019), likely resulting in gravitational wave
events (Abbott et al., 2016). Therefore, understanding the formation of YMCs
is crucial in astronomy. However, their formation processes are still debated,
especially for the physical conditions of GMCs forming YMCs. (e.g., Longmore
et al., 2014).

As a cloud collapses gravitationally, gas converts into stars. Once massive
stars form, they give energy and momentum into the surrounding medium via
radiation, outflow, and stellar wind referred to as "feedback" (e.g., Krumholz
et al., 2019). These feedback processes hamper gravitational contraction and
star formation. Hence the stellar feedback can be a key factor in the
determination of star formation efficiency (SFE) that is defined by the ratio
of final stellar mass to initial cloud one. The SFE needs to be larger than
$15-30$ percent to keep the star cluster gravitationally bounded (e.g., Lada
et al., 1984; Kroupa et al., 2001; Baumgardt & Kroupa, 2007; Goodwin, 2009;
Shukirgaliyev et al., 2017). However, the SFEs are typically less than 10
percent in local galaxies (e.g., Murray, 2011; Vutisalchavakul et al., 2016;
Kruijssen et al., 2019; Chevance et al., 2020b). In the Milky Way, most stars
belong to star clusters initially, and then disperse gradually after
disruption of gas clouds (e.g., Lada & Lada, 2003). Therefore, special
conditions are likely to be required for the formation of YMCs.

The observations of nearby galaxies showed that the lifetimes of the GMCs are
less than $30~{}{\rm Myr}$ (Kawamura et al., 2009; Fukui & Kawamura, 2010).
Kruijssen et al. (2019) suggested that the duration time of star formation
could be as short as the dynamical time of clouds, comparing the spatial
correlation between star-forming molecular gas and Hii regions. This timescale
is shorter than the lifetime of massive stars. Therefore, star-forming clouds
can be dispersed by pre-supernovae feedback, such as the radiative feedback or
the stellar wind (e.g., Kruijssen et al., 2019; Chevance et al., 2020b; Kim et
al., 2020; Chevance et al., 2020a).

Extreme ultraviolet (EUV; $13.6~{}{\rm eV}\leqq h\nu\leqq 1~{}{\rm keV}$)
photons ionize and heat up gas. The high thermal pressure in Hii bubbles
pushes ambient neutral gas and allows rapid expansion of the bubbles (Spitzer,
1987; Hosokawa & Inutsuka, 2005, 2006). Analytical studies showed that the
photoionization feedback suppresses the star formation and makes the SFE lower
(e.g., Williams & McKee, 1997; Matzner, 2002; Krumholz & Matzner, 2009; Fall
et al., 2010; Kim et al., 2016; Inoguchi et al., 2020). If the density field
is inhomogeneous, numerical simulations are required to evaluate the impacts
of the photoionization feedback. In particular, turbulence motion in GMCs
induces the high-density contrast and the star formation in a stochastic
manner (Larson, 1981). Recently, radiation hydrodynamics simulations have been
performed to study star cluster formation and destruction of clouds (e.g.,
Vázquez-Semadeni et al., 2010; Dale et al., 2012, 2013; Howard et al., 2017;
Geen et al., 2017; Gavagnin et al., 2017; He et al., 2019; Decataldo et al.,
2020; Grudić et al., 2018; Grudić et al., 2020; Ali, 2021; Fujii et al.,
2021). Kim et al. (2018) showed that the photoionization feedback played a
dominant role in the suppression of star formation. They indicated that the
SFE increased with the surface density of clouds. In their simulations, the
SFE was less than 10 percent at $\Sigma\sim 100~{}M_{\odot}{\rm pc^{-2}}$,
which is the typical surface density of GMCs in the Milly Way. The
photoionization feedback becomes ineffective if the escape velocity from a
clouds exceeds the sound speed in Hii regions $c_{\rm s}\sim 10~{}{\rm
kms^{-1}}$ (Dale et al., 2012; Bressert et al., 2012). These massive star-
forming clouds typically form in the starburst and merger galaxies (e.g.,
Leroy et al., 2015, 2018; Tsuge et al., 2020). Recently, the numerical
simulations found that theses clouds could form in the high velocity colliding
flow ($\gtrsim 20~{}{\rm km/s}$, Dobbs et al., 2020; Liow & Dobbs, 2020; Maeda
et al., 2020). Such GMCs with high surface density are likely to be needed for
the formation of YMCs against the photoionization feedback.

GMCs are generally optically thick for UV radiation. Therefore, even if the
photoionization feedback is ineffective, the radiation pressure can regulate
the star formation as stellar luminosity increases (Krumholz & Matzner, 2009;
Fall et al., 2010; Murray et al., 2010; Kim et al., 2016). The energy absorbed
by the dust is reprocessed as thermal emission at infrared (IR) wavelengths.
If a cloud is optically thick even for IR photons, the additional radiation
pressure is exerted. This further suppresses the star formation (e.g., Murray
et al., 2010). Skinner & Ostriker (2015) performed RHD simulations with the
radiation pressure from IR photons. They suggested that the IR radiation
pressure could regulate star formation significantly only if the opacity of a
cloud is higher than $\kappa\sim 15~{}{\rm cm^{2}\,g^{-1}}$ that corresponds
to the super-solar environment.

The strength of radiative feedback sensitively depends on the metallicity of
gas as shown in Fukushima et al. (2020b, hereafter Paper I). The turbulent
motion induces a filamentary structure in which the dust column density can
become high enough to attenuate EUV radiation. Therefore, the expansion of Hii
bubbles is hampered because of the filaments in clouds enriched with metal and
dust (see Paper I). Whereas, in cases with the metallicity lower than $\sim
0.1~{}Z_{\odot}$, the filaments are easily disrupted due to the feedback, and
the star formation is quenched rapidly. Besides, the temperature of Hii
regions increases as the metallicity decreases because of the inefficient
metal line cooling such as Oii and Oiii (Draine, 2011). This results in the
faster propagation of ionizing front in the low-metallicity gas. Therefore,
the photoionization feedback is reinforced at the lower-metallicity (He et
al., 2019; Fukushima et al., 2020b).

The initial condition of turbulent motions can be another factor to control
the SFE. The turbulent motion is frequently characterized by the virial
parameter $\alpha_{\rm vir}$ defined by the ratio of kinetic to gravitational
energy (Bertoldi & McKee, 1992). Various values have been reported by
observations of star-forming clouds in the Milky Way and nearby galaxies
(Roman-Duval et al., 2010; Leroy et al., 2016; Sun et al., 2018, 2020).
Recently, Kim et al. (2021) showed that the SFEs decrease as the virial
parameter increases. However, they focused only on the clouds with the surface
density similar to the typical one in the Milky Way. Therefore, the impacts of
turbulent motions is still unclear in more massive compact clouds.

As stated above, the star formation in GMCs is related to various factors as
the radiative feedback, the metallicity, the virial parameter in addition to
the initial cloud mass and radius. Therefore, the formation of YMCs has not
been understood yet because of the complicated processes of them. In this
paper, we study the conditions for the YMC formation, performing the
simulations of clouds with various surface densities, metallicities, and
turbulent strength. We perform 3D RHD simulations with the SFUMATO-M1, the
modified version of a self-gravitational magnetohydrodynamics code with an
Eulerian adaptive mesh refinement (AMR), SFUMATO (Matsumoto, 2007; Matsumoto
et al., 2015). We have newly developed the radiation transfer scheme based on
the moment method with M1-closure. In Paper I, we limited the parameters of
clouds whose escape velocity is less than $10~{}{\rm kms^{-1}}$ to focus on
the effects of the photoionization feedback. We here study star formation in
more compact clouds in the range of surface densities,
$\Sigma=32-3200~{}{M_{\odot}\,{\rm pc^{-2}}}$. We also calculate non-
equilibrium chemical reactions with $\rm H_{2}$ and CO molecules. This allows
us to follow the spatial distribution of the molecules in each evolutionary
stage of a cloud.

We organize the rest of the paper as the following. In Section 2, we describe
the numerical method and the set-up of simulations. We show the results of
simulations in Section 3. Then, we develop a semi-analytical model with
simulation results in Section 4. Section 5 is summary and discussion. In
Appendix A and B, we show the detail of the chemical network and thermal
processes in our simulations. We describe the numerical methods of the moment-
based scheme of radiation transfer in Section C.

##  2 Numerical Method

Table 1: Models considered

model | $M_{\rm cl}\,[\,M_{\odot}\,]$ | $R_{\rm cl}\,[\,{\rm pc}\,]$ | $Z\,[\,Z_{\odot}\,]$ | $\alpha_{\rm vir}$ | $n_{\rm ini}\,[\,{\rm cm^{-3}}\,]$ | $\Sigma_{\rm cl}\,[\,M_{\odot}\,{\rm pc^{-2}}\,]$ | $\sigma_{0}\,[\,{\rm kms^{-1}}\,]$ | $t_{\rm ff}\,[\,{\rm Myr}\,]$ | $v_{\rm esc}\,[\,{\rm km/s}\,]$
---|---|---|---|---|---|---|---|---|---
M5R5Z0A1 | $10^{5}$ | $5$ | $1$ | 1 | $5600$ | $1300$ | $7.2$ | $0.59$ | $13$
M5R5Z0A2 | $10^{5}$ | $5$ | $1$ | 1 | $5600$ | $1300$ | $10$ | $0.59$ | $13$
M5R5Z-1A1 | $10^{5}$ | $5$ | $10^{-1}$ | 1 | $5600$ | $1300$ | $7.2$ | $0.59$ | $13$
M5R5Z-2A1 | $10^{5}$ | $5$ | $10^{-2}$ | 1 | $5600$ | $1300$ | $7.2$ | $0.59$ | $13$
M5R8Z0A1 | $10^{5}$ | $8$ | $1$ | 1 | $1400$ | $500$ | $5.7$ | $1.2$ | $10$
M5R10Z0A1 | $10^{5}$ | $10$ | $1$ | 1 | $700$ | $320$ | $5.1$ | $1.7$ | $9.3$
M5R10Z0A2 | $10^{5}$ | $10$ | $1$ | 2 | $700$ | $320$ | $7.2$ | $1.7$ | $9.3$
M5R10Z-1A1 | $10^{5}$ | $10$ | $10^{-1}$ | 1 | $700$ | $320$ | $5.1$ | $1.7$ | $9.3$
M5R10Z-2A1 | $10^{5}$ | $10$ | $10^{-2}$ | 1 | $700$ | $320$ | $5.1$ | $1.7$ | $9.3$
M5R12Z0A1 | $10^{5}$ | $12$ | $1$ | 1 | $410$ | $220$ | $4.6$ | $2.2$ | $8.5$
M5R20Z0A1 | $10^{5}$ | $20$ | $1$ | 1 | $87$ | $80$ | $3.6$ | $4.7$ | $6.6$
M5R20Z0A2 | $10^{5}$ | $20$ | $1$ | 2 | $87$ | $80$ | $5.1$ | $4.7$ | $6.6$
M5R20Z-1A1 | $10^{5}$ | $20$ | $10^{-1}$ | 1 | $87$ | $80$ | $3.6$ | $4.7$ | $6.6$
M5R20Z-2A1 | $10^{5}$ | $20$ | $10^{-2}$ | 1 | $87$ | $80$ | $3.6$ | $4.7$ | $6.6$
M6R10Z0A1 | $10^{6}$ | $10$ | $1$ | 1 | $7000$ | $3200$ | $16$ | $0.52$ | $29$
M6R10Z0A2 | $10^{6}$ | $10$ | $1$ | 2 | $7000$ | $3200$ | $23$ | $0.52$ | $29$
M6R10Z-1A1 | $10^{6}$ | $10$ | $10^{-1}$ | 1 | $7000$ | $3200$ | $16$ | $0.52$ | $29$
M6R10Z-2A1 | $10^{6}$ | $10$ | $10^{-2}$ | 1 | $7000$ | $3200$ | $16$ | $0.52$ | $29$
M6R17.5Z-2A1 | $10^{6}$ | $17.5$ | $10^{-2}$ | 1 | $1300$ | $1000$ | $12$ | $1.2$ | $22$
M6R20Z0A1 | $10^{6}$ | $20$ | $1$ | 1 | $870$ | $800$ | $11$ | $1.5$ | $20$
M6R20Z0A2 | $10^{6}$ | $20$ | $1$ | 2 | $870$ | $800$ | $16$ | $1.5$ | $20$
M6R20Z-1A1 | $10^{6}$ | $20$ | $10^{-1}$ | 1 | $870$ | $800$ | $11$ | $1.5$ | $20$
M6R20Z-2A1 | $10^{6}$ | $20$ | $10^{-2}$ | 1 | $870$ | $800$ | $11$ | $1.5$ | $20$
M6R25Z0A1 | $10^{6}$ | $25$ | $1$ | $1$ | $450$ | $510$ | $10$ | $2.1$ | $19$
M6R25Z-2A1 | $10^{6}$ | $25$ | $10^{-2}$ | $1$ | $450$ | $510$ | $10$ | $2.1$ | $19$
M6R30Z0A1 | $10^{6}$ | $30$ | $1$ | $1$ | $260$ | $350$ | $9.3$ | $2.7$ | $17$
M6R30Z-2A1 | $10^{6}$ | $30$ | $10^{-2}$ | $1$ | $260$ | $350$ | $9.3$ | $2.7$ | $17$
M6R32.5Z0A1 | $10^{6}$ | $32.5$ | $1$ | $1$ | $200$ | $300$ | $8.9$ | $3.1$ | $16$
M6R35Z0A1 | $10^{6}$ | $35$ | $1$ | $1$ | $160$ | $260$ | $8.6$ | $3.4$ | $16$
M6R35Z-2A1 | $10^{6}$ | $35$ | $10^{-2}$ | $1$ | $160$ | $260$ | $8.6$ | $3.4$ | $16$
M6R40Z0A1 | $10^{6}$ | $40$ | $1$ | 1 | $110$ | $200$ | $8.0$ | $4.2$ | $15$
M6R40Z0A2 | $10^{6}$ | $40$ | $1$ | 2 | $110$ | $200$ | $11$ | $4.2$ | $15$
M6R40Z-1A1 | $10^{6}$ | $40$ | $10^{-1}$ | 1 | $110$ | $200$ | $8.0$ | $4.2$ | $15$
M6R40Z-2A1 | $10^{6}$ | $40$ | $10^{-2}$ | 1 | $110$ | $200$ | $8.0$ | $4.2$ | $15$
M6R60Z0A1 | $10^{6}$ | $60$ | $1$ | 1 | $32$ | $88$ | $6.6$ | $7.7$ | $12$
M6R60Z0A2 | $10^{6}$ | $60$ | $1$ | 2 | $32$ | $88$ | $9.3$ | $7.7$ | $12$
M6R60Z-1A1 | $10^{6}$ | $60$ | $10^{-1}$ | 1 | $32$ | $88$ | $6.6$ | $7.7$ | $12$
M6R60Z-2A1 | $10^{6}$ | $60$ | $10^{-2}$ | 1 | $32$ | $88$ | $6.6$ | $7.7$ | $12$

Notes. Column 1: model names, Column 2: cloud masses, Column 3: cloud radii,
Column 4: metallicity, Column 5: virial parameters, Column 6: initial number
densities, Column 7: surface densities defined as $\Sigma_{\rm cl}=M_{\rm
cl}/(\pi R_{\rm cl}^{2})$, Column 8: three-dimensional velocity dispersions,
Column 9: free fall times, Column 10: escape velocity

To perform RHD simulations, we introduce radiative transfer (RT) into an
adaptive mesh refinement code, SFUMATO (Matsumoto, 2007; Matsumoto et al.,
2015), which is dubbed SFUMATO-M1. We develop the RT module based on the
moment equations with the M1 closure, which has been adopted in previous
studies (e.g., Rosdahl et al., 2013; Skinner & Ostriker, 2013; Takahashi &
Ohsuga, 2013; Rosdahl & Teyssier, 2015; Kannan et al., 2019). In the moment-
based scheme, the RT is calculated from the gradient of radiation energy
densities between adjacent cells. In Paper I, we used the RT module based on
the adaptive ray-tracing method (Abel & Wandelt, 2002) which was developed by
Sugimura et al. (2020). Its calculation amount is proportional to the number
of sources. Therefore, if the number exceeds $\sim 1000$, the computational
cost is quite expensive even with current computational facilities. This
regulates the parameter space with low-mass and low-density clouds. Thus, in
this study, we adopt the moment-based scheme of which the calculation amount
is independent of the number of sources.

In addition, in the case of massive compact clouds, the dust column density
becomes high and can be opaque even to infrared radiation from dust. This
gives the additional force on gas. Therefore, we here take into account the RT
of IR photons. The detailed methodology of RT and test simulations are
described in Appendix C. Also, SFUMATO includes the non-equilibrium chemistry
solver developed by Sugimura et al. (2020), allowing us to calculate the
thermal evolution of low-temperature gas accurately. The details of the
thermal processes and chemical network in our simulations are described in
Appendix A and B.

In this study, we set the size of calculation boxes as three times the cloud
radius $(R_{\rm cl})$ on a side. The maximum refinement level is fixed at
$l_{\rm max}=4$ and the minimum cell size is $\Delta x=0.059~{}{\rm pc}(R_{\rm
cl}/20~{}{\rm pc})$. The simulations end at when $4\times$ free-fall time
$(t_{\rm ff})$ elapses.

###  2.1 Basic Equations of hydrodynamics

We perform the three-dimensional hydrodynamics simulations with Cartesian
coordinate. We solve the following basic equations of compressible
hydrodynamics: the equation of continuity,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\bm{v}\right)=0,$ (1)

the equation of motion,

$\displaystyle\frac{\partial\left(\rho\bm{v}\right)}{\partial
t}+\nabla\cdot\left(\rho\bm{v}\otimes\bm{v}\right)+\nabla
P=\rho\left(\bm{g}+\bm{f}\right),$ (2)

and the energy equation,

$\displaystyle\frac{\partial\left(\rho E\right)}{\partial
t}+\nabla\cdot\left[\left(\rho
E+P\right)\bf{v}\right]=\rho\left({\bf{g}+\bf{f}}\right)\cdot{\bf{v}}+\Gamma-\Lambda,$
(3)

where $E$ is total energy defined as

$\displaystyle E=\frac{|{\bf
v}|^{2}}{2}+\left(\gamma-1\right)^{-1}\frac{P}{\rho},$ (4)

$\rho$, $P$, $\bf{v}$, $\bf{g}$, $\Gamma$ and $\Lambda$ are the density,
pressure, velocity, gravitational force, the heating and cooling functions. We
estimate the adiabatic exponent $\gamma$ as in Omukai & Nishi (1998). In
equations (2) and (3), $\bf{f}$ represents the radiation pressure force. We
consider radiation pressure induced by absorption of EUV photons by Hi and
dust grains, FUV and IR photons absorbed by dust grains.

We take into account the chemical networks of 11 species: $\rm H$, $\rm
H_{2}$, $\rm H^{-}$, $\rm H^{+}$, $\rm H_{2}^{+}$, $\rm e$, $\rm CO$, Cii, Oi,
Oii, and Oiii (see also, Fukushima et al., 2020b). The number density of the
$i$-th specie is calculated as

$\displaystyle\frac{\partial(y_{i}n_{\rm H})}{\partial
t}+\nabla\cdot(y_{i}n_{\rm H}{\bf{v}})=y_{i}n_{\rm H}R_{i},$ (5)

where $n_{\rm H}$ is the number density of hydrogen nuclei,
$y_{i}=n_{i}/n_{\rm H}$ is the fractional abundance of each chemical specie,
and $R_{i}$ is the reaction rate of the $i$-th specie.

We adopt the simple chemical network of Nelson & Langer (1997) for CO
formation, which has been used for the RHD simulations (e.g., Hosokawa &
Inutsuka, 2006; Nakatani et al., 2018). To obtain the relative abundances of
Oi, Oii, and Oiii, we adopt the following procedure (see also, Fukushima et
al., 2020a). We assume that the ionization rate of Oi is equal to that of Hi
because the ionization potential energies of Oi and Hi are similar. The
abundances of doubly ionized oxygen are determined as the chemical equilibrium
between Oii and Oiii.

As the heating/cooling processes, we include (1) the thermal processes related
to the chemical reactions, (2) line cooling of $\rm H_{2}$, Cii, CO, Oi, Oii,
and Oiii, and (3) energy transfer between gas and dust grains. As the
radiative processes, we include heating of Hi photoionization and $\rm H_{2}$
photodissociation. The dust grain temperature is estimated from the energy
balance between the absorption/emission of radiation and energy transfer with
gas. We set the temperature floor at $T=10~{}{\rm K}$ as Fukushima et al.
(2020b). All thermal processes are summarized in Appendix B.

###  2.2 Sink particles and radiation sources

If gas cells reach the maximum refinement level, they form sink particles as
star clusters. The model of the sink particles was developed in Matsumoto et
al. (2015). We summarize their procedures as follows. More details of them
were described in Appendix A of Matsumoto et al. (2015).

The production procedure of sink particles is the same as used in Fukushima et
al. (2020b). The conditions for the production are as follows (Federrath et
al., 2010): (1) the gas density is higher than the threshold value $\rho_{\rm
thr}$; (2) the birthplace is the local minimum of gravitational potential; (3)
the velocity divergence $\nabla\cdot\bf{v}$ and all the eigenvalues of the
symmetric part of the velocity gradient tensor $\nabla\bf{v}$ are negative;
(4) the sum of the thermal, kinetic and gravitational energy is negative. The
density threshold is set as $\rho_{\rm thr}=8.86c_{\rm s}^{2}/(\pi G\Delta
x^{2})=9.3\times 10^{-19}\,(R_{\rm cl}/20\,{\rm pc})^{-2}\,{\rm g\,cm^{-3}}$
(Gong & Ostriker, 2013; Kim et al., 2018) where $c_{\rm s}$ is the sound
speed, and $T=20~{}{\rm K}$ is used here. We set a sink radius as $r_{\rm
sink}=2\Delta x$. We do not consider mergers of sink particles in this work.

In a hydrodynamics step, each sink particle accretes the excess of the gas
mass over the threshold density $\rho_{\rm thr}$ in a sink radius. The sink
particles also obtain momentum from the accreted gas. Therefore, the total
momentum of gas and sink particles is conserved.

The gravitational force among the sink particles is considered for calculating
their motions. Here, we set the softening radius as $r_{\rm soft}=r_{\rm
sink}$. The gravitational force from gas is also taken into account. When a
cell is inside the softening radius of a sink particle, we subdivided each
cell into $8^{3}$ cells, and the gravitational forces from the subcells are
summed (Krumholz et al., 2004). The force from sink particles to gas is
evaluated as the back reactions of them. Thus the conservation of the momentum
is ensured through the interaction between sink particles and gas.

We estimate the luminosity of a sink particle from the stellar isochrone of
Chen et al. (2015) and the Chabrier IMF (Chabrier, 2003) with the stellar mass
range from $0.1~{}M_{\odot}$ to $150~{}M_{\odot}$. We use their isochrone at
$1~{}{\rm Myr}$, which roughly corresponds to star formation duration in the
simulations. In our model, the mass-luminosity relation to all sink particles
is fixed. However, we may overestimate the luminosity if the sink mass is not
massive enough to form a massive star that dominates the luminosity.
Therefore, as stellar radiation sources, we consider only sink particles with
masses higher than $50~{}M_{\odot}$ where the expected number of massive stars
($>10\,M_{\odot}$) is larger than unity. In our models, resolution dependence
of radiative properties appears when the typical sink mass is comparable to
the threshold value. As shown in Appendix E and D, it occurs in the diffuse
and low-mass cloud models (M5R20Z0A1 in Table 1) because the typical sink mass
is comparable to $\sim 50~{}M_{\odot}$. On the other hand, the peaks of the
sink mass distributions appear at $\sim 100~{}M_{\odot}$ in the cases of
$M_{\rm cl}=10^{6}~{}M_{\odot}$. In such a case, total emissivity of sink
particles does not depend on the numerical resolution significantly. Note that
the results of star cluster formation can also change with the models of
radiation sources, such as ones with the IMF-averaged values or stochastic
sampling (Grudić & Hopkins, 2019; Smith, 2021). To reproduce the IMF smoothly,
the total stellar mass of $\gtrsim 10^{4}~{}M_{\odot}$ is required (e.g., Kim
et al., 2016). Thus, our results for star clusters with the mass $\gtrsim
10^{4}~{}M_{\odot}$ are unlikely to be sensitive to the models of radiation
sources.

###  2.3 Initial conditions

As initial conditions, we consider the clouds with the mass $M_{\rm
cl}=10^{5}$ or $10^{6}\,M_{\odot}$, and the radius ranging from $R_{\rm cl}=5$
to $60~{}\rm pc$ as summarized in Table 1. The above parameters include the
typical values of surface density of GMCs in the Milky Way, $\Sigma\sim
100~{}{M_{\odot}\,{\rm pc^{-2}}}$ (e.g., Solomon et al., 1987; Roman-Duval et
al., 2010), and $\sim 3\times 10^{3}~{}M_{\odot}\,{\rm pc^{-2}}$ as in merger
galaxies (e.g., Leroy et al., 2016; Sun et al., 2018). Massive stars can form
only in massive clouds. Kauffmann & Pillai (2010) found that there is the
lower limit for massive star formation as $M_{\rm cl}\gtrsim
870~{}M_{\odot}(r/{\rm pc})^{1.33}$ in the Milky Way (see also Urquhart et
al., 2018). Our models fulfill this condition. In paper I, we focused only on
the clouds whose escape velocities are less than the sound speed of ionized
gas $(v_{\rm esc}<10~{}{\rm km/s})$ and investigated the cloud disruption
process due to the photo-ionization feedback. As the cloud becomes massive and
compact, the deep gravitation potential well can host gas against the
feedback, likely resulting in higher star formation efficiency. Therefore we
here extend the parameter range to take the clouds with $v_{\rm esc}>20~{}{\rm
km/s}$ into account. Also, we change the metallicity from $Z=10^{-2}Z_{\odot}$
to $Z_{\odot}$ and investigate its impacts on the star-formation process. The
initial gas temperature is set as $T_{\rm g}=10~{}{\rm K}$. The abundance of
molecules sensitively depends on the metallicity. We assume that the gas is
fully molecular at $Z=10^{-1}Z_{\odot}$ and $Z_{\odot}$. In the case with
$Z=10^{-2}~{}Z_{\odot}$, we set $y_{\rm H_{2}}=3\times 10^{-3}$ as the initial
abundance because the formation time scale of $\rm H_{2}$ is longer than the
dynamical time of the cloud.

The simulations finish when the elapsed time reaches four times the free-fall
time of the clouds.

As in paper I, we take into account the turbulent velocity field. We assume
that the velocity power spectrum of $P(k)\propto k^{-4}$ where $k$ is the wave
number. The velocity field is generated with the random seeds. We re-scale the
same velocity fields in each simulation to investigate feedback dependence on
the cloud mass and compactness. We note that the results, such as the SFEs,
are affected by choice of this seed (Grudić et al., 2020; Kim et al., 2021).
The amplitude of the turbulent motion is defined by the virial parameter as

$\displaystyle\alpha_{0}=\frac{E_{\rm kin}}{|2E_{\rm
grav}|}=\frac{5\sigma_{0}^{2}R_{\rm cl}}{3GM_{\rm cl}},$ (6)

where $\sigma_{0}$, $E_{\rm kin}$, and $E_{\rm grav}$ are the 3D velocity
dispersion, kinetic, and gravitational energy. We adopt the value of
virialized clouds as $\alpha_{0}=1$ as the fiducial value. Recent observations
showed that the virial parameters somewhat change depending on the formation
sites of the clouds. (e.g., Roman-Duval et al., 2010; Leroy et al., 2016; Sun
et al., 2018). Thus, we also perform the simulations with $\alpha=2$ only with
$Z=Z_{\odot}$ to clarify the effects of turbulent motions. Here, the
turbulence decays freely from the start of the simulations, and it leads to
the rapid collapse of the cloud compared with the cases with the external
turbulent sources (e.g., Krumholz et al., 2011).

##  3 Results

We first study the effects of the cloud compactness by comparing M6R20Z0A1 and
M6R40Z0A1 in Section 3.1. In Section 3.2, we show the SFEs and properties of
star clusters in all models. In Section 3.3, we discuss the threshold surface
density for the formation of gravitationally bounded clusters. In Section 3.4
and 3.5, we describe the dependence of SFEs on the virial parameters and the
escape velocities of the clouds. In Table 2, we summarize the results obtained
in our simulations. In Appendix E, we compare the results with the different
maximum refinement levels and their impacts on our study.

###  3.1 Star cluster formation in massive cloud

####  3.1.1 compact cloud model

Figure 1:  The structure of a cloud with $(M_{\rm cl},R_{\rm
cl},Z)=(10^{6}~{}M_{\odot},20~{}{\rm pc},Z_{\odot})$ at $t=0.6(0.4~{}t_{\rm
ff}),1.9(1.3~{}t_{\rm ff}),2.9(2~{}t_{\rm ff})$ and $4.4~{}{\rm
Myr}(3~{}t_{\rm ff})$. Each panel shows the surface density of gas, the number
column densities of electron, $\rm H_{2}$ and CO molecules from top to bottom.
Stellar particles are shown as white dots.

Figure 2:  Upper panels: The time evolution of the stellar mass normalized by
the initial cloud one and the bound fraction in the models of $(M_{\rm
cl},R_{\rm cl},Z)=(10^{6}~{}M_{\odot},20~{}{\rm pc},Z_{\odot})$ (left) and
$(10^{6}~{}M_{\odot},40~{}{\rm pc},Z_{\odot})$ (right). Solid and dashed black
lines represent the total stellar masses and the bound fractions. Lower
panels: The masses of neutral atomic hydrogen ($M_{\rm n}$, blue), molecular
hydrogen ($M_{\rm H_{2}}$, orange), and ionized gas ($M_{\rm HII}$, green)
normalized by total gas mass. The labels of (1)-(4) represent the four epochs
as shown in Figure 1 and 3.

We first describe the model of $(M_{\rm cl},R_{\rm cl},Z_{\rm
cl})=(10^{6}M_{\odot},20{\rm pc},Z_{\odot})$ as the fiducial one. In Figure 1,
we show the time evolution of this cloud. In the early phase, the turbulent
motion controls the gas dynamics and induces the formation of the filamentary
structures in which stars form. Yet, most new born stars distributed at the
center of the cloud and are tightly bound at $t\sim 1.3\,t_{\rm ff}$. The star
cluster at the center keeps the compactness until the end of the simulation.
In this model, some remains even at $t=3~{}t_{\rm ff}$.

The second column of Figure 1 shows the electron column density, representing
the spatial distributions of Hii regions. The emissivity of ionizing photons
increases with the total stellar mass, and the Hii regions gradually expand.
At $t=1.3\,t_{\rm ff}$, the Hii regions are localized around the stars. Then,
the entire volume is highly ionized at $t>2~{}t_{\rm ff}$.

The spatial distributions of $\rm H_{2}$ and $\rm CO$ molecules also change
with the star formation (see the third and fourth columns of Figure 1). After
the onset of the star formation, FUV photons dissociate these molecules.
Especially, the molecules are efficiently photodissociated in the low-density
regions. On the other hand, the molecules in the high-density filaments can
survive against the radiative feedback for a longer time because of the dust-
shielding effect (Fukushima et al., 2020b). However, most of the molecules
disappear finally because they are consumed by the star formation.

Figure 2 presents the stellar mass and the bound fraction of the star cluster
as a function of time. Here, we adopt a similar procedure in Shukirgaliyev et
al. (2017) to estimate the bound fraction. We calculate the gravitational
binding energy and the kinetic energy of sink particles, and then remove the
unbound particles. We repeat the above processes until all remained sink
particles are gravitationally bounded. The star formation begins when half of
the free-fall time elapses. The total stellar mass reaches 0.1 of the cloud
mass at $t\sim 1.3\,t_{\rm ff}$. At the time, the bound fraction rapidly
increase. The increase of the bound fraction at ${\rm SFE}\gtrsim 0.1$ is
consistent with the previous study (Shukirgaliyev et al., 2017). At $t\sim
1.5\,t_{\rm ff}$, the bound fraction exceeds 0.9, and it is almost constant
until the end of the simulation. The star formation rate (SFR) also starts to
increase at $t\sim 1.3\,t_{\rm ff}$, and it continues until $t\sim 2\,t_{\rm
ff}$. After that, radiative feedback slows down the star formation, but gas
around the star cluster cannot disperse due to the deep gravitational
potential. Thus, the star formation continues for a long time until $t\sim
3.5\,t_{\rm ff}$, finally resulting in the SFE of 0.7.

The bottom panels of Figure 2 show the time evolution of non-ionized (atomic
and molecular hydrogen), ionized, and $\rm H_{2}$ molecule gas. The hydrogen
molecule abundance rapidly decreases as the star formation proceeds because of
the photodissociation in the low-density regions. It becomes $\lesssim 0.2$
when half of the cloud is converted into stars. At $t\gtrsim 3.5\,t_{\rm ff}$,
all hydrogen molecules are dissociated and more than eighty percent of the gas
is ionized. In this phase, the gravitational collapse is prevented due to the
pressure of the ionized gas, resulting in the quenching of the star formation.

####  3.1.2 Diffuse cloud model

Figure 3:  Same as Figure 1, but for the case with $(M_{\rm cl},R_{\rm
cl},Z)=(10^{6}~{}M_{\odot},40~{}{\rm pc},Z_{\odot})$. The snapshots are at
$t=1.7(0.4~{}t_{\rm ff})$, $5.7(1.3~{}t_{\rm ff})$, $8.5(2~{}t_{\rm ff})$, and
$12.4~{}{\rm Myr}(3~{}t_{\rm ff})$.

To investigate the impacts of the cloud compactness on the star formation, we
simulate a diffuse cloud with $(M_{\rm cl},R_{\rm
cl},Z)=(10^{6}~{}M_{\odot},40~{}{\rm pc},Z_{\odot})$ and compare with the
fiducial model. Figure 3 shows the evolution of this cloud as in Figure 1.
Same as the fiducial model, the turbulent motions drive the filamentary
structures in the early phase. The star formation begins at $t\sim
0.5~{}t_{\rm ff}$. In this case, the stellar distribution is more extended,
not concentrated at the center, unlike the fiducial compact case. (see Fig 3-3
and 4). As the stellar mass increases, the gas is evacuated due to the
radiative feedback. The shallower gravitational potential well can not hold
the hot ionized gas. Therefore, once the cloud is ionized, it is dispersed
($t\sim 3~{}t_{\rm ff}$).

In this diffuse case, the dust column densities in the filaments are not
enough to shield FUV photons. Therefore, the photodissociation fronts rapidly
expand. At $t\sim 1.3~{}t_{\rm ff}$, more than ninety per cent of hydrogen
molecules are dissociated as shown in Figure 3-(2). After $t\sim 2.0~{}t_{\rm
ff}$, all hydrogen molecules disappear, resulting in the quenching of the star
formation.

The stellar mass and the bound fraction are presented in the right panels of
Figure 2. The star formation begins at $t\sim 0.5~{}t_{\rm ff}$ and continues
until $t\sim 2~{}t_{\rm ff}$. The final SFE is 0.18 which is lower than the
fiducial case by a factor of $\sim 4$. Unlike the fiducial case, the SFR does
not change with the time for $t\sim 0.5-2.0~{}t_{\rm ff}$. The bound fraction
is much smaller than 0.1 at $t<1.5~{}t_{\rm ff}$, but it rapidly increases
once the SFE exceeds $\sim 0.15$. At the end of star formation, the bound
fraction becomes $\sim 0.6$ and remains constant until the end of the
simulation. This trend is consistent with the previous case.

####  3.1.3 low-metallicity cases

Figure 4:  Same as Figure 1 but for the case with $(M_{\rm cl},R_{\rm
cl},Z)=(10^{6}~{}M_{\odot},20~{}{\rm pc},10^{-2}Z_{\odot})$. The snapshots are
at $t=0.6(0.4~{}t_{\rm ff})$, $1.9(1.3~{}t_{\rm ff})$, $3(2~{}t)$, and
$4.5~{}{\rm Myr}(3.1~{}t_{\rm ff})$.

As shown in Paper I, the SFE sensitively depends on the metallicity of clouds.
If the metallicity is low, the temperature of ionized gas becomes higher and
the dust-shielding of UV radiation is ineffective, resulting in the lower SFE.
However, the previous study did not take massive compact clouds into account.
Therefore, here we show the results of the low-metallicity cloud with $(M_{\rm
cl},R_{\rm cl},Z)=(10^{6}~{}M_{\odot},20~{}{\rm pc},10^{-2}Z_{\odot})$.

As in the fiducial model, the low-metallicity cloud begins to form stars at
$t\sim 1.3~{}t_{\rm ff}$. Then, the dense massive cluster form at the center,
and most of the stars are bound with $f_{\rm bd}=0.89$. However, the SFE
($\epsilon_{*}=0.28$) is smaller than that of the fiducial model
($\epsilon_{*}=0.7$) due to the stronger feedback.

Figure 4 also shows the column densities of ${\rm H_{2}}$ and ${\rm CO}$
molecules. In this case, we assume that the gas is atomic initially. Also,
molecular formation on dust grains is insufficient to make gas fully
molecular. Thus, the molecules only exist in the high-density regions, and the
maximum column density of $\rm H_{2}$ molecules are 10 times smaller than that
of the case with $Z=Z_{\odot}$. After the onset of the star formation, FUV
photons start to dissociate molecules in the low-density regions. In Paper I,
we showed that FUV photons completely photodissociate molecules in the lower-
mass clouds. The star formation occurs from atomic gas at the low-metallicity
(Krumholz, 2012). On the other hand, the molecules remain in the high density
filaments in the cloud with $M_{\rm cl}=10^{6}M_{\odot}$, as shown in Figure
4. In Paper I, we estimated the optical depth of the filaments for FUV
photons, considering the Jeans unstable filaments as

$\displaystyle\tau_{\rm sf}=4.6\left(\frac{\Sigma}{800~{}M_{\odot}{\rm
pc^{-2}}}\right)\left(\frac{\alpha_{0}}{1}\right)\left(\frac{Z}{10^{-2}Z_{\odot}}\right).$
(7)

The star-forming filaments are optically-thick for FUV photons, and dust
shielding prevents photodissociation even at $Z=10^{-2}Z_{\odot}$. Therefore,
the molecular gas remains until the high-density filaments are completely
destroyed.

###  3.2 Dependence of SFEs and SFRs on cloud properties

####  3.2.1 Star formation efficiencies and lifetimes of clouds

Figure 5:  Upper panel: The SFEs of clouds as a function of surface densities.
Each symbol represents the different cloud masses: $M_{\rm
cl}=10^{5}~{}M_{\odot}$ (open) and $10^{6}~{}M_{\odot}$ (filled), and with the
different radii $R_{\rm cl}=5~{}{\rm pc}$ (square), $10~{}{\rm pc}$
(triangle), $20~{}{\rm pc}$ (circle), $40~{}{\rm pc}$ (diamond), and
$60~{}{\rm pc}$ (hexagon). Each color shows the different metallicities:
$Z=10^{-2}Z_{\odot}$ (green), $10^{-1}Z_{\odot}$ (orange), and $Z_{\odot}$
(blue). The black symbols represent the cases with $\alpha_{0}=2$. Lower
panel: The ratio of the cloud lifetimes to the free-fall times.

By using the total stellar mass $M_{*}$ at the end of the simulations, we
estimate the SFEs as $\epsilon_{*}=M_{*}/M_{\rm cl}$. Figure 5 shows the SFEs
as the function of the cloud surface densities. As in Paper I, the SFE
gradually increases with the surface density (see also, Fall et al., 2010;
Raskutti et al., 2016; Geen et al., 2017; Grudić et al., 2018; Kim et al.,
2018; He et al., 2019). In the case with $Z=Z_{\odot}$, the SFE is $0.08$ at
$\Sigma_{\rm cl}=80~{}M_{\rm\odot}{\rm pc^{-2}}$, which is similar with the
results in the previous works of Kim et al. (2018) and He et al. (2019).
Comparing with the results in Paper I, we simulate conditions with higher
turbulent velocity. The higher turbulent motions make the cloud less bound,
and the star formation rate decreases. In such a case, more gas is evaporated
before it is converted into stars. Thus, the SFE is about half as small as
that of our previous work. We will discuss the impacts of the turbulent motion
in detail at Sec. 3.4.

For the higher surface densities, the SFEs can exceed 0.1. For the clouds with
$M_{\rm cl}=10^{5}~{}M_{\odot}$, the SFE achieves 0.47 at $\Sigma_{\rm
cl}=1300~{}M_{\odot}\,{\rm pc^{-2}}$ (M5R5Z0A1). In the cases with the massive
clouds of $M_{\rm cl}=10^{6}\,M_{\odot}$, the SFEs become higher, and 0.7 for
$800~{}M_{\odot}\,{\rm pc^{-2}}$ (M6R20Z0A1) and 0.78 for
$3200~{}M_{\odot}\,{\rm pc^{-2}}$ (M6R10Z0A1). In these two clouds, the cloud
escape velocities exceed $20~{}{\rm km/s}$ that is higher than the sound speed
of the Hii regions. In such cases, the deep gravitational potential can hold
the ionized gas and allows further star formation (Dale et al., 2012, 2013).

In the massive compact clouds, the radiation pressure plays a main role in
suppressing the star formation instead of the photoionization feedback (e.g.,
Krumholz & Matzner, 2009; Fall et al., 2010). The radiation pressure is due to
the direct light from the radiation sources and dust thermal emission (IR
emission). We can roughly estimate the conditions that radiation pressure
suppresses star formation, assuming the spherical shell accelerated by the
radiation pressure. The radiation pressure of the direct light on the shell is
estimated as

$\displaystyle F_{\rm rad}=\frac{L}{c}=\frac{\epsilon_{*}M_{\rm cl}l_{*}}{c},$
(8)

where $\epsilon_{*}$ is the SFE, and $l_{*}$ is the luminosity per unit mass
defined as $L=l_{*}M_{*}$. We estimate the gravitational force on the shell at
the outer boundary of the cloud as (Kim et al., 2016)

$\displaystyle F_{\rm grav}=\frac{GM_{\rm sh}(M_{*}+M_{\rm sh}/2)}{R_{\rm
cl}^{2}}=\frac{GM_{\rm cl}^{2}}{2R_{\rm
cl}^{2}}\left(1-\epsilon_{*}^{2}\right),$ (9)

where the self-gravity of the shell is included. The shell mass is $M_{\rm
sh}=(1-\epsilon_{*})M_{\rm cl}$. By taking the balance between the radiation
pressure (Equation 8) and the gravitational force (Equation 9), the SFE is
estimated as (Kim et al., 2016; Krumholz et al., 2019),

$\displaystyle\frac{\epsilon_{*}}{1-\epsilon_{*}^{2}}=\frac{\Sigma_{\rm
cl}}{\Sigma_{\rm PR,UV}},$ (10)

where $\Sigma_{\rm cl}$ is the cloud’s surface density, and

$\displaystyle\Sigma_{\rm RP,UV}=\frac{2l_{*}}{\pi Gc}\simeq 3.9\times
10^{3}~{}M_{\odot}{\rm pc^{-2}}.$ (11)

The radiation pressure can evacuate the gas if the SFE exceeds the above. Note
that, the SFE can be larger than that given in Equation (11) because there is
a time lag of a free-fall time from the launching to the gas evacuation.
Equation (10) indicates the radiation pressure from the direct light is
effective only if $\Sigma_{\rm cl}<\Sigma_{\rm RP,UV}$. If $\Sigma_{\rm
cl}>\Sigma_{\rm RP,UV}$, the SFE becomes almost unity. All clouds in this
study satisfy the condition of $\Sigma_{\rm cl}<\Sigma_{\rm RP,UV}$. However,
the surface density of the clouds with $(M_{\rm cl},R_{\rm
cl})=(10^{6}~{}M_{\odot},10~{}{\rm pc})$ is close to this limit. Thus, the
star formation hardly stops, and the SFE is close to 80 per cent. The
radiation pressure from dust thermal emission is only important if the optical
depth of the IR light is larger than unity as

$\displaystyle\tau_{\rm IR}=\rho\kappa_{\rm IR}R_{\rm
cl}=0.5\left(\frac{\Sigma}{3.200~{}M_{\odot}\,{\rm
pc^{-2}}}\right)\left(\frac{\kappa_{\rm IR}}{1~{}{\rm
cm^{2}g^{-1}}}\right)>1.$ (12)

In our simulations, the initial states of clouds are optically thin for IR
photons even at the solar metallicity. Thus, the radiation pressure of IR
photons is unlikely to be effective. However, as the cloud collapses, the
surface density can be larger than $10^{4}~{}M_{\odot}{\rm pc^{-2}}$ locally
as shown in Figure 1. In these high-density regions, the IR radiation pressure
works in suppressing star formation. If an individual star is resolved in
simulations, the IR radiation pressure also plays a role in regulating the
accretion rate onto a proto-star (e.g., Wolfire & Cassinelli, 1987; Krumholz
et al., 2009; Kuiper et al., 2010; Fukushima et al., 2018, 2020a). Our current
simulations do not take this process into account. Note that Skinner &
Ostriker (2015) argued that the IR radiation pressure works only for cases
with the opacity $\kappa_{\rm IR}\gtrsim 15~{}{\rm cm^{2}}\,{\rm g^{-1}}$
which corresponds to the metallicity higher than the solar abundance. Hence,
the radiation pressure from stellar light is a dominant feedback in regulating
star formation.

The SFE decreases at lower-metallicity. In the cases of the clouds with
$\Sigma_{\rm cl}<10^{3}~{}M_{\odot}{\rm pc^{-2}}$, the SFEs for
$Z=10^{-2}Z_{\odot}$ are lower than that for $Z=Z_{\odot}$ by a factor of
$\sim 3$. In this case, the dust shielding of ionizing photons is ineffective
(Fukushima et al., 2020b). Therefore, the propagation of the ionization front
is significantly faster, resulting in the early quenching of the star
formation. On the other hand, if $\Sigma>10^{3}~{}M_{\odot}{\rm pc^{-2}}$, the
SFE is not sensitive to the metallicity. In these compact clouds, the main
feedback is the radiation pressure due to the photon absorption by neutral
hydrogen or interstellar dust. Since the optical depth of the dust for UV
photons is larger than unity even at $Z=10^{-2}Z_{\odot}$, the strength of the
radiation pressure does not change significantly irrespective of the
metallicity.

Table 2: Simulation results

model | $\epsilon_{*}$ | $\epsilon_{\rm ff,0}$ | $\epsilon_{\rm ff}$ | $t_{\rm life}\,[\,{\rm Myr}\,]$ | $t_{\rm dr}\,[\,{\rm Myr}\,]$ | $M_{\rm bd}\,[\,{M_{\odot}}\,]$ | $f_{\rm bd}$ | $r_{\rm h}\,[\,{\rm pc}\,]$ | $t_{\rm rh}\,[\,{\rm Myr}\,]$
---|---|---|---|---|---|---|---|---|---
M5R5Z0A1 | $0.47$ | $0.15$ |  $0.30$ a | $1.71(2.91t_{\rm ff})$ | $1.24(2.12t_{\rm ff})$ | $4.38\times 10^{4}$ | $0.92$ | $0.88$ | $1.09\times 10^{2}$
M5R5Z0A2 | $0.33$ | $0.11$ |  $0.19$ a | $1.69(2.88t_{\rm ff})$ | $1.24(2.12t_{\rm ff})$ | $2.91\times 10^{4}$ | $0.87$ | $1.06$ | $1.22\times 10^{2}$
M5R5Z-1A1 | $0.34$ | $0.15$ | $0.22$ | $1.27(2.17t_{\rm ff})$ | $0.83(1.41t_{\rm ff})$ | $2.77\times 10^{4}$ | $0.82$ | $0.87$ | $8.87\times 10$
M5R5Z-2A1 | $0.29$ | $0.11$ |  $0.21$ a | $1.39(2.38t_{\rm ff})$ | $0.92(1.57t_{\rm ff})$ | $2.38\times 10^{4}$ | $0.83$ | $0.99$ | $1.01\times 10^{2}$
M5R8Z0A1 | $0.32$ | $0.13$ |  $0.24$ a | $2.81(2.37t_{\rm ff})$ | $1.88(1.59t_{\rm ff})$ | $2.72\times 10^{4}$ | $0.86$ | $1.75$ | $2.51\times 10^{2}$
M5R10Z0A1 | $0.23$ | $0.11$ | $0.17$ | $3.22(1.94t_{\rm ff})$ | $1.94(1.17t_{\rm ff})$ | $1.54\times 10^{4}$ | $0.68$ | $2.50$ | $3.43\times 10^{2}$
M5R10Z0A2 | $0.13$ | $0.05$ |  $0.10$ a | $3.86(2.33t_{\rm ff})$ | $2.65(1.6t_{\rm ff})$ | $1.81\times 10^{2}$ | $0.01$ | $0.69$ | $1.07\times 10$
M5R10Z-1A1 | $0.13$ | $0.08$ | $0.14$ | $2.61(1.58t_{\rm ff})$ | $1.37(0.83t_{\rm ff})$ | $6.15\times 10^{2}$ | $0.05$ | $1.63$ | $5.64\times 10$
M5R10Z-2A1 | $0.09$ | $0.05$ | $0.08$ | $3.05(1.84t_{\rm ff})$ | $1.77(1.07t_{\rm ff})$ |  $1.43\times 10^{2}$ b |  $0.02$ c | \- d | \- e
M5R12Z0A1 | $0.16$ | $0.09$ | $0.16$ | $3.70(1.70t_{\rm ff})$ | $2.06(0.94t_{\rm ff})$ | $5.95\times 10^{3}$ | $0.37$ | $2.39$ | $2.23\times 10^{2}$
M5R20Z0A1 | $0.08$ | $0.05$ | $0.08$ | $7.43(1.58t_{\rm ff})$ | $3.89(0.83t_{\rm ff})$ | $5.39\times 10^{2}$ | $0.07$ | $3.51$ | $1.70\times 10^{2}$
M5R20Z0A2 | $0.05$ | $0.03$ | $0.05$ | $7.04(1.50t_{\rm ff})$ | $3.70(0.79t_{\rm ff})$ | $1.09\times 10^{2}$ | $0.02$ | $0.19$ | $1.38$
M5R20Z-1A1 | $0.05$ | $0.04$ | $0.09$ | $5.81(1.24t_{\rm ff})$ | $2.19(0.47t_{\rm ff})$ | $4.84\times 10^{2}$ | $0.10$ | $2.29$ | $8.71\times 10$
M5R20Z-2A1 | $0.03$ | $0.02$ | $0.03$ | $7.21(1.54t_{\rm ff})$ | $3.39(0.72t_{\rm ff})$ |  $1.71\times 10^{2}$ b |  $0.06$ c | \- d | \- e
M6R10Z0A1 | $0.78$ | $0.20$ |  $0.38$ a | $1.91(3.65t_{\rm ff})$ | $1.46(2.79t_{\rm ff})$ | $7.51\times 10^{5}$ | $0.96$ | $1.28$ | $6.14\times 10^{2}$
M6R10Z0A2 | $0.53$ | $0.14$ |  $0.18$ a | $1.92(3.66t_{\rm ff})$ | $1.49(2.84t_{\rm ff})$ | $4.87\times 10^{5}$ | $0.92$ | $1.69$ | $7.75\times 10^{2}$
M6R10Z-1A1 | $0.84$ | $0.29$ |  $0.46$ a | $1.42(2.72t_{\rm ff})$ | $0.95(1.82t_{\rm ff})$ | $8.19\times 10^{5}$ | $0.97$ | $0.69$ | $2.52\times 10^{2}$
M6R10Z-2A1 | $0.74$ | $0.28$ |  $0.40$ a | $1.33(2.54t_{\rm ff})$ | $0.86(1.63t_{\rm ff})$ | $7.32\times 10^{5}$ | $0.99$ | $0.66$ | $2.27\times 10^{2}$
M6R175Z-2A1 | $0.41$ | $0.18$ | $0.29$ | $2.56(2.11t_{\rm ff})$ | $1.54(1.27t_{\rm ff})$ | $3.82\times 10^{5}$ | $0.94$ | $1.90$ | $8.36\times 10^{2}$
M6R20Z0A1 | $0.70$ | $0.20$ |  $0.39$ a | $4.82(3.25t_{\rm ff})$ | $3.50(2.36t_{\rm ff})$ | $6.70\times 10^{5}$ | $0.96$ | $1.73$ | $9.17\times 10^{2}$
M6R20Z0A2 | $0.37$ | $0.09$ |  $0.12$ a | $5.88(3.97t_{\rm ff})$ | $4.63(3.12t_{\rm ff})$ | $3.53\times 10^{5}$ | $0.94$ | $2.33$ | $1.10\times 10^{3}$
M6R20Z-1A1 | $0.42$ | $0.16$ |  $0.28$ a | $3.70(2.49t_{\rm ff})$ | $2.48(1.67t_{\rm ff})$ | $4.05\times 10^{5}$ | $0.96$ | $1.98$ | $9.09\times 10^{2}$
M6R20Z-2A1 | $0.28$ | $0.13$ | $0.19$ | $3.15(2.13t_{\rm ff})$ | $1.94(1.31t_{\rm ff})$ | $2.50\times 10^{5}$ | $0.88$ | $3.53$ | $1.78\times 10^{3}$
M6R25Z0A1 | $0.62$ | $0.21$ |  $0.40$ a | $5.73(2.77t_{\rm ff})$ | $3.90(1.88t_{\rm ff})$ | $6.01\times 10^{5}$ | $0.97$ | $2.03$ | $1.12\times 10^{3}$
M6R25Z-2A1 | $0.16$ | $0.08$ | $0.13$ | $3.85(1.86t_{\rm ff})$ | $2.27(1.10t_{\rm ff})$ | $5.54\times 10^{4}$ | $0.34$ | $5.74$ | $1.99\times 10^{3}$
M6R30Z0A1 | $0.53$ | $0.20$ |  $0.35$ a | $6.77(2.49t_{\rm ff})$ | $4.42(1.62t_{\rm ff})$ | $5.15\times 10^{5}$ | $0.97$ | $2.11$ | $1.11\times 10^{3}$
M6R30Z-2A1 | $0.12$ | $0.06$ | $0.09$ | $5.21(1.92t_{\rm ff})$ | $3.16(1.16t_{\rm ff})$ | $2.40\times 10^{3}$ | $0.02$ | $0.87$ | $3.51\times 10$
M6R32.5Z0A1 | $0.41$ | $0.16$ |  $0.24$ a | $7.23(2.35t_{\rm ff})$ | $4.68(1.53t_{\rm ff})$ | $3.90\times 10^{5}$ | $0.96$ | $2.38$ | $1.18\times 10^{3}$
M6R35Z0A1 | $0.27$ | $0.12$ | $0.18$ | $7.23(2.11t_{\rm ff})$ | $4.52(1.32t_{\rm ff})$ | $2.39\times 10^{5}$ | $0.88$ | $5.23$ | $3.14\times 10^{3}$
M6R35Z-2A1 | $0.09$ | $0.05$ | $0.10$ | $5.25(1.53t_{\rm ff})$ | $2.69(0.79t_{\rm ff})$ |  $1.65\times 10^{3}$ b |  $0.02$ c | \- d | \- e
M6R40Z0A1 | $0.18$ | $0.09$ | $0.13$ | $8.17(1.95t_{\rm ff})$ | $5.08(1.21t_{\rm ff})$ | $7.77\times 10^{4}$ | $0.44$ | $10.0$ | $5.26\times 10^{3}$
M6R40Z0A2 | $0.08$ | $0.04$ | $0.06$ | $8.00(1.91t_{\rm ff})$ | $5.30(1.26t_{\rm ff})$ | $7.68\times 10^{2}$ | $0.01$ | $4.27$ | $2.57\times 10^{2}$
M6R40Z-1A1 | $0.08$ | $0.05$ | $0.09$ | $6.48(1.54t_{\rm ff})$ | $3.66(0.87t_{\rm ff})$ | $1.64\times 10^{3}$ | $0.02$ | $3.97$ | $2.99\times 10^{2}$
M6R40Z-2A1 | $0.07$ | $0.04$ | $0.08$ | $6.49(1.55t_{\rm ff})$ | $3.39(0.81t_{\rm ff})$ |  $1.75\times 10^{3}$ b |  $0.03$ c | \- d | \- e
M6R60Z0A1 | $0.08$ | $0.05$ | $0.09$ | $11.2(1.46t_{\rm ff})$ | $5.87(0.76t_{\rm ff})$ |  $1.41\times 10^{3}$ b |  $0.02$ c | \- d | \- e
M6R60Z0A2 | $0.04$ | $0.02$ | $0.04$ | $12.6(1.64t_{\rm ff})$ | $7.80(1.01t_{\rm ff})$ |  $9.21\times 10^{2}$ b |  $0.02$ c | \- d | \- e
M6R60Z-1A1 | $0.04$ | $0.03$ | $0.06$ | $9.72(1.26t_{\rm ff})$ | $4.67(0.61t_{\rm ff})$ |  $1.56\times 10^{3}$ b |  $0.04$ c | \- d | \- e
M6R60Z-2A1 | $0.03$ | $0.02$ | $0.03$ | $11.9(1.55t_{\rm ff})$ | $6.38(0.83t_{\rm ff})$ |  $1.61\times 10^{3}$ b |  $0.05$ c | \- d | \- e

Notes. Column 1: model names, Column 2: star formation efficiency, Column 3:
star formation rate per free fall time, Column 4: specific star formation rate
per free fall time in the duration time of star formation $t_{\rm dr}$,
Column5: lifetime of clouds, Column 6: duration time of star formation, Column
7: gravitationally bounded mass, Column 8: gravitational bound fraction of
star cluster, Column 9: half mass radius, Column 10: half mass relaxation
time.

a We estimate the parameters $\epsilon_{\rm ff}$ from the star formation rates
in the period between the start of the star formation and $1~{}t_{\rm ff}$
later.

b,c All star particles are not gravitationally bounded, and we obtain the
gravitational bounded mass and the bound fractions from the maximum mass of
sink particles.

d,e The half mass radius and the relaxation time are not evaluated if all
stellar particles are gravitationally unbounded.

The bottom panel of Figure 5 shows the cloud lifetimes. The cold gas and
molecules inside it are almost dispersed by radiative feedback at the end of
the star formation. Therefore, the end of the cloud lifetimes is almost the
same as the end of the star formation. Here, we define the cloud lifetimes as
the period from the starting time of the simulations to the time when the
total stellar mass reaches 95 percent of the final stellar mass. The lifetimes
increase slightly in the higher surface density clouds. Also, we estimate the
duration time of the star formation defined as the time required for the total
stellar mass to increase from 5 percent to 95 percent of the final stellar
mass (Table 2). The star formation mainly occurs after one free-fall time in
all cases and continues for one free-fall time in the clouds with $\Sigma_{\rm
cl}<10^{3}~{}M_{\odot}{\rm pc^{-2}}$. In the cases with $\Sigma_{\rm
cl}=3200~{}M_{\odot}{\rm pc^{-2}}$, the clouds are not dispersed, and the star
formation continues until $t\sim 3.3t_{\rm ff}$. In the low-metallicity
environments, the lifetimes are shorter than that of the solar metallicity.

Figure 6:  The star formation parameters $\epsilon_{\rm ff,0}$ and
$\epsilon_{\rm ff}$ as the function of the surface densities. The upper panel
shows the fraction of the cloud mass transformed into stars per free-fall time
averaged in the lifetimes of clouds $(t_{\rm life})$ (Eq. 13). The lower panel
represents the same fraction but averaged in the duration time of the star
formation $(t_{\rm dr})$ (Eq. 14). The styles of symbols are the same as in
Figure 5.

The star formation rate (SFR) is also related to the initial conditions of
clouds. In observations, the SFR is speculated by measuring properties of star
clusters and cloud life times for various clouds. By introducing a conversion
efficiency from the gas into stars $\epsilon_{\rm ff,0}$, the SFR is expressed
as

$\displaystyle\dot{M}_{*,0}$ $\displaystyle=\frac{M_{*}}{t_{\rm
life}}=\epsilon_{\rm ff,0}\frac{M_{\rm cl}}{t_{\rm ff}},$ (13)

where $t_{\rm life}$ is the cloud lifetime and $t_{\rm
ff}=\sqrt{3\pi/(32G\rho)}$ is the free-fall time. This estimate corresponds to
take into account the various evolution stages of clouds over the lifetime.
However, as shown in Figure 2, most stars form within a shorter time scale
than the cloud lifetime. Therefore, we here evaluate the SFR by using the
duration time of the star formation ($t_{\rm dr}$):

$\displaystyle\dot{M}_{*}=\frac{M_{*}}{t_{\rm dr}}=\epsilon_{\rm
ff}\frac{M_{\rm cl}}{t_{\rm ff}}.$ (14)

The above $\epsilon_{\rm ff}$ is related to the star formation efficiency as

$\displaystyle\epsilon_{\rm ff}=\epsilon_{*}\frac{t_{\rm ff}}{t_{\rm dr}}.$
(15)

As shown in Figure 2, the duration time of the star formation is much longer
when the stellar core formation occurs. As discussed in Appendix E, our
current simulations may somewhat overestimate the impact of the radiation
pressure. However, even in such a case, the radiation pressure cannot evacuate
the gas and allows the gas accretion onto the stellar core until its mass
becomes massive enough. Consequently, the overestimated radiation pressure
does not alter the SFE, while it makes the duration time of star formation
longer. Here, we underestimate the SFR of the main phase of the star
formation. Thus, we recalculate $\epsilon_{\rm ff}$ in the periods from the
starting time of the star formation to $1.0~{}t_{\rm ff}$ later if the
duration time is longer than $1.5~{}t_{\rm ff}$. Figure 6 presents
$\epsilon_{\rm ff,0}$ and $\epsilon_{\rm ff}$ in each cloud. Because of
$t_{\rm dr}$ shorter than $t_{\rm life}$, $\epsilon_{\rm ff}$ is higher than
$\epsilon_{\rm ff,0}$ by a factor of $\sim 2$. Both panels show that the SFRs
increase with the higher surface densities. Besides, the SFRs are smaller as
the metallicity decreases. In Section 3.3, we further discuss the increase of
the SFRs in the compact clouds.

###  3.3 Rapid increase of SFE

As in Paper I, if the clouds are disrupted due to the photoionization
feedback, the SFE is expressed as the power-law function of the surface
density $(\propto\Sigma^{-1/2})$. However, in the case of massive compact
clouds, the SFE is unlikely to obey the simple power-law function because the
photoionization does not suppress the star formation. We perform the
additional simulations of the $10^{6}~{}M_{\odot}$ clouds at $Z=Z_{\odot}$
with $R_{\rm cl}=25$, $30$, $32.5$, and $35~{}{\rm pc}$ to investigate the
relation between the SFE and the surface density. We further perform the
simulations of the $10^{5}~{}M_{\odot}$ clouds at $Z=Z_{\odot}$ with $R_{\rm
cl}=8$, $12~{}{\rm pc}$, and the $10^{6}~{}M_{\odot}$ clouds at
$Z=10^{-2}~{}Z_{\odot}$ with $R_{\rm cl}=17.5$, $25$, $30$, $35~{}{\rm pc}$ to
consider the dependency on the cloud masses and metallicity.

Figure 7:  The time evolution of the stellar mass in the case with $(M_{\rm
cl},\alpha_{0})=(10^{6}~{}M_{\odot},1)$. Each line represents the different
initial cloud radius: $R_{\rm cl}=20$, $25$, $30$, $32.5$, $35$, and
$40~{}{\rm pc}$ from top to bottom.  Figure 8:  The SFEs in the cases with
$(M_{\rm cl},Z)=(10^{6}~{}M_{\odot},Z_{\odot})$ (blue filled circles),
$(10^{5}~{}M_{\odot},Z_{\odot})$ (blue open circles), and
$(10^{6}~{}M_{\odot},10^{-2}~{}Z_{\odot})$ (green filled circles) as the
function of the surface densities. The dashed lines show the power law
functions as $\epsilon_{*}\propto\Sigma_{\rm cl}^{1/2}$ with different
normalization factors related to physical properties of clouds.

Figure 7 shows the time evolutions of the total stellar masses in the clouds
of $10^{6}~{}M_{\odot}$ with the different initial radii. All models show the
similar stellar masses until $t\sim 1.3~{}t_{\rm ff}$. Then, the growth curves
of the stellar masses look quite different. In the cases with $R_{\rm
cl}\gtrsim 35~{}{\rm pc}$, the stellar masses increase slowly with almost
constant SFRs until $t\sim 2~{}t_{\rm ff}$. Then the star formation is
quenched and the stellar masses become constant. On the other hand, the SFRs
significantly increase at $t\sim 1.3~{}t_{\rm ff}$ in the more compact clouds.
Besides, the duration times of the star formation are longer, resulting in the
SFEs higher than 0.4. Note that, although the surface density of M6R30Z0A1 is
only 1.3 times larger than that of M6R35Z0A1, the SFE is double.

In Figure 8, we present the dependence of the SFEs on the surface densities.
We find that the SFEs can be fit with the power law function of
$\epsilon_{*}\propto\Sigma_{\rm cl}^{1/2}$ at $\Sigma_{\rm cl}\lesssim
200~{}M_{\odot}{\rm pc^{-2}}$. The SFE jumps at $\Sigma\sim
300~{}M_{\odot}{\rm pc^{-2}}$ in the cases with $(M_{\rm
cl},Z)=(10^{6}~{}M_{\odot},1)$, and achieves $\sim 0.7$ at $\sim
800~{}M_{\odot}{\rm pc^{-2}}$. In the other cases, the SFEs are gradually
apart from the power-law lines.

Figure 9:  The ratio of gas pressure $(f_{\rm p})$ and radiation one ($f_{\rm
r}$) to the gravitational force $(f_{\rm gr})$ on the plane crossing the
center of mass of the star clusters. Each panel shows the number density of
hydrogen $(n_{\rm H})$, $f_{\rm p}/f_{\rm gr}$, and $f_{\rm r}/f_{\rm gr}$
from top to bottom. Each figure shows the cases with models of M6R20Z0A1
(left) and M6R40Z0A1 (right). The snapshots are taken from $t=1.3~{}t_{\rm
ff}$, $2~{}t_{\rm ff}$, and $3~{}t_{\rm ff}$. The white contours in middle and
bottom panels indicate the distributions of number densities of gas. White
dots in top panels represent the potions of sink particles.

At the time when the SFE starts to jump, the stars and the gas concentrate at
the centers of the clouds, resulting in deeper gravitational potential. This
can make the bound fraction of the star cluster larger as shown in Figure 2,
resulting in the formation of stellar cores. In this phase, surrounding gas
accretes to the central regions and form stars efficiently. The rapid increase
of the SFR starts when the SFE achieves $\sim 0.1$ regardless of the cloud
model, and it is consistent with the previous study (e.g., Shukirgaliyev et
al., 2017). In the cases of the diffuse clouds, the radiative feedback can
evacuate the gas from the central region and suppress the formation of the
stellar cores.

The conditions of the stellar core formation are determined by the competition
between the gravitational force and the radiative feedback. In Figure 9, we
show the spatial distributions of the ratio of thermal pressure $(f_{\rm p})$
and radiation pressure $(f_{\rm r})$ to the gravitational force $(f_{\rm gr})$
in the models of M6R20Z0A1 and M6R40Z0A1. The gravitational force includes the
contributions from gas and stars. In the case of M6R20Z0A1, the core forms at
$t\sim 1.3~{}t_{\rm ff}$. At this time, the Hii regions start to expand, and
thermal pressure is dominant in the low-density regions. However, the
gravitational force overcomes thermal pressure in the high-density regions
around the core. Also, the radiation pressure is weaker than the gravitational
force until $t\sim 2~{}t_{\rm ff}$. Thus, the high-density gas remains in the
central region, and the star formation proceeds rapidly as shown in Figure 7.
Then, the radiation pressure increases with the stellar mass and evacuates the
gas at $t\sim 3~{}t_{\rm ff}$, making the hole structure as shown in the
figure. In the case of M6R40Z0A1, the stars are concentrated at
$t=1.3~{}t_{\rm ff}$ as in M6R20Z0A1. However, the force from the thermal
pressure is larger than the gravitational force even in the higher-density
region. Therefore, the gas around the star cluster is dispersed rapidly. As a
result, no stellar core is formed, and stars are widely distributed. To
summarize the above discussion, the SFE jump condition is whether the
gravitational force overcomes the thermal pressure in the high-density region.
We further discuss the condition with the semi-analytical model in Section
4.1.

Figure 10:  Same as Figure 5 but only for the cases with $M_{\rm
cl}=10^{6}~{}M_{\odot}$. The styles of the symbols are the same as in Figure
8. The dashed lines show the fitting function of Eq.(16). The fitting
parameters are tabulated in Table 3.

As discussed in Section 3.2.1, the SFRs is characterized by the parameter
$\epsilon_{\rm ff,0}$ or $\epsilon_{\rm ff}$. Figure 10 shows the dependencies
of $\epsilon_{\rm ff,0}$ and $\epsilon_{\rm ff}$ on the surface densities.
Both $\epsilon_{\rm ff,0}$ and $\epsilon_{\rm ff}$ are like step functions
transiting at the specific surface densities. We make a fitting function as

$\displaystyle\epsilon_{\rm ff}=\epsilon_{\rm ff,1}+\frac{\epsilon_{\rm
ff,2}-\epsilon_{\rm ff,1}}{1+e^{-f(\Sigma_{\rm cl})}},$ (16)

where $\epsilon_{\rm ff,1}$ and $\epsilon_{\rm ff,2}$ are the fitting
parameters for the low and high surface densities. We use the values at the
most compact and most diffuse cases for $\epsilon_{\rm ff,1}$ and
$\epsilon_{\rm ff,2}$. In Equation (16), $f(\Sigma_{\rm cl})$ is defined as

$\displaystyle f(\Sigma_{\rm cl})=a\frac{\Sigma_{\rm cl}-\Sigma_{\rm
th}}{\Sigma_{\rm th}},$ (17)

where we use $a=5$ as the gain of the sigmoid function. In Equation (17),
$\Sigma_{\rm th}$ is the threshold surface density of the stellar core
formation. We further consider the dependence of this threshold surface
density on the cloud mass and metallicities in Section 4.1, and the details of
the fitting parameters are summarized in Table 3.

###  3.4 Dependence on virial parameter

Figure 11:  The time evolution of the stellar mass in the case with $(M_{\rm
cl},R_{\rm cl},Z)=(10^{6}\,M_{\odot},20\,{\rm pc},Z_{\odot})$. Blue and black
lines represent the cases with $\alpha_{0}=1$ and $\alpha_{0}=2$,
respectively.

Turbulent motions prevent the gravitational collapse of clouds. Therefore, the
SFRs are likely to decrease for a larger virial parameter $\alpha_{0}$. In
Figure 11, we show the evolution of the total stellar mass with the different
virial parameters $\alpha_{0}=1$ and $2$ in the clouds with $(M_{\rm
cl},R_{\rm cl},Z)=(10^{6}~{}M_{\odot},20~{}{\rm pc},Z_{\odot})$. In the model
of $\alpha_{0}=1$, the star formation rate rapidly increases at $t\simeq
t_{\rm ff}$ via the core formation. In the case of $\alpha_{0}=2$, the kinetic
energy is comparable to the gravitational one. Accordingly, the turbulent
motions can prevent the cloud collapse, resulting in the lower SFR.

The slow star formation gives time for the expansion of Hii regions, which
disrupts the star-forming clumps. Besides, the shear motions of turbulence
decrease local gas density. In the lower density regions, the Hii regions
expand rapidly. Therefore, the SFEs are smaller for the higher virial
parameter (Kim et al., 2021). Figure 12 shows the SFEs as the function of
surface densities with $\alpha_{0}=1$ and $2$. The SFEs of $\alpha_{0}=2$ are
lower than that of $\alpha_{0}=1$ by a factor of $\sim 1.5-2.0$.

Figure 12:  Same as Figure 5 but for the cases of $(M_{\rm cl},R_{\rm
cl},Z)=(10^{6}\,M_{\odot},20\,{\rm pc},Z_{\odot})$ with $\alpha_{0}=1$ (blue)
and $2$ (black).

We also investigate the cloud lifetimes for the different virial parameters as
shown in Figure 12. At $\Sigma_{\rm cl}\lesssim 200~{}M_{\odot}{\rm pc^{-2}}$,
the cloud lifetime does not change with the virial parameter. At
$200~{}M_{\odot}{\rm pc^{-2}}<\Sigma_{\rm cl}\lesssim 10^{3}~{}M_{\odot}{\rm
pc^{-2}}$, the lifetime of $\alpha=2$ becomes longer than that of $\alpha=1$.
In such a case, the turbulent motions delay the star formation, but the
radiative feedback cannot easily blow away gas from the massive clouds.
Therefore, the star formation continues until most of the gas is exhausted,
resulting in the longer lifetime. At the higher surface density, the cloud
lifetimes are slightly different. In these cases, radiative feedback prevents
gas from converting to stars around the stellar cores. This process mainly
determines the duration time of the star formation.

The turbulent motion also has impacts on the properties of star clusters like
the bound fraction. For example, the bound fraction of M6R40Z0A2 is much lower
than the same mass and radius cloud with a lower virial parameter $\alpha=1$.
In this model, the SFE is lower than $\sim 0.1$ that is the critical value to
induce the formation of a high-density stellar core. Once the surface density
of the clouds exceeds $800\,M_{\odot}{\rm pc^{-2}}$, the SFEs are much higher
than 0.1 even with $\alpha=2$. Therefore, the star clusters with the high
bound fractions can form irrespective of the virial parameter.

###  3.5 Dependency on escape velocity

Previous studies suggested that the escape velocity from a cloud was a key
parameter to understand the impacts of the photoionization feedback on cloud
disruption and star formation. (Dale et al., 2012; Bressert et al., 2012).
According to their simulation results, the photoionization feedback is
ineffective if the escape velocity is higher than the sound speed of ionized
gas, $\sim 10~{}\rm km~{}s^{-1}$. On the other hand, our simulations show that
the SFEs mainly depends on the initial surface densities of the clouds, as
shown in Figure 5 and 12. Although the escape velocity of the clouds with
$(M_{\rm cl},R_{\rm cl})=(10^{6}~{}M_{\odot},40~{}{\rm pc})$ is $15~{}{\rm
km/s}$, the clouds are dispersed rapidly. The SFEs are lower than $0.1$ except
for M6R40Z0A1 of which the value is $0.18$. Thus, the escape velocity is not
only the parameter to understand the suppression of star formation. In section
4, we will further discuss this by developing a semi-analytical model.

###  3.6 Properties of star clusters

Figure 13:  The gravitationally bounded masses and bound fractions of star
clusters at the cloud lifetimes $t=t_{\rm life}$ as the function of the
surface densities. The styles of symbols are the same as Figure 5.

A significant fraction of stars can escape from their natal cluster. Therefore
the final mass of the cluster becomes lower than the initial state. Figure 13
shows the masses of bounded stars at the end of the cloud lifetimes. At
$\Sigma_{\rm cl}\gtrsim 800\,M_{\odot}{\rm pc^{-2}}$, the young massive
clusters (YMCs, $M_{\rm bd}>10^{4}\,M_{\odot}$) form in all cases. On the
other hand, the cluster mass reaches $10^{4}\,M_{\odot}$ only in the two
models of $Z=Z_{\odot}$ (M5R10Z0A1 and M6R40Z0A1) at
$\Sigma<800\,M_{\odot}{\rm pc^{-2}}$. YMCs do not form from the low-
metallicity and low surface density clouds. In particular, the total stellar
mass exceeds $10^{4}\,M_{\odot}$ in the models of M6R40Z-1A1 and M6R40Z-2A1,
but most of the stars are dispersed.

The bound fractions $(f_{\rm bd})$ are related to the compactness of clouds.
Figure 13 shows the bound fractions at the end of the cloud lifetime. At the
solar metallicity, the bound fraction increases from 0.05 to $>0.9$ with the
surface density in the range of $\Sigma_{\rm cl}=80-800~{}M_{\odot}{\rm
pc^{-2}}$. In the cases of the low-metallicty clouds, the bound fractions are
lower than that of the solar metallicity at $\Sigma_{\rm cl}\leqq
320~{}M_{\odot}{\rm pc^{-2}}$, but become almost similar ($\gtrsim 0.9$) if
$\Sigma\geqq 800~{}M_{\odot}{\rm pc^{-2}}$.

Figure 14:  The bound fractions of star clusters at the cloud lifetime
$t=t_{\rm life}$ as the function of the star formation efficiency. The styles
of symbols are the same as Figure 5. The dashed line is the fitting function
given as Equation (18).

Figure 14 presents the bound fraction as the function of the SFEs. In the
cases with the SFEs lower than 0.1, the star clusters can not keep holding
their own stars, resulting in the disperse of them and the low bound
fractions. Once the SFE exceeds 0.1, the bound fraction abruptly increases and
achieves $\gtrsim 0.9$ at ${\rm SFE}\gtrsim 0.3$. These results are consistent
with previous numerical simulations (e.g., Li et al., 2019; Grudić et al.,
2020). We make a fitting function as

$\displaystyle f_{\rm bd}=a_{1}10^{-|\log_{10}(\epsilon_{*})|^{a_{2}}}.$ (18)

The best fit parameters are $(a_{1},a_{2})=(3.9,0.98)$. This fitting formula
is steeper with decreasing SFEs than the analytical model of Li et al. (2019)
that assumed the Maxwellian velocity distribution of stars. The bound fraction
at $\epsilon_{*}\lesssim 0.1$ scatters significantly. This trend was also
reported in Grudić et al. (2020). Note that, our estimate is likely to be
inaccurate if $\epsilon_{*}\ll 0.1$ because a sink particle is modelling
bounded several stars. In this case, the actual bound fraction can be smaller
than our results.

Baumgardt & Kroupa (2007) showed that the relation between the SFE and the
bound fraction depends on the timescale of gas removal. In their results, the
SFE is needed to be higher than 0.33 to form the bounded cluster if the gas is
removed instantaneously. The critical SFE decreases to 0.15 if the removal
timescale is similar to the crossing time of the star cluster. In our
simulations, the removal timescale is comparable to the free-fall time of the
clouds, and the critical SFEs roughly match their results.

Figure 15:  The evolution of half mass radius in the models of $(M_{\rm
cl},R_{\rm cl})=(10^{6}\,M_{\odot},20\,{\rm pc})$ (top) and $(M_{\rm
cl},R_{\rm cl})=(10^{6}\,M_{\odot},10\,{\rm pc})$ (botom). Each line shows the
different metallicity: $Z=Z_{\odot}$ (blue), $10^{-1}Z_{\odot}$ (orange), and
$10^{-2}Z_{\odot}$ (green). The black lines represent the cases for
$\alpha_{\rm vir}=2$ and $Z=Z_{\odot}$.

Figure 15 shows the evolution of half mass radii in the cases with $(M_{\rm
cl},R_{\rm cl})=(10^{6}\,M_{\odot},20\,{\rm pc})$ and
$(10^{6}\,M_{\odot},10\,{\rm pc})$. As shown in Paper I, the radius of a star
cluster rapidly increases in the cases of the low-surface density clouds.
Whereas the bound fraction of the star clusters in the massive compact clouds
is high ($>0.9$) and the radii of them do not increase as shown in the figure.
At $t\lesssim 1.0t_{\rm ff}$, the spatial distributions of stars are
stochastically determined due to the turbulent motions, which induces the
fluctuation of the radii. Once the bound fraction exceeds $\sim 0.1$, the
radii become almost constant and result in the values of $\sim 0.05-0.2$ times
the initial cloud ones.

Figure 16:  The relaxation time $t_{\rm rh}$ of star cluster at $t=t_{\rm
life}$ as the function of the surface densities. The relaxation time $t_{\rm
rh}$ is given as Eq. (19). The styles of symbols are the same as Figure 5.

The structure of a bound star cluster would change via 2-body interaction
between stars. The time-scale is characterized by the relaxation time given as
(Spitzer, 1987; Kumamoto et al., 2019)

$\displaystyle t_{\rm rh}=0.0477\frac{N}{(G\rho_{\rm hm})^{1/2}\log(0.4N)},$
(19)

where $\rho_{\rm hm}$ is the stellar mass density inside the half-mass radius
$r_{\rm h}$. The number of stars $N$ is obtained as $N=M_{\rm bd}/\langle
m\rangle$ where $\langle m\rangle$ is the mean stellar mass. In this work, we
assume the Chabrier IMF $(0.1~{}M_{\odot}\leqq M_{*}\leqq 150~{}M_{\odot})$ to
sink particles. Therefore, the mean stellar mass is estimated as $\langle
m\rangle=0.75~{}M_{\odot}$. Figure 16 shows the relaxation time of each
models. In the cases of $(M_{\rm
cl},Z,\alpha_{0})=(10^{5}~{}M_{\odot},Z_{\odot},2)$ and $(M_{\rm
cl},Z)=(10^{5}~{}M_{\odot},10^{-1}Z_{\odot})$, the relaxation times are
shorter than $10^{2}~{}{\rm Myr}$. As shown in Fujii & Portegies Zwart (2014),
the star cluster can induce the core-collapse within $\sim 0.2\times~{}\rm
t_{\rm rh}$. Due to the core collapse, the stellar density at the center
increases, which induces the formation of stellar-binary systems. In
particular, if the core collapse occurs within a few Myr, massive stars can
migrate to the center and form the binary before the end of their lives. Such
a star cluster can be a potential formation site of the binaries that induce
the gravitational-wave events (Abbott et al., 2016). In our simulations, the
star clusters formed in the clouds with $M=10^{5}~{}M_{\odot}$ and
$\alpha_{0}=2$ or $Z=10^{-1}Z_{\odot}$ can be the formation sites of the
gravitational wave sources.

Note that, since the sink particle masses are much higher than the average
stellar mass of the IMF, the relaxation time of these sink particles might be
underestimated. For further discussions on this problem, higher resolution
simulations are awaited.

##  4 Analytical arguments

###  4.1 Star formation efficiencies

As shown in Section 3, the SFEs mainly depend on the surface densities
($\Sigma$). In Paper I, we constructed the semi-analytical model to describe
the dependence of the SFEs on the initial surface density and metallicity. We
assumed that the duration time of the star formation is equal to the
propagation time of the expanding shell outside the cloud. Then, the total
stellar mass $M_{*}=\epsilon_{*}M_{\rm cl}$ is calculated by multiplying
$t_{\rm HII}$ with the star formation rate $\dot{M}_{*}$ in the clouds.
However, this model is not reasonable for massive-compact clouds where the
photoionization feedback alone cannot quench the star formation. Therefore, we
develop a new semi-analytical model taking the cases of massive-compact clouds
into account.

Here, we use a expanding shell model as in Krumholz & Matzner (2009); Kim et
al. (2016). We assume that a uniform density sphere whose mass and radius are
$M_{\rm cl}$ and $R_{\rm cl}$. The mass of the expanding shell is $M_{\rm
sh}=M_{\rm cl}(1-\epsilon_{*})(r_{\rm sh}/R_{\rm cl})^{3}$, where $r_{\rm sh}$
and $\epsilon_{*}$ are the radius and the SFE. Here, we assume that the SFE is
constant in the entire volume. The radiation pressure and gravitational force
on the shell are estimated from equation (8) and (9). The equation of motion
of the expanding shell is represented as (Krumholz & Matzner, 2009; Kim et
al., 2016)

$\displaystyle\frac{d}{dt}\left(M_{\rm sh}\dot{r}_{\rm sh}\right)=4\pi r_{\rm
sh}^{2}\rho_{\rm i}c_{\rm i}^{2}+\frac{L_{*}}{c}-\frac{GM_{\rm
sh}(M_{*}+M_{\rm sh}/2)}{r_{\rm sh}^{2}},$ (20)

where $\rho_{\rm i}$ and $c_{\rm i}$ are the density and the sound speed in
Hii regions. The second and third terms are radiation force and gravity. We
assume that the shell is optically thick for direct light and thin for thermal
radiation from dust. The first term on the right-hand side is the thermal
pressure on the shell from the Hii regions. The density of the Hii regions is
determined by the balance between the ionization and recombination as

$\displaystyle n_{\rm i}=\left(\frac{\rho_{\rm i}c_{\rm i}^{2}}{k_{\rm
B}T_{\rm i}}\right)=\left(\frac{3f_{\rm ion}S_{\rm ion}}{4\pi r_{\rm
sh}^{3}\alpha_{\rm B}}\right)^{1/2},$ (21)

where $T_{\rm i}$, $S_{\rm ion}$ and $\alpha_{\rm B}$ are the temperature of
ionized gas, the emissivity of ionizing photons and the recombination rate
coefficient $\alpha_{\rm B}=2.6\times 10^{-13}(T_{\rm i}/10^{4}\,{\rm
K})^{-0.8}\,{\rm cm^{3}s^{-1}}$ (Osterbrock, 1989). In the Hii regions, dust
grains absorb a part of ionizing photons. We adopt the absorption rate by
hydrogen atoms as $f_{\rm ion}=1-0.27\,(Z/Z_{\odot})$ (McKee & Williams, 1997;
Krumholz & Matzner, 2009).

In this model, we assume that the gas converts into stars after the shell
passes with the conversion rate $\epsilon_{*,0}$. Here, we can regard
$\epsilon_{*,0}$ as the SFE $\epsilon_{*}$. Therefore, we simply use the
notation $\epsilon_{*}$ for the conversion rate hereafter. The mass of stars
formed in the shell is given as $M_{*}=\epsilon_{*}/(1-\epsilon_{*})M_{\rm
sh}$. The luminosity and the photon emissivity are estimated as

$\displaystyle L_{*}=\frac{\epsilon_{*}}{1-\epsilon_{*}}l_{*}M_{\rm sh},$ (22)
$\displaystyle S_{\rm ion}=\frac{\epsilon_{*}}{1-\epsilon_{*}}s_{*}M_{\rm
sh},$ (23)

where the luminosity and the emissivity per unit mass are given as
$l_{*}=1.3\times 10^{3}\,L_{\odot}\,{M_{\odot}}^{-1}$ and $s_{*}=7.5\times
10^{46}\,{\rm s^{-1}}\,{M_{\odot}^{-1}}$.

We introduce a characteristic radius $r_{\rm ch,r}$ at where the radiation
force is equivalent to the thermal pressure of Hii regions , using equations
(20), (21), (22) and (23):

$\displaystyle r_{\rm ch,r}$ $\displaystyle=\frac{ck_{\rm B}T_{\rm
i}}{l_{*}}\left(\frac{12\pi f_{\rm ion}s_{*}R_{\rm cl}^{3}}{\alpha_{\rm
B}M_{\rm cl}\epsilon_{*}}\right)^{1/2}$ $\displaystyle\simeq 26\,{\rm
pc}~{}\left(\frac{\epsilon_{*}}{0.1}\right)^{-1/2}\left(\frac{\Sigma_{\rm
cl}}{80\,M_{\odot}\,{\rm pc}^{-2}}\right)^{-3/4}\left(\frac{M_{\rm
cl}}{10^{5}\,M_{\odot}}\right)^{1/4}\left(\frac{T_{\rm i}}{8000\,{\rm
K}}\right)^{7/5},$ (24)

where we use $f_{\rm ion}=0.73$ in the second equation. We also define another
characteristic radius $r_{\rm ch,g}$ for the equilibrium between the thermal
pressure and the gravity force as

$\displaystyle r_{\rm ch,g}$ $\displaystyle=\left[\frac{2l_{*}R_{\rm
cl}^{3}}{GcM_{\rm
cl}}\left(\frac{\epsilon_{*}}{1-\epsilon_{*}^{2}}\right)r_{\rm
ch,r}\right]^{1/2}$ $\displaystyle\simeq 51\,{\rm
pc}\left(\frac{\epsilon_{*}}{0.1}\right)^{1/4}\left(\frac{\Sigma_{\rm
cl}}{80\,M_{\odot}{\rm pc}^{-2}}\right)^{-9/8}\left(\frac{M_{\rm
cl}}{10^{5}M_{\odot}}\right)^{3/8}\left(\frac{T_{\rm i}}{8000\,{\rm
K}}\right)^{7/10}.$ (25)

Next, we rewrite the equation (20) by introducing dimensionless parameters:

$\displaystyle x_{\rm sh}=r_{\rm sh}/r_{\rm ch,r},$ (26)

and

$\displaystyle\tau=t/t_{\rm ch},$ (27)

where

$\displaystyle t_{\rm ch}$
$\displaystyle=\left[\frac{c}{l_{*}}\left(\frac{1-\epsilon_{*}}{\epsilon_{*}}\right)r_{\rm
ch,r}\right]^{1/2}$ $\displaystyle\simeq 2.9\,{\rm
Myr}\,\left(\frac{\epsilon_{*}}{0.1}\right)^{-3/4}$ $\displaystyle\hskip
28.45274pt\left(\frac{\Sigma_{\rm cl}}{80\,M_{\odot}{\rm
pc}^{-2}}\right)^{-3/4}\left(\frac{M_{\rm
cl}}{10^{5}M_{\odot}}\right)^{1/2}\left(\frac{T_{\rm i}}{8000{\rm
K}}\right)^{7/10}.$ (28)

Substituting equation (26) and (27) into equation (20), we rewrite the
equation of motion as

$\displaystyle\frac{d}{d\tau}\left(x_{\rm sh}^{3}\dot{x}_{\rm
sh}\right)=x_{\rm sh}^{2}(1+x_{\rm sh}-R_{\rm g}^{2}x_{\rm sh}^{2}),$ (29)

where $R_{\rm g}=r_{\rm ch,r}/r_{\rm ch,g}$. On the right-hand side, each term
represents the contributions from the thermal pressure, the radiation force,
and the gravity force. At $x_{\rm sh}\ll 1$, the thermal pressure plays a
dominant role in the acceleration of the shell. The radiation force becomes
important only if the cloud radius is larger than $r_{\rm ch,r}$.
Additionally, the thermal pressure cannot overcome the gravity in the regions
outside $r_{\rm g}$. Here, we assume that the duration time of the star
formation is equal to the crossing time of the shell over the cloud $(t_{\rm
exp})$ as

$\displaystyle M_{*}=\epsilon_{*}M_{\rm cl}=\dot{M}_{*}t_{\rm exp}.$ (30)

where we use the SFR defined in Equations (14). The expanding time $t_{\rm
exp}$ also depends on the SFE $\epsilon_{*}$. We need to solve Equation (29)
and (30) consistently. Furthermore, the parameter of SFR $\epsilon_{\rm ff}$
is given by Equation (16). Note that, however the dependence of the threshold
surface density $\Sigma_{\rm th}$ on the cloud mass and metallicity is still
uncertain.

If the cloud radius $(R_{\rm cl})$ is smaller than $r_{\rm ch,r}$ and $r_{\rm
ch,g}$, the thermal pressure alone contributes the shell dynamics. Therefore,
as a first step, we consider the case only with the thermal pressure. In this
case, we can estimate the SFEs analytically. Considering only the first term
in the right-hand side of Equation (29), the self-similar solution is derived
as $x_{\rm sh}=\tau/\sqrt{3}$ and $\dot{x}_{\rm sh}=1/\sqrt{3}$. The shell
arrives at the surface of the cloud with the time-scale:

$\displaystyle t_{\rm exp,th}$ $\displaystyle=\frac{\sqrt{3}t_{\rm ch}R_{\rm
cl}}{r_{\rm ch,r}}$ $\displaystyle=3.9\,{\rm
Myr}\,\left(\frac{\epsilon_{*}}{0.1}\right)^{-1/4}\left(\frac{\Sigma_{\rm
cl}}{80\,M_{\odot}{\rm pc}^{-2}}\right)^{1/4}\left(\frac{T_{\rm i}}{8000{\rm
K}}\right)^{-7/10}.$ (31)

Substituting $t_{\rm exp,th}$ into (30), we obtain the SFE as

$\displaystyle\epsilon_{*,{\rm th}}\simeq 0.09\left(\frac{\epsilon_{\rm
ff}}{0.1}\right)^{4/5}\left(\frac{\Sigma_{\rm cl}}{80\,M_{\odot}\,{\rm
pc}^{-2}}\right)^{1/2}\left(\frac{M_{\rm
cl}}{10^{5}\,M_{\odot}}\right)^{1/10}\left(\frac{T_{\rm i}}{8000\,{\rm
K}}\right)^{-14/25}.$ (32)

This is same as the model derived in Paper I where we assumed the central
point sources. The SFEs mainly depends on the surface densities, and it is
slightly changed with the cloud mass. The shell velocity does not depend on
the position of the shell in this solution. We obtain the shell velocity as

$\displaystyle v_{\rm sh,th}$ $\displaystyle=\frac{r_{\rm
ch,r}}{\sqrt{3}t_{\rm ch}}$ $\displaystyle=4.9~{}{\rm
km/s}\left(\frac{\epsilon_{\rm ff}}{0.1}\right)^{1/5}$ $\displaystyle\hskip
14.22636pt\left(\frac{\Sigma_{\rm cl}}{80\,M_{\odot}\,{\rm
pc}^{-2}}\right)^{-1/4}\left(\frac{M_{\rm
cl}}{10^{5}\,M_{\odot}}\right)^{3/20}\left(\frac{T_{\rm i}}{8000\,{\rm
K}}\right)^{14/25}.$ (33)

As discussed in Section 3.3, the parameter $\epsilon_{\rm ff}$ can be fit with
the step function transiting at the threshold density $\Sigma_{\rm th}$. If
the velocity of the expanding shell is smaller than the escape velocity from
the core, the shell is likely to fall back and enhance the star formation at
the core. As shown in Figure 10, the thermal pressure mainly contributes to
accelerating the shell around the core. Assuming that the escaping shell from
the core has the same velocity as the self-similar solution given as Equation
(33), the condition of the SFR enhancement is given as

$\displaystyle v_{\rm th}<v_{\rm esc,core}=\sqrt{\frac{2GM_{\rm core}}{R_{\rm
core}}},$ (34)

where $M_{\rm core}$ and $R_{\rm core}$ are the core mass and radius. As
discussed in Section 3.3, the core formation starts when the total stellar
mass exceeds 0.1 times the cloud mass. As shown in Figures 2 and 15, the bound
fraction $(f_{\rm bd})$ and the half-mass radius $(R_{\rm h})$ are typically
$f_{\rm bd}\sim 0.1$ and $R_{\rm h}\sim 0.1$ at this epoch, regardless of the
cloud model. We adopt this values in Equation (34). Substituting $M_{\rm
core}=10^{-2}~{}M_{\rm cl}$ and $R_{\rm core}=0.1~{}R_{\rm cl}$, we obtain the
condition of the surface density for the SFR enhancement as

$\displaystyle\Sigma_{\rm cl}>\Sigma_{\rm th},$ (35)

where

$\displaystyle\Sigma_{\rm th}\simeq 280~{}M_{\odot}{\rm
pc^{-2}}\,\left(\frac{\epsilon_{\rm ff}}{0.1}\right)^{2/5}\left(\frac{M_{\rm
cl}}{10^{6}M_{\odot}}\right)^{-1/5}\left(\frac{T_{\rm i}}{8000~{}{\rm
K}}\right)^{28/25}.$ (36)

In the cases with $(M_{\rm
cl},Z,\alpha_{0})=(10^{6}~{}M_{\odot},Z_{\odot},1)$, the threshold density is
estimated as $\Sigma_{\rm th}=270~{}M_{\odot}{\rm pc^{-2}}$. The SFR
enhancement occurs around $\sim 300~{}M_{\odot}{\rm pc^{-2}}$ in the numerical
simulations. Thus, the estimate in Equation (36) reproduces the simulation
results well.

Table 3: the parameters used in the analytical estimate $M_{\rm cl}\,[M_{\odot}]$ | $Z\,[Z_{\odot}]$ | $\epsilon_{\rm ff,1}$ | $\epsilon_{\rm ff,2}$ | $T_{\rm i}\,[{\rm K}]$ | $\Sigma_{\rm th}\,[M_{\odot}{\rm pc^{-2}}]$
---|---|---|---|---|---
$10^{5}$ | $1$ | $0.08$ | $0.30$ | $8.0\times 10^{3}$ | $4.1\times 10^{2}$
$10^{5}$ | $10^{-1}$ | $0.09$ | $0.22$ | $1.6\times 10^{4}$ | $9.4\times 10^{2}$
$10^{5}$ | $10^{-2}$ | $0.03$ | $0.21$ | $2.2\times 10^{4}$ | $8.8\times 10^{2}$
$10^{6}$ | $1$ | $0.09$ | $0.38$ | $8.0\times 10^{3}$ | $2.7\times 10^{2}$
$10^{6}$ | $10^{-1}$ | $0.06$ | $0.46$ | $1.6\times 10^{4}$ | $4.9\times 10^{2}$
$10^{6}$ | $10^{-2}$ | $0.03$ | $0.40$ | $2.2\times 10^{4}$ | $5.6\times 10^{2}$

We estimate the shell crossing time and the velocity by integrating equation
(29) from $x_{\rm sh}=0$ to the cloud radius $x_{\rm cl}=R_{\rm cl}/r_{\rm
ch,r}$. Equation (29) has the asymptotic solution at $x_{\rm sh}\ll 0$ as
$x_{\rm sh}=1/\sqrt{3}\tau$ and $\dot{x}_{\rm sh}=1/\sqrt{3}$. We use these
solutions as the inside boundary conditions. We iterate the calculations of
Equation (29) until the SFE $(\epsilon_{*})$ satisfies Equation (30). The SFR
is given by Equations (14) and (15). We adopt Equation (16) to obtain the
parameter $\epsilon_{\rm ff}$. The SFR sensitively depends on various
processes such as thermal pressure, radiative feedback and turbulent motions.
Therefore, it is difficult to predict the SFR analytically. In this work, we
use the SFRs obtained from the numerical simulations. The parameters
$\epsilon_{\rm ff,1}$ and $\epsilon_{\rm ff,2}$ are evaluated from the results
of the most diffuse and most compact clouds for each cloud mass and
metallicity. The threshold surface density of the SFR enhancement is given by
Equation (36). The parameters are summarized in Table 3.

Figure 17:  The analytical estimate of star formation efficiencies for the
clouds with $M_{\rm cl}=10^{6}~{}{M_{\odot}}$ (solid) and
$10^{5}~{}{M_{\odot}}$ (dashed). Each symbol represents the simulation results
with $M_{\rm cl}=10^{6}~{}M_{\odot}$ (circle) and $10^{5}~{}M_{\odot}$
(triangle). We only plot the cases with $\alpha_{0}=1$. The colors of the
symbols are the same as the lines.

Figure 17 shows the SFEs obtained by the semi-analytical model. The semi-
analytical model reproduces the simulation results well. In the high- and low-
surface limits, the SFE is the simple power-law function of the surface
density. The SFE jump occurs at the threshold surface density $(\Sigma_{\rm
th})$. The threshold surface densities increases in the lower mass cloud,
e.g., from $\Sigma_{\rm th}\simeq 270~{}M_{\odot}{\rm pc^{-2}}$ at $M_{\rm
cl}=10^{6}~{}M_{\odot}$ to $\Sigma_{\rm th}\simeq 410~{}M_{\odot}{\rm
pc^{-2}}$ at $M_{\rm cl}=10^{5}~{}M_{\odot}$ in the cases with $Z=Z_{\odot}$.
Also, it is sensitive to the temperature of Hii regions, i.e., metallicity.
The thresholds of the clouds with $Z=10^{-2}~{}\rm Z_{\odot}$ are higher than
that for $Z=~{}\rm Z_{\odot}$ by a factor of $\sim 2$.

###  4.2 Mass of bounded objects

Figure 18:  The semi-analytical estimates of the gravitationally bounded
stellar masses in the clouds with $M_{\rm cl}=10^{5}\,M_{\odot}$ (dashed) and
$10^{6}\,M_{\odot}$ (solid) as the function of the surface densities. The
symbols represent the simulation results, and the styles of symbols are same
as Figure 17. The black dashed line represents $M_{\rm
cl}=10^{4}~{}M_{\odot}$.

The gravitationally bound fraction is tightly related to the SFEs, as shown in
Figure 18. Combining the semi-analytical model of the SFEs and the relation
between the SFEs and the bound fractions given as Equation (18), we obtain the
masses of the gravitationally bounded object as shown in Figure 18. We find
that young massive clusters (YMC, $M_{\rm cl}>10^{4}\,M_{\odot}$) with $Z\sim
Z_{\odot}$ are likely to form in the clouds of $10^{6}~{}M_{\odot}$ with the
surface densities higher than $\sim 100\,M_{\odot}{\rm pc}^{-2}$. The mass and
surface density are similar to the typical GMCs in the Milky Way (Roman-Duval
et al., 2010). On the other hand, at $Z<Z_{\odot}$, YMCs forms from clouds
with $10^{6}~{}M_{\odot}$ only if the surface density is larger than
$\Sigma_{\rm cl}\sim 200~{}M_{\odot}{\rm pc^{-2}}$ ($\sim 350~{}M_{\odot}{\rm
pc^{-2}}$) at $Z=10^{-1}~{}Z_{\odot}$ ($Z=10^{-2}~{}Z_{\odot}$). Besides, the
condition of YMC formation for $10^{5}~{}M_{\odot}$ is strict. Thus, we
suggest that YMCs do not form frequently in low-metallicity environments.

##  5 Summary and Discussion

We have studied the formation processes of young massive star clusters (YMCs,
$>10^{4}~{}M_{\odot}$) by performing three-dimensional radiation hydrodynamics
simulations. Our simulations include photoionization, photodissociation of
molecules, radiation pressure. We have investigated giant-molecular clouds
(GMCs) with the wide parameter ranges: $\Sigma\simeq 80-3200~{}M_{\odot}{\rm
pc^{-2}}$, $Z=10^{-2}Z_{\odot}-Z_{\odot}$ and the virial parameter
$\alpha_{0}=1$ or $2$. Also, we have developed a semi-analytical model that
reproduces the simulation results nicely. Our findings are summarized as
follows:

(i)

In the early phase of the cloud collapse, turbulent motions induce the
formation of high-density filaments in which stars form. Then, after one free-
fall time, the star formation rate rapidly increases at the center of a cloud.

(ii)

The star clusters are gravitationally bounded once the total stellar mass
exceeds 0.1 times the initial cloud mass. Simultaneously, the stellar core
forms if the cloud is compact enough. In such a case, the gravitational force
from the star cluster can pull even ionized hot gas and induces long-lasting
star formation. We obtain the threshold surface density for the core formation
as $\Sigma_{\rm th}\sim 300~{}M_{\odot}{\rm pc^{-2}}$ in the cases with the
cloud mass $M_{\rm cl}=10^{6}~{}M_{\odot}$ at $Z=Z_{\odot}$.

(iii)

The star formation efficiencies (SFEs) sensitively depend on the surface
densities of clouds. If the surface density is lower than the threshold value
of the core formation, the photoionization feedback easily disrupts clouds,
resulting in the SFEs lower than 0.2. In compact clouds, the star formation
can continue against the photoionization feedback. In these cases, the SFEs
are $\epsilon_{*}\geq 0.3$, regardless of metallicities and virial parameters.

(iv)

The SFE sensitively depends on the metallicity in the cases of diffuse clouds.
The temperature of Hii regions of low-metallicity gas is higher than that of
metal-rich gas because of the lower cooling rate. Therefore, the Hii regions
rapidly expand and suppress the star formation. On the other hand, in the
cases of massive compact clouds, the SFE does not change significantly with
the metallicity, because the photoionization feedback is not dominant.

(v)

Turbulent motions delay the gravitational contraction of the clouds.
Therefore, in the cases with higher turbulence energy, the star formation
proceeds slowly and the SFEs become lower. Also, we find that the star
clusters can be unbounded due to the strong turbulent motion.

(vi)

We have constructed the semi-analytical model reproducing the simulation
results. We have used the expanding shell model, considering the thermal
pressure of Hii regions, radiation pressure, and gravity from the formed star
cluster and the shell. The SFEs are estimated based on the assumption in which
the star formation continues only for the shell expanding. We also consider
the enhancement of the star formation rates due to the stellar core formation.
The threshold surface density of the core formation $(\Sigma_{\rm th})$ is
analytically estimated by considering the condition that the velocity of the
expanding shell is comparable to the escape velocity of the stellar core. This
threshold density increases with the lower-metallicity or the smaller cloud,
e.g., $\Sigma_{\rm th}\sim 560~{}M_{\odot}{\rm pc^{-2}}$ ($410~{}M_{\odot}{\rm
pc^{-2}}$) with the cloud mass $M_{\rm cl}=10^{6}~{}M_{\odot}$ and the
metallicity $Z=10^{-2}Z_{\odot}$ ($10^{5}~{}M_{\odot}$ and $Z_{\odot}$).

(vii)

Combining the analytical estimation of SFEs and the bound fraction obtained
from the numerical simulations, we have derived the condition for the YMC
formation. We have shown that massive and higher-metallicity clouds are
favored for the formation of YMCs. In the clouds with $M_{\rm
cl}=10^{6}~{}M_{\odot}$, the YMCs are formed if the surface density is larger
than $100~{}M_{\odot}{\rm pc^{-2}}$.

In the lower-metallicity environments, higher surface densities are needed for
producing YMCs. In particular, we predicted that YMCs are formed only in
clouds with $\Sigma_{\rm cl}\gtrsim 350~{}M_{\odot}{\rm pc^{-2}}$ at
$Z=10^{-2}~{}Z_{\odot}$ even if the cloud mass is $10^{6}~{}M_{\odot}$. These
conditions are related to the minimum metallicity of observed globular
clusters (GCs: Forbes et al., 2018). Forbes et al. (2018) showed the
metallicity floor of GCs in the Milky Way with $Z\sim 10^{-2.5}Z_{\odot}$ (but
see , Wan et al., 2020). Recently Larsen et al. (2020) observed the lower
metallicity GC of $Z\sim 10^{-2.9}Z_{\odot}$ in Andromeda galaxy. Kruijssen
(2019) pointed out that this floor is related to the minimum mass for GCs
formation of the host galaxies and their maximum redshift. We have suggested
additional condition that the higher surface densities $(\Sigma_{\rm
cl}\gtrsim 350~{}M_{\odot}{\rm pc}^{-2})$ is required for the formation of
low-metallicity GCs. Previous studies indicated that the low-metallicity
massive clouds are formed in the early galaxies (e.g., Ricotti et al., 2016;
Arata et al., 2018). We will investigate the star cluster formation in
galaxies modeled by cosmological simulations in future work.

Pfalzner (2009) suggested that there are two types of young star clusters with
$\sim 10^{4}~{}M_{\odot}$ , YMCs, and "leaky clusters." Leaky clusters show
the low-stellar density ($\sim 1-10~{}M_{\odot}{\rm pc^{-3}}$) compared with
YMCs ($\gtrsim 10^{3}~{}M_{\odot}{\rm pc^{-3}}$, Portegies Zwart et al.,
2010). Fujii & Portegies Zwart (2016) showed that these two populations had
different surface densities at birth (also see, Pfalzner, 2011). In our
simulation with $10^{6}~{}M_{\odot}$ clouds at $Z=Z_{\odot}$, the SFEs jump
occurs in the narrow surface density range. This SFE jump is likely to
separate star clusters into the two populations. Note that, in the low-
metallicity cases, it is difficult for the lower surface density clouds to
form star clusters with $10^{4}~{}{M_{\odot}}$ due to the low SFE. Therefore,
we indicate that these two populations do not emerge in the low-metallicity
environments.

In this work, we have focused on radiative feedback. There are other feedback
processes such as stellar wind, outflow, and supernovae. These feedback
processes are also likely to regulate star formation. For example, the stellar
wind pushes out the ambient gas. This makes low-density bubbles where Hii
regions can expand rapidly. Recently, Geen et al. (2020) showed the
combination effects of the stellar wind and the photoionization processes (see
also, Dale et al., 2014; Geen et al., 2020; Decataldo et al., 2020). The
collimated outflow from stars also has impacts on the star formation via the
injection of turbulence energy (Li & Nakamura, 2006; Nakamura & Li, 2007).
Recently, Guszejnov et al. (2021) found that the outflow has a crucial role in
regulating the IMF. They also found that the low-surface density cloud $(\sim
60~{}M_{\odot}{\rm pc^{-2}})$ with $2\times 10^{4}~{}M_{\odot}$ wad disrupted
only by the outflow. The properties of the outflow sensitively depend on the
physical states of the star-forming regions. The direction of the outflow from
the single star depends on both angular momentum and magnetic fields (e.g.,
Hirano et al., 2020), and the outflow is failed in massive and weakly
magnetized massive core (Machida & Hosokawa, 2020). Some cases in our
simulations have cloud lifetimes longer than the lifetimes of OB stars. In
this case, SNe should turn on and evacuate the gas (e.g., Geen et al., 2016).
The magnetic field can affect the gas structures and suppress the star
formation rate by a factor of $\sim 2$ (e.g., Federrath & Klessen, 2012). It
also suppresses shell expansion caused by radiative feedback or stellar wind
(e.g., Offner & Liu, 2018; Kim et al., 2021). In future works, we will
incorporate these feedback effects in our simulations and investigate the
impacts of multiple feedback.

In this work, we have assumed the constant mass-to-light ratio for sink
particles. This assumption is reasonable only for massive star clusters with a
mass larger than $10^{4}~{}M_{\odot}$ (Kim et al., 2016). In practice, in the
cases of clouds with $\lesssim 10^{5}~{}M_{\odot}$, the masses of star
clusters are below $10^{4}~{}M_{\odot}$, which indicates that the expected
number of massive stars can be lower than unity. In such a case, formation
sites and timings of radiation sources depend on the models. Grudić & Hopkins
(2019) showed that the SFEs could change with models of star formation by a
factor of $\sim 3$ for the clouds with $M_{\rm cl}=10^{5}~{}M_{\odot}$ and
$R_{\rm cl}=20~{}{\rm pc}$. Also, total emissivities of star clusters depend
on the resolution in the cases with diffuse clouds of $M_{\rm
cl}=10^{5}~{}M_{\odot}$ (see Appendix E). Thus, we suggest that more
sophisticated modeling of star formation should be developed in studying star
clusters in clouds with $\lesssim 10^{5}~{}M_{\odot}$ in the future.

We have followed the formation and dissociation of CO molecules. We have found
that the spatial distribution and total abundance of CO molecules depended on
the metallicity. Also, the CO distribution changes with the cloud evolution.
Therefore, the surface brightness maps of CO lines can be a tool to understand
the physical properties of star-forming clouds. The CO observations in the
low-metallicity galaxies (such as LMC/SMC) are progressing (e.g., Fukui et
al., 2015; Muraoka et al., 2017; Tokuda et al., 2019). We plan to model the CO
map to compare with the observations in a future study.

## Acknowledgements

The authors wish to express their cordial thanks to Profs. Masayuki Umemura
and Ken Ohsuga for their continual interest, advice, and encouragement. We
would like to thank Michiko Fujii, Tomoaki Matsumoto, Takashi Hosokawa,
Kazuyuki Sugimura, Kazutaka Kimura, Kazuyuki Omukai and Shu-ichiro Inutsuka
for useful discussions and comments. The numerical simulations were performed
on the Cray XC50 (Aterui II) at the Center for Computational Astrophysics of
National Astronomical Observatory of Japan and Yukawa-21 at Yukawa Institute
for Theoretical Physics in Kyoto University. This work is supported in part by
the MEXT/JSPS KAKENHI Grant Number 17H04827, 20H04724, National Astronomical
Observatory of Japan (NAOJ) ALMA Scientific Research Grant Number 2019-11A,
and JST Fusion Oriented REsearch for disruptive Science and Technology
(FOREST) (HY).

## Data Availability

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

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abel & Wandelt (2002) Abel T., Wandelt B. D., 2002, MNRAS, 330, L53
  * Abel et al. (1997) Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astron., 2, 181
  * Ali (2021) Ali A. A., 2021, MNRAS, 501, 4136
  * Arata et al. (2018) Arata S., Yajima H., Nagamine K., 2018, MNRAS, 475, 4252
  * Baumgardt & Kroupa (2007) Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589
  * Bertoldi & McKee (1992) Bertoldi F., McKee C. F., 1992, ApJ, 395, 140
  * Bressert et al. (2012) Bressert E., Ginsburg A., Bally J., Battersby C., Longmore S., Testi L., 2012, ApJ, 758, L28
  * Cen (1992) Cen R., 1992, ApJS, 78, 341
  * 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
  * Chevance et al. (2020a) Chevance M., et al., 2020a, arXiv e-prints, p. arXiv:2010.13788
  * Chevance et al. (2020b) Chevance M., et al., 2020b, MNRAS, 493, 2872
  * Commerçon et al. (2011) Commerçon B., Teyssier R., Audit E., Hennebelle P., Chabrier G., 2011, A&A, 529, A35
  * Dale et al. (2012) Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 427, 2852
  * Dale et al. (2013) Dale J. E., Ercolano B., Bonnell I. A., 2013, MNRAS, 430, 234
  * Dale et al. (2014) Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694
  * Decataldo et al. (2020) Decataldo D., Lupi A., Ferrara A., Pallottini A., Fumagalli M., 2020, MNRAS, 497, 4718
  * Dobbs et al. (2020) Dobbs C. L., Liow K. Y., Rieder S., 2020, MNRAS, 496, L1
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Draine & Bertoldi (1996) Draine B. T., Bertoldi F., 1996, ApJ, 468, 269
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath et al. (2010) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, ApJ, 713, 269
  * Ferland et al. (1992) Ferland G. J., Peterson B. M., Horne K., Welsh W. F., Nahar S. N., 1992, ApJ, 387, 95
  * Forbes et al. (2018) Forbes D. A., et al., 2018, Proceedings of the Royal Society of London Series A, 474, 20170616
  * Forrey (2013) Forrey R. C., 2013, ApJ, 773, L25
  * Fujii & Portegies Zwart (2014) Fujii M. S., Portegies Zwart S., 2014, MNRAS, 439, 1003
  * Fujii & Portegies Zwart (2016) Fujii M. S., Portegies Zwart S., 2016, ApJ, 817, 4
  * Fujii et al. (2017) Fujii M. S., Tanikawa A., Makino J., 2017, PASJ, 69, 94
  * Fujii et al. (2021) Fujii M. S., Saitoh T. R., Hirai Y., Wang L., 2021, arXiv e-prints, p. arXiv:2103.02829
  * Fukui & Kawamura (2010) Fukui Y., Kawamura A., 2010, ARA&A, 48, 547
  * Fukui et al. (2015) Fukui Y., et al., 2015, ApJ, 807, L4
  * Fukushima et al. (2018) Fukushima H., Omukai K., Hosokawa T., 2018, MNRAS, 473, 4754
  * Fukushima et al. (2020a) Fukushima H., Hosokawa T., Chiaki G., Omukai K., Yoshida N., Kuiper R., 2020a, MNRAS, 497, 829
  * Fukushima et al. (2020b) Fukushima H., Yajima H., Sugimura K., Hosokawa T., Omukai K., Matsumoto T., 2020b, MNRAS, 497, 3830
  * Galli & Palla (1998) Galli D., Palla F., 1998, A&A, 335, 403
  * Gavagnin et al. (2017) Gavagnin E., Bleuler A., Rosdahl J., Teyssier R., 2017, MNRAS, 472, 4155
  * Geen et al. (2016) Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2016, MNRAS, 463, 3129
  * Geen et al. (2017) Geen S., Soler J. D., Hennebelle P., 2017, MNRAS, 471, 4844
  * Geen et al. (2020) Geen S., Bieri R., Rosdahl J., de Koter A., 2020, MNRAS,
  * Glover (2015) Glover S. C. O., 2015, MNRAS, 453, 2901
  * Glover & Jappsen (2007) Glover S. C. O., Jappsen A. K., 2007, ApJ, 666, 1
  * Gnedin & Abel (2001) Gnedin N. Y., Abel T., 2001, New Astron., 6, 437
  * Gong & Ostriker (2013) Gong H., Ostriker E. C., 2013, ApJS, 204, 8
  * Goodwin (2009) Goodwin S. P., 2009, Ap&SS, 324, 259
  * Grudić & Hopkins (2019) Grudić M. Y., Hopkins P. F., 2019, MNRAS, 488, 2970
  * Grudić et al. (2018) Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, MNRAS, 475, 3511
  * Grudić et al. (2020) Grudić M. Y., Kruijssen J. M. D., Faucher-Giguère C.-A., Hopkins P. F., Ma X., Quataert E., Boylan-Kolchin M., 2020, arXiv e-prints, p. arXiv:2008.04453
  * Guszejnov et al. (2021) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 502, 3646
  * Harten et al. (1983) Harten A., Lax P. D., Leer B. v., 1983, SIAM Review, 25, 35
  * Hasegawa & Semelin (2013) Hasegawa K., Semelin B., 2013, MNRAS, 428, 154
  * He et al. (2019) He C.-C., Ricotti M., Geen S., 2019, MNRAS, 489, 1880
  * Hirano et al. (2020) Hirano S., Tsukamoto Y., Basu S., Machida M. N., 2020, ApJ, 898, 118
  * Hollenbach & McKee (1979) Hollenbach D., McKee C. F., 1979, ApJS, 41, 555
  * Hosokawa & Inutsuka (2005) Hosokawa T., Inutsuka S.-i., 2005, ApJ, 623, 917
  * Hosokawa & Inutsuka (2006) Hosokawa T., Inutsuka S.-i., 2006, ApJ, 646, 240
  * Hosokawa et al. (2016) Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N., 2016, ApJ, 824, 119
  * Howard et al. (2017) Howard C. S., Pudritz R. E., Harris W. E., 2017, MNRAS, 470, 3346
  * Inoguchi et al. (2020) Inoguchi M., Hosokawa T., Mineshige S., Kim J.-G., 2020, MNRAS, 497, 5061
  * Kannan et al. (2019) Kannan R., Vogelsberger M., Marinacci F., McKinnon R., Pakmor R., Springel V., 2019, MNRAS, 485, 117
  * Kauffmann & Pillai (2010) Kauffmann J., Pillai T., 2010, ApJ, 723, L7
  * Kawamura et al. (2009) Kawamura A., et al., 2009, ApJS, 184, 1
  * 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
  * Kim et al. (2020) Kim J., et al., 2020, arXiv e-prints, p. arXiv:2012.00019
  * Kim et al. (2021) Kim J.-G., Ostriker E. C., Filippova N., 2021, ApJ, 911, 128
  * Klassen et al. (2014) Klassen M., Kuiper R., Pudritz R. E., Peters T., Banerjee R., Buntemeyer L., 2014, ApJ, 797, 4
  * Kreckel et al. (2010) Kreckel H., Bruhns H., Čížek M., Glover S. C. O., Miller K. A., Urbain X., Savin D. W., 2010, Science, 329, 69
  * Kroupa et al. (2001) Kroupa P., Aarseth S., Hurley J., 2001, MNRAS, 321, 699
  * Kruijssen (2014) Kruijssen J. M. D., 2014, Classical and Quantum Gravity, 31, 244006
  * Kruijssen (2019) Kruijssen J. M. D., 2019, MNRAS, 486, L20
  * Kruijssen et al. (2019) Kruijssen J. M. D., et al., 2019, Nature, 569, 519
  * Krumholz (2012) Krumholz M. R., 2012, ApJ, 759, 9
  * Krumholz (2018) Krumholz M. R., 2018, MNRAS, 480, 3468
  * Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352
  * Krumholz et al. (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
  * Krumholz et al. (2009) Krumholz M. R., Klein R. I., McKee C. F., Offner S. S. R., Cunningham A. J., 2009, Science, 323, 754
  * Krumholz et al. (2011) Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 740, 74
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland -Hawthorn J., 2019, ARA&A, 57, 227
  * Kuiper et al. (2010) Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ, 722, 1556
  * Kumamoto et al. (2019) Kumamoto J., Fujii M. S., Tanikawa A., 2019, MNRAS, 486, 3942
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Lada et al. (1984) Lada C. J., Margulis M., Dearborn D., 1984, ApJ, 285, 141
  * Laor & Draine (1993) Laor A., Draine B. T., 1993, ApJ, 402, 441
  * Larsen et al. (2020) Larsen S. S., Romanowsky A. J., Brodie J. P., Wasserman A., 2020, arXiv e-prints, p. arXiv:2010.07395
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lee et al. (1996) Lee H. H., Herbst E., Pineau des Forets G., Roueff E., Le Bourlot J., 1996, A&A, 311, 690
  * Leroy et al. (2015) Leroy A. K., et al., 2015, ApJ, 801, 25
  * Leroy et al. (2016) Leroy A. K., et al., 2016, ApJ, 831, 16
  * 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
  * Li & Nakamura (2006) Li Z.-Y., Nakamura F., 2006, ApJ, 640, L187
  * Li et al. (2019) Li H., Vogelsberger M., Marinacci F., Gnedin O. Y., 2019, MNRAS, 487, 364
  * Liow & Dobbs (2020) Liow K. Y., Dobbs C. L., 2020, MNRAS, 499, 1099
  * 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. p. 291 (arXiv:1401.4175), doi:10.2458/azu_uapress_9780816531240-ch013
  * Machida & Hosokawa (2020) Machida M. N., Hosokawa T., 2020, MNRAS, 499, 4490
  * Maeda et al. (2020) Maeda R., Inoue T., Fukui Y., 2020, arXiv e-prints, p. arXiv:2011.11650
  * Martin et al. (1998) Martin P. G., Keogh W. J., Mandy M. E., 1998, ApJ, 499, 793
  * Matsumoto (2007) Matsumoto T., 2007, PASJ, 59, 905
  * Matsumoto et al. (2015) Matsumoto T., Dobashi K., Shimoikura T., 2015, ApJ, 801, 77
  * Matzner (2002) Matzner C. D., 2002, ApJ, 566, 302
  * McKee & Williams (1997) McKee C. F., Williams J. P., 1997, ApJ, 476, 144
  * Mihalas & Mihalas (1984) Mihalas D., Mihalas B. W., 1984, Foundations of radiation hydrodynamics
  * Millar (1991) Millar T. J., 1991, A&A, 242, 241
  * Muraoka et al. (2017) Muraoka K., et al., 2017, ApJ, 844, 98
  * Murray (2011) Murray N., 2011, ApJ, 729, 133
  * Murray et al. (2010) Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191
  * Nakamura & Li (2007) Nakamura F., Li Z.-Y., 2007, ApJ, 662, 395
  * Nakatani et al. (2018) Nakatani R., Hosokawa T., Yoshida N., Nomura H., Kuiper R., 2018, ApJ, 857, 57
  * Nelson & Langer (1997) Nelson R. P., Langer W. D., 1997, ApJ, 482, 796
  * Offner & Liu (2018) Offner S. S. R., Liu Y., 2018, Nature Astronomy, 2, 896
  * Omukai (2000) Omukai K., 2000, ApJ, 534, 809
  * Omukai & Nishi (1998) Omukai K., Nishi R., 1998, ApJ, 508, 141
  * Omukai et al. (2005) Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627
  * Osterbrock (1989) Osterbrock D. E., 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Palla et al. (1983) Palla F., Salpeter E. E., Stahler S. W., 1983, ApJ, 271, 632
  * Petrosian et al. (1972) Petrosian V., Silk J., Field G. B., 1972, ApJ, 177, L69
  * Pfalzner (2009) Pfalzner S., 2009, A&A, 498, L37
  * Pfalzner (2011) Pfalzner S., 2011, A&A, 536, A90
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Raskutti et al. (2016) Raskutti S., Ostriker E. C., Skinner M. A., 2016, ApJ, 829, 130
  * Ricotti et al. (2016) Ricotti M., Parry O. H., Gnedin N. Y., 2016, ApJ, 831, 204
  * Roman-Duval et al. (2010) Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492
  * Rosdahl & Teyssier (2015) Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380
  * Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188
  * Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  * Shapiro & Kang (1987) Shapiro P. R., Kang H., 1987, ApJ, 318, 32
  * Shukirgaliyev et al. (2017) Shukirgaliyev B., Parmentier G., Berczik P., Just A., 2017, A&A, 605, A119
  * Skinner & Ostriker (2013) Skinner M. A., Ostriker E. C., 2013, ApJS, 206, 21
  * Skinner & Ostriker (2015) Skinner M. A., Ostriker E. C., 2015, ApJ, 809, 187
  * Smith (2021) Smith M. C., 2021, MNRAS, 502, 5417
  * Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters
  * Sugimura et al. (2017) Sugimura K., Hosokawa T., Yajima H., Omukai K., 2017, MNRAS, 469, 62
  * Sugimura et al. (2020) Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ, 892, L14
  * Sun et al. (2018) Sun J., et al., 2018, ApJ, 860, 172
  * Sun et al. (2020) Sun J., et al., 2020, ApJ, 901, L8
  * Takahashi & Ohsuga (2013) Takahashi H. R., Ohsuga K., 2013, ApJ, 772, 127
  * Tielens & Hollenbach (1985) Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 722
  * Tokuda et al. (2019) Tokuda K., et al., 2019, ApJ, 886, 15
  * Trebitsch et al. (2017) Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224
  * Tsuge et al. (2020) Tsuge K., Fukui Y., Tachihara K., Sano H., Tokuda K., Ueda J., Iono D., Finn M. K., 2020, PASJ,
  * Urquhart et al. (2018) Urquhart J. S., et al., 2018, MNRAS, 473, 1059
  * Vázquez-Semadeni et al. (2010) Vázquez-Semadeni E., Colín P., Gómez G. C., Ballesteros-Paredes J., Watson A. W., 2010, ApJ, 715, 1302
  * Vutisalchavakul et al. (2016) Vutisalchavakul N., Evans Neal J. I., Heyer M., 2016, ApJ, 831, 73
  * Wan et al. (2020) Wan Z., et al., 2020, Nature, 583, 768
  * Williams & McKee (1997) Williams J. P., McKee C. F., 1997, ApJ, 476, 166
  * Wise et al. (2012) Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50
  * Wolcott-Green & Haiman (2019) Wolcott-Green J., Haiman Z., 2019, MNRAS, 484, 2467
  * Wolfire & Cassinelli (1987) Wolfire M. G., Cassinelli J. P., 1987, ApJ, 319, 850
  * Wolfire et al. (2010) Wolfire M. G., Hollenbach D., McKee C. F., 2010, ApJ, 716, 1191
  * Yajima et al. (2017) Yajima H., Nagamine K., Zhu Q., Khochfar S., Dalla Vecchia C., 2017, ApJ, 846, 30
  * Yajima et al. (2020) Yajima H., et al., 2020, arXiv e-prints, p. arXiv:2011.11663
  * van Dishoeck & Black (1988) van Dishoeck E. F., Black J. H., 1988, ApJ, 334, 771

##  Appendix A Chemical network

In table 4, we summarize the chemical reaction related with ${\rm H}$, ${\rm
H_{2}}$, ${\rm H}^{+}$, ${\rm H}^{-}$, ${\rm H_{2}}^{+}$, and ${\rm e}$. We
also include $\rm CO$, $\rm C^{+}$, $\rm O$, ${\rm O^{+}}$, and ${\rm O}^{2+}$
as the metal species. We use the chemical network of Nelson & Langer (1997)
for CO formation as in Fukushima et al. (2020b).

Table 4: Chemical Reactions

Number | Reaction | Rate Coefficient | Reference
---|---|---|---
${\rm H0}$ | ${\rm H}+{\rm e}\rightarrow{\rm H}^{+}+2e$ | $k_{\rm H0}=\exp[-32.71396786+13.536556\times\left(\ln T\left({\rm ev}\right)\right)-5.73932875\left(\ln T\left({\rm ev}\right)\right)^{2}$ | 1
|  | $+1.563154998\left(\ln T\left({\rm ev}\right)\right)^{3}-0.2877056\left(\ln T\left({\rm ev}\right)\right)^{4}+3.48255977\times 10^{-2}\left(\ln T\left({\rm ev}\right)\right)^{5}$ |
|  | $-2.63197617\times 10^{-3}\left(\ln T\left({\rm ev}\right)\right)^{6}+1.11954395\times 10^{-4}\left(\ln T\left({\rm ev}\right)\right)^{7}-2.03914985\times 100^{-6}\left(\ln T\left({\rm ev}\right)\right)^{8}$ |
${\rm H1}$ | ${\rm H^{+}}+e\rightarrow{\rm H}+\gamma$ |  $k_{\rm H1,A}=1.269\times 10^{-13}(315614/T)^{1.503}[1+(604625/T)^{0.47}]^{-1.923}$ (case A) | 2, 3
|  |  $k_{\rm H1,B}=2.753\times 10^{-14}(315614/T)^{1.5}[1+(115188/T)^{0.407}]^{-2.242}$ (case B) |
${\rm H2}$ | ${\rm H^{-}}+{\rm H}\rightarrow{\rm H_{2}}+{\rm e}$ | $k_{\rm H2}=1.35\times 10^{-9}\left[T^{0.098493}+0.32852T^{0.5561}+2.771\times 10^{-7}T^{2.1826}\right]$ | 4
|  | $/\left[1+6.191\times 10^{-3}T^{1.0461}+8.9712\times 10^{-11}T^{3.0424}+3.2576\times 10^{-14}T^{3.7741}\right]$ |
${\rm H3}$ | ${\rm H_{2}}+{\rm H}^{+}\rightarrow{\rm H_{2}}^{+}+{\rm H}$ |  $k_{\rm H3}=3\times 10^{-10}\exp(-21050/T)$ $T<10^{4}~{}{\rm K}$ | 5
|  |  $1.5\times 10^{-10}\exp(-14000/T)$ $T>10^{4}~{}{\rm K}$ |
${\rm H4}$ | ${\rm H_{2}}+{\rm e}\rightarrow 2{\rm H}+{\rm e}$ | $k_{\rm H4}=4.4\times 10^{-10}T^{0.35}\exp(-102000/T)$ | 5
${\rm H5}$ | ${\rm H2}+{\rm H}\rightarrow 3{\rm H}$ | see the reference | 6
${\rm H6}$ | $3{\rm H}\rightarrow{\rm H_{2}}+{\rm H}$ | $k_{\rm H6}=6\times 10^{-32}T^{-0.25}+2\times 10^{-31}T^{-0.5}$ | 7
${\rm H7}$ | $2{\rm H}+{\rm H_{2}}\rightarrow 2{\rm H_{2}}$ | $k_{\rm H7}=k_{\rm H6}/8$ | 8
${\rm H8}$ | $2{\rm H_{2}}\rightarrow 2{\rm H}+{\rm H_{2}}$ | $k_{\rm H8}=k_{\rm high}^{1-a}k_{\rm low}^{a}$ | 8
|  |  $k_{\rm low}=1.18\times 10^{-10}\exp(-6.95\times 10^{4}/T)$ |
|  |  $k_{\rm high}=8.125\times 10^{-8}T^{-1/2}\exp(-5.2\times 10^{4}/T)\left[1-\exp(-6000/T)\right]$ |
|  |  $a=\left(1+n/n_{\rm cr}\right)^{-1}$, $\log_{10}(n_{\rm cr})=4.845-1.3\log_{10}\left(T/10^{4}\right)+1.62\left[\log_{10}\left(T/10^{4}\right)\right]^{2}$ |
${\rm H9}$ | ${\rm H}+{\rm e}\rightarrow{\rm H}^{-}+\gamma$ | $k_{\rm H9}=1.4\times 10^{-18}T^{0.928}\exp(-T/16200)$ | 5
${\rm H10}$ | $2{\rm H}\rightarrow{\rm H}^{+}+{\rm e}+{\rm H}$ | $k_{\rm H10}=1.7\times 10^{-4}k_{\rm H0}$ | 5
${\rm H11}$ | ${\rm H}^{-}+{\rm e}\rightarrow{\rm H}+2{\rm e}$ | $k_{\rm H11}=\exp[-18.01849334+2.3608522\times\left(\ln T\left({\rm ev}\right)\right)-0.28274430\left(\ln T\left({\rm ev}\right)\right)^{2}+1.62331664\times 10^{-2}\left(\ln T\left({\rm ev}\right)\right)^{3}$ | 1
|  |  $-3.36501203\times 10^{-2}\left(\ln T\left({\rm ev}\right)\right)^{4}+1.17832978\times 10^{-2}\left(\ln T\left({\rm ev}\right)\right)^{5}-1.65619470\times 10^{-3}\left(\ln T\left({\rm ev}\right)\right)^{6}$ |
|  |  $+1.06827520\times 10^{-4}\left(\ln T\left({\rm ev}\right)\right)^{7}-2.63128581\times 10^{-6}\left(\ln T\left({\rm ev}\right)\right)^{8}$ |
${\rm H12}$ | ${\rm H}^{-}+{\rm H}^{+}\rightarrow{\rm H_{2}}^{+}+{\rm e}$ |  $k_{\rm H12}=6.9\times 10^{-9}/T^{0.35}$ $T<8000~{}{\rm K}$ | 5
|  |  $9.6\times 10^{-7}/T^{0.9}$ $T>8000~{}{\rm K}$ |
${\rm H13}$ | ${\rm H}^{-}+{\rm H}^{+}\rightarrow 2{\rm H}$ | $k_{\rm H13}=6.3\times 10^{-8}+5.7\times 10^{-6}/\sqrt{T}-9.2\times 10^{-11}\sqrt{T}+4.4\times 10^{-13}T$ | 5
${\rm H14}$ | ${\rm H}+{\rm H}^{+}\rightarrow{\rm H_{2}}^{+}+\gamma$ | $k_{\rm H14}=10^{-19.38-1.523\log T+1.118(\log T)^{2}-0.1269(\log T)^{3}}$ | 5
${\rm H15}$ | ${\rm H_{2}}^{+}+{\rm H}\rightarrow{\rm H_{2}}+{\rm H}^{+}$ | $k_{\rm H15}=6.4\times 10^{-10}$ | 5
${\rm H16}$ | ${\rm H_{2}}^{+}+{\rm e}\rightarrow 2{\rm H}$ | $k_{\rm H16}=2\times 10^{-7}/\sqrt{T}$ | 5
${\rm H17}$ | ${\rm H_{2}}^{+}+{\rm H}^{-}\rightarrow{\rm H_{2}}+{\rm H}$ | $k_{\rm H17}=2.3\times 10^{-7}/\sqrt{T/300}$ | 9
${\rm H18}$ | ${\rm 2H}+{\rm grain}\rightarrow{\rm H_{2}}$ | $k_{\rm H18}=6.0\times 10^{-17}\sqrt{T/300}f_{a}\left(Z/Z_{\odot}\right)[1.0+4.0\times 10^{-2}\sqrt{T+T_{\rm gr}}+2.0\times 10^{-3}T+8.0\times 10^{-6}T^{2}]^{-1}$ | 10
|  |  $f_{a}=[1.0+\exp(7.5\times 10^{2}(1/75-T_{\rm gr}^{-1}))]^{-1}$ |
${\rm M1}$ | ${\rm C^{+}}+{\rm O}\rightarrow{\rm CO}$ | see the references | 11,12
${\rm RH1}$ | ${\rm H}+\gamma\rightarrow{\rm H}^{+}+{\rm e}$ | $R_{\rm HI}$ | Eq. (62)
${\rm RH2}$ | ${\rm H_{2}}+\gamma\rightarrow 2{\rm H}$ | $R_{\rm H_{2}}$ | Eq. (66)
${\rm RH3}$ | ${\rm H}^{-}+\gamma\rightarrow{\rm H}+{\rm e}$ | $R_{\rm H^{-}}$ | 13
${\rm RCO}$ | ${\rm CO}+\gamma\rightarrow C^{+}+{\rm O}$ | $R_{\rm CO}$ | Eq. (66)

References. (1) Abel et al. (1997) (2) Ferland et al. (1992) (3) Glover &
Jappsen (2007) (4) Kreckel et al. (2010) (5) Galli & Palla (1998) (6) Martin
et al. (1998) (7) Forrey (2013) (8) Palla et al. (1983) (9) Millar (1991) (10)
Tielens & Hollenbach (1985) (11) Nelson & Langer (1997) (12) Nakatani et al.
(2018) (13) Hosokawa et al. (2016)

##  Appendix B Thermal processes

In table 5, we summarize the heating and cooling processes included in our
simulations. We incorporate line cooling of $\rm H_{2}$ rovibrational
transitions and metal line cooling of Cii, CO, Oi, Oii, and Oiii. We use the
fitting function derived in Glover (2015) for the line cooling from the $\rm
H_{2}$ rovibrational transition, and the escape probability tabulated by
Fukushima et al. (2018). The cooling rates of Cii, Oi, Oii, and Oiii are
estimated by solving the statistical equilibrium of each energy level as
Fukushima et al. (2020a). Here, we assume that the ionization rate of Oi is
the same as Hi because the ionization energies of Oi and Hi are almost the
same. We solve equilibrium state between Oii and Oiii as Fukushima et al.
(2020a). The energy transport between gas and dust grain is formulated as the
function of dust temperature $T_{\rm d}$ (Omukai et al., 2005).

Table 5: Thermal processes.

Number | Process | Rate ($\rm erg\,cm^{-3}\,s^{-1}$) | Reference
---|---|---|---
Heating
1 |  ${\rm H_{2}}$ formation | $\Gamma_{1}=[\,3.73(1+n_{\rm cr}/n_{\rm H})^{-1}k_{\rm H2}n({\rm H})n({\rm H^{-}})+4.48(1+n_{\rm cr}/n_{\rm H})^{-1}k_{\rm H6}n^{3}({\rm H})$ | 1,2
|  |  $+\,\left(0.2+4.2(1+n_{\rm cr}/n_{\rm H})^{-1}\right)k_{\rm H18}n^{2}({\rm H})]~{}\rm{eV}$ |
|  | $n_{\rm cr}=10^{6}T^{-1/2}/\left\\{1.6n({\rm H})/n_{\rm H}\exp[-(400/T)^{2}]+1.4n({\rm H_{2}})/n_{\rm H}\exp[-1200/(T+1200)]\right\\}~{}{\rm cm^{-3}}$ |
2 |  ${\rm H}$ photoionization | $\Gamma_{2}$ | Eq. (63)
3 |  $\rm H_{2}$ photodissocian |  $\Gamma_{3}$, see the reference | 1
Cooling
1 |  ${\rm H_{2}}$ dissociation | $\Lambda_{1}=4.48\,[\,k_{\rm H4}n({\rm H_{2}})n({\rm e})+k_{\rm H5}n({\rm H_{2}})n({\rm H})+k_{\rm H8}n^{2}({\rm H_{2}})\,]\,{\rm eV}$ | 1,2
2 |  ${\rm H}$ ionization | $\Lambda_{2}=13.6k_{\rm H0}n({\rm H})n({\rm e}){\rm eV}$ | 1,2
3 |  ${\rm H}$ recombination | $\Lambda_{3}=\exp[\ln 10\times(-26.02+0.9187\log_{10}T-3.733(\log_{10}T))^{2}+0.1174(\log_{10}T))^{3}$ |
|  |  $-0.01707(\log_{10}T))^{4}+8.119\times 10^{-4}(\log_{10}T))^{5}]n({\rm H^{+}})n({\rm e})$ | 3, 4
4 |  ${\rm H^{-}}$ free-bound | $\Lambda_{4}=0.755k_{\rm H9}n({\rm H})n({\rm e})$ | 5
5 |  ${\rm H}$ excitation | $\Lambda_{5}=7.50\times 10^{-19}[1+(T/10^{5})^{1/2}]^{-1}\exp[-118348/T]n({\rm e})n({\rm H^{+}})$ | 6
6 |  ${\rm He}^{+}$ excitation | $\Lambda_{6}=5.54\times 10^{-17}T^{-0.397}[1+(T/10^{5})^{1/2}]^{-1}\exp[-473638/T]n({\rm e})n({\rm He^{+}})$ | 6
7 | Free-free | $\Lambda_{7}=1.426\times 10^{-27}T^{1/2}g_{\rm ff}(T)n(\rm H^{+})n(\rm e)$ | 7
|  |  $g_{\rm ff}(T)=0.79464+0.1243\log_{10}T$ $T<3.2\times 10^{5}~{}{\rm K}$ |
|  |  $=2.13164-0.1240\log_{10}T$ $T>3.2\times 10^{5}~{}{\rm K}$ |
8 | Compton | $\Lambda_{8}=1.017\times 10^{-37}T_{\rm CMB}^{4}(T-T_{\rm CMB})n({\rm e})$ | 6
9 | Line cooling | $\Lambda_{9}=\Lambda_{\rm H_{2}}+\Lambda_{\rm CII}+\Lambda_{\rm CO}+\Lambda_{\rm OI}+\Lambda_{\rm OII}+\Lambda_{\rm OIII}$ |
10 | Gas-grain heat transfer | $\Lambda_{10}=5.83\times 10^{-8}n_{\rm H}\rho(T/10^{3})^{1/2}[1-0.8\exp(-75/T)](T-T_{\rm d})(Z/Z_{\odot})$ | 1,2

References. (1) Hollenbach & McKee (1979) (2) Omukai (2000) (3) Ferland et al.
(1992) (4) Sugimura et al. (2017) (5) Hosokawa et al. (2016) (6) Cen (1992)
(7) Shapiro & Kang (1987)

##  Appendix C Radiation Transfer

We use the moment-based radiative transfer (RT) with M1-closure (e.g, Rosdahl
et al., 2013; Rosdahl & Teyssier, 2015; Kannan et al., 2019) for extreme
ultraviolet (EUV), far-ultraviolet (FUV), and infrared (IR) photons. In this
method, the computational cost does not depend on the number of radiative
sources. Therefore, it is suitable for calculating massive clouds where a lot
of stars are likely to form. The moment-based technique has been introduced in
recent simulation studies of the star cluster formation (e.g., Skinner &
Ostriker, 2015; Raskutti et al., 2016; Geen et al., 2017; He et al., 2019).

###  C.1 Momentum-based RT equations

In the RT calculations, we solve the specific intensity
$I_{\nu}(\bf{x},t,\nu,\bf{n})$, which is the radiative energy ($dE_{\nu}$)
crossing the point $\bf{x}$ at time $t$, per unit area $dA$, per unit time
$dt$, per unit solid angle $d\Omega$ around the direction $\bf{n}$, and per
unit frequency range $d\nu$ as

$\displaystyle dE_{\nu}=I_{\nu}dAdtd\Omega d\nu.$ (37)

The equation of RT is

$\displaystyle\frac{1}{c}\frac{\partial I_{\nu}}{\partial t}+\bm{n}\cdot\nabla
I_{\nu}=j_{\nu}-\alpha_{\nu}I_{\nu},$ (38)

where $j_{\nu}$ is the emissivity, and $\alpha_{\nu}$ is the absorption
coefficient.

We obtain the zeroth and first moment equations from integration of the RT
equation (38) over all solid angle as

$\displaystyle\frac{\partial E}{\partial t}+\nabla\cdot{\bf{F}}=S-\alpha_{\rm
E}\tilde{c}E,$ (39)
$\displaystyle\frac{1}{\tilde{c}}\frac{\partial{\bf{F}}}{\partial
t}+\tilde{c}\nabla\cdot{\bf{P}}=-\alpha_{\rm F}{\bf{F}},$ (40)

where $E$, $\bf{F}$, and $\bf{P}$ are the radiation energy density, the flux
and the radiation pressure tensor (Rosdahl et al., 2013), and $\tilde{c}$
represents the reduced light speed in the numerical simulations (also see Sec
C.2.4). In Equation (40), $S$ is the source term, $\alpha_{\rm E}$ and
$\alpha_{\rm F}$ are the energy density and flux weighted absorption
coefficients.

The radiation pressure tensor is defined with the Eddington tensor $(\bf{D})$
and the radiation energy density ($E$) as

$\displaystyle{\bf P}=E{\bf D}.$ (41)

We need to calculate the RT equation (38) to obtain the accurate values of the
Eddington tensor (the variable Eddington tensor method, Mihalas & Mihalas,
1984; Fukushima et al., 2018). However, direct calculations of Equation (38)
are computationally expensive. Thus, we adopt the M1 closure scheme in which
the Eddington tensor $\bf{D}$ is approximated as (Levermore, 1984):

$\displaystyle{\bf D}=\frac{1-\chi}{2}{\bf I}+\frac{3\chi-1}{2}{\bf
n}\otimes{\bf n},$ (42)

where

$\displaystyle{\bf n}=\frac{{\bf F}}{|{\bf
F}|},~{}\chi=\frac{3+4f^{2}}{5+2\sqrt{4-3f^{2}}},~{}{\rm and}~{}f=\frac{|{\bf
F}|}{\tilde{c}E}.$ (43)

This approximation is accurate for the optically thick or thin limits
($f\rightarrow 0$, ${\bf{D}}\rightarrow\frac{1}{3}{\bf I}$ or $f\rightarrow
1$, ${\bf{D}}\rightarrow{\bf n}\otimes{\bf n}$). However, it fails to follow
the radiation transport accurately when fluxes from different directions
collide. In this work, a lot of stellar particles distributes at the center of
cloud and most of radiation energy propagates radially. Therefore, the
colliding effect is unlikely to be serious (see also, Trebitsch et al., 2017).

We estimate the photoionization and the photodissociation rates from the
number density of EUV and FUV photons. In this case, it is convenience to
calculate transfer of the photon number density directly instead of Equation
(39) and (40) for EUV and FUV photons as (Rosdahl et al., 2013)

$\displaystyle\frac{\partial N_{\gamma}}{\partial t}+\nabla\cdot{\bf
F_{\gamma}}=\dot{N}_{\gamma,*}-\bar{\alpha}\tilde{c}N_{\gamma},$ (44)
$\displaystyle\frac{\partial{\bf F}_{\gamma}}{\partial
t}+\tilde{c}^{2}\nabla\cdot{\bf P}_{\gamma}=-\bar{\alpha}\tilde{c}{\bf
F}_{\gamma},$ (45)

where $N_{\gamma}$, ${\bf F}_{\gamma}$, and $\dot{N}_{\gamma,*}$ are the
photon number density, the photon number flux, and the photon injection rate
from a radiation source. The frequency-averaged absorption coefficient
$\bar{\alpha}$ is given as

$\displaystyle\bar{\alpha}=n_{\rm i}\bar{\sigma}_{i},$ (46)
$\displaystyle\bar{\sigma}_{i}=\frac{\int^{\nu_{2}}_{\nu_{1}}\frac{4\pi
J_{\nu}}{h\nu}\sigma_{i}(\nu)d\nu}{\int^{\nu_{2}}_{\nu_{1}}\frac{4\pi
J_{\nu}}{h\nu}d\nu},$ (47)

where $\sigma_{i}(\nu)$ is the frequency dependent cross-section, and
$J_{\nu}$ is the mean intensity:

$\displaystyle J_{\nu}=\frac{1}{4\pi}\int I_{\nu}d\Omega.$ (48)

As in Equation (41), the radiation pressure tensor ${\bf P}_{\gamma}$ is
evaluated from the M1-closure relation (see eq. 42).

###  C.2 Computational procedure for RT

We solve the moment equations on Cartesian grid. Each grid has four variables
$(E,{\bf F})$ or $(N_{\gamma},{\bf F_{\gamma}})$. According to Rosdahl et al.
(2013), we adopt the operator-splitting method to advance the time step
$(\Delta t)$. Here, we calculate the moment equations in the following three
steps.

####  C.2.1 The injection step

In this step, we only inject the photons emitted from the radiation sources as

$\displaystyle N_{\gamma}^{n+1}=N_{\gamma}^{n}+\dot{N}_{\gamma,*}\Delta t,$
(49)

where $\dot{N}_{\gamma,*}$ is the local injection rate from sources. We
estimate the injection rate with the photon emissivity $S_{\rm i}$ as

$\displaystyle\dot{N}_{\gamma,*}=\frac{S_{\rm i}}{\Delta V},$ (50)

where $\Delta V$ is the cell volume. As discussed in Section 2.2, we assume
that the photon emissivities of EUV, FUV, and IR components are proportional
to the sink mass in this study.

####  C.2.2 The transport step

In the transport step, we calculate the photon propagation without the source
and absorption terms in Equations (44) and (45) as