# Characterization of Population III Stars with Stellar Atmosphere and
Evolutionary Modeling and Predictions of their Observability with the James
Webb Space Telescope

Mikaela M. Larkin  Center for Astrophysics and Space Sciences, University of
California San Diego, La Jolla, CA 92093, USA  Roman Gerasimov Center for
Astrophysics and Space Sciences, University of California San Diego, La Jolla,
CA 92093, USA  Adam J. Burgasser Center for Astrophysics and Space Sciences,
University of California San Diego, La Jolla, CA 92093, USA

(Received TBD; Accepted TBD)

###### Abstract

Population III stars were the first stars to form after the Big Bang, and are
believed to have made the earliest contribution to the metal content of the
universe beyond the products of the Big Bang Nucleosynthesis. These stars are
theorized to have had extremely short lifespans, and therefore would only be
observable at high redshifts ($z\geq 3-17$) and faint apparent magnitudes
($m_{AB}\gtrsim 40$). The direct detection of Population III stars therefore
remains elusive. However, the recently launched James Webb Space Telescope
(JWST) may be capable of detecting stars in the relevant magnitude range in
the event of favorable gravitational lensing. Theoretical models are required
to interpret these future observations. In this study, new evolutionary models
and non-equilibrium model atmospheres were used to characterize the observable
properties of zero-age main sequence Population III stars. The calculated
models cover a wide range of possible Population III stellar masses, from the
minimum mass predicted by star formation studies to the maximum mass capable
of maintaining hydrostatic equilibrium. Synthetic photometry and theoretical
color-magnitude diagrams were calculated for the bands of the Near-Infrared
Camera (NIRCam) on JWST. The final results are compared to the scales of known
lensing events and JWST magnitude limits. The purpose of this study is to
calculate the observable parameters of Population III stars in the most
optimal JWST bands in order to provide a theoretical foundation for
anticipated future observations of this stellar population.

Population III stars (1285) — Gravitational lensing (670) — Theoretical models
(2107) — Limiting magnitude (923) — Proton-proton cycle (1299) — CNO cycle
(194)

††journal: ApJ††facilities: JWST (NIRCam, MIRI)††software:  Astropy (Astropy
Collaboration et al., 2013, 2018), Matplotlib (Hunter, 2007), NumPy (Harris et
al., 2020), SciPy (Virtanen et al., 2020), ATLAS (Kurucz, 1970; Sbordone et
al., 2004; Castelli, 2005a; Kurucz, 2014), DFSYNTHE (Castelli, 2005b), SYNTHE
(Kurucz & Avrett, 1981), MESA (Paxton et al., 2011, 2013, 2015, 2018, 2019)

tablenum SIXtablenum

##  1 Introduction

Population III stars are an elusive addition to the two traditional stellar
populations identified in Baade (1944). This yet unobserved population
accommodates the earliest and nearly metal-free stars that formed shortly
after the Big Bang. Nearly seven decades ago, Schwarzschild & Spitzer (1953)
recognized that early stars were likely more massive, luminous and short-lived
compared to their present-day Population I and II counterparts. The need for a
non-standard formation mechanism to produce such stars and the relatively high
metallicity measured in the most metal-poor Population II stars known at the
time ($[\mathrm{Fe}/\mathrm{H}]\approx-3$, Bond 1980) led Antonov & Chernin
(1977)111Note that the argument in Antonov & Chernin (1977) is based on the
assumption that dark matter is composed of low-mass stars and brown dwarfs
which has been largely ruled out (Bahcall et al., 1994; Graff & Freese, 1996;
Freese, 2017). and Wagner (1978) to independently identify the earliest stars
as a distinct stellar population. In addition to the earliest chemical
enrichment of the primordial gas, Population III stars may have contributed to
the cosmic microwave background (Rowan-Robinson, 1983), the cosmic infrared
background (Santos et al., 2002; Kashlinsky et al., 2005; Madau & Silk, 2005;
Kashlinsky, 2005), the gravitational wave background (Suwa et al., 2007;
Kowalska et al., 2012; Kinugawa et al., 2014), reionization and reheating of
the universe (Haiman & Loeb, 1997; Cojazzi et al., 2000; Cen, 2003a; Wyithe &
Loeb, 2003; Sokasian et al., 2004; Maio et al., 2016) and likely had a
noticeable feedback effect on the formation of the first galaxies (Greif et
al., 2010; Abe et al., 2021).

Recent searches for extremely metal-poor stars have revealed a sharp cut-off
in the metallicity distribution around $[\mathrm{Fe}/\mathrm{H}]\approx-4$,
with fewer than $50$ known objects below this cut-off (Abohalima & Frebel,
2018). These include the most metal-poor star known
($[\mathrm{Fe}/\mathrm{H}]=-6.2$, Nordlander et al. 2019) as well as at least
one object with an unmeasured iron abundance and the estimated upper limit
$[\mathrm{Fe}/\mathrm{H}]\lesssim-7$ (Keller et al., 2014). However, the
comparatively high abundances of other chemical elements (e.g.
$[\mathrm{C}/\mathrm{H}]\gg-3$ for both stars) challenge the status of these
sources as true Population III representatives, unless the observed abundances
are acquired from the interstellar medium through selective accretion, as
proposed by Johnson (2015).

The lack of identifiable metal-free stars in surveys is consistent with the
long-standing theoretical expectation of a top-heavy initial mass function
(IMF) that precludes primordial stars from having sufficiently long lifespans
to survive until the present day. This expectation is motivated by the lack of
metals in the early universe, which leaves $\mathrm{H}_{2}$ and $\mathrm{HD}$
as the only available cooling agents in molecular clouds, thereby suppressing
fragmentation and producing unusually heavy stars (Larson, 1998; Maio et al.,
2010). Early numerical simulations (e.g. Abel et al. 2000; Bromm et al. 2002;
Yoshida et al. 2006) suggested that Population III stars predominantly formed
with masses in excess of $100\ M_{\odot}$ and collapsed into black holes, with
the exception of a subset of stars with masses between $140\ M_{\odot}$ and
$260\ M_{\odot}$ that produced pair-instability supernovae (Schneider et al.,
2002; Heger & Woosley, 2002). Later studies with a more detailed treatment of
radiative feedback, interactions between stars, turbulence, etc., have
challenged this picture by allowing formation of $\gtrsim 10\ M_{\odot}$
Population III stars (e.g. Hosokawa et al. 2011; Stacy et al. 2012; Hirano et
al. 2015, 2014) or even $<1\ M_{\odot}$ stars (e.g. Stacy et al. 2016; Susa et
al. 2014; Greif et al. 2011; Clark et al. 2011) that may exist in the present-
day Milky Way, likely disguised by metal-enriched mass transfer from heavier
stars (Stacy & Bromm, 2014). The existence of lower-mass primordial stars
enables additional mechanisms of chemical enrichment, potentially explaining
the observed inconsistency of the abundance patterns in metal-poor stars with
predictions of pair-instability supernova yields (Umeda & Nomoto, 2002). A
distinct hypothetical population of supermassive primordial stars with masses
in excess of $10^{3}-10^{6}\ M_{\odot}$ has also been proposed as seeds for
the supermassive black holes in high-redshift quasars (Fuller et al., 1986;
Muñoz et al., 2021; Woods et al., 2021; Herrington et al., 2022).

Since the majority of Population III stars are expected to have masses between
a few tens and a few hundred solar masses, the correspondingly short lifespans
($\lesssim 20\ \mathrm{Myr}$, Windhorst et al. 2018) necessitate direct
observation of such objects at high redshifts. The first stars begin to form
once the primordial molecular clouds, concentrated around growing dark matter
over-densities, cool down sufficiently to become unstable against
gravitational collapse. The star formation rate as a function of redshift can
be traced in simulations of cosmological hydrodynamics (Xu et al., 2013; Wise
et al., 2012b, a; Jaacks et al., 2019). The first Population III stars are
expected to form around $z\approx 30$ ($\sim 0.1\ \mathrm{Gyr}$ after the Big
Bang) and the maximum formation rate density ($\sim 10^{-4}-10^{-3}\
M_{\odot}\ \mathrm{yr}^{-1}\ \mathrm{Mpc}^{-3}$) is attained around
$17\lesssim z\lesssim 10$ ($\sim 0.2-0.5\ \mathrm{Gyr}$ after the Big Bang).
This result is generally consistent with early reionization optical depth
measurements from the Wilkinson Microwave Anisotropy Probe (WMAP) (Spergel et
al., 2003; Kogut, 2003; Cen, 2003b), although later WMAP (Bennett et al.,
2013) and Planck (Planck Collaboration et al., 2016) measurements have cast
doubt on the usability of this parameter as a probe of primordial star
formation (Yung et al., 2020). Population III stars give way to Population II
stars once the metal mass fraction ($Z$) of the interstellar medium reaches a
critical value, $Z_{\mathrm{cr}}\gtrsim 10^{-8}-10^{-6}$
($Z_{\mathrm{cr}}\gtrsim 10^{-6}Z_{\odot}-10^{-4}Z_{\odot}$222Solar
metallicity taken as $Z_{\odot}=0.01$ to the nearest order of magnitude.;
Bromm et al. 2001a; Omukai et al. 2005; Schneider et al. 2006; Clark et al.
2008; Sharda & Krumholz 2022) that allows for efficient cooling and
fragmentation of collapsing molecular clouds. The lowest redshift at which
Population III stars may be observed remains uncertain as isolated metal-free
pockets may last for extended periods of time, producing new Population III
stars at later epochs. Population III star formation is expected to continue
until at least $z=6$ ($\sim 1\ \mathrm{Gyr}$ after the Big Bang, Trenti et al.
2009; Muratov et al. 2013) and possibly much later under special circumstances
(Liu & Bromm, 2020). Searching for metal-free stars at redshifts as low as
$z=3$ ($\sim 2\ \mathrm{Gyr}$ after the Big Bang) is particularly important
given observations of the Lynx arc – a star forming region at $z=3.4$ with a
Population III-consistent ionization source (Fosbury et al., 2003); and
LLS1249 – a dense gas cloud at $z=3.5$ with a Population III remnant
metallicity (Crighton et al., 2016).

A renewed interest in direct detection of individual Population III stars has
developed in anticipation of observations with the next generation of ground-
and space-based facilities, in particular the recently launched James Webb
Space Telescope (JWST). Using the Main Sequence properties of Population III
stars from the evolutionary calculations in Schaerer (2002) and new metal-free
model atmospheres calculated with the TLUSTY code (Hubeny & Lanz, 1995),
Rydberg et al. (2013) estimated the observable properties of isolated
primordial stars at a range of redshifts and stellar masses. Without
gravitational lensing, Population III stars appear far too faint to be
detected even at the lowest redshift ($z=2$) and in extremely long exposures
($100\ \mathrm{hr}$). Rydberg et al. (2013) also considered the case of a
favorable lensed observation through the galaxy cluster MACS J0717.5$+$3745 –
one of the largest gravitational lenses known (Zitrin et al., 2009). Even in
the lensed case, a realistic detection was found to require either an
extremely heavy Population III candidate ($\geq\ 300\ M_{\odot}$) or a
primordial star formation rate, far in excess of theoretical expectation.

However, very high magnifications ($\mu$) may be attained for brief periods of
time in the event of gravitational lensing during a caustic transit. A caustic
is the locus of points in the source plane where the determinant of the
magnification matrix vanishes, i.e. where a true point source would experience
infinite magnification (Narayan & Wallington, 1992). A compact light source
such as an individual star crossing a caustic may experience extreme
magnification up to $\mu\sim 10^{7}$ from a lensing cluster with a continuous
distribution of mass under most favorable conditions (Miralda-Escude, 1991).
In practice, microlenses within the galaxy cluster will distort the lens
caustics, reducing the maximum magnification to $\mu\sim 10^{4}$ (Diego et
al., 2018; Diego, 2019).

Since larger magnifications require more favorable and increasingly less
likely configurations, the true expected magnification in any given survey
will strongly depend on the redshift of interest, the number of observable
targets, and the scope of the survey itself. For example, Zackrisson et al.
(2015) calculate $\mu\gtrsim 700$ as a realistic magnification estimate for
detecting Population III stars in a $100\ \mathrm{deg}^{2}$ ultra-deep survey.
Extreme lensing events with $\mu\gg 10^{3}$ have allowed for recent
discoveries of the most distant individual stars known at $z=1.5$ (“Icarus”,
$\mu>2000$, Kelly et al. 2018), $z=2.7$ (currently unnamed, $\mu\gtrsim
10^{4}$, Chen et al. 2022) and $z=6.2$ (“Earendel”, $\mu>4000$, Welch et al.
2022a, b). It has been suggested that the last source may in fact be a
Population III star (Schauer et al., 2022).

Adopting more optimistic magnifications of $\mu\sim 10^{4}-10^{5}$, Windhorst
et al. (2018) used new metal-free evolutionary tracks calculated with the MESA
(Modules for Experiments in Stellar Astrophysics) code (Paxton et al., 2011,
2013, 2015, 2018, 2019) and assumed blackbody emergent spectra to characterize
Population III stars in the context of future observations with JWST. The
study estimated that a decade-long observational program monitoring up to $30$
candidate lensing clusters will be required for a reliable detection.

In this study, we contribute to the ongoing effort of predicting future
observations of Population III stars using stellar modelling. In particular,
we focus on primordial stars with initial masses between $1\ M_{\odot}$ and
$\sim 10^{3}\ M_{\odot}$ on the Zero Age Main Sequence (ZAMS). ZAMS properties
of Population III stars provide a lower bound on observability since later
evolutionary stages are intrinsically more luminous and subjected to less
interstellar absorption due to redder colors (Schaerer, 2002; Windhorst et
al., 2018). Additionally, the lower surface gravities of post-ZAMS stars may
lead to super-Eddington luminosities in high-mass candidates, requiring
detailed modelling of mechanical motion in the atmosphere that falls beyond
the scope of this study (however, see Liu et al. 2021; Yoon et al. 2012).

We present theoretical color-magnitude diagrams of ZAMS Population III stars
in JWST NIRCam filters for a broad range of redshifts based on new metal-free
evolutionary models and model atmospheres. The observable properties of
primordial stars are analyzed for their dependence on individual opacity
sources and non-equilibrium distribution of the radiation field throughout the
atmosphere. In particular, we demonstrate that even in the absence of non-grey
opacity sources, Population III atmospheres cannot be approximated as
blackbodies and always require detailed modelling. The calculated metal-free
physical parameters on ZAMS are compared to their metal-poor counterparts and
matched to simple, physically-motivated analytic relationships. The models are
also evaluated in the context of the Eddington limit, that is of particular
importance at high initial masses.

In this paper, Section 2 describes our modelling toolkit and presents the new
atmosphere and evolutionary models calculated in this study. The key physical
properties of Population III stars inferred from the models, such as the
dependence of stellar evolution on the dominant energy production mechanism
and the Eddington limit, are discussed in this section as well. Section 3
details our color-magnitude calculations at high redshift and presents the
predicted color-magnitude and mass-magnitude relationships for Population III
stars in the context of future JWST observations with gravitational lensing.
The effect of cosmological parameters on our predictions is also estimated in
Section 3. The key findings and important shortcomings of this investigation
are summarized in Section 4. Population III model parameters are tabulated in
Appendix A.

##  2 Modeling

###  2.1 Overview of the methodology

In this study, predictions of colors and magnitudes of Population III stars
are drawn from synthetic photometry of metal-free stellar models at ZAMS. Each
model is parameterized exclusively by the initial stellar mass and must
account for all relevant physical processes governing the evolution of the
star from its pre-Main Sequence (PMS) phase until steady-state hydrogen
fusion. Synthetic photometry is obtained from the evolved emergent spectrum of
the star, which is, in turn, calculated by solving the radiative transfer
equation at every wavelength throughout the outer layers of the model
comprising the stellar atmosphere.

Stellar atmospheres are particularly challenging to model due to the presence
of neutral and partially ionized species, resulting in complex, wavelength-
dependent opacity from significant non-grey contributions of bound-free and
bound-bound sources. At high effective temperatures and extremely low
metallicities considered in this work, the effect of non-grey atmospheric
opacity on the structure of the stellar interior is expected to be
insignificant, avoiding the need for detailed opacity calculations in the
evolutionary models. However, the atmospheric opacity must be re-introduced
into the model when calculating the final emergent spectrum of the star. We
therefore calculate all models in multiple stages. First, evolutionary
modelling was carried out from PMS to ZAMS with grey atmospheric opacity. The
evolved stellar radii and luminosities were used to derive simple analytic
relationships between stellar mass and the ZAMS surface parameters (effective
temperature, $T_{\mathrm{eff}}$, and surface gravity, $\log_{10}(g)$).
Finally, the analytic relationships were evaluated at a broad range of stellar
masses (from $1\ M_{\odot}$ to $1000\ M_{\odot}$) and the resulting surface
parameters were used as inputs to dedicated model atmosphere calculations with
the full opacity treatment and spectral synthesis.

###  2.2 Evolutionary Modeling

We calculated all evolutionary models with the MESA code (Paxton et al., 2011,
2013, 2015, 2018, 2019), version 15140. Evolutionary calculations in MESA are
carried out in adaptive time steps until the chosen termination condition is
met. By default, structure equations at each step are solved using the grey
atmosphere optical depth-temperature relationship as the surface boundary
condition. As demonstrated in our previous work (Gerasimov et al., 2022), this
approximation is satisfactory at effective temperatures ($T_{\mathrm{eff}}$)
above $5000\ \mathrm{K}$ – a condition met by all models considered in this
study.

In all evolutionary models, we chose to use $Y=0.25$ as the approximation for
the primordial helium mass fraction based on the Planck measurement in Planck
Collaboration et al. (2020), taken to the nearest percent to match the default
precision in MESA.

At initial masses below $90\ M_{\odot}$, models are initialized as PMS with a
central temperature of $5\times 10^{5}\ \mathrm{K}$ (following Choi et al.
2016; Gerasimov et al. 2022), uniform chemical composition, and a density
profile that satisfies the structure equations and results in the desired
stellar mass. At significantly larger initial masses, a convergent PMS
solution may not exist. Instead, the $90\ M_{\odot}$ PMS is used as the
starting point and the mass is slowly increased to the required value using
the mass relaxation routine provided by MESA. As an example, at the initial
mass of $1000\ M_{\odot}$, the relaxation process lasts $\approx 1300\
\mathrm{yr}$ and results in an object with the central temperature of
$72.7\times 10^{6}\ \mathrm{K}$.

We extract the ZAMS luminosity and stellar radius from all evolutionary models
once their nuclear power output begins to exceed $90\%$ of the total
luminosity. The sensitivity of the derived ZAMS parameters on the primordial
helium mass fraction ($Y$) as well as the adopted PMS settings was estimated
by computing multiple models at the representative initial masses of $10\
M_{\odot}$ and $1000\ M_{\odot}$ for a range of $Y$ values from $0.24$ to
$0.26$; a range of initial central temperatures from $3\times 10^{5}\
\mathrm{K}$ to $7\times 10^{5}\ \mathrm{K}$; and a range of maximum (pre-
relaxation) PMS masses from $50\ M_{\odot}$ to $100\ M_{\odot}$ (stored in the
max_mass_to_create variable of the build_pre_ms_model() subroutine). We found
the effect of PMS settings to not exceed $0.002\ \mathrm{dex}$ in both
luminosity and radius at ZAMS. The effect of the chosen $Y$ value was slightly
larger, reaching $0.02\ \mathrm{dex}$ at lower masses. However, neither of the
aforementioned uncertainties exceeds the average accuracy of the calculated
analytic mass-radius and mass-luminosity relationships (to be introduced
below) that were estimated as $0.02\ \mathrm{dex}$ and $0.03\ \mathrm{dex}$
respectively. Therefore, our results are insensitive to the input parameters
within the aforementioned ranges.

We calculated $428$ evolutionary models with initial stellar masses ranging
from $1\ M_{\odot}$ to $1000\ M_{\odot}$ and used the ZAMS radii and
luminosities to derive analytical mass-radius and mass-luminosity
relationships for Population III stars. The input settings (inlist files) for
all evolutionary models are available
https://zenodo.org/record/7145568#.Yzzax33MLAU (catalog
online)333https://zenodo.org/record/7145568#.Yzzax33MLAU. Radii of main
sequence stars are generally well-described by power law relations assuming
that the dominant energy production and transport mechanisms do not vary
significantly (Gimenez & Zamorano, 1985; Lacy, 1979; Demircan & Kahraman,
1991). For stars with solar metallicity, the power law index changes
noticeably around $\sim 1\ M_{\odot}$ due to the dissipation of the outer
convective zone (Kippenhahn & Weigert, 1994, Ch. 22.1) and the transition of
the main hydrogen fusion mechanism from the proton-proton chain to the carbon-
nitrogen-oxygen (CNO) cycle (Salaris & Cassisi, 2005, Ch. 5.2). The latter
effect is particularly important for zero-metallicity stars, as it is expected
that enough carbon will be produced at sufficiently high masses to display a
similar transition, thereby offsetting the power law break into the range of
masses considered in this study. We therefore model the mass-radius
relationship of ZAMS Population III stars as a broken power law with the break
point mass ($M_{\mathrm{bp}}^{(R)}$) treated as a free parameter. The adopted
relationship is

$R\propto\begin{cases}M^{\alpha},&\text{if }M\leq M_{\mathrm{bp}}^{(R)}\\\
M^{\beta},&\text{if }M>M_{\mathrm{bp}}^{(R)}\\\ \end{cases}$ (1)

where $M$ is the initial stellar mass, $R$ is the corresponding stellar
radius, and $\alpha$ and $\beta$ are the proton-proton and CNO power indices
respectively. The mass-luminosity ($M$–$L$) relationship is slightly more
complicated due to the dependence on the dominant pressure support in the
interior. High-mass ($M\gtrsim 10\ M_{\odot}$) stars behave approximately as
Eddington standard models ($n=3$ polytrope; Eddington 1918; Köhler et al.
2015) with the mass-luminosity power index gradually changing from $\geq 3$ at
$M\lesssim 100\ M_{\odot}$ to the asymptotic linear relationship ($L\propto
M$) in the limit of $M\rightarrow\infty$. The transition occurs as the
equation of state in the interior changes from the ideal gas law to radiation
pressure dominance. To accommodate this additional complexity, we allow the
index of the power law to change linearly with $\log_{10}(M)$ below some break
point mass ($M_{\mathrm{bp}}^{(L)}$), which is treated as a free parameter, as
before. The relationship is

$L\propto\begin{cases}M^{\delta\log_{10}\left(M/M_{\odot}\right)+A},&\text{if
}M\leq M_{\mathrm{bp}}^{(L)}\\\ M^{\gamma},&\text{if
}M>M_{\mathrm{bp}}^{(L)}\\\ \end{cases}$ (2)

Here, $L$ is the luminosity, $\gamma$ is the high-mass power law index,
$\delta$ is the rate of change of the index in the low-mass regime, and $A$ is
a constant, fixed by the requirement of the relationship to remain smooth
around the break point,

$A=\gamma-2\delta\log_{10}\left(M_{\mathrm{bp}}^{(L)}/M_{\odot}\right)$ (3)

The best-fit values of all six free parameters, as well as best-fit
normalization factors, are given in Table 1. Note that the standard ZAMS mass-
radius relationship for Population I ($\mathrm{PI}$) stars has
$\alpha_{\mathrm{PI}}\approx 0.8$ and $\beta_{\mathrm{PI}}\approx 0.57$
(Popper 1980, Kippenhahn & Weigert 1994, Ch. 22.1). While the high-mass index
matches the Population III value in Table 1, the low-mass index is discrepant
by a factor of $\sim 4$. This discrepancy arises because the energy production
mechanism transition in Population I stars ($\approx 1.3\ M_{\odot}$; Salaris
& Cassisi 2005, Ch. 5.2) very nearly coincides with the onset of convection in
the envelope ($\sim 1\ M_{\odot}$), while Population III stars maintain
radiative envelopes on either side of $M_{\mathrm{bp}}^{(R)}$. To compare
Population III and Population I mass-luminosity relationships, we consider the
average power law index in the $1\leq M/M_{\odot}\leq 10$ and $1\leq
M/M_{\odot}\leq 40$ ranges. The standard values for Population I are $3.88$
and $3.35$ respectively (Popper 1980, Kippenhahn & Weigert 1994, Ch. 22.1),
while the corresponding Population III values using Eq. 2 and Table 1 are
$3.87$ and $3.50$, suggesting that the Population III relationship is nearly
identical to its Population I counterpart at low masses and marginally steeper
at higher masses.

Figure 1: Mass-radius relationship from the metal-free evolutionary models
calculated in this study alongside the best-fit power law approximation given
in Eq. 1 and using the best-fit parameters in Table 1. Equivalent
relationships for non-zero metallicities were calculated as well and are shown
with dashed lines for comparison. Each dashed curve is parameterized by the
total metal mass fraction, $Z$. The vertical line indicates the break in the
best-fit power law relationship that originates from the onset of the CNO
cycle in the core. The point where the energy production rates of the proton-
proton chain ($L_{\mathrm{pp}}$) and the CNO cycle ($L_{\mathrm{CNO}}$) match
is indicated with blue squares for every mass-radius relationship shown.

The mass-radius relationship is plotted in Fig. 1 alongside the direct output
from selected MESA models. For comparison, we also calculated multiple grids
of MESA models at non-zero metallicities, whose mass-radius relationships are
overplotted in the figure as well. The expected broken power law behavior is
observed at all considered metallicities; however, the break point occurs at
progressively decreasing stellar masses with increasing metallicity due to the
larger initial carbon abundance. The point of equality in the energy
production rates of the proton-proton chain and the CNO cycle is also
indicated in the figure for every mass-radius relationship shown. Note that
while the equality point is strongly correlated with the power law break
point, the latter occurs at a lower stellar mass than the former as even a
subdominant contribution from the CNO cycle is sufficient to influence the
pressure structure within the star. In particular, we calculated the zero-
metallicity equality point as $\approx 20.5\ M_{\odot}$ – over $4\ M_{\odot}$
higher than $M_{\mathrm{bp}}^{(R)}$.

$1$ Parameter Value Error $\alpha$ $0.1982$ $\pm 0.0019$ $\beta$ $0.5527$ $\pm
0.0008$ $M_{\mathrm{bp}}^{(R)}$ $16.03$ $\pm 0.11$ $M_{\odot}$ $R$ at $1\
M_{\odot}$ $0.8792$ $\pm 0.0034$ $R_{\odot}$ $\delta$ $-0.6893$ $\pm 0.0030$
$\gamma$ $1.3137$ $\pm 0.0080$ $M_{\mathrm{bp}}^{(L)}$ $153.0$ $\pm 2.8$
$M_{\odot}$ $L$ at $1\ M_{\odot}$ $1.850$ $\pm 0.023$ $L_{\odot}$

Table 1: Analytic fit parameters

The effective temperature ($T_{\mathrm{eff}}$) and surface gravity
($\log_{10}(g)$) are related to $L$, $R$ and $M$ according to Eqs. 4 (Stefan-
Boltzmann law) and 5, where $\sigma$ and $G$ are the Stefan-Boltzmann and
gravitational constants respectively.

$T_{\mathrm{eff}}=\left(\frac{L}{4\pi\sigma R^{2}}\right)^{1/4}$ (4)

$\log_{10}(g)=\log_{10}\left(\frac{GM}{R^{2}\ [1\ \mathrm{cm}\
\mathrm{s}^{-2}]}\right)$ (5)

###  2.3 Atmosphere Modeling

Model atmospheres were calculated with version 9 of the ATLAS code (Kurucz,
1970; Sbordone et al., 2004; Castelli, 2005a; Kurucz, 2014). The code attains
high efficiency by sampling opacity from pre-tabulated opacity distribution
functions (ODFs), described in Kurucz et al. (1974) and Carbon (1984). The
ATLAS suite also contains the DFSYNTHE program (Castelli, 2005b) that may be
used to calculate ODFs for any given set of abundances, and the SYNTHE code
(Kurucz & Avrett, 1981) that computes the emergent spectrum from converged
ATLAS models by sampling the opacity directly at the wavelengths of interest.
All atmosphere models in this study were calculated at zero metallicity and
with a more precise estimate of the primordial helium mass fraction,
$Y=0.2448$, adopted from Valerdi et al. (2021).

For this project, we developed a universal Python dispatcher that combines all
three codes (originally written in Fortran) in a single user-friendly
pipeline, complete with intermediate consistency checks and comprehensive
documentation. Our dispatcher (Gerasimov & Larkin, 2022) is available
online444https://github.com/Roman-UCSD/BasicATLAS.

ATLAS stratifies the atmosphere into $72$ plane-parallel layers spanning the
range of Rosseland mean optical depths ($\tau$) from $\tau=10^{2}$ at the
bottom to $\tau=10^{-7}$ at the top. For all models with stellar masses below
$20\ M_{\odot}$, the ODFs were calculated following the “new” format (Castelli
& Kurucz, 2003) at $57$ temperatures between $10^{3.3}\ \mathrm{K}$ and
$10^{5.3}\ \mathrm{K}$. At higher masses, the range fails to accommodate the
deepest layers of the atmosphere that may exceed $10^{5.3}\ \mathrm{K}$ in
temperature. As such, a second set of ODFs was calculated with an extended
upper temperature limit of $10^{5.85}\ \mathrm{K}$. For those calculations,
the definition of the temperature grid stored in the variable TABT of the
LINOP() subroutine in the ATLAS source code was modified according to the
altered ODF format. Furthermore, the wavelength grid for opacity and radiation
field sampling in ATLAS, that by default spans from $\approx 9\ \mathrm{nm}$
to $160\text{\,}\mathrm{\SIUnitSymbolMicro m}$, had to be extended, first to
$4\ \mathrm{nm}$ at stellar masses over $6\ M_{\odot}$, and then to $0.1\
\mathrm{nm}$ at stellar masses over $29\ M_{\odot}$. These extensions avoid
errors in flux density and opacity integration due to significant
contributions outside the default wavelength range. For each model, the
adopted wavelength range was validated by ensuring that both the Planck
function ($B_{\nu}(T)$) and its derivative ($dB_{\nu}/dT$) drop below $0.1\%$
of their maximum values at the wavelength range bounds in each layer of the
atmosphere. The changes were applied to the WBIG variable in the BLOCKR()
subroutine of ATLAS.

Atmosphere calculations in ATLAS are carried out through iterative
improvements of an initial guess of the temperature profile throughout the
atmosphere until a new profile is found that meets both the hydrostatic
equilibrium condition for the prescribed surface gravity and the energy
equilibrium condition for the prescribed effective temperature. On each
iteration, the hydrostatic equilibrium condition is applied first to determine
the pressure profile corresponding to the current temperature profile. Energy
equilibrium is then evaluated throughout the atmosphere to determine
corrections for the current temperature profile as well as the current
percentage error in the flux and its derivative in each layer.

Since hydrostatic equilibrium is a hard requirement in ATLAS, no models can be
calculated for stars above the Eddington limit – a critical luminosity (or,
equivalently, critical effective temperature) above which the radiation
pressure gradient begins to exceed gravitational attraction. To determine this
limit, we considered the range of gravities between $\log_{10}(g)=4.6$ and
$\log_{10}(g)=5.2$ and searched for the maximum effective temperature
($T_{\mathrm{eff}}^{\mathrm{max}}$) at which ATLAS is able to find a solution
that is both in hydrostatic equilibrium and has flux and flux derivative
errors below $1\%$ over the course of $50$ iterations using the grey
temperature profile (Mihalas, 1978, Ch. 3),

$T(\tau)=T_{\mathrm{eff}}\left(\frac{3}{4}\tau+\frac{1}{2}\right)^{1/4}$ (6)

as the initial guess, where $T(\tau)$ is the temperature at optical depth
$\tau$. The calculated values of $T_{\mathrm{eff}}^{\mathrm{max}}$ showed a
nearly perfect exponential dependence on $\log_{10}(g)$ with deviations not
exceeding $1\%$ throughout the entire range of considered gravities. The
functional form of the relationship is given by

$\log_{10}(T_{\mathrm{eff}}^{\mathrm{max}})=C_{1}\log_{10}(g)+C_{2}$ (7)

where $C_{1}$ and $C_{2}$ are the best-fit parameters calculated as
$C_{1}=0.2516\pm 0.0013$ and $C_{2}=3.768\pm 0.006$ for CGS units. Eq. 7 is
directly comparable to the commonly adopted “classical” Eddington limit
(Mihalas 1978, Ch. 7-2, Rybicki & Lightman 1986, Ch. 1, Beznogov et al. 2020)
in CGS units, based on grey opacity dominated by Thomson scattering off
electrons and local thermodynamic equilibrium (LTE):

$\log_{10}(T_{\mathrm{eff}}^{\mathrm{max}})=\frac{1}{4}\log_{10}(g)+\frac{1}{4}\log_{10}\left(\frac{c\rho}{\sigma
n_{e}\sigma_{T}}\right)$ (8)

Here, $c$ is the speed of light, $\rho$ is the mass density of the atmospheric
layer, $n_{e}$ is the corresponding electron number density, and $\sigma_{T}$
is the Thomson scattering cross-section for an electron. When the layer is
fully ionized, $n_{e}/\rho$ only depends on the helium mass fraction:

$\frac{n_{e}}{\rho}=2\frac{Y}{m_{\mathrm{He}}}+\frac{1-Y}{m_{\mathrm{H}}}$ (9)

where $m_{\mathrm{He}}$ and $m_{\mathrm{H}}$ are the helium and hydrogen
atomic masses respectively. The “classical” equivalents of $C_{1}$ and
$C_{2}$, denoted $C_{1}^{(T)}$ and $C_{2}^{(T)}$ respectively, can be
evaluated numerically as $C_{1}^{(T)}=0.25$ and $C_{2}^{(T)}=3.795$. Since
$C_{1}^{(T)}\approx C_{1}$ and $C_{2}^{(T)}>C_{2}$, the exact Eddington limit
estimated with ATLAS is slightly lower than its “classical” counterpart at all
considered surface gravities due to additional non-grey opacity sources and
non-LTE effects in the atmosphere in the complete treatment.

Figure 2: Surface parameters of ZAMS Population III stars compared to the
exact calculation of the Eddington limit, including non-LTE effects and non-
grey atmospheric opacity. The dashed curve traces out an alternative locus of
surface parameters in the absence of the CNO break in the mass-radius
relationship at $M_{\mathrm{bp}}^{(R)}\approx 16\ M_{\odot}$. Surface
parameters extracted directly from selected MESA models are shown with red
markers and labelled by the initial masses (in $M_{\odot}$). The maximum
stellar mass with a convergent ATLAS atmosphere is shown with the blue circle
at $820.2\ M_{\odot}$. The “classical” (grey atmosphere, LTE) Eddington limit
from Eq. 8 is shown for reference. Note that the exact Eddington limit results
in a lower maximum mass value.

Eq. 7 is plotted in Fig. 2 alongside the locus of surface parameters of ZAMS
Population III stars derived from our analytic mass-radius and mass-luminosity
relationships in Eqs. 1 and 2. The intersection between the two curves
approximately represents the maximum mass of Population III stars with
atmospheres in hydrostatic equilibrium, which also serves as the maximum
initial stellar mass considered in this study. The maximum mass was calculated
as $820.2\ M_{\odot}$ by gradually increasing the initial stellar mass of the
model in increments of $0.1\ M_{\odot}$ until no convergent atmosphere
solution could be found. Figure 2 emphasizes the importance of the CNO cycle
in Population III stars, as the shape of the surface parameters locus is
clearly dominated by the power law break in the derived mass-radius
relationship at $M\approx 16\ M_{\odot}$.

Figure 3: Synthetic spectra of ZAMS Population III stars for initial stellar
masses between $1\ M_{\odot}$ and the Eddington limit ($820.2\ M_{\odot}$).
For clarity, each spectrum is shown on a separate color-coded vertical scale.
The blackbody spectra at the corresponding effective temperatures are shown in
dashes for reference. Important bound-free absorption breaks are highlighted
with vertical lines and labelled.

Overall, we calculated $59$ ATLAS atmospheres, logarithmically sampling the
range of initial stellar masses between $1\ M_{\odot}$ and $820.2\ M_{\odot}$
and using the derived analytic relationships for surface parameters. For each
model, the number of iterations was incremented in batches of $15$ until the
maximum flux error and the maximum flux derivative error dropped below the
standard convergence requirements of $1\%$ and $10\%$ respectively (Sbordone &
Bonifacio 2005; Mészáros et al. 2012; see Appendix A for details). Synthetic
spectra for each model were then calculated with SYNTHE between $0.5\
\mathrm{nm}$ and $2.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at the resolution
of $\lambda/\delta\lambda=6\times 10^{5}$. The chosen wavelength range ensures
that the flux density falls below $1\%$ of its maximum value at the range
bounds for all calculated model atmospheres. To account for the limited buffer
size in SYNTHE, the spectral synthesis for all atmospheres was carried out in
three batches: between $0.5\ \mathrm{nm}$ and $14\ \mathrm{nm}$; between $14\
\mathrm{nm}$ and $400\ \mathrm{nm}$; and between $400\ \mathrm{nm}$ and
$2.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The calculations were run in
parallel using the Triton Shared Computing Cluster (San Diego Supercomputer
Center, 2022). All calculated models are made public in our online
repository555https://atmos.ucsd.edu/. The key properties of all models as well
as their convergence parameters are tabulated in Appendix A.

Figure 4: The $5^{\mathrm{th}}$ and the $95^{\mathrm{th}}$ percentiles of the
wavelength distribution of the source function at a range of Rosseland mean
optical depths for the solar atmosphere (Bottom) and a near-Eddington
Population III atmosphere with $M=820.2\ M_{\odot}$ (Top). The cases of the
exact solution for the source function and the LTE approximation are shown in
solid red and dashed black respectively. The optical depth where the gas
temperature matches the effective temperature of the star is highlighted for
reference. The background color scheme corresponds to the total atmospheric
opacity as a function of both wavelength and Rosseland optical depth. The
discontinuous changes in opacity are due to bound-free absorption breaks.

A few representative synthetic spectra are shown in Fig. 3. The spectra of
stars with $M\gtrsim 50\ M_{\odot}$ display a considerable flux excess
blueward of the $\mathrm{He\ II}$ ionization break ($\approx 22.8\
\mathrm{nm}$) compared to their corresponding blackbody profiles. Since
ultraviolet emission is heavily attenuated by the interstellar medium in the
early universe, this effect results in an overall reduction of the observed
brightness of Population III stars. The blue excess is primarily caused by
non-LTE scattering of photons from deeper (and hotter) layers of the
atmosphere. The departures from LTE in the radiation field are shown in Fig. 4
for the highest-mass Population III model considered in this study ($820.2\
M_{\odot}$) as well as a solar atmosphere model ($T_{\mathrm{eff}}=5770\
\mathrm{K}$, $\log_{10}(g)=4.44$, abundances from Gerasimov et al. 2022). The
figure shows the $5^{\mathrm{th}}$ and the $95^{\mathrm{th}}$ percentiles of
the wavelength distribution of the source function ($S_{\lambda}$) for the
case of LTE ($S_{\lambda}=B_{\lambda}$, no scattering) and the complete
solution of the integral equation for $S_{\lambda}$ (Mihalas, 1978, Ch. 6-1).
For the solar model, both cases are nearly indistinguishable in all but the
outermost layers of the atmosphere that do not contribute to the emergent
spectrum significantly. On the other hand, the departure from LTE is far more
prominent in the Population III atmosphere, with a noticeable blue excess in
the radiation field at Rosseland mean optical depths shallower than $\sim 10$.

The line features in Fig. 3 diminish at higher masses due to the reduced
populations of neutral species in the atmosphere required for bound-bound
absorption. Selected lines may appear stronger than shown here due to
unaccounted higher-order NLTE effects (e.g. overpopulation of excited levels
as described in Auer & Mihalas 1972), as captured in model sets with a more
detailed treatment of NLTE line profiles (Bromm et al., 2001b; Rauch et al.,
2018; Rydberg et al., 2013). However, the impact of narrow line features on
broadband synthetic photometry is expected to be insignificant, especially in
the JWST bands chosen in this study (see Section 3) that mostly occupy the
comparatively line-free wavelength interval between the Lyman series of
hydrogen and the Fowler series of ionized helium. To verify this claim, we
recalculated our synthetic photometry (introduced in Section 3) with all line
profiles artificially strengthened by the extreme factor of $10$. We found
that at stellar masses over $100\ M_{\odot}$, the predicted magnitudes do not
deviate from their nominal values by more than $0.003\ \mathrm{mag}$ in the
chosen JWST bands.

##  3 Observable parameters

Figure 5: Color-magnitude diagrams of ZAMS Population III stars in the most
optimal JWST bands listed in Table 2 at four different redshifts. Synthetic
photometry is shown for both detailed opacity calculations based on model
atmospheres and the corresponding blackbody profiles. Selected initial stellar
masses are indicated in both cases with red and blue markers respectively. The
red markers are labelled in solar masses. The blue markers correspond to the
same masses as the red markers in the same order along the color-magnitude
curve. Observability limits for JWST are shown in grey and labelled with the
required gravitational lensing magnification, $\mu$. Figure 6: Mass-magnitude
relationships for ZAMS Population III stars in the most optimal JWST bands
listed in Table 2 at four different redshifts, color-coded. Note that a
different band is used at each redshift. The equivalent relationships for the
corresponding blackbody profiles are shown in dashed lines for comparison. The
observability limits for JWST are indicated with black markers for each
relationship shown. The displayed observability limits are grouped by the
required gravitational lensing magnitification, $\mu$.

Color-magnitude diagrams (CMDs) in Fig. 5 and mass-magnitude relationships in
Fig. 6 are provided for Population III stars in the most efficient JWST
transmission bands. All synthetic photometry in this study was carried out in
the ABMAG system (Oke & Gunn, 1983) from the newly calculated synthetic
spectra. For comparison, we also provide equivalent results for blackbody
atmospheres at the corresponding effective temperatures.

###  3.1 High-redshift synthetic photometry

In ABMAG, the apparent magnitude, $m_{\mathrm{AB}}$, is calculated from the
observed spectrum as:

$m_{AB}=-2.5\log_{10}\left(\frac{\int(\nu)^{-1}f_{\nu}(\nu)e(\nu)d\nu}{\int
3631\ \mathrm{Jy}\ (\nu)^{-1}e(\nu)d\nu}\right)$ (10)

Here, $f_{\nu}(\nu)$ is the apparent flux density per unit frequency, $e(\nu)$
is the efficiency of the instrument, and the integrals are evaluated over all
frequencies ($\nu$) in the frame of reference of the observer. The factor of
$(\nu)^{-1}$ is included to adapt the relationship to a photon-counting
instrument (Bohlin et al., 2014). Equivalently, Eq. 10 may be written in terms
of the observed wavelength, $\lambda$, to match the output of SYNTHE:

$m_{AB}=-2.5\log_{10}\left(\frac{\int\lambda
f_{\lambda}(\lambda)e(\lambda)d\lambda}{\int 3631\ \mathrm{Jy}\
(\lambda)^{-1}e(\lambda)cd\lambda}\right)$ (11)

At high redshift, the apparent flux density, $f_{\lambda}$, in Eq. 10 and Eq.
11 is derived from the modelled surface flux density, $F_{\lambda}$:

$f_{\lambda}(\lambda)=F_{\lambda}(\lambda_{\mathrm{e}})T(\lambda,z)\frac{R^{2}}{(1+z)D_{L}^{2}}$
(12)

(e.g. see Blanton & Roweis 2007), where $\lambda_{\mathrm{e}}=\lambda/(1+z)$
is the emitted wavelength at redshift $z$, $R$ is the radius of the star and
$D_{L}$ is the luminosity distance to the star. In the equation,
$T(\lambda,z)$ is the integrated transmissivity of the interstellar medium
across the line of sight to the source. In the wavelength range of interest,
the most significant contributions to interstellar attenuation are the bound-
free and bound-bound absorption by ground state neutral hydrogen at
$\lambda_{\mathrm{e}}\leq$1215.67\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ –
the Lyman $\alpha$ wavelength.

The photons absorbed by the interstellar medium will be re-emitted at longer
wavelengths. Depending on the dynamical evolution of the medium under
radiative feedback and the spatial density of Population III stars, this
reprocessed radiation may contribute significantly to the observed spectrum of
the star or have low impact on direct observations due to dilution over large
surface areas (e.g. compare “type A” and “type C” environments in Zackrisson
et al. 2011; also see Greif et al. 2009). In this study, we focus on the
purely stellar spectra; however, see Rydberg et al. (2013); Sibony et al.
(2022); Tumlinson et al. (2003); Kitayama et al. (2004) for various approaches
to detailed feedback modelling.

At $z\leq 7$, we adopt $T(\lambda,z)$ from the numerical simulation in Meiksin
(2006). At higher redshifts, the absorption by neutral hydrogen is
sufficiently strong to be well-approximated by a hard cut-off:

$T(\lambda,z;z>7)\approx\begin{cases}0,&\text{if }\lambda\leq(1+z)\
$1215.67\text{\,}\mathrm{\SIUnitSymbolAngstrom}$\\\ 1,&\mathrm{otherwise}\\\
\end{cases}$ (13)

As will be demonstrated, even at redshifts below $z=7$ considered in this
study, the most appropriate JWST bands for detecting Population III stars have
blue cut-offs at
$\lambda_{\mathrm{e}}\gtrsim$1215.67\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$,
thereby ensuring that Eq. 13 remains a good approximation at all considered
redshifts.

The luminosity distance, $D_{L}$, is calculated as a function of redshift as:

$D_{L}=(1+z)c\int_{0}^{z}{\frac{1}{H(z^{\prime})}dz^{\prime}}$ (14)

(Dodelson, 2003, Ch. 2.2) where $H(z^{\prime})$ is the Hubble parameter
(Dodelson, 2003, Ch. 2.4):

$H^{2}(z)=H_{0}^{2}\left(\Omega_{r}(1+z)^{4}+\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}\right)$
(15)

In Eq. 15, $H_{0}$ is the Hubble constant and $\Omega_{r}$, $\Omega_{M}$ and
$\Omega_{\Lambda}$ are the fractional present-day contributions of radiation
(including relativistic matter), non-relativistic matter, and dark energy,
respectively. We adopt the Hubble constant value of $H_{0}=67.4\ \mathrm{km}\
\mathrm{s}^{-1}\ \mathrm{Mpc}^{-1}$ and $\Omega_{M}=0.315$ in accordance with
Planck Collaboration et al. (2020). The universe is assumed to be flat
($\Omega_{\Lambda}=1-\Omega_{r}-\Omega_{M}$) and $\Omega_{r}$ is calculated as
(Dodelson, 2003, Ch. 2.4.4):

$\Omega_{r}=\left(1+\frac{7}{8}\left(\frac{4}{11}\right)^{\frac{4}{3}}N_{\mathrm{eff}}\right)\frac{4\sigma}{c^{3}}\frac{T^{4}_{\mathrm{CMB}}}{\rho_{\mathrm{crit}}}$
(16)

We adopt $N_{\mathrm{eff}}=3.04$ (Mangano et al., 2002) as the effective
number of neutrino flavors, $T_{\mathrm{CMB}}=2.725\ \mathrm{K}$ (Fixsen,
2009) as the present day temperature of the cosmic microwave background, and
$\rho_{\mathrm{crit}}=3/(8\pi G)H_{0}^{2}$ as the critical density of the
universe. The calculation implicitly assumes massless neutrinos. The accuracy
of these cosmological assumptions is examined in Section 3.4.

###  3.2 Choice of bands

This study considers detecting Population III stars with JWST at
$z\in(3,7,12,17)$. The chosen range spans between the redshift of the
predicted saturation of the Population III star formation rate in Xu et al.
(2013) and the redshift of the candidate Population III ionization source in
Fosbury et al. (2003). If the aim of the experiment is a simple detection of a
Population III candidate in a single band, the optimal observation band for
each stellar mass and redshift may be chosen by seeking the largest value of
the predicted signal-to-noise ratio. At the lowest stellar masses, this
condition will be met by a wide band, situated closest to the peak wavelength
of the model spectrum. Due to extensive attenuation of flux by the
interstellar medium at wavelengths shorter than the Lyman $\alpha$, the
optimal band at higher masses remains redward of the Lyman $\alpha$ in the
observer’s frame of reference instead of following the peak wavelength. The
transition occurs around $\sim 10\ M_{\odot}$ for $z=3$ and at $\lesssim 3\
M_{\odot}$ for $z\geq 7$. Since high-mass Population III stars are
overwhelmingly more likely to be observable, the same high-mass optimal band
may be safely employed for all Population III candidates at a given redshift.

We calculated the limiting magnitudes in all JWST Near Infrared Camera
(NIRCam) and Mid-Infrared Instrument (MIRI) bands using the JWST Exposure Time
Calculator (Pontoppidan et al., 2016) as the faintest magnitudes resulting in
a signal-to-noise ratio of $3$ in a $\sim 10\ \mathrm{hr}$ exposure (NIRCam:
20 groups, 9 integrations, DEEP2 readout pattern; MIRI: 100 groups, 132
integrations, FASTR1 readout pattern) for a flat frequency continuum. The best
band for each redshift was chosen as the one corresponding to the smallest
difference between the expected apparent magnitude of Population III stars in
the high-mass regime and the limiting magnitude of the band. The chosen bands
are listed in the “Best single” column of Table 2.

$2$ $z$ Best single Best pair $3$ NIRCam F070W NIRCam F070W NIRCam F090W $7$
NIRCam F150W2 NIRCam F115W NIRCam F150W $12$ NIRCam F200W NIRCam F200W NIRCam
F322W2 $17$ NIRCam F322W2 NIRCam F277W NIRCam F356W

Table 2: Optimal bands for Population III detection

$3$ Band Lim. mag Band Lim. mag NIRCam F070W $29.3664$ NIRCam F090W $29.6480$
NIRCam F150W2 $30.5796$ NIRCam F115W $29.8404$ NIRCam F150W $30.0816$ NIRCam
F200W $30.2616$ NIRCam F322W2 $30.5926$ NIRCam F277W $30.2774$ NIRCam F356W
$30.3625$

Table 3: Limiting magnitudes Table 4: Calculated limiting magnitudes are
listed for JWST bands in Table 2, assuming the detection signal-to-noise ratio
of $3$ in a $\sim 10\ \mathrm{hr}$ exposure (20 groups, 9 integrations, DEEP2
readout pattern) of a flat frequency continuum.

A more detailed experiment may be designed with the aim of measuring the
colors of Population III candidates in addition to simple detection, requiring
a choice of two filters without significant overlap in their transmission
profiles. We determine the optimal pairs of JWST filters for each redshift by
considering all possible non-overlapping pairs of bands and choosing the one
with the smallest average difference between the predicted magnitude of
Population III candidates and the limiting magnitude of the band. As before,
the choices are made in the high-mass regime. The resulting optimal pairs of
filters are listed in the “Best pair” column of Table 2.

NIRCam bands were found to be most optimal for both simple detection and color
measurements. The limiting magnitudes of all bands chosen in Table 2 are
listed in Table 4.

###  3.3 Results

Predicted color-magnitude diagrams are presented in Fig. 5 for all four
redshifts in the most optimal JWST band pairs listed in Table 2. Synthetic
photometry for blackbody atmospheres at the same range of effective
temperatures is shown in the figure as well for comparison. The overall trend
of the color-magnitude relationship is nearly unchanged between redshifts, as
the most efficient bands are placed in similar positions with respect to the
redshifted energy density distribution predicted by the model atmospheres. At
low masses ($M\lesssim 100\ M_{\odot}$), predicted colors shift blueward with
increasing effective temperature, with the exception of a brief inversion of
the trend around $6\ M_{\odot}$. We refer to this color-magnitude diagram
feature as the “helium loop”, as the inversion is caused by bound-free
absorption of singly ionized helium in the second excited state (“Fowler
break”, $\lambda_{\mathrm{e}}=205.1\ \mathrm{nm}$, Silberstein 1922). Once
formed, the break disproportionately suppresses flux in the bluer band,
resulting in the redder overall color.

In the high-mass regime ($M\gtrsim 100\ M_{\odot}$), the predicted color
shifts redward with increasing temperature due to progressively decreasing
contributions of free-free and bound-free opacities, both of which vary as
$\lambda^{3}$ and, therefore, redistribute the flux towards shorter
wavelengths. At high masses, Population III stars are predicted to be slightly
fainter than blackbodies with identical effective temperatures due to the non-
LTE distribution of the radiation field illustrated in Fig. 4.

Magnitude predictions for the most efficient single-band observations at each
redshift are shown in Fig. 6 as functions of mass. Both Figs. 5 and 6 also
contain the estimated JWST observability limits for different gravitational
lensing magnifications, from $\mu=10^{3}$ (approximate minimum required
magnification for direct observations of Population III stars) to $\mu=10^{7}$
(maximum theoretical magnification from Miralda-Escude 1991). The
observability limits are based on the calculated limiting magnitudes in each
band listed in Table 4.

###  3.4 Cosmological parameters

$4$ Adopted Value Range Effect $\mathrm{[mag]}$ Min Max Min Max
$H_{0}$aaHubble constant in $\mathrm{km}\ \mathrm{s}^{-1}\ \mathrm{Mpc}^{-1}$
$67.4$ $66.9$ $75.0$ $+0.02$ $-0.23$ $\Omega_{M}$ $0.315$ $0.300$ $0.322$
$+0.04$ $-0.02$ $m_{\nu}$bbNeutrino mass in $\mathrm{eV}\ c^{-2}$, assumed
identical for all neutrino species $0.0$ $0.0$ $0.9$ $0.0$ $-0.15$

Table 5: Effect of cosmological parameters

Predicted colors and magnitudes of Population III stars depend on the
cosmological parameters adopted when calculating the luminosity distance,
$D_{L}$ (Eq. 14). Since $D_{L}$ is independent of the physical properties of
the star, the effect is identical for all initial masses. In this section, we
estimate the magnitude of the effect for individual variations in the Hubble
constant, $H_{0}$; the present-day matter contribution, $\Omega_{M}$ and the
average neutrino mass, $m_{\nu}$. The nominal value of each parameter adopted
in this study, the considered range of variation, and the shift in the
predicted magnitudes at the lower and upper bounds of the considered range are
provided in Table 5. All tests are carried out at the largest considered
redshift, $z=17$, where the effect is expected to be most significant.

For the Hubble constant, the adopted range spans from the lower error bound of
the adopted nominal value ($67.4\pm 0.5\ \mathrm{km}\ \mathrm{s}^{-1}\
\mathrm{Mpc}^{-1}$, Planck Collaboration et al. 2020), based on Planck
observations of the cosmic microwave background (CMB), to the upper error
bound of the local Hubble constant estimate in Riess et al. (2016) ($73.24\pm
1.74\ \mathrm{km}\ \mathrm{s}^{-1}\ \mathrm{Mpc}^{-1}$), based on the updated
distance calibration to type Ia supernovae.

The discrepancy between the CMB and local measurements of $H_{0}$, known as
the Hubble Tension, depends on the adopted $\Omega_{M}$ prior. The two
measurements have been shown to be consistent at $95\%$ confidence for
$\Omega_{M}\lesssim 0.3$ (Wei & Melia, 2022). We therefore adopt
$\Omega_{M}=0.3$ as the lower bound on the variation range and take the error
in the nominal value ($0.315\pm 0.007$, Planck Collaboration et al. 2020) as
the upper bound.

To estimate the effect of massive neutrinos, we replace the relativistic
neutrino density in Eq. 16 with the approximation for neutrinos with identical
masses in Komatsu et al. (2011), implemented in Astropy (Astropy Collaboration
et al. 2013, 2018). We adopted $m_{\nu}=0.9\ \mathrm{eV}\ c^{-2}$ as the
maximum neutrino mass, measured in the Karlsruhe Tritium Neutrino (KATRIN)
experiment (Katrin Collaboration et al., 2022).

As demonstrated in Table 5, $H_{0}$ and $m_{\nu}$ have the largest effect on
the predicted photometry. Larger values for both parameters lead to shorter
lookback times to a given redshift and, therefore, brighter apparent
magnitudes of Population III stars. However, the gain in magnitude for both
parameters was calculated to fall below $0.25$ at the most extreme, which is
expected to remain within the measurement uncertainty at the adopted limiting
magnitude signal-to-noise ratio of $3$.

##  4 Conclusion

In this study, we calculated new evolutionary models and model atmospheres for
ZAMS Population III stars in hydrostatic equilibrium. The new models were used
to investigate the physical properties of the first stars in the universe as
well as to produce predictions of their colors and magnitudes as may be
observed in the near future with JWST under strong gravitational lensing. The
analysis was carried out at a broad range of plausible redshifts for
Population III stars from $z=3$ to $z=17$. Our predictions of Population III
color-magnitude diagrams and mass-magnitude relations are provided in Figs. 5
and 6 respectively. All predictions are given for the optimal JWST bands
listed in Table 2. Our other findings are listed below:

  * •

The mass-radius relationship of ZAMS Population III stars is well approximated
by a broken power law, similar to their Population I counterparts. However,
the break in the power law occurs at a much higher mass ($\approx 16\
M_{\odot}$ for Population III stars as opposed to $\approx 1\ M_{\odot}$ for
Population I) due to the suppressed CNO cycle.

  * •

Despite the initial absence of metals in Population III stars, the required
amount of carbon to sustain the CNO cycle is produced at sufficiently high
masses. The CNO cycle becomes the dominant energy production mechanism in
Population III stars around $M\approx 20.5\ M_{\odot}$, in agreement with Yoon
et al. (2012).

  * •

The mass-luminosity relationship of ZAMS Population III stars may be
approximated as a power law with a variable power index that decreases at
higher masses. This behavior is observed in Population I stars as well and is
approximately consistent with the Eddington standard model.

  * •

The evolution of true metal-free stars is nearly indistinguishable from the
evolution of extremely metal-poor stars with $Z\lesssim 10^{-9}$. This result
is more conservative but consistent with the $Z=10^{-8}$ limit derived in
Windhorst et al. (2018). Furthermore, both values agree with the lower bound
of the expected threshold of the Population III / Population II transition
($Z_{\mathrm{cr}}\gtrsim 10^{-8}$; Omukai et al. 2005; Schneider et al. 2006;
Clark et al. 2008).

  * •

The maximum mass of ZAMS Population III stars, at which hydrostatic
equilibrium is possible in the atmosphere (the Eddington limit) was calculated
as $\approx 820\ M_{\odot}$. This value is well above the commonly considered
range of initial masses for primordial stars, suggesting that hydrostatic
equilibrium may be an adequate approximation in Population III models. The
exact Eddington limit was found to be slightly lower than predicted by the
classical formula (Eq. 8) due to non-grey opacity sources in the atmosphere as
well as non-LTE effects. The influence of the CNO cycle on the internal
structure of Population III stars was determined to be a key factor in setting
the maximum mass.

  * •

Atmospheres of high-mass Population III stars host strongly non-LTE radiation
fields, resulting in significant excess in the UV flux compared to the
corresponding blackbody profiles.

  * •

The color-magnitude diagrams of Population III stars depend strongly on the
non-grey opacity sources in the atmosphere with notable features including the
“helium loop” at $M\sim 6\ M_{\odot}$ and the color-temperature inversion at
$M\gtrsim 100\ M_{\odot}$. In general, ZAMS Population III stars are fainter
than expected from blackbody profiles.

  * •

At the lowest redshift ($z=3$), the highest-mass Population III stars
considered in this study ($M\gtrsim 700\ M_{\odot}$) are just observable with
a gravitational lensing magnification of $\mu\sim 10^{3}$. A more plausible
range of stellar masses ($M\gtrsim 100\ M_{\odot}$) would likely require
$\mu\sim 10^{4}$. Such magnification is comparable to that inferred from
previous detections of the most distant individual stars known (e.g. Welch et
al. 2022a). At higher redshifts, the required magnification for an equivalent
detection increases to $\mu\sim 10^{5}$.

  * •

Our predictions of Population III observability do not depend significantly on
the adopted cosmological parameters; however, the maximum calculated effect of
$\sim 0.25\ \mathrm{mag}$ is comparable to the adopted signal-to-noise ratio
of JWST observations and may therefore be measurable under more generous
gravitational lensing conditions than the minimum detection requirement
considered in this study.

This study is limited to ZAMS Population III stars and may therefore be
considered a lower limit on true JWST observability since later evolution
stages are generally expected to be more luminous and less attenuated by the
interstellar medium. The reduction in surface gravity during the post-main
sequence evolution of Population III stars may drive the highest-mass stars
considered in this study above the calculated Eddington limit, requiring a
more detailed modelling approach allowing for mechanical motion in the
atmosphere as well as mass loss.

In our predictions of the observational signatures of Population III stars,
only radiation emitted directly from the stellar photospheres was considered.
Realistic regions of Population III formation are likely to display
significant flux contributions from the surrounding interstellar nebula. A
follow-up study could deploy an analytic ionization model, as in Sibony et al.
(2022), or a numerical simulation of radiative feedback, to derive the
necessary corrections. In this context, the predictions drawn here may once
again be interpreted as lower limits of the true observability of individual
Population III sources.

While the detailed modelling of atmospheric opacity was shown to produce
noticeably different results from the commonly adopted blackbody approximation
(e.g. Windhorst et al. 2018; Fosbury et al. 2003), the quantitative difference
in predictions of the two approaches will likely remain within the measurement
uncertainty for Population III candidates at the observability threshold (e.g.
at $z=3$, the magnitude difference between the two approaches is $\approx 0.2\
\mathrm{mag}$ at the highest considered mass). However, the discrepancy may be
detectable under marginally stronger gravitational lensing and should
therefore be taken into account in more detailed observational studies of
Population III stars.

Our overall result generally agrees with previous studies of Population III
observability (e.g. Windhorst et al. 2018; Rydberg et al. 2013) that detection
of the first stars in the universe may be possible with JWST under strong but
realistic gravitational lensing, assuming sufficiently high stellar mass.
Placing more specific constraints on the expected rate of detection remains
challenging due to the highly debated initial mass function of Population III
stars.

Finally, we note that Population III stars likely formed in clusters or
galaxies rather than in isolation (e.g. Jaacks et al. 2019; Zackrisson et al.
2011; Johnson 2010; Johnson et al. 2009; DeSouza & Basu 2015; Visbal et al.
2017). The combined luminosity of such objects makes them more accessible
targets, requiring less extreme gravitational lensing. However, modelling the
spectral energy distributions of Population III clusters is further
complicated by the dependency on the highly uncertain initial mass function
and the concurrent formation of Population III and Population II stars (Xu et
al., 2013; Wise et al., 2012b, a; Jaacks et al., 2019) at later epochs.

M. Larkin acknowledges funding support from the University of California at
San Diego (UCSD) Department of Physical Sciences Summer Research Award and the
UCSD Triton Research & Experiential Learning Scholars (TRELS) program. R.
Gerasimov acknowledges funding support from HST Program GO-15096, provided by
NASA through a grant from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy,
Incorporated, under NASA contract NAS5-26555. This work was conducted at UCSD,
which was built on the unceded territory of the Kumeyaay Nation. Today, the
Kumeyaay people continue to maintain their political sovereignty and cultural
traditions as vital members of the San Diego community. We acknowledge their
tremendous contributions to our region and thank them for their stewardship.

##  Appendix A Model parameters

Table 6 lists the defining parameters and convergence criteria of all ATLAS-9
model atmospheres calculated in this study. Initial masses are sampled
logarithmically between $1\ M_{\odot}$ and the estimated Eddington limit
($820\ M_{\odot}$). The corresponding stellar radii, luminosities, effective
temperatures and surface gravities are calculated using Eqs. 1, 2, 4 and 5 as
well as the best-fit parameters in Table 1. Only the effective temperature and
surface gravity are used as inputs to model atmospheres. The final convergence
is parameterized in terms of the maximum absolute flux error and the maximum
absolute flux derivative error with respect to the depth-integrated mass
density (see Kurucz 1970). Flux errors and flux derivative errors are used to
calculate the temperature corrections between iterations using the Avrett &
Krook (1963) scheme at large optical depths and the $\Lambda$-iteration scheme
(Böhm-Vitense 1964, Mihalas 1978, Ch. 3-3) at shallow optical depths,
respectively. All models calculated in this study meet the standard
convergence target of flux error below $1\%$ and flux derivative error below
$10\%$ (Sbordone & Bonifacio, 2005; Mészáros et al., 2012).

$5$ Initial Mass Effective Temperature Surface gravity Radius Luminosity Max
Flux Max Flux $M\ [M_{\odot}]$ $T_{\mathrm{eff}}\ [\mathrm{K}]$ $\log_{10}(g)$
$R\ [R_{\odot}]$ $\log_{10}(L/L_{\odot})$ Error $[\%]$ Derivative Error $[\%]$
$1.000$ $7180$ $4.550$ $0.879$ $0.267$ $0.24$ $6.58$ $1.124$ $8047$ $4.581$
$0.900$ $0.485$ $0.62$ $1.54$ $1.264$ $9000$ $4.612$ $0.921$ $0.700$ $0.61$
$3.11$ $1.421$ $10045$ $4.642$ $0.942$ $0.911$ $0.71$ $7.49$ $1.597$ $11189$
$4.673$ $0.964$ $1.118$ $0.76$ $8.30$ $1.796$ $12437$ $4.704$ $0.987$ $1.322$
$0.78$ $7.88$ $2.019$ $13796$ $4.734$ $1.010$ $1.523$ $0.85$ $6.84$ $2.270$
$15273$ $4.765$ $1.034$ $1.719$ $0.79$ $5.00$ $2.551$ $16873$ $4.796$ $1.058$
$1.912$ $0.84$ $3.89$ $2.868$ $18602$ $4.827$ $1.083$ $2.102$ $0.74$ $4.94$
$3.225$ $20466$ $4.857$ $1.108$ $2.288$ $0.58$ $5.76$ $3.625$ $22472$ $4.888$
$1.134$ $2.471$ $0.53$ $5.73$ $4.075$ $24623$ $4.919$ $1.161$ $2.650$ $0.45$
$5.17$ $4.582$ $26924$ $4.949$ $1.188$ $2.825$ $0.35$ $4.50$ $5.151$ $29381$
$4.980$ $1.216$ $2.997$ $0.34$ $3.23$ $5.790$ $31996$ $5.011$ $1.245$ $3.165$
$0.31$ $1.63$ $6.510$ $34772$ $5.042$ $1.274$ $3.330$ $0.23$ $1.31$ $7.318$
$37712$ $5.072$ $1.304$ $3.491$ $0.24$ $0.93$ $8.227$ $40817$ $5.103$ $1.334$
$3.648$ $0.28$ $0.60$ $9.249$ $44087$ $5.134$ $1.365$ $3.802$ $0.31$ $0.65$
$10.398$ $47521$ $5.164$ $1.397$ $3.953$ $0.32$ $0.47$ $11.690$ $51118$
$5.195$ $1.430$ $4.100$ $0.32$ $0.36$ $13.141$ $54874$ $5.226$ $1.464$ $4.243$
$0.27$ $0.27$ $14.774$ $58785$ $5.256$ $1.498$ $4.383$ $0.24$ $0.23$ $16.609$
$62432$ $5.276$ $1.554$ $4.519$ $0.20$ $0.26$ $18.672$ $65238$ $5.270$ $1.658$
$4.652$ $0.18$ $0.31$ $20.991$ $68031$ $5.265$ $1.769$ $4.781$ $0.17$ $0.46$
$23.598$ $70799$ $5.260$ $1.887$ $4.906$ $0.15$ $0.49$ $26.529$ $73528$
$5.254$ $2.013$ $5.028$ $0.12$ $0.56$ $29.825$ $76205$ $5.249$ $2.148$ $5.147$
$0.14$ $0.59$ $33.529$ $78819$ $5.243$ $2.291$ $5.261$ $0.15$ $0.70$ $37.694$
$81355$ $5.238$ $2.445$ $5.373$ $0.21$ $0.74$ $42.376$ $83800$ $5.233$ $2.608$
$5.480$ $0.26$ $0.75$ $47.639$ $86143$ $5.227$ $2.783$ $5.584$ $0.20$ $0.84$
$53.557$ $88369$ $5.222$ $2.969$ $5.685$ $0.22$ $1.09$ $60.209$ $90467$
$5.216$ $3.167$ $5.782$ $0.17$ $1.15$ $67.688$ $92425$ $5.211$ $3.379$ $5.875$
$0.21$ $1.08$ $76.095$ $94232$ $5.206$ $3.605$ $5.965$ $0.22$ $0.93$ $85.547$
$95878$ $5.200$ $3.846$ $6.052$ $0.23$ $0.69$ $96.172$ $97352$ $5.195$ $4.103$
$6.134$ $0.23$ $0.64$ $108.118$ $98647$ $5.189$ $4.378$ $6.214$ $0.23$ $0.65$
$121.547$ $99754$ $5.184$ $4.671$ $6.289$ $0.26$ $0.53$ $136.645$ $100666$
$5.179$ $4.983$ $6.361$ $0.31$ $0.45$ $153.617$ $101379$ $5.173$ $5.317$
$6.430$ $0.37$ $0.51$ $172.698$ $102000$ $5.168$ $5.672$ $6.497$ $0.44$ $0.59$
$194.149$ $102624$ $5.162$ $6.052$ $6.563$ $0.47$ $0.60$ $218.264$ $103253$
$5.157$ $6.456$ $6.630$ $0.40$ $0.49$ $245.375$ $103885$ $5.152$ $6.888$
$6.697$ $0.46$ $0.54$ $275.853$ $104521$ $5.146$ $7.349$ $6.764$ $0.32$ $0.38$
$310.117$ $105160$ $5.141$ $7.840$ $6.831$ $0.19$ $0.26$ $348.637$ $105804$
$5.136$ $8.365$ $6.898$ $0.13$ $0.26$ $391.941$ $106452$ $5.130$ $8.924$
$6.964$ $0.15$ $0.39$ $440.624$ $107103$ $5.125$ $9.521$ $7.031$ $0.18$ $0.43$
$495.354$ $107759$ $5.119$ $10.158$ $7.098$ $0.19$ $0.41$ $556.881$ $108419$
$5.114$ $10.838$ $7.165$ $0.22$ $0.44$ $626.052$ $109082$ $5.109$ $11.562$
$7.232$ $0.25$ $0.48$ $703.814$ $109750$ $5.103$ $12.336$ $7.299$ $0.30$
$6.53$ $791.234$ $110422$ $5.098$ $13.161$ $7.365$ $0.37$ $7.67$ $820.200$
$110629$ $5.096$ $13.425$ $7.386$ $0.60$ $4.59$

Table 6: Parameters of Population III models calculated in this study

## References

  * Abe et al. (2021) Abe, M., Yajima, H., Khochfar, S., Dalla Vecchia, C., & Omukai, K. 2021, MNRAS, 508, 3226, doi: 10.1093/mnras/stab2637
  * Abel et al. (2000) Abel, T., Bryan, G. L., & Norman, M. L. 2000, ApJ, 540, 39, doi: 10.1086/309295
  * Abohalima & Frebel (2018) Abohalima, A., & Frebel, A. 2018, ApJS, 238, 36, doi: 10.3847/1538-4365/aadfe9
  * Antonov & Chernin (1977) Antonov, V. A., & Chernin, A. D. 1977, Astrofizika, 13, 271
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Auer & Mihalas (1972) Auer, L. H., & Mihalas, D. 1972, ApJS, 24, 193, doi: 10.1086/190253
  * Avrett & Krook (1963) Avrett, E. H., & Krook, M. 1963, ApJ, 137, 874, doi: 10.1086/147564
  * Baade (1944) Baade, W. 1944, ApJ, 100, 137, doi: 10.1086/144650
  * Bahcall et al. (1994) Bahcall, J. N., Flynn, C., Gould, A., & Kirhakos, S. 1994, ApJ, 435, L51, doi: 10.1086/187592
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Beznogov et al. (2020) Beznogov, M. V., Page, D., & Ramirez-Ruiz, E. 2020, ApJ, 888, 97, doi: 10.3847/1538-4357/ab5fd6
  * Blanton & Roweis (2007) Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734, doi: 10.1086/510127
  * Bohlin et al. (2014) Bohlin, R. C., Gordon, K. D., & Tremblay, P. E. 2014, PASP, 126, 711, doi: 10.1086/677655
  * Böhm-Vitense (1964) Böhm-Vitense, E. 1964, SAO Special Report, 167, 99
  * Bond (1980) Bond, H. E. 1980, in Bulletin of the American Astronomical Society, Vol. 12, 831
  * Bromm et al. (2002) Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23, doi: 10.1086/323947
  * Bromm et al. (2001a) Bromm, V., Ferrara, A., Coppi, P. S., & Larson, R. B. 2001a, MNRAS, 328, 969, doi: 10.1046/j.1365-8711.2001.04915.x
  * Bromm et al. (2001b) Bromm, V., Kudritzki, R. P., & Loeb, A. 2001b, ApJ, 552, 464, doi: 10.1086/320549
  * Carbon (1984) Carbon, D. F. 1984, Line blanketing. (Cambridge University Press), 395–426
  * Castelli (2005a) Castelli, F. 2005a, Memorie della Societa Astronomica Italiana Supplementi, 8, 25
  * Castelli (2005b) —. 2005b, Memorie della Societa Astronomica Italiana Supplementi, 8, 34
  * Castelli & Kurucz (2003) Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20.  https://arxiv.org/abs/astro-ph/0405087
  * Cen (2003a) Cen, R. 2003a, ApJ, 591, 12, doi: 10.1086/375217
  * Cen (2003b) —. 2003b, ApJ, 591, L5, doi: 10.1086/377068
  * Chen et al. (2022) Chen, W., Kelly, P. L., Treu, T., et al. 2022, arXiv e-prints, arXiv:2207.11658.  https://arxiv.org/abs/2207.11658
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, The Astrophysical Journal, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Clark et al. (2008) Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2008, ApJ, 672, 757, doi: 10.1086/524187
  * Clark et al. (2011) Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011, ApJ, 727, 110, doi: 10.1088/0004-637X/727/2/110
  * Cojazzi et al. (2000) Cojazzi, P., Bressan, A., Lucchin, F., Pantano, O., & Chavez, M. 2000, MNRAS, 315, L51, doi: 10.1046/j.1365-8711.2000.03670.x
  * Crighton et al. (2016) Crighton, N. H. M., O’Meara, J. M., & Murphy, M. T. 2016, MNRAS, 457, L44, doi: 10.1093/mnrasl/slv191
  * Demircan & Kahraman (1991) Demircan, O., & Kahraman, G. 1991, Ap&SS, 181, 313, doi: 10.1007/BF00639097
  * DeSouza & Basu (2015) DeSouza, A. L., & Basu, S. 2015, MNRAS, 450, 295, doi: 10.1093/mnras/stv523
  * Diego (2019) Diego, J. M. 2019, A&A, 625, A84, doi: 10.1051/0004-6361/201833670
  * Diego et al. (2018) Diego, J. M., Kaiser, N., Broadhurst, T., et al. 2018, ApJ, 857, 25, doi: 10.3847/1538-4357/aab617
  * Dodelson (2003) Dodelson, S. 2003, Modern Cosmology (Academic Press)
  * Eddington (1918) Eddington, A. S. 1918, ApJ, 48, 205, doi: 10.1086/142427
  * Fixsen (2009) Fixsen, D. J. 2009, ApJ, 707, 916, doi: 10.1088/0004-637X/707/2/916
  * Fosbury et al. (2003) Fosbury, R. A. E., Villar-Martín, M., Humphrey, A., et al. 2003, ApJ, 596, 797, doi: 10.1086/378228
  * Freese (2017) Freese, K. 2017, International Journal of Modern Physics D, 26, 1730012, doi: 10.1142/S0218271817300129
  * Fuller et al. (1986) Fuller, G. M., Woosley, S. E., & Weaver, T. A. 1986, ApJ, 307, 675, doi: 10.1086/164452
  * Gerasimov et al. (2022) Gerasimov, R., Burgasser, A. J., Homeier, D., et al. 2022, ApJ, 930, 24, doi: 10.3847/1538-4357/ac61e5
  * Gerasimov & Larkin (2022) Gerasimov, R., & Larkin, M. 2022, Roman-UCSD/BasicATLAS: v1.0.0, v1.0.0, Zenodo, doi: 10.5281/zenodo.7145514
  * Gimenez & Zamorano (1985) Gimenez, A., & Zamorano, J. 1985, Ap&SS, 114, 259, doi: 10.1007/BF00653969
  * Graff & Freese (1996) Graff, D. S., & Freese, K. 1996, ApJ, 467, L65, doi: 10.1086/310195
  * Greif et al. (2010) Greif, T. H., Glover, S. C. O., Bromm, V., & Klessen, R. S. 2010, ApJ, 716, 510, doi: 10.1088/0004-637X/716/1/510
  * Greif et al. (2009) Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2009, MNRAS, 399, 639, doi: 10.1111/j.1365-2966.2009.15336.x
  * Greif et al. (2011) Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75, doi: 10.1088/0004-637X/737/2/75
  * Haiman & Loeb (1997) Haiman, Z., & Loeb, A. 1997, ApJ, 483, 21, doi: 10.1086/304238
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Heger & Woosley (2002) Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532, doi: 10.1086/338487
  * Herrington et al. (2022) Herrington, N. P., Whalen, D. J., & Woods, T. E. 2022, arXiv e-prints, arXiv:2208.00008.  https://arxiv.org/abs/2208.00008
  * Hirano et al. (2015) Hirano, S., Hosokawa, T., Yoshida, N., Omukai, K., & Yorke, H. W. 2015, MNRAS, 448, 568, doi: 10.1093/mnras/stv044
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60, doi: 10.1088/0004-637X/781/2/60
  * Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250, doi: 10.1126/science.1207433
  * Hubeny & Lanz (1995) Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875, doi: 10.1086/175226
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jaacks et al. (2019) Jaacks, J., Finkelstein, S. L., & Bromm, V. 2019, MNRAS, 488, 2202, doi: 10.1093/mnras/stz1529
  * Johnson (2010) Johnson, J. L. 2010, MNRAS, 404, 1425, doi: 10.1111/j.1365-2966.2010.16351.x
  * Johnson (2015) —. 2015, MNRAS, 453, 2771, doi: 10.1093/mnras/stv1815
  * Johnson et al. (2009) Johnson, J. L., Greif, T. H., Bromm, V., Klessen, R. S., & Ippolito, J. 2009, MNRAS, 399, 37, doi: 10.1111/j.1365-2966.2009.15158.x
  * Kashlinsky (2005) Kashlinsky, A. 2005, ApJ, 633, L5, doi: 10.1086/498243
  * Kashlinsky et al. (2005) Kashlinsky, A., Arendt, R. G., Mather, J., & Moseley, S. H. 2005, in American Astronomical Society Meeting Abstracts, Vol. 207, American Astronomical Society Meeting Abstracts, 100.02
  * Katrin Collaboration et al. (2022) Katrin Collaboration, Aker, M., Beglarian, A., et al. 2022, Nature Physics, 18, 160, doi: 10.1038/s41567-021-01463-1
  * Keller et al. (2014) Keller, S. C., Bessell, M. S., Frebel, A., et al. 2014, Nature, 506, 463, doi: 10.1038/nature12990
  * Kelly et al. (2018) Kelly, P. L., Diego, J. M., Rodney, S., et al. 2018, Nature Astronomy, 2, 334, doi: 10.1038/s41550-018-0430-3
  * Kinugawa et al. (2014) Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963, doi: 10.1093/mnras/stu1022
  * Kippenhahn & Weigert (1994) Kippenhahn, R., & Weigert, A. 1994, Stellar Structure and Evolution (Springer Berlin / Heidelberg)
  * Kitayama et al. (2004) Kitayama, T., Yoshida, N., Susa, H., & Umemura, M. 2004, ApJ, 613, 631, doi: 10.1086/423313
  * Kogut (2003) Kogut, A. 2003, New A Rev., 47, 977, doi: 10.1016/j.newar.2003.09.029
  * Köhler et al. (2015) Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71, doi: 10.1051/0004-6361/201424356
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18
  * Kowalska et al. (2012) Kowalska, I., Bulik, T., & Belczynski, K. 2012, A&A, 541, A120, doi: 10.1051/0004-6361/201118604
  * Kurucz (1970) Kurucz, R. L. 1970, SAO Special Report, 309
  * Kurucz (2014) —. 2014, in Determination of Atmospheric Parameters of B-, A-, F- and G-Type Stars (Springer International Publishing), 39–51, doi: 10.1007/978-3-319-06956-2_4
  * Kurucz & Avrett (1981) Kurucz, R. L., & Avrett, E. H. 1981, SAO Special Report, 391
  * Kurucz et al. (1974) Kurucz, R. L., Peytremann, E., & Avrett, E. H. 1974, Blanketed model atmospheres for early-type stars (Smithsonian Astrophysical Observatory)
  * Lacy (1979) Lacy, C. H. 1979, ApJ, 228, 817, doi: 10.1086/156908
  * Larson (1998) Larson, R. B. 1998, MNRAS, 301, 569, doi: 10.1046/j.1365-8711.1998.02045.x
  * Liu & Bromm (2020) Liu, B., & Bromm, V. 2020, MNRAS, 497, 2839, doi: 10.1093/mnras/staa2143
  * Liu et al. (2021) Liu, B., Sibony, Y., Meynet, G., & Bromm, V. 2021, MNRAS, 506, 5247, doi: 10.1093/mnras/stab2057
  * Madau & Silk (2005) Madau, P., & Silk, J. 2005, MNRAS, 359, L37, doi: 10.1111/j.1745-3933.2005.00031.x
  * Maio et al. (2010) Maio, U., Ciardi, B., Dolag, K., Tornatore, L., & Khochfar, S. 2010, MNRAS, 407, 1003, doi: 10.1111/j.1365-2966.2010.17003.x
  * Maio et al. (2016) Maio, U., Petkova, M., De Lucia, G., & Borgani, S. 2016, MNRAS, 460, 3733, doi: 10.1093/mnras/stw1196
  * Mangano et al. (2002) Mangano, G., Miele, G., Pastor, S., & Peloso, M. 2002, Physics Letters B, 534, 8, doi: 10.1016/S0370-2693(02)01622-2
  * Meiksin (2006) Meiksin, A. 2006, MNRAS, 365, 807, doi: 10.1111/j.1365-2966.2005.09756.x
  * Mészáros et al. (2012) Mészáros, S., Allende Prieto, C., Edvardsson, B., et al. 2012, AJ, 144, 120, doi: 10.1088/0004-6256/144/4/120
  * Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres (W. H. Freeman)
  * Miralda-Escude (1991) Miralda-Escude, J. 1991, ApJ, 379, 94, doi: 10.1086/170486
  * Muñoz et al. (2021) Muñoz, V., Takhistov, V., Witte, S. J., & Fuller, G. M. 2021, J. Cosmology Astropart. Phys, 2021, 020, doi: 10.1088/1475-7516/2021/11/020
  * Muratov et al. (2013) Muratov, A. L., Gnedin, O. Y., Gnedin, N. Y., & Zemp, M. 2013, ApJ, 773, 19, doi: 10.1088/0004-637X/773/1/19
  * Narayan & Wallington (1992) Narayan, R., & Wallington, S. 1992, in Gravitational Lenses, ed. R. Kayser, T. Schramm, & L. Nieser, Vol. 406 (Springer, Berlin, Heidelberg), 12, doi: 10.1007/3-540-55797-0_76
  * Nordlander et al. (2019) Nordlander, T., Bessell, M. S., Da Costa, G. S., et al. 2019, MNRAS, 488, L109, doi: 10.1093/mnrasl/slz109
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627, doi: 10.1086/429955
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, 991016, doi: 10.1117/12.2231768
  * Popper (1980) Popper, D. M. 1980, ARA&A, 18, 115, doi: 10.1146/annurev.aa.18.090180.000555
  * Rauch et al. (2018) Rauch, T., Demleitner, M., Hoyer, D., & Werner, K. 2018, MNRAS, 475, 3896, doi: 10.1093/mnras/sty056
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Rowan-Robinson (1983) Rowan-Robinson, M. 1983, in Early Evolution of the Universe and its Present Structure, ed. G. O. Abell & G. Chincarini, Vol. 104, 109–112
  * Rybicki & Lightman (1986) Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics (Wiley-VCH)
  * Rydberg et al. (2013) Rydberg, C.-E., Zackrisson, E., Lundqvist, P., & Scott, P. 2013, MNRAS, 429, 3658, doi: 10.1093/mnras/sts653
  * Salaris & Cassisi (2005) Salaris, M., & Cassisi, S. 2005, Evolution of Stars and Stellar Populations (John Wiley and Sons)
  * San Diego Supercomputer Center (2022) San Diego Supercomputer Center. 2022, doi: 10.57873/T34W2R
  * Santos et al. (2002) Santos, M. R., Bromm, V., & Kamionkowski, M. 2002, MNRAS, 336, 1082, doi: 10.1046/j.1365-8711.2002.05895.x
  * Sbordone & Bonifacio (2005) Sbordone, L., & Bonifacio, P. 2005, Atlas Cookbook.  http://atmos.obspm.fr/index.php/documentation/7
  * Sbordone et al. (2004) Sbordone, L., Bonifacio, P., Castelli, F., & Kurucz, R. L. 2004, Memorie della Societa Astronomica Italiana Supplementi, 5, 93.  https://arxiv.org/abs/astro-ph/0406268
  * Schaerer (2002) Schaerer, D. 2002, A&A, 382, 28, doi: 10.1051/0004-6361:20011619
  * Schauer et al. (2022) Schauer, A. T. P., Bromm, V., Drory, N., & Boylan-Kolchin, M. 2022, arXiv e-prints, arXiv:2207.02863.  https://arxiv.org/abs/2207.02863
  * Schneider et al. (2002) Schneider, R., Ferrara, A., Natarajan, P., & Omukai, K. 2002, ApJ, 571, 30, doi: 10.1086/339917
  * Schneider et al. (2006) Schneider, R., Omukai, K., Inoue, A. K., & Ferrara, A. 2006, MNRAS, 369, 1437, doi: 10.1111/j.1365-2966.2006.10391.x
  * Schwarzschild & Spitzer (1953) Schwarzschild, M., & Spitzer, L. 1953, The Observatory, 73, 77
  * Sharda & Krumholz (2022) Sharda, P., & Krumholz, M. R. 2022, MNRAS, 509, 1959, doi: 10.1093/mnras/stab2921
  * Sibony et al. (2022) Sibony, Y., Liu, B., Simmonds, C., Meynet, G., & Bromm, V. 2022, arXiv e-prints, arXiv:2205.15125.  https://arxiv.org/abs/2205.15125
  * Silberstein (1922) Silberstein, L. 1922, Nature, 110, 247, doi: 10.1038/110247b0
  * Sokasian et al. (2004) Sokasian, A., Yoshida, N., Abel, T., Hernquist, L., & Springel, V. 2004, MNRAS, 350, 47, doi: 10.1111/j.1365-2966.2004.07636.x
  * Spergel et al. (2003) Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175, doi: 10.1086/377226
  * Stacy & Bromm (2014) Stacy, A., & Bromm, V. 2014, ApJ, 785, 73, doi: 10.1088/0004-637X/785/1/73
  * Stacy et al. (2016) Stacy, A., Bromm, V., & Lee, A. T. 2016, MNRAS, 462, 1307, doi: 10.1093/mnras/stw1728
  * Stacy et al. (2012) Stacy, A., Greif, T. H., & Bromm, V. 2012, MNRAS, 422, 290, doi: 10.1111/j.1365-2966.2012.20605.x
  * Susa et al. (2014) Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32, doi: 10.1088/0004-637X/792/1/32
  * Suwa et al. (2007) Suwa, Y., Takiwaki, T., Kotake, K., & Sato, K. 2007, ApJ, 665, L43, doi: 10.1086/521078
  * Trenti et al. (2009) Trenti, M., Stiavelli, M., & Shull, J. M. 2009, ApJ, 700, 1672, doi: 10.1088/0004-637X/700/2/1672
  * Tumlinson et al. (2003) Tumlinson, J., Shull, J. M., & Venkatesan, A. 2003, ApJ, 584, 608, doi: 10.1086/345737
  * Umeda & Nomoto (2002) Umeda, H., & Nomoto, K. 2002, ApJ, 565, 385, doi: 10.1086/323946
  * Valerdi et al. (2021) Valerdi, M., Peimbert, A., & Peimbert, M. 2021, MNRAS, 505, 3624, doi: 10.1093/mnras/stab1543
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Visbal et al. (2017) Visbal, E., Bryan, G. L., & Haiman, Z. 2017, MNRAS, 469, 1456, doi: 10.1093/mnras/stx909
  * Wagner (1978) Wagner, R. L. 1978, A&A, 62, 9
  * Wei & Melia (2022) Wei, J.-J., & Melia, F. 2022, ApJ, 928, 165, doi: 10.3847/1538-4357/ac562c
  * Welch et al. (2022a) Welch, B., Coe, D., Diego, J. M., et al. 2022a, Nature, 603, 815, doi: 10.1038/s41586-022-04449-y
  * Welch et al. (2022b) Welch, B., Coe, D., Zackrisson, E., et al. 2022b, arXiv e-prints, arXiv:2208.09007.  https://arxiv.org/abs/2208.09007
  * Windhorst et al. (2018) Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ApJS, 234, 41, doi: 10.3847/1538-4365/aaa760
  * Wise et al. (2012a) Wise, J. H., Abel, T., Turk, M. J., Norman, M. L., & Smith, B. D. 2012a, MNRAS, 427, 311, doi: 10.1111/j.1365-2966.2012.21809.x
  * Wise et al. (2012b) Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012b, ApJ, 745, 50, doi: 10.1088/0004-637X/745/1/50
  * Woods et al. (2021) Woods, T. E., Patrick, S., Elford, J. S., Whalen, D. J., & Heger, A. 2021, ApJ, 915, 110, doi: 10.3847/1538-4357/abfaf9
  * Wyithe & Loeb (2003) Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 586, 693, doi: 10.1086/367721
  * Xu et al. (2013) Xu, H., Wise, J. H., & Norman, M. L. 2013, ApJ, 773, 83, doi: 10.1088/0004-637X/773/2/83
  * Yoon et al. (2012) Yoon, S. C., Dierks, A., & Langer, N. 2012, A&A, 542, A113, doi: 10.1051/0004-6361/201117769
  * Yoshida et al. (2006) Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6, doi: 10.1086/507978
  * Yung et al. (2020) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2020, MNRAS, 496, 4574, doi: 10.1093/mnras/staa1800
  * Zackrisson et al. (2015) Zackrisson, E., González, J., Eriksson, S., et al. 2015, MNRAS, 449, 3057, doi: 10.1093/mnras/stv492
  * Zackrisson et al. (2011) Zackrisson, E., Rydberg, C.-E., Schaerer, D., Östlin, G., & Tuli, M. 2011, ApJ, 740, 13, doi: 10.1088/0004-637X/740/1/13
  * Zitrin et al. (2009) Zitrin, A., Broadhurst, T., Rephaeli, Y., & Sadeh, S. 2009, ApJ, 707, L102, doi: 10.1088/0004-637X/707/1/L102