JWST has brought us new insights into Cosmic Dawn with tentative detection of the unique signatures of metalfree Population III Pop III stars such as strong He II emission extremely blue ultraviolet spectrum and enhanced nitrogen abundance. Selfconsistent theoretical predictions of the formation rates sites and masses of Pop III stars are crucial for interpreting the observations but are challenging due to complex physical processes operating over the large range of lengthscales involved. One solution is to combine analytical models for the smallscale star formation process with cosmological simulations that capture the largescale physics such as structure formation radiation backgrounds and baryondark matter streaming motion that regulate the conditions of Pop III star formation. We build an analytical model to predict the final masses of Pop III starsclusters from the properties of starforming clouds based on the key results of smallscale star formation simulations and stellar evolution models. Our model for the first time considers the interplay between feedback and fragmentation and covers different modes of Pop III star formation ranging from ordinary small inlineformulatexmath idTM0001 notationLaTeXsim102000 rm Modottexmathinlineformula clusters in molecularcooling clouds to massive inlineformulatexmath idTM0002 notationLaTeXgtrsim104 rm Modottexmathinlineformula clusters containing supermassive inlineformulatexmath idTM0003 notationLaTeXsim1043times 105 rm Modottexmathinlineformula stars under violent collapse of atomiccooling clouds with large gas accretion rates of inlineformulatexmath idTM0004 notationLaTeXgtrsim0.1 rm Modot  yr1texmathinlineformula. As an example the model is applied to the Pop III starforming clouds in the progenitors of typical haloes hosting highz luminous quasars inlineformulatexmath idTM0006 notationLaTeXMrm hsim 1012 rm Modottexmathinlineformula at inlineformulatexmath idTM0007 notationLaTeXzsim 6texmathinlineformula which shows that formation of Pop III massive clusters is common inlineformulatexmath idTM0008 notationLaTeXsim2070 rm per centtexmathinlineformula in such biased inlineformulatexmath idTM0009 notationLaTeXsim4sigmatexmathinlineformula regions and the resulting heavy black hole seeds from supermassive stars can account for a significant fraction of observed luminous inlineformulatexmath idTM0010 notationLaTeXgtrsim1046 rm erg s1texmathinlineformula quasars at inlineformulatexmath idTM0011 notationLaTeXzsim 6texmathinlineformula.
['Astrophysics - Astrophysics of Galaxies']
Towards a universal analytical model for Population III star formation interplay between feedback and fragmentation
{'Towards a universal analytical model for Population III star formation: interplay between feedback and fragmentation': "Boyuan Liu ★ 1 , 2 , James Gurian 3 , Kohei Inayoshi 4 , Shingo Hirano 5 , Takashi Hosokawa 6  \n, Volker Bromm 7 , 8 , and Naoki Yoshida 9 , 10 , 11  \n- 1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 2 Institut für Theoretische Astrophysik, Zentrum für Astronomie, Universität Heidelberg, Albert Ueberle Straße 2, D-69120 Heidelberg, Germany\n- 3 Perimeter Institute for Theoretical Physics, Waterloo, Ontario, N2L 2Y5, Canada\n- 4 Kavli Institute for Astronomy and Astrophysics, Peking University, Beĳing 100871, People's Republic of China\n- 5 Department of Astronomy, School of Science, University of Tokyo, Tokyo 113-0033, Japan\n- 6 Department of Physics, Kyoto University, Sakyo, Kyoto 606-8502, Japan\n- 7 Department of Astronomy, University of Texas, Austin, TX 78712, USA\n- 8 Weinberg Institute for Theoretical Physics, University of Texas, Austin, TX 78712, USA\n- 9 Department of Physics, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan\n- 10 Research Center for the Early Universe, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan\n- 11 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'JWSThasbrought us new insights into Cosmic Dawn with tentative detection of the unique signatures of metal-free Population III (Pop III) stars, such as strong HeII emission, extremely blue UV spectrum, and enhanced nitrogen abundance. Self-consistent theoretical predictions of the formation rates, sites, and masses of Pop III stars are crucial for interpreting the observations, but are challenging due to complex physical processes operating over the large range of length scales involved. One solution is to combine analytical models for the small-scale star formation process with cosmological simulations that capture the large-scale physics such as structure formation, radiation backgrounds, and baryon-dark matter streaming motion that regulate the conditions of Pop III star formation. We build an analytical model to predict the final masses of Pop III stars/clusters from the properties of star-forming clouds, based on the key results of small-scale star formation simulations and stellar evolution models. Our model for the first time considers the interplay between feedback and fragmentation and covers different modes of Pop III star formation ranging from ordinary small ( ∼ 10 -2000 M ⊙ ) clusters in molecular-cooling clouds to massive ( ≳ 10 4 M ⊙ ) clusters containing supermassive ( ∼ 10 4 -3 × 10 5 M ⊙ ) stars under violent collapse of atomic-cooling clouds with large gas accretion rates of ≳ 0 . 1 M ⊙ yr -1 . As an example, the model is applied to the Pop III star-forming clouds in the progenitors of typical haloes hosting high𝑧 luminous quasars ( 𝑀 h ∼ 10 12 M ⊙ at 𝑧 ∼ 6), which shows that formation of Pop III massive clusters is common ( ∼ 20 -70%) in such biased ( ∼ 4 𝜎 ) regions, and the resulting heavy black hole seeds from supermassive stars can account for a significant fraction of observed luminous ( ≳ 10 46 erg s -1 ) quasars at 𝑧 ∼ 6.  \nKey words: stars: Population III - stars: black holes - dark ages, reionization, first stars', '1 INTRODUCTION': "The first generation of so-called Population III (Pop III) stars, formed in metal-free primordial gas, are expected to be fundamentally different from present-day metal-enriched Population I/II (Pop I/II) stars (reviewed by, e.g., Bromm et al. 2009; Bromm 2013; Loeb & Furlanetto 2013; Haemmerlé et al. 2020; Klessen & Glover 2023). Due to their pristine formation environments, Pop III stars tend to form in small clusters with a broad top-heavy initial mass function (IMF) extending up to a few ∼ 10 5 M ⊙ according to small-scale (radiative) (magneto-) hydrodynamic simulations (e.g. Clark et al. 2011; Greif et al. 2011, 2012; Susa et al. 2014; Hirano et al. 2014, 2015, 2018;  \nStacy et al. 2013; Stacy & Bromm 2014; Stacy et al. 2016; Hosokawa et al. 2016; Hirano & Bromm 2017; Hirano et al. 2017; Suazo et al. 2019; Susa 2019; Sugimura et al. 2020, 2023; McKee et al. 2020; Wollenberg et al. 2020; Park et al. 2021, 2023, 2024; Chon et al. 2018, 2021, 2022; Chon & Omukai 2020; Sakurai et al. 2020; Regan et al. 2014, 2017, 2020a; Sharda et al. 2020, 2021; Latif et al. 2021, 2022b; Woods et al. 2021, 2024; Riaz et al. 2018, 2023a,b; Prole et al. 2022a,b, 2023, 2024; Regan 2023; Reinoso et al. 2023; Toyouchi et al. 2023; Kiyuna et al. 2024; Sharda & Menon 2024; Sadanari et al. 2024). Further, stellar evolution models predict that Pop III stars are more compact and hotter than Pop I/II stars, likely fast rotating, have negligible mass loss from metal line-driven winds (e.g. Schaerer 2002; Meynet et al. 2006; Ekström et al. 2008; Heger & Woosley 2010; Yoon et al. 2012; Tanikawa et al. 2020; Farrell  \net al. 2021; Murphy et al. 2021; Aryan et al. 2023; Martinet et al. 2023; Nandal et al. 2023; Volpato et al. 2023), and end their lives in diverse supernova events with distinct metal yields or collapse directly into massive back holes (BHs) (e.g., Heger & Woosley 2002; Heger et al. 2003; Maeda & Nomoto 2003; Umeda & Nomoto 2003, 2005; Iwamoto et al. 2005; Kobayashi et al. 2006; Tominaga 2009; Heger & Woosley 2010; Nomoto et al. 2013).  \nWith these peculiar properties, Pop III stars play important roles in the first billion years of cosmic history (i.e., Cosmic Dawn) through their radiation, metal enrichment, and by seeding BHs, which produce unique signatures in direct observations (e.g., Windhorst et al. 2018; Grisdale et al. 2021; Nakajima & Maiolino 2022; Vikaeus et al. 2022; Trussler et al. 2023; Katz et al. 2023; Larkin et al. 2023; Venditti et al. 2024a; Zackrisson et al. 2011, 2012; Zackrisson et al. 2024), imprints in cosmic chemical and thermal evolution (see, e.g., Karlsson et al. 2013; Nomoto et al. 2013; Frebel & Norris 2015; Barkana 2016; Dayal & Ferrara 2018, for reviews), and impact on the formation and evolution of subsequent populations of stars, galaxies and supermassive BHs (SMBHs) (see, e.g., Bromm & Yoshida 2011; Johnson et al. 2013; Pawlik et al. 2013; Jeon et al. 2015; Smith et al. 2015; Regan et al. 2017; Sakurai et al. 2017; Haemmerlé et al. 2020; Inayoshi et al. 2020; Schauer et al. 2021; Volonteri et al. 2021; Chon et al. 2022; Sarmento & Scannapieco 2022; Chiaki et al. 2023; Sanati et al. 2023; Reinoso et al. 2023; Wise 2023; Kiyuna et al. 2024; Regan & Volonteri 2024; Rossi et al. 2024), allowing us to constrain their properties through direct and indirect observations.  \nRecently, JWST has revealed the Cosmic Dawn with unprecedented discoveries of puzzling phenomena unexpected from conventional theoretical models (Adamo et al. 2024), some of which can be explained by the unique signatures of Pop III stars and their BH remnants, such as metal-poor strong HeII-emission systems (Wang et al. 2024; Vanzella et al. 2023; Maiolino et al. 2024b), galaxies with extremely blue UV spectra (e.g., Nanayakkara et al. 2023; Austin et al. 2024; Cullen et al. 2024), strong carbon and nitrogen enrichment (D'Eugenio et al. 2023; Cameron et al. 2023; Senchyna et al. 2024; Ji et al. 2024; Schaerer et al. 2024; Topping et al. 2024) likely from fast-rotating or very massive Pop III stars (Liu et al. 2021b; Nandal et al. 2024d,c,b; Tsiatsiou et al. 2024), highly magnified individual massive stars (Welch et al. 2022; Schauer et al. 2022; Zackrisson et al. 2024), a slowly evolving UV luminosity function hinting at high UV variability, non-constant star formation efficiency (SFE), extreme nebular emission powered by hot ionizing sources, and top-heavy IMF (e.g., Inayoshi et al. 2022b; Finkelstein et al. 2024; Shen et al. 2023; Conselice et al. 2024; Cueto et al. 2024; Harikane et al. 2024; Katz et al. 2024; Trinca et al. 2024; Ventura et al. 2024), and overabundant luminous quasars with overmassive (rapidly spinning) BHs at 𝑧 ≳ 5 (e.g., Goulding et al. 2023; Harikane et al. 2023; Maiolino et al. 2023, 2024a; Kokorev et al. 2023; Übler et al. 2023; Akins et al. 2024; Durodola et al. 2024; Greene et al. 2024; Juodžbalis et al. 2024; Kocevski et al. 2024; Matthee et al. 2024; Inayoshi & Ichikawa 2024) likely originating from heavy BH seeds ( ≳ 10 4 -5 M ⊙ ). Although there is still no conclusive direct detection of Pop III-dominated galaxies, which is likely challenging for JWST (Gardner et al. 2006; Riaz et al. 2022; Nakajima & Maiolino 2022; Katz et al. 2023; Bovill et al. 2024) and may require large extraterrestrial telescopes (Angel et al. 2008; Rhodes et al. 2020; Schauer et al. 2020), these observations highlight the importance of the Pop III components in high𝑧 galaxies. As JWST continues its revolutionary observations, other facilities (e.g., Roman Space Telescope and Vera C. Rubin Observatory) will come online in the near future, which may bring us new detection of Pop III signatures, such as the transient signals of Pop III gamma-ray bursts and supernova  \nexplosions (e.g., Fryer et al. 2022; Lazar & Bromm 2022; Hartwig et al. 2023; Venditti et al. 2024b; Wiggins et al. 2024).  \nGravitational wave observations of massive binary BH (BBH) mergers have also tentatively discovered a non-negligible contribution from Pop III BH remnants (up to 10% of the detected events, see, e.g., Kinugawa et al. 2020, 2021b; Iwaya et al. 2023), especially for the most massive mergers such as GW190521 with unusual BH masses (Abbott et al. 2020b; Abbott et al. 2020a) in the range ∼ 50 -130 M ⊙ mostly forbidden for BHs from Pop I/II stars by standard pair-instability supernova (PISN) models (e.g., Heger et al. 2003; Belczynski et al. 2016; Yoshida et al. 2016; Woosley 2017; Spera & Mapelli 2017; Marchant et al. 2019; Marchant & Moriya 2020; Mapelli et al. 2020). Thanks to their massive, compact nature and lack of mass loss, Pop III stars can produce BHs up to ∼ 100 M ⊙ before triggering PISNe and more massive ( ≳ 130 M ⊙ ) BHs via direct collapse (e.g., Farrell et al. 2021; Volpato et al. 2023; Santoliquido et al. 2023; Mestichelli et al. 2024), which can end up in BBH mergers like GW190521 via isolated binary stellar evolution under proper conditions (Kinugawa et al. 2021a; Tanikawa et al. 2021, 2022; Santoliquido et al. 2023; Tanikawa 2024) or through dynamical interactions in massive ( ≳ 10 4 M ⊙ ) dense clusters of Pop III stars/BHs themselves (Wang et al. 2022; Liu et al. 2023; Mestichelli et al. 2024) and (Pop I/II) nuclear star clusters (Liu & Bromm 2020c, 2021; Liu et al. 2024). In the next decades, the 3rd-generation GW detectors such as the Einstein Telescope (ET, Punturo et al. 2010; Maggiore et al. 2020) and the Cosmic Explorer (Reitze et al. 2019; Evans et al. 2023) will be able to detect BBH mergers up to 𝑧 ∼ 30, among which the contribution of Pop III BH remnants are expected to be larger and better characterized (e.g., Tanikawa et al. 2022; Franciolini et al. 2024; Santoliquido et al. 2024).  \nTo fully understand the roles played by Pop III stars at Cosmic Dawn and the underlying physics (of early star formation, stellar feedback, BH formation and growth) through observations, selfconsistent theoretical predictions are required not only for  \n- (1) the formation rates and sites of Pop III stars, but also for\n- (2) the masses of Pop III star clusters and the stellar IMF.  \nEvidently, the total mass of Pop III stars in a galaxy is vital for the detectability of Pop III signatures embedded in the signals from coexisting Pop I/II stars and AGN, which are common in the galaxies detectable by JWST (e.g., Sarmento et al. 2018; Liu & Bromm 2020b; Riaz et al. 2022; Venditti et al. 2023, 2024a). The cluster mass also determines whether dynamical interactions can efficiently produce (distinguishable) massive, eccentric Pop III BBH mergers (Wang et al. 2022; Liu et al. 2023; Mestichelli et al. 2024), which may not be possible for isolated binary stellar evolution (Costa et al. 2023; Santoliquido et al. 2023). The Pop III IMF affects the UV spectra of Pop III clusters (e.g., Zackrisson et al. 2011, 2012; Bovill et al. 2024) and the X-ray outputs of Pop III X-ray binaries (Sartorio et al. 2023), shaping the imprints of Pop III stars on the 21-cm signal and reionization (e.g., Fialkov 2014; Gessey-Jones et al. 2022, 2024; Salvador-Solé et al. 2017, 2022; Fialkov et al. 2023). Moreover, the shape of the Pop III IMF in the high-mass regime ( ≳ 100 M ⊙ ) has essential implications on the rates/abundances of strongly lensed massive Pop III stars detectable by JWST (Zackrisson et al. 2024), Pop III PISNe (e.g., Lazar & Bromm 2022; Wiggins et al. 2024; Venditti et al. 2024b), very massive stars (VMSs, ∼ 10 3 -10 4 M ⊙ ) with peculiar metal enrichment (e.g., Nandal et al. 2024c,b), and supermassive stars (SMSs, ≳ 10 4 M ⊙ ) as progenitors of heavy BH seeds (Smith & Bromm 2019; Inayoshi et al. 2020; Haemmerlé et al. 2020; Volonteri et al. 2021; Wise 2023; Regan & Volonteri 2024) that are likely needed to explain the SMBHs in high𝑧 luminous  \nquasars (e.g., Fan et al. 2001; Willott et al. 2010; Matsuoka et al. 2016; Onoue et al. 2019; Shen et al. 2020; Yang et al. 2021; Fan et al. 2023; Akins et al. 2024; Durodola et al. 2024; Greene et al. 2024; Matthee et al. 2024; Kocevski et al. 2024).  \nThese two aspects are expected to be closely related to each other in the rapidly evolving environments of Cosmic Dawn. However, they are usually modelled separately in current theoretical studies, as it is challenging to combine them in one theoretical framework due to the large range of scales involved (from over-Mpc to sub-AU). The first aspect is mostly investigated with large-scale cosmological hydrodynamic simulations (e.g., Johnson et al. 2013; Smith et al. 2015; Xu et al. 2016; Sarmento et al. 2017, 2018; El-Badry et al. 2018; Liu et al. 2019; Liu & Bromm 2020a,b; Liu et al. 2022; Skinner & Wise 2020; Kulkarni et al. 2021; Schauer et al. 2019a, 2021; Kulkarni et al. 2022; Yajima et al. 2022, 2023; Kiyuna et al. 2023; Garcia et al. 2023; Venditti et al. 2023; Incatasciato et al. 2023; Correa Magnus et al. 2024; Lenoble et al. 2024; Sugimura et al. 2024; Smith et al. 2024) and semi-analytical models (e.g., Manrique et al. 2015; SalvadorSolé et al. 2017; Griffen et al. 2018; Dayal et al. 2020; Visbal et al. 2020; Li et al. 2021; Lupi et al. 2021; Hartwig et al. 2022; Trinca et al. 2022, 2024; Hegde & Furlanetto 2023; Nebrin et al. 2023; Bovill et al. 2024; Ventura et al. 2024; Feathers et al. 2024), which show that the formation rates and sites of Pop III star formation are regulated by cosmic/halo-scale processes, such as cosmic structure formation (for different dark matter models), dissociation of H 2 by Lyman-Werner (LW) radiation, streaming motion between baryons and dark matter, ionization and heating by X-rays, dynamical heating by halo mergers, and metal enrichment. For instance, LW radiation, streaming motion, and halo mergers can delay Pop III star formation, even shifting the formation sites to more massive ( ≳ 10 8 M ⊙ ) atomic-cooling haloes (e.g. Yoshida et al. 2003; Tseliakhovich & Hirata 2010; Tseliakhovich et al. 2011; Fialkov et al. 2012; Tanaka & Li 2014; Fernandez et al. 2014; Schauer et al. 2017; Inayoshi et al. 2018; Regan et al. 2020b; Kulkarni et al. 2021; Schauer et al. 2021). Such large-scale models lack the resolution to follow the sub-pc scale star formation process in detail (which is computationally prohibitive in a cosmologically representative volume). As a result, they estimate the final outcomes of Pop III star-formation with simple assumptions on the SFE and IMF treated as independent global parameters. On the other hand, for the second aspect, small-scale high-resolution (idealized/zoomin) simulations of star-forming clouds (see above for references) have revealed that the final mass of Pop III stars formed and their IMF are determined by the complex interplay between fragmentation, radiative feedback, and stellar evolution, which have non-trivial correlations with the (initial) properties of the star-forming cloud. For instance, the standard pathway of Pop III star formation in molecularcooling minihaloes ( 𝑀 h ∼ 10 5 -6 M ⊙ at 𝑧 ∼ 20 -30 with gas infall rates ∼ 10 -4 -0 . 01 M ⊙ yr -1 ) typically produces Pop III stars with masses ∼ 10 -1000 M ⊙ regulated by radiative feedback (e.g., Hirano et al. 2014, 2015; Park et al. 2021; Sugimura et al. 2023), while more violent collapse ( ≳ 0 . 1 M ⊙ yr -1 ) in atomic-cooling haloes can overcome/suppress radiative feedback and produce SMSs up to a few 10 5 M ⊙ (e.g., Bromm & Loeb 2003; Chon et al. 2018; Chon & Omukai 2020; Sakurai et al. 2020; Toyouchi et al. 2023; Reinoso et al. 2023; Regan 2023; Kiyuna et al. 2024).  \nTo bridge the aforementioned insights from theoretical efforts on the (1) large- and (2) small-scale aspects and achieve self-consistent modelling of Pop III star formation across all scales, we build a universal analytical model to predict the outcomes of Pop III star formation from the properties of star-forming clouds, as described in Sec. 2. Our model inherits the basic elements from previous analytical studies (e.g., Omukai & Inutsuka 2002; McKee & Tan 2002,  \n2003, 2008; Tan & McKee 2004; Hosokawa & Omukai 2009; Johnson et al. 2012; Yajima & Khochfar 2017; Fukushima et al. 2018; Li et al. 2021; Toyouchi et al. 2023) such as the balance between mass inflow from the collapsing cloud and mass loss by feedback, the reaction of protostars to accretion, as well as their finite lifetimes and maximum masses set by stellar evolution (Sec. 2.2). It captures the key trends in the evolution of Pop III protostar systems inferred from small-scale simulations with simple scaling relations governed by physically motivated parameters (Sec. 2.1). Furthermore, our model can track the growth and feedback of multiple protostars, and is used to explore the interplay between feedback and fragmentation in determining the formation efficiency and IMF of Pop III stars (Sec.2.3, Appendices A and B). As an example, the model is applied to the Pop III star-forming clouds in the progenitor haloes of high𝑧 luminous quasar host galaxies (Li et al. 2021) to evaluate the mass distributions of Pop III stars/clusters formed in such over-dense regions (Sec. 3.1). We also discuss the caveats in our model and possible ways to improve it (Sec. 3.2). Given its simplicity and flexibility, our model can be easily incorporated into large-scale semi-analytical models and cosmological simulations.", '2 ANALYTICAL MODEL': 'In this section, we describe our analytical model for Pop III star formation. We first introduce the scaling relations learned from simulations to capture the general evolution of Pop III star-forming discs and protostar systems governed by gravity and hydrodynamics (Sec. 2.1). Next, we consider the roles played by feedback and stellar evolution, which are combined with the scaling relations to calculate the final mass under the simple assumption that only a single star forms per cloud (Sec. 2.2). Finally, we generalize the model to consider multiple protostars and explore the effects of fragmentation/multiplicity on the final stellar mass (Sec. 2.3). For conciseness, we only show the results for select examples below. To demonstrate the flexibility of our model, we provide more detailed results for a broader exploration of parameter space and an alternative parametrisation of multiplicity in Appendices A and B. The qualitative conclusions shown in this section remain valid for other choices of parameters/assumptions, although the predicted final stellar mass (for the same gas infall rate) can vary up to an order of magnitude.', '2.1 Scaling relations for Pop III star-forming discs and protostar systems': "It is shown in Liu et al. (2021a, see their sec. 2.1) that the total masses 𝑀 and radii 𝑅 of Pop III protostar systems in hydrodynamic simulations of typical primordial star-forming clouds in minihaloes ( 𝑀 h ∼ 10 5 -10 6 M ⊙ , 𝑧 ∼ 20 -30) can be reproduced with an analytical solution within a factor of 3 at least in the early stage (see also fig. 9 in Susa 2019):  \n𝑀 ≃ 400 M ⊙ h 𝑡 / GLYPH<16> 10 5 yr GLYPH<17> i 𝛽 , (1) 𝑅 ≃ AU ( 𝑡 / yr ) 𝛿 , (2)  \nwhere 𝛽 = 4 -3 𝛾 eff = 0 . 73, 𝛿 = 2 -𝛾 eff = 0 . 91, given 𝛾 eff = 1 . 09 as the polytropic index in the effective equation of state (EoS) 𝑃 ∝ 𝑛 𝛾 eff of collapsing primordial gas under canonical H 2 cooling 1 in the  \ntemperature range 𝑇 ∼ 300 -1000 K (Omukai & Nishi 1998; Omukai et al. 2005). To generalize this solution for Pop III star formation in smaller/larger haloes with weaker/stronger gas inflows, we re-write Eq. 1 with a normalisation factor 𝑀 0 and a characteristic timescale 𝑡 0 , and express 𝑅 with 𝑀 :  \n𝑀 = 𝑀 0 ( 𝑡 / 𝑡 0 ) 𝛽 , (3)  \n𝑅 = h 10 5 𝛽 𝑀 /( 400 M ⊙ ) i 𝛿 / 𝛽 AU ∝ 𝑀 𝛿 / 𝛽 ≃ 𝑀 1 . 25 . (4)  \nTaking time derivative of Eq. 3, the accretion rate / 𝑀 can be expressed as a function of mass 𝑀 :  \n/ 𝑀 ( 𝑀 ) = ( 𝛽𝑀 / 𝑡 0 )( 𝑀 / 𝑀 0 ) -1 / 𝛽 ∝ 𝑀 1 -1 / 𝛽 ≃ 𝑀 -0 . 37 . (5)  \nThis scaling is also predicted by (spherically symmetric) similarity solutions for polytropic gas with 𝑃 ∝ 𝑛 𝛾 eff given 𝛾 eff = 1 . 09 (e.g., Yahil 1983; Suto & Silk 1988; Omukai & Nishi 1998; Tan & McKee 2004), as well as cosmological zoom-in simulations of the accretiondriven growth of a Pop III protostar (Bromm & Loeb 2004).  \nNext, we set the characteristic timescale 𝑡 0 using the relation between time and density normalisation parameters in self-similar transformations of the hydrodynamic equations (see eq. 1 in Liu et al. 2021a)  \n𝑡 0 = 1 / √︁ 4 𝜋𝐺𝜌 0 ≃ 3 × 10 4 yr , (6)  \nwhere 𝜌 0 = 𝜇𝑚 H 𝑛 0 , given the proton mass 𝑚 H , the mean molecular weight of primordial gas 𝜇 ≃ 1 . 22, and the typical density 𝑛 0 = 10 6 cm -3 around which primordial gas in run-away collapse under H 2 cooling settles to the effective EoS 𝑃 ∝ 𝑛 1 . 09 . Here, 𝑛 0 = 10 6 cm -3 is also the typical density at the edge of a Pop III star-forming disc (see, e.g., fig. 8 in Toyouchi et al. 2023). Finally, to determine 𝑀 0 , we equate the accretion rate of the protostar at 𝑡 0 , / 𝑀 | 𝑡 = 𝑡 0 = 𝛽𝑀 0 / 𝑡 0 , to the cloud-scale initial gas inflow rate / 𝑀 in onto the disc, multiplied by a factor 𝜂 ∼ 0 . 4 -0 . 6 (Sakurai et al. 2016; Toyouchi et al. 2023) that captures the rotational support against collapse, which gives  \n𝑀 0 = 𝜂 / 𝑀 in 𝑡 0 / 𝛽 , (7)  \n𝑀 ( 𝑡 ) = ( 𝜂 / 𝛽 ) / 𝑀 in 𝑡 1 -𝛽 0 𝑡 𝛽 ∝ ( 𝜂 / 𝛽 ) / 𝑀 in 𝑛 ( 𝛽 -1 )/ 2 0 𝑡 𝛽 . (8)  \nThe effects of other aspects, such as ejections of protostars by dynamical interactions and magnetic/radiative outflows, can also be absorbed into the parameter 𝜂 . Thanks to the simple power-law scaling, once the slope 𝛽 = 4 -3 𝛾 eff is known, the mass evolution is completely governed by one parameter, 𝐴 ≡ ( 𝜂 / 𝛽 ) / 𝑀 in 𝑡 1 -𝛽 0 . This also indicates that our results are not very sensitive to the disc edge density 𝑛 0 given the weak dependence of 𝐴 ∝ 𝑡 1 -𝛽 0 ∝ 𝑛 ( 𝛽 -1 )/ 2 0 = 𝑛 -0 . 135 0 on 𝑛 0 . The cloud-scale initial gas inflow rate / 𝑀 in and 𝜂 play more important roles, which capture the dynamics of the collapsing cloud evaluated at the density (time) scale 𝑛 0 ( 𝑡 0 ).  \nEq. 4 provides a good estimate of the radius of the star-forming disc in the early stage. A similar scaling relation 𝑅 ∝ 𝐾 -1 / 𝛽 ( 𝑀 / 𝜂 ) 𝛿 / 𝛽 is obtained in Tan & McKee (2004, see their eq. 17) based on angular momentum conservation and analytical spherical collapse solutions (McKee & Tan 2002, 2003), where 𝐾 ≡ 𝑃 / 𝑛 𝛾 eff ∝ 𝑇𝑛 1 -𝛾 eff is the entropy parameter. The normalisation of their scaling law is consistent with that in Eq. 4 learned from simulations if the disc mass is similar to the mass of protostars ( 𝜂 ∼ 0 . 5) and the typical value of  \nHD or atomic-cooling (e.g., Hirano et al. 2014; Regan 2023; Gurian et al. 2024b). For simplicity, we adopt the standard value 𝛾 eff = 1 . 09 for H 2 cooling by default and explore the dependence of our results on 𝛾 eff in Appendix A2.  \n𝐾 is adopted for infalling gas with 𝑇 ∼ 300 K and 𝑛 ∼ 10 6 cm -3 . However, this simple power-law scaling can produce unphysically large radii ( ≳ 10 pc) when 𝑀 becomes very large ( ≳ 10 4 M ⊙ ) in the late stage if stellar feedback is inefficient (see Sec. 2.2). Therefore, we further impose an upper limit 𝑅 max on the disc radius that can be related to the (initial) cloud size (see Sec. 2.2 below) and introduce the dependence on 𝜂 and 𝐾 :  \n𝑅 ( 𝑀 ) = GLYPH<26> 𝑅 max , ( 𝐾 / 𝐾 ' ) -1 / 𝛽 h 10 5 𝛽 𝑀 /( 800 𝜂 M ⊙ ) i 𝛿 / 𝛽 AU GLYPH<27> . (9)  \nFor simplicity, 𝐾 is fixed throughout this paper as 𝐾 = 𝐾 ' ≡ 𝑘 B ( 300 K )( 10 6 cm -3 ) 1 -𝛾 eff , where 𝑘 B is the Boltzmann constant, although in reality 𝐾 depends on the detailed thermo-chemical evolution of the collapsing gas. In fact, as discussed below, if 𝑀 becomes very large ( ≳ 10 4 M ⊙ ) before star formation is terminated by feedback, the cloud must undergo violent collapse with high inflow rates / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 . This usually happens in hotter, atomic-cooling clouds, where 𝐾 can be higher than that adopted here, leading to smaller discs. In a companion paper (Gurian et al. 2024a), we develop a semi-analytical model for the thermo-chemistry and collapse dynamics, which can be used to calculate / 𝑀 in and 𝐾 under different conditions of Pop III star formation (see also Li et al. 2021; Sharda & Krumholz 2022; Smith et al. 2024). We plan to explore the correlation between 𝐾 and / 𝑀 in in further work using this model 2 .", '2.2 Maximum/final mass of a single star per cloud': 'We first consider the simple case that only one protostar exists in the cloud, so Eqs. 8 and 9 delineate the (smoothed) median mass and radius evolution of the protostar and circumstellar disc, respectively, governed only by gravity and hydrodynamics . For simplicity, we assume that Eqs. 8 and 9 hold for 𝑡 ≳ 𝑡 0 until star formation is terminated quasi-instantaneously by factors other than gravity and hydrodynamics, such as stellar feedback, stellar collapse/explosion, and gas supply, although in reality they may cause gradual reductions of accretion with long-term deviations from the scaling in Eqs. 5 and 8 (McKee & Tan 2008). Each terminating effect corresponds to a characteristic maximum/final mass. Below, we consider five characteristic masses, which are combined to derive a universal formula for the final mass ˆ 𝑀 ★ of the central Pop III star, as illustrated in Fig. 1.', '(i) Final mass regulated by photo-ionization feedback ( 𝑀 ★, f ):': "In our previous work (Gurian et al. 2024b), we derive the final mass regulated by photo-ionization feedback assuming that a D-type shock is launched into the cloud once the protostar is thermally relaxed to produce ionizing photons efficiently (Alvarez et al. 2006), and star formation is shut down when the homogenized downstream medium is fully ionized (which marks the evaporation of the cloud). However, the final masses estimated in this way are lower than those found in 2D radiative hydrodynamic simulations (Hirano et al. 2014, 2015) by a factor of few. This is mainly caused by the disregard of the circumstellar disc, which can still feed the protostar even if bipolar ionized bubbles have expanded significantly. In fact, it has long been argued that photo-evaporation of the primordial protostellar disc may ultimately set the upper mass scale for Pop III stars (e.g., McKee & Tan 2008; Hosokawa et al. 2011; Stacy et al. 2012; Fukushima et al.  \n2020). As shown in Toyouchi et al. (2023), the feedback-regulated final mass can be estimated from the mass loss rate of the circumstellar disc by photo-evaporation (Tanaka et al. 2013)  \n/ 𝑀 pe ≃ 0 . 015 M ⊙ yr -1 GLYPH<18> / 𝑄 ion 10 52 s -1 GLYPH<19> 1 / 2 GLYPH<18> 𝑅 pe 10 4 AU GLYPH<19> 1 / 2 , (10)  \nwhere / 𝑄 ion is the production rate of ionizing photons (as a function of 𝑀 ), and 𝑅 pe is the physical scale within which the disc suffers from photo-evaporation. Assuming that the accretion rate of the protostar is comparable to the gas inflow rate onto the circumstellar disc, the final mass 𝑀 ★, f can be obtained by solving / 𝑀 pe = / 𝑀 , above which the disc is rapidly evaporated by ionizing radiation and can no longer feed the protostar. In Toyouchi et al. (2023), both / 𝑀 and 𝑅 pe are fixed as / 𝑀 = 𝜂 / 𝑀 in and 𝑅 pe = 10 4 AU. Here, we further consider their time evolution, i.e., mass dependence, using the scaling relations in Sec. 2.1. We set 𝑅 pe to the circumstellar disc radius given by Eq. 9 (i.e., 𝑅 pe ∝ 𝑀 1 . 25 ) and use Eq. 5 to derive / 𝑀 ( 𝑀 ) . Finally, to calculate / 𝑄 ion , we use the fitting formula from Schaerer (2002) for non-accreting, non-rotating 3 Pop III stars with 𝑀 ∼ 5 -500 M ⊙ combined with the linear relation for massive stars radiating at the Eddington limit (Bromm et al. 2001; Johnson et al. 2012):  \n/ 𝑄 ion ( 𝑀 ) ≃ ( 10 𝑎 0 + 𝑎 1 𝑥 + 𝑎 2 𝑥 2 s -1 , 𝑀 ≤ 224 M ⊙ , 1 . 55 × 10 48 + 𝑥 s -1 , 𝑀 > 224 M ⊙ , (11)  \nwhere 𝑥 ≡ log ( 𝑀 / M ⊙ ) , 𝑎 0 = 43 . 61, 𝑎 1 = 4 . 9, and 𝑎 2 = -0 . 83. We assume that the effects of accretion on stellar evolution, such as initial swelling and expansion by radiation pressure enhanced by the accretion shock (Hosokawa & Omukai 2009) become unimportant when the disc is about to diminish under photo-ionization feedback. This assumption is expected to be valid in most cases with / 𝑀 in ≳ 10 -4 M ⊙ yr -1 , where the protostar can reach 𝑀 ≳ 10 M ⊙ within an accretion timescale 𝑡 acc ( 𝑀 ) = 𝑡 0 ( 𝑀 / 𝑀 0 ) 1 / 𝛽 ≳ 0 . 3 Myr that is longer than the Kelvin-Helmholtz (KH) timescale 4 𝑡 KH . Note that disc photo-evaporation only happens after the ionized bubble breaks out into the cloud, while it is shown in the simulations by Jaura et al. (2022) that the bubble can be trapped at the center in the early stage for an extended period 𝑡 trap ≳ 2 × 10 4 yr. Here, we assume 𝑡 acc > 𝑡 trap such that the initial confinement of ionizing photons does not affect / 𝑀 ★, f . We find that this assumption is valid as long as 𝑡 trap ≲ 𝑡 0 = 3 × 10 4 yr.  \n(ii) Maximum mass gained during the feedback-free bloating phase ( 𝑀 ★, b ): In the calculation of 𝑀 ★, f , we assume that the protostar is always able to thermally relax and contract, which is only valid for low accretion rates. In fact, the protostar will expand significantly to enter a bloating phase (Omukai & Palla 2001, 2003), when the accretion rate is above a critical value / 𝑀 ★, crit ∼ 0 . 01 -0 . 04 M ⊙ yr -1 (Hosokawa et al. 2013; Haemmerlé et al. 2018; Herrington et al.  \n3 Pop III stars are likely born as fast rotators due to rapid accretion of gas with high angular momentum and inefficient magnetic braking during protostar growth (e.g., Stacy et al. 2011, 2013; Hirano & Bromm 2018; Kimura et al. 2023). In the extreme case, Pop III stars can undergo chemically-homogeneous evolution under efficient mixing triggered by fast rotation. As a result, the star becomes more compact and hotter, such that the production efficiency of ionizing photons 𝐵 ≡ / 𝑄 ion / 𝑀 is boosted by a factor of ∼ 2 -8 (Sibony et al. 2022). The final mass regulated by photo-ionization feedback will then be reduced slightly (by up to a factor of 2), as will be shown below that 𝑀 ★, f ∝ 𝐵 -0 . 33 when 𝐵 is independent of 𝑀 for massive stars.  \n2023; Nandal et al. 2023). In this case, the surface temperature becomes as low as 𝑇 eff ∼ 5000 -8000 K (see fig. 1 in Toyouchi et al. 2023), and the emission of ionizing photons is negligible (with a rate lower than that from a contracted star of the same mass by a factor of ∼ 10 10 ), which leads to effectively feedbackfree evolution. According to Eq. 5, the accretion rate of the protostar slowly decreases with time and will drop below / 𝑀 ★, crit at 𝑡 > 𝑡 b = 𝑡 0 GLYPH<0> 𝜂 / 𝑀 in / / 𝑀 ★, crit GLYPH<1> 1 /( 1 -𝛽 ) so the maximum mass accreted during the bloating phase is 𝑀 ★, b = 𝑀 0 ( 𝑡 b / 𝑡 0 ) 𝛽 ∝ / 𝑀 𝛽 /( 1 -𝛽 ) in .  \n(iii) Maximum mass regulated by the stellar lifetime ( 𝑀 ★, l ): Beyond stellar feedback, star formation also terminates when the central protostar runs out of fuel for nuclear fusion in its core and then collapses/explodes. To the first order, this means that the time taken for the protostar to reach a mass 𝑀 , i.e., the accretion timescale 𝑡 acc ( 𝑀 ) , cannot exceed the lifetime 𝑡 ★ ( 𝑀 ) of a star with the same mass. The maximum mass regulated by the stellar lifetime 𝑀 ★, l satisfies 𝑡 acc ( 𝑀 ★, l ) = 𝑡 ★ ( 𝑀 ★, l ) . It is non-trivial to calculate the lifetimes of accreting massive stars (e.g., Begelman 2010; Johnson et al. 2012; Haemmerlé 2021; Nandal et al. 2024a; Saio et al. 2024; Shibata et al. 2024). Here, we simply estimate 𝑡 ★ with (extrapolation of) the fitting formula from Schaerer (2002) for non-accreting Pop III stars (in the mass range ∼ 5 -500 M ⊙ ) bound by a lower limit 𝑡 ★, min to be calibrated (see below):  \n𝑡 ★ ( 𝑀 ) = max GLYPH<16> 𝑡 ★, min , 10 𝑏 0 + 𝑏 1 𝑥 + 𝑏 2 𝑥 2 + 𝑏 3 𝑥 3 yr GLYPH<17> , (12)  \nwhere 𝑥 ≡ log ( 𝑀 / M ⊙ ) , 𝑏 0 = 9 . 785, 𝑏 1 = -3 . 759, 𝑏 2 = 1 . 413, and 𝑏 3 = -0 . 186.  \n(iv) Maximum mass regulated by general-relativity instability (GRI, 𝑀 ★, g ): When the accretion rate is very high ( / 𝑀 > / 𝑀 ★, crit , GR ≃ 0 . 1 M ⊙ yr -1 ), the star will collapse by GRI before running out of fuel, and the final stellar mass can be estimated with the fitting formula (Woods et al. 2017; Li et al. 2021)  \n𝑀 ★, GR = GLYPH<20> 0 . 83 log GLYPH<18> / 𝑀 M ⊙ yr -1 GLYPH<19> + 2 . 48 GLYPH<21> × 10 5 M ⊙ , (13)  \nassuming a constant accretion rate (see also Haemmerlé 2021, 2024; Nandal et al. 2024a; Saio et al. 2024; Shibata et al. 2024). In our case, the accretion rate / 𝑀 decreases slowly as / 𝑀 ( 𝑀 ) ∝ 𝑀 -0 . 37 ∝ 𝑡 -0 . 27 (Eq. 5), and the protostar with a current mass 𝑀 spends most of its time at accretion rates similar to the current rate / 𝑀 ( 𝑀 ) . We use the current rate to evaluate Eq. 13 and find the final mass regulated by GRI 𝑀 ★, g by solving the equation 𝑀 = 𝑀 ★, GR ( / 𝑀 ( 𝑀 )) for 𝑀 .  \n(v) Maximum mass regulated by the gas supply ( 𝑀 c ): Finally, the stellar mass cannot exceed the total mass of gas 𝑀 c available in the collapsing cloud.  \nAmong the five maximum masses, the first two are closely related to each other. If 𝑀 ★, b > 𝑀 ★, f , the circumstellar disc will evaporate immediately after the bloating phase ends, and the final mass will be close to 𝑀 ★, b , which we call the bloating mode. On the other hand, if 𝑀 ★, f > 𝑀 ★, b , the protostar can further grow after the bloating phase until it reaches 𝑀 ★, f to destroy the disc, which we call the feedback-regulated mode. Therefore, if we only consider gravity, hydrodynamics and photo-ionization feedback, the final mass can be well estimated by ˆ 𝑀 ★ = max GLYPH<0> 𝑀 ★, b , 𝑀 ★, f GLYPH<1> . Considering the other three masses that serve as upper limits, we have  \nˆ 𝑀 ★ = min GLYPH<2> max GLYPH<0> 𝑀 ★, b , 𝑀 ★, f GLYPH<1> , 𝑀 ★, l , 𝑀 ★, g , 𝑀 c GLYPH<3> . (14)  \nThere are nine parameters in our single-star model: / 𝑀 in , 𝑛 0 , 𝑀 c , 𝐾 , 𝛾 eff , 𝜂 , 𝑅 max , / 𝑀 ★, crit , and 𝑡 ★, min . The first three reflect the properties of star-forming clouds (and host haloes) on large scales ( ≳ 1 pc).  \nThelast two parameters are relevant for primordial stellar evolution at small (sub-AU) scales. These two aspects are connected to each other through our model for star-forming discs (Eqs. 3-9) at intermediate scales governed by 𝐾 , 𝛾 eff , 𝜂 , and 𝑅 max .  \nFor illustration, below we treat / 𝑀 in at 𝑛 0 = 10 6 cm -3 as the independent variable and use physically-motivated values for the other parameters. We adopt 𝛾 eff = 1 . 09 for typical primordial collapsing gas under H 2 cooling with 𝐾 = 𝐾 ' . The cloud mass 𝑀 c is associated with / 𝑀 in through the collapse timescale 𝑡 col , i.e., 𝑀 c = / 𝑀 in 𝑡 col . As an optimistic estimate, we set the collapse timescale to the free-fall timescale 5 𝑡 ff = √︁ 3 𝜋 /( 32 𝐺𝑛 crit 𝜇𝑚 H ) ≃ 1 . 5Myrforacritical density 𝑛 crit = 10 3 cm -3 , beyond which H 2 cooling saturates in minihaloes (Abel et al. 2002; Bromm et al. 2002; Gurian et al. 2024b) and H 2 formation becomes efficient in atomic-cooling haloes (Sugimura et al. 2024). This density also marks the onset of collapse by runaway cooling when a core with a enclosed mass larger than the BonnorEbert mass first emerges (Gurian et al. 2024a). Given the cloud mass and the critical density, we estimate the initial cloud size as 𝑅 c = [ 3 𝑀 c /( 4 𝜋𝜇𝑛 crit 𝑚 H )] 1 / 3 , which then sets the upper limit of the circumstellar disc radius 𝑅 max = 𝑅 c in Eq. 9. We further fix 𝑡 ★, min using the critical accretion rate / 𝑀 ★, crit , GR ≃ 0 . 1 M ⊙ yr -1 above which GRI collapse is expected to occur (Woods et al. 2017, see their fig. 4). To be specific, we require that / 𝑀 = / 𝑀 ★, crit , GR when 𝑀 ★, l = 𝑀 ★, g , which gives 𝑡 ★, min ≃ 1 . 2 Myr. We adopt an intermediate value 𝜂 = 0 . 5 based on the finding 𝜂 ∼ 0 . 4 -0 . 6 from simulations (Sakurai et al. 2016; Toyouchi et al. 2023). This means that half of the infalling gas goes into the protostar system, while the rest constitutes the star-forming disc (or leaves the star-forming cloud via ejected protostars and outflows). Following Toyouchi et al. (2023), we use / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 by default to obtain conservative estimates of ˆ 𝑀 ★ . In reality, these input parameters can have non-trivial correlations with / 𝑀 in and each other and show a strong dependence on environmental factors (e.g., LW radiation and streaming motion between baryons and dark matter). For simplicity, we ignore such correlations in this section and explore their possible impact on the results in Appendix A.  \nFig. 1 shows the characteristic masses 𝑀 ★, b , 𝑀 ★, f , 𝑀 ★, l , 𝑀 ★, g , 𝑀 c , and the final mass ˆ 𝑀 ★ as functions of / 𝑀 in . The transition from the feedback-regulated mode ( ˆ 𝑀 ★ = 𝑀 ★, f ) to the bloating mode ( ˆ 𝑀 ★ = 𝑀 ★, b ) occurs almost exactly at / 𝑀 in = / 𝑀 ★, crit / 𝜂 = 0 . 08 M ⊙ yr -1 , corresponding to ˆ 𝑀 ★ ∼ 2000 M ⊙ . Approaching this transition point, whenthe production rate of ionizing photons is proportional to stellar mass as / 𝑄 ion ( 𝑀 ) = 𝐵𝑀 for massive stars with 𝑀 ≳ 200 M ⊙ , 𝑀 ★, f can be derived analytically for 𝑅 < 𝑅 max :  \n/ 𝑀 ( 𝑀 ★, f ) = / 𝑀 pe ( 𝑀 ★, f ) , (15)  \n( 𝜂 / 𝑀 in ) 1 𝛽 ( 𝛽𝑀 ★, f / 𝑡 0 ) 𝛽 -1 𝛽 = 𝐶 ( 𝐵𝑀 ★, f ) 1 2 𝐾 -1 2 𝛽 ( 𝑀 ★, f / 𝜂 ) 𝛿 2 𝛽 , (16)  \n( 𝜂 / 𝑀 in ) 1 𝛽 ( 𝛽 / 𝑡 0 ) 𝛽 -1 𝛽 𝑀 𝛽 -1 𝛽 ★, f = 𝐶𝐵 1 2 𝐾 -1 2 𝛽 𝜂 -𝛿 2 𝛽 𝑀 𝛽 + 𝛿 2 𝛽 ★, f , (17)  \n𝑀 𝛽 -2 -𝛿 2 𝛽 ★, f = ( 𝛽 / 𝑡 0 ) 1 -𝛽 𝛽 𝐶𝐵 1 2 𝐾 -1 2 𝛽 𝜂 -𝛿 + 2 2 𝛽 / 𝑀 -1 𝛽 in , (18)  \n𝑀 ★, f = ( 𝛽 / 𝑡 0 ) 2 ( 𝛽 -1 ) 𝜖 𝐶 -2 𝛽𝜖 𝐵 -𝛽𝜖 𝐾 𝜖 𝜂 ( 𝛿 + 2 ) 𝜖 / 𝑀 2 𝜖 in . (19)  \nFigure 1. Final stellar mass ( ˆ 𝑀 ★ , thick solid), and the characteristic masses regulated by photo-ionization feedback ( 𝑀 ★, f , solid), bloating phase ( 𝑀 ★, b , dashed), stellar lifetime ( 𝑀 ★, l , dash-dotted), GRI ( 𝑀 ★, g , dotted), and gas supply ( 𝑀 c , thick dashed) as functions of the initial gas inflow rate / 𝑀 in . The fitting formula of final stellar mass (Eq. 20) based on 2D simulations from Hirano et al. (2014, see their fig. 14) is shown with the long-dashed line, and the dots show the underlying data. The results of the 3D AMR simulations in Sugimura et al. (2023) with rescaled / 𝑀 in (see the main text) are denoted by the squares, while those of the 3D spherical simulations in Toyouchi et al. (2023) are shown with the diamonds. The curve for 𝑀 ★, g is only meaningful when / 𝑀 ( 𝑀 ★, g ) ≳ 0 . 1 M ⊙ yr -1 , corresponding to 𝜂 / 𝑀 in ≳ 0 . 2 M ⊙ yr -1 .  \n<!-- image -->  \nFigure 2. Same as Fig. 1 but assuming a fixed scale 𝑅 pe = 10 4 AU for disc photo-evaporation following Toyouchi et al. (2023). As expected, the simulation results from Toyouchi et al. (2023, denoted by the diamonds) are better reproduced here compared with the case in Fig. 1 with disc size evolution 𝑅 pe ∝ 𝑀 1 . 25 .  \n<!-- image -->  \nHere we have simplified the final expression by defining 𝜖 ≡ 1 /( 2 + 𝛿 -𝛽 ) ≃ 0 . 46 and used Eqs. 5, 7, 9, and 10 to obtain the second line given the constant 𝐶 ≡ 0 . 015 M ⊙ yr -1 h 10 -56 𝐾 ' 1 / 𝛽 10 5 𝛿 ( 800 M ⊙ ) -𝛿 / 𝛽 s i 1 / 2 . A quasilinear relation between ˆ 𝑀 ★ = 𝑀 ★, f and / 𝑀 in is found as 𝑀 ★, f ∝ 𝐵 -𝛽𝜖 𝐾 𝜖 𝜂 ( 𝛿 + 2 ) 𝜖 / 𝑀 2 𝜖 in ≃ 𝐵 -0 . 33 𝐾 0 . 46 𝜂 1 . 33 / 𝑀 0 . 92 in . On the other hand, in the bloating mode, the final stellar mass increases rapidly with / 𝑀 in as ˆ 𝑀 ★ = 𝑀 ★, b ∝ / 𝑀 𝛽 /( 1 -𝛽 ) in ≃ / 𝑀 2 . 7 in before it reaches the limits placed by the stellar lifetime and GRI. The gas supply ( 𝑀 c ) never becomes a limiting factor, and we always have 𝑡 acc > 𝑡 0 .  \nIn the feedback-regulated regime, our results ( 𝑀 ★, f ) agree perfectly with the fitting formula based on 2D simulations (for / 𝑀 in ∼ 10 -4 -0 . 1 M ⊙ yr -1 ) from Hirano et al. (2014, see their fig. 14)  \n𝑀 ★ = 100 M ⊙ GLYPH<18> / 𝑀 in 2 . 8 × 10 -3 M ⊙ yr -1 GLYPH<19> 0 . 8 , (20)  \nand the results of the 3D simulations with adaptive mesh refinement (AMR) over Cartesian coordinates in Sugimura et al. (2023), which produce 2 to 4 stars per cloud in the end rather than a single star. Here, we have applied a correction factor ≃ 3 . 7 to the initial gas inflow rates in Sugimura et al. (2023, see their table 1) to be consistent with the definition of / 𝑀 in in Hirano et al. (2014) 6 . The final mass is underestimated by our model at / 𝑀 in ≲ 10 -4 M ⊙ compared with the results of two simulations from Hirano et al. (2014) in this regime where the initial KH contraction phase likely becomes important. This so-called low-mass mode of Pop III star formation (Stacy & Bromm 2014) is expected to be rare.  \nWealso compare our predictions with the results of the 3D singlestar simulations with spherical coordinates in Toyouchi et al. (2023, see their table 3), finding a good agreement for / 𝑀 in ≲ 0 . 02 M ⊙ yr -1 . This shows that our model accurately captures the median evolution of Pop III protostars and circumstellar discs in the feedback-regulated regime at least for / 𝑀 in ≲ 0 . 02 M ⊙ yr -1 . However, two simulations in Toyouchi et al. (2023) with / 𝑀 in ∼ 0 . 06 -0 . 08 M ⊙ yr -1 obtain higher masses ( ∼ 10 4 M ⊙ ) than our results ( ˆ 𝑀 ★ ∼ 2000 M ⊙ ) and those from the 2D simulations in Hirano et al. (2014) under similar conditions. This implies that the final mass can vary greatly from system to system around the critical inflow rate ∼ / 𝑀 ★, crit / 𝜂 required to trigger the bloating phase, likely due to the deviation from the scaling relations (Eqs. 3-9) in the (stochastic) evolution of individual systems, which we will discuss in Sec. 3.2 and Appendix A2. In the lifetime/GRI-regulated regime, our model reproduces the result of Toyouchi et al. (2023) for / 𝑀 in ∼ 0 . 6 M ⊙ yr -1 within a factor of 2, showing a slight overestimation 7 .  \nFor comparison, Fig. 2 shows the results assuming a fixed scale 𝑅 pe = 10 4 AU for disc photo-evaporation following Toyouchi et al. (2023). As expected, the predicted ˆ 𝑀 ★ -/ 𝑀 in relation is generally consistent with the simulation data from Toyouchi et al. (2023). The final mass regulated by photo-ionization feedback 𝑀 ★, f is overestimated (by up to a factor of ∼ 10) at / 𝑀 in ∼ 0 . 02 -0 . 1 M ⊙ yr -1 compared with the fit from Hirano et al. (2014) 8 . In this way, the feedback-  \n6 Underlying the fitting formula Eq. 20 from Hirano et al. (2014), the initial gas inflow rate / 𝑀 in is measured at the moment the central hydrogen density reaches 𝑛 H , cen = 10 12 cm -3 and at the radius where the ratio between enclosed mass and Bonnor-Ebert mass reaches the maximum, corresponding to a density 𝑛 ∼ 10 6 cm -3 (see their figs. 11 and 13) consistent with our definition of / 𝑀 in for 𝑛 0 = 10 6 cm -3 . When / 𝑀 in is evaluated at an earlier stage with a smaller central density 𝑛 H , cen = 10 7 cm -3 in the 2D simulations of Hirano et al. (2015) as well as in the 3D simulations of Sugimura et al. (2023), a best-fit 𝑀 ★ = 250 M ⊙[ / 𝑀 in /( 2 . 8 × 10 -3 M ⊙ yr -1 ) ] 0 . 7 is obtained with a similar slope but a higher normalisation, because / 𝑀 in is lower compared with the 𝑛 H , cen = 10 12 cm -3 case. The two fits are consistent with each other given the relation / 𝑀 in ( 𝑛 H , cen = 10 12 cm -3 ) = 3 . 7 / 𝑀 in ( 𝑛 H , cen = 10 7 cm -3 ) , which is used to correct the inflow rates reported in Sugimura et al. (2023). 7 The final mass for / 𝑀 in ∼ 0 . 6 M ⊙ yr -1 in Toyouchi et al. (2023, see their LWH-10 model) is obtained by assuming constant mass growth for a duration of 0 . 5 Myr using the average accretion rate in the first 0.2 Myr covered by the simulation. Therefore, it should be regarded as a lower limit if the gas inflow can be sustained beyond 0.5 Myr.  \n- 8 Similarly, 𝑀 ★, f is overestimated by a factor of ∼ 10 in Li et al. (2021, see their eq. 27), where smaller values of 𝑅 pe are adopted based on the gravitational influence radii of protostars.  \nregulated mode is also important for / 𝑀 in ∼ 0 . 08 -0 . 2 M ⊙ yr -1 , and the bloating mode becomes rare. However, for / 𝑀 in ≲ 0 . 002M ⊙ yr -1 , 𝑀 ★, f is slightly lower for 𝑅 pe = 10 4 AU compared with the fit from Hirano et al. (2014) and our results with 𝑅 pe ∝ 𝑀 1 . 25 . This indicates that the evolution of 𝑅 pe is crucial for the accurate modelling of photo-ionization feedback and calculation of ˆ 𝑀 ★ .", '2.3 Effects of fragmentation/multiplicity': 'In the previous subsection, we derive the final stellar mass under the assumption that only one protostar forms per cloud. However, in the last two decades, small-scale simulations have shown that fragmentation is common in primordial star-forming clouds/discs (Bromm 2013; Haemmerlé et al. 2020; Klessen & Glover 2023), producing multiple protostars that undergo complex evolution and interactions (accretion, migration, gravitational scatters, and mergers). There is also evidence for multiplicity in Pop III star formation from stellar archaeology data (Hartwig et al. 2023). However, the final fates of the protostars are still in debate due to computational limitations (e.g., resolution, time coverage, volume, and treatment of feedback). For instance, in the dense discs formed under strong gas inflows ( / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 ), gaseous dynamical friction can drive rapid inward migrations of fragments/protostars to merge with one central object that dominates the total stellar mass, converging to the singlestar scenario, or feed a small number of SMSs and many low-mass stars (Inayoshi & Haiman 2014; Boekholt et al. 2018; Suazo et al. 2019; Chon & Omukai 2020; Latif et al. 2021; Woods et al. 2021, 2024; Prole et al. 2024; Reinoso et al. 2023; Schleicher et al. 2023; Kiyuna et al. 2024). In general, the interplay between fragmentation and feedback has profound implications for the formation efficiency and IMF of Pop III stars, similar to the case of Pop I/II star formation (e.g., Menon et al. 2024).  \nIn light of this, we generalize the single-star model in Sec. 2.2 to follow the growth and feedback of multiple protostars per cloud and explore the effects of multiplicity on the final total stellar mass. We assume that the cloud has modest turbulence at the onset of collapse due to inefficient cooling (Chon et al. 2021) and only forms one star-forming disc that fragments by disc instability. We only follow the median evolution of the mass distribution of surviving protostars and ignore transient protostars that are promptly lost in mergers or ejections 9 . As a heuristic attempt, we consider a simple phenomenological multiplicity model under the following two assumptions.  \n- (i) The mass distribution of protostars follows a power law 𝑑𝑁 / 𝑑𝑚 ∝ 𝑚 -𝛼 .  \n(ii) The fraction of the total accretion rate / 𝑚 / / 𝑀 that feeds a protostar is identical to the fraction of total stellar mass 𝑚 / 𝑀 occupied by the protostar. Therefore, if the number of protostars is constant, the mass distribution of protostars undergoes self-similar evolution: The shape (i.e., power-law index) remains unchanged and the boundaries increase proportionally with the total stellar mass 𝑀 .  \nHere, the second assumption is motivated by the quasi-linear relation between / 𝑚 / / 𝑀 and 𝑚 / 𝑀 found in the (magneto-) hydrodynamic  \n- 9 Although ejections are commonly seen in simulations (e.g., Prole et al. 2022a; Raghuvanshi & Dutta 2023; Reinoso et al. 2023), their effects on the final total stellar mass are still unclear. Since protostars can hardly grow once ejected from the star-forming cloud, the mass carried by ejected stars is expected to be negligible compared with the final mass of the remaining star cluster. However, the radiative feedback from ejected stars may also regulate star formation. We defer the inclusion of ejected stars to future work.  \nN  \n=3  \nN  \n=10  \nFigure 3. Relation between 𝑓 1 and 𝑓 2 for 𝛼 = 2 . 5 (solid), 2 (dashed), 1.5 (dash-dotted), 1 (dotted), 0.5 (long-dashed), and 0 (dash-dot-dotted) given 𝑁 = 3 (top-left), 10 (top-right), 30 (bottom-left), and 100 (bottom-right).  \n<!-- image -->  \nsimulations from Sharda et al. (2020) when the accretion rate is measured over a timescale ≳ 0 . 01 𝑡 KH in the early stage with 𝑀 ≲ 50 M ⊙ . Interestingly, this quasi-linear relation between accretion rate and mass appears to be an intermediate case compared with the extreme regimes of tidal lobe accretion ( / 𝑚 ∝ 𝑚 2 / 3 , Bonnell et al. 2001b) and Bondi-Hoyle accretion ( / 𝑚 ∝ 𝑚 2 , Bondi & Hoyle 1944). According to the theory of competitive accretion in star clusters (Bonnell et al. 2001b), tidal lobe accretion is expected to happen in the early stage when the cluster mass/potential is dominated by gas (Bonnell et al. 2001a), while Bondi-Hoyle accretion takes over once stars dominate the potential and are virialized. For Pop III star formation, the masses of stars and gas are comparable in the star-forming disc ( 𝜂 ∼ 0 . 5), so it is reasonable that the accretion rate-mass relation has a power-law slope ∼ 1 between 2 / 3 and 2.  \nWith self similarity, the protostar mass distribution is determined by three parameters 𝛼 , 𝑓 1 , and 𝑓 2 , where the last two parameters define the lower and upper bounds as 𝑚 min = 𝑓 1 𝑀 and 𝑚 max = 𝑓 2 𝑀 . These parameters are related to the number of protostars 𝑁 through  \n1 𝑁 = 1 𝑀 ∫ 𝑓 2 𝑀 𝑓 1 𝑀 𝑚 1 -𝛼 𝑑𝑚 ∫ 𝑓 2 𝑀 𝑓 1 𝑀 𝑚 -𝛼 𝑑𝑚 = ∫ 𝑓 2 𝑓 1 𝑓 1 -𝛼 𝑑𝑓 ∫ 𝑓 2 𝑓 1 𝑓 -𝛼 𝑑𝑓 , (21)  \nwhere 𝑓 ≡ 𝑚 / 𝑀 given 𝑚 the mass of individual protostars. We treat 𝛼 and 𝑓 2 as free parameters and associate 𝑓 1 to 𝑁 , such that 𝑁 becomes the third parameter of our model. Given 𝛼 , 𝑓 2 , and 𝑁 , we first solve Eq. 21 to obtain 𝑓 1 ( 𝑁, 𝑓 2 , 𝛼 ) , as shown in Fig. 3. Then we distribute 𝑁 bins 𝑓 𝑖 in the range 𝑓 ∈ ( 𝑓 1 , 𝑓 2 ) according to a power law with index 𝛼 and apply small corrections to ensure " 𝑁 𝑖 𝑓 𝑖 = 1. The masses and accretion rates of individual protostars  \nare simply 𝑚 𝑖 = 𝑓 𝑖 𝑀 and / 𝑚 𝑖 = 𝑓 𝑖 / 𝑀 given the total stellar mass 𝑀 and accretion rate / 𝑀 ( 𝑀 ) of the protostar cluster (Eqs. 5 and 7). Here, we assume that the evolution of / 𝑀 with 𝑀 (governed by cloud-scale collapse) is completely determined by / 𝑀 in and not affected by the fragmentation process before feedback kicks in, since recent simulations have shown that the growth of total stellar mass is insensitive to the number of protostars at least in the early stage without feedback (see, e.g., Wollenberg et al. 2020; Liu et al. 2021a; Prole et al. 2022a). In this way, the effects of fragmentation on the final total stellar mass only manifest via its interplay with feedback , as discussed below.  \nOnce the relation between 𝑀 and the masses and accretion rates of individual protostars is known, we generalize the analysis in Sec. 2.2 to multiple protostars. We first absorb the effects of the bloating phase into the calculation of ionizing photon production rates, in which we reduce the value from Eq. 11 by a factor of 10 10 when the accretion rate is above / 𝑀 ★, crit (see fig. 4 in Toyouchi et al. 2023). We derive the total disc photo-evaporation rate under the assumption that all protostars reside in a central region much smaller than the starforming disc and contribute to the evaporation of the disc effectively as one point source. This central concentration may occur by inward migrations from dynamical friction (Riaz et al. 2023a). In this case, we simply replace / 𝑄 ion in Eq. 10 with the summation of the rates / 𝑞 ion ,𝑖 from individual protostars and still use the radius of the entire disc from Eq. 9 with 𝑅 max = 𝑅 c to estimate the disc evaporation scale  \n/ 𝑀 pe ≃ 0 . 015 M ⊙ yr -1 " 𝑁 𝑖 / 𝑞 ion ,𝑖 10 52 s -1 ! 1 / 2 GLYPH<20> 𝑅 ( 𝑀 ) 10 4 AU GLYPH<21> 1 / 2 . (22)  \nAn alternative scenario with distributed feedback is explored in Appendix A1.  \nNext, we consider the limits on the total stellar mass placed by the lifetime and GRI of the most massive protostar in the cluster (with the index 𝑖 = 𝑁 ), assuming that once this star ends its life, gas inflows will be shut down either by a supernova explosion or by black hole accretion feedback. To be specific, we now use 𝑚 𝑁 = 𝑓 𝑁 𝑀 rather than 𝑀 to evaluate the stellar lifetime 𝑡 ★ formula (Eq. 12) and compare 𝑚 𝑁 with the GRI mass from Eq. 13 (using / 𝑚 𝑁 = 𝑓 𝑁 / 𝑀 instead of / 𝑀 ). The final total stellar mass ˆ 𝑀 ★ is given by the maximum value of 𝑀 that satisfies  \n/ 𝑀 pe ( 𝑀, 𝑓 𝑖 ) ≤ / 𝑀 , 𝑡 acc ( 𝑀 ) = 𝑡 0 ( 𝑀 / 𝑀 0 ) 1 / 𝛽 ≤ 𝑡 ★ ( 𝑓 𝑁 𝑀 ) , 𝑓 𝑁 𝑀 ≤ 𝑀 ★, GR ( 𝑓 𝑁 / 𝑀 ) , (23)  \nand again capped by the cloud mass 𝑀 c .  \nTo demonstrate the effects of multiplicity on ˆ 𝑀 ★ , we explore the parameter space 10 of 𝛼 , 𝑓 2 , and 𝑁 with 3 × 4 = 12 examples while fixing / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 and 𝜂 = 0 . 5. Here, we assume that 𝑁 is constant (i.e., independent of 𝑀 ) for simplicity and consider the concentrated feedback scenario (Eq. 22). In practice, we impose an upper limit 𝑁 max on 𝑁 corresponding to 𝑓 1 ( 𝑁 max , 𝑓 2 , 𝛼 ) = 0 . 1 M ⊙ / 𝑀 when evaluating Eq. 23 to avoid the production of objects below the minimum mass ∼ 0 . 1 M ⊙ of Pop III protostars (Greif et al. 2011, see their fig. 7). Therefore, 𝑁 is not strictly constant in the model, and the input value of 𝑁 as a parameter should be regarded as the optimistic target value. These choices ( / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 , 𝜂 = 0 . 5, constant 𝑁 , and concentrated feedback) together with the assumptions of 𝑅 pe ∝ 𝑀 1 . 25 and 𝛾 eff = 1 . 09 are adopted as the default setup unless otherwise noted. More detailed results, including those with different prescriptions/assumptions and choices of parameters for fragmentation, mass-size scaling relation, and photo-ionization feedback are shown in Appendix A.  \nFig. 4 shows the relation between ˆ 𝑀 ★ and / 𝑀 in for 𝑁 = 3, 10, 30, and 100 given 𝑓 2 = 0 . 8 and 𝛼 = 2. Note that 𝑓 1 decreases with 𝑁 (Fig. 3), low-mass stars make up a larger fraction of the total stellar mass given a larger 𝑁 , which are less efficient at producing ionizing photons compared with massive stars if we do not consider the bloating phase. Therefore, as 𝑁 increases, photo-ionization feedback is weaker in the low-accretion rate regime 𝜂 / 𝑀 in ≲ / 𝑀 ★, crit without bloating, leading to higher ˆ 𝑀 ★ . The effect is stronger for smaller / 𝑀 in . On the contrary, for 𝜂 / 𝑀 in ≳ / 𝑀 ★, crit , the transition from the nobloating feedback-regulated regime ( ˆ 𝑀 ★ ≲ 2000 M ⊙ , where all stars have contracted before star formation is terminated by feedback) to the lifetime/GRI-regulated regime ( ˆ 𝑀 ★ ≳ 10 5 M ⊙ ) is smoother and delayed as 𝑁 increases, leading to smaller ˆ 𝑀 ★ in the intermediate stage (corresponding to the bloating mode in the single-star model). This is because the full bloating phase, where all protostars have accretion rates above / 𝑀 ★, crit , is more difficult to trigger when we distribute the total accretion rate among multiple protostars. Even if the most massive stars have entered the bloating phase, low-mass stars can still contract and produce ionizing photons efficiently under low accretion rates. The fraction of stars that experience bloating when star formation terminates increases with / 𝑀 in until it reaches unity. However, ˆ 𝑀 ★ increases again with 𝑁 in the GRI-regulated regime, because the GRI mass limit (Eq. 13) is applied to the most massive star, whose mass fraction in the cluster is lower with higher 𝑁 .  \nFigure 4. Final total stellar mass as a function of initial gas inflow rate for 𝑁 = 3 (thin solid), 10 (dashed), 30 (dash-dotted), and 100 (dotted) with 𝑓 2 = 0 . 8 and 𝛼 = 2. The results for the single-star model ( 𝑁 = 1) is shown with the thick solid curve for comparison. The fitting formula of final stellar mass (Eq. 20) based on 2D (single-star) simulations from Hirano et al. (2014, see their fig. 14) is shown with the long-dashed line, and the dots show the underlying data. The results of the 3D AMR simulations in Sugimura et al. (2023) are denoted by the squares, while those of the 3D spherical simulations in Toyouchi et al. (2023) are shown with the diamonds.  \n<!-- image -->  \n1  \nFigure 5. Same as Fig. 4 but for 𝑓 2 = 0 . 8 (thin solid), 0.6 (dashed), 0.4 (dash-dotted), and 0.2 (dotted) with 𝑁 = 10 and 𝛼 = 1.  \n<!-- image -->  \nFig. 5 shows the results for 𝑓 2 = 0 . 8, 0.6, 0.4, and 0.2 given 𝑁 = 10 and 𝛼 = 1. Similarly to trends with 𝑁 , ˆ 𝑀 ★ decreases with 𝑓 2 for 𝜂 / 𝑀 in ≲ / 𝑀 ★, crit and in the GRI-regulated regime, but increases with 𝑓 2 in the intermediate regime. This is because the importance of massive stars simply increases with 𝑓 2 , resulting in stronger feedback if there is no bloating but meanwhile making the bloating phase and GRI easier to occur. Besides, since 𝑓 1 increases when 𝑓 2 is smaller (Fig. 3), the transition to the GRI-regulated regime is generally sharper.  \nFig. 6 shows the results for 𝛼 = 2 . 5, 1.5, 1, and 0.5 given 𝑁 = 3 and 𝑓 2 = 0 . 8. Here, the dependence on 𝛼 is generally weaker compared to the trends in 𝑁 and 𝑓 2 . The biggest effect is that the full bloating phase occurs at lower / 𝑀 in when the protostar mass distribution is more bottom-heavy with higher 𝛼 because 𝑓 1 decreases with 𝛼 when the input values of 𝑁 and 𝑓 2 are fixed (Fig. 3).  \nFigure 6. Same as Fig. 4 but for 𝛼 = 2 . 5 (thin solid), 1.5 (dashed), 1 (dashdotted), and 0.5 (dotted) with 𝑁 = 3 and 𝑓 2 = 0 . 8.  \n<!-- image -->', '3.1 Implications for the cluster mass distribution and IMF of Pop III stars in quasar progenitor haloes': "As an example, we apply our analytical model for Pop III star formation to the primordial star-forming clouds in the (main) progenitor haloes of high𝑧 luminous quasar host galaxies and derive the mass distributions of Pop III clusters ( 𝑀 ★ ) and stars ( 𝑚 ★ ). The (cloudscale initial) gas inflow rates / 𝑀 in of such clouds are calculated by Li et al. (2021) with merger trees targeting typical luminous quasar host haloes with masses 𝑀 h = 10 12 M ⊙ and a co-moving number density 11 of 𝑛 h ∼ 2 . 3 × 10 -6 cMpc -3 at 𝑧 = 6, corresponding to ∼ 4 𝜎 over-dense peaks that only occupy a small fraction ≲ 6 × 10 -5 of cosmic volume in the standard Λ CDMcosmology (Planck Collaboration et al. 2016). They use a one-zone model to follow the cloud collapse until 𝑛 = 10 6 cm -3 to predict the gas inflow rate, consistent with our definition of / 𝑀 in for 𝑛 0 = 10 6 cm -3 . In addition to standard thermo-chemistry of collapsing primordial gas, they also model the effects of Lyman-Werner (LW) radiation from nearby star-forming galaxies, dynamical heating by halo mergers, and streaming motion between baryons and dark matter, which can delay Pop III star formation to more massive (atomic-cooling) clouds/haloes with higher gas inflow rates, favoring the formation of SMSs and heavy BH seeds (see, e.g., Inayoshi et al. 2020; Haemmerlé et al. 2020; Volonteri et al. 2021; Wise 2023; Regan & Volonteri 2024, for reviews).  \nHere we consider the / 𝑀 in distributions predicted by Li et al. (2021) for four values of baryon-dark matter streaming velocity 𝑣 bc = 0, 1, 2, and 3 𝜎 rms (where 𝜎 rms is the cosmic root-mean-squared velocity), as shown in Fig. 7. The / 𝑀 in distribution for 𝑣 bc = 0 covers a broad range / 𝑀 in ∼ 0 . 003 -1 M ⊙ yr -1 with multiple peaks corresponding to haloes of different assembly histories and LW background intensities (Li et al. 2021, see their fig. 6), and the distribution is narrower and dominated by higher / 𝑀 in as 𝑣 bc increases, reaching up to / 𝑀 in ∼ 5 M ⊙ yr -1 for 𝑣 bc = 3. We further construct the / 𝑀 in distribution in the cosmic (volume) average case by numerically integrating the  \nFigure 7. Gas inflow rate distribution of primordial star-forming clouds in high𝑧 luminous quasar progenitor haloes predicted by Li et al. (2021, see their fig. 6) for 𝑣 bc = 0 (solid), 1 (dashed), 2 (dash-dotted) and 3 𝜎 rms (dotted), each with 10000 realisations. The shaded region shows the cosmic average distribution estimated by numerical integration of the distributions for the four 𝑣 bc values using the trapezium rule, given the probability density function 𝑝 ( 𝑣 bc ) of 𝑣 bc (Eq. 24), following the same normalisation, i.e., ∫ ( 𝑑𝑁 / 𝑑 / 𝑀 in ) 𝑑 / 𝑀 in = 10000. Here, the number counts for different 𝑣 bc models are not directly comparable, as they are not normalized to reflect 𝑝 ( 𝑣 bc ) but to have the same total number of haloes/clouds.  \n<!-- image -->  \nresults for the four 𝑣 bc values over the probability density function of 𝑣 bc (Tseliakhovich et al. 2011; Fialkov 2014; Schauer et al. 2019b)  \n𝑝 ( 𝑣 bc ) = GLYPH<18> 3 2 𝜋𝜎 rms GLYPH<19> 3 / 2 4 𝜋𝑣 2 bc exp -3 𝑣 2 bc 2 𝜎 2 rms ! , (24)  \nusing the trapezium rule. To be specific, the results for 𝑣 bc / 𝜎 rms = 0, 1, 2, and 3 are combined with the weights 𝑤 0 = 0 . 5 𝐹 ( 0 , 1 ) ≃ 0 . 304, 𝑤 1 = 0 . 5 [ 𝐹 ( 0 , 1 ) + 𝐹 ( 1 , 2 )] ≃ 0 . 496, 𝑤 2 = 0 . 5 [ 𝐹 ( 1 , 2 ) + 𝐹 ( 2 , 3 )] ≃ 0 . 196, and 𝑤 3 = 0 . 5 𝐹 ( 2 , 3 ) + 𝐹 ( 3 , ∞) ≃ 0 . 004, respectively, where the function 𝐹 is defined as 𝐹 ( 𝑥, 𝑦 ) ≡ ∫ 𝑦𝜎 rms 𝑥𝜎 rms 𝑝 ( 𝑣 bc ) 𝑑𝑣 bc . Since we only have 4 bins of 𝑣 bc , the numerical integration is only meant to capture the general shape of the cosmic average distribution, while the detailed features in the integrated distribution may be numerical artifacts, which can be smoothed out given a finer grid of 𝑣 bc .  \nFig. 8 shows the resulting mass distributions of Pop III stars in the single-star model ( 𝑁 = 1, 𝑚 ★ = 𝑀 ★ ) under the default setup for 𝑣 bc = 0, 1, 2, and 3 𝜎 rms as well as the cosmic average. For comparison, we also plot the mass distribution of Pop III stars formed in typical minihaloes in unbiased regions with 𝑣 bc = 0 from the 2D (single-star) simulations in Hirano et al. (2015). Some statistical properties of these mass distributions are summarized in Table 1. In the simple case with no streaming motion ( 𝑣 bc = 0), the mass distribution spans a wide range 𝑚 ★ ∼ 100 -2 × 10 5 M ⊙ with a complex shape caused by the multimodal distribution of / 𝑀 in . There are two low-mass peaks around 𝑚 ★ ∼ 160 and 500 M ⊙ and a highmass peak at 𝑚 ★ ∼ 10 5 M ⊙ that contains a small fraction ( ∼ 7%) of stars. In between the second low-mass peak and the high mass peak, the distribution approximately follows a broken power-law form with 𝑑𝑁 / 𝑑 log 𝑚 ★ ∝ 𝑚 -1 . 7 ★ at 𝑚 ★ ∼ 500 -2000 M ⊙ and 𝑑𝑁 / 𝑑 log 𝑚 ★ ∝ 𝑚 -0 . 4 ★ at 𝑚 ★ ∼ 5000 -7 × 10 4 M ⊙ . Although the majority ( ∼ 80%) of primordial star-forming clouds in high𝑧 luminous quasar progenitor haloes form ordinary Pop III stars with 𝑚 ★ ∼ 100 -2000 M ⊙ (in the feedback-regulated regime), a non-negligible fraction ( 𝑁 SMS ∼ 14%) of clouds produce SMSs with 𝑚 ★ > 10 4 M ⊙ that can become  \n5  \nFigure 8. Mass distribution of Pop III stars formed in high𝑧 luminous quasar progenitor haloes based on the / 𝑀 in data (Fig. 7) from Li et al. (2021) for 𝑣 bc = 0 (solid), 1 (dashed), 2 (dash-dotted), 3 𝜎 rms (dotted), and the cosmic average (shaded), assuming a single star per cloud/halo. The distribution approximately follows a (broken) power-law form in the range 𝑚 ★ ∼ 500 -7 × 10 4 M ⊙ for 𝑣 bc = 0 and the cosmic average, as indicated by the longdashed lines. The mass distribution of typical Pop III stars formed in unbiased regions with 𝑣 bc = 0 from Hirano et al. (2015) is shown with the thick dashed curve for comparison. These distributions are normalized such that the total number of haloes is 10000 in each case.  \n<!-- image -->  \nheavy BH seeds. Besides, a small fraction ( ∼ 8%) of clouds produce VMSswith 𝑚 ★ ∼ 2000 -9000 M ⊙ that can explain the high nitrogen abundances observed in several high𝑧 galaxies by JWST (Nandal et al. 2024b).  \nSimilarly to our case with 𝑣 bc = 0, the mass distribution predicted by Hirano et al. (2015) is also dominated by ordinary Pop III stars ( 𝑚 ★ ∼ 100 -1000 M ⊙ ) with a peak at 𝑚 ★ ∼ 300 M ⊙ . However, it has a sharp cutoff at ∼ 1000 M ⊙ in the lack of enhanced environmental effects (halo mergers and LW radiation) that can produce high / 𝑀 in and more massive Pop III stars in high𝑧 luminous quasar progenitor haloes. Hirano et al. (2015) also obtain a small ( ∼ 12%) fraction of Pop III stars below 100 M ⊙ , mostly formed in HD cooling haloes at 𝑧 ≲ 15, while such cases are extremely rare in high𝑧 luminous quasar progenitor haloes, since Pop III star formation mostly occur at 𝑧 ∼ 15 -50 in the relevant biased regions (Li et al. 2021). The masses of Pop III stars are significantly boosted by higher streaming velocities between baryons and dark matter, such that most ( ≳ 90%) stars have 𝑚 ★ > 10 4 ( 5 ) M ⊙ for 𝑣 bc = 1 ( 2 ) 𝜎 rms , while the number of VMSs is reduced. With the contributions of such high𝑣 bc regions, the cosmic average mass function of Pop III stars is dominated ( ∼ 70%) by SMSs, showing a very top-heavy power law 𝑑𝑁 / 𝑑 log 𝑚 ★ ∝ 𝑚 0 . 5 ★ at 𝑚 ★ ∼ 3000 -7 × 10 4 M ⊙ followed by two higher-mass peaks that account for half of the entire population. These differences indicate that the environmental effects at halo and cosmic scales strongly regulate the masses of Pop III stars. In general, the broad mass distributions predicted by our analytical model are consistent with the picture emerging from recent numerical simulations that the mass spectrum of BH seeds from the first stars is a continuum (Li et al. 2023; Regan & Volonteri 2024).  \nGiven the number density 𝑛 h ∼ 2 . 3 × 10 -6 cMpc -3 of typical high𝑧 luminous quasar host haloes (with 𝑀 h ∼ 10 12 M ⊙ at 𝑧 ∼ 6), the number density of heavy ( ≳ 10 4 M ⊙ ) BH seeds from Pop III SMSs formed in their main progenitors is 𝑛 seed = 𝑁 SMS 𝑛 h ∼ 1 . 6 × 10 -6 cMpc -3 (averaging over regions with different 𝑣 bc ). This should be regarded as a lower limit of the abun-  \nTable 1. Statistics of Pop III clusters and stars in high𝑧 luminous quasar progenitor haloes. Column 1: streaming motion velocity, where 'average' denotes the cosmic average case constructed by numerical integration based on Eq. 24. Columns 2-4: multiplicity parameters (Sec. 2.3). Column 5: fraction of haloes hosting massive clusters with 𝑀 ★ > 10 4 M ⊙ . Column 6: average number per halo of VMSs in the mass range 𝑚 ★ ∼ 2000 -9000 M ⊙ as sources of large nitrogen yields (Nandal et al. 2024b). Column 7: average number per halo of SMSs with 𝑚 ★ > 10 4 M ⊙ as progenitors of heavy BH seeds. Column 8: average number per halo of extremely massive stars above 10 5 M ⊙ whose BH remnants are able to grow rapidly (Regan & Volonteri 2024).  \nluminous quasar progenitor haloes, as Pop III formation can occur in multiple branches of a merger tree beyond the main branch. Nevertheless, this lower limit already accounts for ∼ 11% of the number density 𝑛 quasar ∼ 1 . 4 × 10 -5 cMpc -3 of luminous quasars with bolometric luminosities 𝐿 bol ≳ 10 46 erg -1 at 𝑧 ∼ 6 inferred from JWST observations 12 (Akins et al. 2024, see their fig. 8). Given Eddington ratios and duty cycles of the order of unity (Fragione & Pacucci 2023), such heavy BH seeds from Pop III SMSs (formed at 𝑧 ≳ 15) can grow to 𝑚 BH ≳ 10 8 M ⊙ by 𝑧 ∼ 6, powering luminous quasars with 𝐿 bol ∼ 𝐿 Edd ≃ 1 . 2 × 10 38 erg s -1 ( 𝑚 BH / M ⊙ ) ≳ 10 46 erg -1 . In fact, about half of the Pop III SMSs (and their BH remnants) have masses above 10 5 M ⊙ , which is the minimum BH mass required to activate efficient growth found in recent simulations (Regan & Volonteri 2024). A detailed investigation for the formation (e.g., Haemmerlé 2021, 2024; Cammelli et al. 2024; Coughlin & Begelman 2024; Nagele & Umeda 2024; Nandal et al. 2024a; Saio et al. 2024; Shibata et al. 2024) and growth (e.g., Johnson & Haardt 2016; Latif & Khochfar 2020; Toyouchi et al. 2021; Inayoshi et al. 2022a; Hu et al. 2022; Park et al. 2022; Li et al. 2023; Lupi et al. 2024; Massonneau et al. 2023; Sassano et al. 2023; Bhowmick et al. 2024; Jeon et al. 2023, 2024) of BH seeds from Pop III SMSs is beyond the scope of  \nthis paper. In general, Pop III SMSs can make a significant contribution to the (heavy) BH seed population underlying high𝑧 luminous quasars, at least under the single-star assumption.  \nNext, we explore the effects of multiplicity at the smaller scales of star-forming clouds. The 12 multiplicity models shown in Sec. 2.3 and listed in Table 1 (with 𝑁 > 1) are applied to the / 𝑀 in datasets from Li et al. (2021). For conciseness, we focus on the dependence on 𝑁 as the most important multiplicity parameter, considering the cosmic average case for 𝑣 bc . The mass distributions of Pop III clusters and stars with 𝑁 = 3, 10, 30, and 100 given 𝑓 2 = 0 . 8 and 𝛼 = 2 are shown in Fig. 9. The corresponding statistics are summarized in Table 1, which also include the results for the other 8 models varying 𝑓 2 and 𝛼 . As expected, fragmentation generally suppresses the formation of massive Pop III clusters with 𝑀 ★ ≳ 10 4 M ⊙ , reducing their number fraction from 𝑓 MS ∼ 70% in the single-star model to 𝑓 MS ∼ 25 -32% for 𝑁 ∼ 3 -100, as the full bloating phase becomes rare in the presence of low-mass slowly-accreting stars. However, in the 𝑁 = 3 case, the maximum total stellar mass increases to 𝑀 ★ ∼ 4 × 10 5 M ⊙ because the collapse of the most massive star by GRI is delayed. Interestingly, the Pop III cluster masses predicted by our model remain below the upper limit 𝑀 ★, max ∼ 4 × 10 5 -7 × 10 5 M ⊙ (Liu & Bromm 2020b) estimated from the non-detection of Pop III-dominated galaxies in the Hubble Frontier Fields at 𝑧 ∼ 6 -9 with a limiting rest-frame UV magnitude of -13.5 (Bhatawdekar & Conselice 2021) using the stellar population synthesis code /y.pc/g.pc/g.pc/d.pc/r.pc/a.pc/s.pc/i.pc/l.pc (Zackrisson et al. 2011). The decrease of 𝑓 MS with 𝑁 slows down at 𝑁 ≳ 10 and is even reversed for 𝑁 ≳ 30, because the ionizing power is lower when there are more (low-mass) stars for a fixed total stellar mass in the no-bloating feedback-regulated regime. In fact, there is always a significant fraction 𝑓 MC ∼ 22 -34% of haloes hosting massive clusters for the 12 models considered here covering 𝑁 ∼ 3 -100, 𝑓 2 ∼ 0 . 2 -0 . 8, and 𝛼 ∼ 2 . 5 -0 . 5. This indicates that massive star clusters are common outcomes of primordial star formation in over-dense regions, whose 𝑁 -body dynamics has crucial implications on the formation and evolution of merging Pop III BBHs and their gravitational wave signatures (Wang et al. 2022; Liu et al. 2023; Mestichelli et al. 2024).  \nThe number of SMSs is also reduced as 𝑁 increases, showing a sharp drop from 𝑁 = 3 to 𝑁 = 10, down to 𝑁 SMS ∼ 0 . 02 per halo for 𝑁 = 100. The corresponding number density of heavy BH seeds is 𝑛 seed ≃ 9 . 8 × 10 -7 ( 4 . 5 × 10 -8 ) cMpc -3 for 𝑁 = 3 ( 100 ) , lower than the observed number density of 𝑧 ∼ 6 luminous ( 𝐿 bol ≳ 10 46 erg s -1 ) quasars by a factor of ∼ 15 ( 300 ) . Besides, given 𝑁 ≳ 3, most ( ≳ 80%) SMSs have 𝑚 ★ ≲ 10 5 M ⊙ , which may not be massive enough to produce rapidly growing BH seeds (Regan & Volonteri 2024). This indicates that fragmentation can significantly reduce the efficiency of heavy BH seed formation in Pop III starforming clouds. On the other hand, formation of VMSs with 𝑚 ★ ∼ 2000 -9000 M ⊙ as potential sources of strong nitrogen enrichment (Nandal et al. 2024b) is significantly enhanced by fragmentation. The number of such VMSs is 𝑁 VMS ∼ 0 . 5 -0 . 9 per halo for 𝑁 ∼ 3 -100 compared with 𝑁 VMS ∼ 0 . 06 for 𝑁 = 1. This implies that understanding the effects of fragmentation/multiplicity on the final products of primordial star formation is essential for interpreting the peculiar observational signatures of high𝑧 galaxies/quasars unveiled by JWST (Adamo et al. 2024). Note that here we focus on the natal population of stars/BHs, while VMSs, SMSs, and heavy BH seeds can further form and grow via run-away collisions in massive clusters (e.g., Sakurai et al. 2017; Boekholt et al. 2018; Wang et al. 2022). Last but not least, when 𝑁 is larger, the shape of the overall IMF of Pop III stars is increasingly influenced by the protostar mass distribution assumed in the multiplicity model and extends to lower masses (down to 𝑚 ★ ∼ 1 M ⊙ for 𝑁 ≳ 30). The mass distribution of  \n<!-- image -->  \nFigure 9. Mass distributions of Pop III clusters (top) and stars (bottom) formed in high𝑧 luminous quasar progenitor haloes for 𝑁 = 3 (thin solid), 10 (dashed), 30 (dash-dotted), and 100 (dotted) with 𝑓 2 = 0 . 8 and 𝛼 = 2. Here, the cosmic average distributions are constructed from the results with different values of 𝑣 bc using Eq. 24. The mass distribution in the single-star model ( 𝑁 = 1) is shown with the shaded region. For 𝑁 ≥ 10, the mass distribution of stars has a power-law shape similar to that assumed in the cloud/disc-scale multiplicity model ( 𝛼 = 2) in a mass range that is broader for a larger 𝑁 , as illustrated by the long-dashed line. The mass distribution of typical Pop III stars formed in unbiased regions with 𝑣 bc = 0 from Hirano et al. (2015) is shown with the thick dashed curve, which is normalized to have a total number of haloes of 10000 consistent with the other distributions.  \n<!-- image -->  \nstars follows a power-law form 𝑑𝑁 / 𝑑 log 𝑚 ★ ∝ 𝑚 -1 ★ consistent with the input IMF at the cloud/disc scales ( 𝛼 = 2) from 𝑚 ★ ∼ 5000 M ⊙ to 𝑚 ★ ∼ 500 / 100 / 50 M ⊙ for 𝑁 = 10 / 30 / 100, while in the 𝑁 = 3 case, the mass distribution of stars resembles more the cluster mass distribution regulated by larger-scale physics.  \nCompared with 𝑁 , the effects of 𝑓 2 and 𝛼 are weaker. The most obvious trend is that 𝑁 SMS decreases when 𝑓 2 is lower (see the third section in Table 1). Indeed, in our default multiplicity model, beyond 𝑁 , it is mainly 𝑓 2 (rather than 𝛼 ) that characterizes the oligarchy of the most massive objects from competitive accretion, which is important for SMS formation. It will be shown in Appendix B that 𝛼 can have stronger effects under an alternative parametrisation of multiplicity.  \nWe end this section with some brief comments on the Pop III SFE in high𝑧 luminous quasar progenitor haloes, defined as 𝑓 ★ ≡ 𝑀 ★ / 𝑀 gas given the total gas mass 𝑀 gas in the halo at the moment of star formation. A broad range 𝑓 ★ ∼ 0 . 001 -0 . 1 is obtained regardless of the choices of multiplicity parameters. In the single-star model, the  \nefficiency is generally higher in more massive haloes with larger 𝑀 ★ . The typical values are 𝑓 ★ ∼ 0 . 004 for ordinary Pop III star formation ( 𝑀 ★ < 2000 M ⊙ ) in minihaloes with 𝑀 h ∼ 2 × 10 5 -3 × 10 6 M ⊙ and 𝑓 ★ ∼ 0 . 02 for SMS formation ( 𝑀 ★ > 10 4 M ⊙ ) in atomic-cooling haloes ( 𝑀 h ≳ 10 7 M ⊙ ). When fragmentation is included, 𝑓 ★ hardly changes in minihaloes but is significantly reduced for massive Pop III clusters ( 𝑀 ★ > 10 4 M ⊙ ) given 𝑁 ≥ 10, where the cosmic average efficiency is ∼ 0 . 008. We plan to conduct more detailed analyses on Pop III SFE and its correlation with IMF in future work by applying our model to a cosmologically representative sample of Pop III starforming haloes.", '3.2 Caveats and future improvements': "The goal of this paper is to provide a flexible analytical framework that can predict the outcomes of Pop III star formation from the properties of star-forming clouds and clarify the key physical processes that determine the final masses of Pop III stars. As an initial version, our model is kept as simple as possible, while reproducing the key results of small-scale simulations of primordial star-forming clouds (Hirano et al. 2014; Sugimura et al. 2023; Toyouchi et al. 2023). It does not necessarily take into account all possible factors governing Pop III star formation, and some model elements can be improved to better reflect the complexity in simulations and reality, as we discuss below.  \n(i) The scaling relations used in our model (Sec. 2.1) describe the (smoothed) median evolution for a large sample of Pop III protostar systems focusing on the surviving protostars. However, the evolution of individual systems shows strong stochasticity and diversity in simulations. For instance, the accretion rates of protostars fluctuate significantly at very short timescales for Pop III star-forming clouds across a broad range of initial collapse rates / 𝑀 in ∼ 10 -4 -1 M ⊙ yr -1 , modulated by the interactions between protostars and the clumpy structure of the unstable star-forming disc (Sakurai et al. 2016; Hosokawa et al. 2016; Chon & Omukai 2020; Park et al. 2023, 2024; Sugimura et al. 2020, 2023; Toyouchi et al. 2023; Prole et al. 2024; Reinoso et al. 2023). This has interesting implications on the bloating mode of Pop III star formation (Sec. 2.2). For typical protostars in this regime with masses ∼ 1000 -10 5 M ⊙ and KH timescales 𝑡 KH ∼ 100 -1000 yr, the bloating phase will terminate if the accretion rate stays below the critical rate / 𝑀 ★, crit ∼ 0 . 01 -0 . 04 M ⊙ yr -1 for a duration longer than ∼ 10 -100 𝑡 KH ∼ 10 3 -10 5 yr (Schleicher et al. 2013; Sakurai et al. 2015; Toyouchi et al. 2023; Reinoso et al. 2023) in which the star contracts to produce intense ionizing radiation. Such contraction is more likely to occur at later stages as the accretion rate generally decreases with time. In the current model, we simply assume that the contraction occurs exactly when the (smoothed) median accretion rate drops below / 𝑀 ★, crit , while in reality, the timing and the resulting final mass may have a large scatter due to the stochastic nature of this process. We plan to consider such fluctuating accretion rates as well as protostar mergers and ejections in future work based on statistical analyses of the protostar accretion/merger/ejection histories in simulations and/or physically motivated models for fragmentation via disc instability (e.g., Tsukamoto et al. 2015; Inoue & Yoshida 2020; Kimura et al. 2021). In this way, we can explicitly follow the mass evolution of individual protostars and their circumstellar discs in a fully dynamical manner.  \n(ii) Weassume that Pop III star formation happens via monolithic collapse with weak turbulence based on the simulations by Chon et al. (2021) for isolated primordial star-forming clouds with masses 𝑀 c ∼ 2000 M ⊙ (see also Higashi et al. 2022). However, this as-  \nn may not always hold in realistic environments of Pop III star formation, especially for more massive clouds with violent collapse ( / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 ). The presence of strong turbulence in star-forming clouds has two consequences: First, gas inflows onto the star-forming disc can be highly dynamical, deviating from the steady declining trend described by Eq. 5 (Regan et al. 2020a, see their fig. 6). Second, core fragmentation can break the cloud into multiple clumps each hosting their own star-forming discs and star clusters, which is particularly relevant when strong streaming motion between baryons and dark matter produces filamentary clouds (Hirano et al. 2018, 2023; Latif et al. 2022a). These effects are expected to reduce the final masses of Pop III stars formed. It is feasible to include turbulent fragmentation in our framework like what has been done for present-day star formation (e.g., Larson 1973; Hennebelle & Chabrier 2008; Hopkins 2012; Guszejnov & Hopkins 2015; Thomasson et al. 2024), but a larger sample of simulations is required to better understand the statistical properties of turbulent fragmentation under various conditions of primordial star formation.  \n(iii) Weadopt a simple power-law form with an invariant slope for the mass function of protostars, in the hope that it is flexible enough to qualitatively demonstrate the effects of multiplicity. Our model is not meant to reflect the reality, where the time evolution and shape of protostar mass function can be more complex, especially for SMS formation with high gas inflow rates ( / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 ). Since high accretion rates ( > / 𝑀 ★, crit ∼ 0 . 01 -0 . 04 M ⊙ yr -1 , Hosokawa et al. 2013; Haemmerlé et al. 2018; Herrington et al. 2023; Nandal et al. 2023) cause significant expansion of some protostars during such violent collapse, super-competitive accretion and stellar collisions can cause a 'class bifurcation' in the stellar population, where the total stellar mass is dominated by a small number of SMSs accompanied by a large number of low-mass stars (e.g., Chon & Omukai 2020; Prole et al. 2024; Reinoso et al. 2023; Schleicher et al. 2023). The exact form of the mass function of Pop III (proto)stars (at the cloud scale) is still in debate, as different simulations produce highly divergent results which often suffer from poor statistics with a small number of realisations. We defer the consideration of more complex mass distributions to future work.  \n(iv) Related to the previous point (iii), the dependence of fragmentation outcomes on environmental factors (e.g., 𝑣 bc , LW radiation, and halo dynamics) and cloud/disc properties (e.g., / 𝑀 in , 𝑀 c , 𝐾 , 𝛾 eff , and 𝜂 , or even the underlying temperature, density, and velocity profiles) regulated by larger (cosmic/halo)-scale physics is also ignored in our calculation of the ˆ 𝑀 ★ -/ 𝑀 in relation. Recent simulations start to find hints for the existence of a non-trivial dependence (see, e.g., fig. 4 in Regan & Volonteri 2024). For instance, it is shown by Saavedra-Bastidas et al. (2024, see their fig. 7) that 𝑓 2 can be accurately predicted from the (cloud-scale) SFE 𝜖 ★ ≡ ˆ 𝑀 ★ / 𝑀 c using a machine learning model trained on existing simulation results, which implies that 𝑓 2 strongly depends on / 𝑀 in . The parameters for cloud/disc properties may also correlate with each other, while they are treated as independent variables in our case. One needs to analyse a large set of (multi-scale) simulations systematically covering the diverse environments of Pop III star formation to fully characterize the correlations between cloud/disc properties and multiplicity parameters and their dependence on the larger-scale condition of primordial star formation. Here, we simply estimate the range of possible effects of such correlations on the final masses of Pop III star clusters with controlled numerical experiments (see Sec. 2.3 and Appendix A).  \n(v) Our model does not include the effects of magnetic fields, which have been intensively investigated in recent numerical studies (e.g., McKee et al. 2020; Sharda et al. 2020, 2021; Hirano et al. 2021; Prole et al. 2022b; Saad et al. 2022; Stacy et al. 2022; Hirano &  \nMachida2022;Higashi et al. 2024; Sadanari et al. 2023, 2024; Sharda &Menon 2024). Currently, (radiative) magneto-hydrodynamic simulations have converged on the picture that magnetic fields are significantly amplified by a turbulent dynamo and rotation to near equipartition with kinetic energy in the central region of primordial starforming clouds, regardless of the initial field strength. Nevertheless, there is no census on the roles played by such strong magnetic fields during the star formation process, which may depend on the initial cloud properties. For instance, Sharda & Menon (2024, see their fig. 7) find that dynamically strong magnetic fields slow down accretion onto protostars in typical primordial star-forming clouds with 𝑀 c ∼ 1000 M ⊙ . However, the simulations by Hirano et al. (2021) targeting more massive ( 𝑀 c ∼ 2 × 10 6 M ⊙ ), atomic-cooling clouds show that efficient extraction of angular momentum by magnetic fields increases the accretion rate and enhances fragmentation and protostar mergers, facilitating the growth of the primary protostar. Besides, it is shown by several groups (Sharda et al. 2020; Saad et al. 2022; Stacy et al. 2022; Hirano & Machida 2022; Sadanari et al. 2024) that magnetic pressure and torques can efficiently suppress and even eliminate fragmentation, while Prole et al. (2022b) find that the number and total mass of protostars are unaffected by saturated magnetic fields at equipartition. It is non-trivial to model magnetic fields explicitly in a simple analytical manner. However, their effects can possibly be absorbed into some existing parameters (e.g., / 𝑀 in , 𝐾 , 𝛾 eff , 𝜂 , and 𝑁 ) if clearer trends can be deduced from simulations, which is an interesting topic for follow-up studies.", '4 SUMMARY AND OUTLOOK': 'We build an analytical model for Pop III star formation to predict the final masses of Pop III clusters and stars from the (initial) properties of star-forming clouds, using physically motivated prescriptions for the following physical processes and quantities:  \n- (1) gas inflow onto the star-forming disc,\n- (2) disc geometry and fragmentation, the resulting spatial, mass, and accretion rate distributions of protostars,\n- (3) destruction of the star-forming disc/cloud by stellar feedback,\n- (4) evolution of protostars regulated by their accretion histories.  \nThe first two processes are expected to be closely related to the cloud-scale (initial) condition of star formation, while the last two processes lead to the termination of star formation (via depletion of gas supply or stellar collapse/explosion). Any self-consistent model of star formation must consider these four aspects. The basic idea of our model is to express them as functions of the total mass of protostars, so that the final mass can be easily obtained in a rootfinding process.  \nTo be specific, given the initial gas inflow rate / 𝑀 in onto the starforming disc, we use power-law scaling relations between mass, size, and time to describe the (smoothed) median evolution of Pop III protostar systems (Sec. 2.1), which is combined with a phenomenological fragmentation model (see Sec. 2.3 and Appendix A) to handle points (i) and (ii), based on the results of (magneto-) hydrodynamic simulations (Liu et al. 2021b; Sharda et al. 2020). For point (iii), we calculate the disc photo-evaporation rate as a function of the production rates of ionizing photons from individual protostars (Tanaka et al. 2013) considering different spatial configurations of protostars and their circumstellar discs from point (ii). Besides, for point (iv), we include the bloating phase of protostar evolution with negligible ionizing power, triggered by high accretion rates ( ≳ / 𝑀 ★, crit ∼ 0 . 01 -0 . 04 M ⊙ yr -1 , Hosokawa et al. 2013; Haem-  \nHerrington et al. 2023; Nandal et al. 2023), and consider the finite lifetimes and maximum masses of stars before collapse/explosion (see Sec. 2.2).  \nDespite its simplicity, our model covers the full range of outcomes of Pop III star formation known to date, from ordinary small ( 𝑀 ★ ∼ 10 -2000 M ⊙ ) clusters in molecular-cooling clouds to massive ( 𝑀 ★ ≳ 10 4 M ⊙ ) clusters containing supermassive ( 𝑚 ★ ∼ 10 4 -3 × 10 5 M ⊙ ) stars (SMSs) under violent collapse ( / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 ) of atomic-cooling clouds, reproducing relevant simulation results (Hirano et al. 2014; Sugimura et al. 2023; Toyouchi et al. 2023). Moreover, our model specifically considers the interplay between feedback and fragmentation. We find that enhanced fragmentation tends to increase the final mass of the star cluster regulated by photo-ionization feedback under low gas inflow rates ( / 𝑀 in ≲ 0 . 1 M ⊙ yr -1 ). However, it suppresses the formation of SMSsinhigher-/ 𝑀 in clouds not only by internal starvation (Prole et al. 2022a) but also due to the decrease of final cluster mass caused by enhanced photo-ionization feedback from low-mass slowly-accreting ( ≲ / 𝑀 ★, crit ) stars.  \nAs an example, we apply our model to the Pop III star-forming clouds in the main progenitors of typical haloes hosting high𝑧 luminous quasars with 𝑀 h ∼ 10 12 M ⊙ at 𝑧 ∼ 6 (Li et al. 2021), which predicts broad mass distributions for Pop III clusters and stars extending to a few 10 5 M ⊙ (Sec. 3.1). This shows that the formation of Pop III massive clusters is common (among ∼ 20 -70%ofthe Pop III star-forming haloes) in such biased ( ∼ 4 𝜎 ) regions, which favours the production of dynamical Pop III BBH mergers with distinct gravitational wave signatures (Wang et al. 2022; Liu et al. 2023; Mestichelli et al. 2024). The corresponding SFE (ratio of the total stellar and gas masses in the halo) covers a broad range 𝑓 ★ ∼ 0 . 001 -0 . 1, and is generally higher for more massive clusters. As the number of stars 𝑁 per cloud increases, the shape of the overall stellar mass distribution becomes more similar to that of the protostar mass distribution assumed in the fragmentation model (for individual clouds), while it closely follows the cluster mass distribution when fragmentation is inefficient ( 𝑁 ≲ 3). This indicates that the overall IMF of Pop III stars is regulated by both large-scale physics (embodied by the gas inflow rate distribution for Pop III star-forming clouds) and small-scale fragmentation/feedback processes. When fragmentation is inefficient ( 𝑁 ≲ 3), the resulting Pop III SMSs as progenitors of heavy ( ≳ 10 4 M ⊙ ) BH seeds can account for a significant fraction ( ≳ 7%) of observed luminous ( ≳ 10 46 erg s -1 ) quasars at 𝑧 ∼ 6 (e.g., Greene et al. 2024; Matthee et al. 2024; Kokorev et al. 2023; Akins et al. 2024) powered by SMBHs with 𝑚 BH ≳ 10 8 M ⊙ . On the other hand, if we extend this trend to present-day star formation in a metalenriched medium (by increasing 𝑁 and the slope 𝛼 of the protostar mass distribution), where cooling and fragmentation are much more efficient compared with the Pop III case, the overall IMF is mostly shaped by small-scale processes rather than large-scale physics, such that it becomes nearly universal as long as the IMF in individual star clusters is relatively insensitive to cloud/cluster properties 13 .  \nFinally, we also discuss the potential missing pieces and directions to improve the model, such as stochastic evolution, core fragmenta-  \n13 To what extent (or at which scale) the IMF of present-day stars remains universal is still in debate. Recent observations and simulations have discovered hints for a non-trivial dependence of the IMF on various properties of the star-forming environment (e.g., metallicity, temperature, density, star formation rate, and cloud mass, Marks et al. 2012; Jeřábková et al. 2018; Lacchin et al. 2020; Gunawardhana et al. 2011; Cheng et al. 2023; Dib 2023; Grudic et al. 2023; Rusakov et al. 2023; Yan et al. 2023; Tanvir & Krumholz 2024).  \ntion, complex forms of protostar mass distribution, correlations between cloud/disc properties and multiplicity parameters, their dependence on the larger-scale condition of primordial star formation, and the effects of magnetic fields (Sec. 3.2). In conclusion, our cloud/discscale analytical model for Pop III star formation can be incorporated into semi-analytical models (e.g., Manrique et al. 2015; Griffen et al. 2018; Dayal et al. 2020; Visbal et al. 2020; Li et al. 2021; Lupi et al. 2021; Hartwig et al. 2022; Hegde & Furlanetto 2023; Nebrin et al. 2023; Bovill et al. 2024; Ventura et al. 2024; Feathers et al. 2024; Trinca et al. 2024) and cosmological hydrodynamic simulations (e.g., Johnson et al. 2013; Smith et al. 2015; Skinner & Wise 2020; Kulkarni et al. 2021; Schauer et al. 2021; Kiyuna et al. 2023; Yajima et al. 2022; Garcia et al. 2023; Lenoble et al. 2024; Sugimura et al. 2024) that focus on larger (cosmic/halo)-scale physics (e.g., structure formation, radiation backgrounds, baryons-dark matter streaming motion, and metal enrichment) to self-consistently follow the formation and feedback of Pop III stars across all scales in a cosmologically representative volume. Such self-consistent modelling is necessary to make comprehensive predictions on the observational signatures and imprints of Pop III stars and their BH remnants, which can be directly compared with multi-messenger observations of Cosmic Dawn to validate and advance our theories of early star/galaxy/BH formation.', 'ACKNOWLEDGEMENTS': 'The authors thank the anonymous referee for insightful comments and helpful suggestions. We thank Wenxiu Li for providing the data underlying Li et al. (2021). We thank Piyush Sharda for providing the data for protostar accretion histories from the simulations shown in Sharda et al. (2020). We thank Daisuke Toyouchi for useful discussions and Anastasia Fialkov for helpful comments on the manuscript. BL gratefully acknowledges the support of the Royal Society University Research Fellowship. Research at the Perimeter Institute for Theoretical Physics is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. The development of this paper benefited from in-person discussions at the Dark Matter, First Light workshop hosted by the Perimeter Institute for Theoretical Physics and the First Star VII conference organized by the Flatiron Institute and Columbia University. We thank the organizers of these events for providing us platforms to exchange ideas. This work excessively used the public packages numpy (van der Walt et al. 2011), matplotlib (Hunter 2007), and scipy (Virtanen et al. 2020). The authors wish to express their gratitude to the developers of these packages and to those who maintain them.', 'DATA AVAILABILITY': 'The code and data underlying this paper will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'Abbott R., et al., 2020a, Phys. Rev. Lett., 125, 101102 Abbott R., et al., 2020b, ApJ, 900, L13 Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93  \nAdamo A., et al., 2024, arXiv e-prints, p. arXiv:2405.21054  \n1756  \nf  \n2 = 0.8,  \n=2 (concentrated)  \nFigure A1. Map of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in space with other parameters fixed as in Fig. 4. The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \nYajima H., Abe M., Fukushima H., Ono Y., Harikane Y., Ouchi M., Hashimoto T., Khochfar S., 2023, MNRAS, 525, 4832  \nYan Z., Jerabkova T., Kroupa P., 2023, A&A, 670, A151  \nYang J., et al., 2021, ApJ, 923, 262  \nYoon S. C., Dierks A., Langer N., 2012, A&A, 542, A113  \nYoshida N., Abel T., Hernquist L., Sugiyama N., 2003, ApJ, 592, 645  \nYoshida T., Umeda H., Maeda K., Ishii T., 2016, MNRAS, 457, 351  \nZackrisson E., Rydberg C.-E., Schaerer D., Östlin G., Tuli M., 2011, ApJ, 740, 13  \nZackrisson E., et al., 2012, MNRAS, 427, 2212  \nZackrisson E., et al., 2024, MNRAS, 533, 2727  \nvan der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science  \nand Engineering, 13, 22', 'APPENDIX A: DETAILED RESULTS IN A LARGER PARAMETER SPACE': 'In Sec. 2.3, we demonstrate the dependence of our results on mutiplicity parameters with select examples in the fiducial case where 𝑁 is independent of 𝑀 in the concentrated feedback scenario (Eq. 22) given / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 , 𝑅 pe ∝ 𝑀 1 . 25 , and 𝛾 eff = 1 . 09. In this Appendix, we first provide more detailed results in the fiducial case as maps of ˆ 𝑀 ★ in Figs. A1-A3 covering the curves shown in Figs. 4-6. Then we vary the feedback and fragmentation prescriptions (Appendix A1) and the parameters governing the mass-size scaling relation and feedback strength (Appendix A2) to evaluate their impact on our results.', 'A1 Feedback and fragmentation prescriptions': 'In this section, we relax the assumptions of concentrated feedback (Eq. 22) and constant 𝑁 . As an extreme opposite to concentrated feedback, we consider the scenario where the distances between protostars are comparable to the disc size caused by outward migrations from accretion of angular momentum, as seen in recent radiative hydrodynamic simulations (e.g., Sugimura et al. 2020, 2023; Park et al. 2023, 2024), so the ionizing flux from a protostar mostly affects its own circumstellar disc, and the total mass loss rate of photoevaporation is the summation of the mass loss rates from the sub-discs  \nN  \n=10,  =1 (concentrated)  \nFigure A2. Map of ˆ 𝑀 ★ in the 𝑓 2 -/ 𝑀 in space with other parameters fixed as in Fig. 5. The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \nFigure A3. Map of ˆ 𝑀 ★ in the 𝛼 -/ 𝑀 in space with other parameters fixed as in Fig. 6. The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \naround individual protostars:  \n/ 𝑀 pe ≃ 0 . 015 M ⊙ yr -1 𝑁 ∑︁ 𝑖 GLYPH<26> GLYPH<18> / 𝑞 ion ,𝑖 10 52 s -1 GLYPH<19> GLYPH<20> 𝑅 ( 𝑚 𝑖 ) 10 4 AU GLYPH<21> GLYPH<27> 1 / 2 . (A1)  \nHere, we assume that the sub-discs of individual protostars are always connected to each other (by spiral arms) so that the termination of protostar growth happens spontaneously when / 𝑀 pe exceeds / 𝑀 . The maximum sub-disc radius 𝑅 max ,𝑖 of each protostar 𝑖 is chosen to ensure that the total disc area does not exceed 𝜋𝑅 2 c :  \n𝑅 max ,𝑖 = 𝑅 c 𝑓 𝛿 / 𝛽 𝑖 / " 𝑁 ∑︁ 𝑖 𝑓 2 𝛿 / 𝛽 𝑖 # 1 / 2 , (A2)  \nassuming that the size-mass scaling relation (Eq. 4) holds for each protostar 14 . The photo-evaporation rate / 𝑀 pe is smaller in this scenario of distributed feedback due to the super-linear scaling relation between mass and size (Eq. 9).  \nBesides, it is shown by Susa (2019) and Liu et al. (2021a) that  \nFigure A4. Same as Fig. 4 but for four cases combining the treatments of 𝑁 and / 𝑀 pe given 𝑁 ( 𝑁 0 ) = 3, 𝑓 2 = 0 . 8, and 𝛼 = 1. The solid (dash-dotted) and dashed (dotted) curves show the results without (with) time evolution of 𝑁 under concentrated and distributed feedback, respectively.  \n<!-- image -->  \nthe (median) number of surviving protostars increases with time following a power-law scaling 𝑁 ∝ 𝑡 0 . 3 at least in the early stage 15 in most simulations of Pop III star-forming clouds. It is unclear whether this trend can hold till the end of star formation. In our default setup (Sec. 2.3), we assume that the number of surviving protostars has saturated at an input value (e.g., Shima & Hosokawa 2021; Chon et al. 2021) independent of 𝑀 before the accretion of protostars is turned down by feedback (i.e., 𝑡 frag < 𝑡 acc in the formalism of Liu et al. 2021a). Now, as an alternative model, we assume that 𝑁 keeps growing in time and introduce the dependence of 𝑁 on 𝑀 as  \n𝑁 = max h 1 , 𝑁 0 ( 𝑀 / 𝑀 0 ) 0 . 3 / 𝛽 i , (A3)  \nusing Eq. 8, where 𝑁 0 is the new parameter defined as the number of protostars at 𝑡 = 𝑡 0 . We keep 𝑓 2 and 𝑁 0 fixed to the input value, and let 𝑓 1 ( 𝑁, 𝑓 2 , 𝛼 ) evolve with 𝑀 through the dependence on 𝑁 , as fragmentation is expected to be important mainly for the smallest objects in the system.  \nFig. A4 shows the effects of the treatments for / 𝑀 pe and 𝑁 for the exemplar multiplicity model with 𝑁 ( 𝑁 0 ) = 3, 𝑓 2 = 0 . 8, and 𝛼 = 1. As expected, ˆ 𝑀 ★ increases when we consider distributed feedback, in particular for the intermediate regime, where the circumstellar discs of the most massive stars that undergo bloating are intact from photo-evaporation. On the other hand, allowing 𝑁 to increase with 𝑀 reduces ˆ 𝑀 ★ in the intermediate regime under concentrated feedback. However, the effect is minor for distributed feedback, because when the disc sizes approach their upper limits ( 𝑀 ≳ 2 × 10 4 M ⊙ ), enhanced fragmentation increases the ionizing flux from small protostars that do not enter the bloating phase but meanwhile reduce the disc sizes, and the two effects cancel out each other in the disc evaporation rate.  \nMore detailed results with varying 𝑁 ( 𝑁 0 ) in terms of maps of ˆ 𝑀 ★ in the 𝑁 ( 𝑁 0 ) -/ 𝑀 in space for the three additional combinations of fragmentation and feedback prescriptions considered here are shown in Figs. A5-A7, to be compared with the fiducial case (Fig. A1).  \nFigure A5. Map of ˆ 𝑀 ★ in the 𝑁 0 -/ 𝑀 in space for the 𝑁 evolution model in Eq. A3, with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \nf  \n2 = 0.8,  \n=2 (distributed)  \nFigure A6. Map of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in under distributed feedback (Eq. A1), with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->', 'A2 Mass-size scaling relation and feedback parameters': 'Next, we explore the effects of the mass-size scaling relation and feedback parameters on our results by varying the polytropic index 𝛾 eff of the effective EoS, the assumption on the disc photonevaporation scale 𝑅 pe , and / 𝑀 ★, crit for the single-star case ( 𝑁 = 1) and the exemplar multiplicity model with 𝑁 = 3 (constant), 𝑓 2 = 0 . 8, 𝛼 = 1, and concentrated feedback, as shown in Fig. A8. Here we consider three new models with 𝛾 eff = 1, 𝑅 pe = 10 4 AU, and / 𝑀 ★, crit = 0 . 01 M ⊙ yr -1 , varying one aspect each time with respect to the fiducial model ( 𝛾 eff = 1 . 09, 𝑅 pe ∝ 𝑀 𝛿 / 𝛽 , and / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 ). The 𝛾 eff = 1 model is meant to capture Pop III star formation under extremely strong LW radiation (with an intensity 𝐽 LW ≳ 1000 in units of 10 -21 erg s -1 Hz -1 , Sugimura et al. 2014) or the cosmic microwave background at very high redshifts ( 𝑧 ≳ 500, Ito & Omukai 2024), where H 2 formation is significantly suppressed and the gas remains quasi-isothermal during collapse by atomic cooling (see fig. 4 in Li et al. 2021).  \nIf the star-forming gas is isothermal with 𝛾 eff = 1, the accretion rate does not evolve with time and 𝑀 , so in the single-star case, there is a sharp transition from the feedback-regulated regime to the  \nf  \n2 = 0.8,  \n=2 (distributed)  \nFigure A7. Map of ˆ 𝑀 ★ in the 𝑁 0 -/ 𝑀 in space for the 𝑁 evolution model in Eq. A3 under distributed feedback (Eq. A1), with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \nFigure A8. Same as Fig. 4 but for four models with different parameters for the mass-size scaling relation and feedback strength. The solid curves show the results for the fiducial model with 𝛾 eff = 1 . 09, 𝑅 pe ∝ 𝑀 𝛿 / 𝛽 , and / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 , while we adopt 𝛾 eff = 1, 𝑅 pe = 10 4 AU, and / 𝑀 ★, crit = 0 . 01 M ⊙ yr -1 for the dashed, dash-dotted, and dotted curves, respectively, with the other parameters set to the fiducial values. In each case, the thin and solid curves show the results for the exemplar multiplicity model (with 𝑁 = 3, 𝑓 2 = 0 . 8, and 𝛼 = 1), and the single-star model ( 𝑁 = 1), respectively.  \n<!-- image -->  \nlifetime/GRI-regulated regime (where the bloating phase lasts until the the star collapses) at 𝜂 / 𝑀 in = / 𝑀 ★, crit . Therefore, ˆ 𝑀 ★ is higher for 𝛾 eff = 1 when the bloating phase is involved. However, ˆ 𝑀 ★ becomes slightly smaller in the feedback-regulated regime because the disc size is smaller when 𝛽 increases (Eq. 9) even though the accretion rate become higher (see Eq. 5). The trends are similar when multiple protostars are considered. Now there is a jump in the ˆ 𝑀 ★ -/ 𝑀 in curve each time the accretion rate of a protostar crosses / 𝑀 ★, crit such that it is in the bloating phase right before the end of star formation, which also occurs at a smaller / 𝑀 in compared with the 𝛾 eff = 1 . 09 case.  \nIn the single-star case and the exemplar multiplicity model, using 𝑅 pe = 10 4 AU instead of 𝑅 pe ∝ 𝑀 𝛿 / 𝛽 increases (reduces) ˆ 𝑀 ★ for ˆ 𝑀 in ≳ ( ≲ ) 0 . 006 M ⊙ yr -1 , while reducing / 𝑀 ★, crit facilitates the transition to the bloating or lifetime/GRI-regulated regime. With  \nFigure A9. Map of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in space for 𝛾 eff = 1 (instead of 𝛾 eff = 1 . 09) with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \nFigure A10. Map of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in space assuming a fixed disc photoevaporation scale 𝑅 pe = 10 4 AU (instead of the evolution model in Eq. 9) with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->  \n𝑅 pe = 10 4 AU, the effects of multiplicity on ˆ 𝑀 ★ vanish for ˆ 𝑀 in ∼ 0 . 01 -0 . 2 M ⊙ yr -1 in the feedback-regulated regime. These clouds are dominated by massive ( ≳ 200 M ⊙ ) stars that follow the linear relation between / 𝑄 ion and stellar mass, so the total production rate of ionizing photons and disc evaporation rate at a given total stellar mass 𝑀 hardly vary with 𝑁 at least for 𝑁 ≤ 100. Interestingly, the models with 𝑅 pe = 10 4 AUand / 𝑀 ★, crit = 0 . 01 M ⊙ yr -1 have similar results for / 𝑀 in ≳ 0 . 02 M ⊙ yr -1 in the single-star case, consistent with the simulation results in Toyouchi et al. (2023), where SMS formation starts at slightly lower / 𝑀 in compared with the fiducial model. This degeneracy is broken in the exemplar multiplicity model.  \nTo further illustrate the dependence on 𝑁 , Figs. A9, A10, and A11 show the maps of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in space for the three new models with 𝛾 eff = 1, 𝑅 pe = 10 4 AU, and / 𝑀 ★, crit = 0 . 01 M ⊙ yr -1 , respectively, in comparison with the results of the fiducial model in Fig. A1.  \nFigure A11. Map of ˆ 𝑀 ★ in the 𝑁 -/ 𝑀 in space for 𝑀 ★, crit = 0 . 01 M ⊙ yr -1 (instead of 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 ) with other parameters fixed as in Fig. 4 (to be compared with Fig. A1). The dotted contours denote 5 specific levels of ˆ 𝑀 ★ in units of M ⊙ .  \n<!-- image -->', 'APPENDIX B: ALTERNATIVE PARAMETRISATION OF MULTIPLICITY': 'In Sec. 2.3, we build a phenomenological multiplicity model with four variables 𝛼 , 𝑁 , 𝑓 1 , and 𝑓 2 , which are correlated with each other via Eq. 21, so only three of them can be chosen as free model parameters. Here, 𝛼 and 𝑁 have clear physical meanings as the slope of mass distribution and target number of stars and therefore, they should be adopted as model parameters. 𝑓 1 and 𝑓 2 are related to the least and most massive star in the cluster, respectively. As we are mainly concerned with the formation of SMSs as the most massive objects in the clusters, we use 𝑓 2 as the last free parameter by default.  \nMotivated by possible constraints on the minimum mass of metalfree stars by stellar archaeology (e.g., Salvadori et al. 2007; Frebel & Norris 2015; Hartwig et al. 2015; Komiya et al. 2016; Ishiyama et al. 2016; Magg et al. 2018, 2019; Dutta et al. 2020; Rossi et al. 2021), here, we consider an alternative parametrisation of multiplicity by adopting 𝑚 min as the third free parameter (in replacement of 𝑓 2 ) with 𝑓 1 = 𝑚 min / 𝑀 given the total mass 𝑀 . Once 𝑓 1 is known according to 𝑚 min , 𝑓 2 = 𝑓 2 ( 𝑁, 𝑓 1 , 𝛼 ) is derived using Eq. 21. Since the solution of 𝑓 2 may not exist for the input (target) value of 𝑁 , we constrain 𝑁 with an upper limit 𝑁 max = 𝑀 / 𝑚 min and a lower limit 𝑁 min corresponding to 𝑓 2 ( 𝑁 min , 𝑓 1 , 𝛼 ) = 1. Below, we briefly explore the effects of multiplicity with the new parametrisation under the default setup (constant 𝑁 , concentrated feedback, 𝛾 eff = 1 . 09, 𝑅 pe ∝ 𝑀 1 . 25 , / 𝑀 ★, crit = 0 . 04 M ⊙ yr -1 , and 𝜂 = 0 . 5).  \nFig. B1 shows the evolution of ˆ 𝑀 ★ with / 𝑀 in for 𝑁 = 3, 10, 30, and 100 given 𝛼 = 1 and 𝑚 min = 1 M ⊙ . We adopt 𝑚 min = 1 M ⊙ as the fiducial choice based on observations (Salvadori et al. 2007; Hartwig et al. 2015; Magg et al. 2019; Rossi et al. 2021). The trends are similar to the case where 𝑓 2 is fixed (Fig. 4), as ˆ 𝑀 ★ becomes higher in the feedback-regulated regime ( / 𝑀 in ≲ 0 . 1 M ⊙ yr -1 ) when 𝑁 increases, showing a slower transition to the lifetime/GRI-regulated regime. However, the effects of varying 𝑁 are stronger with the new parametrisation for 𝑁 ≳ 10 because 𝑓 2 rapidly decreases with 𝑁 given a fixed 𝑚 min , which reduces the ionization flux at / 𝑀 in ≲ 0 . 1 M ⊙ yr -1 and the mass fraction of stars in the bloating phase at / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 . The final mass is even regulated by the gas supply for 𝑁 ≳ 10 at the low / 𝑀 in end. Interestingly, the results for  \nFigure B1. Final total stellar mass as a function of initial gas inflow rate for 𝑁 = 3 (thin solid), 10 (dashed), 30 (dash-dotted), and 100 (dotted) with 𝑚 min = 1 M ⊙ and 𝛼 = 1. The results for the single-star model ( 𝑁 = 1) is shown with the thick solid curve for comparison. The fitting formula of final stellar mass (Eq. 20) based on 2D (single-star) simulations from Hirano et al. (2014, see their fig. 14) is shown with the long-dashed line, and the dots show the underlying data. The results of the 3D AMR simulations in Sugimura et al. (2023) are denoted by the squares, while those of the 3D spherical simulations in Toyouchi et al. (2023) are shown with the diamonds.  \n<!-- image -->  \nFigure B2. Same as Fig. B1 but for 𝑚 min = 1 (thin solid), 3 (dashed), and 10 M ⊙ (dash-dotted), with 𝑁 = 30 and 𝛼 = 1. We also show the results for the default parametrisation with a fixed 𝑓 2 ∼ 0 . 26 as the dotted curve, which agree well with the results for the new parametrisation at / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 where ˆ 𝑀 ★ ≳ 2000 M ⊙ and 𝑓 1 is close to zero.  \n<!-- image -->  \n𝑁 = 3 and 10 are identical at / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 because 𝑁 min ≃ 10 in this regime ( 𝑀 ≳ 2000 M ⊙ ).  \nThe dependence on 𝑚 min is illustrated in in Fig. B2 for 𝑚 min = 1, 3, and 10 M ⊙ given 𝑁 = 30 and 𝛼 = 1. Here the evolution of ˆ 𝑀 ★ with / 𝑀 in has two stages. When / 𝑀 in is low, 𝑁 max is below the input target value of 𝑁 and increases with / 𝑀 in , which explains the rapid increase of ˆ 𝑀 ★ with / 𝑀 in . For 𝑚 min = 10 M ⊙ , ˆ 𝑀 ★ = 𝑚 min is achieved at / 𝑀 in ≲ 10 -4 M ⊙ where 𝑁 max ∼ 1. Once 𝑁 max excesses the target value of 𝑁 with higher / 𝑀 in , the increase of ˆ 𝑀 ★ first slows down (within the no-bloating feedback-regulated regime) before accelerating again when the bloating effect kicks in. In this stage, the dependence on 𝑚 min is minor since 𝑓 1 is close to 0 for all the cases considered here, whose results are well reproduced by the  \n7  \nFigure B3. Same as Fig. B1 but for 𝛼 = 2 . 5 (thin solid), 1.5 (dashed), 1 (dash-dotted), and 0.5 (dotted) with 𝑁 = 3 and 𝑚 min = 1 M ⊙ .  \n<!-- image -->  \nmodel with a fixed 𝑓 2 ≃ 0 . 26 from the default parametrisation, given lim 𝑓 1 → 0 [ 𝑓 2 ( 𝑁 = 30 , 𝑓 1 , 𝛼 = 1 )] ≃ 0 . 26.  \nFinally, we consider the dependence on 𝛼 , as shown in Fig. B3 with 𝑁 = 3 and 𝑚 min = 1 M ⊙ for 𝛼 ∼ 0 . 5 -2 . 5. At / 𝑀 in ≲ 0 . 1 M ⊙ yr -1 , ˆ 𝑀 ★ is higher when the mass distribution is more bottom-heavy, showing a consistent but stronger trend compared with the default case with fixed 𝑓 2 (Fig. 6), which reflects the weaker ionizing power of smaller stars. The trend will be even stronger if 𝑁 is larger. However, for / 𝑀 in ≳ 0 . 1 M ⊙ yr -1 , increasing 𝛼 with a fixed 𝑚 min has a similar effect as increasing 𝑁 while fixing 𝑓 2 (and 𝛼 , see Fig. 4) for 𝛼 ∼ 0 . 5 -2. The reason is that the input value 𝑁 = 3 is no longer effective here as 𝑁 min > 3 and 𝑁 min increases with 𝛼 . In the extreme case of 𝛼 = 2 . 5 strongly dominated by low-mass stars, the final mass is always regulated by photo-evaporation feedback for / 𝑀 in ≤ 10 M ⊙ yr -1 , which is weaker when 𝛼 is larger, so ˆ 𝑀 ★ is higher than that for 𝛼 = 1 . 5, holding the trend at / 𝑀 in ≲ 0 . 1 M ⊙ yr -1 .  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}