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

James Gurian  Institute for Gravitation and the Cosmos, The Pennsylvania State
University, University Park, PA 16802, USA  Department of Astronomy and
Astrophysics, The Pennsylvania State University, University Park, PA, 16802,
USA  Perimeter Institute for Theoretical Physics, Waterloo, Ontario, N2L 2Y5,
Canada  Donghui Jeong  Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Department of
Astronomy and Astrophysics, The Pennsylvania State University, University
Park, PA, 16802, USA,  School of Physics, Korea Institute for Advanced Study
(KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul, 02455, Republic of Korea  Boyuan Liu
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK

###### Abstract

We develop an analytic model for the mass of the first stars forming in the
center 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
dependences 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.

Cosmology(343), Population III stars(1285) , Star formation(1569)

##  1 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\sim 30$ and in halos
of masses $\sim 10^{5}$–$10^{6}\,\rm M_{\odot}$ (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 $\Lambda$CDM cosmology, the characteristic abundance of the
mini-halos which could host Pop. III stars is $\sim 100\,\rm 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 $\sim\rm
cMpc^{3}$, while the characteristic scale of the proto-stellar accretion disk
is $\sim 100\rm AU$. 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 (Susa et al., 2014; Hirano et
al., 2015; Stacy et al., 2016; Susa, 2019) 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, 2012; Hirano et al., 2014;
Hosokawa et al., 2016; Sharda et al., 2020, 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\,{\rm Mpc}/h)^{3}$, resolving
$\sim 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 ${\rm
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 ${\rm 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 information 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 well-
understood 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 star-forming 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.

##  2 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_{ff}=\sqrt{\frac{3\pi}{32G\rho}}$, where $\rho$ 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}=\sqrt{\frac{\gamma k_{B}T}{\mu m_{P}}}$ is the adiabatic
sound speed with the adiabatic index $\gamma$, temperature $T$ and mean
molecular weight $\mu$ of the cloud, which is the time for a pressure wave to
propagate across the cloud.

  * •

The cooling timescale $t_{{\cal C}}={\cal E}_{th}/{\cal C}$, where ${\cal
E}_{th}$ is the thermal energy density of the cloud
($\rm[energy][volume]^{-1}$) and ${\cal C}$, the volumetric cooling rate
($\rm[energy][time]^{-1}[volume]^{-1}$), which is the time for a gas cloud to
loose its thermal energy by cooling. Note that the cooling rate depends on the
density, composition, and temperature of the cloud.

  * •

The ${\rm H_{2}}$ formation timescale $t_{{\rm H_{2}}}$, which is the time for
enough molecular hydrogen to be generated. Although we use $t_{{\rm H_{2}}}$
only in qualitative discussions in this section, we will provide a more
quantitative definition of $t_{\rm H_{2}}$ in Sec. 4.1.

  * •

The collapse timescale $t_{\rm 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_{ff}=t_{s}$) of
the Jeans stability condition $t_{ff}>t_{s}$. This stability condition can
also be translated to the Jean’s mass,

$\displaystyle M_{J}$ $\displaystyle=\frac{4\pi\rho c_{s}^{3}}{3t_{ff}^{3}}$
(2.1) $\displaystyle\approx 1.44\left(\frac{k_{B}T}{\mu
m_{P}G}\right)^{3/2}\rho^{-1/2},$ (2.2)

with $k_{B}$ Boltzmann’s constant, $T$ the gas temperature, $\mu$ the mean
molecular weight, $m_{P}$ the proton mass, $G$ Newton’s constant, and $\rho$
the gas density. In the second line, we have taken $\gamma=5/3$ for a
monatomic ideal gas. Note that different derivations of the Jean’s 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 $\rho\simeq 178\bar{\rho}_{\rm
m}(z)$ , where $\bar{\rho}_{\rm 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 $\sim 10^{4}$ K, is
molecular hydrogen. The cosmological molecular hydrogen fraction $x_{{\rm
H_{2}},0}\equiv n_{{\rm H_{2}}}/n_{\rm 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 ${\rm H_{2}}$ production timescale
$t_{\rm H_{2}}$.

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
overcomed either by cooling (reducing $M_{J}$) or by growing more massive
($M_{\rm 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_{\rm crit}\sim 10^{3}\,\rm cm^{-3}$, collisional de-excitation
reduces the efficiency of ${\rm 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 $\sim 200\,\rm K$ at $n_{\rm crit}$ (see,
for example, Fig. 1). The Jean’s mass at this temperature and $n_{\rm crit}$
is $\sim 1000\,M_{\odot}$, which simulations confirm is the approximate mass
scale of the gravitationally unstable clouds, when the cloud is cooled by
${\rm H_{2}}$ alone (Bromm et al., 2002).

However, if appreciable ${\rm HD}$ is formed, the gas can reach temperatures
as low as $\sim 50\,\rm K$, which suppresses the Jeans mass at the critical
density and hence the mass of the star-forming cloud. At temperatures above
$\sim 500\,\rm K$, ${\rm HD}$ is efficiently converted into ${\rm H_{2}}$, so
the ${\rm HD}$ formation begins at $T=500\ \rm K$. When collapse occurs
rapidly (on the free-fall timescale $t_{ff}$, for example) the time between
$T=500\,\rm K$ and the end of the loitering phase at $T=200\,\rm K$ is too
brief to form significant ${\rm HD}$. However, if the collapse occurs more
slowly, with a long enough ${\rm H_{2}}$ production timescale $t_{{\rm
H_{2}}}$, then appreciable ${\rm 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 ${\rm HD}$
cooling (Magnus et al., 2023; Johnson & Bromm, 2006). In this work we consider
only the case where ${\rm 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.

##  3 Gas Chemistry

Our analytical estimation is based on a so-called “one-zone” 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 $\rho(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 time-
scales after a fixed elapsed time, and one can think of this inner part as the
later evolutionary stage of the one-zone 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:

$\frac{dT}{dt}=(\gamma-1)\left(\frac{\dot{n}}{n}T-\frac{{\cal
C}(T,{\vec{n}})}{k_{B}n}\right),$ (3.1)

where $T$ is the temperature, $\gamma$ the adiabatic index, $k_{B}$
Boltzmann’s constant, and $\vec{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 $TV^{\gamma-1}=$const.

One-zone calculations do not self-consistently solve for gravity: the density
evolution $\dot{n}$ must be independently specified, and we use a generic
parameterization in terms of the collapse timescale $t_{\rm col}$ as

$\displaystyle\frac{dn}{dt}$
$\displaystyle=-3n\frac{\dot{r}}{r}\approx\frac{n}{t_{{\rm col}}(n)},$ (3.2)

which comports with the definition of $t_{\rm col}$. Then, the temperature
evolution can be rewritten in terms of the relevant timescales, the collapse
timescale $t_{\rm col}$ and the cooling timescale $t_{\cal C}$.

$\frac{d\log T}{d\log n}=(\gamma-1)\left[1-\frac{t_{{\rm
col}}(n)}{(\gamma-1)t_{{\cal C}}({\vec{n}},T)}\right].$ (3.3)

As expected, for $t_{{\rm col}}\ll t_{{\cal C}}$ we regain the adiabatic
heating $T\propto n^{(\gamma-1)}$, while for $t_{{\rm col}}\gg t_{\cal C}$ the
gas cools.

The temperature and density evolutions described in Eqs. (3.1)–(3.2) are self
regulatory: if $t_{\rm col}$ differs greatly from $t_{\cal C}$ at any point,
the system will evolve towards $t_{{\rm col}}\sim t_{\cal C}$ and vice versa.
Specifically, the solution

$t_{{\rm col}}(n)=(\gamma-1)t_{{\cal C}}({\vec{n}},T)$ (3.4)

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_{\rm col}(n)$
and the chemical abundances for the relevant species: free electrons and the
primordial coolants ${\rm H_{2}}$ and ${\rm HD}$ as functions of temperature
and density.

We estimate the collapse timescale by parameterizing $t_{{\rm
col}}(n)=ft_{ff}(n)$. 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 ${\rm H_{2}}$ formation timescale $t_{\rm H_{2}}$ (see Sec. 4 for a
quantitative discussion).

$k_{\rm H,1}$ | $\rm p+e\rightarrow H$
---|---
$k_{\rm H,3}$ | $\rm H+e\rightarrow H^{-}$
$k_{\rm D,3}$ | $\rm D+H^{+}\rightarrow D^{+}+H$
$k_{\rm D,4}$ | $\rm D^{+}+H\rightarrow D+H^{+}$
$k_{\rm D,8}$ | $\rm D^{+}+H_{2}\rightarrow HD+H^{+}$
$k_{\rm D,{10}}$ | $\rm HD+H^{+}\rightarrow H_{2}+D^{+}$
Table 1: The minimal reaction network, which includes only the dominant formation pathways for ${\rm H_{2}}$ and ${\rm HD}$, as well as ${\rm H_{2}}$-${\rm HD}$ interconversion.  Species | Initial Abundance | Source
---|---|---
$x_{e}$ | $2.5\times 10^{-4}$ | Recfast (Seager et al., 1999)
$x_{{\rm H_{2}}}$ | $7\times 10^{-7}$ | Hirata & Padmanabhan (2006)
$x_{\rm D}$ | $2.5\times 10^{-5}$ | Cooke et al. (2018)
$x_{\rm D^{+}}$ | $6.3\times 10^{-9}$ | $x_{\rm D+}/x_{\rm D}\equiv x_{\rm H+}/x_{\rm H}$
$x_{\rm HD}$ | $1.8\times 10^{-11}$ | $x_{\rm HD}/x_{\rm D}\equiv x_{\rm H_{2}}/x_{\rm H}$
Table 2: The initial fractional abundances and their sources.

Following the argument in Tegmark et al. (1997), we estimate the chemical
abundances by adopting a minimal 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_{\rm H}$,
where $n_{\rm H}$ is the total number density of hydrogen (ionized and
neutral),

$\displaystyle\frac{dx_{e}}{dn_{{\rm H}}}$ $\displaystyle=-k_{{\rm
H},1}x_{e}^{2}n_{\rm H}{t_{{\rm col}}(n_{{\rm H}})}\frac{dt}{dn_{{\rm H}}},$
(3.5)

and approximate the right-hand side as

$\displaystyle\frac{dx_{e}}{dn_{{\rm H}}}$ $\displaystyle=-k_{{\rm
H},1}x_{e}^{2}n_{\rm H}\left(\frac{t_{\rm col}(n_{{\rm H}})}{\sqrt{n_{{\rm
H}}/n_{{\rm H},0}}}\frac{1}{n_{{\rm H}}}\right),$ (3.6)

to find the analytic solution. As a result, we obtain the following solutions
for abundances given their initial values (with the subscript 0):

$\displaystyle x_{{\rm e}}(n_{{\rm H}})$ $\displaystyle=\frac{x_{{\rm
e},0}}{1+2k_{\rm H,1}t_{{\rm col},0}x_{{\rm e},0}(\sqrt{n_{{\rm H}}n_{{\rm
H},0}}-n_{{\rm H},0})}$ (3.7) $\displaystyle x_{\rm H_{2}}(n_{{\rm H}})$
$\displaystyle=x_{{\rm H_{2}},0}$ $\displaystyle\quad+\frac{k_{\rm
H,3}}{k_{\rm H,1}}\log{\left[1+2k_{\rm H,1}t_{{\rm col},0}x_{{\rm
e},0}(\sqrt{n_{{\rm H}}n_{{\rm H},0}}-n_{{\rm H},0})\right]}$ (3.8)
$\displaystyle x_{{\rm HD}}(n_{{\rm H}})$ $\displaystyle=x_{{\rm
HD},0}\left(\frac{x_{\rm e}(n_{{\rm H}})}{x_{{\rm e},0}}\right)^{k_{\rm
D,10}/k_{\rm H,1}}$ $\displaystyle\quad+\frac{k_{\rm{\rm HD},eff}}{k_{\rm
D,10}}x_{\rm D,0}\left(1-\left(\frac{x_{\rm e}(n_{{\rm H}})}{x_{{\rm
e},0}}\right)^{k_{\rm D,10}/k_{\rm H,1}}\right),$ (3.9)

where we have assumed that neutral $\rm D$ and $\rm H$ are not depleted and
that $x_{e}$ is reduced only by the recombination of neutral hydrogen. Here,
we have defined the effective ${\rm HD}$ formation rate

$\displaystyle k_{\rm{\rm HD},eff}$ $\displaystyle=k_{\rm D,3}\frac{k_{\rm
D,8}}{k_{\rm D,8}+k_{\rm D,4}/x_{\rm H_{2}}}.$ (3.10)

In deriving Eq. (3.9), we treat $x_{{\rm H_{2}}}$ as constant because the
${\rm H_{2}}$ fraction is typically nearing its asymptotic value before ${\rm
HD}$ production becomes significant. Tab. 2 summarizes the initial abundances
that we use for the computation. The rate coefficients $k_{X}$ are drawn from
Table C1 of Grassi et al. (2014).

Finally, to compute the cooling rate ${\cal C}$ for given abundances, we use
the ${\rm H_{2}}$ cooling rates of Hollenbach & McKee (1979) and the ${\rm
HD}$ cooling rate of Lipovka et al. (2005). More recent calculations of the
${\rm H_{2}}$ cooling rates differ only marginally from the results of
Hollenbach & McKee (1979), which have a simpler functional form (Ryan et al.,
2022).

We solve Eq. (3.4) for a halo with mass $M_{\rm 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 $\rho_{V}=178\bar{\rho}_{\rm m}(z)$, and
estimate the virial temperature by

$T_{V}=\left(\frac{4\pi}{3}\rho_{V}\right)^{1/3}\frac{G\mu m_{H}}{k_{B}}M_{\rm
Halo}^{2/3}.$ (3.11)

We start the computation by setting the cosmological abundances (shown in Tab.
2) as initial values, at $\rho_{i}=\rho_{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_{{\cal C}}>t_{{\rm
col}}$) under these initial conditions, we initially assume adiabatic heating
$T\propto 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: $t_{{\rm
col}}(n)=(\gamma-1)t_{{\cal C}}({\vec{n}},T)$. Computing abundances using Eqs.
(3.7)–(3.9) requires a temperature, which during this heating phase is
supplied as the geometric mean of the initial and current temperature,
$\bar{T}=\sqrt{TT_{i}}$.

From this intercept with the $t_{{\rm col}}=(\gamma-1)t_{{\cal C}}$ curve
onward, the trajectory is determined by solving $t_{{\rm
col}}=(\gamma-1)t_{{\cal C}}$ for the temperature.

Figure 1: Evolution of the one-zone parcel in phase space for the $t_{{\rm
col}}=t_{ff}$ 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.7-3.9) are always evaluated at a single temperature
$\bar{T}=\sqrt{200T_{\rm int}}$ 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 $n_{{\rm H}}\in[n_{\rm int},10^{8}\ \rm cm^{-3}]$, where to
solve for the temperature $T_{n}$ at a step $n$, the reaction rates are
evaluated at a temperature $\bar{T}=\sqrt{T_{n_{\rm prev}}T_{n}}$ given the
temperature $T_{n_{\rm prev}}$ from the previous step $n_{\rm prev}$. The
initial density and temperature are those of a halo with $M_{\rm Halo}=5\times
10^{5}M_{\odot}$ at $z=15$.

In the case where $f=1$ (that is, $t_{{\rm col}}=t_{ff}$), collapse occurs too
quickly for ${\rm HD}$ to contribute to the cooling, and only the first two
reactions in Tab. 1 are necessary. Moreover, the ${\rm H_{2}}$ production rate
depends weakly enough on temperature that evaluating the temperature at each
density as $\bar{T}=\sqrt{200T_{\rm int}}$, where $\sim 200\,\rm K$ is the
molecular cooling limit and $T_{\rm int}$ is the temperature when for the
first time $t_{{\rm col}}=(\gamma-1)t_{{\cal C}}$, can reproduce the one-zone
results calculated using KROME (Grassi et al., 2014), as shown in Fig. 1.

Figure 2: Evolution of the one-zone parcel in phase space for five $f\equiv
t_{{\rm col}}/t_{ff}$ 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 $t_{{\cal C}}=(\gamma-1)t_{{\rm 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
$M_{J}$: The ${\rm HD}$ formation lowers the minimum temperature, which
reduces the Jeans mass at the critical density $n_{\rm crit}\sim 10^{3}\,{\rm
cm}^{-3}$. The initial density and temperature are those of a halo with
$M_{\rm Halo}=5\times 10^{5}M_{\odot}$ at $z=15$.

However, for $f=t_{{\rm col}}/t_{ff}\gtrsim 3$, ${\rm HD}$ production becomes
important. The relevant reaction rates, for example, $k_{\rm 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_{{\rm H}}=10^{8}\,{\rm cm}^{-3}$. That is, we solve Eq. (3.4) on
20 logarithmically spaced steps of $n_{{\rm H}}$ beginning at $T_{\rm int}$
and $\vec{n}_{\rm int}$. In each subsequent step $n_{{\rm H}}$, to solve for
the temperature $T_{n}$, the abundances with the subscript 0 in Eqs. (3.7),
(3.8) and (3.9) are supplied by the previous step, and the reaction rates are
evaluated at a temperature $\bar{T}=\sqrt{T_{n_{\rm prev}}T_{n}}$ given the
temperature $T_{n_{\rm prev}}$ from the previous step $n_{\rm prev}$. Fig. 2
demonstrates an excellent match between the resulting phase-space diagram over
a range of values of $f\in[1,10]$. Here, the dashed lines are the result from
the algebraic approximation 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
${\rm HD}$ dissociation reaction ${\rm HD}+\rm H\rightarrow{\rm H_{2}}+\rm D$,
our model will require only the temperature at the critical density $n_{\rm
crit}\sim 10^{3}\,\rm 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.

##  4 Stellar Mass

We are now in a position to estimate the final stellar mass. We will define
the collapse timescale $t_{\rm 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 given in Hirano et al. (2014).

###  4.1 Collapse Timescale and Mass of the Cloud

As discussed in Sec. 2, the collapse timescale is determined by the ${\rm
H_{2}}$ production timescale. If the collapse timescale is long, significant
${\rm 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 ${\rm H_{2}}$ efficiently within
their initial (virial) free-fall timescale. On the other hand, low-mass halos
($M_{\rm halo}\sim 10^{5}\,M_{\odot}$ at $z\sim 20$, for example) are
initially too cold to rapidly form ${\rm H_{2}}$ and cool further.

To be specific, we define the collapse timescale as

$t_{{\rm col}}\approx\min{(t_{ff},t_{{\rm H_{2}}})},$ (4.1)

where $t_{\rm H_{2}}$ is the ${\rm H_{2}}$ production timescale that we define
in Sec. 2. Specifically, $t_{\rm H_{2}}$ is the timescale required to produce
enough ${\rm H_{2}}$ to satisfy $t_{{\cal C}}=t_{ff}$ at the virial density
and temperature. That is, the critical ${\rm H_{2}}$ abundance for cooling
$x_{{\rm H_{2}},{\cal C}}$ is defined by

$t_{{\cal C}}(n_{V},x_{{\rm H_{2}},{\cal C}},T_{V})-t_{{\cal C}}(n_{V},x_{{\rm
H_{2}},{\cal C}},T_{{\rm H_{2}}})=t_{ff}(n_{V}),$ (4.2)

where the subscript $V$ stands for quantities in the initial, virial
equilibrium and $T_{{\rm H_{2}}}\approx 200~{}{\rm K}$ is the minimum
temperature achievable by ${\rm 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_{\rm H_{2}}$. Then, with the
${\rm H_{2}}$ production rate

$\dot{x}_{{\rm H_{2}}}=k_{\rm H,3}(T_{V})n_{{\rm H},V}x_{e,V},$ (4.3)

we have

$t_{{\rm H_{2}}}=\frac{x_{{\rm H_{2}},{\cal C}}}{\dot{x}_{{\rm H_{2}}}}.$
(4.4)

In Fig. 2, we find that the ${\rm HD}$ production indeed lowers the minimum
achievable temperature. Fig. 2 also shows that the ${\rm HD}$ production
saturates before $t_{{\rm col}}\approx 10t_{ff}$, which motivates us to set
$10t_{ff}$ as the maximum of $t_{{\rm col}}$.

We then extract the mass of the cloud from the $n$–$T$ trajectory as
$M_{c}=M_{J}(n_{\rm crit},T(n_{\rm crit})),$ where $n_{\rm crit}=1.75\times
10^{3}\,\rm cm^{-3}$ characterizes the loitering phase, the density at which
the ${\rm H_{2}}$ collisional de-excitation rates equal the spontaneous
radiative decay rates. For both large and small values of $f$ (corresponding
to efficient and inefficient ${\rm HD}$ formation, respectively) the critical
density also approximately corresponds to the minimum temperature. However,
for the transitional values of $f\sim 3$ the minimum temperature occurs at
higher densities as ${\rm HD}$ continues to form later in the collapse. Even
for these cases, we persist in extracting Jean’s mass from the temperature at
$n_{\rm crit}$ because, in order to affect the mass of the cloud, the ${\rm
HD}$ must be formed before ${\rm H_{2}}$ cooling becomes inefficient at
$n_{\rm crit}$.

Note that while Hirano et al. (2014) have discussed rotation as an important
determinant of $t_{\rm 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 proto-star once the cloud
mass is fixed.

###  4.2 Accretion Rate

The initial infall rate, as assumed in the KROME calculation [Eq. (3.2)], is

$\dot{M}_{\star}=M_{c}/t_{{\rm col}}.$ (4.5)

However, the accretion rate in later stages and on small scales is typically
limited by angular momentum transport. To model this, we adopt the
$\alpha$-disk parameterization (Shakura & Sunyaev, 1973) for the viscosity
$\nu=\alpha hc_{s}$ with $h$ being the scale height. For the thin disk,

$\frac{h}{R}\simeq\frac{c_{s}}{v_{c}}\,\to\,\nu=\alpha c_{s}^{2}/\Omega,$
(4.6)

where $R$ is disk radius, $v_{c}=\Omega R$ (with $\Omega$ the angular
velocity) is the circular velocity, the viscous timescale becomes

$t_{\nu}=\frac{R}{\nu/R}=\frac{R^{2}\Omega}{\alpha c_{s}^{2}}.$ (4.7)

Meanwhile, the angular momentum can be characterized by the spin parameter

$\beta\equiv\frac{v_{c}^{2}}{3|\Phi|}=\frac{\Omega^{2}R^{3}}{3GM_{c}},$ (4.8)

which is the ratio of rotational to gravitational energy in the cloud. In
terms of $\beta$, the viscous timescale is

$t_{\nu}=\frac{\sqrt{3\beta_{\rm crit}GM_{c}c_{s,\rm crit}t_{ff,\rm
crit}}}{\alpha c_{s}^{2}}$ (4.9)

where the subscript “${\rm crit}$” indicates the value of the quantity
evaluated at the critical density. Initially, $t_{{\rm col}}>t_{\nu}$. As the
density increases, $t_{{\rm col}}$ decreases as $1/\sqrt{n}$ while $t_{\nu}$
varies only due to the factor-of-few change in sound speed as the gas cools
and then heats. Hence, eventually $t_{\nu}>t_{\rm col}$. Therefore, the
accretion rate onto the protostar is ultimately set by the viscous timescale.
We estimate the accretion rate as

$\dot{M}_{\star}=M_{c}/t_{\nu}.$ (4.10)

In calculating $t_{\nu}$ we assume $T=1000\,\rm K$, typical of the molecular
disk. As the fiducial case we adopt $\beta=0.3$, which is a typical value in
Hirano et al. (2014). Here, we take $\alpha=1$. The typical values of $\alpha$
in Hirano et al. (2014) are a factor of a few lower, but we find that applying
$\alpha=1$ in Eq. (4.10) better matches the accretion rates in that work (see
App. A). We assume $\dot{M}_{\star}$ is constant, which is equivalent to
replacing $\dot{M}_{\star}(t)$ with its average value. Note that Liu et al.
(2020) gives a semi-analytic universal solution for Pop. III growth $M\propto
t^{4-3\gamma_{eff}}\approx t^{0.7}$ (with $\gamma_{eff}=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.

###  4.3 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:

$R_{\star}=26R_{\odot}\left(\frac{M_{\star}}{M_{\odot}}\right)^{0.27}\left(\frac{\dot{M}_{\star}}{10^{-3}\,\rm
M_{\odot}\,yr^{-1}}\right)^{0.41},$ (4.11)

which is surrounded by optically thick radiative precursor with a photospheric
radius of $R_{\rm ph}=1.4R_{\star}$. 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_{\star}^{3}$, and Hosokawa et al. (2012) find the following
approximate relationship:

$L_{\star}\simeq 10L_{\odot}\left(\frac{M_{\star}}{M_{\odot}}\right)^{3}\,.$
(4.12)

Two timescales are relevant in determining the contracting mass scale of the
protostellar core: the accretion timescale

$t_{\rm acc}=\frac{M_{\star}}{\dot{M}_{\star}},$ (4.13)

and the Kelvin-Helmholtz timescale,

$t_{\rm KH}=\frac{GM_{\star}^{2}}{R_{\star}L_{\star}}.$ (4.14)

Initially, the accretion timescale is short compared to the Kelvin-Helmholtz
timescale, and the protostar expands according to Eq. (4.11). Eventually,
however, the mass and luminosity growth of the protostar cause Kelvin-
Helmholtz contraction to dominate over accretion: $t_{KH}<t_{\rm acc}$, which
happens at $M_{\rm eq}$ (Hosokawa et al., 2012):

$M_{\rm eq}\simeq 15\,{\rm M_{\odot}}\left(\frac{\dot{M}_{\star}}{10^{-2}\rm
M_{\odot}yr^{-1}}\right)^{0.26}\,,$ (4.15)

at which point the protostellar cores with luminosity less than the Eddington
luminosity, $L(M_{\rm eq})<L_{\rm Edd}$, begin to contract. The following
energy balance equation can model the contraction:

$\frac{d}{dt}\left(\frac{W}{2}\right)=\frac{3}{2(5-n)}\frac{d}{dt}\left(\frac{GM_{\star}^{2}}{R}\right)=L(M_{\star}),$
(4.16)

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\propto\rho^{(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_{\star}\propto M_{\star}^{3}$,
we solve the energy balance equation (Eq. (4.16)) to find the radius-mass
relationship:

$R(M_{\star})=R_{\rm eq}\left[\frac{(M_{\star}/M_{\rm
eq})^{2}}{((M_{\star}/M_{\rm eq})^{4}-1)/3+1}\right]\,.$ (4.17)

Eventually, the luminosity of the collapsing proto-star reaches the Eddington
limit (Hosokawa et al., 2012),

$L_{\rm Edd}=3.8\times
10^{6}\left(\frac{M_{\star}}{100M_{\odot}}\right)L_{\odot}\,,$ (4.18)

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.16) 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_{\star}\propto M_{\star}^{0.5}$, which we adopt here. Note that Hosokawa et
al. (2012) find that for very high accretion rates $\gtrsim 0.1\,\rm
M_{\odot}/yr$ the proto-star begins to expand even before $L_{\rm 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 $\beta=0.3$. This omission will lead to an underestimate of the
masses of the most massive, slowly rotating stars.

Figure 3: The radius and effective temperature implied by Eq. (4.11) (radius
at $t_{eq}$), Eq. (4.17) (evolution to the Eddington luminosity), and
$R_{\star}\propto M_{\star}^{0.5}$ (upon reaching the Eddington limit).
Compare with Fig. 12 of Hirano et al. (2014) and Fig. 5 of Hosokawa et al.
(2012).

We show the evolution of radius and effective temperature of the protostars in
Fig. 3 for three different accretion rates.

###  4.4 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_{\rm 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 at $t_{\rm eq}$ (i.e. the time
when the mass reaches $M_{\rm eq}$) which is the characteristic epoch where
the protostar begins to emit significant UV radiation. 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 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 $\rho\propto r^{-2}$ as the initial
condition. There, $t_{B}$ is calculated by equating the ionizing photon flux
to the recombination rate behind the shock front:

$S_{\rm
EUV}(t_{B})=4\pi\alpha_{B}\int_{0}^{r_{sh}(t_{B})}drr^{2}n(r,t_{B})^{2},$
(4.19)

where $\alpha_{B}$ is the case-B recombination rate coefficient for hydrogen
at the temperature characteristic of the photo-heated gas $\sim 10^{4}\,\rm
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_{\rm SIS}$ of the initial SIS, which
we take to be the temperature at the loitering phase, i.e. $T_{\rm
SIS}=T(n_{\rm crit})$. When the shocked gas is much hotter than the
surrounding isothermal sphere (which holds in all cases we encounter), we have
$x_{s}\approx 2.56$.

Alvarez et al. (2006) provide the following approximate relationship for the
breakout time, which we find matches the exact result to within $\sim 10\%$:

$\begin{split}t_{B}&=6.5\times 10^{4}{\,\rm yr}\left(\frac{c_{s}x_{s}}{40\,\rm
km\,s^{-1}}\right)^{-1}\left(\frac{T_{\rm SIS}}{300\,\rm K}\right)^{2}\\\
&\times\left[\frac{S_{\rm EUV}(t_{B})}{3\times 10^{50}\,\rm
s^{-1}}\right]^{-1}.\end{split}$ (4.20)

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.20 by $x_{\rm
H}^{-2}=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_{\rm H}=0.927$. We finally solve Eq. (4.20) for $t_{B}$,
and then compute the final stellar mass as

$M_{\star}=M_{\rm eq}+\dot{M}_{\star}t_{B}.$ (4.21)

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) fit the redshift and
mass dependence of the stellar mass as

$M_{\star}=100{\,\rm
M_{\odot}\,}\left(\frac{1+z}{20}\right)^{3}\left(\frac{M_{\rm halo}}{3\times
10^{5}\,\rm M_{\odot}}\right)^{2}.$ (4.22)

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 $\beta=0.3$ as a function of halo mass and redshift. We have also plotted
the stellar sample of Hirano et al. (2014) (dots), their fit to the redshift
and halo mass dependence of the sample (dashed white lines) 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.

Figure 4: The predicted stellar mass in our model compared with the primordial
stars found in the simulations of Hirano et al. (2014) (colored circles) and
the fit of Hirano et al. (2014) to these stellar masses (dashed lines). 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.
Our model is accurate to within a factor of a few where Hirano et al. (2014)
have robust data (between $\sim 30$ and $\sim 300$ solar mass).

Our model succeeds quite well in predicting the transition between ${\rm
H_{2}}$ cooled clouds (final stellar mass of hundreds of solar masses) and
${\rm HD}$ cooled clouds (final stellar masses of order tens of solar masses).
Over the range where Hirano et al. (2014) have reasonable statistics ($\sim
30-300\,\rm M_{\odot}$), 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 $\beta=0.3$ in Fig. 4 (see Fig. 5). There are
several physically reasonable explanations which may contribute to this
discrepancy, which are explored in App. A. First, the rapidly accreting stars
of Hirano et al. (2014) can grow by around $100\,\rm M_{\odot}$ 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_{\nu}$ 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 $\sim 1000\,M_{\odot}$. 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 $\beta$ on a low mass, ${\rm HD}$
cooled cloud and a high mass, ${\rm H_{2}}$ cooled cloud at redshift 25. In
our model, the largest stellar mass attainable with a minimal realistic
rotation parameter $\beta\approx 0.05$ is $M\sim 200\,\rm M_{\odot}$, with an
accretion rate of $\sim 0.01\,\rm M_{\odot}/yr$. Rotation has a larger effect
at higher accretion rates (lower $\beta$), 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).

Figure 5: The role of the rotation parameter $\beta$ in the final stellar mass
for an ${\rm H_{2}}$ cooled ($10^{6}\,\rm M_{\odot}$, gold) cloud and ${\rm
HD}$ cooled ($10^{5}\,\rm M_{\odot}$, purple) cloud, both at $z=25$, showing
the sronger dependence on $\beta$ at higher cloud masses.

##  5 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_{{\rm
col}}=(\gamma-1)t_{{\cal 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 to within a factor of a few for stellar masses
between $\sim 30-300\,\rm M_{\odot}$. 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).

The model makes liberal use of the $\approx$ 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 $\alpha$, the
characteristic temperature of the disk used to evaluate the viscous timescale
$t_{\nu}$, and the characteristic temperature of the singular isothermal
sphere used in Eq. (4.19). The overall trends are robust to these choices.

We comment briefly that based only on the chemical-thermal evolution of the
gas (Sec. 3) two alternative estimates of the stellar mass are already
possible. First, the Jean’s mass at the loitering phase can be multiplied by
some global efficiency factor $M_{\star}=\epsilon M_{c}$, where to match the
characteristic $\sim 100\,\rm M_{\odot}$ Pop. III stellar mass $\epsilon\sim
0.1$. However, in our model the efficiency factor depends on the cloud mass.
We find that $\epsilon$ is a factor of a few larger for smaller, ${\rm HD}$
cooled clouds than for larger, ${\rm 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_{\star}=NM_{l}$, where
$M_{l}\approx
1.4M_{\odot}\mu^{-9/4}\left(\frac{k_{B}T}{m_{p}c^{2}}\right)^{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\sim
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_{f}$ on $T$ means that the temperature difference
between ${\rm H_{2}}$ and ${\rm HD}$ cooled clouds amounts to less than a
factor of two difference in $M_{l}$ and hence in $M_{\star}$. Additionally,
neither of these simpler estimates include the role of angular momentum.

Our model correctly determines the order-of-magnitude 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., 2020, 2023; Park et
al., 2023a, b), which is an intriguing direction for future research.

We thank Naoki Yoshida for helpful comments and for tracking down and sharing
the data from Hirano et al. (2014). This work was supported at Pennsylvania
State University by NASA ATP Program No. 80NSSC22K0819. DJ is also supported
by KIAS Individual Grant PG088301 at Korea Institute for Advanced Study. BL is
supported by the Royal Society University Research Fellowship. Research at
Perimeter Institute 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.

## References

  * Alvarez et al. (2006) Alvarez, M. A., Bromm, V., & Shapiro, P. R. 2006, The Astrophysical Journal, 639, 621, doi: 10.1086/499578
  * Bromm (2013) Bromm, V. 2013, Reports on Progress in Physics, 76, 112901, doi: 10.1088/0034-4885/76/11/112901
  * Bromm et al. (2002) Bromm, V., Coppi, P. S., & Larson, R. B. 2002, Astrophys. J., 564, 23, doi: 10.1086/323947
  * Bromm & Larson (2004) Bromm, V., & Larson, R. B. 2004, Annual Review of Astronomy and Astrophysics, 42, 79, doi: 10.1146/annurev.astro.42.053102.134034
  * Cooke et al. (2018) Cooke, R. J., Pettini, M., & Steidel, C. C. 2018, The Astrophysical Journal, 855, 102, doi: 10.3847/1538-4357/aaab53
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1, doi: 10.1016/j.physrep.2017.12.002
  * Eddington (1926) Eddington, A. S. 1926, The Internal Constitution of the Stars
  * Grassi et al. (2014) Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, Monthly Notices of the Royal Astronomical Society, 439, 2386, doi: 10.1093/mnras/stu114
  * Gurian et al. (2022) Gurian, J., Ryan, M., Schon, S., Jeong, D., & Shandera, S. 2022, The Astrophysical Journal Letters, 939, L12, doi: 10.3847/2041-8213/ac997c
  * Haemmerlé et al. (2020) Haemmerlé, L., Mayer, L., Klessen, R. S., et al. 2020, Space Science Reviews, 216, doi: 10.1007/s11214-020-00673-y
  * Hippert et al. (2022) Hippert, M., Setford, J., Tan, H., et al. 2022, Physical Review D, 106, doi: 10.1103/physrevd.106.035025
  * Hirano et al. (2015) Hirano, S., Hosokawa, T., Yoshida, N., Omukai, K., & Yorke, H. W. 2015, Mon. Not. R. Astron. Soc., 448, 568, doi: 10.1093/mnras/stv044
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, The Astrophysical Journal, 781, 60, doi: 10.1088/0004-637X/781/2/60
  * Hirata & Padmanabhan (2006) Hirata, C. M., & Padmanabhan, N. 2006, Monthly Notices of the Royal Astronomical Society, 372, 1175–1186, doi: 10.1111/j.1365-2966.2006.10924.x
  * Hollenbach & McKee (1979) Hollenbach, D., & McKee, C. F. 1979, Astrophys. J. Supp., 41, 555, doi: 10.1086/190631
  * Hosokawa et al. (2016) Hosokawa, T., Hirano, S., Kuiper, R., et al. 2016, Astrophys. J., 824, 119, doi: 10.3847/0004-637X/824/2/119
  * Hosokawa et al. (2012) Hosokawa, T., Omukai, K., & Yorke, H. W. 2012, The Astrophysical Journal, 756, 93, doi: 10.1088/0004-637X/756/1/93
  * Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250, doi: 10.1126/science.1207433
  * Johnson & Bromm (2006) Johnson, J. L., & Bromm, V. 2006, Mon. Not. R. Astron. Soc., 366, 247, doi: 10.1111/j.1365-2966.2005.09846.x
  * Klessen & Glover (2023) Klessen, R. S., & Glover, S. C. O. 2023, The first stars: formation, properties, and impact.  https://arxiv.org/abs/2303.12500
  * Latif et al. (2022) Latif, M. A., Whalen, D., & Khochfar, S. 2022, The Astrophysical Journal, 925, 28, doi: 10.3847/1538-4357/ac3916
  * Lipovka et al. (2005) Lipovka, A., Nú ñez-López, R., & Avila-Reese, V. 2005, Monthly Notices of the Royal Astronomical Society, 361, 850, doi: 10.1111/j.1365-2966.2005.09226.x
  * Liu et al. (2020) Liu, B., Meynet, G., & Bromm, V. 2020, Monthly Notices of the Royal Astronomical Society, 501, 643, doi: 10.1093/mnras/staa3671
  * Magnus et al. (2023) Magnus, L. C., Smith, B. D., Khochfar, S., et al. 2023, Formation Channels for Population III Stars at Cosmic Dawn.  https://arxiv.org/abs/2307.03521
  * McKee & Tan (2008) McKee, C. F., & Tan, J. C. 2008, The Astrophysical Journal, 681, 771, doi: 10.1086/587434
  * Nakazato et al. (2022) Nakazato, Y., Chiaki, G., Yoshida, N., et al. 2022, The Astrophysical Journal Letters, 927, L12, doi: 10.3847/2041-8213/ac573e
  * Nebrin et al. (2023) Nebrin, O., Giri, S. K., & Mellema, G. 2023, Monthly Notices of the Royal Astronomical Society, 524, 2290, doi: 10.1093/mnras/stad1852
  * Omukai & Nishi (1998) Omukai, K., & Nishi, R. 1998, The Astrophysical Journal, 508, 141, doi: 10.1086/306395
  * Park et al. (2023a) Park, J., Ricotti, M., & Sugimura, K. 2023a, Mon. Not. R. Astron. Soc., 521, 5334, doi: 10.1093/mnras/stad895
  * Park et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2307.14562, doi: 10.48550/arXiv.2307.14562
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, Astrophys. J., 187, 425, doi: 10.1086/152650
  * Rees (1976) Rees, M. J. 1976, Mon. Not. R. Astron. Soc., 176, 483, doi: 10.1093/mnras/176.3.483
  * Ripamonti (2007) Ripamonti, E. 2007, Monthly Notices of the Royal Astronomical Society, 376, 709, doi: 10.1111/j.1365-2966.2007.11460.x
  * Ryan et al. (2022) Ryan, M., Gurian, J., Shandera, S., & Jeong, D. 2022, Astrophys. J., 934, 120, doi: 10.3847/1538-4357/ac75ef
  * Ryan & Radice (2022) Ryan, M., & Radice, D. 2022, Phys. Rev. D, 105, 115034, doi: 10.1103/PhysRevD.105.115034
  * Seager et al. (1999) Seager, S., Sasselov, D. D., & Scott, D. 1999, The Astrophysical Journal, 523, L1–L5, doi: 10.1086/312250
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astron. Astrophys., 24, 337
  * Shandera et al. (2018) Shandera, S., Jeong, D., & Gebhardt, H. S. G. 2018, Physical Review Letters, 120, doi: 10.1103/physrevlett.120.241102
  * Sharda et al. (2020) Sharda, P., Federrath, C., & Krumholz, M. R. 2020, Monthly Notices of the Royal Astronomical Society, 497, 336, doi: 10.1093/mnras/staa1926
  * Sharda et al. (2021) Sharda, P., Federrath, C., Krumholz, M. R., & Schleicher, D. R. G. 2021, Monthly Notices of the Royal Astronomical Society, 503, 2014, doi: 10.1093/mnras/stab531
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, Monthly Notices of the Royal Astronomical Society, 323, 1, doi: 10.1046/j.1365-8711.2001.04006.x
  * Shu et al. (2002) Shu, F. H., Lizano, S., Galli, D., Canto, J., & Laughlin, G. 2002, The Astrophysical Journal, 580, 969, doi: 10.1086/343859
  * Singh et al. (2021) Singh, D., Ryan, M., Magee, R., et al. 2021, Physical Review D, 104, doi: 10.1103/physrevd.104.044015
  * Stacy & Bromm (2013) Stacy, A., & Bromm, V. 2013, Monthly Notices of the Royal Astronomical Society, 433, 1094, doi: 10.1093/mnras/stt789
  * Stacy et al. (2016) Stacy, A., Bromm, V., & Lee, A. T. 2016, Monthly Notices of the Royal Astronomical Society, 462, 1307, doi: 10.1093/mnras/stw1728
  * Stahler et al. (1986) Stahler, S. W., Palla, F., & Salpeter, E. E. 1986, Astrophys. J., 308, 697, doi: 10.1086/164542
  * Steigman (2007) Steigman, G. 2007, Annual Review of Nuclear and Particle Science, 57, 463, doi: 10.1146/annurev.nucl.56.080805.140437
  * Sugimura et al. (2020) Sugimura, K., Matsumoto, T., Hosokawa, T., Hirano, S., & Omukai, K. 2020, Astrophys. J. Lett., 892, L14, doi: 10.3847/2041-8213/ab7d37
  * Sugimura et al. (2023) —. 2023, arXiv e-prints, arXiv:2307.15108, doi: 10.48550/arXiv.2307.15108
  * Susa (2019) Susa, H. 2019, The Astrophysical Journal, 877, 99, doi: 10.3847/1538-4357/ab1b6f
  * Susa et al. (2014) Susa, H., Hasegawa, K., & Tominaga, N. 2014, Astrophys. J., 792, 32, doi: 10.1088/0004-637X/792/1/32
  * Tegmark et al. (1997) Tegmark, M., Silk, J., Rees, M. J., et al. 1997, The Astrophysical Journal, 474, 1, doi: 10.1086/303434
  * Yoshida et al. (2006) Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, Astrophys. J., 652, 6, doi: 10.1086/507978

##  Appendix A Detailed Comparison with Simulations

In addition to the fit to the final stellar mass [Eq. (4.22)], Hirano et al.
(2014) provide fits to the accretion rate given the cloud mass and $\beta$ and
to the final stellar mass given the accretion rate. These are:

$\displaystyle\dot{M_{\star}}$ $\displaystyle=2.8\times 10^{-3}\,{\rm
M_{\odot}\,yr^{-1}}\left(\frac{M_{c}}{350\,M_{\odot}}\cdot\frac{0.3}{\beta}\right)$
(A1) $\displaystyle M_{\star}$ $\displaystyle=100\,{\rm
M_{\odot}}\left(\frac{\dot{M}_{\star}}{2.8\times
10^{-3}\,M_{\odot}\,yr^{-1}}\right)^{0.8}.$ (A2)

We can apply these estimates to better understand the discrepancies between
our model and the results of Hirano et al. (2014), by replacing the relevant
pieces of our analytic model with these fits. This comparison will also
demonstrate the flexible, modular nature of our model. The fiducial model aims
for physical intuitiveness throughout the star formation process, but
depending on the application it can easily be modified to provide greater
fidelity to simulations.

First, we compare the accretion rate calculated in our model for $\beta=0.3$
with the sample of Hirano et al. (2014) and the corresponding fit Eq. (A1), in
Fig. 6. Compared to the simulation results, the analytic calculation
overestimates the accretion rate by a factor of a few for the least massive
clouds, and underestimates the accretion rate by a similar factor for the
heaviest.

Figure 6: The accretion rate from Eq. 4.10 (black) and from Eq. A1 (red,
dashed) and the sample of Hirano 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.

On the other hand, if we use our original accretion model Eq. (4.10) but
replace our analytic model of the stellar growth and ionizing feedback with
Eq. (A2), the largest effect is a factor of a few increase in the stellar mass
from the heavy, ${\rm H_{2}}$ cooled clouds (Fig. 8.)

Figure 7: The stellar mass calculated in our model using the accretion rate
Eq. (A1), compared to our fiducial result (Fig. 4). The range of stellar
masses is increased by a factor of a few.

Figure 8: The stellar mass calculated from Eq. (A2) with the accretion rate
supplied by Eq. (4.10). For most of the parameter space, the stellar mass is
enhanced compared to the fiducial model. This enhancement is strongest for the
most massive halos (and hence the most massive stars).

Similarly, we can use the accretion rate calculated in our model, Eq. (4.10),
but replace our treatment of the stellar evolution with the fit Eq. (A2). This
scenario is compared to the fiducial result in Fig. 8.

These comparisons show that both the accretion rate calculation and the
treatment of the subsequent stellar evolution and feedback contribute to the
discrepancies between our model and Hirano et al. (2014). The largest factor
arises from the stellar growth and feedback: Hirano et al. (2014) find that
rapidly accreting stars grow to higher masses than predicted by our model.
This is consistent with their finding that such stars can grow significantly
even after UV breakout, likely due to self-shielding in aspherical accretion
flows and early expansion of the protostar at very high accretion rates
($\gtrsim 0.1\ \rm M_{\odot}/yr$, Hosokawa et al., 2012), which are not
considered in our model.


end result of simplify still being a symbol, albeit simpler in form. Lastly,
inspired by a SymPy library function free_sym, we introduce another equivalent
function in our library, which we call symIdentifier; this function returns
all symbols found in expression entries of a given multidimensional array.

###  C.2 New functions for the computer library

In this appendix we discuss the computer library of our new functions which
augments the set of known functions in public libraries discussed in §C.1. We
briefly describe each new built-in function that is used in our algorithms in
§IV.2, and then we present the detailed information about these functions,
including types of input and output. These functions are presented as their
name in camelCase followed by a plain description of the function, then their
inputs and outputs.

  1. 1.

Oracle for group/algebra properties: Our library includes one oracular
function, accepting a binary input that determines the type of the group,
which is either a finite discrete group or a compact connected Lie group. Next
the input of this function is a name of a finite discrete group or a compact
connected Lie group. This oracular function also accepts an integer input
corresponding to Hilbert-space dimension $d$. This function returns pertinent
information about the properties of the input finite discrete group or the
input compact connected Lie group, depending on the type of the group
determined by the first input of the function. This oracular function yields
six outputs. Depending on the type of the group in the input, the first output
is the number of inequivalent irreps of the finite discrete group or the
number of inequivalent irreps of Lie algebra corresponding to the compact
connected Lie group with dimension less than or equal to $d$. The second
output is the number of inequivalent representations of the finite discrete
group or the number of inequivalent representations of the Lie algebra
corresponding to the compact connected Lie group, with respect to the group
type in the input, for given $d$. The third output is the rank of the finite
discrete group or the number of generators of the Lie algebra corresponding to
the compact connected Lie group. The fourth output is the dimension of the
inequivalent irreps of the finite discrete group or the dimension of all
inequivalent irreps of the compact connected Lie algebra corresponding to the
compact connected Lie group, with dimension less than or equal to $d$. Next we
have the fifth output, which is all inequivalent unitary irreps of the finite
discrete group or all inequivalent unitary irreps of the Lie algebra
corresponding to the compact connected Lie group with dimension less than or
equal to $d$. Finally, the sixth output is all unitary representations of the
finite discrete group or all unitary representations of the Lie algebra
corresponding to the compact connected Lie group, for the given $d$.

     * •

props $\triangleright$ Properties and representations of a discrete group or a
compact connected Lie group gName

       * INPUT:

         * binary[ ] gType $\triangleright$ Flag: $0$ for finite discrete groups and $1$ for compact connected Lie groups.

         * character[ ] gName $\triangleright$ Name of finite discrete group or compact connected Lie group.

         * posinteger hDim $\triangleright$ Hilbert-space dimension.

       * OUTPUT:

         * posinteger numIrrep $\triangleright$ Number of inequivalent irreps of the gName for gType=0 and number of inequivalent irreps of algebra generating gName for gType=1 with dimension less than or equal to hDim.

         * posinteger numRep $\triangleright$ Number of inequivalent reps of the gName for gType=0 and number of inequivalent irreps of algebra generating gName for gType=1 with dimension hDim.

         * posinteger numGen $\triangleright$ rank of gName for gType=0 and number of generators of the Lie algebra generating gName for gType=1.

         * posinteger[ ] dim $\triangleright$ Dimension of numIrrep inequivalent irreps of gName with dimension less than or equal to hDim for gType=0 and dimension of numIrrep inequivalent irreps of Lie algebra generating gName with dimension less than or equal to hDim.

         * symjagged[ ][ ][ ][ ] irrep $\triangleright$ For gType=0 entries of the first hyperrow label the numIrrep inequivalent irreps of gName with dimension less than or equal to hDim. Entries of the second hyperrow label numGen generators of gName. The third and fourth hyperrows are matrix elements for each irrep of gName with dimension less than or equal to hDim. For gType=1 entries of the first hyperrow label the numIrrep inequivalent irreps of Lie algebra generating gName with dimension less than or equal to hDim. Entries of the second hyperrow label numGen generators of the Lie algebra generating gName. The third and fourth hyperrows are matrix elements for each irrep of Lie algebra generating gName with dimension less than or equal to hDim.

         * symbol[ ][ ][ ][ ] rep $\triangleright$ For $\textsc{gType}=0/1$ entries of the first hyperrow label the numRep inequivalent reps of gName/Lie algebra generating gName with dimension hDim. Entries of the second hyperrow label numGen generators of gName/Lie algebra generating gName. The third and fourth hyperrows are matrix elements for each rep of gName/Lie algebra generating gName with dimension hDim.

  2. 2.

Imposing a trace preserving condition on completely positive maps: Our library
requires a function for accepting $N$ number of $d\times d$ Kraus operators of
a completely positive map and returning $N$ number of $d\times d$ Kraus
operators of the same completely positive map that satisfy the trace
preserving condition in Eq. (4). Given the set of Kraus operators, this
function constructs $\Xi$ in Eq. (5) and solves the matrix equation in Eq. (4)
for parameters in Kraus operators of the completely positive map. If the
solution exists, it employs the solution to simplify the Kraus operators and
returns 1 for TP, indicating that the trace-preserving condition is satisfied
in addition to the set of Kraus operators satisfying the trace-preserving
condition (4). If the solution does not exist, the function returns 0 for TP,
indicating that the solution does not exist.

     * •

solveChannel $\triangleright$ Apply the trace-preserving condition to the
input CP map to yield either a channel or no solution.

       * INPUT:

         * symbol[n][d][d] kraus $\triangleright$ n $\textsc{d}\times\textsc{d}$ Kraus matrices for CP map.

       * OUTPUT:

         * binary[ ] TP $\triangleright$ Flag: $0$ for no solution and $1$ for existence of the solution.

         * symbol[n][d][d] kraus $\triangleright$ n $\textsc{d}\times\textsc{d}$ Kraus matrices for channel.

       * PROCEDURE:

         1. (a)

symbol[d][d] xi $\triangleright$ Represents $\Xi$ (5).

         2. (b)

symbol[$\frac{\textsc{d}(\textsc{d}-1)}{2}$] eqs $\triangleright$ Constrained
expressions due to the trace-preserving condition.

         3. (c)

symbol[ ] symbols $\triangleright$ Yields parameters obtained from kraus.

         4. (d)

symbol[ ] sol $\triangleright$ Relations between parameters in kraus due to
the trace-preserving constraint; same size as symbols.

         5. (e)

binary flag $\triangleright$ TRUE if solutions exist.

         6. (f)

Compute xi $\triangleright$ From Eq. (5)

         7. (g)

For all $\textsc{i}\in[\textsc{d}]$,
eqs[i]$\leftarrow\textsc{xi}[\textsc{i}][\textsc{i}]-1$ $\triangleright$ Apply
trace-preserving condition for diagonal elements.

         8. (h)

For i from $\textsc{d}+1$ to $\frac{\textsc{d}(\textsc{d}-1)}{2}$, eqs gets
upper-diagonal elements of xi $\triangleright$ Apply trace-preserving
condition for upper-diagonal elements.

         9. (i)

symbols$\leftarrow$symIdentifier(kraus) $\triangleright$ Extracts symbols from
kraus.

         10. (j)

sol$\leftarrow$solve(eqs,params) $\triangleright$ Assigns eqs=0, then solves
for params; flag$\leftarrow$FALSE if a solution is not found.

         11. (k)

If flag, proceed to the next step.

         12. (l)

RETURN kraus$\leftarrow$ simplify(kraus,sol) $\triangleright$ Employ
expressions in sol to simplify kraus

  3. 3.

Solving systems of linear equations: The next function solves systems of
homogeneous linear equations symbolically; this system is $A_{n}\bm{x}=\bm{0}$
for $(d_{1}\times d_{2})$-dimensional symbolic matrices $A_{n}$s with
$n\in[N]$. Solving this system of equations is accomplished by computing the
intersection between kernels of $\left\\{A_{n}\right\\}$, which is
accomplished by solving a system of linear homogeneous equations. The matrix
of coefficients in this system of coupled linear equations is an $Nd_{1}\times
d_{2}$ matrix constructed by stacking $N$ instances of $A_{n}$ matrices. This
system of homogeneous linear equations is solved by standard methods such as,
singular value decomposition.

     * •

symSolve $\triangleright$ Solves symbolic $\bm{x}$ for
$A_{n}\bm{x}=\bm{0}\forall n\in[N]$, with each $A_{n}$ a
$\textsc{k}\times\textsc{l}$ symbolic matrix, and simplify employed to
simplify all algebraic expressions.

       * INPUT: symbol[n][k][l] a

       * OUTPUT: symbol[l] x

  4. 4.

Reshape: Our library incorporates a function that reshapes the input symbolic
vector into a square matrix, where the input vector has a length that is a
squared integer $d^{2}$ and the matrix has size $d\times d$. Reshaping is
accomplished by writing in order each element of the vector into each element
of the first row of the matrix until that row is full. Then we continue by
writing the next elements of the vector into the next row of the matrix until
that row is filled. This procedure is complete when all rows of the vector are
written into all elements of the square matrix, and we have ensured that our
matrix is exactly the right size for this transcription from a length $d^{2}$
vector to work properly.

     * •

reshape $\triangleright$ Reshape symbolic vector §C.1 but here only for a
vector of squared-integer length to a square matrix, converted according to
the rule that the first row fills the matrix, then the second row and so on.

       * INPUT: symbol$\left[\textsc{Dim}^{2}\right]$ vector

       * OUTPUT: symbol[Dim][Dim] matrix

  5. 5.

Convert a vector to multiple matrices: Our library furthermore requires a
function that converts a given vector to a set of square matrices, which
generalizes the previous function that maps a vector to a single square
matrix. The input is a symbolic vector of size $Kd^{2}\times 1$ yielding $K$
number of $d\times d$ matrices at the output. This function chops an input
vector of size $Kd^{2}$ into $K$ length $d^{2}$ vectors and then reshapes each
of these vectors into a square $d\times d$ matrix.

     * •

vecToMatr $\triangleright$ Chops a length-$KD^{2}$ symbolic vector into $K$
length $D^{2}$ vectors and reshape each vector to a square $D$-dimensional
symbolic.

       * INPUT:

         * symbol$\left[\textsc{k}*\textsc{d}^{2}\right]$ vector,

         * posinteger k

       * OUTPUT: symbol[k][d][d] matrixSeq

## References

  * Ludwig (1968) G. Ludwig, “Attempt of an axiomatic foundation of quantum mechanics and more general theories. iii.” Commun. Math. Phys 9, 1–12 (1968).
  * Hellwig and Kraus (1969) K. Hellwig and K. Kraus, “Pure operations and measurements,” Commun. Math. Phys 11, 214–220 (1969).
  * Doplicher  _et al._ (1971) S. Doplicher, R. Hagg, and J. E. Roberts, “Local observables and particle statistics i,” Commun.Math. Phys. 23, 199–230 (1971).
  * Holevo (2001) A. S. Holevo,  _Statistical Structure of Quantum Theory_ , Lecture Notes in Physics Monographs, Vol. 67 (Springer-Verlag, Berlin, 2001).
  * Shannon (1948a) C. E. Shannon, “A mathematical theory of communication,” Bell System Technical Journal 27, 379–423 (1948a).
  * Shannon (1948b) C. E. Shannon, “A mathematical theory of communication,” Bell System Technical Journal 27, 623–656 (1948b).
  * Landau and Streater (1993) L. J. Landau and R. F. Streater, “On birkhoff’s theorem for doubly stochastic completely positive maps of matrix algebras,” Linear Algebra Its Appl. 193, 107–127 (1993).
  * Choi (1975) Man-Duen Choi, “Completely positive linear maps on complex matrices,” Linear Algebra Its Appl. 10, 285–290 (1975).
  * Fujiwara and Algoet (1998) A. Fujiwara and P. Algoet, “Affine parameterization of quantum channels,” in  _Proc. 1998 IEEE International Symposium on Information Theory_ (1998) p. 87.
  * Fujiwara and Algoet (1999) Akio Fujiwara and Paul Algoet, “One-to-one parametrization of quantum channels,” Phys. Rev. A 59, 3290–3294 (1999).
  * Ruskai  _et al._ (2002) M. B. Ruskai, S. Szarek, and E. Werner, “An analysis of completely positive trace-preserving maps on $m_{2}$,” Linear Algebra Its Appl. 347, 159–187 (2002).
  * Braun  _et al._ (2014) Daniel Braun, Olivier Giraud, Ion Nechita, Clément Pellegrini, and Marko Žnidarič, “A universal set of qubit quantum channels,” J. Phys. A 47, 135302 (2014).
  * Mendl and Wolf (2009) Christian B. Mendl and Michael M. Wolf, “Unital quantum channels – convex structure and revivals of birkhoff’s theorem,” Commun. Math. Phys. 289, 1057–1086 (2009).
  * Haagerup  _et al._ (2020) U. Haagerup, M. Musat, and M. B. Ruskai, “Extreme points and factorizability for new classes of unital quantum channels,” arXiv:2006.03414 (2020).
  * Pérez-García  _et al._ (2006) David Pérez-García, Michael M. Wolf, Denes Petz, and Mary Beth Ruskai, “Contractivity of positive and trace-preserving maps under lp norms,” J. Math. Phys. 47, 083506 (2006).
  * King (2002) Christopher King, “Additivity for unital qubit channels,” J. Math. Phys. 43, 4641–4653 (2002), https://aip.scitation.org/doi/pdf/10.1063/1.1500791 .
  * Fukuda (2007) M. Fukuda, “Simplification of additivity conjecture in quantum information theory,” Quantum Inf. Process. 6 (2007).
  * Feynman (1982) R. P. Feynman, “Simulating physics with computers,” Int. J. Theor. Phys. 21, 467–488 (1982).
  * Lloyd (1996) Seth Lloyd, “Universal quantum simulators,” Science 273, 1073–1078 (1996).
  * Buluta and Nori (2009) I. Buluta and F. Nori, “Quantum simulators,” Science 326, 108–111 (2009).
  * Aharonov and Ta-Shma (2003) Dorit Aharonov and Amnon Ta-Shma, “Adiabatic quantum state generation and statistical zero knowledge,” in  _Proc. Thirty-Fifth Annual ACM Symposium on Theory of Computing_, STOC ’03 (Association for Computing Machinery, New York, USA, 2003) pp. 20–29.
  * Berry  _et al._ (2007) D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, “Efficient quantum algorithms for simulating sparse hamiltonians,” Commun. Math. Phys. 270, 359–371 (2007).
  * Childs (2010) A. M. Childs, “On the relationship between continuous- and discrete-time quantum walk,” Commun. Math. Phys 294, 581–603 (2010).
  * Wiebe  _et al._ (2011) Nathan Wiebe, Dominic W Berry, Peter Høyer, and Barry C Sanders, “Simulating quantum dynamics on a quantum computer,” J. Phys. A 44, 445308 (2011).
  * Sanders (2013) B. C. Sanders, “Efficient algorithms for universal quantum simulation,” in  _Reversible Computation. RC 2013_, Lecture Notes in Computer Science: Programming and Software Engineering, Vol. 7948, edited by Gerhard W. Dueck and D. Michael Miller (Springer, 2013) pp. 1–10.
  * Dawson and Nielsen (2006) Christopher M. Dawson and Michael A. Nielsen, “The solovay-kitaev algorithm,” Quantum Info. Comput. 6, 81–95 (2006).
  * Nielsen and Chuang (2010) Michael A. Nielsen and Isaac L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, Cambridge, 2010).
  * Wang  _et al._ (2013) Dong-Sheng Wang, Dominic W. Berry, Marcos C. de Oliveira, and Barry C Sanders, “Solovay-kitaev decomposition strategy for single-qubit channels,” Phys. Rev. Lett. 111, 130504 (2013).
  * Lu  _et al._ (2017) He Lu, Chang Liu, Dong-Sheng Wang, Luo-Kan Chen, Zheng-Da Li, Xing-Can Yao, Li Li, Nai-Le Liu, Cheng-Zhi Peng, Barry C. Sanders, Yu-Ao Chen, and Jian-Wei Pan, “Experimental quantum channel simulation,” Phys. Rev. A 95, 042310 (2017).
  * Wang and Sanders (2015) Dong-Sheng Wang and Barry C. Sanders, “Quantum circuit design for accurate simulation of qudit channels,” New J. Phys. 17, 043004 (2015).
  * Wang (2016) Dong-Sheng Wang, “Convex decomposition of dimension-altering quantum channels,” Int. J. Quantum Inform. 14, 1650045 (2016).
  * Iten  _et al._ (2017) Raban Iten, Roger Colbeck, and Matthias Christandl, “Quantum circuits for quantum channels,” Phys. Rev. A 95, 052316 (2017).
  * Scutaru (1979) H. Scutaru, “Some remarks on covariant completely positive linear maps on c∗-algebras,” Rep. Math. Phys. 16, 79–87 (1979).
  * Ruskai (2000) M. B Ruskai, “Some open problems in quantum information theory,” arXiv:0708.1902 (2000).
  * Holevo (2002) A. S. Holevo, “Remarks on the classical capacity of quantum channel,” arXiv:quant-ph/0212025 (2002).
  * König and Wehner (2009) Robert König and Stephanie Wehner, “A strong converse for classical channel coding using entangled inputs,” Phys. Rev. Lett. 103, 070504 (2009).
  * Datta  _et al._ (2016) N. Datta, M. Tomamichel, and M. M. Wilde, “On the second-order asymptotics for entanglement-assisted communication,” Quantum Inf. Process. 15, 2569–2591 (2016).
  * Wilde  _et al._ (2017) M. M. Wilde, M. Tomamichel, and M. Berta, “Converse bounds for private communication over quantum channels,” IEEE Trans. Inf. Theory 63, 1792–1817 (2017).
  * Das  _et al._ (2020) Siddhartha Das, Stefan Bäuml, and Mark M. Wilde, “Entanglement and secret-key-agreement capacities of bipartite quantum interactions and read-only memory devices,” Phys. Rev. A 101, 012344 (2020).
  * Karimipour  _et al._ (2011) Vahid Karimipour, Azam Mani, and Laleh Memarzadeh, “Characterization of qutrit channels in terms of their covariance and symmetry properties,” Phys. Rev. A 84, 012321 (2011).
  * Mozrzymas  _et al._ (2017) Marek Mozrzymas, Michał Studziński, and Nilanjana Datta, “Structure of irreducibly covariant quantum channels for finite groups,” J. Math. Phys. 58, 052204 (2017), https://doi.org/10.1063/1.4983710 .
  * Siudzińska and Chruściński (2018) Katarzyna Siudzińska and Dariusz Chruściński, “Quantum channels irreducibly covariant with respect to the finite group generated by the weyl operators,” J. Math. Phys. 59, 033508 (2018).
  * Datta  _et al._ (2006) N. Datta, M. Fukuda, and A. S. Holevo, “Complementarity and additivity for covariant channels,” Quantum Inf. Process. 5, 179–207 (2006).
  * Arora and Barak (2009) S. Arora and B. Barak,  _Computational Complexity: A Modern Approach_ (Cambridge University Press, Cambridge, 2009).
  * Coxeter and Moser (1972) H. S. M. Coxeter and W. O. J. Moser,  _Generators and Relations for Discrete Groups_, Vol. 14 (Springer-Verlag Berlin Heidelberg, 1972).
  * Barut and Raczka (1986) A. Barut and R. Raczka,  _Theory of Group Representations and Applications_ (WORLD SCIENTIFIC, 1986).
  * Herstein (1975) I. N. Herstein,  _Topics in Algebra_ (John Wiley & Sons, 1975).
  * Young (1900) A. Young, “On quantitative substitutional analysis,” Proc. London Mathematical Society 33, 97–145 (1900).
  * Karimipour and Memarzadeh (2008) Vahid Karimipour and Laleh Memarzadeh, “Matrix product representations for all valence bond states,” Phys. Rev. B 77, 094416 (2008).
  * Duan and Guo (1998) L. M. Duan and G. C. Guo, “Reducing decoherence in quantum-computer memory with all quantum bits coupling to the same environment,” Phys. Rev. A 57, 737–741 (1998).
  * Macchiavello and Palma (2002) C. Macchiavello and G. M. Palma, “Entanglement-enhanced information transmission over a quantum channel with correlated noise,” Phys. Rev. A 65, 050301(R) (2002).
  * Holevo (2006) A. S. Holevo, “The additivity problem in quantum information theory,” in  _Proc. International Congress of Mathematicians_ (European Mathematical Society, 2006) pp. 999–1018.
  * Hastings (2009) M. B. Hastings, “Superadditivity of communication capacity using entangled inputs,” Nat. Phys. 5, 255–257 (2009).
  * Blum  _et al._ (1989) L. Blum, M. Shub, and S. Smale, “On a theory of computation and complexity over the real numbers: Np-completeness, recursive functions and universal machines,” Bull. Am. Math. Soc. 21, 1–46 (1989).
  * Jeong  _et al._ (2013) Youn-Chang Jeong, Jong-Chan Lee, and Yoon-Ho Kim, “Experimental implementation of a fully controllable depolarizing quantum operation,” Phys. Rev. A 87, 014301 (2013).
  * Shaham and Eisenberg (2012) Assaf Shaham and Hagai S. Eisenberg, “Realizing a variable isotropic depolarizer,” Opt. Lett. 37, 2643–2645 (2012).
  * Sciarrino  _et al._ (2004) F. Sciarrino, C. Sias, M. Ricci, and F. De Martini, “Realization of universal optimal quantum machines by projective operators and stochastic maps,” Phys. Rev. A 70, 052305 (2004).
  * Turing (1937) A. M. Turing, “On computable numbers, with an application to the entscheidungsproblem,” in  _Proc. London Mathematical Society_, 2, Vol. 42 (1937) pp. 230–265.
  * Yao (1993) A. Chi-Chih Yao, “Quantum circuit complexity,” in  _Proc. 1993 IEEE 34th Annual Foundations of Computer Science_, SFCS ’93 (IEEE Computer Society, USA, 1993) pp. 352–361.
  * Molina and Watrous (2019) A. Molina and J. Watrous, Proc. R. Soc. A 475, 20180767 (2019).
  * Parnas  _et al._ (1976) D. L. Parnas, John E. Shore, and David Weiss, “Abstract types defined as classes of variables,” in  _Proc. 1976 Conference on Data: Abstraction, Definition and Structure_ (Association for Computing Machinery, New York, USA, 1976) pp. 149–154.
  * Jafarzadeh  _et al._ (2020) M. Jafarzadeh, Ya-Dong Wu, Y. R. Sanders, and B. C. Sanders, “Randomized benchmarking for qudit clifford gates,” New J. Phys. 22, 063014 (2020).
  * Knuth (1997) Donald Ervin Knuth,  _The Art of Computer Programming_ , 3rd ed., Vol. 1 (Addison Wesley Longman, Reading, MA, 1997).
  * Wilkes  _et al._ (1951) Maurice Wilkes, David. Wheeler, and Stanley Gill,  _The preparation of programs for an electronic digital computer: with special reference to the EDSAC and the use of a library of subroutines_ (Addison-Wesley Press, Boston, 1951).
  * Meurer  _et al._ (2017) Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz, “SymPy: Symbolic Computing in Python,” PeerJ Computer Science 3, e103 (2017).
  * Oliphant (2015) Travis Oliphant,  _Guide to NumPy_ (CreateSpace Independent Publishing Platform, 2015).
  * (67) The NAG Library, The Numerical Algorithms Group (NAG), Oxford, United Kingdom, www.nag.com .