Making a Supermassive Star by Stellar Bombardment

Approximately 200 supermassive black holes (SMBHs) have been discovered within the first ∼gigayear after the Big Bang. One pathway for the formation of SMBHs is through the collapse of supermassive stars (SMSs). A possible obstacle to this scenario is that the collapsing gas fragments and forms a cluster of main-sequence stars. Here, we raise the possibility that stellar collisions may be sufficiently frequent and energetic to inhibit the contraction of the massive protostar, avoiding strong UV radiation driven outflows, and allowing it to continue growing into an SMS. We investigate this scenario with semianalytic models incorporating star formation; gas accretion; dynamical friction from stars and gas; stellar collisions; and gas ejection. We find that when the collapsing gas fragments at a density of ≲3 × 1010 cm−3, the central protostar contracts due to infrequent stellar mergers, and in turn photoevaporates the remaining collapsing gas, resulting in the formation of a ≲104 M⊙ object. On the other hand, when the collapsing gas fragments at higher densities (expected for a metal-poor cloud with Z ≲ 10−5 Z⊙ with suppressed H2 abundance) the central protostar avoids contraction and keeps growing via frequent stellar mergers, reaching masses as high as ∼105–106 M⊙. We conclude that frequent stellar mergers represent a possible pathway to form massive BHs in the early universe.

Growth from stellar-mass BH remnants of PopulationIII stars (e.g., Madau & Rees 2001;Abel et al. 2002;Heger & Woosley 2002;Tan & McKee 2004;McKee & Tan 2008;Yoshida et al. 2008;Clark et al. 2011b;Greif et al. 2011;Hirano et al. 2014;Susa et al. 2014;Stacy et al. 2016) to SMBHs is difficult because the gas accretion rate is suppressed by radiative and kinetic feedback processes (Whalen et al. 2004;Alvarez et al. 2009;Milosavljevic et al. 2009;Tanaka & Haiman 2009;Tanaka et al. 2012;Regan et al. 2019) and growth by mergers is made inefficient by large recoil induced by gravitational wave emission during mergers, which unbinds the merger remnant BHs from the shallow potential wells of their early hosts (Haiman 2004). These difficulties have motivated several alternative pathways.
One pathway is the direct collapse of supermassive stars (SMSs); (e.g., Omukai 2001;Oh & Haiman 2002;Bromm & Loeb 2003;Begelman et al. 2006;Spaans & Silk 2006;Shang et al. 2010;Agarwal et al. 2012;Hosokawa et al. 2012Hosokawa et al. , 2016Latif et al. 2013;Ferrara et al. 2014;Inayoshi et al. 2014;Sugimura et al. 2014;Tanaka & Li 2014;Becerra et al. 2015;Chon et al. 2016;Umeda et al. 2016;Haemmerlé et al. 2018). If a gas cloud in a massive halo with virial temperature T8000 K has no metals or H 2 molecules, the gas cloud can collapse without fragmentation and grow to become an SMS (Oh & Haiman 2002;Volonteri & Rees 2005). However, the background UV radiation flux required to prevent H 2 molecule formation is as high as a few times 10 4 in units of J 21 (see, e.g., Wolcott-Green & Haiman 2019, and references therein) because of the high density reached via atomic cooling (Omukai 2001;Oh & Haiman 2002) and self-shielding of H 2 for realistic UV spectra produced by PopulationII stars (Wolcott-Green et al. 2011Sugimura et al. 2014;Agarwal & Khochfar 2015). The condition of such a strong background radiation is satisfied only in rare cases, in collapsing halos that have bright nearby neighbors (Dijkstra et al. 2008). While this is a rare special configuration, it appears feasible for a sufficient number of such pairs of halos to form nearly simultaneously (Visbal et al. 2014), while avoiding metal pollution (Dijkstra et al. 2014), tidal disruption (Chon et al. 2016), and photoevaporation (Regan et al. 2017). For gas in halos located in regions of unusually high baryonic streaming motions , and/or in halos with unusually rapid merger histories experiencing compressional heating (Yoshida et al. 2003;Fernandez et al. 2014;Inayoshi et al. 2018a), the UV flux required to avoid H 2 cooling can be significantly reduced (Wise et al. 2019).
A second possible pathway is hyper-Eddington accretion onto a stellar-mass BH (Begelman 1979;Volonteri & Rees 2005;Pacucci et al. 2015Pacucci et al. , 2017Inayoshi et al. 2016;Sakurai et al. 2016;Sugimura et al. 2018;Takeo et al. 2018;Toyouchi et al. 2019). Here the problem is that inefficient angular momentum transfer is estimated to reduce the accretion rate (Inayoshi et al. 2018b;Sugimura et al. 2018, but see Alexander & Natarajan 2014 for a possible solution if the seed BH is surrounded by a massive and dense star cluster). Then the accretion of an isothermal rotating disk (Oh & Haiman 2002) may not be rapid enough to increase the mass of a BH by several orders of magnitude (Sugimura et al. 2018). Also kinetic feedback may limit the growth rate of BHs .
In this study, we focus on environments similar to the direct collapse scenario. A significant caveat of this scenario is that the collapsing gas may fragment efficiently, resulting in the formation of a cluster of stars, preventing the gas from fueling the formation of a central SMS. This would then lead to a third pathway, which is expected to produce BH remnants with masses ≈10 3-4 M e . On the other hand, we propose here that if stars themselves continue to be accreted efficiently, a more massive SMS may form despite the fragmentation of the parent cloud. In order for this to occur, incoming stars must collide with the central star in sufficiently rapid succession such that the central star never has time to cool and contract and settle on the main sequence. This scenario is similar to the runaway mergers above, but differs in detail. Stars are brought to the central region of the halo by both gas and stellar dynamical friction. The central SMS bloats up to astronomical unit size, facilitating the continued accretion of other stars. To distinguish this from the usual "runaway merger" case, we refer to this variant as "stellar bombardment." To investigate the feasibility of such a scenario, we have performed numerical modeling, incorporating star formation, dynamical friction by gas and stars, gas accretion, stellar collisions, and gas ejection.
After the submission of our paper, Chon & Omukai (2020) presented results investigating similar scenarios (collapsing gas with a small amount of metal pollution without H 2 molecules) using three-dimensional hydrodynamical simulations. They find that for metallicities up to 10 −4 Z e , a central star can keep growing to ∼10 4  M over ∼10 4 yr with high growth rate due to gas accretion and stellar accretion. Their results confirm the expectations in this paper.

Physical Picture
The SMBHs of ∼10 9  M at z∼6 are rare objects with an abundance ∼1 Gpc −3 , and thus rare conditions may be required to explain their formation (e.g., Buchner et al. 2019). The situation we consider is similar to the usual direct collapse scenario (e.g., Bromm & Loeb 2003;Shang et al. 2010). In this scenario, H 2 molecules are disrupted in the collapsing cloud by strong background radiation from nearby galaxies, the host halo is massive, and the collapsing cloud is not polluted by metals. These conditions keep the cloud at a high temperature, and so enable the cloud to collapse into an SMS without fragmentation. Recent studies have suggested that it may be difficult to satisfy these conditions (Latif et al. 2015), particularly because a large H 2 -dissociating flux may be required for an extended period, prior to reaching the "atomic cooling" threshold (Regan et al. 2017). This may be alleviated only in rare overdense regions via dynamical heating accompanying unusually rapid merger histories (Wise et al. 2019).
Here, instead, we relax the assumption of (the lack of) metal pollution. We consider a massive host halo, a moderate amount of metal pollution of ∼10 −5 Z e , and no H 2 molecules in a collapsing cloud. In such environments, fragmentation only occurs in high-density regions of ∼10 11 cm −3 due to weak cooling by a small amount of dust grains Latif et al. 2016). After an ultrahigh-density star cluster forms via gas fragmentation, runaway mergers can proceed. In this process, the final mass of the central object is constrained by radiation and supernova (SN) feedback onto a collapsing cloud from newly formed stars, since if gas was ejected by feedback, the central object could grow at most to some fraction of the masses of stars (and compact objects) in the cluster.
The main feedback processes from stars are photoionizing UV radiation and/or SN explosions, which can eject gas from the host halo (Kitayama et al. 2004;Whalen et al. 2004;Kitayama & Yoshida 2005). Photoionization feedback from a star influences gas on large scales when the Strömgren radius ion, gas 2 rec,B 1 3 exceeds the effective Bondi Here R i eff,B, is the radius within which ambient gas is bound to the star and is modified from the standard Bondi radius to incorporate radiation pressure to ionized gas (McKee & Tan 2008), a rec,B is the case-B recombination coefficient for H (evaluated at T=10 4 K), n gas is the gas number density, G is the gravitational constant, c s is the sound velocity of gas, Q i ion, is the ionizing photon number flux emitted from a star, m i is the mass of a star, L i is the luminosity of a star, L i E, is the Eddington luminosity, and subscript i represents the ith star in the cluster.
For main-sequence stars of 30 M e in high gas density environments of ∼10 11 cm −3 , the Bondi radius always exceeds the Strömgren radius (Tan & McKee 2004;McKee & Tan 2008;Hosokawa et al. 2012). Thus, photoionization feedback from the low-mass stars of ∼0.1-1 M e , expected to be born from metal-poor gas Dopcke et al. 2013) cannot quench accretion and star formation unless these stars grow to ∼30  M by gas accretion or mergers. Furthermore, photoionization feedback from a massive star becomes efficient only after the star contracts (Hosokawa et al. 2012). Stars typically contract on the Kelvin-Helmholtz (KH) timescale, which is the timescale for a star to radiate away its gravitational binding energy. On the other hand, Hosokawa et al. (2012) and Haemmerlé et al. (2018) have shown that when the accretion rate onto a protostar exceeds a critical rate of ( -)   m M 0.006 0.03 yr cri 1 , the protostar continues expanding, because the heating rate of its envelope due to gas accretion exceeds the radiative cooling rate. The production rate of ionizing photons emitted by the soft spectrum of the bloated star is so low that the gas dynamics is not influenced by photoionization feedback (Kitayama et al. 2004). Sakurai et al. (2015) have found that even when there are quiescent phases of accretion onto protostars, they keep expanding if the timeaveraged accretion rate within the KH timescale (evaluated at the stellar surface) exceeds  m cri . The above suggests that if the growth rate of massive stars by mergers with other stars, averaged on the KH timescale, exceeds  m cri , massive stars would continue expanding for the same reason. The growth of massive stars may remain efficient in this way, until gas is ejected by an SN explosion of one of the massive stars or by accretion feedback from a collapsed massive BH. Thus, there is a possibility that efficient stellar accretion may help to keep the stellar envelope expanding and so inhibit strong feedback from a contracting massive star, thereby leading to the formation of an SMS. In the rest of this paper, unless specified otherwise, the expression "stellar accretion" refers to the central protostar colliding and merging with other stars in the core of the halo.
In this paper, we calculate the evolution of stars that form in high gas density environments as predicted in Omukai et al. (2008). Sakurai et al. (2017) calculated the evolution of stars formed in a massive halo. While they assumed that some fraction of gas is converted to stars at the beginning of the simulation and at the same time gas is ejected, in this paper we consider the evolution of stars including the effects of continuous star formation. Our pathway is similar to the situation in Boekholt et al. (2018), who calculated collisions of accreting stars. In their model, stars are assumed to be kept in the expanded phase due to high gas accretion rates of 0.03  M /yr per star, and dynamical interactions are restricted to the two-body relaxation among stars, while hydrodynamical interactions with gas are neglected. Boekholt et al. (2018) and Reinoso et al. (2018) find that the efficiency of collisions increases as the radii of the stars grow. On the other hand, we find that even when the gas accretion rate is limited by the Bondi-Hoyle-Lyttleton accretion rate, the central star can keep expanding due to accretion of stars.
We find that the central star can grow efficiently due to the following feedback loop. First, stars are captured by the central star by efficient migration due to stellar dynamical friction. The radius of the central star then grows, both because of its increase in mass, and because of the heating of its envelope by stellar accretion. Due to the larger stellar radius, more stars can be captured by the central star. Thus stellar accretion is facilitated by both the mass segregation due to dynamical relaxation processes and the growth of the stellar radius due to stellar accretion. As mentioned above, we refer to this growth process by stellar accretion as stellar bombardment to distinguish it from the usual runaway collisions.
In the usual runaway collisions, only a small fraction of stars in the cluster forms a core and the core collapses. The core is maintained due to the heating by hard binaries, whose binding energy can be a large fraction of the binding energy of the cluster (e.g., Portegies Zwart et al. 1999). On the other hand, during stellar bombardment, the binding energy of hard binaries is a tiny fraction of the binding energy of the cluster, since the central star can be expanding during the evolution, so stars can accrete onto the central star almost without being heated by hard binaries. Thus both the dynamical evolution of surrounding stars and the final outcome (i.e., the mass of the central BH remnant) from the runaway collision and the stellar bombardment scenarios are qualitatively different.

Method
To investigate how stars form, migrate inward, and crash into the central star, and how they are affected by feedback, we use a semianalytical model incorporating the effects of star formation, protostellar evolution, gas accretion, dynamical friction by stars and gas, collisions, and gas ejection (Figure 1). In this section, we provide an overview of our simulations.

Setup and Initial Conditions
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. 3 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 m cent 0,max . As a fiducial model, we set , which is roughly the Jeans mass at which fragmentation occurs at the density -10 cm 10 3 . We assume that there are no surrounding stars initially. When stars are expanding, we set their radius to 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 0.58 (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.

Density Profile
We set the number density profile for gas n gas (r) as 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 11 cm −3 , and the temperature of inflowing gas to T=10 4 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 kT T 10 km s 10 K s 1 2 1 4 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 (e.g., Stahler et al. 1980;Begelman et al. 2006). The core density = n 10 cm c 11 3 roughly matches the density at which gas with a metallicity of ∼10 −5 Z e and a suppressed H 2 fraction (by strong background radiation) begins to fragment due to the decrease of its temperature . In the fiducial model, the initial value of the core radius for the collapsing gas is chosen to be ´r 3 , measured in high-resolution cosmological hydrodynamical simulations of metal-and H 2 -free gas in atomic-cooling halos (e.g., Regan et al. 2014) is also found to be ∼10 −4 -10 −3 pc. We also checked that our results are not significantly influenced by changing the value of r c,ini to 10 −3 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. 4 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).

Gravitational Potential
We adopt a four-component gravitational potential, ( ) F r gas , and ( ) F r cent are, respectively, the gravitational potential at the position r of the dark matter halo, surrounding stars, collapsing gas, and the central star.
We set is the scale radius of the halo, 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 7  M 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 9  M at z∼6 via gas accretion, the seed BH mass is required to be ∼10 5  M at z∼10 (Di Matteo et al. 2012).
The gravitational potential of the gas is derived from Equation where m H is the hydrogen mass. 4 The main effects of the background gas distribution in the simulation are (i) to generate a potential that influences the stellar collision probability (Section 3.5), and (ii) to drive stellar bombardment through gas dynamical friction (Section 3.4.2). The growth rate of the central star is not directly influenced by n gas since we limit the stellar mass accretion rate and star formation by the inflowing gas supply rate from the outer regions  M in (Sections 3.6.2 and 3.2).
Similarly, for the gravitational potential of the surrounding stars: 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 −8 pc to R max =10 pc, spaced uniformly in logr. The gravitational potential by the central star is 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.

Star Formation
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 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 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 where Ω is the orbital frequency of the gas disk, and ( ) ( ( ) ) W = r GM r r Kep 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, W = W  Kep 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 r Q 1 c ,i n i , unless the total accretion rate is as high as  M 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 r µ --r r 1 2 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 r Q 1 c ,i n i , while surrounding stars form uniformly from r Q=1 to r c,ini when < = r r Q 1 c ,i n i . Accordingly the core radius is set to r c =r Q=1 in each time step using Equation (11). From Equation (11), 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 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 ). In our model, gas is heated by gas dynamical friction, with the total heating rate given by We assume that even when gas dynamical friction does not reduce the velocity of surrounding stars at (Section 3.4.2), gas is heated by gas dynamical friction, and we substitute a i GDF, 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 Since the core region is optically thin to dust emission , the net cooling rate by dust grains within the core region is 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).

Expanding and Contracting Stars
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 . 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, , or the time from its formation is shorter than the KH timescale for zero-age main-sequence ), star i keeps expanding, and otherwise it contracts. We calculate the surface KH timescale t i surf,KH, 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 , 6  M m i <50  M , and m i 50  M , respectively. Note that the critical condition (  m 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.

Radial Motion
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 2 is the specific kinetic energy, Dk 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 where a i is the acceleration of the ith particle. We assume Kep is the Keplerian orbital velocity at r. The acceleration is given as where a i SDF, , a i GDF, , and a i acc, 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 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  M 4.1 10 3 ,  M 3.8 10 3 , and  M 3.7 10 3 , 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.

Stellar Dynamical Friction
Stellar dynamical friction is modeled using the analytic formula (Binney & Tremaine 2008) where s * is the velocity dispersion of background stars, r * is the stellar density, ln max min is the Coulomb logarithm, and b max and b min are the maximum and minimum impact parameters for weak stellar encounters. We set (21) assumes an isotropic and Maxwellian velocity distribution for background stars. We set , which sets the value in the square parenthesis in Equation (21) to 0.86 since we assume . 5 In practice we define  á ñ m i t , the smoothed accretion rate of star i at time t, recursively using the instantaneous accretion rate  The Astrophysical Journal,892:36 (19pp), 2020 March 20 Tagawa, Haiman, & Kocsis To obtain the background density r * at each time step, we compute an average stellar density in each radial cell [ ] Î l 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 , 3.7×10 3 , and 3.6×10 3 , 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 (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 m i 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 (Equation (20)) and Δp l =m i v i , we approximate this equation as 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 4 yr. Hence the number of surrounding stars is mostly large enough to validate Equation (21) in our models.

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 where L¢ ln is a Coulomb logarithm for the gas distribution, and v i rel, is the relative velocity between the ith surrounding star and the background gas. Referring to the result of numerical simulations by Chapon et al. (2013), we adopt L¢ = ln 3.1. We 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 ( ) 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.

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 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.

Stellar Collisions among Surrounding Stars
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 s l , * , respectively, the expected rate of collisions within the time step Δt in a cell l is given by (Equation (7.194 2 coll

* * *
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, m l is the average stellar mass in cell l, t l coll, 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 . For describing collisions between two surrounding stars i and j, we assume that the relative velocity v rel, * is drawn from a Maxwellian distribution with the dispersion of s 2 * 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 rel, * exceeds the escape velocity from the stars, 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 v rel, 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  M 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   M m in cent (see Section 3.6.2). For collisions with < v v rel, 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 . 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  where R p is the pericenter distance at collision. Referring to Adams et al. (2004), we set = f 0.16 loss,max 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 loss,max on our results, we compared the final mass of the central star in one of the models ("Model 2" below) with We set the distribution of R p so that b 2 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 R p coll ). The fraction of mass f loss (Equation (28)) is subtracted from the mass of the collided surrounding stars, and added to the inflow rate  M in . Following Bailey & Davies (1999), we also set the condition for merger into a single star to 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. Stellar and Gaseous Accretion onto the Central Star 3.6.1. Stellar Accretion Surrounding stars collide with and accrete onto the central star when the distance from the central star to some surrounding star r i becomes smaller than the sum of the radii After a star accretes onto the central star, we add the mass of the accreted surrounding star to the mass of the central star, and the radius of the central star increases according to Equations (1) or (2). We assume no mass loss during this collision/accretion event. Freitag & Benz (2005) show that when the collision velocity (v coll ) is smaller than the escape velocity from the surface of the collided star (v esc ), the mass loss is small. If the accreted surrounding star orbits in a gravitational potential dominated by the central star, the collision velocity becomes smaller than the escape velocity from the central star. This can be violated and some fraction of the envelope of the central star will be lost if surrounding stars accrete on highly eccentric orbits, which cannot be accounted for in our present model. After accreting a surrounding star, we assume that the envelope of the central star is heated since the orbital energy of the accreted surrounding star is converted to thermal energy in the envelope of the central star. The accreted surrounding star then sinks to the core of the central star, and the central star is expected to expand, similar to the case for gas accretion (Sakurai et al. 2015). We determine the expansion rate of the central star according to the averaged mass accretion rate (Section 3.3). 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  m i acc, onto the central object (in our case the ith star) is given by . We set a = 0.01 SS 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

Gas Accretion
BHL, 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 ith surrounding star is sufficiently high, R i BHL, may become smaller than R i . In this case, the gas accretion rate is determined by direct collision with the stellar surface, 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 3 gas H , i.e., the gas mass within the Bondi-Hoyle-Lyttleton radius (R i BHL, ) of the ith 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. , we reduce the accretion rate  m i by a constant factor f red . We assume =f 10 red 3 as a fiducial value. Even when fragmentation is expected within r i BHL, , we assume that surrounding stars form at r c (Section 3.2).
Thus, in summary we calculate the gas accretion rate of stars as

Feedback Effects
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)  1 2 5 , the gas density is not affected by photoionization feedback. On the other hand, when Q ion exceeds Q cri , gas is blown away from the halo. We adopt this criterion as the gas ejection condition. Q i ion, strongly depends on whether the ith star is in the expanding phase or the contracting phase, with ion cri is ever satisfied, all gas is assumed to be ejected from the system.
Although we add up the ionizing photons emitted by lowmass stars (1  M ), 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 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~10 yr our total calculation timescale of 3Myr, which justifies the assumption of instant gas ejection.

Central Star Evolution
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.
In the early stages, the gas accretion rate onto the central star  m i acc, is as low as ( )  (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 3 yr, the Bondi mass of the central star exceeds its own mass, m cent =31  M (blue and black lines in panel (a)), and the gas accretion rate is reduced by f red =10 −3 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 10 yr 5 , gas accretion begins to dominate the central star's growth rate. In this phase, most of the gas flowing in from large scales  M 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   M 10 5 without any contraction. During the evolution,6.1 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  M 3.5 10 2 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 . 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 . 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  M 6.7 10 5 and6.6 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 880 au c at = t 10 yr 4 is  M 9.9 10 2 , while that of stars is  M 2.1 10 3 . 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) . 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 3 10 cm c 10 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 3 10 KH,ZAMS 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 -10 10 yr 2 3 (e.g., Sakurai et al. 2015), which justifies the assumption of abrupt contraction in our Table 1 The Results of Our Simulations in the Fiducial Model (Model 1) and 21 Variants

Input Output
Model 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 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 ( ( g = G L max HC GDF dust )).
calculations. The central star continues to grow by accreting surrounding stars, but at a more moderate rate, reaching the mass of  M 3.7 10 3 at 3Myr. Therefore a massive BH may still be produced in Model2, but the mass of this massive BH is »100 times below that of the BH remnant in Model 1.
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 selfgravitating 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 (  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)).

Evolution of Surrounding Stars
At2.5 10 yr 3 , 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.
At3.2 10 yr 4 , 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 The growth rate of the central star (black), the rate of stellar accretion onto the central star (red), the gas accretion rate onto the central star (blue), the total star formation rate (orange), and the critical accretion rate below which the central star contracts when the age of the central star exceeds the KH timescale t KH (gray). The black, red, and blue lines are smoothed on a timescale of t surf,KH since the behavior (contraction or expansion) of the central star depends on the growth rate averaged over this timescale. (d) Black and gray lines are the total production rate of ionizing photons and the critical production rate of ionizing photons at which gas is ejected from the halo, respectively. The cyan and blue lines show the cooling rate by dust grains and the heating rate due to gas dynamical friction by surrounding stars, respectively. In this model, the high growth rate of the central star enables it to continue expanding and growing into a supermassive star with 6.7×10 5  M at the end of the simulation at 3Myr. 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. At2 10 yr 5 , 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 10 yr 6 , this star orbits at = r 380 au i , and the mass is 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.

Parameter Dependence
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.
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 (  µ M T in 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 10 cm c 11 3 with T 3 10 3 (Models 1, 3, and 4) the stellar accretion rate onto the central star exceeds the critical rate for contraction before t KH,ZAMS for the central star, allowing it to grow to   M 10 5 . 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 m Bondi cent , which is satisfied for  Table 1), illustrating a case when the central star contracts. The slow decline of the gas accretion rate onto the central star (blue line in panel (c)) after gas dispersal is due to smoothing in time calculated from Equation (16). , and then the gas accretion rate onto the central star exceeds the critical rate  m 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 ej 3Myr), gas accretion contributes very little to the final mass (see the values of m fin and m 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 (  r µ~m r SF c 3 * ). Due to the scaling relations  µm T SF 3 2 and µr n c c 1 3 (Equation (12)), we can expect r µ~T n 3 2 c * . Figure 6 shows the relation between the final  The binding energy of the star (black) and the cumulative change in the kinetic energy due to the stellar dynamical friction (orange and blue), accretion drag (cyan), and gas dynamical friction (brown). Decrease and increase of the kinetic energy by stellar dynamical friction are shown separately by orange and blue lines, respectively. This star decreases its mass due to frequent collisions, and ends up as a less massive star orbiting at ∼380 au at 3Myr. mass for the central star and the product T n 3 2 c . We can indeed see the rough correlation between the final mass and T n 3 2 c as expected in the cases in which the central star contracts (dashed line and circles in Figure 6). Also, the expected relation r µ~T n 3 2 c * is roughly confirmed in Figure 7, which shows the maximum stellar density within the core radius as a function of T n 3 2 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.
As discussed above, if T is low, r * remains low, and therefore stellar dynamical friction is inefficient. Additionally, since the mass of the stellar cluster in the core is approximately limited by  M t in , if T is low, then  M 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 n 3 2 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 3 Myr fin 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 G > L GDF 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 G L  GDF dust .

Discussion
In this section, we discuss assumptions in our models.

Inconsistency between Assumptions
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 3 to 10 5 yr in Model 1, and never reaches very small values. Second, we expect that an isotropic distribution is established by relaxation   Figure 6, but the y-axis represents the maximum stellar density within the core radius.
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~W h c s 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.

Star Formation Efficiency
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 yr 4 in our fiducial Model 1 (see below), and it increases with time. Observationally, some massive molecular clouds are found to have an SFE of >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 ´M 2 10 6 , and the mass of the stellar cluster is at most  M 10 4 (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 4 K, the star formation rates for T=5×10 3 , 3×10 3 , and 10 3 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.

The Eccentric Orbit
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 e de 2 i i i (e.g., Jeans 1919;Heggie 1975 ). 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(Hayasaki et al. , 2016Bonnerot & 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.

Evolution Following the Formation of the Massive Black Hole
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 f t l l evap, evap relax, (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), * * is the relaxation timescale (Binney & Tremaine 2008), and we set the Coulomb logarithm to be 10. Although we set = f 300 evap , 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. 7 Figure 4 shows that the collision and evaporation timescales in the outer regions of the cluster are longer (at 2000 au where - r~-M 10 pc 7 8 3 * ) and comparable to the Hubble time of 10 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  M 10 3 . 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 (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).

Comparison with Hydrodynamical Simulations
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 (  - M 1 yr 1 ). They find that the power-law index of the initial mass function is ∼−1, the maximum mass of surrounding stars is , β=−1, and =T 3 10 K 4 . Since they use a barotropic equation of state, we allow the star formation even when G > L GDF dust . Since f red is highly uncertain, we varied f red between 10 −3 and 1.
The evolution of = f 0.1 red is shown in Figure 8. The central star grows to  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 0.001 red , 0.01, 0.1, and 1, respectively, the mass of the central star at 10 kyr is 5.5×10 3 , 5.8×10 3 , 1.0×10 4 , and 1.1×10 4  M , 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 4  M , 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 ) for Z∼10 −3 -10 −5 compared to the value adopted in our models (0.08  M ). 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 −3 -1 is as high as 3.5×10 6 due to the high inflow rate (   -M M 1.1 yr in 1 ).

Conclusions
In this paper, we propose a process for forming supermassive stars via stellar collisions and accretion in high-redshift protogalaxies. The scenario envisioned here shares some aspects of both the popular "direct collapse" and the "runaway collision" scenarios. We focus on environments in which a gas cloud is polluted only by a moderate amount of metals, and its H 2 abundance is suppressed. In such environments, a gas cloud fragments only at very high density, producing a high-density stellar cluster ). If gas is ejected soon after stars form, the final mass of a central star becomes ∼10 3  M (Sakurai et al. 2017).
The novel aspect proposed here is that if subsequent frequent capture and accretion of surrounding stars onto a central star efficiently heats the envelope of the central star, the central star continues expanding, and gas will be retained in the system due to the lack of strong UV radiation and weak photoionization feedback from the bloated central star. The central star can therefore keep growing until the supply of surrounding stars and gas runs out due to gas ejection by SN explosions or by accretion feedback from a collapsed massive BH. We call such a rapid stellar accretion process "stellar bombardment," which could be caused by efficient stellar migration via relaxation processes, the increase of the stellar radius by the mass increase, and most importantly, the heating and bloating of the stellar envelope due to the frequent stellar accretion itself.
To investigate the viability of this "stellar bombardment" scenario, we have performed numerical modeling using a semianalytic toy model. The model includes dynamical friction by stars and gas, star formation, gas accretion, collisions, and gas ejection. Our main results can be summarized as follows: 1. When the central core density exceeds 10 11 cm −3 and the gas temperature is 3×10 3 K, the central star continues growing without contracting until it reaches a mass of ∼10 5 -10 6 M e at 3Myr. The central star grows mainly by stellar bombardment early on, and by gas accretion in the later phases. 2. When the central core density is below 3×10 10 cm −3 , the central star contracts due to the subcritical rate of accretion and heating by surrounding stars. After the contraction, photoionization feedback ejects gas from the system, reducing the final central star mass by about two orders of magnitude, to 10 4 M e . 3. The final mass of the central star depends strongly on the gas temperature and the core density of the gas, in addition to whether the central star contracts (Figure 6). This is because the efficient growth of the central star by stellar accretion requires a high-density cluster. Highdensity star clusters can be realized for high star formation rates and/or compact core sizes, which in turn are produced for high gas temperature and core gas density, respectively. In a cosmological setting, these conditions can arise in metal-poor atomic-cooling halos, in which the H 2 abundance has been suppressed, leading to inefficient cooling until very high densities are reached.
In this paper we have used a simple toy model to illustrate the possibility of this new evolutionary process. To understand this pathway in more detail, including its viability, future N-body and hydrodynamical simulations will be required, which are able to follow stellar evolution and radiation feedback onto the collapsing cloud. Chon & Omukai (2020) have recently presented results from three-dimensional hydrodynamical simulations, albeit with a prescribed barotropic equation of state, and without radiation. They find that supermassive stars likely form from gas with metallicities up to ∼10 −4 Z e due to stellar accretion and gas accretion. Thus our predictions are confirmed by more realistic numerical simulations.