Supermassive stars SMSs and heavy seed black holes as their remnants are promising candidates for supermassive black hole SMBH progenitors especially for ones observed in the early universe inlineformulatexmath idTM0001 notationLaTeXzsimeq 8.510texmathinlineformula by recent JWST observations. Expected cradles of SMSs are the atomic cooling haloes inlineformulatexmath idTM0002 notationLaTeXMrm halosimeq 107 mathrmModot texmathinlineformula where cold accretion emerges and possibly forms SMSs. We perform a suit of cosmological radiation hydrodynamics simulations and investigate star formation after the emergence of cold accretion solving radiative feedback from stars inside the halo. We follow the mass growth of the protostars for inlineformulatexmath idTM0003 notationLaTeXsim 3 mathrmMyrtexmathinlineformula resolving the gas inflow down to inlineformulatexmath idTM0004 notationLaTeXsim 0.1 mathrmpctexmathinlineformula scales. We discover that after cold accretion emerges multiple SMSs of inlineformulatexmath idTM0005 notationLaTeXmstar gtrsim 105 mathrmModot texmathinlineformula form at the halo centre with the accretion rates maintained at inlineformulatexmath idTM0006 notationLaTeXdotmstar simeq 0.04 mathrmModot mathrmyr1texmathinlineformula for inlineformulatexmath idTM0007 notationLaTeXlesssim 3 mathrmMyrtexmathinlineformula. Cold accretion supplies gas at a rate of inlineformulatexmath idTM0008 notationLaTeXdotMrm gasgtrsim 0.010.1 mathrmModot mathrmyr1texmathinlineformula from outside the halo virial radius to the central gas disc. Gravitational torques from spiral arms transport gas further inwards which feeds the SMSs. Radiative feedback from stars suppresses inlineformulatexmath idTM0009 notationLaTeXmathrmH2texmathinlineformula cooling and disc fragmentation while photoevaporation is prevented by a dense envelope which attenuates ionizing radiation. Our results suggest that cold accretion can bring efficient BH mass growth after seed formation in the later universe. Moreover cold accretion and gas migration inside the central disc increase the mass concentration and provide a promising formation site for the extremely compact stellar clusters observed by JWST.
['Astrophysics - Astrophysics of Galaxies']
Sequential formation of supermassive stars and heavy seed BHs through the interplay of cosmological cold accretion and stellar radiative feedback
{'Sequential formation of supermassive stars and heavy seed BHs through the interplay of cosmological cold accretion and stellar radiative feedback': 'Masaki Kiyuna, 1 ★ Takashi Hosokawa, 1 and Sunmyon Chon 2  \n1 Department of Physics, Graduate School of Science, Kyoto University, Sakyo, Kyoto 606-8502, Japan 2 Max-Planck-Institut f ¥ 𝑢 r Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany  \nAccepted 2024 October 15. Received 2024 October 15; in original form 2024 August 29', 'ABSTRACT': 'Supermassive stars (SMSs) and heavy seed black holes, as their remnants, are promising candidates for Supermassive Black Hole (SMBH) progenitors, especially for ones observed in the early universe 𝑧 ≃ 8 . 5 -10 by recent JWST observations. Expected cradles of SMSs are the atomic cooling halos ( 𝑀 halo ≃ 10 7 M ⊙ ), where "cold accretion" emerges and possibly forms SMSs. We perform a suit of cosmological radiation hydrodynamics simulations and investigate star formation after the emergence of cold accretion, solving radiative feedback from stars inside the halo. We follow the mass growth of the protostars for ∼ 3 Myr, resolving the gas inflow down to ∼ 0 . 1 pc scales. We discover that, after cold accretion emerges, multiple SMSs of 𝑚 ★ ≳ 10 5 M ⊙ form at the halo centre with the accretion rates maintained at / 𝑚 ★ ≃ 0 . 04 M ⊙ yr -1 for ≲ 3 Myr. Cold accretion supplies gas at a rate of / 𝑀 gas ≳ 0 . 01 -0 . 1 M ⊙ yr -1 from outside the halo virial radius to the central gas disc. Gravitational torques from spiral arms transport gas further inward, which feeds the SMSs. Radiative feedback from stars suppresses H2 cooling and disc fragmentation, while photoevaporation is prevented by a dense envelope, which attenuates ionising radiation. Our results suggest that cold accretion can bring efficient BH mass growth after seed formation in the later universe. Moreover, cold accretion and gas migration inside the central disc increase the mass concentration and provide a promising formation site for the extremely compact stellar clusters observed by JWST.  \nKey words: quasars: supermassive black holes - stars: Population III - galaxies: formation.', '1 INTRODUCTION': 'Recent studies reveal the presence of over 200 active galactic nuclei (AGNs) at redshifts 𝑧 = 6 -10, suggesting that supermassive black holes (SMBHs) with masses ranging from 10 6 to 10 10 M ⊙ already formed at the time when the universe was a few × 0 . 1 Gyr years old (e.g. Mortlock et al. 2011; Bañados et al. 2018; Inayoshi et al. 2020; Wang et al. 2021; Volonteri et al. 2021). The formation of these SMBHs is one of the greatest mysteries in modern astrophysics, as it appears to be highly challenging for them to form in such a short duration. Remarkable discoveries of high-z SMBHs at 𝑧 ≃ 8 . 5 -10 by deep observations with the JWST (Kokorev et al. 2023), combined with Chandra (Bogdán et al. 2024; Kovács et al. 2024), make the situation more serious.  \nOne possible formation channel for those high-z SMBHs is mass growth from the remnant BHs of Population (Pop) III stars with 𝑚 · ∼ 10 1 -3 M ⊙ formed in mini-halos with 𝑀 halo ∼ 10 5 -6 M ⊙ (light seed model; Bromm et al. 1999, 2002; Abel et al. 2002; Yoshida et al. 2003, 2006; Hosokawa et al. 2011; Hirano et al. 2014; Susa et al. 2014; Sugimura et al. 2020, 2023). If the mass accretion rate onto such BHs is maintained continuously at the Eddington rate, this scenario cannot explain the SMBHs observed at 𝑧 ≳ 8 . 5. Furthermore, the mass accretion rate decreases below Eddington rate  \ndue to radiation feedback (e.g. Alvarez et al. 2009; Milosavljević et al. 2009; Park & Ricotti 2011, 2012, 2013; Aykutalp et al. 2014; Smith et al. 2018).  \nAn alternative pathway for the SMBH formation is the heavy seed model, where supermassive stars (SMSs) collapse and generate seed BHs with 𝑚 · ∼ 10 4 -6 M ⊙ (Bromm & Loeb 2003). They potentially emerge out of the pristine gas cloud within atomic cooling halos (ACHs) with masses of 𝑀 halo ∼ 10 7 -8 M ⊙ . SMS formation begins with the collapse of a cloud triggered by effective Ly𝛼 cooling, when H 2 cooling is inhibited. In this scenario, extremely rapid accretion onto a protostar at rates / 𝑀 ∼ 0 . 1 M ⊙ yr -1 facilitates SMS formation within the stellar lifetime, ∼ Myr(Latif et al. 2013; Chon et al. 2018). The SMS formation in ACHs benefits the subsequent mass growth of the heavy seed BHs. In ACHs, the gas dissipates its internal energy and concentrates to the halo centre due to Ly𝛼 cooling. Given that these seed BHs are formed at the very centre of the halo, they can grow efficiently in dense gas. Cold accretion will accumulate the gas to the halo centre and help efficient growth (Di Matteo et al. 2012).  \nDuring standard Pop III star formation, H 2 molecules act as effective coolants, lowering the gas temperature to around 𝑇 ≃ 200 K and thus hindering SMS formation. Various authors have proposed different physical mechanisms to inhibit H 2 cooling. These include intense Lyman-Werner (LW; 11 . 2 eV ≤ ℎ𝜈 ≤ 13 . 6 eV) background radiation originating from nearby galaxies (Omukai 2001; Dĳkstra et al. 2008; Shang et al. 2010; Sugimura et al. 2014; Chon et al.  \nFigure 1. A schematic picture showing the evolution outline followed by our cosmological radiation hydrodynamics simulations.  \n<!-- image -->  \n2016), turbulence caused by halo mergers (Wise et al. 2019; Regan et al. 2020; Latif et al. 2022; Toyouchi et al. 2023), and the relative streaming motion between baryons and dark matter (Hirano et al. 2015; Schauer et al. 2017).  \nInayoshi & Omukai (2012) proposed a model in which SMSs can form in a high-temperature cloud, which is shock-heated by cold accretion. When a gas cloud is compressed by shocks, the gas becomes both dense ( 𝑛 ≳ 10 4 cm -3 ) and hot ( 𝑇 ≳ 8000 K). In such case, H 2 is collisionally dissociated and the cloud collapses without H 2 cooling. Cold accretion is a phenomenon found in studies of mature galaxy formation ( 𝑀 halo ∼ 10 10 -12 M ⊙ ) at lower redshifts (Birnboim & Dekel 2003; Kereš et al. 2005; Dekel & Birnboim 2006; Dekel et al. 2009), which involves the gas stream reaching the halo centre before experiencing virialisation shocks due to efficient radiative cooling. Dense and strong shocks are generated only near the halo centre with post-shock temperatures around 𝑇 ∼ 𝑇 vir ≃ 10 4 K, which is potentially available for SMS formation. Inayoshi & Omukai (2012) investigated the thermal evolution when shock heating occurs in the ACH and the possibility of the formation of SMSs by the one-zone model, including shock heating caused by the cold accretion in ACHs (Wise & Abel 2008; Greif et al. 2008). However, Fernandez et al. (2014) performed cosmological simulations and showed that cold accretion does not occur, at least until their simulation end with the first cloud collapse in ACHs.  \nIn our previous work Kiyuna et al. (2023) (hereafter K23), we have studied the emergence of the cold accretion in ACHs and whether it causes SMS formation. We have found that cold accretion does emerge in ACHs when their mass exceeds 𝑀 halo ∼ 10 7 M ⊙ , while normal Pop III stars form slightly before its emergence. Cold accretion transports filamentary gas flows into the halo centre, forming a dense rotation-supported disc. We have shown the presence of hot and dense gas within the disc, where collisional dissociation of H 2 molecules occurs efficiently. Instead of following SMS formation by simulations, we post-processed the snapshot to estimate the maximum amount of the hot and dense gas potentially available for SMS formation. This limitation arose because we did not include the radiation emitted by Pop III stars within the same halo, which prevented us from quantitatively evaluating the amount of hot and dense gas. In reality, these stars serve as sources of radiative feedback, potentially photodissociating H 2 molecules in the disc to promote the formation of SMSs. Concurrently, these stars could also ionise the gas, leading to disc photoevaporation. We have to update our simulations, considering such internal stellar radiative feedback.  \nIn this study, we perform cosmological radiation hydrodynamics simulations to explore the onset of cold accretion and subsequent star formation influenced by radiative feedback from stars within the same  \nhalo. Fig. 1 illustrates the evolutionary outline we follow. Previous studies, including K23, have shown that with a moderate background LW radiation, normal Population III star formation begins when 𝑇 vir ≃ 10 4 K ( 𝑀 halo ∼ 10 7 M ⊙ ). We refer to these Population III stars as the "preceding star", as its formation precedes the emergence of the cold accretion. The onset of the cold accretion is delayed since the radiation by the preceding star evacuates the halo gas. Once the cold accretion channels a substantial amount of gas into the halo centre, it triggers further star formation. These "subsequent stars" also cause radiative feedback, which affects the thermal and chemical conditions of the gas in the dense gas disc. We study the entire evolution, specifically to investigate whether SMSs emerge among the subsequent stars following the advent of cold accretion.  \nThe organisation of this paper is as follows. Section 2 outlines our simulation methodology. Our simulation results are presented in Section 3. We focus on the development of cold accretion with radiative feedback in Section 3.2. Section 3.3 discusses the emergence of SMSsfollowing the onset of cold accretion. In Section 4, we examine the combined effects of cold accretion and radiative feedback on star formation. Discussions and summaries are provided in Sections 5 and 6, respectively.', '2 METHODS': 'We perform a series of cosmological N-body + Smoothed Particle Hydrodynamics (SPH) simulations using the code Gadget-3 (Springel 2005), as in K23 (see Sections 2.1 and 2.2). We improve the spatial resolution to follow the star formation more in detail (Section 2.3). We also consider stellar radiative feedback, solving the transfer of radiation emitted by massive protostars (Section 2.4). The simulation extends for ≳ 3 Myr after the onset of cold accretion, tracking the long-term evolution of star formation considering the stellar radiative feedback.  \nThroughout the paper, we adopt the cosmological parameters of PLANCK13 with the following values: Ω m = 0 . 3086, Ω Λ = 0 . 6914, Ω b = 0 . 045, ℎ = 0 . 6777, 𝜎 8 = 0 . 8288, 𝑛 spec = 0 . 9611 (Planck Collaboration et al. 2014).', '2.1 N-body + zoom-in SPH simulation': 'Wegenerate cosmological initial conditions at 𝑧 = 99 with a volume of 𝑉 = ( 1 ℎ -1 cMpc ) 3 using the code MUSIC (Hahn & Abel 2013). We use the same random seed as "halo-A" in K23, where the cold accretion emerges at the earliest epoch among the three halo samples, 𝑧 ≃ 18.  \nWe first perform a dark matter (DM) only N-body simulation with 256 3 particles. The mass of each DM particle is 𝑚 DM = 4360 ℎ -1 M ⊙ . We identify the position of the most massive halo in the simulation box. We set a zoom-in region around the specified halo with a volume of ( 0 . 4 ℎ -1 cMpc ) 3 and regenerate the initial conditions at 𝑧 = 99 consistent with the random field at the base resolution. Inside the zoom-in region, we effectively place ( 1024 ) 3 gas and DM particles. The mass resolutions inside the zoom-in region are 𝑚 DM = 68 ℎ -1 M ⊙ for DM and 𝑚 gas = 11 . 6 ℎ -1 M ⊙ for gas.  \nWe perform a N-body + SPH + Ray-tracing simulation from this initial condition. The calculation proceeds in the same way as in K23 up to the epoch of the normal Pop III star formation before the onset of cold accretion (the formation of "preceding star" in Fig. 1). Following this phase, we continue to examine star formation, incorporating radiative feedback from stars.', '2.2 Chemical network and thermal processes': 'Weassumezero metallicity (Z = 0) for all the gas particles and solve a non-equilibrium chemistry network with H , H + , e -, H -, and H 2 using an implicit scheme. We follow the non-equilibrium evolution of the thermal energy considering radiative processes by Ly𝛼 emission, ro-vibrational transition of H 2 , and continuum processes related to H atom following Matsukoba et al. (2021).  \nWe include the LW background radiation to suppress the star formation in the mini-halos with 𝑀 halo ≲ 10 6 M ⊙ (Haiman et al. 2000). We set the black-body spectrum with effective temperature 𝑇 eff = 10 4 Kandtheintensity 𝐽 21 = 10, which is typical at 𝑧 ≃ 20 -10 (Dĳkstra et al. 2008; Holzbauer & Furlanetto 2012), where 𝐽 21 is the specific intensity at the LW band ( ℎ𝜈 = 11 . 2 -13 . 6 eV ) normalised by 10 -21 erg cm -2 s -1 Hz -1 str -1 . This value is much lower than the critical intensity for which H 2 cooling is completely disabled, i.e. 𝐽 21 , crit ≃ 10 2 for 𝑇 eff = 10 4 K (Omukai 2001; Shang et al. 2010; Sugimura et al. 2014). We consider the self-shielding effect against LW background (Wolcott-Green et al. 2011), measuring the column density as 𝑁 H2 = 𝑛 H 𝑦 H2 𝜆 Jeans , where 𝑦 H2 is H 2 fraction and 𝜆 Jeans is the local Jeans length. We note that this effect is neglected in K23, where ignoring radiative feedback makes it difficult to assess the strength of the LW field. We consider the photoionisation of H, photodissociation of H 2 , and photo-detachment of H -by the radiation emitted by stars forming within the halo, which will be described in Section 2.4.', '2.3 Sink prescription and particle splitting': "We introduce a sink particle once the density exceeds 𝑛 H, th , largely following the method described in K23. We differently set the threshold density for the 'preceding star' and 'subsequent stars', which form before and after the onset of cold accretion (Fig. 1).  \nWe set 𝑛 H , th = 2 × 10 6 cm -3 for the preceding star, which forms at 𝑧 = 18 . 9 before the emergence of cold accretion (see Fig. 1 and Section 3.2). We set the radius of the sink particles to be three times the smoothing length of the SPH particle at the sink formation. The resulting sink radius is 𝑟 sink ≃ 1 pc for the preceding star. Note that this sink radius is not small enough to follow the formation of the individual Pop III stars. For this reason, we use a particular treatment to model the radiative feedback from the preceding star (see Section 2.4.1).  \nTobetter resolve the subsequent star formation, we perform particle splitting just after the preceding star formation for the particles within a radius of 40 ℎ -1 ckpc ≃ 3 kpc around it. Following Kitsionas &  \nWhitworth (2002), each gas particle is split into 13 daughter particles (Chon et al. 2021, for detail). The mass of the split gas particle is 𝑚 gas = 0 . 89 ℎ -1 M ⊙ . With the improved resolution, we set a higher threshold density to create sink particles 𝑛 H , th = 2 × 10 8 cm -3 and a smaller sink radius 0 . 1 pc. We also shrink the sink radius of the preceding star particle from 1 pc to 0 . 1 pc, according to the resolution enhancement. This length corresponds to the Bondi radius of a star with mass 𝑚 ★ = 10 3 M ⊙ in the ambient gas with a temperature of 𝑇 = 10 4 K. This indicates that we can safely resolve gas accretion onto the sink particles for stars with masses larger than ∼ 10 3 M ⊙ .  \nWe allow the merger of sink particles once the separation of a pair of sink particles becomes smaller than the sum of sink radii. To suppress spurious fragmentation before the sink formation, we turn off the cooling at the density 𝑛 H > 𝑛 H , ad ≡ 0 . 5 × 𝑛 H, th (Chon et al. 2018; Susa 2019).  \nWe can only resolve the gravitational collapse of the gas with 𝑇 ≳ 10 3 K, where the Jeans mass at the sink formation is resolved by more than 80 gas particles (Bate & Burkert 1997). In our simulation, a small portion of the gas at 𝑛 H ∼ 10 8 cm -3 is efficiently cooled via H 2 cooling below 𝑇 ∼ 10 3 K. That cold gas component can form Pop III stars with 𝑚 ★ ≲ 10 2 M ⊙ , which will be unresolved in our simulation. Since our goal is to study the formation of SMSs, we follow the formation of Pop III stars with 𝑚 ★ ≳ 10 3 M ⊙ including SMSs, neglecting the formation of such less massive stars.", '2.4 Radiative feedback models and radiation transfer': 'We incorporate radiative feedback from protostars, assuming their spectra to be black-body with luminosity 𝐿 and effective temperature 𝑇 eff . We differently treat stars forming before and after the onset of the cold accretion (i.e., the preceding and subsequent stars as referred to in Fig. 1).', '2.4.1 Before the emergence of cold accretion': 'In our simulation, the preceding star forms before the cold accretion emerges, at which we have yet to perform particle splitting. The luminosity of the preceding star 𝐿 p is given by  \n𝐿 p = 𝜂𝐿 Edd ( 𝑚 ★ ) , (1)  \nwhere 𝜂 is the non-dimensional parameter and 𝑚 ★ is the mass of the star particle. We take 𝜂 = 0 . 03 as a fiducial value for the following reasons. In the preceding star formation, our mass resolution is insufficient to estimate the stellar mass (see Section 2.3). We assume 10% of the mass of the sink particle 3000 M ⊙ is converted into stars, following the result of radiation hydrodynamic simulations by Hirano et al. (2014). This assumption gives a single 300 M ⊙ star, whose luminosity is ≃ 30% of the Eddington value, following one-dimensional stellar evolution calculations (e.g. Hosokawa et al. 2013). Since the Eddington luminosity is proportional to the stellar mass, the luminosity in the fiducial model is 𝐿 p ( 𝑚 ★ ) = 0 . 3 × 𝐿 Edd ( 0 . 1 𝑚 ★ ) = 0 . 03 𝐿 Edd ( 𝑚 ★ ) , yielding 𝜂 = 0 . 03. We assume that the effective temperature is 𝑇 eff = 10 5 K, typical to the Pop III zero-age main-sequence (ZAMS) stars (Bromm et al. 2001; Hosokawa et al. 2013). In Appendix A, we explain how varying the feedback efficiency 𝜂 affects our results. We discuss the expected effect of the lifetime and the supernova feedback in Sections 5.2, 5.3.', '2.4.2 After the emergence of cold accretion': 'We regard sink particles forming after the emergence of cold accretion as individual stars, for which we employ the higher spatial resolution. We assume the stellar luminosity to be the Eddington luminosity, 𝐿 = 𝐿 Edd ( 𝑚 ★ ) = 1 . 2 × 10 41 erg s -1 ( 𝑚 ★ / 10 3 M ⊙ ) , which gives a good approximation of the massive stars with 𝑚 ★ ≳ 10 3 M ⊙ .  \nWemodelthestellar effective temperature as a function of the mass accretion history, based on the stellar evolution calculations (e.g. Hosokawa et al. 2013). We consider that the effective temperature immediately falls to 𝑇 eff = 5000 K once the accretion rate exceeds / 𝑚 ★ crit = 0 . 04 M ⊙ yr -1 , assuming an inflated stellar envelope with very rapid accretion. When the accretion rate drops below / 𝑚 ★ crit , the star experiences KH contraction, where the effective temperature is kept 5000 K until the time interval from the last epoch of / 𝑚 ★ > / 𝑚 ★ crit exceeds the surface KH timescale given by Sakurai et al. (2016),  \n𝑡 KH = 1 . 0 × 10 4 yr GLYPH<18> 𝑚 ★ 5 × 10 4 M ⊙ GLYPH<19> . (2)  \nAfter the KH contraction stage, the effective temperature depends on the accretion rate and is given by 𝑇 eff = 3 × 10 4 K for 0 . 004 < / 𝑚 ★ < 0 . 04 M ⊙ yr -1 , and 𝑇 eff = 10 5 Kfor / 𝑚 ★ < 0 . 004 M ⊙ yr -1 . Our assumption of 𝑇 eff = 10 5 K applies to ZAMS stars, as low accretion rates allow the KH contraction to persist until hydrogen burning begins. The transitional regime where 𝑇 eff = 3 × 10 4 Kismotivated by the "oscillatory" evolution of the stellar radius under these conditions (Omukai et al. 2008; Hirano et al. 2014). The mass accretion rate is measured and averaged for every 10 4 yr, the dynamical time at the sink surface.  \nWe also investigate an additional shielding of stellar radiation by the gas supposedly existing inside the sink particle, though it ultimately turns out to be a minor effect. The sink particle method cannot resolve the dense gas that will be present very close to a star in reality. We assume that the stellar UV radiation is completely absorbed when the following condition is satisfied. We consider the mass loss rate by photoevaporation of an "unresolved" disc within the sink as  \n/ 𝑚 PE = 1 . 5 × 10 -2 M ⊙ yr -1 GLYPH<18> Φ EUV 10 52 s -1 GLYPH<19> 1 / 2 GLYPH<18> 𝑟 sink 10 4 au GLYPH<19> 1 / 2 (3)  \n(Tanaka et al. 2013), where Φ EUV ( 𝐿, 𝑇 eff ) is the UV emissivity at ℎ𝜈 ≥ 13 . 6 eV, and 𝑟 sink is the sink radius which is assumed to be the size of the unresolved disc. We assume that no UV photons escape from a sink particle when the mass accretion rate onto the sink exceeds the above photoevaporation rate, supposing that the unresolved disc blocks the UV radiation. We do not consider possible directional dependence of the escape fraction, although it may be more realistic. In addition, we also assume that H -dissociating photon ( ℎ𝜈 ≃ 2 eV ) is always optically thin and unshielded by the unresolved disc.', '2.4.3 Radiation transfer': 'We solve the transfer of the UV radiation and H 2 column density by using a ray-tracing method based on Susa (2006), which has already been implemented and validated in the previous studies (Chon & Latif 2017; Chon et al. 2024). In this scheme, we calculate optical depths for ionising photon ( ℎ𝜈 ≥ 13 . 6 eV ) 𝜏 ion and for H 2 dissociating photons ( 11 . 2 eV ≤ ℎ𝜈 ≤ 13 . 6 eV ) 𝜏 LW , as well as the column density of H 2 . We attenuate the UV intensity by exp (-𝜏 ) and determine the effect of the self-shielding of H 2 based on the column density of H 2 (Wolcott-Green et al. 2011).', '3.1 Overall evolution': 'Our simulation shows that the cold accretion emerges during the formation of the atomic cooling halo, and the massive gas accretion due to cold accretion induces the formation of SMSs. Fig. 2 summarises the overall evolution of the formation of SMSs, showing the density distributions for the different time epochs and spatial scales. The top panel shows the gas distribution at 7 kpc scale, and the gas along the filamentary structures brings a large amount of the gas into the atomic cooling halo. While the normal Pop III star forms preceding the onset of the cold accretion and the radiation feedback from it evacuates the gas from the halo (Fig. 1), the cold accretion emerges and the gas flows deep inside the virial radius of the halo at 𝑧 = 16 . 9. The bottom panels show the density distribution around the halo centre for three different epochs. The cold streams carry a significant amount of the gas and accumulate the gas within 10 -2 𝑟 vir around the halo centre. More gas accumulates at the halo centre as time goes on, forming a dense gas disc 𝑛 H ≳ 10 4 cm -3 . The disc is gravitationally unstable and fragments into multiple protostars. The protostars efficiently accrete the disc gas and finally grow into SMSs with 𝑚 ★ ∼ 10 5 M ⊙ .  \nThe top panel of Fig. 3 shows the time evolution of the masses of the stars formed during our simulation. Four of the stars rapidly grow in mass and evolve into SMSs with 𝑚 ★ ∼ 10 5 M ⊙ within the initial ∼ 3-4 Myr. The bottom panel shows the time evolution of the accretion rate onto the most massive protostar. The protostar keeps the accretion rate of 0 . 01-0 . 1 M ⊙ yr -1 , which is expected for the direct collapse model. The accretion rate decreases as the ionised region expands at Δ 𝑡 ∼ 5 Myr. Four SMSs form the gravitationally bound multiple system at the end of our simulation.', '3.2 Emergence of cold accretion with radiative feedback': '3.2.1 Before cold accretion; a normal Pop III star in atomic cooling halo  \nIn our simulation, no star forms until the halo mass exceeds ∼ 10 7 M ⊙ since the background LW radiation destroys hydrogen molecules and suppresses molecular cooling. At 𝑧 = 18 . 9, when 𝑀 halo ≃ 10 7 M ⊙ and 𝑇 vir ≃ 10 4 K, Ly𝛼 cooling becomes efficient, leading to cloud collapse within the halo. This induces star formation before the emergence of cold accretion in the same manner as in K23. At an early stage of the collapse with 𝑛 H < 10 3 cm -3 , the gas temperature remains at 𝑇 ≃ 8000 K by Ly𝛼 cooling. After the density exceeds 10 3 cm -3 , a sufficient amount of H 2 forms, blocking the LW radiation by the self-shielding effect, thereby cooling the gas to 𝑇 ≲ 1000 K. We expect that the normal Pop III stars form in this situation, as H 2 molecular cooling regulates the cloud temperature (Hirano et al. 2014). As the cloud collapse proceeds and the density reaches 10 6 cm -3 , we introduce a sink particle. The mass of the sink particle grows to ≃ 3000 M ⊙ within the local free-fall time. Since we do not have enough spatial resolution to resolve individual protostars before applying the particle splitting, we assume that 10% of the accreted mass contributes to the stellar mass (Sections 2.3 and 2.4.1). We assume that this sink particle emits radiation with an effective temperature of 10 5 K and luminosity of ∼ 10 40 erg s -1 , which are typical values of a Pop III star with ≃ 300 M ⊙ .  \nFigure 2. Hierarchical view of our simulation showing the sequential formation of supermassive stars induced by the cold accretion. The top panel presents 2D-projected gas density maps, showing a large-scale filamentary structure surrounding the focused halo at the epoch of 𝑧 = 18 . 9 ( 𝑡 = 189 Myr ) . The central black solid circle denotes the virial radius of the halo. The panel embedded in the top panel presents a 2D-projected map of the gas accretion rate at 𝑡 = 223 . 1 Myr, indicating the cold accretion that feeds the halo. The white solid circle represents the virial radius of the halo. The same panel is described in detail in Fig. 4. The bottom three panels illustrate the time evolution of 2D-projected gas density maps, showing the formation of SMS within the dense disc at the centre of the halo. In each panel, the white points represent individual SMSs, and the white dashed circle denotes 0 . 01 𝑟 vir .  \n<!-- image -->', '3.2.2 Accretion flows v.s. radiative feedback': 'In K23, we performed similar cosmological simulations without radiative feedback from stars within the halo and demonstrated that cold accretion arises approximately 10 Myr after the formation of the stars. In this paper, we newly incorporate the radiation from the stars formed preceding the emergence of cold accretion. We have similarly observed the emergence of cold accretion, but the epoch of the emergence is delayed by the feedback as we will see later.  \nFig. 4 shows the evolution of the halo-scale gas structure until the emergence of cold accretion. At 𝑧 = 18 . 9 (first row of panels), immediately after the initial cloud collapse, a small HII bubble begins to expand around the preceding Pop III star. Filamentary streams feed  \nthe gas at a rate of / 𝑀 ≳ 10 -2 M ⊙ yr -1 , but they are stalled around 20% of the viral radius and do not feed the gas inside this radius. At 𝑧 = 18 . 2 (second row), the spherical HII bubble expands and the dense region with 𝑛 H ≳ 10 cm -3 near the halo centre also expands. The accretion flow along the filament is terminated by the expansion of the HII bubble and still stalled at ≃ 0 . 2 𝑟 vir . At 𝑧 = 17 . 7 (third row), the HII bubble keeps expanding but becomes asymmetrical in shape. One stream with a high-mass accretion rate reaches the halo centre. As the stream brings the gas with high density, it shields the ionising radiation and causes the HII bubble to shrink in the direction of the stream. At 𝑧 = 16 . 9 (fourth row), the filamentary streams from all directions penetrate deep into the halo centre and feed the mass to  \nFigs. 7 and 9 show that the star formation at this stage occurs through H 2 molecular cooling, despite the LW radiation field produced by the pre-existing (or preceding) star. We estimate the intensity of the LW intensity at the positions of the cores as 𝐽 21 ≃ 2000, assuming the distance from the star as ≃ 30 pc ≃ 0 . 1 𝑟 vir when H 2 cooling becomes effective with 𝑛 H ≳ 10 3 cm -3 . This value is  \n<!-- image -->  \n∆  \n∆  \nFigure 3. (Top) History of stellar mass growth following the onset of cold accretion. The horizontal axis illustrates the age of the universe 𝑡 , with time interval Δ 𝑡 measured from the star formation event at 𝑡 = 223 Myr. The red line represents the mass evolution of the star that ultimately becomes the most massive by the end of the simulation. The black lines depict the mass evolution of the other stars. (Bottom) The mass accretion history of the star highlighted by the red line in the top panel. The red (blue) parts of the line denote the periods when the accretion rate is higher (lower) than 0 . 04M ⊙ yr -1 , for which 𝑇 eff = 5000 K ( 𝑇 eff = 3 × 10 4 K or 10 5 K). The orange dashed line denotes the disc photoevaporation rate with 𝑇 eff = 3 × 10 4 K given by Eq. (3). As described in Section 2.4.2, we assume that stellar UV radiation is shielded by an "unresolved" disc within a sink particle when the accretion rate / 𝑚 ★ exceeds the photoevaporation rate / 𝑚 PE, as for / 𝑚 ★ < 0 . 04 M ⊙ yr -1 . When / 𝑚 ★ > 0 . 04 M ⊙ yr -1 , EUV emissivity is low enough to always satisfy / 𝑚 ★ ≫ / 𝑚 PE with 𝑇 eff = 5000 K.  \n<!-- image -->  \nthe halo centre at a rate of / 𝑀 ≳ 10 -2 M ⊙ yr -1 . The dense gas carried by the accretion flows confines the HII bubble into a small region ( < 0 . 1 𝑟 vir ). This marks the period of the onset of cold accretion.  \nFig. 5 shows the time evolution of 𝑟 shock , the distance from the halo centre where the accreting flow is terminated by the shock. To derive this radius, we spherically averaged the mass inflow rate and identified the position where the velocity of the flow transitions from supersonic to subsonic as in K23. The figure indicates that, when radiative feedback is considered, the onset of cold accretion is delayed by ≃ 20 Myr compared to our simulation without radiative feedback. This additional delay is attributed to the outward expansion  \ncaused by the thermal pressure of the HII bubble, which opposes the inward ram pressure of the accretion streams.  \nK23 derived the condition for the first emergence of the cold accretion in the 𝑧 -𝑀 halo plane, by the spherical accretion model as in Birnboim & Dekel (2003), which well describes the onset of the cold accretion found in the simulation without radiation feedback. In Fig. 6, we plot the time evolution of the halo mass and highlight the epoch when the cold accretion emerges by stellar symbols. The figure shows that the condition given by K23 well explains when cold accretion emerges when we include the radiation feedback, while it delays the onset of the emergence.', '3.3.1 Initial cloud collapse induced by cold accretion': 'Cold accretion supplies a substantial mass of the gas and accumulates gas near the centre of the halo, triggering subsequent star formation. The left panels of Fig. 7 present the mass histograms in the 2D plane of the density-temperature (top) and density-H 2 fraction (bottom panels) at 𝑧 = 16 . 9, at the moment when cold accretion emerges ( 𝑡 = 223 . 1 Myr). The presence of extremely dense gas, surpassing 10 6 cm -3 , indicates that the gravitational collapse of the cloud occurs by this epoch. The bottom left panel of Fig. 7 shows that the gas with 𝑛 H ≲ 10 3 cm -3 has a low H 2 fraction, while the gas with 𝑛 H ≳ 10 3 cm -3 lies the region with high H 2 fraction. This is above a critical H 2 fraction 𝑦 H2 > 10 -5 ≡ 𝑦 H2 , crit , where H 2 cooling timescale is shorter than dynamical timescale (Yoshida et al. 2003). This indicates that cloud collapse is triggered by Ly𝛼 cooling until ≲ 10 3 cm -3 , while H 2 cooling becomes effective and becomes the dominant cooling process after the density exceeds ∼ 10 3 cm -3 .  \nFig. 8 presents a 2D-projected gas structure at the halo centre within 𝑟 ∼ 10 -2 𝑟 vir at the epoch of 𝑡 = 223 . 1 Myr, just before subsequent star formation. The accreting filament fragments into two dense cores and stars form inside the cores afterwards. The fragmentation is caused by the rapid temperature decrease caused by H 2 molecular cooling, as the temperature of the filament gas decreases to ∼ 1000 K. This cold and dense gas corresponds to the dense component with 𝑛 H ≳ 10 7 cm -3 in Fig. 7.  \nThe upper panel of Fig. 9 illustrates the radial distribution of the enclosed gas mass 𝑀 gas. The red dashed line denotes the critical Bonnor-Ebert mass (Ebert 1955; Bonnor 1956), above which the cloud cannot be in a hydro-static equilibrium and becomes gravitationally unstable,  \n𝑀 BE ( 𝑟 ) = 1 . 18 𝑐 4 s ( 𝑟 ) 𝑃 1 / 2 ( 𝑟 ) 𝐺 3 / 2 , (4)  \nwhere 𝑐 s ( 𝑟 ) and 𝑃 ( 𝑟 ) are sound velocity and gas pressure, calculated as mass-weighted averages within each radial bin for 𝑟 < 𝑟 \' < 𝑟 + 𝛿𝑟 . The enclosed gas mass becomes larger than the Bonnor-Ebert mass at the radius ≃ 1 . 5 pc ≃ 4 × 10 -3 𝑟 vir . This indicates that the cloud becomes gravitationally unstable at this scale and fragments into clumps with a typical mass of 𝑀 BE ∼ 10 4 M ⊙ . The blue point shows the mass and distance from the halo centre of the preceding star. The star formation after the emergence of cold accretion occurs about 15 pc ∼ 0 . 04 𝑟 vir away from the position of the preceding star.  \n10  \nFigure 4. Mass-weighted projection maps of the halo-scale gas structures at epochs of 𝑧 = 18 . 9 , 18 . 2 , 17 . 7 , and 16 . 9 in descending order. The epoch of the first row of the panels corresponds to the initial cloud collapse and the resulting normal Pop III ("preceding") star formation. In each panel, the central white point denotes the preceding star particle, and the white circle denotes the virial radius of the halo. The epoch of the bottom row corresponds to that of the emergence of the cold accretion, i.e., when the filamentary accretion streams penetrate the halo. The left, middle, and right columns represent the gas density, the accretion rate, and the degree of ionisation, respectively. The gas accretion rate at a given radius 𝑟 is defined as / 𝑀 ≡ 𝜌𝑟 2 𝑣 inf , where 𝜌 is the gas density, and 𝑣 inf the radial infalling velocity.  \n<!-- image -->  \ncomparable to the critical intensity 𝐽 21 , crit ≃ 2000 for the radiation temperature 𝑇 eff = 10 5 K, below which H 2 cooling operates to cool the gas to a few 100 K (Sugimura et al. 2014). The intensity at the fragmented cores is comparable to the critical density estimated by the one-zone model by Sugimura et al. (2014), while the threedimensional simulations show that the turbulence can promote the  \nformation of H 2 and increase the required intensity to disable H 2 cooling (Latif et al. 2018).', '3.3.2 Protostellar accretion in dense disc fed by cold accretion': 'Fig. 3 shows the stellar mass growth histories after the onset of cold accretion. The upper panel shows that stars rapidly increase their', '𝑛 H [cm -3 ]': "Figure A1. Mass-weighted 2D-projection maps of the halo at epochs of 𝑧 = 18 . 9 , 18 . 7 , and 18 . 2 in descending order, in the case with the feedback luminosity 𝐿 p = 0 . 1 𝐿 Edd . The epoch of the first row corresponds to when the preceding star (not SMS) forms. The left, middle, and right columns represent the gas density, accretion rate, and degree of ionisation, as the tracer of virialised gas, accretion flows, and HII bubble, respectively. The white point in each panel denotes the preceding star particle. The white circle in each panel denotes the halo's virial radius.  \n<!-- image -->  \n10  \nሶ  \n-1", '3.3.3 Photoevaporation': 'Fig. 11 shows the density and temperature distributions around the central stellar system at the later evolutionary stage of the mass accretion phase of the SMSs. The dense gas disc surrounding the SMSs at 𝑡 = 226 . 7 yr, is finally photoevaporated by the radiation from the SMSs at 227 . 8 Myr. The photoevaporation is triggered by the radiation from the ejected star with 𝑚 ★ ≃ 5000 M ⊙ . This star is ejected by the three-body encounter around the epoch of 𝑡 = 226 . 7 Myr and wanders around the low-density region. The accretion  \nrate onto the ejected star becomes lower due to the small gas density. It falls below 4 × 10 -3 M ⊙ yr -1 and the stellar radius shrinks and reaches the ZAMS stage due to the KH contraction. This increases the effective temperature of the stars to 10 5 K, and an HII region around the star begins to expand. The disc is exposed to strong ionising radiation afterwards and gradually photoevaporated. The bottom panels show that the disc gas is completely cleared by the strong radiation by 𝑡 = 227 . 8 Myr. Mass accretion onto all the stars completely ceases by this period (Fig. 3). We note that the photoevaporation of the disc is accelerated once the HII region starts to expand in the following manner. Once the feedback from one star expels the gas that surrounds other stars, they also transition to the ZAMS phase due to the low accretion rate. It develops HII bubbles with high ionising emissivities, further enhancing the feedback.  \nFig. 12 shows the distributions of mass accretion rate and ionisation degree projected onto the face-on and edge-on view of the central gas disc, showing how radiative feedback from the SMSs affects the accretion flow at the halo scale. The right panels indicate that a large bipolar HII region develops and eventually extends beyond the virial radius at 𝑡 = 227 . 1 Myr. Meanwhile, the left panels reveal that the cold accretion flows with / 𝑀 ≳ 10 -2 M ⊙ yr -1 are progressively interrupted from the inside outward. The expansion of H II region finally halts the mass accretion onto the central halo region by 𝑡 = 227 . 8 Myr.  \nFigure 8. 2D-projected maps of the gas density (top) and temperature (bottom) around the halo centre at the epoch of 𝑡 = 223 . 1 Myr, when the star formation is about to resume after the emergence of the cold accretion. The white dashed circles denote 𝑟 = 10 -2 𝑟 vir .  \n<!-- image -->', '4 WHAT DETERMINES ACCRETION RATE?': 'Our simulation has shown that the protostars efficiently grow after the emergence of cold accretion, maintaining a high mass accretion rate with 0 . 04 M ⊙ yr -1 for ≃ 3 Myr. Since this timescale is longer than the typical dynamical time of the disc, 𝑟 disc / 𝑐 s ∼ 0 . 02 𝑟 vir / 𝑐 s ≃ 0 . 5 Myr, the high mass accretion rate should be caused by the physics associated with larger scales. This section clarifies the underlying physics that sustains the high accretion rate: the cold accretion and gravitational torque in the disc.', '4.1 Cold accretion': 'Cold accretion is a cosmological phenomenon in which gas falls from scales larger than the virial radius directly to the centre of the halo (Birnboim & Dekel 2003). Unlike the classical accretion model (Rees & Ostriker 1977), the accreting gas is never stalled by shock fronts around the virial radius of the halo. This is evident in Fig. 2 that cold accretion keeps the gas accretion rate radially constant from virial radius 𝑟 vir to the rotation-supported disc radius 𝑟 disc ∼ 10 -2 𝑟 vir , resulting in / 𝑀 gas ( 𝑟 vir ) ≃ / 𝑀 gas ( 𝑟 disc ) . This allows us to relate the gas accretion rate onto the central disc / 𝑀 CA to the mass accretion rate of the halo / 𝑀 halo as / 𝑀 CA ≡ / 𝑀 gas ( 𝑟 disc ) ≃ 𝑓 br / 𝑀 halo , where 𝑓 br = 0 . 15 is the baryon fraction. By estimating / 𝑀 halo with the spherical accretion model (Gunn & Gott 1972; Lahav et al. 1991), we analytically determine / 𝑀 CA as follows, where we approximate Λ CDM cosmology we assume for our simulation with Einsten-de Sitter universe for simplicity.  \nThe virial radius of a halo with mass 𝑀 halo at redshift 𝑧 vir is given by 𝑟 vir ≃ ( 𝑀 halo / 24 𝜋 3 ¯ 𝜌 0 ) 1 / 3 ( 1 + 𝑧 vir ) -1 , where ¯ 𝜌 0 is the mean cosmic matter density at 𝑧 = 0. We consider a mass shell falling into the halo which encloses 𝑀 halo at 𝑟 vir . The infall velocity of the shell  \nFigure 9. The comparison between the radial profiles of enclosed gas mass (black solid lines) and estimated Bonnor-Ebert mass (red dashed lines). The panels show the snapshots at 𝑡 = 223 . 1 , 223 . 9, and 224 . 6 Myr in descending order. The horizontal axes represent the distance from the centre normalised by the virial radius. The blue squares mark the positions and masses of the star particles. Gas is considered gravitationally bound at a given radius if the enclosed mass exceeds the Bonnor-Ebert mass, 𝑀 gas ( < 𝑟 ) ≳ 𝑀 BE ( 𝑟 ) . The central position is defined as the peak of gas density for the initial snapshot at 𝑡 = 223 . 1 Myr, before star formation starts in the disc, and as the most massive star for 𝑡 ≥ 223 . 9 Myr.  \n<!-- image -->  \n𝑣 can be written as  \n𝑣 ( 𝑧 vir ) = √︄ 𝐺𝑀 halo 𝑟 vir ≡ 𝑣 vir . (5)  \nAssuming the radial profile of the DM density to follow 𝜌 halo ( 𝑟 ) ∝  \nFigure 10. The 2D-projected maps of the dense star-forming gas disc around the halo centre at the epoch of 𝑡 = 224 . 6 Myr. The upper, middle and lower rows of the panels represent the gas density, degree of ionisation, and H2 fraction. The left and right columns represent the different viewing angles of the face-on and edge-on of the disc, respectively. In each panel, the white dashed circle denotes 𝑟 = 10 -2 𝑟 vir . The star particles are denoted by white points in the top and middle rows and black points in the bottom rows. The coordinate origin is set at the position of the most massive star.  \n<!-- image -->  \n𝑟 -2 , we obtain the matter density at 𝑟 = 𝑟 vir as  \n𝑀 halo = ∫ 𝑟 vir 4 𝜋𝑟 2 𝜌 halo ( 𝑟 ) 𝑑𝑟 = 4 𝜋𝑟 3 vir 𝜌 halo ( 𝑟 vir )  \n𝜌 halo ( 𝑧 vir , 𝑟 vir ) = 𝑀 halo 4 𝜋𝑟 3 vir .  \n0 (6) (7)  \nNote that the actual density profile of the halos follows the NFW profile (Navarro et al. 1997), introducing a factor of order unity into Eq. (7), which we neglect here for convenience. Using the above equations, we can rewrite / 𝑀 CA ( 𝑡 ) as  \n/ 𝑀 halo ( 𝑡 ) = 4 𝜋𝑟 2 vir 𝑣 vir 𝜌 halo ( 𝑧 vir , 𝑟 vir ) ≃ 𝑀 halo 𝑣 vir 𝑟 vir = 2 𝜋 𝑀 halo 𝑡 (8) / 𝑀 CA ( 𝑡 ) ≃ 𝑓 br / 𝑀 halo ( 𝑡 ) = 0 . 055 M ⊙ yr -1 GLYPH<18> 𝑀 halo 2 × 10 7 M ⊙ GLYPH<19> GLYPH<18> 1 + 𝑧 vir 18 GLYPH<19> 3 / 2 , (9)  \nwhere we use 𝑡 ≃ GLYPH<2> 6 𝜋𝐺 ¯ 𝜌 0 ( 1 + 𝑧 vir ) 3 GLYPH<3> -1 / 2 , which holds in matterdominant epoch ( 𝑧 vir > 1 ) , and  \n𝑣 vir 𝑟 vir = √︄ 𝐺𝑀 halo 𝑟 3 vir = √︃ 24 𝜋 3 𝐺 ¯ 𝜌 0 ( 1 + 𝑧 vir ) 3 = 2 𝜋 𝑡 (10)  \nThe accretion rate / 𝑀 CA can also be represented by virial temperature  \nFigure 11. Destruction of the central gas disc by the UV radiative feedback associated with the expansion of an HII bubble. The left and right columns of the panels show 2D-projected maps of the gas density and temperature, respectively. Snapshots at 𝑡 = 226 . 7, 227.1, and 227 . 8 Myr are shown in the upper, middle, and lower rows, respectively, with a face-on view of the disc structure in all panels. The white dashed circles in the left column and the black dashed circles in the right column denote 𝑟 = 10 -2 𝑟 vir . The star particles are denoted by white points in the left column and black points in the right column.  \n<!-- image -->  \nof the halo  \n𝑇 vir ≡ 𝜇𝑚 H 𝑣 2 vir 2 𝑘 B = 1 . 3 × 10 4 K GLYPH<18> 𝑀 halo 2 × 10 7 M ⊙ GLYPH<19> 2 / 3 GLYPH<18> 1 + 𝑧 vir 18 GLYPH<19> (11)  \nas  \n/ 𝑀 CA ≃ 0 . 040 M ⊙ yr -1 GLYPH<18> 𝑇 vir 10 4 K GLYPH<19> 3 / 2 . (12)  \nSince the cold accretion occurs when 𝑇 vir ≳ 10 4 K, the gas accretion rate by cold accretion is always larger than / 𝑀 CA ≃ 0 . 04 M ⊙ yr -1 .  \nFig. 13 presents the time evolution of the radial profiles of the gas accretion rate within the ACH in our simulation. The top panel shows the snapshot before cold accretion appears, indicating that the accretion rate at 𝑟 ≳ 0 . 2 𝑟 vir is comparable to / 𝑀 CA (Eq. 12). We see that there is an outflow component at 𝑟 ≲ 0 . 2 𝑟 vir , which is caused by the expansion of the HII bubble (see also the second row of Fig. 4). The middle and bottom panels show the evolution after the emergence of cold accretion. The accretion rate follows / 𝑀 CA and is almost constant from 𝑟 ∼ 𝑟 vir down to the scales smaller than the  \nሶ  \n𝑀[M  \n-3  \n-1  \n]  \n10  \n-1  \n10  \n⊙  \n10  \nyr  \n-2  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 12. The halo-scale photoevaporation, the UV stellar feedback caused by the SMSs formed at the centre of the halo. The left and right columns of panels display the 2D-projected maps of the accretion rate and degree of ionisation, respectively. Snapshots at 𝑡 = 226 . 7, 227.1, and 227 . 8 Myr are shown in the upper, middle, and lower rows, respectively, with a edge-on view of the central disc in all panels. The white or black solid circle in each panel denotes the virial radius. The star particles are denoted with the white or the black points.  \n<!-- image -->  \ncentral disc, 10 -3 𝑟 vir . These panels indicate that once cold accretion begins at 𝑡 ≃ 223 Myr, the cosmological halo-scale accretion rate governs the accretion rate onto the central disc for ≃ 3 Myr.', '4.2 Gravitational torque in the disc': 'The cold accretion brings a large mass accretion rate from the cosmological scales to the halo centre. It is tempting to infer that the cold accretion directly provides the gas to the SMSs. However, the cold accretion provides the gas only to the disc scale 𝑟 disc ≳ 10 pc, where the centrifugal barrier becomes effective. The Bondi scales of the accreting SMSs are much smaller, 𝑟 Bondi ∼ 1 pc, indicating that we need additional physics that migrates the gas further inward. In this section, we show that the gravitational instability excites the spiral arms and resulting torque enables efficient transfer of the mass and angular momentum across the disc.  \nThe Toomre 𝑄 parameter (Toomre 1964), which measures the  \n10  \n-3  \n-1  \n𝑥  \ne  \n-2  \n<!-- image -->  \n10  \n10  \n1  \n0  \n.  \n8  \nkpc  \ndegree of the gravitational instability, is given by  \n𝑄 ≡ 𝑐 s 𝜅 𝜋𝐺 Σ gas , (13)  \nwhere Σ gas is the surface gas density of the disc, 𝜅 is epicyclic frequency 𝜅 2 ≡ 𝑟 -3 𝑑 ( 𝑟 4 Ω ( 𝑟 ) 2 )/ 𝑑𝑟 , and Ω is the angular velocity. Fig. 14 shows radial profiles of the physical quantities involving Toomre 𝑄 parameter after cold accretion emerges. We calculate surface gas density in each radial bin Σ gas ( 𝑟 ) using enclosed gas mass 𝑀 enc ( < 𝑟 ) by  \nΣ gas ( 𝑟 ) ≡ 1 2 𝜋𝑟 𝑑𝑀 enc ( < 𝑟 ) 𝑑𝑟 , (14)  \nwhich holds inside the disc radius 𝑟 ≲ 10 -2 𝑟 vir ≃ 𝑟 disc where most of the gas is concentrated in the disc plane. We calculate 𝜅 ( 𝑟 ) by differentiating Ω ( 𝑟 ) , which is evaluated as  \nΩ ( 𝑟 ) 2 ≡ 𝐺 ( 𝑀 enc ( < 𝑟 ) + 𝑀 ★ ( < 𝑟 )) 𝑟 3 , (15)  \nwhere 𝑀 ★ ( < 𝑟 ) is a total stellar mass within the radius 𝑟 . This figure shows that the Toomre 𝑄 parameters is an order of unity for 10 -3 𝑟 vir ≲ 𝑟 ≲ 𝑟 disc at the onset ( 𝑡 = 223 . 1 Myr) and in the middle (223 . 9 Myr) of the SMS formation, indicating that the gravitational instability operates across the disc. We have confirmed that 𝑄 ∼ 1 is maintained for ≃ 3 Myr after the onset of SMS formation.  \nAccording to Shakura & Sunyaev (1973), the accretion rate within a disc can be described as  \n/ 𝑀 disc = 𝛼 Σ gas 𝑐 2 s Ω , (16)  \nwhere 𝛼 is a non-dimensional parameter which depends on the process of angular-momentum transportation.  \nCombining Eqs. 13 and 16 and assuming 𝜅 ≃ Ω , which is true for the disc with Keplerian rotation, we can evaluate the accretion rate as  \n/ 𝑀 disc ∼ GLYPH<18> 𝛼 𝜋𝑄 GLYPH<19> 𝑐 3 s 𝐺 ≃ 0 . 068 M ⊙ yr -1 GLYPH<18> 𝑇 8000 K GLYPH<19> 3 / 2 𝛼 𝑄 . (17)  \nEq. (17) indicates that the gas accretion rate can be described by its temperature 𝑇 , 𝑄 , and 𝛼 . Assuming 𝛼 ∼ O( 1 ) for the self-gravitating disc (e.g. Boley et al. 2006), the accretion rate / 𝑀 disc depends on only the temperature 𝑇 , with comparable rate to / 𝑀 disc ≃ / 𝑀 CA at 𝑇 ≃ 𝑇 vir 1 .  \nThe red lines in Fig. 13 represent the accretion rate / 𝑀 disc when we assume 𝛼 = 4. The accretion rate is well estimated by our / 𝑀 disc ( 𝑟 ) in the range of 10 -3 𝑟 vir ≲ 𝑟 ≲ 𝑟 disc at the onset ( 𝑡 = 223 . 1 Myr) and in the middle (223 . 9 Myr) of the SMS formation, which supports the idea that gravitational instability drives the rapid mass accretion through the disc.  \nWehaveseen that / 𝑀 disc ∼ / 𝑀 CA is maintained for scales more than three orders of magnitude below the virial radius. Assuming the gas disc is H 2 deficient and hot, keeping 𝑄 as low as 𝑄 ∼ 1 is required for maintaining a high accretion rate / 𝑀 disc . If / 𝑀 disc ≪ / 𝑀 CA with 𝑄 ≫ 1, cold accretion supplies more gas than the disc can transport inward. It accumulates more gas on the disc and increases the surface density, decreasing 𝑄 . This makes the disc more unstable and causes a higher accretion rate. If / 𝑀 disc > / 𝑀 CA with 𝑄 < 1, more disc gas is transported inward than the gas supplied by cold accretion.  \nFigure 13. The radial profiles of gas inflow (blue) and outflow (green lines) rates in the halo at 𝑡 = 199 . 4 , 223 . 1 , and 223 . 9 Myr. The horizontal axes represent the distance from the centre normalised by the virial radius. We take the position of the most massive star as the centre at 𝑡 = 119 . 4 and 223 . 9 Myr. For the snapshot at 𝑡 = 223 . 1Myr,wechoosethepositionofthemaximumgas density as the centre. The black and red lines denote the analytically estimated accretion rate, where / 𝑀 CA = 𝑓 br / 𝑀 halo = 0 . 04 M ⊙ yr -1 ( 𝑇 vir / 10 4 K ) 3 / 2 and / 𝑀 disc = 𝛼 Σ gas 𝑐 2 s / Ω , where 𝛼 = 4, respectively (see Section 4.2).  \n<!-- image -->  \nThis decreases the surface density and increases 𝑄 , thus reducing the accretion rate inside the disc. These processes bring 𝑄 to the equilibrium value 𝑄 eq where / 𝑀 disc ∼ / 𝑀 CA satisfied as  \n𝑄 eq ≃ 2 𝛼 GLYPH<18> 𝑇 𝑇 vir GLYPH<19> 3 / 2 ∼ 1 . 4 𝛼 GLYPH<18> 𝑇 8000 K GLYPH<19> 3 / 2 GLYPH<18> 𝑇 vir 10 4 K GLYPH<19> -3 / 2 . (18)  \nNote that in our simulation, the disc and spiral arms sometimes fragment. In general, the spiral arm in the disc will fragment by selfgravity if 𝑄 ≲ 0 . 6 is satisfied (Takahashi et al. 2016). This suggests that 𝑄 ≲ 0 . 6 is locally satisfied at some epochs. Even though, as long as the equilibrium value does satisfy 𝑄 eq > 0 . 6, fragmentation should be globally suppressed and inefficient. In summary, the self-gravity of the central disc regulates the accretion rate radially constant for scales at 10 -3 ≲ 𝑟 / 𝑟 vir ≲ 10 -2 . The accretion rate is ≃ 0 . 04 M ⊙ yr -1 , which is determined by the accretion rate brought by the cold accretion as in Eq. 12.  \nFigure 14. The radial profiles of gas surface density Σ gas (blue), temperature 𝑇 (black), and Toomre Q parameter (red lines) at two different epochs, 𝑡 = 223 . 1 , and 223 . 9 Myr. The horizontal axes represent the distance from the centre normalised by the virial radius.  \n<!-- image -->', '5.1 Formation channels of molecular hydrogen': 'We have incorporated H -channel for H 2 formation and do not include H 2 formation via three-body reactions nor H + 2 channel. The contribution from those formation channels to the total H 2 abundance is small and negligible in environments we consider in our simulation, as we will discuss below. Since three-body reaction becomes dominant to produce H 2 only in a density of 𝑛 H ≳ 10 10 cm -3 (Inayoshi et al. 2014), we can neglect this channel as we replace gas with a sink particle above the density 10 8 cm -3 . We can also neglect H + 2 channel for the following reasons. This channel can be dominant only in the environments with high redshift where CMB temperature is high (Tegmark et al. 1997) or the background radiation with loweffective temperature with 𝑇 crit ≲ 7000 K (Sugimura et al. 2016). In our calculation, the LW background has a higher effective temperature than 𝑇 crit , so the formation rate of H 2 is dominated by the H -channel. After the SMSs form and they accrete mass at a high rate, we model their spectra with blackbody radiation with an effective temperature of 𝑇 eff = 5000 K (Section 2.4). When their radiation dominates, the formation rate of H 2 by H + 2 channel becomes higher than that by H -channel (Sugimura et al. 2016). Even in such a case, the dissociation rate by the FUV radiation exceeds the formation rate via H + 2 channel since the luminosity of the accreting SMSs is very high. We thus expect the contribution to H 2 production via H + 2 channel does not affect our results, as the intense radiation from the  \nSMSs always decreases the H 2 abundance small enough and makes that contribution to the thermal evolution negligible.', '5.2 On our simple modelling of Pop III star formation induced by initial cloud collapse': 'In our simulation, we give an assumption on the formation of the normal Pop III star formation induced in the ACH before the onset of cold accretion. Instead of following the detailed evolution of Pop III stars, we model the resulting stellar system to have a blackbody spectrum with a luminosity 𝐿 and effective temperature 𝑇 eff and ignore a stellar lifetime. As a result, this star continues to emit radiation for ≳ 30 Myr. In this section, we discuss the uncertainty in modelling this Pop III star and how it affects our results of forming SMSs.  \nDue to the lack of spatial resolution, we do not follow the possible fragmentation, which may occur inside the sink particle introduced in our simulation. Previous studies, which track Pop III star formation in ACHs, indicate that fragmentation should occur and a star cluster should form rather than a single massive star (Regan et al. 2020; Latif et al. 2021). Regan et al. (2020) show the stellar mass function in ACHs, which is similar to that among different mini-halos (e.g. Hirano et al. 2014, 2015). Since the stellar lifetimes and fates depend on their mass, the luminosity and the spectra of the Pop III cluster should depend on the initial mass function (IMF) and evolve with time. The star formation efficiency (SFE), or the ratio of the total stellar mass to the cloud mass, is another uncertainty, which is not considered in our simulation. One can obtain more realistic IMF and SFE by solving star formation with higher resolution (Chon et al. 2024) and stellar evolution consistently. Liu et al. (2024) have compiled the results of Pop III star formation of gas clouds in halos, including ACHs, and built an analytic model of IMF and SFE. They have shown that IMF and SFE in ACHs can vary among different halos, depending on their properties, such as the cloud mass, accretion rate depending on temperature, and degree of fragmentation in the clouds. To cover the possible uncertainties in SFE and IMF, we have performed an experimental simulation changing our model of stellar radiation, which is shown in Appendix A.  \nIn our simulation, the star continues to emit radiation after the onset of cold accretion. We define the delay Δ 𝑡 as the duration from the initial collapse of a cloud in the ACH to the onset of cold accretion and 𝑡 life as the lifetime of the most massive star in the cluster. When Δ 𝑡 < 𝑡 life is satisfied, our model qualitatively approximates the situation as the stars emit radiation until the cold accretion emerges. To justify the assumption in the modelling of the Pop III stars, the stellar lifetime should be larger than 30 Myr. This condition can be rephrased as the stellar mass of the individual stars, which should be smaller than 10 M ⊙ (Schaerer 2002). Although this mass range may be lower than the typical mass of Pop III stars, our simulation possibly approximates some cases described below. While we investigate the evolution in one specific halo in this study, the previous work without radiative feedback (K23) shows the variety in delay Δ 𝑡 in three halos, suggesting the variety in Δ 𝑡 even with radiative feedback. As for the case with 10 times shorter Δ 𝑡 ∼ 3 Myr, Pop III stars with ∼ 100 M ⊙ , consistent with recent numerical simulation and have the lifetimes of 2-3 Myr (Sugimura et al. 2023), can survive until the cold accretion emerges. We expect that some halos from a larger sample will satisfy the sufficiently shorter delay time and the massive Pop III stars emit radiation until the cold accretion emerges, which allows the formation of supermassive stars. We also provide statistical discussion about this in Section 5.6 later.  \nIn reality, the situation with Δ 𝑡 > 𝑡 life might be typical among  \nACHs. In this case, halos would follow different evolutions, such as the SNe explosion, before the emergence of cold accretion and SMS formation. If the delay is very long Δ 𝑡 ≫ 𝑡 life , it is possible to have episodes of normal Pop III star formation ( 𝑚 ★ ≲ 10 3 M ⊙ ) before the cold accretion appears. We discuss the consequence in such cases in Section 5.3.', '5.3 Potential effects of supernova feedback and metal enrichment': 'In our simulations, we do not include the SNe feedback before the onset of cold accretion. However, if the emergence is much later than Δ 𝑡 ∼ 3 -10 Myr, the Pop III stars will cause SNe explosions depending on their stellar masses (Heger & Woosley 2002). This will change the evolution of star formation as it injects thermal and kinetic energy and enriches the surrounding material with heavy elements. We discuss the possibility of SMS formation following the emergence of cold accretion in the presence of SNe and associated metal pollution. It has been pointed out that, for SMS formation, the star-forming gas should be metal-poor (Omukai et al. 2008; Chon & Omukai 2020). Otherwise, SMS formation is inhibited by rapid cooling due to heavy elements, leading to substantial fragmentation. The critical metallicity threshold for SMS formation is around 𝑍 crit ∼ 10 -4 - 10 -3 Z ⊙ .  \nIt is important to note the various uncertainties linked to the evolution of SNe within ACHs. Firstly, there is uncertainty regarding the frequency of core-collapse SNe (CCSNe; 𝑚 ★ ∼ 10 M ⊙ ) and pair-instability SNe (PISNe; 𝑚 ★ ∼ 10 2 M ⊙ ) in ACHs due to the unknown nature of the Pop III IMF and SFE (see Section 5.2). For Pop III stars that form in mini-halos, the predicted mass distribution shows two peaks around ∼ 10 M ⊙ and ∼ 10 2 M ⊙ (Hirano et al. 2014, 2015), which implies a frequent occurrence of SNe explosions per unit of total stellar mass. Regan et al. (2020) illustrate that the Pop III IMF in ACHs exhibits double peaks that are somewhat outside the mass range typically associated with SNe on the larger mass side, indicating a lower occurrence rate.  \nFurthermore, the development of SNe bubbles and metal spreading in ACHs is uncertain. Recent numerical models of halos with 𝑀 halo ∼ 10 6 M ⊙ reveal varying evolutions based on the type of SNe involved. In halos of 𝑀 halo ∼ 10 6 M ⊙ , an individual core-collapse supernova (CCSN) releases only a minimal portion of gas from the halo (Chiaki &Wise2019). Metal enrichment occurs locally within the supernova remnant (SNR), characterized by a mass and size of ≳ 2 × 10 4 M ⊙ and 1 -10 pc, respectively. This size is less than the virial radius of mini-halos, approximately ∼ 10 2 pc. The resultant metallicity is estimated to be 𝑍 ≲ 10 -2 Z ⊙ (Magg et al. 2020) with a large scatter (Ritter et al. 2015; Tarumi et al. 2020). In contrast, an individual PISN expels most of the gas from halos with 𝑀 halo ∼ 10 6 M ⊙ (Magg et al. 2022). In this case, the expelled gas returns to the halos as they expand. The time taken for the gas to fall back varies significantly, ranging from ∼ 5 -100 Myr, probably influenced by the mass of the SN progenitor and the density structure surrounding the SNe. The enrichment occurs relatively uniformly across the virial radius. The resulting metallicity within the SNR is estimated to be ≲ 3 × 10 -3 Z ⊙ . Note that the uncertainty in the size of SNR potentially reduces the metallicity by four orders of magnitude, such as 𝑍 ≳ 10 -6 Z ⊙ for the CCSNand ≳ 3 × 10 -7 Z ⊙ for the PISN (Magg et al. 2020). The above evolutions are the results for halos with 𝑀 halo ∼ 10 6 M ⊙ , which may not be applicable to the ACHs. We also note that if multiple SNe associate, the development of the explosion could differ.  \nThese uncertainties present significant challenges to our understanding of metal enrichment in ACHs. We here discuss possible  \nSMS formation channels induced by the cold accretion, under the assumption that either CCSNe or PSINe dominate.  \nIf CCSNe are dominant, the majority of gas will remain within the ACHs. Consequently, following the first generation of Pop III stars, the halos will experience several episodes of star formation until the emergence of cold accretion. When the total stellar mass becomes significantly large, the radiative feedback from the stars inhibits H 2 cooling, thereby halting further PopIII star formation during their lifetime. Conversely, if the stellar mass is relatively small, the weaker feedback permits star formation to proceed more quickly. Therefore, star formation can proceed at a steady rate, at which the feedback modulates subsequent star formation. Our simple modelling of the preceding Pop III star formation with constant luminosity may represent this situation, except for the effects of the metal enrichment. Although modelling of advancing metal enrichment remains incomplete, as mentioned above, we expect that the earlier onset of cold accretion is more advantageous for SMS formation.  \nIf PISNe dominate, most of the gas will be expelled from the ACHs by one or a few PISNe, halting star formation until the gas falls back into the halos. If the fallback time is long enough, the halos will grow and cold accretion will emerge around the time of the fallback. When star formation resumes with cold accretion, a situation similar to our simulation can arise, where SMS can form. The metallicity at this time is set by the initial PISN event at 𝑍 ∼ 3 × 10 -7 -3 × 10 -3 Z ⊙ . This value may be below the critical metallicity for SMS formation, 𝑍 crit ∼ 10 -4 -10 -3 Z ⊙ .  \nThe above suggests that SMS formation may be possible in different channels, depending on how effective CCSNe and PISNe are in modifying the evolution. Investigating the possibility of SMS formation under the influence of SNe is one of our important future works.  \nIf the previous episodes of SNe pollute the gas above the critical metallicity and cold accretion emerges later on, fragmentation may occur, and the star cluster will form rather than the SMSs. We expect the star cluster to be very compact as cold accretion and the mass transfer inside the central disc concentrate the stellar distribution. Our simulation has shown that after the emergence of the cold accretion, the gas is transferred to the central region of ≲ 10 -3 𝑟 vir at a rate of 0 . 1 M ⊙ yr -1 . This accumulates the mass of 10 5 M ⊙ inside the ∼ 0 . 1 pc region within a few Myr (Fig. 3). If a similar mechanism works for the metal-enriched gas and concentrates the mass in the halo centre, the stellar surface density becomes Σ ★ = 𝜖 × 10 5 M ⊙ /( 0 . 1 pc ) 2 = 10 6 M ⊙ pc -2 ( 𝜖 / 0 . 1 ) , where 𝜖 is the conversion efficiency from the gas to stars. Assuming the fiducial value 𝜖 = 0 . 1, we expect the star cluster with a high mass density of 10 6 M ⊙ pc -2 . EUV and SN feedback will be inefficient in halting star formation during the entire star formation. The gas surface density at the onset of star formation will exceed Σ gas ∼ 300 -10 3 M ⊙ pc -2 , above which the stellar EUV feedback and its radiative pressure on gas from a stellar cluster is weaker than its gravitational force and cannot halt the star formation (Kim et al. 2018; Grudić et al. 2018; Fukushima & Yajima 2021; Menon et al. 2023). Since the mass accretion is completed within a few times free-fall time 𝑡 ff ≃ 0 . 5 Myr ( 𝑛 H / 10 4 cm -3 ) -1 / 2 , shorter than the stellar lifetime of Myr, the star cluster forms before the onset of SNe, similar to the feedback-free star formation (Dekel et al. 2023). This stellar density is comparable to the compact star cluster observed by JWST (Vanzella et al. 2022, 2023; Fujimoto et al. 2024; Adamo et al. 2024; Harikane et al. 2024). While the compactness of the star cluster should depend on where the fragmentation occurs, which may change the mass transfer rate inside the central disc as the disc gas is converted into the stars, our simulation result may offer a promising  \nenvironment for the formation of the compact star cluster observed in the high-z universe.  \nIf metal-enriched ( 𝑍 ≳ 10 -3 Z ⊙ ) dense star clusters form instead of SMSs after the onset of cold accretion, heavy seed BHs may form through a different channel. Our simulation predicts that the mass and radius of the cluster-forming clouds will be 𝑀 cloud ∼ 10 5 M ⊙ and 𝑅 cloud ∼ 0 . 1 pc. As discussed above, the star formation within such clouds is expected not to be impeded by both radiative feedback and SNe feedback, presumably resulting in very compact star clusters. Since Fukushima & Yajima (2021) predict that the SFE can reach 0.5 and the cluster radius is about 10% of the original cloud radius, the cluster mass and radius are estimated as 𝑀 cluster ≲ 0 . 5 𝑀 cloud ∼ 5 × 10 4 M ⊙ and 𝑅 cluster ∼ 0 . 1 𝑅 cloud ∼ 0 . 01 pc.  \nThe "runaway collision" resulting from the core collapse of a cluster is among the potential mechanisms for the formation of heavy seed BHs. This process takes place when the relaxation timescale 𝑡 relax is shorter than the stellar lifetime 𝑡 life (e.g. Portegies Zwart & McMillan 2002; Sakurai et al. 2017). Based on the previously mentioned estimates of the cluster mass and radius, we estimate the relaxation timescale as 𝑡 relax ∼ 3 × 10 -4 -3 × 10 -2 Myr, assuming that the mean stellar mass is in the range of 1 -10 2 M ⊙ . With 𝑡 relax ≪ 𝑡 life ∼ 10 Myr, runaway collision possibly occurs in the clusters we consider. In this case, however, the mass of the central object is believed to be limited to 0 . 1 -1% of the cluster (Portegies Zwart &McMillan 2002; Devecchi & Volonteri 2009; Sakurai et al. 2017; Reinoso et al. 2018). As a result, runaway collisions will produce BHs with masses around 𝑚 · ∼ 10 2 -10 3 M ⊙ , which are only slightly larger than light seed BHs originating from typical PopIII stars. Consequently, the runaway collision followed by the formation of star clusters is not an efficient mechanism for producing massive seed BHs.  \nEscala (2021) has proposed another process with dense stellar clusters, where a massive star can form by rapid stellar collision without going through two-body relaxation. If the cluster is sufficiently dense that the stellar collision timescale is shorter than the stellar lifetime 𝑡 coll < 𝑡 life , recurrent stellar mergers lead to the formation of a single massive star. Vergara et al. (2023, 2024) investigated this process by calculating the mass of the final massive object from various initial cluster masses and sizes using their toy model. Using their results, we expect that stellar clusters we consider marginally satisfy 𝑡 coll < 𝑡 life ∼ 10 Myr and estimate the mass of the most massive star as ∼ 0 . 1 -0 . 2 × 𝑀 cluster . This suggests that relatively heavy seed BHs with 𝑚 · ∼ 10 4 M ⊙ possibly form after the onset of cold accretion, even if metal enrichment induces the formation of star clusters instead of SMSs.  \nDavies et al. (2011) and Kroupa et al. (2020) have also proposed another channel, where the BH cluster, formed by the massive compact stellar cluster, collapses via efficient gravitational wave emission to provide a seed BH. In this channel, it is crucial that the velocity dispersion of the BH cluster achieves relativistic values such as ∼ 10% of the speed of light. This requires the radius of a cluster smaller than ∼ 10 -3 pc with a limited available gas mass of ≲ 10 5 M ⊙ around the halo centre. This radius is smaller than that of the clusters we consider, 𝑅 cluster ∼ 0 . 01 pc. If some process removes energy and angular momentum from the cluster, shrinking its size by a factor of ∼ 10, this channel can be realised. Gaete et al. (2024) predict that a BH with ∼ 10 3 M ⊙ will be produced from the BH cluster with ≲ 10 5 -6 M ⊙ . Since the BH cluster mass will be lower than or, at most, comparable to the gas cloud mass ∼ 10 5 M ⊙ , the seed BHs with their masses of ≲ 10 3 M ⊙ can be produced via this channel.', '5.4 Expected fate and evolution of supermassive stars': 'Previous stellar evolution calculations show that the fates of SMSs depend on the accretion rates (Umeda et al. 2016; Woods et al. 2017, 2020; Haemmerlé et al. 2018; Haemmerlé 2021; Herrington et al. 2023; Saio et al. 2024). Umeda et al. (2016) show that, in the case with / 𝑚 ★ ≥ 0 . 1 M ⊙ yr -1 , general relativistic (GR) pulsational instability induces stellar collapse (Chandrasekhar 1964) during or before Hburning. In the case with / 𝑚 ★ ≲ 10 -1 . 5 M ⊙ yr -1 , the collapse occurs before the GR instability becomes effective; immediately on core hydrogen exhaustion for / 𝑚 ★ = 10 -1 . 5 M ⊙ yr -1 , and after silicon burning for / 𝑚 ★ = 0 . 01 M ⊙ yr -1 . A resulting BH mass is, for example, 𝑚 · = 8 . 4 × 10 4 M ⊙ for / 𝑚 ★ = 10 -1 . 5 M ⊙ yr -1 (Woods et al. 2017; Haemmerlé et al. 2018).  \nIn our simulation, the mean accretion rates on SMSs are / 𝑚 ★ ≃ 0 . 04 M ⊙ yr -1 . According to the previous studies mentioned above, these stars collapse into BHs when they exceed 𝑚 · ≃ 8 × 10 4 M ⊙ . Fig. 3 shows that the most massive star exceeds 8 × 10 4 M ⊙ in mass at 𝑡 ≃ 225 Myr, with its accretion rate / 𝑚 ★ ≃ 0 . 04 M ⊙ yr -1 . The other three SMSs exceed this mass later at 𝑡 ≃ 226 . 5 Myr, before the gas disc is dispersed by photoevaporation. The seed BHs may take over the efficient accretion via cold accretion and the disc, if these gas structures persist for a while (also see Section 5.5). In that case, the SMS formation channel we studied provides a favourable environment for the subsequent growth of heavy seed BHs.  \nNagele et al. (2022) show that SMSs explode as very energetic SNe aided by the GR instability in a narrow mass range of 2 . 63 . 0 × 10 4 M ⊙ with no remnant BH, while their model does not include mass accretion onto the star (see also Chen et al. 2014). Most of the SMSs in our simulation exceed this mass range, and it is unlikely that such GR SNe occur unless some mechanism regulates the stellar mass into the narrow mass window. Note that, in our simulation, we merge two star particles when get closer than the sum of sink radii, thus reducing the number of stars. However, with higher resolution, star particles may form binary or multiple systems instead of merging (Greif et al. 2011; Prole et al. 2022; Kirihara et al. 2023), resulting in a typically lower stellar mass due to sharing of the disc mass. If this is the case, there may be more chances that GR SNe will occur.', '5.5 Succeeding BH growth in atomic cooling halos': 'As discussed in Section 5.4, some of the SMSs that form in our simulation can collapse into BHs in the middle of the evolution in reality. They are expected to remain within the dense disc, characterised by 𝑛 H ≳ 10 4 cm -3 and 𝑇 ∼ 10 4 K. We here discuss the possible subsequent evolution of gas accretion onto the seed BHs and whether the typical accretion rate for SMSs, / 𝑚 ≃ 0 . 04 M ⊙ yr -1 , is also available for seed BHs. If such rapid mass accretion is maintained for another ≳ 10 Myr, the BH grows to 𝑚 · ≳ 4 × 10 5 M ⊙ , which is favourable for later mass growth.  \nAccretion rates to BHs will be related to the mass transfer rate through the disc, which depends on the typical gas temperature 𝑇 (Eq. 17). In the case of the SMSs, the rapid accretion is achieved with a high disc temperature at 𝑇 ∼ 10 4 K. If H 2 molecular cooling is effective, the gas temperature decreases to 𝑇 < 10 3 K, resulting in an accretion rate that is too low for SMS formation. However, this is not the case because the stellar radiative feedback destroys H 2 molecules (Section 3.3). We investigate whether the radiative feedback from BHs also achieves a similar situation, assuming the luminosity of a seed BH as 𝐿 = min ( 𝐿 Edd , 0 . 1 / 𝑚 · 𝑐 2 ) and SED with two components of multi-color black body and non-thermal X-ray  \nwith power law (Kato et al. 1998; Kuhlen & Madau 2005; Jeon et al. 2014). If we assume the BH mass of 𝑚 ★ = 10 5 M ⊙ and accretion rate of / 𝑚 · = 10 -2 M ⊙ yr -1 , the resulting LW intensity is estimated as 𝐽 21 ∼ 2 × 10 3 ( 𝑟 / 𝑟 disc ) -2 . While this value is five orders of magnitude lower than the LW intensity for a SMS with 𝑇 eff = 3 × 10 4 K, it is still comparable to the critical intensity above which efficient H 2 molecular cooling is prevented. Therefore, we expect that the BH feedback also contributes to maintaining the high temperature of the disc, 𝑇 ∼ 10 4 K, with which rapid mass accretion should continue.  \nAs illustrated in Section 3.3.3, the growth of stellar mass is ultimately halted by the stellar ionising feedback, namely, the expansion of an HII bubble and the photoevaporation of the disc, in our fiducial simulation run. We also investigate whether similar processes limit the mass increase of seed BHs, if these are present in the disc. Using the same BH feedback model, we calculate the emissivity of ionising photons from an accreting seed BH as Φ = 8 × 10 52 s -1 , which is somewhat lower than that for SMSs with 𝑇 eff = 3 × 10 4 or 10 5 K. Given the comparable strength of ionising feedback from BHs, they would be capable of accreting about ∼ 10 5 M ⊙ of gas, similar to many of the SMSs in our fiducial run.  \nIn Section 3.3.3, we have also shown that the ionising feedback of SMSs becomes significant once initiated. When an HII bubble expands within the disc, driven by gravitational three-body interactions and the scattering of a ≃ 5000 M ⊙ star, it diminishes the accretion rates of other SMSs. These SMSs, which experience slow accretion, then become powerful ionising sources with 𝑇 eff = 10 5 K. The ionising feedback intensifies as it continues. We expect that such a drastic evolution is mitigated if some SMSs are replaced by seed BHs. When an accretion rate to a BH drops, its ionising emissivity correspondingly reduces. The destruction of the disc might occur somewhat slowly, even if it occurs. Therefore, the seed BHs might have the potential to grow in mass over several Myr following their formation, which will be confirmed in future simulations.', '5.6 Statistics and seed BH density': "This section discusses the number density of seed BH possibly formed by the channel investigated in this work. In Section 5.3 in K23, we already estimated the seed BH number density as  \n𝑛 BH = 𝑛 ACH 𝑓 p d 𝑁 merger d 𝑡 Δ 𝑡 (19) ∼ 10 -9 cMpc -3 GLYPH<18> 𝑛 ACH 10 -3 cMpc -3 GLYPH<19> GLYPH<18> 𝑓 p 2 × 10 -4 GLYPH<19> GLYPH<18> Δ 𝑡 3 Myr GLYPH<19> ,  \nwhere 𝑛 ACH is the number density of ACHs, 𝑓 p the metal-pristine fraction, i.e. the fraction of ACHs which have experienced no prior star formation, d 𝑁 merger / d 𝑡 is the merger rate of 𝑀 halo ∼ 10 7 M ⊙ with 𝑀 ' halo ≳ 10 6 M ⊙ (Fakhouri et al. 2010), and Δ 𝑡 the time delay until the onset of cold accretion from the first cloud collapse in a given ACH. Since Li et al. (2021) show that ACHs at 𝑧 ≃ 20 -40 gorw to become host halos of SMBH at 𝑧 ≳ 6, we evaluate 𝑛 ACH and d 𝑁 merger / d 𝑡 at 𝑧 = 30 in Equation (19). Since the metal-pristine fraction at this epoch is unknown, we instead use the value at 𝑧 ≃ 10 provided by Fernandez et al. (2014). Equation (19) predicts 𝑛 BH comparable to the SMBH number density at 𝑧 ≳ 6, ∼ 10 -9 cMpc -3 .  \nThe term of ( d 𝑁 merger / d 𝑡 ) Δ 𝑡 considers the condition as follows; if cold accretion occurs too late compared to the stellar lifetime, typically 𝑡 life = 3 Myr, the gas in the halo will be metal-enriched by SNe before the emergence of cold accretion, inhibiting SMS formation. In Section 5.2, we have also discussed that our current simulations where SNe feedback is ignored will approximate such  \ncases with Δ 𝑡 < 3 Myr. Equation (19) suggests that even if these cases are rare, their number density will be large enough to explain the origin of SMBHs at 𝑧 ≳ 6. We also point out that the radiative feedback effect somewhat prolongs Δ 𝑡 , as show in Section 3.3.1 (see Figures 4 - 6). In this paper, we only study a specific ACH where Δ 𝑡 ∼ 10 Myr without incorporating the radiative feedback. It is uncertain how the feedback effect extends the delay for rare ACHs which originally had much shorter durations, such as Δ 𝑡 ∼ Myr. This remains a task for future work.  \nNote that our estimate with Equation (19) only counts the number density of ACHs without considering the possible mass growth histories of seed BHs up to the epoch of 𝑧 ≃ 6. Since we have shown that the sequential formation of multiple SMSs can be induced by the onset of cold accretion (Section 3.3.2), a very massive seed BH with ∼ 10 6 M ⊙ may be provided if the BHs with ∼ 10 5 M ⊙ produced by the SMSs merge efficiently. If it occurs, this is advantageous for subsequent growth of the BH mass. With such a very massive seed BH, we may consider the ACHs at a relatively late starting epoch of BH growth at 𝑧 ≃ 20 instead of 𝑧 = 30 in Equation (19). Note that the number density of ACHs increases as 𝑛 ACH ∼ 0 . 1 cMpc -3 at 𝑧 = 20, which is two orders of magnitude higher than the value used in Equation (19) (Barkana & Loeb 2001).  \nConsidering the lower redshifts may be disfavored in terms of the metal-pristine fraction 𝑓 p. However, recall that the value referred to in Equation (19) is an estimate for 𝑧 ≃ 10, which will be a lower limit for the higher redshift 𝑧 ≃ 20. Furthermore, as discussed in Section 5.3, the formation of SMS may still take place with non-zero metallicities provided 𝑍 < 𝑍 crit ∼ 10 -4 -10 -3 Z ⊙ . Allowing prior star formation within this condition, we further amplify 𝑛 BH by some orders of magnitude.", '6 CONCLUSIONS': "We have investigated the formation of Pop III stars within an ACH ( 𝑀 halo ∼ 10 7 M ⊙ ) at a redshift of 𝑧 ≃ 17, when the cold accretion first appears in the early universe, performing cosmological radiation hydrodynamics simulations. Following up our previous work K23, we have explored the radiative feedback effects caused by stars forming within the same halo. Specifically, we examined the potential formation of SMS triggered by cold accretion under the influence of stellar radiative feedback. We have used zoom-in and particlesplitting techniques to achieve high spatial resolution. We have followed a long-term (several Myr) evolution of star formation in a dense gas disc fed by cold accretion. Our findings are summarised as follows.  \n- · As studied in K23, the formation of Pop III stars first occurs at the epoch of 𝑧 ≃ 18 . 9, when the virial temperature of the halo is 𝑇 vir ≃ 10 4 K, before the emergence of cold accretion. The radiative feedback from such a 'preceding star' affects the evolution even in this early stage. An HII region develops at the halo centre and temporarily halts further star formation. This effect delays the emergence of the cold accretion by ≃ 20 Myr. However, cold accretion emerges when the virial mass of the halo exceeds 𝑀 halo ∼ 10 7 M ⊙ at 𝑧 ∼ 16 . 9. The large ram pressure of the cold accretion streams overcomes the thermal pressure of the photoionised gas, which finally quenches the HII bubble.\n- · Immediately after cold accretion emerges, the star formation is induced owing to the rapid mass supply to the halo centre. This begins with cloud collapse induced by H 2 molecular cooling. Consequently, massive Pop III stars with ∼ 10 3 M ⊙ emerge, and cold accretion persists. A dense gas disc forms around these stars in the  \ncentral region within 𝑟 ≲ 10 -2 𝑟 vir , rapidly accreting gas to increase its mass. The average accretion rates for individual stars are a few × 0 . 01 M ⊙ yr -1 , allowing them to grow into SMSs of 𝑚 ★ ≃ 10 5 M ⊙ within a few Myr. As these rapidly accreting stars develop, a significant portion of the dense disc remains atomic. While the disc undergoes fragmentation to create multiple stars, the efficient mass growth continues. Only several stars form, but most of them grow into SMSs over the timescale of ∼ Myr. The stellar photoionising feedback is ineffective in halting the mass accretion during that.  \n- · The substantial accretion rate observed for each star can be attributed to the following two effects. Firstly, cold accretion delivers gas directly to the halo centre at rates ranging from 0 . 04 -0 . 1 M ⊙ yr -1 , forming a central dense disc which sustains the stars with a nearly constant accretion rate. Secondly, radiative feedback efficiently destroys H 2 molecules and H -ions, maintaining the atomic nature of the disc. The feedback from ionising (EUV) photons is relatively weak because rapidly accreting stars possess large radii and low effective temperatures of 𝑇 eff = 5 × 10 3 K or 3 × 10 4 K. Additionally, the dense disc formed by cold accretion shields the stellar ionising photons, preventing the expansion of an HII region. The atomic disc does not fragment efficiently, resulting in the formation of only a few stars. Thus, the gas provided by cold accretion is not competed for by many stars, allowing a high accretion rate for each individual star.\n- · The ionising feedback eventually operates after ≃ 3 Myr from the SMS formation induced by the cold accretion. Multiple SMSs emerge due to the fragmentation of the central gas disc by this epoch and gravitational interactions among them eject one star outside the disc. The ejected star cannot sustain a high accretion rate since the surrounding density outside the disc is very low. The accretion rate falls below 0 . 004 M ⊙ yr -1 and the star contracts to reach a high effective temperature of 𝑇 eff ≃ 10 5 K. The stellar emissivity of the ionising photons increases significantly, leading to the expansion of an HII region. The strong feedback disrupts the central disc and halts the further growth of stellar mass. The HII region continues to grow, primarily along the polar directions, forming a large bipolar bubble. This bubble expands beyond the virial radius over ∼ Myr, causing the photoevaporation of the entire halo gas. The cold streams toward the halo centre are finally disturbed by this effect.  \nIn summary, we conclude that the interplay between cold accretion and stellar radiative feedback facilitates the sequential formation of SMSs. The heavy seed BHs formed by this process have the great advantage for further growth that they are embedded in the dense, massive gas disc at the halo centre, and the disc is continuously fed by cold accretion. The next key to understanding SMBH formation is to study how efficiently they subsequently grow in mass over the cosmological timescale.", 'ACKNOWLEDGEMENTS': 'The authors express their cordial gratitude to Prof. Takahiro Tanaka for his continuous interest and encouragement. We sincerely appreciate Kazuyuki Omukai, Naoki Yoshida, Zoltán Haiman, Avishai Dekel, Daniel Schaerer, Kohei Inayoshi, Koutarou Kyutoku, Kazuyuki Sugimura, Shingo Hirano, Gen Chiaki, Daisuke Toyouchi, Ryoki Matsukoba, Devesh Nandal, Kazutaka Kimura, Yohsuke Enomoto, and Tomoya Suzuguchi for the fruitful discussions and comments. We sincerely appreciate Volker Springel for the development of the simulation code GADGET-3 we make use of, which is essential for our calculations. The numerical simulations were carried out on XC50 Aterui II at the Center for Computational Astro-  \nhysics (CfCA) of the National Astronomical Observatory of Japan. This research could never be accomplished without the support by Grants-in-Aid for Scientific Research (TH:19H01934, 19KK0353, 22H00149) from the Japan Society for the Promotion of Science and JST SPRING, Grant Number JPMJSP2110. We use the SPH visualisation tool SPLASH (Price 2007, 2011) in Figs 2, 4, 8, 10, 11, and 12.', 'DATA AVAILABILITY': 'The data underlying this article will be shared on a reasonable request of the corresponding author.', 'REFERENCES': 'record ascl:1103.004 (ascl:1103.004)  \n- Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., 2022, MNRAS, 510, 4019', 'APPENDIX A: SIMULATION WITH DIFFERENT FEEDBACK EFFICIENCY 𝐿 P = 0 . 1 𝐿 Edd': 'In our fiducial simulation run, we have placed a star particle whose mass is ≃ 3000 M ⊙ as a result of the initial cloud collapse prior to the onset of the cold accretion. This mass should be considered an upper limit due to our limited spatial resolution, as described in Section 2.3. We have used 𝜂 = 0 . 03 in Eq. 1), which is approximated by, if a single star is assumed, a star with 𝑚 ★ ≃ 300 M ⊙ . While this may be a  \ntypical stellar mass of normal Pop III stars, it is important to consider another case with different feedback efficiency 𝜂 . For comparison, here we describe the case with 𝜂 = 0 . 1, assuming a more massive star or cluster. Note that we allow the star to emit radiation for more than ∼ 10 Myr, which may be beyond the lifetime of very massive stars (see Section 5.2 for a relevant discussion).  \nFig. A1 shows the time evolution of the gas structures at the halo scale, in the same manner as in Fig. 4. In this case, an HII bubble breaks out within ≃ 5 -10 Myr, and the gas within the halo is fully ionised and heated at 𝑇 > 10 4 K. A large amount of the gas is expelled from the halo due to photoevaporation on a halo scale. The lower panels indicate that filamentary streams remain within the virial radius, although these structures will eventually disperse as well. While not explored here, there will be considerable delay before cold accretion begins, which will occur once the halo mass substantially surpasses the ACH range. Even in such a situation, the formation of SMS induced by cold accretion might occur anyway, which should be confirmed in future simulations.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}