© The Authors 2025

1. Introduction

More than 200 supermassive black holes (SMBHs) have been detected at z > 5.7 (Fan et al. 2006; Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018, 2021; Reed et al. 2019; Onoue et al. 2019; Wang et al. 2021), with masses higher than 10⁹ M_⊙. This also includes SMBHs as massive as 10¹⁰ M_⊙ for SDSS J010013.02+280225.8 (Wu et al. 2015) or from the recent new observations of the James Webb Space Telescope (JWST) in collaboration with the Chandra X-ray Observatory, which is the currently highest-redshift quasar UHZ1 at z = 10.3 with a mass of about M_bh ∼ 10⁷ − 10⁸ M_⊙, as reported by Bogdán et al. (2024). The possible origin of these objects was discussed early by Rees (1984), who also discussed a variety of possible formation pathways, including the possibility of a direct collapse, in which a massive gas cloud collapses directly into one or a small number of objects (e.g., Bromm & Loeb 2003; Koushiappas et al. 2004; Wise et al. 2008; Begelman & Shlosman 2009; Schleicher et al. 2010; Latif et al. 2013). A viable alternative for the formation of SMBHs is the formation of very massive objects through collisions in dense stellar clusters (Omukai et al. 2008; Devecchi & Volonteri 2009; Katz et al. 2015; Sakurai et al. 2017, 2019; Reinoso et al. 2018, 2020; Vergara et al. 2021, 2022, 2023; Escala 2021) or black hole mergers in dense black hole clusters (Davies et al. 2011; Lupi et al. 2014).

Another frequently discussed possibility is that the first SMBHs may originate from the first stars in the Universe (Abel et al. 2002; Heger & Woosley 2002; Heger et al. 2003; Klessen 2019) through their top-heavy initial mass function (Clark et al. 2011; Greif et al. 2011; Susa et al. 2014; Stacy et al. 2016; Fraser et al. 2017; Riaz et al. 2018; Sharda et al. 2020). In comparison to the other scenarios these seeds would be relatively light, however. Another challenge is that these stars very effectively removing gas from the low-mass halos through radiative and supernova feedback. This impedes subsequent growth (Johnson & Bromm 2007; Smith et al. 2018).

More recent investigations have also suggested that a realistic pathway for the formation of a very massive object may not only consist of the idealized scenarios originally proposed by Rees (1984), but that several variations are also possible, including hybrid scenarios that combine gas and stellar dynamical effects. For example, while cosmological simulations found that fragmentation occurs in gas dynamical studies, the resulting clumps may also merge again and do not necessarily inhibit the formation of a very massive object (see, e.g., Grete et al. 2019; Suazo et al. 2019; Regan 2023; Regan & Downes 2018a,b; Prole et al. 2024). It is well known that gravitational torques in self-gravitating disks may produce exactly this type of behavior (Latif & Schleicher 2015). Similarly, the very high accretion rates that were often found in gas dynamical simulations are expected to enhance the radii of the protostars, which in turn enhances their probability for collisions (e.g., Hosokawa et al. 2013; Schleicher et al. 2013; Haemmerlé et al. 2018, 2019, 2021).

A more systematic exploration of these hybrid scenarios has been pursued by Boekholt et al. (2018), who presented toy models to study the interplay of accretion and collisions. The authors showed that the latter could lead to the formation of very massive central objects in dense compact clusters. Alister Seguel et al. (2020) included the mass loss during collisions of protostars in these simulations. Tagawa et al. (2020) developed a semianalytic model considering stellar bombardment and accretion, while Das et al. (2021) showed that the collision probability is enhanced in the presence of accretion through the absence of mass and linear momentum conservation.

The question regarding the origin remains very important today. The original formation needs to be taken into account together with subsequent accretion and other processes, such as the growth via merger events (see, e.g., the reviews by Volonteri 2010; Woods et al. 2019). Semianalytic models that include these processes are very helpful because they provide an assessment of the black hole population that results from different assumptions on the different seeds. This provided some tentative indications that hybrid models that might work at very low but nonzero metallicities may be an important ingredient to consider (Habouzit et al. 2016; Sassano et al. 2021; Trinca et al. 2022). Chon & Omukai (2020) presented simulations that incorporated dust cooling and showed that the resulting fragmentation does not have to impede the formation of a massive object. An analytic scenario showing that collisions and accretion can be expected to interact and form massive objects was presented by Schleicher et al. (2022). Reinoso et al. (2023) presented a suite of numerical simulations that consistently modeled the gas and N-body dynamics of the protostars in a compact cluster, demonstrating in detail how these processes interact and lead to the formation of very massive objects.

We follow up on the recent study by Reinoso et al. (2023), with the objective to study even more unstable systems, that is, systems with a higher ratio of the gas mass to the Jeans mass, because we expect this regime to be particularly relevant and important for the formation of very massive objects. For this purpose, we varied the gas temperature from 8000 K, where atomic hydrogen cooling dominates, to 500 K, where H₂ cooling and dust thermal emission become efficient and where gravitational instability might be increased. Our numerical approach is similar to that of Reinoso et al. (2023) and is briefly presented in Section 2 for completeness. Our results are given in Section 4, and we discuss our main conclusions in Section 5.

2. Numerical method

To explore how the gas-dynamics and protostellar dynamics affect the formation of a SMBH seed we model a protostellar cluster inside a compact gas cloud in virial equilibrium. The methodology adopted here is based on the Astrophysical MUlti-purpose Software Environment (AMUSE¹, see Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013; Portegies Zwart & McMillan 2018), a PYTHON interface designed to couple existing numerical codes. We closely follow the methodology presented by Reinoso et al. (2023), including mass accretion onto protostars, sink particle creation, a treatment for stellar collisions, and a mass-radius relation for the protostars.

The N-body code PH4 (McMillan & Hut 1996) and the smoothed-particle hydrodynamics (SPH) code FI (Hernquist & Katz 1989; Gerritsen & Icke 1997; Pelupessy et al. 2004) are coupled via the BRIDGE method (Fujii et al. 2007) in AMUSE. The gravitational acceleration at the position of the N-body particles is calculated using the SPH particles and vice versa. To ensure that the coupling does not violate Newton’s third law, we employ the code FastKick with a constant gravitational smoothing length of 0.5 au corresponding to the smallest smoothing length in the simulations. The particles in the N-body code are solved with a smaller smoothing length of 1 R_⊙ for a very accurate treatment of the gravitational interactions between the protostars. An external pressure floor is included in the SPH code using the treatment of Benz (1990), Clark et al. (2011). The external pressure is set equal to the pressure of the cloud at the cut-off radius, corresponding to 1.78 × 10⁻⁷ g cm⁻¹ s⁻² for the hottest simulation and 1.13 × 10⁻⁸ g cm⁻¹ s⁻² for the coldest one.

We employ a modified equation of state of the form

$$\begin{array}{c}\hfill T=T_0[1+{\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)}^{\gamma -1}],\end{array}$$(1)

such that the gas behaves isothermally, with temperatures from T₀ = 500 K − 8000 K at low densities, but becomes adiabatic at densities above ρ_c = 10¹⁵ cm⁻³, as found in 1D and 3D models including detailed chemical networks with γ = 5/3 (Omukai et al. 2008; Becerra et al. 2015).

The time integration is based on the Kick-Drift-Kick (KDK) integration with bridge including protostellar collisions and sink particle creation using the algorithms described below. We employ a maximum timestep of 5 yr for bridge along with an adaptive time-stepping algorithm, where the timestep is reduced by factors of 2 until becoming smaller than the shortest free-fall time among all the SPH particles. The N-body particles in our simulations represent protostars. Using the particle sets in AMUSE, we assign additional properties such as: stage and luminosity. The stage property indicates the phase of stellar evolution, that is, when the object has already entered the main sequence phase or if it is still in the protostellar phase, allowing to take into account the difference in the protostellar radii, which is highly relevant for collisions. We use a mass-radius relationship according to Reinoso et al. (2023), based on the works by Hosokawa & Omukai (2009), Hosokawa et al. (2012, 2013). The N-body particles evolve through three different stages: protostars, star and supermassive star (SMS). When a particle is in the star stage, it grows as a main sequence star according to the following mass-radius relation:

$$\begin{array}{c}\hfill R_{\mathrm{star}}=0.97{\left(\frac{M_{\mathrm{star}}}{{\mathrm{M}}_{\odot }}\right)}^{0.57}{\mathrm{R}}_{\odot }.\end{array}$$(2)

Particles in a SMS stage have an accretion rate greater than the critical accretion rate ṁ_crit = 0.04 M_⊙ yr⁻¹ (Hosokawa et al. 2013). The mass-radius relation of a particle in a protostar stage is determined by the accretion rate ṁ, where the swelling or contraction of the protostar will bring it closer to either the star or SMS stage, depending on the accretion rate. A protostar can contract to a star stage if its accretion rate falls below the critical accretion rate for a period of time longer than a factor 100 of the Kelvin-Helmholtz timescale (100t_KH). For a detailed description of the mass-radius relation, see Appendix A in Reinoso et al. (2023).

After each KDK step, we calculate the gas accretion onto the protostars employing the algorithm from Hubber et al. (2013). We consider a spherical volume called the interaction zone with radius R_I.Z.. Inside this radius the weighted average gas flux onto the central point mass is computed using a cubic spline kernel function. The interaction radius is adjusted iteratively with a maximum of 50 iterations to maintain a constant gas mass of M_int, max = 50 M_gas/N_SPH (the mass corresponding to 50 SPH particles). We further impose that the interaction zone cannot be smaller than the protostars. We employ angular momentum conservation for the accreted gas and include the angular momentum feedback described by Hubber et al. (2013), including their prescriptions for accretion. We use their NewSink algorithm to create sink particles when the following conditions are fulfilled: (i) an SPH particle reaches a density higher than 10¹⁶ cm⁻³ (Becerra et al. 2015; ii) it sits in a minimum of the gravitational potential among its neighbors; (iii) it does not overlap with existing sinks; (iv) it fulfills the density criterion:

$$\begin{array}{c}\hfill {\rho }_i>{\rho }_{\mathrm{Hill}}\equiv \frac{3X_{\mathrm{Hill}}(-\mathrm{\Delta }{r}_{is\prime }·\mathrm{\Delta }{a}_{is\prime })}{4\pi G|\mathrm{\Delta }{r}_{is\prime }{|}^2},\end{array}$$(3)

for all existing sinks s′ for a given SPH particle i. Here, Δr_is′ and Δa_is′ are the relative position and acceleration of sink candidate i with respect to an existing sink s′, with X_Hill = 4. The last criterion is based on the Hill sphere and ensures that an SPH particle turns into a sink particle in the vicinity of another sink only if the density peak dominates the local gravitational potential.

A collision between two particles occurs when their radii overlap during the N-body integration, that is,

$$\begin{array}{c}\hfill d\le R_1+R_2,\end{array}$$

where d is the separation between the particles and R₁ and R₂ are their radii. We then replace the overlapping particles with a new particle placed at the center of mass where the new velocity is calculated following momentum conservation. The new mass of the particle is the sum of the masses of the colliding particles, that is,

$$\begin{array}{c}\hfill M_{\mathrm{new}}=M_1+M_2.\end{array}$$

3. Initial conditions

The gas is modeled through SPH particles, with a total of 2¹⁸ SPH particles per simulation. In all simulations the mass of the gas cloud is M_gas = 3 × 10⁴ M_⊙, with a virial radius of 0.14 pc. The protostars were simulated by sink particles that interact only through gravity with themselves and with the gas. The initial number of sink particles is 256, each with a mass of 0.1 M_⊙, which are compatible with protostellar masses that emerged from atomic cooling halos (Becerra et al. 2015). Both the protostars and the gas particles follow a plummer distribution (Plummer 1911), with a half-mass radius of 0.1 pc and a cut-off radius of 5 Plummer radii for each model. The initial velocity of the gas particles is zero and for the protostars the velocities follow a virial equilibrium condition. To be consistent with the results found by Latif et al. (2013), we use the same recipe as Reinoso et al. (2023) to inject a spectrum of noncompressive Kolmogorov turbulence with Mach number ℳ = 1. We simulate isothermal gas clouds with different initial temperatures and we assume a primordial composition of 75% atomic hydrogen and 25% helium. In Table 1 we summarize the initial conditions for all simulations performed in this study.

4. Results

In this section we present the main results of our simulations, adopting the initial conditions from Table 1, where we vary the initial temperature to explore its impact on the formation of the most massive object (MMO) in the cluster. The simulations are stopped at approximately ∼30 000 yr, because at this point most of the gas has already been accreted and a supermassive object at the center of the cluster is clearly established.

4.1. General evolution of the cluster

The shape of the collapsing gas cloud after the initial collapse is largely sensitive to the cooling and heating processes, and thus by the gas temperature, which is likely to have a significant impact on fragmentation. Fig. 1 shows the face-on view of our simulations, where we present the morphology at different times of our simulations with temperatures from T = 8000 K to T = 500 K (edge-on view can be found in Appendix A). The maximum gas density reaches ρ ≈ 5.08 × 10⁻⁶ g cm⁻³ for the hottest simulations and ρ ≈ 5.31 × 10⁻³ g cm⁻³ for the coldest simulations. In all temperatures we observe a similar initial evolution with a rapid collapse of the central gas cloud, where 25% of the total mass collapsed between ∼3000 − 4000 yr, thereby producing a larger potential gap in the center of the cluster leading to high accretion rates, runaway collisions and star formation. In nine of our simulations, the MMO was a new sink formed in this central overdense region, whereas in the other run, a sink particle from the initial configuration became the MMO.

Fig. 1. Density projections of the face-on view of a primordial gas cloud with an embedded protostar cluster for simulations with an initial temperature of 500 K − 8000 K at different times. The white dots are the protostars, black dots are the formed stars and the white star is the MMO. The clusters evolve in time from left to right and from the warmer to colder simulations from the top to the bottom panels.

After the formation of the MMO (marked as a white star), we observe a comparable overall morphology of the collapsing clouds, maintaining a spherical shape at the beginning and later forming a compact filamentary structure in the center. This elongated overdensity was produced from a disk structure around the MMO, formed as a result of angular momentum conservation.

Radial profiles of several quantities related to the general properties of the system are presented in Fig. 2. These profiles correspond to the third snapshots shown in Fig. 1, starting with the gas density in the top panel, stellar mass density excluding the MMO in the second panel, radial velocity of the gas in the third panel, and an analytic accretion rate given by $\dot{m}=4\pi r^2\rho (r)v_{\mathrm{rad}}$ in the fourth panel. In the top panel, we observe a general trend in the shape of the density profiles for the different temperatures. At the center, the gas reaches similar high densities, while at larger radii, they become flatter as the temperature decreases. In the case of T = 500 K, the density decreases strongly because most of the gas has collapsed. In the second panel, we see that the stellar mass density is higher for warmer temperatures. This is because higher temperatures reduce fragmentation and higher accretion rates from companion stars can be achieved, allowing the formation of more massive companion stars. The third panel in Fig. 2 indicates that the gas is not in equilibrium around the MMO, because the gas is collapsing in some regions while expanding in others. Even if some gas evaporates, however, it does not reach the escape velocity of the MMO, keeping the same gas mass or increasing it, if the gas continues to flow toward the MMO. In the bottom panel, we show the radial accretion rate, considering only regions when the radial velocity is toward the center (V_rad < 0). Close to the MMO, the accretion rate drops to approximately ∼10⁻² − 1 M_⊙ yr⁻¹ and remains higher at larger radii for colder temperatures, following the trend of the radial density profile.

Fig. 2. Radial profiles relative to the MMO for the gas density, stellar mass density without the MMO, radial velocity, and accretion at approximately 15 000 yr, corresponding to the third snapshot of each temperature shown in Fig. 1.

In our simulations, we implicitly assume the presence of metals in the gas, which can cool the gas to temperatures on the order of 100 K. Omukai et al. (2008) studied the thermal and chemical evolution of low-metallicity gas clouds exposed to UV radiation. Their one-zone model results indicate that strong UV radiation can suppress H₂ formation in a low-metallicity gas cloud, leading to a high-temperature evolution. In the presence of higher metal abundances, the temperature can decrease to ∼200 K, primarily due to dust thermal emission. Gas evolution with lower temperatures can trigger an early fragmentation of the gas cloud that could generate more sink formation and convert the gas mass into protostellar mass more quickly as we can see in Fig. 3, where the blue, green and orange regions are the mass from the gas, protostars and escaped protostars of the system. In Fig. 3 we observe a faster reduction of the gas mass for simulations with lower temperatures, with steeper slope and a displacement of the black dashed lines to the left, indicating that more gas mass was accreted at a given time due to gas accretion or star formation. We find in general that the escaped mass from the system can be ignored, since even in the most extreme case the cluster just loses ∼0.007% of the total mass.

Fig. 3. Time evolution of the mass fraction for the simulations with an initial temperature of 8000 K, 5000 K, 3000 K, 1000 K, and 500 K, corresponding to the gas mass (blue zone), protostars in the cluster (green zone) and the escaped protostellar mass (orange zone). The vertical dotted black lines indicate the times when the protostars accreted 25%, 50%, and 75% of the gas mass.

In all simulations, an MMO is formed and most of the gas is accreted by the stars by the end of the simulations. In all our simulations we find a similar general behavior that can be divided into three phases: (i) rapid initial contraction, (ii) accelerated formation and predominance of the MMO and (iii) a small cluster with the MMO at its center.

4.1.1. Initial contraction

The first phase is marked by a rapid initial contraction of the cluster, triggered by the initial instability with a ratio of M_gas/M_Jeans = 10^0.73 for the hottest simulations and M_gas/M_Jeans = 10^2.54 for the simulations with an initial temperature of 500 K. These more unstable configurations for lower temperatures provide a favorable scenario for high accretion rates and runaway collisions.

This is illustrated in Fig. 4, where we present the Lagrangian radii of 10%, 25%, 50%, 75% and 90% of the total mass of the system relative to the center of mass for different temperatures. For all temperatures, we observe a rapid core collapse before ∼5000 yr, during which around 50% of the gas mass was accreted by the stars, as we shown in Fig. 3. This is caused because the gas cloud departs from equilibrium due to Jeans instability where M_gas > M_Jeans, and the collapse is isothermal, implying that the pressure forces cannot significantly retard the gas collapse. Considering an initial density ρ as a function of the initial radius r₀, we have mass shells within the initial radius r₀ with a mean density of ρ_m(r₀, t₀) in a time t₀. Assuming that the different shells do not cross each other and as the free-fall timescales as ρ^−1/2, we will have different free-fall times for different mass shells becoming shorter for the inner regions, since ρ_m(r, t) decreases with increasing r. Thus, the inner regions will collapse first, with the other mass shells falling in progressively later. This behavior is best reflected for colder temperatures, where the initial instability is strong enough to collapse most of the gas cloud. In contrast, for warmer temperatures, we observe a higher probability that the MMO may move away from the center of the system, making the accretion less efficient.

Fig. 4. Time evolution of the 90%, 75%, 50%, 25%, and 10% Lagrangian radii for our simulations with temperatures of 8000 K, 3000 K, and 500 K. The panels show the Lagrangian radii considering the total mass of the systems (gas and stars) relative to their centers of mass.

4.1.2. Formation of a central object

After the initial collapse, we observe a high density in the central region of the gas cloud, where the sink forms, which later evolves into the MMO. In nine of our simulations, the MMO was a new sink formed during the first collapse, when the mean density in the warmest simulations within the 25% Lagrangian radius was ∼1.11 × 10⁻¹⁰ g cm⁻³ and 1.89 × 10⁻¹⁰ g cm⁻³ in the coldest simulations. The distance from the center of the system to the position where the MMO was formed is ∼0.00348 ± 0.00166 pc, so it is in the inner region of the gas cloud with enough gas to accrete and increase in mass. After the birth of the protostar, it finds itself in a rich environment where it will reach high accretion rates and grows rapidly in mass, from 9.65 ± 4.13 M_⊙ to 5636.53 ± 858.14 M_⊙ in only a thousand years, with a maximum peak in the accretion rate of 20.56 ± 12.67 M_⊙ yr⁻¹, due to the accretion of a massive clump formed by jeans instability (Hosokawa et al. 2016).

4.1.3. High accretion rates

While initially the mass of the MMO is smaller than the total mass in the rest of the protostars, eventually the MMO becomes more massive in all the simulations. This transition occurs when the MMO reaches a mass of approximately 645.69 ± 303.81 M_⊙. The second phase is marked by very high accretion rates in the MMO, with a mean accretion rate of 1.253 ± 0.864 M_⊙ yr⁻¹ for warmer simulations and 2.928 ± 0.271 M_⊙ yr⁻¹ for colder ones. These high accretion rates were reached before ∼6000 yr, due to the collapse of the gas toward the center of the cluster (see Table 2).

Fig. 5 shows the gas mass fractions and the protostellar mass fractions (top panels), the mass of the MMO and the total mass of the remaining protostars (second panels), accretion rate (third panel), radius of MMO and protostars (fourth panel) and the total collision rate (bottom panels). In the first panels we observe that in the evolution of the total mass of the gas and the protostars, the intersection of the curves occurs at roughly 5219.58 ± 705.46 yr, which corresponds to the moment when the dynamics of the system is dominated by the protostars rather than the gas. This point occurs ∼1855.27 yr earlier for simulations with 8000 K compared to simulations with 500 K. Later we see a difference in the mass exchange between the gas and stars, which is faster for more unstable configurations or colder temperatures. It could be argued that the formation of the MMO is not very sensitive to the initial gas temperature, but the growth could be highly sensitive to the initial instability, driving a faster growth for more unstable configurations. The relatively rapid occurrence of the mass transition also suggests that for our configuration the gas plays a role mostly at the beginning, while subsequently the collision dynamics is more relevant. In the second panels of Fig. 5, we find the formation of a single massive object (blue line) which consumes most of the available gas mass, while the total mass of the remaining protostars (orange line) reaches a few thousand solar masses. We also note that in most of our simulations, the MMO is a protostar formed due to the gas collapse and not from the protostars that were part of the initial condition. This is because the first component to collapse is the gas and it will concentrate a lot of mass in the center of the cluster, allowing the formation of some new sink particles within a high density region, where the overlap criterion for new sink formation should not be a problem due to the high number of SPH particles in the region.

Fig. 5. Evolution of the gas and protostellar mass in the simulations with different temperatures (top panel). Comparison of the time evolution of the MMO mass and the total mass of the protostars without the MMO (second panel), accretion rate of the MMO with the mean accretion rate of the protostars without the MMO (third panel), where the black dashed line is the critical accretion rate ṁcrit = 0.04 M⊙ yr−1 and marks the point where the protostars pass to the SMS Track. We further show the radius of the MMO with the mean radii of the protostars (fourth panel) and the evolution of the total collision rates over time (bottom panel).

We observe in the third panels of Fig. 5 that the initial phase of the evolution is marked by high accretion rates of the MMOs, on average with 1.32 ± 3.35 M_⊙ yr⁻¹ at roughly 5000 yr and 1.88 × 10⁻³ ± 0.01 M_⊙ yr⁻¹ for the other protostars at a similar time. The high accretion rates of protostars are caused by the strong gas inflow to the center of the cluster, triggered by the initial gas instability. These high accretion rates rapidly make one of the protostars the MMO in the cluster, leading to a swelling of its radius to ∼10⁴ R_⊙, increasing its probability of collisions very early in the simulation. As a result, the system tends to form a single massive object, with a large mass difference between the MMO and the other protostars in the cluster.

This swelling of the MMO may also explain the high number of collisions observed in the fifth panels of Fig. 5, where the runaway collision period is, at least partly, a result of the increased radius of the central sink particles and the MMO. The initial high number of collisions is due to the increase in the stellar radius of the MMO, triggered by the initial high accretion rates, which increases its cross section as σ ∝ r². The high collision rate is maintained due to the star formation as a result of disk fragmentation around the MMO.

4.1.4. MMO and protostellar cluster

After the growth phase the cluster reaches its final state, where the typical outcome is a small cluster of stars with minimal gas remaining. The typical star masses in the stellar cluster are in the range of 7.35 ± 20.23 M_⊙, with 81.20 ± 35.78 members surrounding the MMO. Up to that time most of the gas was accreted by the stars and the MMO, reaching an efficiency of ∼0.62 for the hottest simulations and ∼0.96 for the coldest one. We have not considered a binary outcome in our results because all second MMOs are less than 400 M_⊙. In Fig. 6 we present the combined mass distribution of the stars bound to the MMO in the final state of our simulations. We note that for the simulations with 8000 K there are two trends. The first with a prominent peak that corresponds to stars with masses around ∼10 M_⊙ and with less presence, stars in the lowest limit with ∼10⁻¹ M_⊙. For colder simulations, we observe a reduction of the most massive peak and there are more stars close to the lowest mass peak. The combined final number of stars is also smaller for colder simulations. This is because if we reduce the gas temperature we will have a more unstable configuration, making the gas collapse faster toward the center of the cluster, leaving the stars in a gas-poor environment for accretion. These poor conditions could inhibit the mass growth of protostars, even though they were accreting from the beginning of the gas collapse. When we compare the ratio between the mean accretion time of the protostars that reach the final time with the free-fall time, we find that in the warmest simulation, approximatly ∼50% of the stars have t_acc/t_ff < 1, and for lower temperatures, t_acc/t_ff ≫ 1 for all final stars. This is shown in Fig. 7, indicating that the gas collapses faster than the protostars could accrete it.

Fig. 6. Combined mass distribution at the end of our simulations for each temperature.

Fig. 7. Mean accretion time over the free-fall time for each star in the cluster at the end of a simulation with each temperature.

As a result of the reduced temperature, the Jeans mass is also reduced, where $M_{\mathrm{Jeans}}\propto T^{\frac{3}{2}}$. This will produce more fragmentation with smaller masses, reducing the accretion for less massive protostars. Hence, due to the reduced accretion in a gas-poor environment, the stars did not grow large enough in radius to increase their cross section and grow through collisions, so most of the available gas will collapse and be consumed by the MMO.

Our simulations agree with the supercompetitive accretion mechanism proposed by Chon & Omukai (2020) and not with a fragmentation-induced starvation because t_acc/t_ff ≫ 1, and the gas collapses toward the center of gas cloud, allowing the formation and evolution of a single MMO.

In Fig. 8 we show the mass of the MMO over the total stellar mass without the MMO, where according to theoretical work by Schleicher et al. (2023), the ratio M_MMO/M_stars must be a factor of 10⁻² to limit the accretion of the MMO. We only see this regime when the MMO was formed, but quickly grows above 1 and in most of the evolution for all temperatures it is greater than 1, therefore the stellar mass is not large enough to limit the accretion of the MMO and reach the fragmentation-induced starvation regime.

Fig. 8. MMO mass over the total stellar mass without the MMO for each temperature. The horizontal dashed black line shows when Mstars/MMMO = 1.

In Fig. 9 we show the combined mass distribution for all ejected stars at each temperature. For all temperatures we have a similar quantity of ejected stars mainly due to the three-body interactions. The ejected stars show a similar behavior in the peak movement as in Fig. 6, where the peak moves from more massive stars to lighter stars, with 9.38 ± 4.77 M_⊙ for the simulations with 8000 K and 0.87 ± 2.10 M_⊙ for 500 K simulations.

Fig. 9. Combined mass distributions of the escaped stars for simulations with all temperatures.

4.2. Sink creation mechanism and disk fragmentation

As we mentioned in the previous Section 4.1.3, one possible mechanism for the creation of new sink particles could be the disk fragmentation. In Figs. 1 and A.1 we observe a disk-like morphology, with a long overdensity around the MMO and both protostars and new sinks are rotating around it. The disk fragmentation implies that the star formation was mainly due to overdensities produced in a turbulent disk structure around the MMO. There are previous numerical studies of disk fragmentation around a massive star in the context of star formation (Oliva & Kuiper 2020; Girichidis et al. 2012; Peters et al. 2010) and observational evidence of a disk formation around massive stars (Maud et al. 2019; Ginsburg et al. 2018). Making the disk fragmentation a possible mechanism to explain the high fragmentation and protostar collisions with the MMO.

The simplest configuration to form a disk structure is a rotating and axisymmetric cloud without magnetic field. Therefore, the gas needs to reach the following conditions:

• Axisymmetric gas cloud, with one component of the angular momentum must predominate.

• Rotating gas cloud, where the rotational kinetic energy must be higher than the gravitational potential energy.

If the gas reaches these conditions, it could form a disk structure around the MMO. In our simulations, we initially impose noncompressive Kolmogorov-turbulence, which induces a net angular momentum of the gas going into gravitational collapse, forming a disk structure around the MMO as a result of angular momentum conservation. After some time, due to the overdensities in the rotating turbulent gas, we observe the formation of sink particles around the MMO.

The conditions of axisymmetric and rotational gas cloud are shown in Fig. 10, with the ratio between the rotational kinetic energy and potential energy of the gas in a region of 0.01 pc around the MMO and the angular momentum in spherical coordinates. In the top panel we observe that the gas needs some time to increase its rotational energy and to be above the potential energy in order to form a disk structure, needing more time for lower temperatures because the radial velocity toward the center of the gas cloud for colder simulations is higher due to the higher initial instability. For all simulations with temperatures over 500 K the kinetic rotational energy begins to be greater than the potential energy near ∼11 000 yr, and for a temperature of 500 K the rotational energy is above the gravitational energy for a short time, in a time around ∼20 000 yr. Making it more likely that the gas will rotate earlier around the MMO for higher temperatures.

Fig. 10. Ratio of the rotational kinetic energy and potential energy of the gas in a region of 0.01 pc around the MMO for each temperature (top panel). The angles ϕ and θ of the total angular momentum are plotted in spherical coordinates (middle and bottom panels).

Another necessary condition for disk formation is a predominant angular momentum component. In the middle and bottom panel we show the ϕ and θ angles of the total angular momentum of the gas in a region of 0.01 pc. In the ϕ component we do not observe a lot of variation, being constant most of the time and therefore no azimuthal rotation is observed. On the other side, the θ component shows more instability for simulations of 1000 K and 500 K, needing more time to reach a stable configuration. In general, we observe a predominant direction on time of the total angular momentum of the gas around the MMO.

We also calculated the components of the inertia tensor relative to the MMO in the following way:

$$\begin{array}{cc}\hfill I_{\mathit{xx}}& ={\displaystyle \sum _j}{m}_j(\mathrm{\Delta }{y}_{j,\mathrm{MMO}}^2+\mathrm{\Delta }{z}_{j,\mathrm{MMO}}^2),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill I_{\mathit{xy}}& ={\displaystyle \sum _j}{m}_j\mathrm{\Delta }{x}_{j,\mathrm{MMO}}\mathrm{\Delta }{y}_{j,\mathrm{MMO}},\hfill \end{array}$$(5)

with m_j, the mass of the SPH particle j and Δy_j, MMO the relative distance in the y component between the gas particle j and the MMO, it is the same for each component of position. The other components of the inertia tensor can be calculated in a similar way. With the inertia tensor we calculate the principal moment of inertia or eigenvalues λ₁, λ₂ and λ₃, which are listed in Table 3 for each temperature. We calculate the mean eigenvalues for a time over 10 000 yr, when approximately the rotational energy is over the gravitational energy. For temperatures of 8000 K, 5000 K and 3000 K, the eigenvalues of two components are almost equal and the value for the perpendicular component for each temperature are 1.84, 1.87 and 1.96 times higher than the other two. For a thin disk, the expected value for the perpendicular component is twice the values in the disk plane in agreement with our result for warmer temperatures. Similarly, for temperatures of 1000 K and 500 K the eigenvalues of two components are similar, but the perpendicular component are 1.16 and 1.28 times the values in the disk plane. This difference at colder temperatures could be because the energy and axisymmetrical conditions are reached at later times. Calculating the eigenvalues for later times we note that from ∼20 000 yr the mean eigenvalues for 1000 K are λ₁ = 5.60, λ₂ = 2.92 and λ₃ = 2.89, for 500 K with a time over ∼35 000 yr, the mean eigenvalues are λ₁ = 5.61, λ₂ = 3.36 and λ₃ = 3.56. This suggests that colder gas clouds need more time to form a structure like a stable thin disk at the same time due to their higher instability.

In Fig. 11, we present the angular velocity and the Toomre-Q parameter as a function of the radius for the last snapshop shown in Fig. 1. In the top panel, the rotational curve for all temperatures tends to be below the Keplerian rotational curve near the MMO, due to the sub-Keplerian correction when considering the contribution of a negative pressure gradient. At larger radii, velocities higher than the Keplerian prediction can be produced by the turbulent behavior triggered by the external torque caused by the stars. In order to quantify gravitational instabilities, we use the Toomre-Q parameter (Toomre 1964), which measures the gravitational stability against small perturbations in a rotating, self-gravitating disc. It is defined as

$$\begin{array}{c}\hfill Q=\frac{c_s\kappa }{\pi G\sum },\end{array}$$(6)

Fig. 11. Radial profiles of the rotation velocity (top panel) and Toomre-Q parameter (bottom panel) around the MMO for all temperatures in the last snapshot shown in Fig. 1.

where ∑ is the surface density of the disc, c_s is the sound speed, and κ is the epicyclic frequency. For warmer temperatures, the general trend shows a stable disk (Q > 1) around the MMO, which then moves into the unstable region at larger radii. In contrast, for simulations with colder temperatures, the disk is Toomre-unstable at all radii, implying higher fragmentation and it can occur anywhere on the disk.

Although we cannot say with certainty that a disk formation exists, the gas tends to rotate around the MMO due to angular momentum conservation, forming a elongated turbulent structure that creates a large number of sink particles that will later collide with the MMO. Most of the protostars that collide with the MMO were formed in this quasi-disk, where the average distance of the created sink particles relative to the MMO for all simulations is 0.0029 ± 0.0027 pc, with 1985.90 ± 514.011 stars with a distance shorter than 0.01 pc, corresponding to 98 ± 2% of the created stars. The combined average lifetime of the new sink particles before colliding with the MMO are 207.32 ± 1101.30 yr for a temperature of 8000 K, 284.56 ± 1434.36 yr for a temperature of 5000 K, 316.25 ± 1463.48 yr for a temperature of 3000 K, 127.80 ± 761.55 yr for a temperature of 1000 K, and finally, 47.24 ± 133.44 yr for a temperature of 500 K. Colliding the protostars faster in colder gas clouds, having less time to grow in mass.

To understand the relationship between collisions and protostar formation over time, we compare in Fig. 12 the rates of new protostars (new stars), total collisions (total coll), stars colliding with the MMO (MMO), stars that are part of other collision products (collision prod.) and escaping stars (escaped) and the number of stars in the system (N stars). Our results show that most of the protostars collide with the MMO for colder simulations and the rate of other collision products increases with the temperature. Additionally, most of these other collision products occur in the most chaotic time, between the birth of the MMO (dotted line) and the moment when it becomes the MMO in the system (dashed line), moment when it dominates collisions due to its large cross section.

Fig. 12. Time evolution of five rate categories: (1) Formation rate (blue; new stars), (2) total collision rate (orange; total coll.), (3) rate of stars that collided with the MMO (green; MMO), (4) rate of stars that are part of other collision products (red; collision prod.), and (5) rate of stars that escaped from the cluster (purple; escaped). The solid black line shows the number of stars. The vertical dotted and dashed black lines indicate the birth of the MMO and the moment at which the system is dominated by the stellar mass Mstars > Mgas.

In our study, we consider that a protostar escaped from the cluster when its distance is greater than 10 R_gas and its kinetic energy is greater than its potential energy. We observe that in general for all simulations, the number of escaping stars is very small compared to their creation rate. Additionally, the escaped stars are not concentrated in a specific period but are instead distributed throughout the entire evolution. Finally, although there are continuous creations and collisions, the total number of protostars tends to be more or less constant, with a more important decrease in protostars for the coldest simulation.

4.3. MMO for different temperatures

Effects related to the change in gas temperature can greatly affect the evolution of the MMO in a collapsing gas cloud. In this section we will show the main differences in the evolution of the MMO for different gas temperatures.

In all our simulations we reach the formation of a single massive object with a mean mass of 2.10 ± 0.20 × 10⁴ M_⊙ for the simulations with T = 8000 K, 5000 K, 3000 K and 2.70 ± 0.10 × 10⁴ M_⊙ for the simulations with T = 1000, 500 K. As we mentioned in Section 4.1.3, the high accretion rates of the central object are driven by a strong gas inflow, triggered by the core contraction. As shown in the third panel of Fig. 5, in all simulations the MMO accretion rates exceed the critical accretion rate ṁ_crit = 0.04 M_⊙ yr⁻¹ determined by the dashed black line for a long period of time. At the same time, the MMOs enter the SMS track, reaching radii of the order of ∼10⁴ R_⊙, and since we assume the most optimistic case analyzed by Reinoso et al. (2023), with t_KH, surf = 100t_KH, the SMS track is maintained until the end of the simulation, when only 0.18 ± 0.11 of the gas mass remains. Through the swelling and the increase in radius, there is an increase in the cross section of the MMO and, therefore it reaches the highest probability of collision.

In Fig. 13 we show the collision rate and radius of the MMO for a simulation of each temperature. We observe that the highest number of collisions occurs when the MMO increases and reaches the largest stellar radius, leading to a maximum collision rate of ∼0.60 ± 0.20 collisions per year.

Fig. 13. Evolution of the MMO collision rate (dashed lines) and the radius R of the MMO (solid lines) for each temperature.

In colder simulations the accretion rates of the MMO tend to be higher and remain for a longer time above the critical accretion rate. The average accretion rate during this phase, for the simulations with T = 8000 K, when it exceeds ṁ_crit is 0.56 M_⊙ yr⁻¹, and it is maintained for 21671.67 yr. In case of T = 500 K, we have 1.20 M_⊙ yr⁻¹ for 24106.43 yr, making the average accretion rate a factor 2.14 higher. The high accretion rate triggered by the falling gas leads to a faster formation of a massive object in colder environments.

The MMO evolves as a SMS due to the high accretion rates and, therefore, the stellar radius leads to an increase of the number of collisions. In our simulations, the only two mechanisms by which MMOs can obtain mass are through gas accretion and/or runaway collisions. In Fig. 14 we show the mass fraction obtained by collisions and the fraction obtained by accretion for the MMO at the end of the simulations. We clearly note a predominance of the collision mechanism as the main way of obtaining mass for the MMO in simulations with higher temperatures, contributing 75.16 ± 4.60% of the total mass of the MMO via collisions. On the other hand, the MMOs in the simulations with lower temperatures obtain their mass mainly through accretion, reaching 60.50 ± 9.55% through this mechanism.

Fig. 14. Final mass fraction of the MMOs for all simulations from 8000 K − 500 K at ∼0.3 Myr. We plot the mass acquired through accretion (blue) and the mass obtained via collisions (orange).

The difference in how the MMO obtains mass in simulations with different temperatures is because the stars that collide with the MMO have more mass in the simulations with higher temperatures, as shown in Fig. 15 where we present the combined mass distributions of the stars that collide with the MMO for each temperature. On average the stars that collide with the MMO in the cluster for a temperature of 8000 K have more time to grow in mass and radius through accretion compared to those in colder simulations, with average values of 13.9 M_⊙, 8.5 M_⊙, 7.1 M_⊙, 5.2 M_⊙, 4.7 M_⊙ for 8000 K, 5000 K, 3000 K, 1000 K, 500 K. The average lifetime of the stars that collide with the MMO for the gas cloud with 8000 K is ∼1340.37 ± 3835.66 yr, while it is ∼869.86 ± 3141.58 yr for the simulation with 500 K.

Fig. 15. Combined mass distribution of the stars that collided with the MMO in the simulations with different temperatures.

In our simulations we obtain an efficiency (M_MMO/M_Total) for the formation of the MMO of about 0.71 ± 0.08 in warmer simulations and 0.93 ± 0.04 in colder ones, which is mainly due to the increased gravitational instability as a result of the high ratio between the gas mass and the Jeans mass. As the collapse of the system is initially driven by the gravitational instability of the gas, it falls toward the center of the cluster faster and higher accretion rates are reached for the central protostars in colder simulations, directly impacting the efficiency of MMO formation, because in less unstable configurations the corresponding timescales will be longer and a fragmentation-induced starvation regime (e.g., Peters et al. 2010) could be possible. In Fig. 16 we show the efficiencies obtained in simulations with different temperatures as a function of the ratio between the total gas mass and Jeans mass, which we further compare with the results of Chon & Omukai (2020), Solar et al. (2022) and Reinoso et al. (2023). We observe a correlation between higher efficiencies and more unstable configurations, reaching efficiencies ≳80% in the simulations with colder temperatures where the ratio between gas mass and Jeans mass is higher.

Fig. 16. Efficiency of the MMO formation as a function of the gravitational instability measured via Mgas/MJeans for simulations with temperatures from 8000 K − 500 K, including the results presented by Reinoso et al. (2023) (light blue triangles), Solar et al. (2022) (gray square), and Chon & Omukai (2020) (pink diamond). The error bars were taken into account according to the 1σ error bars by Reinoso et al. (2023).

At the end of the simulations, we determine the final mass of the MMO, its radius, mean accretion rate and the total number of collisions with the MMO. The results are presented in Fig. 17. The maximum mass reached by the MMO at the end of each simulation increases for lower temperatures with very high efficiencies due to the high ratio of the total gas mass over Jeans mass, reaching close to ϵ ≈ 0.95 for simulations with 500 K after ∼30 000 yr and ϵ ≈ 0.60 for simulations with 8000 K, in agreement with the results reported by Reinoso et al. (2023), where efficiencies of ϵ = 0.90 − 0.68 have been reached by evolving the cluster for 200 kyr.

Fig. 17. Final parameters of the MMO (blue dots) as a function of the initial temperature in each simulation. The mass reached at the end of the simulation (upper left), the mean radii (upper right), the mean accretion rates (lower left), and the total number of collisions with the MMO at the end of the simulation (lower right).

Observing the radii of the MMOs in the upper right panel of Fig. 17, we notice that they increase monotonically with the mass of the MMO. The mean final radius in simulations of lower temperatures is ∼44026.07 R_⊙, a factor ∼1.25 higher compared to warmer simulations. After the initial growth of the radius, the MMOs maintain their radius on a similar level due to high accretion rates, independently of their final masses. The mean accretion rate of the MMOs decreases with increasing temperature, averaging 0.14 M_⊙ yr⁻¹ for the warmest simulations and 0.57 M_⊙ yr⁻¹ for the coldest one, which could explain the rapid evolution of the MMO in colder simulations. Finally, in the bottom right panel we show the total number of collisions for the different temperatures, with a tendency to have more collisions in configurations with lower temperatures. In colder simulations, we also have a larger variation due to random overdensities leading to disk fragmentation.

4.4. Accretion regime

Our results indicate that the supercompetitive accretion scenario can be reached even in the presence of fragmentation. This scenario corresponds to a situation in which only a few central objects grow supermassive with significant accretion inflows of ≳0.1 M_⊙ yr⁻¹, while a larger number of stars compete for the gas reservoir (Chon & Omukai 2020). The conditions where this scenario develops were studied through analytical arguments by Schleicher et al. (2023), taking into account the growth of a massive object via stellar collisions and accretion due to self-gravity and Bondi–Hoyle accretion. They included corrections for angular momentum effects due to a suppression factor (Inayoshi et al. 2018) and the velocity of the massive objects itself.

They found that at the beginning of the mass growth of the MMO, it will be regulated by self-gravity induced accretion because the density and mass of the central object are sufficiently low, regardless of the number of stars. The accretion rate due to self-gravity depends strongly on the temperature, reaching higher accretion rates for warmer gas clouds, with ∼0.2 M_⊙ yr⁻¹ for gas clouds with 8000 K, and ∼10⁻³ M_⊙ yr⁻¹ for temperatures of 300 K. Only when the central object becomes sufficiently massive the gas accretion can be dominated by Bondi–Hoyle accretion due to the reduction of the mean density in a constant spherical region. Since self-gravity induced accretion depends on temperature as $\propto T^{\frac{3}{2}}$, while Bondi-Hoyle accretion scales as $\propto T^{\frac{-3}{2}}$, colder temperatures allow Bondi–Hoyle accretion to dominate at lower masses of the MMO. Because the cooling and heating processes are complicated, however, there are uncertainties of how the temperature will evolve around the MMO in a realistic scenario.

Calculating the accretion rates of the MMO and the growth via collisions, we can compare our results with the theoretical analysis performed by Schleicher et al. (2023). In the top panels in Fig. 18, we compare the accretion rates of the MMO for three of our simulations with temperatures of 8000 K, 3000 K and 500 K, with their analytical predictions for both the accretion rate driven by self-gravity and Bondi–Hoyle accretion rate as a function of the mass of the central massive object. For the self-gravity induced accretion rate we assume an isothermal collapse of 8000 K, 3000 K and 500 K. While the Bondi–Hoyle accretion rate is formally defined at infinity (ρ_∞), we evaluate it through the average density of the cloud. Since Bondi–Hoyle accretion must be evaluated outside the Bondi radius (R_Bondi), it is only valid when R_Bondi < R_gas.

Fig. 18. Comparison of the mass growth rate (top panel) and the growth rate due to stellar collisions (bottom panel) as a function of the mass of the MMO for temperatures of 8000 K, 3000 K, and 500 K with the theoretical accretion rates of Bondi-Hoyle (red line) and accretion regulated by self-gravity (orange line) from Schleicher et al. (2023).

For a temperature of 8000 K and 3000 K, we observe that the self-gravity inflow generally provides a good approximation for the accretion rate, despite some differences, such as self-gravity underestimates the general accretion flow for 3000 K and between 200 M_⊙ − 3000 M_⊙ for 8000 K. At later times, when the mass of the MMO exceeds 10⁴ M_⊙, it becomes difficult to determine whether accretion is dominated by Bondi–Hoyle or self-gravity accretion in the 8000 K simulation. In the case of 3000 K, since R_Bondi is comparable to R_gas, the configuration breaks down the Bondi-Hoyle approximation during the growth of the MMO. In masses larger than ∼7000 M_⊙, self-gravity effects are more relevant. For the simulation of 500 K, the Bondi-Hoyle accretion can be approximated until masses of ∼700 M_⊙, for the central massive object. For larger masses, the accretion exceeds the Bondi accretion due to it is boosted by the self-gravitational effects of the gas, triggered by the high instability and increased density in the center of the cloud. As a result, the self-gravity drives the gas toward the center of the system, in addition to the gravity of the central massive object. The variability in the accretion rates from our result is also due to fluctuations and substructure within the gas.

In the bottom panels in Fig. 18, we compare the growth rate due to stellar collisions with our simulations at three different temperatures inside a region of 0.1 pc, where most of the fragments are distributed around the MMO. To calculate the theoretical contribution from stellar collisions, we used the approximation from Schleicher et al. (2023), where they consider the average mass of the stars over the timescale for collisions from Leigh et al. (2017). Under the conditions of supercompetitive accretion, it is possible to assume $m_{\ast }\sim \overline{m}$, R_* ≪ R_MMO and R_cl ∼ R, leading to the following equation to estimate the contribution from collisions:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{coll}}\sim \frac{12\sqrt{2}{R}_{\mathrm{MMO}}^2{m}_{\ast }^{3/2}N}{G^3{M}^{1/2}{f}_{\ast }^3}{\left(\frac{G}{2R}\right)}^{7/2}\sqrt{M_{\mathrm{tot}}},\end{array}$$(7)

where f_* = M_cl/M_tot, M_cl is the mass of the cluster, R_MMO is the radius of the MMO, m_* is the average stellar mass, N is the number of stars, M is the mass of the MMO, and R the radius of the cluster. To compare the theoretical estimation to our data, we assume an average mass of the fragments that collide with the MMO of m_* = 15.0 M_⊙, 7.4 M_⊙ and 4.8 M_⊙ for gas clouds with 8000 K, 3000 K and 500 K. For our simulations, the typical number of stars is N ≈ 145, with a combined mean mass of m_mean = 5.14 M_⊙, resulting in a ratio f_* = 0.024. In general, we find good agreement between the analytical estimate from Schleicher et al. (2023) and the mass growth of the MMO due to stellar collisions of our results. As the analytical equation predicts, the contribution from collisions increases with the mass of the MMO, because R_MMO increases with the mass following an SMS stage, leading to a larger cross section. Additionally, the collision contribution is sensitive to the mass of the fragments that collide with the MMO, being more massive for higher temperatures and lower for colder temperatures due to fragmentation effects, where more stars merge with the central massive object but with lower masses. We note that the analytic rate has a better match for the warmest simulations and underestimates the mass growth for colder temperatures when the MMO has a mass lower than 10⁴ M_⊙. The formation of a structure like a disk around the MMO (see Section 4.2) will produce uncertainty due to the generation of substructures that later will collapse with the MMO, possibly increasing the collision rates. This could be more relevant for colder temperatures where the gas will produce more fragments and higher star formation rates can be reached. The differences between the analytic equation and our results probably are due to the simplifying assumptions such as an uniform star distribution, a constant number of stars and a constant mass of the fragments within the cluster. Even with the uncertainties associated with the simplifications for the collision rate, we conclude that in general our results are in agreement with the theoretical equations and provide a good approximation to predict the mass of the central massive object obtained via collisions.

5. Discussion and conclusions

When a gas cloud collapses, fragmentation can occur even in the presence of a minimum amount of metals. This leads to the formation of a protostellar cluster that is embedded in a massive gas cloud. Different environmental factors such as the inhibition of H₂ formation through ionizing photons from the first generation of stars or a first quasar as a source of UV radiation can influence the evolution of the gas temperature and its collapse behavior. Additionally, metal contamination from halos hosting this first generation of stars can significantly affect its cooling properties (Schneider et al. 2006, 2003; Omukai et al. 2005). Omukai (2001) studied the effect of different intensities of far-ultraviolet (FUV) radiation on the primordial star formation and reported that under intense FUV radiation, gas collapse is approximately isothermal and regulated by atomic cooling lines. In contrast, a cloud irradiated with weaker FUV radiation collapses dynamically via H₂ cooling and reaches a lower minimum Jeans mass. This leads to increased fragmentation. Omukai et al. (2008) studied the impact of the gas metallicity on the temperature evolution, where gas fragmentation is inevitable above a critical metallicity between Z_crit ≈ 5 × 10⁻⁶ Z_⊙ − 3 × 10⁻⁴ Z_⊙, and dense clusters of low-mass stars can form. These might still form a remnant intermediate-mass BH through stellar collisions.

We explored the formation of a massive object during the collapse of a protostellar cluster that was embedded in a massive gas cloud following fragmentation, accretion, and collisions on subparsec scales. We investigated the impact of different temperatures on the formation of a massive object and its potential to form a SMBH seed through a simplified equation of state, in which the gas evolved isothermically according to Eq. (1). We performed a series of hydrodynamical and N-body simulations using the AMUSE framework, including the routines developed by Reinoso et al. (2023) for the creation of new sink particles, gas accretion, pseudo- (proto-) stellar evolution and stellar collisions. In our simulations, the gas collapsed nearly spherically as a result of the high gravitational instability before it forms an elongated gas cloud around the MMO. Even with this elongated structure, 90% of the particle formation on subparsec scales occurs close to the MMO and does not prevent the formation of a SMS with ∼10⁴ M_⊙.

The formation of this elongated structure was triggered by the conservation of angular momentum, which led to the formation of a disk and, subsequently, to fragmentation of the gas. Effects of disk fragmentation in massive star formation were explored in previous work and under different conditions (Oliva & Kuiper 2020; Girichidis et al. 2012; Peters et al. 2010). The disk was formed from low initial angular momentum in a collapsing gas cloud that grew and formed a disk with spiral arms that fragmented in eccentric and chaotic orbits. Oliva & Kuiper (2020) explored the convergence of their results for different spatial resolutions and reported a minimum resolution for a convergence of the physical properties. Effects such as magnetic fields and supersonic and subsonic turbulence can also affect the homogeneity in the density profile and produce more fragmentation or reduce it. While the fragments are formed, they tend to migrate to the MMO, and most of the sinks in the disk finally merge with the MMO (Latif & Schleicher 2015; Sakurai et al. 2016; Chon & Hosokawa 2019). More work is required to understand whether disk fragmentation is a possible mechanism in the formation of a SMBH seed.

Even with fragmentation, the MMO grows efficiently in mass and reaches masses of ∼10⁴ M_⊙, which agrees with the results of Chon & Omukai (2020), who investigated the impact of metal enrichment in massive gas clouds from a low-metallicity regime of Z/Z_⊙ = 10⁻⁶ to a higher metallicity with Z/Z_⊙ = 10⁻³. They found that for very low and slight metal enrichment, the cloud collapsed and formed a massive object. Even when the cloud presented some fragmentation, most of it either merged with the central star or was ejected from the system. Only when the metallicity was higher than 10⁻³ Z_⊙ did the central star not become supermassive because the accretion rates were lower. We did not consider the chemical evolution, but employed a highly simplified treatment to investigate the effects of different thermodynamic conditions by varying the gas temperature. We found that in the context of highly Jeans-unstable clouds, a massive object with a mass ∼10⁴ M_⊙ can form at different initial temperatures of the gas cloud. Cooling and heating might have an important impact on the evolution of the MMO by decreasing the accretion rate and limiting the growth in mass of the central star, as was shown in the higher-metallicity case by Chon & Omukai (2020). For more stable gas clouds, the cooling could lead to the formation of star clusters and might greatly inhibit the accretion of a massive star. We focused on the regime of strong gravitational instability and explored its effects on the formation of a central massive object. We found a correlation between the efficiency in the formation of the MMO and the Jeans instability of the cloud.

Even for higher protostar formation rates, our results agree with the supercompetitive accretion regime proposed by Chon & Omukai (2020). We note that most of the created protostars collide with the MMO because they were formed by gravitational instabilities in the quasi-disk structure around the MMO. In our simulations, the protostars collide very quickly with the MMO and do not have enough time to increase in mass and reach the theoretical condition proposed by Schleicher et al. (2023) that could lead to fragmentation-induced starvation (Peters et al. 2010) and restrict the accretion of the MMO. The protostars in the cluster have t_acc ≫ t_ff on average, so that the gas collapses faster than it can be accreted by the protostars.

In the supercompetitive accretion scenario, the accretion inflow could be regulated by two mechanisms, self-gravity and Bondi-Hoyle accretion. Our results agree with previous estimates from Schleicher et al. (2023), who reported higher accretion rates in earlier states of the MMO and provided a good approximation for accretion rates in gas clouds with higher temperatures. The MMO initially grows rapidly in mass due to high accretion rates driven by self-gravity. At later times, Bondi-Hoyle accretion is boosted by the increasing gas density in the center of the cloud, driven by the self-gravity of the gas. We conclude that in an atomic cooling halo, the theoretical equation provides a good approximation of the accretion rate, which is mainly gas accretion driven by self-gravity.

Additionally, our results showed a lower contribution through collisions for colder temperatures, according to Schleicher et al. (2023), because at these temperatures the Jeans mass is lower. This results in less massive fragments, and thus, in a smaller contribution to the MMO. Furthermore, at higher temperatures, the average mass of the fragments that collide with the MMO is higher, which leads to a stronger impact on the growth via collisions. The variability in the growth rate via collisions is due to the large difference in the mass of the protostars that formed in the disk around the MMO.

In the simulations we presented, we did not include mass loss due to collisions. This effect was explored by different authors in different contexts, such as local star clusters (Sills & Bailyn 1999; Sills et al. 2000), and the mass ratios of the collision stars (Lombardi et al. 2002) and the effect of stellar structure were considered (Glebbeek et al. 2008). Alister Seguel et al. (2020) explored the impact of mass loss due to collisions and included both constant mass-loss fractions and analytical models (Lombardi et al. 2002). They found that the mass loss can strongly impact the final mass of the SMBH seeds, which in some cases lost between 60%−80% of the final mass and 20%−40% for the similar context we studied. Nevertheless masses of about 10⁴ M_⊙ were obtained, and this is sufficient for them to be viable candidates for SMBH seeds. These effects might be important for warmer temperatures as well, where the contribution to the final mass of the MMO is ≳70% and ≳25% for colder simulations.

We included a recipe for protostellar evolution to evolve the radius with the mass and accretion rate following Reinoso et al. (2023). This recipe allowed the protostar to expand and contract to the main sequence according to the accretion rate, and this also affected ejections and stellar collisions, which are sensitive to their cross sections. The object with the highest collision probability is the MMO because it has the highest accretion rates, which increases its stellar radius and thereby also its cross section. Due to the high inflow of gas and the growth in radii of the MMO, most of the protostars merge with the MMO, and only some are ejected due to three-body interactions (∼0.007% of the total mass for the most extreme case). The stellar contraction to the main sequence triggers UV radiation from the stars, which ionizes the surrounding gas and creates an ionization front. This might affect the accretion rate of the MMO by the surrounding stars that enter the main sequence. Different studies such as Schleicher et al. (2022), Reinoso et al. (2023) concluded that the stellar feedback from surrounding stars that enter the main sequence does not appreciably impact the final mass of the MMO. These authors considered a very metal-poor gas environment that collapses approximately isothermally, however, and they did not resolve the chemical evolution of the gas in detail. More work is required to reduce the uncertainty in these estimates and to determine the final masses of the objects that form under these conditions.

In all our simulations, a central massive object with a mass higher than ∼10⁴ M_⊙ formed through runaway collisions and high accretion rates. The MMO grew faster in mass and became more massive for colder temperatures through more unstable initial conditions even with vigorous fragmentation. This is similar to previous studies that modeled low-metallicity conditions (Chon & Omukai 2020; Solar et al. 2022; Reinoso et al. 2023). For all our simulations, efficiencies higher than 60% were reached and exceeded 90% for colder temperatures due to the strong gravitational instability. We observed a trend for the formation of a massive central object even with some fragmentation when the gas cloud was sufficiently unstable. The mechanism by which the MMO grows in mass is primarily driven by collisions at higher temperatures, while accretion becomes the dominant process in simulations with lower temperatures. This is a result of the combination of the timescale for the formation of the quasi-disk structure, which is the primary mechanism for the creation of new sinks in our simulations, and the lower Jeans mass at lower temperatures. The typical final outcome in our simulations is a central SMS with a small cluster of stars with masses of about 0.1 − 100 M_⊙ for warmer simulations and 0.1 − 1 M_⊙ for colder temperatures.

Acknowledgments

PS acknowledges support through ANID/Doctorado en el Extranjero convocatoria 2022 (funding number 72220198), the Federal Ministry of Education and Research and the state governments for supporting this project as part of the joint funding of National High Performance Computing (NHR) and the Kultrun cluster hosted at the Departamento de Astronomía, Universidad de Concepción. We gratefully acknowledge support by the ANID BASAL project FB21003, as well as via Fondecyt Regular (project code 1201280) and ANID QUIMAL220002. BR acknowledges funding through ANID (CONICYT-PFCHA/Doctorado acuerdo bilateral DAAD/62180013), DAAD (funding program number 57451854), and the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). BR acknowledges the support by the European Research Council via ERC Consolidator grant KETJU (no. 818930). DRGS thanks for funding via the Alexander von Humboldt – Foundation, Bonn, Germany. RSK acknowledges financial support from the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). RSK also thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. RB acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306.

References

Appendix A: Edge-on view

In this appendix, we present the density projection of the edge-on view, similar to Fig. 1.

Fig. A.1. Density projections of the edge-on view of a primordial gas cloud with an embedded protostar cluster for simulations with an initial temperature of 500 K − 8000 K at different times. The white dots are the protostars, black dots are the formed stars and the white star is the MMO. The clusters evolve in time from left to right and from the warmer to colder simulations from the top to the bottom panels.

All Tables

All Figures

Fig. 1. Density projections of the face-on view of a primordial gas cloud with an embedded protostar cluster for simulations with an initial temperature of 500 K − 8000 K at different times. The white dots are the protostars, black dots are the formed stars and the white star is the MMO. The clusters evolve in time from left to right and from the warmer to colder simulations from the top to the bottom panels.

In the text

	Fig. 2. Radial profiles relative to the MMO for the gas density, stellar mass density without the MMO, radial velocity, and accretion at approximately 15 000 yr, corresponding to the third snapshot of each temperature shown in Fig. 1.
In the text

	Fig. 3. Time evolution of the mass fraction for the simulations with an initial temperature of 8000 K, 5000 K, 3000 K, 1000 K, and 500 K, corresponding to the gas mass (blue zone), protostars in the cluster (green zone) and the escaped protostellar mass (orange zone). The vertical dotted black lines indicate the times when the protostars accreted 25%, 50%, and 75% of the gas mass.
In the text

	Fig. 4. Time evolution of the 90%, 75%, 50%, 25%, and 10% Lagrangian radii for our simulations with temperatures of 8000 K, 3000 K, and 500 K. The panels show the Lagrangian radii considering the total mass of the systems (gas and stars) relative to their centers of mass.
In the text

Fig. 5. Evolution of the gas and protostellar mass in the simulations with different temperatures (top panel). Comparison of the time evolution of the MMO mass and the total mass of the protostars without the MMO (second panel), accretion rate of the MMO with the mean accretion rate of the protostars without the MMO (third panel), where the black dashed line is the critical accretion rate ṁcrit = 0.04 M⊙ yr−1 and marks the point where the protostars pass to the SMS Track. We further show the radius of the MMO with the mean radii of the protostars (fourth panel) and the evolution of the total collision rates over time (bottom panel).

In the text

	Fig. 6. Combined mass distribution at the end of our simulations for each temperature.
In the text

	Fig. 7. Mean accretion time over the free-fall time for each star in the cluster at the end of a simulation with each temperature.
In the text

	Fig. 8. MMO mass over the total stellar mass without the MMO for each temperature. The horizontal dashed black line shows when Mstars/MMMO = 1.
In the text

	Fig. 9. Combined mass distributions of the escaped stars for simulations with all temperatures.
In the text

	Fig. 10. Ratio of the rotational kinetic energy and potential energy of the gas in a region of 0.01 pc around the MMO for each temperature (top panel). The angles ϕ and θ of the total angular momentum are plotted in spherical coordinates (middle and bottom panels).
In the text

	Fig. 11. Radial profiles of the rotation velocity (top panel) and Toomre-Q parameter (bottom panel) around the MMO for all temperatures in the last snapshot shown in Fig. 1.
In the text

Fig. 12. Time evolution of five rate categories: (1) Formation rate (blue; new stars), (2) total collision rate (orange; total coll.), (3) rate of stars that collided with the MMO (green; MMO), (4) rate of stars that are part of other collision products (red; collision prod.), and (5) rate of stars that escaped from the cluster (purple; escaped). The solid black line shows the number of stars. The vertical dotted and dashed black lines indicate the birth of the MMO and the moment at which the system is dominated by the stellar mass Mstars > Mgas.

In the text

	Fig. 13. Evolution of the MMO collision rate (dashed lines) and the radius R of the MMO (solid lines) for each temperature.
In the text

	Fig. 14. Final mass fraction of the MMOs for all simulations from 8000 K − 500 K at ∼0.3 Myr. We plot the mass acquired through accretion (blue) and the mass obtained via collisions (orange).
In the text

	Fig. 15. Combined mass distribution of the stars that collided with the MMO in the simulations with different temperatures.
In the text

Fig. 16. Efficiency of the MMO formation as a function of the gravitational instability measured via Mgas/MJeans for simulations with temperatures from 8000 K − 500 K, including the results presented by Reinoso et al. (2023) (light blue triangles), Solar et al. (2022) (gray square), and Chon & Omukai (2020) (pink diamond). The error bars were taken into account according to the 1σ error bars by Reinoso et al. (2023).

In the text

	Fig. 17. Final parameters of the MMO (blue dots) as a function of the initial temperature in each simulation. The mass reached at the end of the simulation (upper left), the mean radii (upper right), the mean accretion rates (lower left), and the total number of collisions with the MMO at the end of the simulation (lower right).
In the text

	Fig. 18. Comparison of the mass growth rate (top panel) and the growth rate due to stellar collisions (bottom panel) as a function of the mass of the MMO for temperatures of 8000 K, 3000 K, and 500 K with the theoretical accretion rates of Bondi-Hoyle (red line) and accretion regulated by self-gravity (orange line) from Schleicher et al. (2023).
In the text

Fig. A.1. Density projections of the edge-on view of a primordial gas cloud with an embedded protostar cluster for simulations with an initial temperature of 500 K − 8000 K at different times. The white dots are the protostars, black dots are the formed stars and the white star is the MMO. The clusters evolve in time from left to right and from the warmer to colder simulations from the top to the bottom panels.

In the text