Zero Metallicity with Zero CPU Hours: Masses of the First Stars on the Laptop

We develop an analytic model for the mass of the first stars forming in the centers of primordial gas clouds as a function of host halo mass, redshift, and degree of rotation. The model is based on the estimation of key timescales determining the following three processes: the collapse of the gas cloud, the accretion onto the protostellar core, and the radiative feedback of the protostellar core. The final stellar mass is determined by the total mass accreted until the radiative feedback halts the accretion. The analytic estimation, motivated by the result of the full numerical simulations, leads to algebraic expressions allowing an extremely fast execution. Despite its simplicity, the model reproduces the stellar mass scale and its parameter dependencies observed in state-of-the-art cosmological zoom-in simulations. This work clarifies the basic physical principles undergirding such numerical treatments and provides a path to efficiently calibrating numerical predictions against eventual observations of the first stars.


INTRODUCTION
The primordial universe contained only trace metals from big bang nucleosynthesis (BBN; Steigman 2007).The first (Population III, hereafter Pop.III) stars must have formed from this pristine gas, which was cooled principally by atomic and molecular hydrogen.Still, no one has observed a Pop.III star, so our knowledge of this process comes entirely from simulations and analytic estimates (see Bromm & Larson 2004;Bromm 2013;Haemmerlé et al. 2020;Klessen & Glover 2023, for reviews).
Due to the lack of observational comparison and complex physics involved, the Pop.III star formation process is still an uncertain and active area of research.In the hierarchical structure formation scheme, these stars first formed in the universe when halos massive enough to cool and collapse their gas began to virialize.This process began at redshift z ∼ 30 and in halos of masses ∼ 10 5 -10 6 M ⊙ (Tegmark et al. 1997;Bromm & Larson 2004).The Pop. III stars ended the dark ages and began the reionization and metal enrichment of the universe, and the turn-on of Pop.III stars marks the epoch when, for the first time since BBN, nuclear physics becomes relevant in the universe.
In the standard ΛCDM cosmology, the characteristic abundance of the mini-halos which could host Pop.III stars is ∼ 100 cMpc −3 at redshift 20 − 30 (Yoshida et al. 2006).Here, cMpc is the comoving megaparsec.Hence, to develop a statistical sample of tens to hundreds of primordial star-forming clouds, simulations must consider a volume of ∼ cMpc 3 , while the characteristic scale of the proto-stellar accretion disk is ∼ 100AU.This corresponds to a dynamic range of some ten orders of magnitude, which requires zoom-in simulations.The zoom-in simulations beginning from cosmological initial conditions (Greif et al. 2011;Hirano et al. 2015;Stacy et al. 2016;Susa et al. 2014) reveal that Pop.III stars form in small star clusters, with relatively more massive stars (between a few tens and a few thousand solar mass) forming in the center of the collapsing cloud and several lower mass companions which originate from fragmentation in the star-forming disk, as illustrated in Figure 1 of Liu et al. (2020).
Studying the primordial star-cluster systems' formation and evolution demands three-dimensional radiation hydrodynamics simulations in a dense, optically thick environment.Due to computational limits, simulations often do not evolve these primordial clusters to their full maturity, when accretion is shut off by (proto-)stellar feedback, and can only report on the trends observed during the simulation time span.Moreover, since these simulations typically select the first star-forming cloud to form in the cosmological volume as the zoom-in region, they are subject to sampling bias and variance (Stacy & Bromm 2013).
On the other hand, tracking the evolution of a single proto-star in a collapsing cloud including accretion, feedback, and radiative transfer is numerically tractable (Hosokawa et al. 2011(Hosokawa et al. , 2012;;Hirano et al. 2014;Hosokawa et al. 2016;Sharda et al. 2020Sharda et al. , 2021;;Latif et al. 2022).Hirano et al. (2014) attempted to characterize the population of primordial stars by beginning with a cosmological volume of (2 Mpc/h) 3 , resolving ∼ 100 individual star-forming clouds.After its formation, the density of each cloud was azimuthally averaged, then used to initialize 2D radiation hydrodynamics simulations, sidestepping the most computationally demanding part of the problem.Note that the azimuthal averaging smears out any possible small-scale fragmentation and does not resolve the formation and evolution of stellar multiples.However, it allows the central, most massive star in each halo to be evolved until accretion is shut off by feedback.Hirano et al. (2014) found that in their simulations the mass of the central star in each gravitationally unstable cloud is tightly correlated with the collapse timescale of the star-forming cloud.They also found that a long collapse timescale (associated with a low halo mass) allows for the production of HD molecules, permitting cooling to lower temperatures and promoting the formation of low-mass stars.The similar study of Susa et al. (2014) did not include HD chemistry, and found no such correlations.
Here, we expand upon Hirano et al. (2014) by deriving the relationships revealed by the simulations.We show that a simple analytic model based on algebraic timescale arguments can capture the formation of the most massive Pop.III star in the center of a collapsing gas cloud.In particular, this model reproduces the stellar-mass distribution of the sophisticated numerical treatment of Hirano et al. (2014).In the process, we clarify the most important physics underlying Pop.III star formation.
Such an analytic method provides a theoretical "handle" on the Pop.III stellar-mass as a function of the mass and formation redshift of the hosting halo and the rotation parameter.As observations of stellar populations at high redshift become available using, for example, the JWST satellite, such a handle will be necessary to efficiently connect the simulated universe to the physical Universe and to extract cosmological infor-mation from the Pop.III observables.We draw upon previous analytic studies (McKee & Tan 2008;Stahler et al. 1986), generalizing those results to include the effect of the host halo mass and formation redshift on the final stellar mass while wherever possible simplifying the arguments to simple, algebraic relations.
Finally, we contrast the Pop.III star modeling presented here to the modeling of the later generations (Population I and II).The project of deriving the mass of Pop.III stars ab initio is saved from hopelessness only by the simplicity of their environment.These stars form from the first clouds to become gravitationally unstable as halos which exceed the cosmological Jeans mass begin to form.Because the first stars arise from the first baryonic structures to depart from the underlying dark matter distribution, there is reason to believe that their properties can be inferred from the wellunderstood evolution of the dark matter distribution and basic timescale arguments.For this reason, environmental effects can be described in terms of a small number of parameters at the scale of the halo or starforming cloud.That is, deriving the population statistics of the first stars does not require resolving the detailed gas physics, feedback, and radiative transfer in the star-forming clouds over a cosmological volume because the initial conditions can be accurately described using only a few parameters.We defer a detailed study of the enviromental effects such as the Lyman-Werner radiation (e.g.Nebrin et al. 2023) and the relative velocity between the baryons and dark matter (e.g.Nakazato et al. 2022) to a future study.
The paper is organized as follows.In Sec. 2 we present a basic description of the dynamics of the collapsing gas, which informs our determination of the mass of the gravitationally unstable cloud.In Sec. 3 we describe the chemical-thermal network for the primordial gas and compare our results with standard one-zone calculations.In Sec. 4 we develop relations governing the growth of the central star and evaluate the final stellar mass over a range of environmental conditions.We conclude with a discussion of the context of this work and possible future avenues for research.

COLLAPSE DYNAMICS AND CLOUD MASS
The evolution of a primordial gas cloud can be broadly understood through the relations between the following timescales: • The free-fall timescale t f f = 3π 32Gρ , where ρ is the total, dark matter and baryon, density, which is the time for a test particle accelerated by gravity to fall to the center of the cloud.
• The sound crossing timescale t s = r/c s where r is the radius of the cloud and c s = γk B T µm P is the adiabatic sound speed with the adiabatic index γ, temperature T and mean molecular weight µ of the cloud, which is the time for a pressure wave to propagate across the cloud.
• The cooling timescale t C = E th /C, where E th is the thermal energy density of the cloud ([energy][volume] −1 ) and C, the volumetric cooling rate ([energy][time] −1 [volume] −1 ), which is the time for a gas cloud to lose its thermal energy by cooling.Note that the cooling rate depends on the density, composition, and temperature of the cloud.
• The H 2 formation timescale t H2 , which is the time for enough molecular hydrogen to be produced to cool the cloud.Although we use t H2 only in qualitative discussions in this section, we will provide a more quantitative definition of t H2 in Sec.4.1.
• The collapse timescale t col is an e-folding time for the cloud's density, which will be determined from the preceding timescales.
Here's the qualitative sketch of how these timescales interplay to determine the collapsing gas cloud's mass.
The gas in a primordial halo is initially near virial equilibrium, which is equivalent (up to factors of order unity) to the boundary (t f f = t s ) of the Jeans stability condition t f f > t s .This stability condition can also be translated to the Jeans mass, with k B Boltzmann's constant, T the gas temperature, µ the mean molecular weight, m P the proton mass, G Newton's constant, and ρ the gas density.In the second line, we have taken γ = 5/3 for a monatomic ideal gas.Note that different derivations of the Jeans criterion yield different values for the order-unity prefactor, here 1.44.Setting M J approximately equal to M H , the halo mass, determines the virial temperature T of the halo at its formation ρ ≃ 178ρ m (z) , where ρm(z) is the background matter (dark matter and baryon) density at redshift z.
Collapse from this marginally stable configuration at the virial equilibrium requires cooling, which increases the sound crossing time by reducing the temperature.
The dominant coolant of primordial clouds, which consist of pristine gas and whose virial temperature is less than ∼ 10 4 K, is molecular hydrogen.The cosmological molecular hydrogen fraction x H2,0 ≡ n H2 /n H that sets the cloud's initial molecular hydrogen abundance is insufficient to cool the cloud.Therefore, cooling occurs only after enough molecular hydrogen piles up, and the collapse timescale is ultimately determined by the H 2 production timescale t H2 .
The ability of pressure to inhibit collapse depends not only on the temperature but also on the mass of the gas cloud: pressure support can be overcome either by cooling (reducing M J ) or by growing more massive (M cloud exceeding fixed M J ).This fact is responsible for the "loitering phase" (Bromm et al. 2002), which is a pause in the condensation of the gas at the density where cooling becomes inefficient.At a critical density n crit ∼ 103 cm −3 , collisional de-excitation reduces the efficiency of H 2 cooling.Beyond this density, collapse can continue only once enough gas particles have condensed into a cloud at this "loitering" density to exceed the corresponding Jeans mass.Molecular hydrogen alone can cool the gas to ∼ 200 K at n crit (see, for example, Fig. 1).The Jeans mass at this temperature and n crit is ∼ 1000 M ⊙ , which simulations confirm is the approximate mass scale of the gravitationally unstable clouds, when the cloud is cooled by H 2 alone (Bromm et al. 2002).
However, if appreciable HD is formed, the gas can reach temperatures as low as ∼ 50 K, which suppresses the Jeans mass at the critical density and hence the mass of the star-forming cloud.At temperatures above ∼ 500 K, HD is efficiently converted into H 2 , so chemical fractionation of HD begins at T = 500 K.When collapse occurs rapidly (on the free-fall timescale t f f , for example) the time between T = 500 K and the end of the loitering phase at T = 200 K is too brief to form significant HD.However, if the collapse occurs more slowly, with a long enough H 2 production timescale t H2 , then appreciable HD can form, allowing for further cooling and hence a lower cloud mass.The mass of the collapsing cloud in turn determines the mass scale of the central protostar.Environmental effects including mergers and shocks can influence the efficiency of HD cooling (Magnus et al. 2023;Johnson & Bromm 2006).In this work we consider only the case where HD is formed due to a delay in the collapse according to the intrinsic properties of the halo: its mass and its redshift, which are assumed to capture all environmental effects.
Our analytical estimation is based on a so-called "onezone" calculation, which follows the chemical-thermal evolution of a uniform density parcel of gas evolving in a free-fall timescale.For the typical density profile of primordial gas clouds, where the density profile ρ(r) decreases as a function of radius, the free-fall timescale of the inner part of the cloud is much shorter than the outskirts.The inner part therefore undergoes more timescales after a fixed elapsed time, and one can think of this inner part as the later evolutionary stage of the onezone calculation.We therefore expect that the time evolution of one-zone calculation describes the radial profile of the cloud at a fixed time, and indeed the radial profile of full three-dimensional simulations shows a clear correspondence with the results of one-zone calculations (Yoshida et al. 2006).
The thermal evolution of such a gas parcel subject to radiative cooling and adiabatic heating is described by: where T is the temperature, γ the adiabatic index, k B Boltzmann's constant, and ⃗ n the number densities of the various species and n the total nucleon density (i.e.including helium).The first term in the parentheses describes compressional heating due to adiabatic collapse T V γ−1 =const.One-zone calculations do not self-consistently solve for gravity: the density evolution ṅ must be independently specified, and we use a generic parameterization in terms of the collapse timescale t col as which comports with the definition of t col .It is possible to numerically integrate Eq. 3.1 along with the equations for the chemical abundances.Here, we will simplify the problem further by reducing Eq. 3.1 to an algebraic equation with a clear physical interpretation.
Because one-zone calculations are already inexpensive (less than a few seconds on modern consumer hardware), the computational speedup is unlikely to be important.However, this derivation emphasizes the simple physics which determine the density-temperature relationship of the collapsing gas.We begin by rewriting the temperature evolution in terms of the relevant timescales: the collapse timescale t col and the cooling timescale t C .
As expected, for t col ≪ t C we regain the adiabatic heating T ∝ n (γ−1) , while for t col ≫ t C the gas cools.
Table 1.The minimal reaction network, which includes only the dominant formation pathways for H2 and HD, as well as H2-HD interconversion, along with their sources.The temperature and density evolutions described in Eqs.(3.1)-(3.2) are self regulatory: if t col differs greatly from t C at any point, the system will evolve towards t col ∼ t C and vice versa.Specifically, the solution is an attractor.In more physical language, the collapsing gas will evolve towards the equilibrium between heating and cooling described by Eq. (3.4).
Using the attractor solution in Eq. (3.4), therefore, we can find the phase-space diagram T (n) from the algebraic equation without solving the full differential equation of Eq. (3.3).That is the approach that we take in this paper.To do so, we need to estimate the collapse timescale t col (n) and the chemical abundances for the relevant species: free electrons and the primordial coolants H 2 and HD as functions of temperature and density.
We estimate the collapse timescale by parameterizing Here, f is a factor that is approximated to be independent of temperature and density (Hirano et al. 2014).The approximation is reasonable because f accounts for the factor by which the condensation to the loitering phase is extended, which is a well-defined epoch with a single characteristic timescale.In our analysis, we determine f from the H 2 formation timescale t H2 (see Sec. 4 for a quantitative discussion).
Following the argument in Tegmark et al. (1997), we estimate the chemical abundances by adopting a min-imal reaction network (Tab. 1) and analytically solve for the abundances at a fixed temperature, instead of solving for the full coupled differential equations for the reaction rates.For example, we take the equation for the free electron abundance, x e = n e /n H , where n H is the total number density of hydrogen (ionized and neutral), and approximate the right-hand side as Similarly for H 2 and HD we have: where we have defined the effective HD formation rate In these equations, we have made several assumptions.First, that neutral D and H are not depleted and that x e ≡ x p is reduced only by the recombination of neutral hydrogen.Second, that the only important channel for H 2 production is through H − .Third, that x H2 can be treated as constant in Eq. 3.9, because the H 2 fraction is typically nearing its asymptotic value before HD production becomes significant.Finally, to derive the effective HD production rate, Eq. 3.9, we have assumed that the D + produced by reaction k D,3 is promptly converted either back into D or into HD.
As a result, we obtain the following solutions for abundances given their initial values (with the subscript 0): kD,10/kH,1 . (3.12) Tab. 2 summarizes the initial abundances that we use for the computation.The rate coefficients k X are drawn from the "primordial" network of KROME (Grassi et al. 2014), except for k H3 , which is sourced from Galli & Palla (1998).
Finally, to compute the cooling rate C for given abundances, we use the H 2 cooling rates of Hollenbach & McKee (1979) and the HD cooling rate of Lipovka et al. (2005).More recent calculations of the H 2 cooling rates (i.e., Coppola et al. 2019) have more complicated func-tional forms and differ by less than a factor of two from the rates of Hollenbach & McKee (1979) in the relevant temperature regime, which is acceptable in the context of our approximate model.
We solve Eq. (3.4) for a halo with mass M Halo at redshift z.First, we use the spherical collapse (section 5.1 of Desjacques et al. 2018) value of the virial density of the halo ρ V = 178ρ m (z), and estimate the virial tem- perature by Halo . (3.13) We start the computation by setting the cosmological abundances (shown in Tab. 2) as initial values, at ρ i = ρ V /3 (that is, somewhat before the virialization) to allow for any changes in the chemical compositions prior to virialization.Because the cooling is inefficient (t C > t col ) under these initial conditions, we initially assume adiabatic heating T ∝ n 2/3 , and correspondingly initialize the temperature at T i = T V (1/3) 2/3 .The density and temperature then evolve along the adiabatic track until the density n reaches the attractor track: Computing abundances using Eqs.(3.10)-(3.12)requires a temperature, which during this heating phase is supplied as the geometric mean of the initial and current temperature, T = √ T T i .From this intercept with the t col = (γ − 1)t C curve onward, the trajectory is determined by solving t col = (γ − 1)t C for the temperature.In the case where f = 1 (that is, t col = t f f ), collapse occurs too quickly for HD to contribute to the cooling, and only the first two reactions in Tab. 1 are necessary.Moreover, the H 2 production rate depends weakly enough on temperature that evaluating the temperature at each density as T = √ 200T int , where ∼ 200 K is the molecular cooling limit and T int is the temperature when for the first time t col = (γ − 1)t C , can reproduce the one-zone results calculated using KROME (Grassi et al. 2014), as shown in Fig. 1.
However, for f = t col /t f f ≳ 3, HD production becomes important.The relevant reaction rates, for example, k D,10 , are strongly temperature dependent.To capture the temperature-dependent interaction rates, we discretize the attractor curve to 20 pieces up to a density of n H = 10 8 cm −3 .That is, we solve Eq. (3.4) on 20 logarithmically spaced steps of n H beginning at T int and ⃗ n int .In each subsequent step n H , to solve for the temperature T n , the abundances with the subscript 0 in Eqs.(3.10), (3.11) and (3.12) are supplied by the previous step, and the reaction rates are evaluated at a temperature T = T nprev T n given the temperature T nprev from the previous step n prev .Fig. 2 demonstrates an excellent match between the resulting phase-space diagram over a range of values of f ∈ [1, 10].Here, the dashed lines are the result from the algebraic approx-imation in this paper, and solid lines are the solution of the full ODE using the KROME (Grassi et al. 2014) package.Note that although the results disagree somewhat for f = 2 due to our omission of the subdominant HD dissociation reaction HD + H → H 2 + D, our model will require only the temperature at the critical density n crit ∼ 10 3 cm −3 , (Sec.4.1) at which point the curves still agree very closely.With Eq. (3.4) we have transposed the problem of solving a system of coupled ordinary differential equations (already a great simplification of the full, three dimensional problem) to a root finding problem involving inexpensive function calls at a small number of grid points.

STELLAR MASS
We are now in a position to estimate the final stellar mass.Here for simplicity, we assume that only a single protostar forms at the center and grows by accretion via a disk until stellar radiation evaporates the cloud.We will define the collapse timescale t col (that is, fixing f ) to estimate the mass of the collapsing gas cloud.We then determine the final stellar mass in two steps.First, we compute the mass of the initial stellar core that forms before Kelvin-Helmholtz contraction dominates over accretion.We then compute the mass accreted onto this stellar core: the accretion rate is given by dividing the cloud mass with the viscous timescale, and the accretion shut-off time is defined as when the star-forming cloud is completely ionized by the radiation from the protostar.Finally, we compare our result with the fitting formula for the Pop.III star mass derived from 2D radiative hydrodynamic simulations in Hirano et al. (2014), which also assume a single protostar in each star-forming disk.
Actually, recent hydrodynamic simulations of primordial star formation have converged on the picture that Pop.III star-forming disks typically undergo fragmentation and produce multiple protostars (Turk et al. 2009;Stacy et al. 2010;Clark et al. 2011;Greif et al. 2011;Susa 2019;Liu et al. 2020;Haemmerlé et al. 2020;Klessen & Glover 2023) unless stabilized by strong magnetic fields (e.g., Sharda et al. 2020Sharda et al. , 2021;;Hirano & Machida 2022).The formation and evolution of protostars during disk fragmentation (e.g., growth by competing accretion, mergers, and ejections by gravitational scatters) are still poorly understood due to the limited resolution and/or time coverage of current simulations.It is likely that multiple protostars will survive and grow to comparable masses by the moment of cloud evaporation (Sugimura et al. 2020(Sugimura et al. , 2023;;Park et al. 2023a,b).In this case, the strength of ionization feedback will be overestimated by the assumption that the star-forming disk only feeds one central protostar, since the efficiency of producing ionizing photons increases with stellar mass.Therefore, the final stellar mass estimated by our idealized model should be regarded as a lower limit.We defer a more general analysis that includes the situation of multiple protostars to future work.

Collapse Timescale and Mass of the Cloud
As discussed in Sec. 2, the collapse timescale is determined by the H 2 production timescale.If the collapse timescale is long, significant HD can form leading to a lower minimum temperature (Ripamonti 2007).It is this minimum temperature that ultimately sets the mass of the Jeans-unstable cloud.Relatively massive halos form H 2 efficiently within their initial (virial) free-fall timescale.On the other hand, low-mass halos (M halo ∼ 10 5 M ⊙ at z ∼ 20, for example) are initially too cold to rapidly form H 2 and cool further.
To be specific, we define the collapse timescale as where t H2 is the H 2 production timescale that we defined in Sec. 2. The timescales t H2 and t f f are evaluated at the virial density and temperature to determine the value of f , which is assumed to be constant during collapse (Hirano et al. 2014).Quantitatively, t H2 is the timescale required to produce enough H 2 to satisfy t C = t f f at the virial density and temperature.That is, the critical H 2 abundance for cooling x H2,C is defined by where the subscript V stands for quantities in the initial, virial equilibrium and T H2 ≈ 200 K is the minimum temperature achievable by H 2 cooling at the loitering phase.The second term on the left-hand side signifies that we are interested in the time until arrival at the loitering phase at T = T H2 .Note that formally the left hand side becomes negative for T V < T H2 .However, Pop.III star formation in such low mass, low redshift halos is likely very uncommon, and we do not consider this case.Then, with the H 2 production rate we have Once f is determined, we solve the chemical-thermal network as described in Sec. 3, with t col,0 = f t f f,V .In Fig. 2, we find that the HD production indeed lowers the minimum achievable temperature.
We then extract the mass of the cloud from the n-T trajectory as M c = M J (n crit , T (n crit )), where n crit = 1.75 × 10 3 cm −3 characterizes the loitering phase.In a two level system, the critical density is the density at which the H 2 collisional de-excitation rate equals the spontaneous radiative decay rate, which remains a useful heuristic for the multi-level system.For both large and small values of f (corresponding to efficient and inefficient HD formation, respectively) the critical density also approximately corresponds to the minimum temperature.However, for the transitional values of f ∼ 3 the minimum temperature occurs at higher densities as HD continues to form later in the collapse.Even for these cases, we persist in extracting Jeans mass from the temperature at n crit because, in order to affect the mass of the cloud, the HD must be formed before H 2 cooling becomes inefficient at n crit .
Note that while Hirano et al. (2014) have discussed rotation as an important determinant of t col , their results reveal little correlation between the cloud mass and the rotational parameter.This is because, typically, rotation becomes an important stabilizing force only after gravitational instability sets in.An unusually high angular momentum is required for rotational support to develop before the loitering phase.Instead, we argue that rotation determines the infall rate onto the protostar once the cloud mass is fixed.

Accretion Rate
The initial infall rate, as assumed in the KROME calculation [Eq.(3.2)], is Ṁ⋆ = M c /t col . (4.6) However, the accretion rate in later stages and on small scales is typically limited by angular momentum transport.To model this, we adopt the α-disk parameterization (Shakura & Sunyaev 1973) for the viscosity ν = αhc s with h being the scale height.For the thin disk where R is disk radius, v c = ΩR (with Ω the angular velocity) is the circular velocity, the viscous timescale becomes We now make the assumption that the specific angular momentum, J ∝ R 2 Ω is conserved, which is reasonable because in a hot, thick Pop III star-forming disk in the absence of strong (regular) magnetic fields, there are no processes (e.g., resonant torques and outflows) that can efficiently transport angular momentum away.We can thus evaluate the viscous timescale of the disk at any time using the value of R 2 Ω at the critical density n crit as with the subscript "crit" indicating the value of the quantity evaluated at the critical density.Here, we have assumed that at n crit the radius of the cloud is of order the Jeans length, R crit ≈ c s,crit t f f,crit , and introduced the spin parameter as the ratio of rotational to gravitational energy in the cloud at n crit to characterize the degree of rotation.Initially, t col > t ν .As the density increases, t col decreases as 1/ √ n while t ν varies only due to the factorof-few change in sound speed as the gas cools and then heats.Hence, eventually t ν > t col .Therefore, the accretion rate onto the protostar is ultimately set by the viscous timescale.We estimate the accretion rate as Ṁ⋆ = M c /t ν . (4.11) In calculating t ν we assume T = 1000 K, typical of the molecular disk.As the fiducial case we adopt β = 0.3, which is a typical value in Hirano et al. (2014).Here, we take α = 1.The characteristic values of α in Hirano et al. (2014) are a factor of a few lower, but we find that applying α = 1 in Eq. (4.11) better matches the accretion rates in that work (see App. A).We assume Ṁ⋆ is constant, which is equivalent to replacing Ṁ⋆ (t) with its average value.Note that Liu et al. (2020) gives a semi-analytic universal solution for Pop.III growth M ∝ t 4−3γ ef f ≈ t 0.7 (with γ ef f = 1.09 as the effective polytropic index of gas in primordial star-forming disks (Omukai & Nishi 1998)) which is not far from the constant growth rate.

Mass, Radius, and Luminosity
Once the proto-stellar core is formed, the radiation field from the core competes against the accretion flow to determine the evolution during the protostar phase.
For Pop. III stars, Stahler et al. (1986) have studied the evolution of protostars while accretion dominates the dynamics.We adopt their analytical estimate for the protostar core radius: , (4.12) which is surrounded by optically thick radiative precursor with a photospheric radius of R ph = 1.4R ⋆ .Here, the star symbol denotes protostellar quantities.Also, when the opacity is dominated by electron scattering, the luminosity of the hydrostatic equilibrium object is proportional to M 3 ⋆ , and Hosokawa et al. (2012) find the following approximate relationship: Two timescales are relevant in determining the contracting mass scale of the protostellar core: the accretion timescale and the Kelvin-Helmholtz timescale, Initially, the accretion timescale is short compared to the Kelvin-Helmholtz timescale, and the protostar expands according to Eq. (4.12).Eventually, however, the mass and luminosity growth of the protostar cause Kelvin-Helmholtz contraction to dominate over accretion: t KH < t acc , which happens at M eq (Hosokawa et al. 2012): at which point the protostellar cores with luminosity less than the Eddington luminosity, L(M eq ) < L Edd , begin to contract.The following energy balance equation can model the contraction: where W is the gravitational binding energy of the protostar and we assume that each step of contraction plus accretion maintains a new virial equilibrium by radiating the energy difference.Here, n is the polytropic index P ∝ ρ (n+1)/n and we choose n = 3, consistent with the Eddington-beta model (Eddington 1926), which is reasonably accurate for massive stars.Again, assuming the hydrostatic equilibrium protostellar core with constant opacity, dominated by electron scattering, L ⋆ ∝ M 3 ⋆ , we solve the energy balance equation (Eq.(4.17)) to find the radius-mass relationship: Eventually, the luminosity of the collapsing proto-star reaches the Eddington limit (Hosokawa et al. 2012), from which point radiation pressure prevents the contraction of the envelope (Hosokawa et al. 2012;Hirano et al. 2014).Then, the protostellar radius must either oscillate or grow.Substituting the Eddington luminosity into Eq.(4.17) gives a nearly constant radius.However, Hosokawa et al. (2012) find that once the protostars reach the Eddington luminosity, the gravothermal evolution contracts the core while expanding the envelope.The result of these more complicated dynamics is to grow the characteristic radius as R ⋆ ∝ M 0.5 ⋆ , which we adopt here.Note that Hosokawa et al. (2012) find that for very high accretion rates ≳ 0.1 M ⊙ /yr the proto-star begins to expand even before L Edd , an effect which we neglect, both because the underlying physics are complex and because this mechanism principally operates for rotation rates lower than our fiducial β = 0.3.This omission will lead to an underestimate of the masses of the most massive, slowly rotating stars.
We show the evolution of radius and effective temperature of the protostars in Fig. 3 for three different accretion rates.

Ionizing Feedback and Stellar Mass
Equipped with the effective temperature and radius that we estimated in the previous section (see Fig. 3 for the result), we can also estimate the ionizing photon flux, S EUV , by integrating the black body spectrum.
The interaction between the UV flux and the hydrodynamical gas launches a shock.We assume that the shock is launched into the cloud at t eq (i.e. the time when the mass reaches M eq ) which is the characteristic epoch where the protostar begins to emit significant UV radiation.Once the shock is launched, it homogenizes the downstream medium and reduces its density.Initially, the UV radiation is trapped behind the shock front by recombinations.Eventually, the density behind the shock (that is, towards the center of the cloud) falls low enough that recombination is no longer efficient.Then, the ionization front escapes the shock front, rapidly ionizing the cloud and shutting down accretion.This cloud breakout time t B is estimated in Alvarez et al. (2006) using the similarity solution of Shu et al. (2002) for a singular isothermal sphere (SIS) with a density profile of ρ ∝ r −2 as the initial condition.There, t B is calculated by equating the ionizing photon flux at t B after the shock onset to the recombination rate behind the shock front: where α B is the case-B recombination rate coefficient for hydrogen at the temperature characteristic of the photoheated gas ∼ 10 4 K, and r sh = x s c s t (with c s the sound speed in the shocked gas) is the radius of the shock front, and the normalization of the density profile n is proportional to the temperature T SIS of the initial SIS, which we take to be the temperature at the loitering phase, i.e.
T SIS = T (n crit ).When the shocked gas is much hotter than the surrounding isothermal sphere (which holds in all cases we encounter), we have x s ≈ 2.56.Alvarez et al. (2006) provide the following approximate relationship for the breakout time, which we find matches the exact result to within ∼ 10%: Note that helium at the cosmological mass fraction Y = 0.24 enters the calculation of the sound speed via the mean molecular weight.We assume the first ionization of helium is coupled with that of hydrogen, and correspondingly have multiplied the prefactor in Eq. 4.21 by x −2 H = 1.16 compared to Alvarez et al. (2006) (who neglected the consumption of ionizing photons by helium) given the primordial number fraction of hydrogen nuclei x H = 0.927.We finally solve Eq. (4.21) for t B , and then compute the final stellar mass as (4.22) The solution in Shu et al. (2002) ignores the gravity from the central star and the fact that the density profile in the central region is shallower than ρ ∝ r −2 in the presence of the star-forming disk.These effects tend to slow down the propagation of the ionization front in the early stage so that the ionized region can be trapped in the center for an extended period t trap ≳ 20 kyr before breaking into the cloud (Jaura et al. 2022).Since this initial trapped phase is not considered in our model, the final stellar mass can be underestimated1 .
Having relied on the timescales and analytical estimate, our estimation of the final stellar mass only takes a few seconds.On the other hand, with a series of 2D radiative hydrodynamic simulations of protostar accretion from initial conditions of primordial star-forming clouds produced by cosmological hydrodynamic (zoom- in) simulations, Hirano et al. (2014) developed a sample of ∼ 100 Population III stars.Fig. 4 shows that our estimation very closely resembles the result of the state of the art simulations!The figure shows the predictions of our model for β = 0.3 as a function of halo mass and redshift.We have also plotted the stellar sample of Hirano et al. (2014) (dots) and the line corresponding to a halo abundance of one in the simulation volume of Hirano et al. (2014) in the Press-Schechter Formalism (Press & Schechter 1974;Sheth et al. 2001).Halos become common enough to typically appear in the simulation volume of Hirano et al. (2014) only to the right of this line.However, Hirano et al. (2014) employ additional selection criteria to ensure that each star-forming cloud is pristine.In the bottom right corner of the figure, T V < T H2 , and our calculation of f does not apply.As suggested above, primordial star formation in this parameter space is not important.
Our model succeeds quite well in predicting the transition between H 2 cooled clouds (final stellar mass of hundreds of solar masses) and HD cooled clouds (final stellar masses of order tens of solar masses).Over the range where Hirano et al. (2014) have reasonable statistics (∼ 30 − 300 M ⊙ ), our model agrees with the simulation results to within a factor of a few.However, our model fails to produce the most massive stars observed in Hirano et al. (2014), a fact which cannot be entirely explained by our choice of fixed β = 0.3 in Fig. 4 (see Fig. 5).There are several physically reasonable explanations which may contribute to this discrepancy, which are explored in more detail in App. A. First, the rapidly accreting stars of Hirano et al. (2014) can grow by around 100 M ⊙ after HII breakout due to the geometry of the accretion disk which is not captured by our spherically-symmetric breakout model.Relatedly, our assumption of a thin disk in calculating t ν breaks down for the slowly rotating clouds which produce the most massive stars, leading to an underestimate of the accretion rate.Finally, as mentioned above our model fails to account for the finding of Hosokawa et al. (2012) that for very high accretion rates the proto-star never undergoes a period of contraction.The low surface temperature associated with this large stellar radius is necessary to attain stellar masses ∼ 1000 M ⊙ .These facts can largely explain the systematically lower masses predicted by our model at high accretion rates.
We also show in Fig. 5 the effect of varying β on a low mass, HD cooled cloud and a high mass, H 2 cooled cloud at redshift 25.In our model, the largest stellar mass attainable with a minimal realistic rotation parameter β ≈ 0.05 is M ∼ 200 M ⊙ , with an accretion rate of ∼ 0.01 M ⊙ /yr.Rotation has a larger effect at higher accretion rates (lower β), because the efficiency of the UV feedback has a strong dependence on the temperature at which the surface of the Eddington radiating stars is "frozen" (Fig. 3).

CONCLUSION & DISCUSSION
We have developed a simplistic model of the formation of Pop.III stars in the center of collapsing primordial gas clouds.The model consists of two parts.The first determines the chemical-thermal evolution of the collapsing gas cloud using the dynamical-thermal equilibrium relation t col = (γ − 1)t C .The second relates the mass and spin parameter of the collapsing cloud to the final mass of the star.For a typical value of the spin parameter, this model agrees with the masses predicted by sophisticated simulations (Hirano et al. 2014) to within a factor of a few for stellar masses between ∼ 30 − 300 M ⊙ .However, the model struggles to produce the most massive stars seen in simulations, for which the aspherical accretion geometry strongly affects the final mass (see App. A).Moreover, simulations such as those of Hirano et al. (2014) which do not resolve the trapped phase of HII region (Jaura et al. 2022) and produce only one star per cloud may underestimate the star formation efficiency.Since our model also assumes one star per cloud and ignores the trapped phase, the predicted final stellar mass is best understood as a lower bound.
The model makes liberal use of the ≈ symbol.Although it contains no explicit free parameters, there are implicit choices involving various order unity factors which can bring the model into better or worse agreement with simulations.These include the numerical prefactors in the sound crossing and free-fall timescales, the value of the viscosity parameter α, the characteristic temperature of the disk used to evaluate the viscous timescale t ν , and the characteristic temperature of the singular isothermal sphere used in Eq. (4.20).The overall trends are robust to these choices.
We comment briefly that based only on the chemicalthermal evolution of the gas (Sec.3) two alternative estimates of the stellar mass are already possible.First, the Jeans mass at the loitering phase can be multiplied by some global efficiency factor M ⋆ = ϵM c , where to match the characteristic ∼ 100 M ⊙ Pop.III stellar mass ϵ ∼ 0.1.However, in our model the efficiency factor depends on the cloud mass.We find that ϵ is a factor of a few larger for smaller, HD cooled clouds than for larger, H 2 cooled clouds.This is because while smaller clouds lead to slower accretion, the slower accretion in turn leads to a longer period of growth before the breakout of the ionization front.Alternatively, one could estimate M ⋆ = N M l , where M l ≈ 1.4M ⊙ µ −9/4 k B T mpc 2 1/4 is the opacity limited minimum fragment mass (Rees 1976), evaluated at the minimum temperature of the gas (Shandera et al. 2018;Singh et al. 2021), and N ∼ 10 4 (as the effective number of fragments that merge to make the central star) to produce reasonable stellar masses.In this case, if N is constant, the weak dependence of M l on T means that the temperature difference between H 2 and HD cooled clouds amounts to less than a factor of two difference in M l and hence in M ⋆ .Additionally, neither of these simpler estimates include the role of angular momentum.
Our model correctly determines the order-ofmagnitude mass scale of the first stars, and the dependence of that mass scale on redshift, halo mass, and the rotation parameter.These dependencies, which were initially emergent properties in radiation-hydrodynamic simulations, are here distilled to simple timescale arguments.If Pop.III stars are observed, this kind of understanding will provide a powerful lever to calibrate simulations against observations.In the meanwhile, this model can be inserted in simulations of cosmological volumes at a low computational cost, allowing a novel treatment of the metal enrichment of the universe and subsequent reionization history.
Another possible application is to the study of dissipative dark matter, which can itself cool to eventually form compact objects (Shandera et al. 2018;Hippert et al. 2022;Gurian et al. 2022;Ryan & Radice 2022).These objects could have masses and compactnesses impossible under ordinary stellar astrophysics, leading to distinctive gravitational wave signatures.However, given the large model space there is a need for simple, inexpensive, and reasonably accurate models of the compact object formation process.The core ideas present in this model could be generalized to other dissipative physics, providing just such a tool.
Finally, although this work only considers the formation of the massive central star in each cloud, Liu et al. (2020) found that the formation of Pop.III star clusters by disk fragmentation can be described by simple scaling laws that capture the key trends in 3D hydrodynamic simulations of primordial star-forming clouds.Our model can be generalized and improved to consider Pop.III star clusters and self-shielding in aspherical accretion flows using the results of 3D radiative hydrodynamic simulations that follow the growth and feedback of multiple protostars (e.g.Sugimura et al. 2020Sugimura et al. , 2023;;Park et al. 2023b,a), which is an intriguing direction for future research.A1), compared to our fiducial result (Fig. 4).Also shown are contours of the fit (Eq.A3) of Hirano et al. (2014) to their sample.The range of stellar masses is increased by a factor of a few.A2) with the accretion rate supplied by Eq. (4.11).Also shown are contours of the fit (Eq.A3) of Hirano et al. (2014).For most of the parameter space, the stellar mass is enhanced compared to our fiducial model.This enhancement is strongest for the most massive halos (and hence the most massive stars).

Figure 1 .
Figure 1.Evolution of the one-zone parcel in phase space for the t col = t f f case.The solid black line shows the result of solving the chemical network using KROME, and the other two lines show two analytic approximations.For the blue dot-dashed line, the reaction rates in Eqs.(3.10-3.12)are always evaluated at a single temperature T = √ 200Tint when Eq. (3.4) is solved at any density.For the green dashed line, we solve Eq. (3.4) at twenty logarithmically spaced steps of density nH ∈ [nint, 10 8 cm −3 ], where to solve for the temperature Tn at a step n, the reaction rates are evaluated at a temperature T = Tn prev Tn given the temperature Tn prev from the previous step nprev.The initial density and temperature are those of a halo with M Halo = 5 × 10 5 M⊙ at z = 15.

Figure 2 .
Figure 2. Evolution of the one-zone parcel in phase space for five f ≡ t col /t f f values.The solid lines show the results of integrating the chemical network using KROME, and the dashed lines show the solution of the algebraic relationship tC = (γ − 1)t col , with 20 discretized density values (see the main text).The algebraic solution matches well with the KROME result.Also shown are lines of constant Jeans mass MJ : The HD formation lowers the minimum temperature, which reduces the Jeans mass at the critical density ncrit ∼ 10 3 cm −3 .The initial density and temperature are those of a halo with M Halo = 5 × 10 5 M⊙ at z = 15.

Figure 4 .
Figure 4.The predicted stellar mass in our model (background color, contours) compared with the primordial stars found in the simulations of Hirano et al. (2014) (colored circles).The gray shaded region indicates halos which will not typically appear in the simulation volume of Hirano et al. (2014) in the Press-Schechter formalism.In the bottom right, TV < TH 2 and we do not calculate a stellar mass.Our model is accurate to within a factor of a few where Hirano et al. (2014) have robust data (between ∼ 30 and ∼ 300 solar mass).

Figure 5 .
Figure 5.The role of the rotation parameter β in the final stellar mass for an H2 cooled (10 6 M⊙, gold) cloud and HD cooled (10 5 M⊙, purple) cloud, both at z = 25, showing the stronger dependence on β at higher cloud masses.

Figure 6 .
Figure6.The accretion rate from Eq. 4.11 (black) and from Eq. A1 (red, dashed) and the sample ofHirano et al. (2014), plotted against the mass of the star-forming cloud.Although the analytically calculated accretion rate is broadly consistent with the data, the power-law index is shallower than the best fit.

Figure 7 .
Figure7.The stellar mass calculated in our model using the accretion rate Eq.(A1), compared to our fiducial result (Fig.4).Also shown are contours of the fit (Eq.A3) ofHirano et al. (2014) to their sample.The range of stellar masses is increased by a factor of a few.

Figure 8 .
Figure 8.The stellar mass calculated from Eq. (A2) with the accretion rate supplied by Eq. (4.11).Also shown are contours of the fit (Eq.A3) ofHirano et al. (2014).For most of the parameter space, the stellar mass is enhanced compared to our fiducial model.This enhancement is strongest for the most massive halos (and hence the most massive stars).

Table 2 .
The initial fractional abundances and their sources.