Making a Supermassive Star by Stellar Bombardment

We consider the following components: a central star, surrounding stars, a gas cloud, and a dark matter (DM) halo. The label "surrounding stars" refers to all stars other than the central star. We reserve the term "stars" throughout the paper to include both surrounding stars and the central star. We follow the evolution of the entire system for 3 Myr, which is roughly the time when either the first surrounding star may be expected to explode or the central star collapses to a BH.

In our semianalytical model, N-body particles represent surrounding stars. Surrounding stars form, migrate, and accrete onto the central star, while the central star is pinned to the center of the system neglecting both gas driven migration and wandering due to dynamical interactions with other stars.³ We investigate several values for the maximum initial mass of surrounding stars (${m}_{0,\max }$), and also set the initial mass of the central star to be ${m}_{\mathrm{cent}}={m}_{0,\max }$. As a fiducial model, we set ${m}_{0,\max }=1\,{M}_{\odot }$, which is roughly the Jeans mass at which fragmentation occurs at the density ${10}^{10}\,{\mathrm{cm}}^{-3}$ (Omukai et al. 2008). We assume that there are no surrounding stars initially.

When stars are expanding, we set their radius to

$$\begin{eqnarray}&&{R}_{i}=2.6\times {10}^{3}\,{R}_{\odot }{\left({m}_{i}/100{M}_{\odot }\right)}^{0.5},\end{eqnarray} \tag{ 1 }$$

following Hosokawa et al. (2012), while after stars are contracted (the condition of these stars depends on their accretion rate and the KH timescale as described in Section 3.3 below), their radius is assumed to be

$$\begin{eqnarray}&&{R}_{i}={R}_{\mathrm{ZAMS},i}=4.6\,{R}_{\odot }{\left({m}_{i}/100{M}_{\odot }\right)}^{0.58}\end{eqnarray} \tag{ 2 }$$

(Hirano & Bromm 2017). Throughout this paper, we refer to a star in the expanding (pre-main-sequence) and the contracting (main-sequence) phases as an "expanding star" and a "contracted star," respectively.

3.1.1. Density Profile

We set the number density profile for gas n_gas(r) as

$$\begin{eqnarray}{n}_{\mathrm{gas}}(r)=\left\{\begin{array}{cl}{n}_{{\rm{c}}} & \mathrm{if}\,r\leqslant {r}_{{\rm{c}}}\\ {n}_{{\rm{c}}}{\left(r/{r}_{{\rm{c}}}\right)}^{-2} & \mathrm{if}\,{r}_{{\rm{c}}}\leqslant r\leqslant {r}_{\mathrm{vir}}\\ 0 & \mathrm{if}\,r\gt {r}_{\mathrm{vir}}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where r_c is the core radius of the collapsing gas and n_c is the core gas density. Outside the core radius, gas is assumed to be collapsing under its self-gravity, while in the dense core, gas is assumed to cool efficiently, fragment, and form stars. As a fiducial model, we set n_c = 10¹¹ cm⁻³, and the temperature of inflowing gas to T = 10⁴ K (e.g., Oh & Haiman 2002; Omukai et al. 2008; Shang et al. 2010). Assuming an isothermal equation of state, the sonic velocity of inflowing gas is ${c}_{{\rm{s}}}={({kT}/\mu )}^{1/2}\simeq 10\,\mathrm{km}\,{{\rm{s}}}^{-1}{(T/{10}^{4}{\rm{K}})}^{1/2}$, where k is the Boltzmann constant and μ = 1.22 is the mean mass per particle, and the accretion rate from large scales is set to

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=\displaystyle \frac{{c}_{{\rm{s}}}^{3}}{G}\simeq 0.22\,\displaystyle \frac{{M}_{\odot }}{\mathrm{yr}}{\left(\displaystyle \frac{T}{{10}^{4}{\rm{K}}}\right)}^{3/2}\end{eqnarray} \tag{ 4 }$$

(e.g., Stahler et al. 1980; Begelman et al. 2006). The core density ${n}_{{\rm{c}}}={10}^{11}\,{\mathrm{cm}}^{-3}$ roughly matches the density at which gas with a metallicity of ∼10⁻⁵ Z_⊙ and a suppressed H₂ fraction (by strong background radiation) begins to fragment due to the decrease of its temperature (Omukai et al. 2008). In the fiducial model, the initial value of the core radius for the collapsing gas is chosen to be ${r}_{{\rm{c}},\mathrm{ini}}\simeq 4\times {10}^{-4}$ pc by matching the cosmological baryon-to-DM mass ratio inside the virial radius ${r}_{\mathrm{vir}}={(3{M}_{\mathrm{halo}}/800\pi {\rho }_{\mathrm{cri}})}^{1/3}$, where M_halo is the halo mass within the virial radius, ${\rho }_{\mathrm{cri}}=3{H}^{2}/8\pi G$ is the cosmological critical density, $H\simeq {H}_{0}{[{{\rm{\Omega }}}_{{\rm{m}}0}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}0}]}^{1/2}$ is the Hubble parameter, ${H}_{0}\simeq 70\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ is the Hubble constant, Ω_m0 = 0.24 is the matter density today, and Ω _Λ0 = 0.76 is the cosmological constant today (Planck Collaboration et al. 2016). We assume that the halo mass within the virial radius M_halo is 10⁷ ${M}_{\odot }$. The radius at which ${n}_{\mathrm{gas}}(r)={10}^{11}\,{\mathrm{cm}}^{-3}$, measured in high-resolution cosmological hydrodynamical simulations of metal- and H₂-free gas in atomic-cooling halos (e.g., Regan et al. 2014) is also found to be ∼10⁻⁴–10⁻³ pc. We also checked that our results are not significantly influenced by changing the value of ${r}_{{\rm{c}},\mathrm{ini}}$ to 10⁻³ pc, which is because r_c quickly evolves due to our assumption of setting r_c to the place where the gas becomes unstable to fragmentation (see Section 3.2). Thus, we assume that the core gas density is fixed while r_c evolves with time (Section 3.2). In our model, the final results depend on the position of star formation at r_c, while they are less affected by other effects related to the gas density distribution.⁴ Therefore we expect that the core gas density fixed in time is not a critical assumption.

When we calculate the acceleration due to gas dynamical friction and accretion (Sections 3.4.2 and 3.6.2), we assume that the gas mean velocity is zero for simplicity. This assumption gives an optimal rate for migration toward the center for surrounding stars due to gas dynamical friction and accretion. Nevertheless, the migration by gas dynamical friction and accretion is found to give small contributions to the final SMS mass (Section 4.1).

3.1.2. Gravitational Potential

We adopt a four-component gravitational potential,

$$\begin{eqnarray}&&{\rm{\Phi }}(r)={{\rm{\Phi }}}_{\mathrm{DM}}(r)+{{\rm{\Phi }}}_{\mathrm{star}}(r)+{{\rm{\Phi }}}_{\mathrm{gas}}(r)+{{\rm{\Phi }}}_{\mathrm{cent}}(r),\end{eqnarray} \tag{ 5 }$$

where ${{\rm{\Phi }}}_{\mathrm{DM}}(r)$, ${{\rm{\Phi }}}_{\mathrm{star}}(r)$, ${{\rm{\Phi }}}_{\mathrm{gas}}(r)$, and ${{\rm{\Phi }}}_{\mathrm{cent}}(r)$ are, respectively, the gravitational potential at the position r of the dark matter halo, surrounding stars, collapsing gas, and the central star.

We set ${{\rm{\Phi }}}_{\mathrm{DM}}(r)$ by the Navarro–Frenk–White (NFW) profile as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{DM}}(r)=\displaystyle \frac{-4\pi G{\rho }_{{\rm{h}}}{r}_{{\rm{h}}}^{3}}{r}\mathrm{ln}\left(\displaystyle \frac{{r}_{{\rm{h}}}+r}{{r}_{{\rm{h}}}}\right),\end{eqnarray} \tag{ 6 }$$

where ${r}_{{\rm{h}}}={r}_{\mathrm{vir}}/C$ is the scale radius of the halo, ${\rho }_{{\rm{h}}}=200{\rho }_{\mathrm{cri}}{C}^{3}/(3\mathrm{ln}(1+C)-C/(1+C))$ is the density parameter of the NFW profile, and C is the concentration parameter (Navarro et al. 1997). We assume that C = 9. We assume a redshift z = 15, since atomic-cooling halos (whose masses are ≈10⁷ ${M}_{\odot }$ at this redshift) start to appear from around this epoch (e.g., Tanaka & Haiman 2009). For reference, we also note that in the context of trying to grow to an SMBH of ∼10⁹ ${M}_{\odot }$ at z ∼ 6 via gas accretion, the seed BH mass is required to be ∼10⁵ ${M}_{\odot }$ at z ∼ 10 (Di Matteo et al. 2012).

The gravitational potential of the gas is derived from Equation (3) using Equation (2.28) in Binney & Tremaine (2008) as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{gas}}(r)=4\pi {{Gn}}_{{\rm{c}}}\mu {m}_{{\rm{H}}}\times \left\{\begin{array}{l}\displaystyle \frac{{r}^{2}}{6}-\tfrac{{r}_{{\rm{c}}}^{2}}{2}-{r}_{{\rm{c}}}^{2}\mathrm{ln}\tfrac{{r}_{\mathrm{vir}}}{{r}_{{\rm{c}}}}\\ \qquad \qquad \mathrm{for}\,r\lt {r}_{{\rm{c}}},\\ \tfrac{2{r}_{{\rm{c}}}^{3}}{3r}-{r}_{{\rm{c}}}^{2}-{r}_{{\rm{c}}}^{2}\mathrm{ln}\tfrac{{r}_{\mathrm{vir}}}{r}\\ \qquad \qquad \mathrm{for}\,r\gt {r}_{{\rm{c}}},\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where m_H is the hydrogen mass.

Similarly, for the gravitational potential of the surrounding stars:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{star}}(r)=-\displaystyle \frac{{{GM}}_{* }(\lt r)}{r}-{\int }_{r}^{{R}_{\max }}4\pi {rG}{\rho }_{* }(r){dr},\end{eqnarray} \tag{ 8 }$$

where ρ_*(r) is the stellar density at r, and M_*(<r) is the stellar mass within r. We integrate and derive the gravitational potential at the center of the spherical radial cell l in each time step using a linear interpolation of Φ_star(r). We use 60 radial cells, covering the range from 10⁻⁸ pc to R_max = 10 pc, spaced uniformly in log r.

The gravitational potential by the central star is ${{\rm{\Phi }}}_{\mathrm{cent}}(r)=-{{Gm}}_{\mathrm{cent}}/r$.

Note that the profiles of Φ_DM(r), Φ_gas(r), and Φ_cent(r) are assumed to be algebraically fixed as given above, while r_c and the normalizations of Φ_gas(r) and Φ_cent(r) are followed in time. Φ_star(r) is assumed to evolve, and its full radial profile is followed during our calculations.

We assume that any gas flowing in from large scales that is not accreted onto stars is converted into new surrounding stars, and thus the star formation rate is

$$\begin{eqnarray}&&{\dot{m}}_{\mathrm{SF}}={\dot{M}}_{\mathrm{in}}-\sum _{i}{\dot{m}}_{i}\end{eqnarray} \tag{ 9 }$$

where the sum goes over all stars, including the central star and surrounding stars. This assumption is not obvious, since the inflowing gas may accumulate in the central region and increase the core density and core radius. Since it is difficult to follow the time evolution of the core density due to this gas accumulation, we assume a constant core density and a constant inflow rate as input parameters. Due to the temperature dependence of the inflow rate given by Equation (4), cases with low star formation rates are investigated in models with low T below.

Following simulations of Population II star formation by Dopcke et al. (2013), we set the initial mass function to

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dm}}\propto {m}^{-\beta }\,\,\mathrm{where}\,\,0.08{M}_{\odot }\leqslant m\leqslant {m}_{0,\max }\end{eqnarray} \tag{ 10 }$$

with a flat logarithmic distribution, β = 0, as a fiducial model. In practice, in each time step we form new stars randomly from this mass function in succession as long as their total mass is less than ${m}_{\mathrm{SF}}={\dot{m}}_{\mathrm{SF}}{\rm{\Delta }}t$. Let us label with n the corresponding largest number of new stars where this criterion holds. Finally we form the last star with probability ${m}_{\mathrm{new},n+1}/({m}_{\mathrm{SF}}-{{\rm{\Sigma }}}_{i=1}^{n}{m}_{\mathrm{new},i})$, where the mass of the n + 1th newly formed star ${m}_{\mathrm{new},n+1}$ is drawn randomly from the mass function.

We assign cells to the newly formed stars based on the following arguments. Although we assume that the collapsing gas cloud is overall spherically symmetric, gas disks may form around the central star or around surrounding stars forming and growing by accreting gas with nonzero angular momentum. Since a gas disk is stabilized by rapid rotation in a steep gravitational potential, surrounding stars form outside the radius where the Toomre parameter

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{{\rm{\Sigma }}}_{\mathrm{gas}}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{Kep}}{\rm{\Omega }}}{2\pi {{Gn}}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}}\end{eqnarray} \tag{ 11 }$$

becomes 1, where Ω is the orbital frequency of the gas disk, and ${{\rm{\Omega }}}_{\mathrm{Kep}}{(r)=({GM}(r)/{r}^{3})}^{1/2}$ is the Keplerian orbital frequency, where M(r) is the enclosed mass, and Σ_gas is the surface density of the gaseous disk. Since gas disks are partially supported also via turbulent and thermal pressure, ${\rm{\Omega }}={\epsilon }_{\mathrm{Kep}}{{\rm{\Omega }}}_{\mathrm{Kep}}$, where _Kep describes the deviation from a fully rotationally supported disk. Referring to the results of simulations for primordial disks (e.g., Greif et al. 2012; Latif et al. 2013; Hirano et al. 2014), we adopt _Kep = 0.5. In the case of ${r}_{Q=1}\gt {r}_{{\rm{c}},\mathrm{ini}}$, unless the total accretion rate is as high as $\sim {\dot{M}}_{\mathrm{in}}$, gas accumulates within the central region. Unlike our simplified assumption of a homogeneous n_gas, a more realistic density distribution for Keplerian rotating gas is $\rho \propto {r}^{-1/2}-{r}^{-3/2}$ (e.g., Inayoshi et al. 2018b). However, even for this density profile within the core, gas is most unstable in the outer region where n_gas ≥ n_c. Thus we assume that surrounding stars form at r_Q=1 when ${r}_{Q=1}\gt {r}_{{\rm{c}},\mathrm{ini}}$, while surrounding stars form uniformly from r_Q=1 to ${r}_{{\rm{c}},\mathrm{ini}}$ when ${r}_{Q=1}\lt {r}_{{\rm{c}},\mathrm{ini}}$. Accordingly the core radius is set to r_c = r_Q=1 in each time step using Equation (11). From Equation (11), r_c = r_Q=1 is satisfied at

$$\begin{eqnarray}&&{r}_{{\rm{c}}}={\left[\displaystyle \frac{{\epsilon }_{\mathrm{Kep}}}{\left(\tfrac{3}{2}-{\epsilon }_{\mathrm{Kep}}\right)}\displaystyle \frac{[{m}_{\mathrm{cent}}+{M}_{\mathrm{stars}}({r}_{{\rm{c}}})+{M}_{\mathrm{DM}}({r}_{{\rm{c}}})]}{\tfrac{4}{3}\pi {n}_{{\rm{c}}}\mu {m}_{{\rm{H}}}}\right]}^{1/3}.\end{eqnarray} \tag{ 12 }$$

Let us introduce M_DM(r) and M_stars(r) to label the enclosed mass of the dark matter and surrounding stars within r, respectively. The dark matter mass is typically subdominant in this expression. In the early phases, the total stellar mass of surrounding stars and the central star is limited by the inflow rate from large scales and ${m}_{\mathrm{cent}}+{M}_{\mathrm{stars}}(r)\sim {\dot{M}}_{\mathrm{in}}t$, which implies that

$$\begin{eqnarray}&&{r}_{{\rm{c}}}\approx \displaystyle \frac{{\epsilon }_{\mathrm{Kep}}^{1/3}}{{\left(\tfrac{3}{2}-{\epsilon }_{\mathrm{Kep}}\right)}^{1/3}}\displaystyle \frac{{\dot{M}}_{\mathrm{in}}^{1/3}{t}^{1/3}}{{\left[\tfrac{4}{3}\pi {n}_{{\rm{c}}}\mu {m}_{{\rm{H}}}\right]}^{1/3}}.\end{eqnarray} \tag{ 13 }$$

We ignore the dynamical effect on the surrounding stars due to the deepening of the gas gravitational potential if r_c moves outwards, which tightens the orbit of surrounding stars outside of r_c.

In order to form stars, cooling from dust grains needs to be stronger than heating, since gas fragmentation is caused by the decrease of the gas temperature due to cooling by dust grains (Omukai et al. 2008). In our model, gas is heated by gas dynamical friction, with the total heating rate given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{GDF}}=\displaystyle \frac{{{dE}}_{\mathrm{GDF}}}{{dt}}=\sum _{i}{m}_{i}{v}_{i}{a}_{\mathrm{GDF},i}.\end{eqnarray} \tag{ 14 }$$

We assume that even when gas dynamical friction does not reduce the velocity of surrounding stars at ${{\rm{\Sigma }}}_{i\in l}4\pi {r}_{\mathrm{Bondi},i}^{3}/3\gt {V}_{l}$ (Section 3.4.2), gas is heated by gas dynamical friction, and we substitute ${a}_{\mathrm{GDF},i}$ calculated by Equation (24) into Equation (14). Thus Equation (14) represents the upper limit for the heating rate due to gas dynamical friction.

Referring to Omukai et al. (2008), the specific cooling rate from dust grains is

$$\begin{eqnarray}&&{\lambda }_{\mathrm{dust}}\sim \left\{\begin{array}{l}100\,\tfrac{\mathrm{erg}}{{\rm{s}}\,{\rm{g}}}\tfrac{{n}_{\mathrm{gas}}}{{10}^{10}\,{\mathrm{cm}}^{-3}}\\ \qquad \mathrm{for}\,{10}^{6}\,{\mathrm{cm}}^{-3}\lt {n}_{\mathrm{gas}}\lt {10}^{10}\,{\mathrm{cm}}^{-3},\\ 100\,\tfrac{\mathrm{erg}}{{\rm{s}}\,{\rm{g}}}{\left(\tfrac{{n}_{\mathrm{gas}}}{{10}^{10}{\mathrm{cm}}^{-3}}\right)}^{0.2}\\ \qquad \mathrm{for}\,{10}^{10}\,{\mathrm{cm}}^{-3}\lt {n}_{\mathrm{gas}}\lt {10}^{12}\,{\mathrm{cm}}^{-3}.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

Since the core region is optically thin to dust emission (Omukai et al. 2008), the net cooling rate by dust grains within the core region is ${{\rm{\Lambda }}}_{\mathrm{dust}}={\lambda }_{\mathrm{dust}}{M}_{{\rm{c}}}$, where ${M}_{{\rm{c}}}=\tfrac{4\pi }{3}{n}_{{\rm{c}}}\mu {m}_{{\rm{H}}}{r}_{{\rm{c}}}^{3}$ is the gas mass within the core radius. Whenever the heating rate exceeds the cooling rate, Γ_GDF > Λ_dust, star formation is assumed to be quenched. When the cooling dominates the heating, the gas temperature is expected to decrease and gas fragments as found in Omukai et al. (2008).

For each star, we specify whether they are in the expanding (pre-main-sequence) or contracting (main-sequence) phase in each simulation time step as follows. In isolation, a protostar contracts on the KH timescale ${t}_{\mathrm{KH},i}={{Gm}}_{i}^{2}/({R}_{i}{L}_{i})$. On the other hand, Hosokawa et al. (2012) have shown that if the mass accretion rate (see Section 3.6.2) exceeds the critical rate, ${\dot{m}}_{i}\gt {\dot{m}}_{\mathrm{cri}}\sim 0.006\mbox{--}0.03\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, stars can keep expanding (see also Haemmerlé et al. 2018 for similar results). Furthermore, Sakurai et al. (2015) have shown that if the mass accretion rate time-averaged over a KH timescale evaluated on the stellar surface, ${t}_{\mathrm{surf},\mathrm{KH},i}\sim 10\,{t}_{\mathrm{KH},i}\sim 10\,{{Gm}}_{i}^{2}/{R}_{i}{L}_{i}$, exceeds the critical rate $\langle {\dot{m}}_{i}\rangle \gt {\dot{m}}_{\mathrm{cri}}$, the star can keep expanding. Following these results, we assume that if the accretion rate ${\dot{m}}_{i}$ smoothed on the surface KH timescale⁵ exceeds ${\dot{m}}_{\mathrm{cri}}=0.01\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, or the time from its formation is shorter than the KH timescale for zero-age main-sequence (ZAMS) stars (${t}_{\mathrm{KH},\mathrm{ZAMS},i}={{Gm}}_{i}^{2}/({R}_{\mathrm{ZAMS},i}{L}_{i})$), star i keeps expanding, and otherwise it contracts. We calculate the surface KH timescale ${t}_{\mathrm{surf},\mathrm{KH},i}$ using the stellar luminosity L_i given by Equations (3), (4), and (5) in Hosokawa et al. (2012) for stars with masses of m_i < 6 ${M}_{\odot }$, 6 ${M}_{\odot }$ ≤ m_i < 50 ${M}_{\odot }$, and m_i ≥ 50 ${M}_{\odot }$, respectively. Note that the critical condition (${\dot{m}}_{\mathrm{cri}}$) could be lower for accretion of stars than that for accretion of gas. This is because high-velocity accretion of highly eccentric stars and/or the internal energy of accreted stars may heat the envelope of the central star more, per unit infalling mass, compared to accreting cold gas at the same rate.

A particle i (i.e., one of the surrounding stars) is described by its mass m_i and its radial distance from the central star r_i. For simplicity, particles are assumed to follow circular orbits, but they are allowed to migrate radially. After each time step Δt, we update the position of particle i to r_i + Δr_i, satisfying

$$\begin{eqnarray}&&E({r}_{i}+{\rm{\Delta }}{r}_{i})=E({r}_{i})+{\rm{\Delta }}{k}_{i}\end{eqnarray} \tag{ 17 }$$

where $E({r}_{i})={\rm{\Phi }}({r}_{i})+k({r}_{i})$ is the total specific energy, $k({r}_{i})=\tfrac{1}{2}{v}_{i}^{2}$ is the specific kinetic energy, ${\rm{\Delta }}{k}_{i}$ is the change in the specific kinetic energy within Δt, and v_i is the orbital velocity of the ith particle. The change in the specific kinetic energy is given as

$$\begin{eqnarray}&&{\rm{\Delta }}{k}_{i}={v}_{i}{a}_{i}{\rm{\Delta }}t\end{eqnarray} \tag{ 18 }$$

where a_i is the acceleration of the ith particle. We assume ${v}_{i}={v}_{\mathrm{Kep}}({r}_{i})$, where ${v}_{\mathrm{Kep}}(r)$ is the Keplerian orbital velocity at r. The acceleration is given as

$$\begin{eqnarray}&&{a}_{i}={a}_{\mathrm{SDF},i}+{a}_{\mathrm{GDF},i}+{a}_{\mathrm{acc},i}\end{eqnarray} \tag{ 19 }$$

where ${a}_{\mathrm{SDF},i}$, ${a}_{\mathrm{GDF},i}$, and ${a}_{\mathrm{acc},i}$ are the acceleration of the ith particle due to stellar dynamical friction, gas dynamical friction, and accretion, respectively (see Sections 3.4.1, 3.4.2, and 3.6.2 below).

For simplicity, in the calculation of the migration rate in Equation (17), we assume that the eccentricity of a surrounding star does not evolve with time, and remains zero. We consider the effects of nonzero eccentricities in Section 4.3.

To follow the migration, we use a shared time step of

$$\begin{eqnarray}&&{\rm{\Delta }}t=\eta \,{\min }_{i}\left[\displaystyle \frac{{v}_{i}}{{a}_{\mathrm{SDF},i}},\displaystyle \frac{{v}_{i}}{{a}_{\mathrm{GDF},i}},\displaystyle \frac{{v}_{i}}{{a}_{\mathrm{acc},i}},\displaystyle \frac{1}{\sqrt{{{Gn}}_{{\rm{c}}}\mu {m}_{{\rm{H}}}}}\right]\end{eqnarray} \tag{ 20 }$$

where the constant η is a time step parameter. On the right-hand side, the four terms are the timescales for stellar dynamical friction, gas dynamical friction, and accretion torque, and the dynamical time within the core radius, respectively.

We set the time step parameter to be η = 0.1. To validate this choice, we compared the final mass of the central star in one of the models ("Model 2" below) with η = 0.4, 0.2, and 0.1. The mass was found to be $4.1\times {10}^{3}\,{M}_{\odot }$, $3.8\times {10}^{3}\,{M}_{\odot }$, and $3.7\times {10}^{3}\,{M}_{\odot }$, respectively. The final mass changes only by <2% between the last two cases (η = 0.2 and η = 0.1), giving us confidence that our results have nearly converged at η = 0.1.

3.4.1. Stellar Dynamical Friction

Stellar dynamical friction is modeled using the analytic formula (Binney & Tremaine 2008) of

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{SDF},i} & = & -\displaystyle \frac{4\pi {G}^{2}{m}_{i}{\rho }_{* }\mathrm{ln}{\rm{\Lambda }}}{{v}_{i}^{2}}\left[\mathrm{erf}\left(\displaystyle \frac{{v}_{i}}{\sqrt{2}{\sigma }_{* }}\right)\right.\\ & & -\,\left.\displaystyle \frac{\sqrt{2}{v}_{i}}{\sqrt{\pi }{\sigma }_{* }}{e}^{-{v}_{i}^{2}/2{\sigma }_{* }^{2}}\right],\end{array}\end{eqnarray} \tag{ 21 }$$

where ${\sigma }_{* }$ is the velocity dispersion of background stars, ${\rho }_{* }$ is the stellar density, $\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}({b}_{\max }/{b}_{\min })$ is the Coulomb logarithm, and b_max and b_min are the maximum and minimum impact parameters for weak stellar encounters. We set ${b}_{\min }={{Gm}}_{i}/{v}_{i}^{2}$, ${b}_{\max }=0.1$ pc. Equation (21) assumes an isotropic and Maxwellian velocity distribution for background stars. We set ${\sigma }_{* }={v}_{\mathrm{Kep}}({r}_{i})/\sqrt{3}$, which sets the value in the square parenthesis in Equation (21) to 0.86 since we assume ${v}_{i}={v}_{\mathrm{Kep}}({r}_{i})$.

To obtain the background density ${\rho }_{* }$ at each time step, we compute an average stellar density in each radial cell $l\in [1,60]$, obtained from the total number of surrounding stars found in each cell assuming spherical symmetry. To check the effect of the number of cells N_cell, we compared the results for Model 2 for N_cell = 40, 60, and 80. The final mass of the central star in these three cases was found to be 4.1 × 10³, 3.7 × 10³, and 3.6 × 10³, respectively. The small difference (<3%) between the latter two cases gives us confidence that our results nearly converge for N_cell = 60.

When m_i is larger than the average mass (${\overline{m}}_{l}$) in some cell l hosting the ith surrounding star, the ith surrounding star migrates inward by the acceleration in Equation (21). On the other hand, when ${m}_{i}\lt {\overline{m}}_{l}$, the ith surrounding star gains kinetic energy from the encounter and migrates outward, which is not accounted for by Equation (21). Due to energy conservation, the total kinetic energy change for all surrounding stars in each cell by stellar dynamical friction is zero. To reduce computational time, we assign equal momentum change (Δp_l) to every below-average surrounding star in each cell. Here Δp_l is determined from energy conservation by solving the following equation in each cell

$$\begin{eqnarray}&&\begin{array}{l}\sum _{{m}_{i}\lt {\overline{m}}_{l}}\displaystyle \frac{1}{2}{m}_{i}\left[{\left({v}_{i}+\displaystyle \frac{{\rm{\Delta }}{p}_{l}}{{m}_{i}}\right)}^{2}-{v}_{i}^{2}\right]\\ \quad =\,-\sum _{{m}_{i}\gt {\overline{m}}_{l}}\displaystyle \frac{1}{2}{m}_{i}\left[{\left({v}_{i}-| {a}_{\mathrm{SDF},i}| {\rm{\Delta }}t\right)}^{2}-{v}_{i}^{2}\right].\end{array}\end{eqnarray} \tag{ 22 }$$

Since ${v}_{i}\gg | {a}_{\mathrm{SDF},i}| {\rm{\Delta }}t$ (Equation (20)) and Δp_l ≪ m_iv_i, we approximate this equation as

$$\begin{eqnarray}&&\sum _{{m}_{i}\lt {\overline{m}}_{l}}{v}_{i}{\rm{\Delta }}{p}_{l}=\sum _{{m}_{i}\gt {\overline{m}}_{l}}{m}_{i}{v}_{i}| {a}_{\mathrm{SDF},i}| {\rm{\Delta }}t.\end{eqnarray} \tag{ 23 }$$

We assign the new radial location to the surrounding stars to match the updated velocity to the circular velocity at that radius (Section 3.4). This procedure ensures that the cluster is in local virial equilibrium everywhere and accounts for two-body relaxation for the stellar cluster in an approximate way.

We assume that stellar dynamical friction operates when the number of surrounding stars within a cell is more than one. In the fiducial model, we verify that the number of surrounding stars within 10 R_cent is more than a hundred at t = 10⁴ yr. Hence the number of surrounding stars is mostly large enough to validate Equation (21) in our models.

3.4.2. Gas Dynamical Friction

When a particle has a nonzero velocity relative to the background gas, it suffers additional dynamical friction from the gas component. Due to this mechanism, surrounding stars may migrate toward the center. We use the gas dynamical friction formulation derived by Ostriker (1999) as

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{GDF},i} & = & -\displaystyle \frac{4\pi {G}^{2}{m}_{i}{n}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}}{{v}_{\mathrm{rel},i}^{2}}f({v}_{\mathrm{rel},i}{/c}_{s}),\\ f(x) & = & \left\{\begin{array}{ll}\tfrac{1}{2}\mathrm{ln}\left(\tfrac{1+x}{1-x}\right)-x & \mathrm{for}\,\,0\lt x\lt 1,\\ \tfrac{1}{2}\,\mathrm{ln}\left({x}^{2}-1\right)+\mathrm{ln}{\rm{\Lambda }}^{\prime} & \mathrm{for}\,\,x\gt 1,\end{array}\right.\end{array}\end{eqnarray} \tag{ 24 }$$

where $\mathrm{ln}\,{\rm{\Lambda }}^{\prime}$ is a Coulomb logarithm for the gas distribution, and ${v}_{\mathrm{rel},i}$ is the relative velocity between the $i\mathrm{th}$ surrounding star and the background gas. Referring to the result of numerical simulations by Chapon et al. (2013), we adopt $\mathrm{ln}{\rm{\Lambda }}^{\prime} =3.1$. We set ${v}_{\mathrm{rel},i}={v}_{i}$ assuming a static background gas distribution.

In the usual formulation of dynamical friction, a body is assumed to be moving on a straight line (but see Kim & Kim 2007; Chapon et al. 2013), relative to an unperturbed background. Since each star disturbs the gas inside its Bondi–Hoyle–Lyttleton sphere, the formulation of gas dynamical friction is not valid within another star's Bondi–Hoyle–Lyttleton sphere. We assume that when the sum of the volumes of the Bondi–Hoyle–Lyttleton spheres of surrounding stars within the spherical cell l (${{\rm{\Sigma }}}_{i\in l}4\pi {R}_{\mathrm{BHL},i}^{3}/3$) exceeds the volume of the cell V_l, gas dynamical friction does not operate in that cell. We likewise neglect gas dynamical friction inside the Bondi radius of the central star.

3.4.3. Accretion Torque

Due to gas accretion (Section 3.6.2), accreted objects receive momentum to satisfy momentum conservation. In this study, we set the acceleration due to gas accretion as

$$\begin{eqnarray}&&{a}_{\mathrm{acc},i}=-\displaystyle \frac{{\dot{m}}_{i}{v}_{i}}{{m}_{i}}.\end{eqnarray} \tag{ 25 }$$

For simplicity we assume that gas is static, and the relative velocity between surrounding stars and gas is always given by the velocity of surrounding stars, i.e., the angular momentum is always reduced, which leads to radially inward migration. Since the collapsing gas may instead have angular momentum in the same sense as the surrounding stars, this prescription gives an upper limit for the deceleration and the resulting radial migration rate for surrounding stars. In Section 4.1, we find that the deceleration by gas accretion has a minor effect on the migration of surrounding stars, even at this upper limit.

We also consider collisions among surrounding stars. Assuming that the surrounding stars' motion is isotropic, the number, the number density, and the velocity dispersion in cell l are N_l, n_l, and ${\sigma }_{* ,l}$, respectively, the expected rate of collisions within the time step Δt in a cell l is given by (Equation (7.194) in Binney & Tremaine 2008)

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{coll},l} & = & \displaystyle \frac{1}{2}{N}_{l}\displaystyle \frac{{\rm{\Delta }}t}{{t}_{\mathrm{coll},l}}\\ & = & \displaystyle \frac{1}{2}{N}_{l}{n}_{* ,l}{\sigma }_{* ,l}\left({R}_{\mathrm{coll}}^{2}+\displaystyle \frac{G{\overline{m}}_{l}}{{\sigma }_{* ,l}^{2}}{R}_{\mathrm{coll}}\right){\rm{\Delta }}t,\end{array}\end{eqnarray} \tag{ 26 }$$

where R_coll is the pericenter distance between the center of mass of two stars needed for a collision, i.e., the sum of the radii of the colliding stars, ${\overline{m}}_{l}$ is the average stellar mass in cell l, ${t}_{\mathrm{coll},l}$ is the collision timescale in cell l, and the factor 1/2 is introduced to prevent double counting due to the fact that two stars participate in the collisions. In practice, we assume that the surrounding star i collides in the simulation with probability

$$\begin{eqnarray}&&{P}_{\mathrm{coll},i}=2\sqrt{\pi }{n}_{* ,l}{\sigma }_{* ,l}\left({R}_{\mathrm{coll},i}^{2}+\displaystyle \frac{{{Gm}}_{i}}{{\sigma }_{* ,l}^{2}}{R}_{\mathrm{coll},i}\right){\rm{\Delta }}t\end{eqnarray} \tag{ 27 }$$

during a time step, where the collision radius is approximated by ${R}_{\mathrm{coll},i}=2{R}_{i}$. For describing collisions between two surrounding stars i and j, we assume that the relative velocity ${v}_{\mathrm{rel},* }$ is drawn from a Maxwellian distribution with the dispersion of $\sqrt{2}{\sigma }_{* }$ as given in Equation (8.45) in Binney & Tremaine (2008). Collisions may occur when the number of surrounding stars within a cell is more than one.

For collisions among contracted stars, when this relative velocity ${v}_{\mathrm{rel},* }$ exceeds the escape velocity from the stars, ${v}_{\mathrm{esc}}={[2G({m}_{i}+{m}_{j})/({R}_{i}+{R}_{j})]}^{1/2}$, contracted stars lose a significant amount of mass at collision instead of simply coalescing into one remnant star (Freitag & Benz 2005). For simplicity, we assume that when ${v}_{\mathrm{rel},* }\gt {v}_{\mathrm{esc}}$ for contracted stars, the colliding surrounding stars are completely disrupted. However, the fraction of the released gas mass that accretes onto the central star and that is converted to form new stars is not well understood. In this study, the mass released during collisions is added to the inflowing gas ${\dot{M}}_{\mathrm{in}}$ from large scales (Equation (4)). The inflowing gas is mostly converted to new surrounding stars during the early phase of the evolution (see Section 3.2) and it is mostly accreted onto the central star when ${\dot{M}}_{\mathrm{in}}\sim {\dot{m}}_{\mathrm{cent}}$ (see Section 3.6.2). For collisions with ${v}_{\mathrm{rel},* }\lt {v}_{\mathrm{esc}}$ between contracted surrounding stars, we assume that the stars coalesce without any mass loss. When two stars i and j coalesce, we assume that m_j accretes onto m_i. The merger remnant star may either become an expanding star or a contracted star depending on the time-smoothed accretion rate as defined in Section 3.3.

When surrounding stars are in an expanding phase (conditions specified in Section 3.2), collisions become more frequent (Boekholt et al. 2018). The mass loss during such collisions is also significantly different from that of contracted stars. Alister Seguel et al. (2020) have recently investigated the effect of mass loss during collisions on mass growth, relevant to collisions between contracted stars, in which high mass-loss rates are predicted.

We used the fraction of total mass lost during collisions among expanding stars from Figure 8 of Adams et al. (2004) as

$$\begin{eqnarray}&&{f}_{\mathrm{loss}}={f}_{\mathrm{loss},\max }\,{10}^{-1.2\left(\tfrac{{R}_{{\rm{p}}}}{{R}_{\mathrm{coll},i}}\right)},\end{eqnarray} \tag{ 28 }$$

where R_p is the pericenter distance at collision. Referring to Adams et al. (2004), we set ${f}_{\mathrm{loss},\max }\,=\,0.16$ as a fiducial value, which is roughly consistent with the results by Bailey & Davies (1999). Note that since Adams et al. (2004) simulated collisions between an expanding star and a contracted star, f_loss for collisions between expanding stars may become lower than that in Equation (28). To see the effect of ${f}_{\mathrm{loss},\max }$ on our results, we compared the final mass of the central star in one of the models ("Model 2" below) with ${f}_{\mathrm{loss},\max }\,=\,0.16$, 0.3, and 1. The final mass of the central star was found to be $3.3\times {10}^{3}\,{M}_{\odot }$, $2.6\times {10}^{3}\,{M}_{\odot }$, and $1.8\times {10}^{3}\,{M}_{\odot }$, respectively, and the total mass lost at collisions was $1.6\times {10}^{3}\,{M}_{\odot }$, $2.9\times {10}^{3}\,{M}_{\odot }$, and $4.0\times {10}^{3}\,{M}_{\odot }$, respectively. Thus the final mass of the central star is affected by at most a factor ∼2 due to ${f}_{\mathrm{loss},\max }$.

The pericenter distance is related to the impact parameter b through (e.g., O'Leary et al. 2009)

$$\begin{eqnarray}&&{R}_{{\rm{p}}}={\left(\sqrt{\displaystyle \frac{1}{{b}^{2}}+\displaystyle \frac{{G}^{2}{\left({m}_{i}+{m}_{j}\right)}^{2}}{{b}^{4}{v}_{\mathrm{rel},* }^{4}}}+\displaystyle \frac{G({m}_{i}+{m}_{j})}{{b}^{2}{v}_{\mathrm{rel},* }^{2}}\right)}^{-1}.\end{eqnarray} \tag{ 29 }$$

We set the distribution of ${R}_{{\rm{p}}}$ so that b² is uniformly distributed between 0 and ${b}_{\max }^{2}$, where b_max is the maximum impact parameter at which collision occurs (b = b_max at ${R}_{{\rm{p}}}={R}_{\mathrm{coll}}$). The fraction of mass ${f}_{\mathrm{loss}}$ (Equation (28)) is subtracted from the mass of the collided surrounding stars, and added to the inflow rate ${\dot{M}}_{\mathrm{in}}$. Following Bailey & Davies (1999), we also set the condition for merger into a single star to

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{coll}}}\lt -0.75\left(\displaystyle \frac{{v}_{\mathrm{rel},* }}{{v}_{\mathrm{esc}}}\right)+1.3.\end{eqnarray} \tag{ 30 }$$

Even when the colliding surrounding stars merge into one single star, the collision leads to some mass loss in the case where at least one of the two colliding stars is expanding, according to Equation (28).

3.6.1. Stellar Accretion

3.6.2. Gas Accretion

Inayoshi et al. (2018b) have considered radiatively inefficient accretion onto a compact object and generalized Bondi accretion to a case with angular momentum. They have found that when the angular momentum is low, so that the centrifugal radius is well inside the Bondi radius and a compact accretion disk forms around the central object, the accretion rate ${\dot{m}}_{\mathrm{acc},i}$ onto the central object (in our case the ith star) is given by

$$\begin{eqnarray}&&{\dot{m}}_{\mathrm{acc},i}={f}_{\sup }{\dot{m}}_{\mathrm{BHL},i}\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&{f}_{\sup }=\min \left\{1,\max \left[{\left(\displaystyle \frac{{\alpha }_{\mathrm{SS}}}{0.01}\right)}^{0.62}\displaystyle \frac{{r}_{\mathrm{in},i}}{{R}_{\mathrm{BHL},i}},{f}_{\sup ,\min }\right]\right\}\end{eqnarray} \tag{ 32 }$$

is the suppression rate from the Bondi accretion rate, α_SS is the viscosity parameter in the standard thin α-disk model (Shakura & Sunyaev 1973), ${\dot{m}}_{\mathrm{BHL},i}=4\pi {G}^{2}{n}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}{m}_{i}^{2}/{({c}_{s}^{2}+{v}_{i}^{2})}^{3/2}$ is the Bondi–Hoyle–Lyttleton accretion rate, ${R}_{\mathrm{BHL},i}={{Gm}}_{i}/({c}_{s}^{2}+{v}_{i}^{2})$ is the Bondi–Hoyle–Lyttleton radius, ${r}_{\mathrm{in},i}$ is the inner radius, and ${f}_{\sup ,\min }$ is the minimum of the suppression rate. Inayoshi et al. (2018b) found that ${f}_{\sup ,\min }\,\sim \,{10}^{-2}\mbox{--}{10}^{-3}$. We set ${f}_{\sup ,\min }\,=\,0.003$. The inner radius is the inner boundary of the calculation introduced in Inayoshi et al. (2018b) due to the computational limit. The inner radius is considered to correspond to the stellar radius ${r}_{\mathrm{in},i}={R}_{i}$. We set ${\alpha }_{\mathrm{SS}}=0.01$ as a fiducial value. This value is motivated by the results for a weak vertical magnetic field by Bai & Stone (2013), which simulates the magnetorotational instability turbulence (see also King et al. 2007). We limit the maximum accretion rate to ${\dot{m}}_{\mathrm{acc},i}={\dot{m}}_{\mathrm{BHL},i}$, if ${\dot{m}}_{\mathrm{acc},i}\gt {\dot{m}}_{\mathrm{BHL},i}$ given by Equation (31), since in this case Equation (31) becomes invalid, which describes the reduction in the accretion rate relative to the Bondi–Hoyle–Lyttleton rate due to rotation. We do not consider the enhancement of the gas density due to the N-body accretion (Kaaz et al. 2019), since the upper limit on the density of gas outside of the Bondi–Hoyle–Lyttleton radius is given by n_c.

If the velocity of the $i\mathrm{th}$ surrounding star is sufficiently high, ${R}_{\mathrm{BHL},i}$ may become smaller than R_i. In this case, the gas accretion rate is determined by direct collision with the stellar surface, ${\dot{m}}_{\mathrm{coll},i}=\pi {R}_{i}^{2}{n}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}{v}_{i}$. We set the accretion rate to $\max [{\dot{m}}_{\mathrm{acc},i},{\dot{m}}_{\mathrm{coll},i}]$. Furthermore, when the total accretion rate ${{\rm{\Sigma }}}_{i}{\dot{m}}_{i}$ onto all stars exceeds the inflow rate from large scales ${\dot{M}}_{\mathrm{in}}$, we normalize the respective accretion rate of each star by the inflow rate by multiplying it by ${\dot{M}}_{\mathrm{in}}/{{\rm{\Sigma }}}_{i}{\dot{m}}_{i}$. When ${{\rm{\Sigma }}}_{i}{\dot{m}}_{i}\sim {\dot{M}}_{\mathrm{in}}$, the gas density should be depleted. However, for simplicity, we assume that the gas density distribution is unchanged; this is justified since whenever this condition is satisfied, the dynamical evolution of surrounding stars is hardly affected by the presence of gas because the gravitational potential is dominated by stars in later phases and star formation ceases.

In cases in which the Bondi mass ${M}_{\mathrm{Bondi},i}\,=\tfrac{4}{3}\pi {R}_{\mathrm{BHL},i}^{3}{n}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}$, i.e., the gas mass within the Bondi–Hoyle–Lyttleton radius (${R}_{\mathrm{BHL},i}$) of the $i\mathrm{th}$ star, is larger than the stellar mass m_i, the gas within the Bondi–Hoyle–Lyttleton radius can be unstable to fragmentation since the Jeans instability can be significant in weak shear regions (e.g., Elmegreen 1994; Kim & Ostriker 2001; Kim et al. 2002). If fragmentation is significant, it is not obvious what fraction of the gas can accrete onto the star. The fraction depends on cooling, turbulence (e.g., Clark et al. 2011a; Elmegreen 2011; Greif et al. 2011; Dopcke et al. 2013), and the efficiency of angular momentum transfer (e.g., Thompson et al. 2005). Following the prescription in Ryu et al. (2016), when the Bondi mass exceeds the stellar mass ${M}_{\mathrm{Bondi},i}\gt {m}_{i}$, we reduce the accretion rate ${\dot{m}}_{i}$ by a constant factor ${f}_{\mathrm{red}}$. We assume ${f}_{\mathrm{red}}={10}^{-3}$ as a fiducial value. Even when fragmentation is expected within ${r}_{\mathrm{BHL},i}$, we assume that surrounding stars form at r_c (Section 3.2).

Thus, in summary we calculate the gas accretion rate of stars as

$$\begin{eqnarray}{\dot{m}}_{i}=\left\{\begin{array}{cc}{\dot{m}}_{i,0} & \mathrm{if}\,{M}_{\mathrm{Bondi},i}\lt {m}_{i},\\ {f}_{\mathrm{red}}{\dot{m}}_{i,0} & \mathrm{if}\,{M}_{\mathrm{Bondi},i}\geqslant {m}_{i},\end{array}\right.\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}&&{\dot{m}}_{i,0}=\min \left({\dot{m}}_{i,1},\displaystyle \frac{{\dot{m}}_{i,1}{\dot{M}}_{\mathrm{in}}}{{\sum }_{i}{\dot{m}}_{i,1}}\right),\end{eqnarray} \tag{ 34 }$$

and

$$\begin{eqnarray}&&{\dot{m}}_{i,1}=\max ({\dot{m}}_{\mathrm{acc},i},{\dot{m}}_{\mathrm{coll},i}).\end{eqnarray} \tag{ 35 }$$

Photoionization and supernova feedback play key roles in ejecting gas from pregalactic halos (Kitayama et al. 2004; Whalen et al. 2004; Kitayama & Yoshida 2005). We did not take into account feedback from supernova explosions since our simulations are limited to the time until a first supernova explosion at 3 Myr. Kitayama et al. (2004) have shown that when the production rate of ionizing photons ${Q}_{\mathrm{ion}}={\sum }_{i}{Q}_{\mathrm{ion},i}$ is below the critical value ${Q}_{\mathrm{cri}}\sim {10}^{51}\,{{\rm{s}}}^{-1}{({M}_{\mathrm{halo}}/{10}^{7}{M}_{\odot })}^{8/5}{[(1+z)/15]}^{12/5}$, the gas density is not affected by photoionization feedback. On the other hand, when ${Q}_{\mathrm{ion}}$ exceeds ${Q}_{\mathrm{cri}}$, gas is blown away from the halo. We adopt this criterion as the gas ejection condition. ${Q}_{\mathrm{ion},i}$ strongly depends on whether the $i\mathrm{th}$ star is in the expanding phase or the contracting phase, with ${Q}_{\mathrm{ion},i}\sim {10}^{36}\,{{\rm{s}}}^{-1}{({m}_{i}/10{M}_{\odot })}^{2}$ and $\sim {10}^{48}\,{{\rm{s}}}^{-1}{({m}_{i}/10{M}_{\odot })}^{2}$ in these phases, respectively (Hosokawa et al. 2012). If ${Q}_{\mathrm{ion}}\gt {Q}_{\mathrm{cri}}$ is ever satisfied, all gas is assumed to be ejected from the system.

Although we add up the ionizing photons emitted by low-mass stars (≲1 ${M}_{\odot }$), these photons do not affect bulk gas dispersal due to their low numbers and because they are trapped within their parent stars' Bondi radii. Furthermore, the contraction timescale for low-mass stars exceeds the calculation time (3 Myr), and we expect that low-mass stars contribute negligibly to the total photon emission rate. Indeed, we find below that gas dispersal (when it occurs) is caused by the contraction of the central star in our models.

After the gas is ejected, gas accretion, star formation, and gas dynamical friction are all assumed to stop operating (${\dot{m}}_{i}={\dot{m}}_{\mathrm{SF}}={n}_{\mathrm{gas}}(r)=0$ for all i and r) during the rest of the simulation.⁶ The radial position of surrounding star i increases to ${r}_{i}+{\rm{\Delta }}r$ due to the decrease of the potential energy as

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}({r}_{i}+{\rm{\Delta }}r)-{{\rm{\Phi }}}_{\mathrm{gas}}({r}_{i}+{\rm{\Delta }}r)+{k}_{\mathrm{ej}}({r}_{i}+{\rm{\Delta }}r)\\ \qquad =\,{\rm{\Phi }}({r}_{i})-{{\rm{\Phi }}}_{\mathrm{gas}}({r}_{i})+{k}_{\mathrm{be}}({r}_{i}),\end{array}\end{eqnarray} \tag{ 36 }$$

where ${k}_{\mathrm{be}}(r)$ and ${k}_{\mathrm{ej}}(r)$ are the specific kinetic energies of an object at radius r with and without gas, respectively. As in Equation (17), we use the zero-eccentricity approximation when we calculate the change in the radial position of surrounding stars.

Although we assume that gas is ejected instantly, the ejection timescale is roughly given by the size of the gas cloud over the ejection speed. In our models, the gas distribution affects the dynamical evolution of surrounding stars, and most surrounding stars are distributed within 0.1 pc. The ejection timescale for gas within 0.1 pc is $\sim {10}^{4}\,\mathrm{yr}$ when the ejection speed is $\sim 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which is set to a rough value of the sound speed of ionized gas. Thus the ejection timescale is much smaller than our total calculation timescale of 3 Myr, which justifies the assumption of instant gas ejection.

We have performed several numerical calculations using the above semianalytical model. We first present the evolution of the central star in the fiducial model (labeled as Model 1 in Table 1). Figure 2 shows the evolution of several other quantities in this model.

Zoom In Zoom Out Reset image size

Table 1.  The Results of Our Simulations in the Fiducial Model (Model 1) and 21 Variants

Input						Output
Model	nc	$T/{10}^{4}$	fred	β	${m}_{0,\max }$	mfin	${m}_{\mathrm{acc},* }$	Mloss	NSF	Nacc	Ncoll	mej	tej	${\gamma }_{\mathrm{HC}}$
	$({\mathrm{cm}}^{-3})$	(K)			$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$				$({M}_{\odot })$	$(\mathrm{yr})$
1	1011	1	10−3	0	1	$6.7\times {10}^{5}$	$9.6\times {10}^{3}$	$3.5\times {10}^{2}$	$8.7\times {10}^{3}$	$8.7\times {10}^{3}$	$6.1\times {10}^{3}$	⋯	⋯	0.50
2	$3\times {10}^{10}$	1	10−3	0	1	3.7 × 103	3.7 × 103	1.6 × 103	1.2 × 104	5.1 × 103	5.7 × 104	2.4 × 102	3.2 × 104	0.19
3	1011	0.5	10−3	0	1	2.4 × 105	3.8 × 103	1.6 × 102	4.1 × 103	4.1 × 103	3.3 × 103	⋯	⋯	0.77
4	1011	0.3	10−3	0	1	1.1 × 105	1.8 × 103	89	2.2 × 103	2.1 × 103	2.0 × 103	⋯	⋯	0.67
5	1011	0.1	10−3	0	1	2.8 × 102	2.4 × 102	57	5.0 × 102	3.2 × 102	3.3 × 103	2.4 × 102	4.8 × 104	0.67
6	1011	1	0	0	1	6.7 × 105	6.7 × 105	4.3 × 102	6.7 × 105	6.7 × 105	7.6 × 103	⋯	⋯	0.50
7	1011	1	10−3	2.35	1	6.7 × 105	1.1 × 104	5.7 × 102	2.6 × 104	2.6 × 104	2.8 × 104	⋯	⋯	0.18
8	3 × 1010	0.5	10−3	0	1	1.7 × 103	1.7 × 103	6.2 × 102	4.8 × 103	2.4 × 103	2.4 × 104	2.4 × 102	3.5 × 104	0.23
9	3 × 1010	0.3	10−3	0	1	9.3 × 102	9.3 × 103	35	2.6 × 103	1.3 × 103	1.5 × 104	2.4 × 102	4.1 × 104	0.32
10	3 × 1010	1	10−3	0	0.1	3.7 × 102	3.5 × 102	1.2 × 102	8.7 × 102	4.8 × 102	5.6 × 103	2.4 × 102	7.5 × 104	0.29
11	3 × 1010	1	10−2	0	1	3.6 × 103	3.6 × 103	1.6 × 103	1.2 × 104	5.0 × 103	5.4 × 104	2.4 × 102	3.1 × 104	0.19
12	3 × 1010	1	10−1	0	1	6.7 × 105	7.1 × 103	1.4 × 103	1.5 × 104	9.2 × 103	5.1 × 104	⋯	⋯	0.19
13	3 × 1010	1	10−3	2.35	1	2.7 × 103	2.7 × 103	2.2 × 103	4.4 × 104	6.4 × 103	1.8 × 105	2.4 × 102	4.1 × 104	0.073
14	1010	1	10−3	0	1	4.1 × 103	4.1 × 103	2.7 × 103	2.3 × 104	5.1 × 103	7.6 × 104	2.4 × 102	5.9 × 104	0.076
15	1010	1	10−2	0	1	4.1 × 103	4.1 × 103	2.7 × 103	2.3 × 104	5.1 × 103	7.6 × 104	2.4 × 102	5.9 × 104	0.076
16	1010	1	10−1	0	1	4.1 × 103	4.0 × 103	2.6 × 103	2.2 × 104	5.0 × 103	7.4 × 104	2.4 × 102	5.7 × 104	0.076
17	1010	1	1	0	1	6.6 × 105	1.3 × 103	3.8 × 102	1.1 × 104	1.5 × 103	1.3 × 104	⋯	⋯	0.076
18	109	1	10−3	0	1	1.2 × 103	1.1 × 103	5.2 × 102	9.6 × 104	1.0 × 103	9.7 × 103	2.4 × 102	2.4 × 105	0.085
19	109	1	1	0	1	2.6 × 102	1.7 × 102	1.5 × 102	9.0 × 104	1.6 × 102	3.2 × 103	2.4 × 102	2.5 × 105	0.085
20	108	1	10−3	0	1	3.2 × 102	2.8 × 102	1.3 × 102	3.6 × 105	2.2 × 102	2.5 × 103	2.4 × 102	8.8 × 105	0.081
21	107	1	10−3	0	1	2.0 × 102	1.9 × 102	56	1.2 × 106	1.6 × 102	1.1 × 103	⋯	⋯	0.082
22	106	1	10−3	0	1	22	21	2.4	1.2 × 106	17	34	⋯	⋯	0.079
23	1011	1	10−3	0	10	6.7 × 105	8.4 × 103	3.1 × 102	1.2 × 103	1.2 × 103	7.7 × 102	⋯	⋯	3.1
24	1010	1	10−3	0	10	6.6 × 105	5.2 × 104	4.8 × 103	6.4 × 103	6.3 × 103	1.1 × 104	⋯	⋯	0.74
25	109	1	10−3	0	10	4.2 × 103	4.2 × 103	9.0 × 102	1.4 × 103	5.3 × 102	2.7 × 103	2.5 × 102	3.2 × 104	0.75
26	1011	1	10−3	0	100	6.6 × 105	8.0 × 103	2.0 × 102	1.6 × 102	1.4 × 102	94	⋯	⋯	8.8
27	1010	1	10−3	0	100	9.0 × 103	8.8 × 103	3.0 × 102	2.1 × 102	1.5 × 102	83	4.3 × 103	5.4 × 104	3.6
28	109	1	10−3	0	100	2.8 × 102	1.3 × 102	0	14	2	0	2.7 × 102	8.1 × 104	3.9

Note. The columns show several input and output parameters in each case, as follows: the model number, the core gas number density (n_c), the gas temperature (T), the reduction factor for the gas accretion rate when the Bondi mass exceeds the central mass (f_red), the power-law index of the stellar initial mass function (β), the maximum initial mass of stars (${m}_{0,\max }$), the mass of the central star at the end of the simulation at 3 Myr (m_fin), the total mass of surrounding stars accreted onto the central star (${m}_{\mathrm{acc},* }$), the total mass lost during collisions (M_loss), the number of newly formed surrounding stars (N_SF), the number of surrounding stars accreted onto the central star (N_acc), the number of collisions between surrounding stars (N_coll), the mass of the central star and the time at the ejection of gas from the system (m_ej and t_ej) for models in which such ejection occurs, and the maximum value for the ratio of the heating rate by gas dynamical friction to the cooling rate by dust grains (${\gamma }_{\mathrm{HC}}=\max ({{\rm{\Gamma }}}_{\mathrm{GDF}}/{{\rm{\Lambda }}}_{\mathrm{dust}}$)).

Download table as:  ASCIITypeset image

In the early stages, the gas accretion rate onto the central star ${\dot{m}}_{\mathrm{acc},i}$ is as low as $\sim 6\times {10}^{-5}{({m}_{i}/1{M}_{\odot })}^{1.5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (blue line in panel (c) of Figure 2), and almost all of the gas flowing in from large scales is converted into new surrounding stars, at the rate of ${\dot{m}}_{\mathrm{SF}}\sim {\dot{M}}_{\mathrm{in}}\simeq 0.22\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (orange line in panel (c)), and so a dense $\sim {10}^{3}\,{M}_{\odot }$ stellar cluster forms in $\sim {10}^{4}$ yr (orange line in panel (a)).

At $2.6\times {10}^{2}\,\mathrm{yr}$, the central star accretes the first surrounding star (red line in panel (c)). The accretion rate of surrounding stars subsequently gradually increases, due to the increasing radius of the central star (eventually to $\sim {10}^{3}$ au) as well as due to the increase in the number of surrounding stars (red line in panel (c), black line in panel (b), and orange line in panel (a)). At the same time, the surface KH time ${t}_{\mathrm{surf},\mathrm{KH}}$ for the central star decreases significantly (blue line in panel (b)) due to its increase in mass, which significantly raises the luminosity up to ${m}_{\mathrm{cent}}\sim 6\,{M}_{\odot }$ (Hosokawa et al. 2012).

Before the radius of the central star increases to 10 au, 77 collisions occur between surrounding stars (cyan line in panel (b)). All of these collisions occur between expanding surrounding stars; 4.7 ${M}_{\odot }$ is lost during stellar collisions; and five pairs of them merge as a result of collisions. The mass lost during collisions is added to the gas inflow rate ${\dot{M}}_{\mathrm{in}}$.

At $1.7\times {10}^{3}\,\mathrm{yr}$ the central star's mass is ${m}_{\mathrm{cent}}=19\,{M}_{\odot }$, and its growth rate exceeds the critical rate required to inhibit contraction (${\dot{m}}_{\mathrm{cri}}=0.01\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$, see Section 3.3; black and gray lines in panel (c)). The accretion rate of surrounding stars onto the central star subsequently increases further, due to the increasing radius of the central star, the increased total number of surrounding stars in the cluster, and the increased masses of the surrounding stars (black line in panel (b), orange and cyan lines in panel (a)). Thus the mass of the central star rapidly increases by stellar bombardment in this phase.

At t = 2.2 × 10³ yr, the Bondi mass of the central star exceeds its own mass, m_cent = 31 ${M}_{\odot }$ (blue and black lines in panel (a)), and the gas accretion rate is reduced by f_red = 10⁻³ due to the fragmentation of gas within the Bondi–Hoyle–Lyttleton radius (see Section 3.6.2).

Since the cooling rate due to emission by dust grains (Equation (15)) always exceeds the heating rate due to gas dynamical friction (blue and cyan lines in panel (d)), surrounding stars continue forming at the core radius (Section 3.2).

At $\gtrsim {10}^{5}\,\mathrm{yr}$, gas accretion begins to dominate the central star's growth rate. In this phase, most of the gas flowing in from large scales ${\dot{M}}_{\mathrm{in}}$ is accreted onto the central star. Star formation ceases, and a large number of surrounding stars are absorbed by the central star due to the radial growth of the central star (orange lines in panels (c) and (a) and black line in panel (b)). Stellar bombardment keeps the central star in the bloated state. Hence, the central star continues growing, and reaches $\gtrsim {10}^{5}\,{M}_{\odot }$ without any contraction.

During the evolution, $6.1\times {10}^{3}$ collisions occur (cyan line in panel (b)); all are between expanding surrounding stars; and 31 pairs of them merge as a result of collisions. In total $3.5\times {10}^{2}\,{M}_{\odot }$ is lost by surrounding stars during collisions (Section 3.5). Thus most collisions between surrounding stars result in a relatively small amount of mass loss in the fiducial model.

To clarify the importance of each mechanism for the growth of the central star, we repeat the above calculation in several variants of the fiducial model, in which parameter settings remain unchanged but some mechanism is turned off and does not operate.

To check whether the stellar bombardment plays an important role, we first run the model with the fiducial settings but no migrating motion for surrounding stars. In this model, the final mass of the central star is found to be 32 ${M}_{\odot }$. Thus via gas accretion alone, we find that the central star contracts and cannot grow into an SMS.

We next investigate the importance of dynamical friction. With stellar dynamical friction turned off, the final mass of the central star is $750\,{M}_{\odot }$. On the other hand, in the model without gas dynamical friction or without gas accretion drag, respectively, the final masses of the central stars are $6.7\times {10}^{5}\,{M}_{\odot }$ and $6.6\times {10}^{5}$ M_☉. We conclude that the migration of surrounding stars is dominated by stellar dynamical friction rather than gas dynamical friction and gas accretion drag. This is essentially because the density of stars dominates the density of gas. For example, in Model 1, the gas mass within the core radius of ${r}_{{\rm{c}}}=880\,\mathrm{au}$ at $t={10}^{4}\,\mathrm{yr}$ is $9.9\times {10}^{2}\,{M}_{\odot }$, while that of stars is $2.1\times {10}^{3}\,{M}_{\odot }$. The stellar density increases closer to the SMBH, while the gas density is set to be constant within the core (e.g., upper and middle panels of Figure 4).

We further simulate a model with the fiducial settings, but without allowing the stars to expand, and instead always setting their radii to the value in Equation (2). In this model, the rate of stellar accretion onto the central star does not increase beyond $\sim 0.005\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, and the final mass is found to be $1.7\times {10}^{3}\,{M}_{\odot }$. Thus the bloating of stars is required to facilitate the growth of the central star due to stellar accretion.

We next present a case in which the central star contracts before it collapses to a BH. Figure 3 shows the results in Model 2, which differs from Model 1 only by a modified value of the gas density (reduced by a factor of 3 to ${n}_{{\rm{c}}}=3\times {10}^{10}\,{\mathrm{cm}}^{-3}$). Initially the radial position at which surrounding stars form is about a factor of 1.5 larger than in Model 1 (orange lines in panels (b) in Figures 3 versus 2). Stellar dynamical friction becomes less efficient due to the lower stellar density, and the average accretion rate of surrounding stars onto the central star becomes lower than in Model 1 (red line in panel (c) of Figure 3). At ${t}_{\mathrm{KH},\mathrm{ZAMS}}(\sim 3\times {10}^{4}$ yr) for the central star, the central star contracts, and then the production rate of ionizing photons exceeds the critical value for gas ejection from the system (black and gray lines in panel (d)). After the ejection, gas accretion and star formation cease (blue and orange lines in panel (c)), and the rate of accretion of surrounding stars decreases (red line in panel (c)). The radius of a star is predicted to contract in $\sim {10}^{2}\mbox{--}{10}^{3}\,\mathrm{yr}$ (e.g., Sakurai et al. 2015), which justifies the assumption of abrupt contraction in our calculations. The central star continues to grow by accreting surrounding stars, but at a more moderate rate, reaching the mass of $3.7\times {10}^{3}\,{M}_{\odot }$ at 3 Myr. Therefore a massive BH may still be produced in Model 2, but the mass of this massive BH is $\approx 100$ times below that of the BH remnant in Model 1.

Zoom In Zoom Out Reset image size

The top panel of Figure 4 shows the stellar and gas density profiles for Model 2. The power-law slope of the stellar density is almost unchanged during the evolution (black, orange, and red curves in the top panel). Coincidentally, such self-similar evolution is also expected for the core collapse of a self-gravitating system driven by two-body relaxation (Binney & Tremaine 2008). The evolution of the stellar density in our model is driven by the combination of gas and stellar dynamical friction, Bondi accretion, and star formation. On the other hand, the outer cutoff of the stellar density distribution slowly evolves from 0.1 to 3 Myr (orange and red curves in the top panel of Figure 4). This is because the timescale for stellar dynamical friction (${t}_{\mathrm{SDF}}\equiv {v}_{\mathrm{Kep}}(r)/{a}_{\mathrm{SDF},i}$ with ${m}_{i}={\overline{m}}_{l}$) at this position (orange curve in the bottom panel) exceeds the calculation timescale of 3 Myr. At 3 Myr, the density profile contains zero surrounding stars within 10 au, and ≈50 surrounding stars within 100 au.

Zoom In Zoom Out Reset image size

Figure 5 shows the evolution of one of the surrounding stars in Model 2. This star is born at $1.2\times {10}^{3}\,\mathrm{yr}$ with mass ${m}_{i}=0.50\,{M}_{\odot }$ at ${r}_{i}=660\,\mathrm{au}$. In the early phase, due to gas dynamical friction, accretion drag and stellar dynamical friction, the star migrates inwards very slightly (orange line in panel (a) and brown, cyan, and orange lines in panel (c)).

Zoom In Zoom Out Reset image size

At $\gtrsim 2.5\times {10}^{3}\,\mathrm{yr}$, the average mass in the cell hosting this star becomes more massive than the mass of the star due to the formation of new stars within the cell (orange and black lines in panel (b)). The star therefore begins to migrate outward due to mass segregation.

At $3.2\times {10}^{4}\,\mathrm{yr}$, gas is ejected from the system due to photoionization feedback (as was shown by the gray and black lines in panel (d) of Figure 3), and so the binding energy of this star decreases abruptly (black line in panel (c) in Figure 5). Gas dynamical friction and gas accretion stop operating due to the lack of gas around the star. Since star formation also stops operating and massive surrounding stars migrate inward, the average masses in cells at ∼100–1000 au begin to decrease. At $\sim 2\times {10}^{5}\,\mathrm{yr}$, this star also begins to migrate inward. In the inner regions, the star and other surrounding stars lose some fraction of their mass due to frequent collisions (black and orange lines in panel (b)). Since the direction of migration due to stellar dynamical friction is influenced by the mass of the star compared to that of surrounding stars, the star wanders around ∼500 au. When the central star collapses into a massive BH at $3\times {10}^{6}\,\mathrm{yr}$, this star orbits at ${r}_{i}=380\,\mathrm{au}$, and the mass is ${m}_{i}=0.32\,{M}_{\odot }$. Hence surrounding stars are redistributed mainly by mass segregation driven by stellar dynamical friction, and only more massive stars can migrate toward the central star. In later phases, frequent collisions reduce the masses of surrounding stars, and they prevent accretion of surrounding stars onto the central star.

The dependence of the results on the input parameters of our model is illustrated through a range of model variants listed in Table 1. The final mass of the central star (m_fin) is most strongly influenced by whether the central star contracts or not, which in turn depends on the parameters we investigated. This is illustrated by the masses shown in Figure 6.

Zoom In Zoom Out Reset image size

We find that for efficient growth via stellar bombardment, the formation of a high-density star cluster is required in order to enhance the inward acceleration by stellar dynamical friction. In cases with high core gas density n_c, the core radius r_c is small, and since surrounding stars form at the core radius, the growth rate of the density of stars in early phases is high (Equation (12)). The growth rate of the stellar density is also high for the high T cases, since the star formation rate during the early stages is mostly given by the gas inflow rate (${\dot{M}}_{\mathrm{in}}\propto {T}^{3/2}$, Equation (4)). In high stellar density environments, the migration time due to stellar dynamical friction is short and the rate of stellar bombardment is high. When the growth rate by stellar bombardment exceeds the critical rate, the central star continues growing without ejecting gas, as seen in the evolution for Model 1 in Figure 2. From Table 1, for ${n}_{{\rm{c}}}={10}^{11}\,{\mathrm{cm}}^{-3}$ with $T\geqslant 3\times {10}^{3}$ (Models 1, 3, and 4) the stellar accretion rate onto the central star exceeds the critical rate for contraction before ${t}_{\mathrm{KH},\mathrm{ZAMS}}$ for the central star, allowing it to grow to $\gtrsim {10}^{5}\,{M}_{\odot }$.

Whether the central star contracts is also influenced by the value of f_red (see Sections 3.3 and 3.6.2). This is the factor by which the gas accretion rate is assumed to be reduced by both fragmentation and by the removal of gas that is captured by the fragmented clumps, when the Bondi mass becomes larger than the star's own mass. A high f_red value increases the gas accretion rate for ${M}_{\mathrm{Bondi}}\gt {m}_{\mathrm{cent}}$, which is satisfied for ${m}_{\mathrm{cent}}\gtrsim 30\,{M}_{\odot }$ in Models 1 and 2 (blue and black lines in panel (a) of Figures 2 and 3). So if the central star can grow to ${m}_{\mathrm{cent}}\gtrsim 30\,{M}_{\odot }$ within ${t}_{\mathrm{KH},\mathrm{ZAMS}}$, the high value for f_red can aid to enhance the growth rate of the central star. From Table 1, we see that for ${n}_{{\rm{c}}}=3\times {10}^{10}\,{\mathrm{cm}}^{-3}$ with ${f}_{\mathrm{red}}\geqslant 0.1$ (Model 12) or ${n}_{{\rm{c}}}={10}^{10}\,{\mathrm{cm}}^{-3}$ with ${f}_{\mathrm{red}}\sim 1$ (Model 17), the central star keeps expanding until 3 Myr when the SMS collapses to a massive BH or when any of the surrounding stars explode as a supernova and blow away all of the gas from the vicinity. In Models 12 and 17, the central star grows mainly via stellar accretion until ${m}_{\mathrm{cent}}\sim 100\,{M}_{\odot }$, and then the gas accretion rate onto the central star exceeds the critical rate ${\dot{m}}_{\mathrm{cri}}$. Thus even for high f_red, stellar accretion is important to enhance the gas accretion rate. Unfortunately, the relevant value for f_red is highly uncertain. To assess it, we need to consider fragmentation of gas inside the Bondi radius, and the evolution of any accretion disk around the central star. These issues are beyond the scope of the present paper and will be investigated elsewhere in the future.

In those cases in which the central star contracts and gas is ejected before 3 Myr (${t}_{\mathrm{ej}}\,\lt$ 3 Myr), gas accretion contributes very little to the final mass (see the values of m_fin and ${m}_{\mathrm{acc},* }$ in Table 1). In these cases, since the growth rate should correlate with the efficiency of stellar dynamical friction, which depends on the stellar density, the final mass of the central star should correlate with the stellar density. We further assume that the stellar density is proportional to the star formation rate over the core radius cubed (${\rho }_{* }\propto \sim {\dot{m}}_{\mathrm{SF}}/{r}_{{\rm{c}}}^{3}$). Due to the scaling relations ${\dot{m}}_{\mathrm{SF}}\propto \sim {T}^{3/2}$ and ${r}_{{\rm{c}}}\propto {n}_{{\rm{c}}}^{-1/3}$ (Equation (12)), we can expect ${\rho }_{* }\propto \sim {T}^{3/2}{n}_{{\rm{c}}}$. Figure 6 shows the relation between the final mass for the central star and the product ${T}^{3/2}{n}_{{\rm{c}}}$. We can indeed see the rough correlation between the final mass and ${T}^{3/2}{n}_{{\rm{c}}}$ as expected in the cases in which the central star contracts (dashed line and circles in Figure 6). Also, the expected relation ${\rho }_{* }\propto \sim {T}^{3/2}{n}_{{\rm{c}}}$ is roughly confirmed in Figure 7, which shows the maximum stellar density within the core radius as a function of ${T}^{3/2}{n}_{{\rm{c}}}$. The offsets for high ${m}_{\max ,0}$ or low T cases (large empty or red circles/squares in Figure 7) are presumably due to the differences in the efficiency of migration, which changes the accretion rate of surrounding stars onto the central star and so affects the stellar density within the core radius.

Zoom In Zoom Out Reset image size

As discussed above, if T is low, ${\rho }_{* }$ remains low, and therefore stellar dynamical friction is inefficient. Additionally, since the mass of the stellar cluster in the core is approximately limited by ${\dot{M}}_{\mathrm{in}}t$, if T is low, then ${\dot{M}}_{\mathrm{in}}$ is low (Equation (4)) and the cluster mass grows slowly. If the cluster mass remains low, the number of surrounding stars that bombard the central star is reduced. This is plausibly the reason why the final mass of the central star at some fixed values of the combination ${T}^{3/2}{n}_{{\rm{c}}}$ in low T models is lower than those for high T models (Figure 6). Also in those cases when the central star keeps expanding until it collapses into a massive BH, the growth rate is determined primarily by the gas accretion rate in the final phase, the final mass depends on the inflow rate and accordingly the gas temperature (square plots in Figure 6). This dependence explains why the final masses in Models 1, 6, and 7, which have the same temperature, are the same. Note that, in the cases in which the central star keeps expanding for 3 Myr, almost all of the gas that fell in from large scales is converted to the central star (${m}_{\mathrm{fin}}\sim 3\,\mathrm{Myr}\times {\dot{M}}_{\mathrm{in}}$).

On the other hand, the power-law slope β of the IMF has only a small effect on the final mass (empty circle and square in Figure 6). This is because β has almost no effect on the density and mass of the stellar cluster, which are the critical factors for the efficiency of migration by stellar dynamical friction.

For high ${m}_{0,\max }$, the final mass is higher than that for low ${m}_{0,\max }$ except in Model 28 (large empty symbols in Figure 6). When the masses of surrounding stars are high, stellar dynamical friction operates prominently, which facilitates the growth of the central star. For high-mass stars, gas dynamical friction also operates efficiently, which enhances the heating of gas. In our models, when ${{\rm{\Gamma }}}_{\mathrm{GDF}}\gt {{\rm{\Lambda }}}_{\mathrm{dust}}$, star formation is assumed to stop operating due to the gas heating. In Model 28, star formation becomes inefficient due to this effect, resulting in a low final mass of the central star. However, our models cannot predict the evolution in this case, since the gas distribution and accretion processes will be affected by the increased gas pressure. Additional studies using three-dimensional hydrodynamical simulations are required to estimate the evolution in these cases, i.e., when ${{\rm{\Gamma }}}_{\mathrm{GDF}}\gtrsim {{\rm{\Lambda }}}_{\mathrm{dust}}$.

To calculate the acceleration rate due to stellar dynamical friction, the velocity dispersion of surrounding stars is assumed to be isotropic. Such a thermalized distribution for the surrounding stars is realized during the evolution due to the nonresonant and resonant relaxation processes (Kocsis & Tremaine 2011). We note that, technically, this isotropic distribution is inconsistent with the assumption that the surrounding stars follow circular orbits (in the calculation of the migration rate in Equation (17)). However, since migration rates for nonzero- and zero-eccentricity surrounding stars are the same when the binding energy is dissipated by the same amount, this inconsistency should have a negligible impact on our results (i.e., on the evolution of the central star).

There is an inconsistency between the star formation prescription, assuming surrounding stars form in a rotating gas disk, and Equation (21), which assumes that surrounding stars are isotropically distributed. We expect that this does not significantly affect our conclusions. First, the gas disk thickness h/r roughly evolves from ∼0.4 to ∼0.08 from 10³ to 10⁵ yr in Model 1, and never reaches very small values. Second, we expect that an isotropic distribution is established by relaxation processes (e.g., Kocsis & Tremaine 2011). Finally, even if relaxation processes are inefficient, stellar dynamical friction would operate more strongly in a disk configuration, due to the higher stellar density and the low relative velocity between surrounding stars, which would facilitate stellar accretion. Thus the isotropic distribution of surrounding stars is a conservative choice for the growth rate of the central star.

Although the disk around the central star is thick, we used the approximation $h\sim {c}_{s}/{\rm{\Omega }}$ for the scale height of the disk. At $h/r=0.4$, the Toomre Q parameter is overestimated by ∼10% due to the approximation. Since Q depends linearly on n_c, this assumption may affect the dependence of the results on n_c by the same factor of 10%, which is well within other uncertainties of our simplified model.

In our models, we allow a high star formation efficiency (SFE), defined as the ratio of the total mass in newly formed stars to the initial gas mass. For example, the SFE within the core radius is ∼0.7 at $t={10}^{4}\,\mathrm{yr}$ in our fiducial Model 1 (see below), and it increases with time. Observationally, some massive molecular clouds are found to have an SFE of $\gt 0.5$ (Turner et al. 2015), though the SFEs of most molecular clouds in the Milky Way are ∼0.002–0.3 (Murray 2011). On the other hand, theoretically the SFE is determined by radiation pressure from ionizing ultraviolet (UV) photons, nonionizing UV photons, and infrared (IR) photons (e.g., Kim et al. 2018). Radiation pressure from nonionizing UV photons does not halt gas collapse when the gas surface density exceeds a critical value (Raskutti et al. 2016; Thompson & Krumholz 2016), and likewise IR photons do not halt collapse unless the IR opacity is very high (Skinner & Ostriker 2015). In our models, the gas surface density within the core radius (Equation (3)) is much higher than the critical value (Raskutti et al. 2016; Thompson & Krumholz 2016), and the IR opacity is extremely low because the gas is metal-poor.

We also estimate whether ionizing UV photons are confined within the Bondi–Hoyle–Lyttleton radius of each star (in which case they do not halt gas collapse; Section 3.7). According to numerical simulations (Skinner & Ostriker 2015; Raskutti et al. 2016; Thompson & Krumholz 2016; Kim et al. 2018), when these feedback effects are inefficient, the SFE is close to unity (but not exactly 1 in their simulations due to the initial turbulent motion). Thus we considered the SFE of ∼1 to be justified in our case. On the other hand, although the SFE within the core radius becomes close to 1, the SFE within the rest of the halo is still low in our models since the baryon mass within the halo is $\sim 2\times {10}^{6}\,{M}_{\odot }$, and the mass of the stellar cluster is at most $\sim {10}^{4}\,{M}_{\odot }$ (see orange line in panel (a) of Figure 2 below), so the SFE might not be so extreme compared to the SFE observed in molecular clouds (0.002–0.5). Also as mentioned earlier, the rate of star formation is sensitive to the gas temperature in our models. Compared to the fiducial case of T = 10⁴ K, the star formation rates for T = 5 × 10³, 3 × 10³, and 10³ K are lower by factor of 2.8, 6.1, and 32, respectively. These lower-T models may be considered as proxies for lower SFE cases.

In this study, surrounding stars are allowed to migrate in- or outward, but are assumed to remain on circular orbits. If angular momentum exchange dominates the accretion of surrounding stars, stellar accretion becomes more efficient than in our model, since the binding energy of a surrounding star required to accrete onto the central star decreases by a factor of 1/(1 − e_i) where e_i is the eccentricity of the ith surrounding star. To investigate the impact of nonzero eccentricity, we examine a case in which the eccentricity distribution for surrounding stars is assumed to be thermalized (e.g., due to two-body relaxation), and has a distribution function of $f({e}_{i})=2{e}_{i}{{de}}_{i}$ (e.g., Jeans 1919; Heggie 1975). In this case, the central star captures surrounding stars from the larger distance ${r}_{i}={R}_{\mathrm{cent}}/(1-{e}_{i})$ (this is the only difference from the models above). Simulating this prescription with the parameter set of Model 2, we find the final mass of the central star to be ${m}_{\mathrm{fin}}=4.6\times {10}^{3}\,{M}_{\odot }$, which is almost unchanged from the final mass in Model 2 (${m}_{\mathrm{fin}}=3.7\times {10}^{3}\,{M}_{\odot }$). However, this neglects other possible effects. For example, a surrounding star with a high eccentricity interacts with stars and gas orbiting over a wider ranges of r, and mass loss should increase when a surrounding star with extremely high eccentricity is captured.

If most stellar accretion onto the central star is highly eccentric, and the mass lost at stellar accretion is typically a large fraction of the mass of the accreted surrounding star, the results of our models may be largely influenced. We intend to explore these issues in a follow-up study, based on direct N-body and hydrodynamical simulations. Here we only briefly consider the possible fate of the lost gas. If the launch velocity of this gas is similar in magnitude to the collision velocity between the stars, then the gas is kicked out to at most the apocenter of the colliding star's orbit before the collision. On the other hand, due to the low specific angular momentum of the ejected gas, it would be circularized (presumably by shocks it encounters) near the central star, similar to the expectation in the context of tidal disruption of stars (Hayasaki et al. 2013, 2016; Bonnerot & Lu 2019). In the vicinity of the central star, the viscous timescale is very short. Thus the gas ejected in high-eccentricity collisions may end up promptly accreted onto the central star, leaving our results largely unchanged.

Finally, let us consider the evolution of the stellar cluster after the central star collapses to a massive BH. Since collisions, relaxation, and evaporation are important mechanisms for cluster evolution, we show the collision (black curve in the bottom panel of Figure 4), stellar dynamical friction (orange curve), and evaporation (red curve) timescales for the stellar cluster at 3 Myr in Model 2. For the collision timescale (Equation (26)), the collision radius is assumed to be twice the radius of stars with the average mass, and stars are assumed to be in the contracted phase (Equation (2)). We adopt the evaporation timescale to be ${t}_{\mathrm{evap},l}={f}_{\mathrm{evap}}{t}_{\mathrm{relax},l}$ (Binney & Tremaine 2008), where the factor f_evap is ∼300 for clusters with a single stellar mass, and without a massive black hole and gas (Spitzer 1987), ${t}_{\mathrm{relax},l}=0.34{\sigma }_{* ,l}^{3}/({G}^{2}{\overline{m}}_{l}{\rho }_{* ,l}\mathrm{ln}{\rm{\Lambda }})$ is the relaxation timescale (Binney & Tremaine 2008), and we set the Coulomb logarithm to be 10. Although we set ${f}_{\mathrm{evap}}=300$, this value may be significantly increased for the cluster with a central massive BH which may help to retain objects from dynamical ejections both by increasing the cluster's escape velocity and by inhibiting binary formation.⁷ Figure 4 shows that the collision and evaporation timescales in the outer regions of the cluster are longer (at ≳2000 au where ${\rho }_{* }\sim {10}^{7\mbox{--}8}\,{M}_{\odot }\,{\mathrm{pc}}^{-3}$) and comparable to the Hubble time of $\sim 10\,\mathrm{Gyr}$, respectively. Thus these clusters could possibly survive to low-z epochs.

If such high-density clusters sink to the centers of massive local galaxies, the relics of such high-density clusters formed at high z may be observationally confused with stellar systems formed at lower redshift, such as infalling dense clusters and in situ formed stars, if those produce similarly high stellar density environments. The stellar density of nuclear star clusters may also be reduced by a supermassive black hole binary following galaxy collisions (Merritt 2006). On the other hand, if such clusters remain isolated, their relics may in principle be clearly identified in the local universe. Such clusters contain low-mass and extremely low-metallicity stars, and an intermediate-mass BH with the mass of $\sim {10}^{3}\,{M}_{\odot }$. Stellar densities within ∼2000 au of galactic nuclei have not been resolved to date (e.g., Nguyen et al. 2018). Extrapolating the observed density profile from diffuse light in the center of the Milky Way, the stellar mass within ∼2000 au from Sgr A* is estimated to be ∼600–800 ${M}_{\odot }$ (e.g., Schödel et al. 2018), which is about a factor ∼3 smaller than that for high-density clusters formed at high z (middle panel of Figure 4). If such high-density nuclear star clusters are identified with low-mass stars in the future, they might represent the fossils of high-z clusters.

In the stellar cluster in Model 2, the accretion of stars will continue after the BH formation. This may contribute to the rate of high-z tidal disruption events (Kashiyama & Inayoshi 2016) or to gravitational wave events observed by the Laser Interferometer Space Antenna (LISA; e.g., Amaro-Seoane et al. 2007; Hartwig et al. 2018).

After this work was submitted and posted on arXiv, Chon & Omukai (2020) presented results for three-dimensional hydrodynamical simulations focusing on similar scenarios. In this section, we briefly discuss the consistency between our predictions and their simulation results.

Chon & Omukai (2020) chose a relatively massive halo formed from cosmological initial conditions, in which the gas inflow rate is very high ($\gtrsim 1\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$). They find that the power-law index of the initial mass function is ∼−1, the maximum mass of surrounding stars is $\sim 100\,{M}_{\odot }$, and stars form at densities $\gtrsim {10}^{11}\,{\mathrm{cm}}^{-3}$. Referring to these findings, we perform models with ${n}_{{\rm{c}}}={10}^{11}\,{\mathrm{cm}}^{-3}$, ${m}_{0,\max }=100\,{M}_{\odot }$, β = −1, and $T=3\times {10}^{4}\,{\rm{K}}$. Since they use a barotropic equation of state, we allow the star formation even when ${{\rm{\Gamma }}}_{\mathrm{GDF}}\gt {{\rm{\Lambda }}}_{\mathrm{dust}}$. Since f_red is highly uncertain, we varied f_red between 10⁻³ and 1.

The evolution of ${f}_{\mathrm{red}}=0.1$ is shown in Figure 8. The central star grows to $\sim {10}^{4}\,{M}_{\odot }$ by 10 kyr (black line in the upper panel) mainly due to mergers in the early phases (red line in the lower panel) and gas accretion in the later phases (blue line in the lower panel).

Zoom In Zoom Out Reset image size

The total stellar mass stops increasing at ∼5 kyr for ${f}_{\mathrm{red}}=0.1$ (orange line in the upper panel) while it keeps increasing for ${f}_{\mathrm{red}}\leqslant 0.01$. Similar trends are seen in panel (d) of Figure 4 in Chon & Omukai (2020); namely the number of stars keeps increasing for $Z={10}^{-4}$ and $Z={10}^{-5}$, while it stops increasing at ∼5 kyr for other models. We conclude from our models that quenching of star formation in their simulations is related to whether the accretion rate onto the central star comes close to the gas inflow rate.

For ${f}_{\mathrm{red}}=0.001$, 0.01, 0.1, and 1, respectively, the mass of the central star at 10 kyr is 5.5 × 10³, 5.8 × 10³, 1.0 × 10⁴, and 1.1 × 10⁴ ${M}_{\odot }$, the number of surviving stars is 896, 743, 426, and 203, and the contribution of mergers to the total final mass of the central star is 98%, 87%, 33%, and 6.8%. Chon & Omukai (2020) find that the contribution of mergers is ∼30%–70%, the central mass is ∼10⁴ ${M}_{\odot }$, and the number of surviving stars is ∼500–4000 at 10 kyr. Thus our models with f_red ∼ 0.01–0.1 can reproduce their results remarkably well. The larger number of surrounding stars in Chon & Omukai (2020) presumably reflects the lower minimum mass of surrounding stars (∼0.01 ${M}_{\odot }$) for Z ∼ 10⁻³–10⁻⁵ compared to the value adopted in our models (0.08 ${M}_{\odot }$).

When we continue the above models beyond the 10 kyr at which Chon & Omukai (2020) stopped their simulation, we find that the final mass (at t = 3 Myr) of the central star for f_red = 10⁻³–1 is as high as 3.5 × 10⁶ due to the high inflow rate (${\dot{M}}_{\mathrm{in}}\sim 1.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$).