# The detection and characterization of highly magnified stars with JWST: prospects of finding Population III

Erik Zackrisson, Adam Hultquist, Aron Kordt, Jose Diego, Armin Nabizadeh, Anton Vikaeus, Ashish Kumar Meena, Adi Zitrin, Guglielmo Volpato, Emma Lundqvist, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Erik Zackrisson, Adam Hultquist, Aron Kordt, Jose Diego, Armin Nabizadeh, et al.. The detection and characterization of highly magnified stars with JWST: prospects of finding Population III. Monthly Notices of the Royal Astronomical Society, 2024, 533 (3), pp.2727-2746. 10.1093/mnras/stae1881. hal-05016709

HAL Id: hal-05016709

https://hal.science/hal-05016709v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The detection and characterization of highly magnified stars with JWST: prospects of finding Population III

Erik Zackrisson<sup>1,2,★</sup> Adam Hultquist,<sup>1</sup> Aron Kordt,<sup>1</sup> Jose M. Diego,<sup>3</sup> Armin Nabizadeh,<sup>1</sup>

Anton Vikaeus<sup>1,†</sup>, <fnref n="1" /> Ashish Kumar Meena<sup>2,§</sup>, <fnref n="4" /> Adi Zitrin<sup>2,§</sup>, <fnref n="4" /> Guglielmo Volpato,<sup>5,6</sup> Emma Lundqvist,<sup>1</sup>

Brian Welch,<sup>7,8</sup> Guglielmo Costa<sup>1,6,9</sup> and Rogier A. Windhorst<sup>10</sup>

Accepted 2024 July 29. Received 2024 July 29; in original form 2023 December 12

## ABSTRACT

Gravitational lensing may render individual high-mass stars detectable out to cosmological distances, and several extremely magnified stars have in recent years been detected out to redshifts $z \approx 6$. Here, we present Muspelheim, a model for the evolving spectral energy distributions of both metal-enriched and metal-free stars at high redshifts. Using this model, we argue that lensed stars will form a highly biased sample of the intrinsic distribution of stars across the Hertzsprung–Russell diagram, and that this bias will typically tend to favour the detection of lensed stars in evolved stages characterized by low effective temperatures, even though stars only spend a minor fraction of their lifetimes in such states. We also explore the prospects of detecting individual, lensed metal-free (Population III) stars at high redshifts using the *James Webb Space Telescope (JWST)*. We find that very massive ($\gtrsim 100 M_\odot$) Population III stars at $z \gtrsim 6$ may potentially be detected by *JWST* in surveys covering large numbers of strong-lensing clusters, provided that the Population III stellar initial mass function is sufficiently top-heavy, that these stars evolve to effective temperatures $\le 15000$ K, and that the cosmic star formation rate density of Pop III stars reaches $\gtrsim 10^{-4} M_\odot \text{cMpc}^{-3} \text{yr}^{-1}$ at $z \approx 6–10$. Various ways to distinguish metal-free lensed stars from metal-enriched ones are also discussed.

**Key words:** gravitational lensing: micro – gravitational lensing: strong – stars: Population III – dark ages, reionization, first stars.

## 1 INTRODUCTION

The very first stars, created from gas with primordial composition, likely started forming at a time when the Universe was $\approx 50–100$ Myr old (cosmic redshift $z \approx 30–50$), thereby bringing the cosmic ‘dark ages’ to an end. Due to the lack of efficient coolants and fragmentation in the chemically unenriched gas at these early epochs, the resulting metal-free (a.k.a. Population III, hereafter Pop III) stars are expected to exhibit a stellar initial mass function (IMF) that differs significantly from that of metal-enriched stars. Current simulations predict this IMF to be top-heavy and possibly logarithmically flat across a very wide mass range (from $\lesssim 1 M_\odot$ and up to several hundred $M_\odot$; for a review, see Klessen & Glover 2023). A substantial Pop III mass fraction may hence be locked up in the form of massive stars ($\gtrsim 10 M_\odot$), with short lifetimes ($\approx 2–40$ Myr for stars in the

$10–1000 M_\odot$ mass range; Yoon, Dierks & Langer 2012). As some of these stars exploded as supernovae, the ambient gas was enriched with heavy elements that initiated the transition to the normal mode of star formation (Population II and I), with much lower characteristic stellar masses as a result. While the Pop III star formation mode likely became cosmologically subdominant already at $z \approx 15–25$ (less than 200 Myr after the big bang), these stars may have continued forming at a low rate for several billion years within rare pockets of gas that remained unpolluted by metals (Magg et al. 2018; Mebane, Mirocha & Furlanetto 2018; Liu & Bromm 2020). However, no massive Pop III stars have ever been directly detected, and – with the exception of supermassive Pop III stars ($\sim 10^5 M_\odot$; e.g. Haemmerlé et al. 2018; Surace et al. 2018, 2019) or Pop III stars at the lowest redshifts ($z \lesssim 1.5$ for $M \le 1000 M_\odot$) – their intrinsic fluxes place them far below the detection limits of both current and upcoming telescopes (Rydberg et al. 2013; Schauer, Drory & Bromm 2020).

Gravitational lensing can boost the fluxes of individual stars beyond the local Universe to detectable levels, and dozens of

<footnote n="1"><sup>1</sup>Observational Astrophysics, Department of Physics and Astronomy, Uppsala University, Box 516, SE-751 20 Uppsala, Sweden</footnote>

<footnote n="2"><sup>2</sup>Swedish Collegium for Advanced Study, Linneanum, Thunbergsvägen 2, SE-752 38 Uppsala, Sweden</footnote>

<footnote n="3"><sup>3</sup>Instituto de Física de Cantabria, (CSIC-UC), Avda. Los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="4"><sup>4</sup>Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, 8410501 Be’er-Sheva, Israel</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica e Astronomia Galileo Galilei, Università di Padova, Vicolo dell’Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="6"><sup>6</sup>INAF – Padova, Vicolo dell’Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="7"><sup>7</sup>Observational Cosmology Lab, NASA Goddard Space Flight Center, Greenbelt, 20771 MD, USA</footnote>

<footnote n="8"><sup>8</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="9"><sup>9</sup>CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, Univ Lyon, Univ Lyon, Ens de Lyon, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="10"><sup>10</sup>School of Earth and Space Exploration, Arizona State University, Tempe, 85287-1404 AZ, USA</footnote>

<footnote>*E-mail: erik.zackrisson@physics.uu.se</footnote>
such stars have already been detected at $z \approx 1-6$ (e.g. Kelly et al. 2018; Welch et al. 2022a; Fudamoto et al. 2024). These highly magnified stars have all been found through observations targeting foreground strong-lensing galaxy cluster fields. In this situation, a background, high-redshift star drifts into the region close to the cluster macrolensing caustic, where it may be subjected to magnifications $\mu > 1000$ due to the combined effects of the macrolensing from the galaxy cluster itself, microlensing from intracluster stars, and potentially also millilensing from globular clusters and dark matter substructures (e.g. Diego et al. 2024a). The magnification is time dependent, both because of the transit of the star across the macrocaustic (time-scale of months to years) and the movement of the star though the corrugated network of microcaustics (time-scale of days to weeks). Most of the lensed stars detected so far have been detected because of their lensing-induced variability (e.g. Kelly et al. 2022; Yan, Jerabkova & Kroupa 2023; Fudamoto et al. 2024), but Earendel, the highest redshift lensed-star candidate currently known, has so far not exhibited any substantial variability (Welch et al. 2022b), possibly because of a very high microlensing optical depth that would result in small magnification fluctuations (Welch et al. 2022a).

Very high levels of gravitational magnification ($\mu \gtrsim 10^3-10^4$) would also be able to lift massive Pop III stars at redshifts as high as $z \approx 20$ into the detection range of the *James Webb Space Telescope (JWST)*, but the probability for such detections is sensitive to the detailed cosmic star formation history and stellar initial mass function of the Pop III stars (Rydberg et al. 2013; Windhorst et al. 2018).

The question of how to distinguish between a lensed, metal-enriched star and a bona fide Pop III object using *JWST* photometry or spectroscopy, and what stellar properties one can realistically hope to derive from such data, remains open.

While many papers have been devoted to the spectral signatures of Pop III galaxies or star clusters (e.g. Tumlinson & Shull 2000; Tumlinson, Giroux & Shull 2001; Schaerer 2002, 2003; Tumlinson, Shull & Venkatesan 2003; Inoue 2011; Zackrisson et al. 2011; Mas-Ribas, Dijkstra & Forero-Romero 2016; Nakajima & Maiolino 2022; Trussler et al. 2023), not much work has been devoted to the spectral signatures of individual lensed stars. Here, we present the Muspelheim<fnref n="1" /> models for the spectral energy distributions of individual, lensed stars at high redshifts, derive the required conditions for detecting highly magnified Pop III stars at $z \gtrsim 6$ and explore the prospects of distinguishing such stars from metal-enriched counterparts using *JWST*.

Section 2 includes an introduction to the computational machinery of the Muspelheim models. In Section 3, we use Muspelheim to explain why lensed, high-$z$ samples of stars obtained using *JWST* photometry will typically result in a highly biased sample of high-mass stars, favouring specific regions in the Hertzsprung-Russell (HR) diagram over others. This turns out to have a large impact on our forecasts for lensed Pop III star detections in Section 4, as some evolutionary pathways predicted for Pop III stars can significantly boost the detection probabilities of such stars in lensing surveys. In Section 5, we explore a number of features that can potentially be used to distinguish Pop III stars from metal-enriched lensed stars. Various uncertainties in the modelling of lensed stars are discussed in Section 6, and Section 7 summarizes our findings.

## 2 THE MUSPELHEIM MODELS

The Muspelheim model predicts the spectral energy distributions (SEDs) of stars as a function of age and initial [Zero Age Main Sequence, (ZAMS)] mass $M_{ZAMS}$ at different metallicities. The model start from a set of stellar evolutionary tracks (Section 2.1) which describe the bolometric luminosity $L$, effective temperature $T_{eff}$, surface gravity log(g), and radius $R$ of stars as a function of age $t$, $M_{ZAMS}$, metallicity ([M/H]) and rotation. SEDs are generated for positions along the tracks separated by $\Delta(t) \approx 10^5$ yr$^2$, or whenever $T_{eff}$ has changed by $\Delta \log_{10} T_{eff} \ge 0.05$ compared to the last saved SED, by matching the $T_{eff}$ and log(g) to a suitable stellar atmosphere spectrum from one of several pre-calculated grids (Section 2.2). This SED is then integrated and rescaled to match the bolometric luminosity of the star. As small mismatches between the $T_{eff}$ and surface gravity log(g) of the track data point and the available pre-calculated stellar atmospheres are inevitable, this rescaling is required to minimize the impact of such defects on the predicted brightness of the star. This results, for every ZAMS mass featured in the stellar evolutionary tracks, in an age sequence of stellar spectra with relatively low spectral resolution (SEDs) that cover the rest-frame ultraviolet to mid-infrared wavelength range. Using these spectra, we derive the ionizing luminosities and redshift-dependent AB magnitudes in a set of *JWST* broad- and medium band filters of these stars, assuming a $\Omega_{\Lambda} = 0.7$, $\Omega_{M} = 0.3$, $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$ cosmology. The Muspelheim model grids (spectra, ionizing fluxes, AB magnitudes in various filter sets) are publicly available (see Data Availability section for access link).

### 2.1 Stellar evolutionary tracks

In this paper, we use the Bonn Optimized Stellar Tracks (BoOST, v. 1.3.; Szécsi et al. 2022) for low-metallicity stars (SMC metallicity; $Z \sim 0.1 Z_{\odot}$) in the $M_{ZAMS} \approx 9-575 M_{\odot}$ range. For Pop III stars, several sets of alternative tracks are used: tracks in the $M_{ZAMS} = 10-1000 M_{\odot}$ range from Yoon et al. (2012), $M_{ZAMS} = 10-100 M_{\odot}$ tracks by Windhorst et al. (2018), $M_{ZAMS} = 9-120 M_{\odot}$ tracks by Murphy et al. (2021) and $M_{ZAMS} = 100-1000 M_{\odot}$ tracks by Volpato et al. (2023). These Pop III tracks differ in their assumptions on e.g. stellar rotation, magnetic torques, mass-loss, and overshooting, and consequently result in different predictions for mass, surface temperature, luminosity as a function on stellar age even for Pop III stars with identical ZAMS masses. We refer to the original papers for detailed descriptions on the input physics of these models, but note that whereas the Windhorst et al. (2018) and Volpato et al. (2023) models are for non-rotating stars only, Yoon et al. (2012) and Murphy et al. (2021) also consider models with various levels of rotation. The Yoon et al. (2012) models cover initial rotational velocities $v_{init} = 0.0-0.6 v_{K}$ for 10–250 $M_{\odot}$ and $v_{init} = 0.0-0.4 v_{K}$ at 300–1000 $M_{\odot}$, where $v_{K}$ is the Keplerian velocity at the equatorial surface. The Murphy et al. (2021) models instead quantify the initial rotation of their models as the fraction of the critical break-up velocity $v_{crit}$, covering $v_{init} = 0.0-0.4 v_{crit}$ throughout the 9–120 $M_{\odot}$ mass interval.

### 2.2 Stellar atmospheres

To produce stellar atmosphere spectra for metal-poor stars, we combine the O- and B-star TLUSTY grids of Lanz & Hubeny

<footnote n="1">¹In Norse mythology, the stars of the night sky are interpreted as glowing sparks from Muspelheim, the realm of fire.</footnote>

<footnote n="2">²Provided that this is possible given the time sampling of the original track data points.</footnote>
(2003, 2007) for $T_{\text{eff}} = 15000-55000$ K with the grid of Lejeune, Cuisinier & Buser (1997) at $T_{\text{eff}} < 15000$ K, both at metallicity [$\text{M/H}$] = -1 and scaled solar abundances. For our Pop III models, we use the TLUSTY stellar atmosphere code (Hubeny & Lanz 1995) to generate a grid of stellar atmosphere models with primordial chemical composition for effective surface temperatures $T_{\text{eff}}$ in the range from 15000 to $\approx 3 \times 10^5$ K, and augment this with $T_{\text{eff}} < 15000$ K with the lowest metallicity models ($\text{[M/H]} = -5$) models from Lejeune et al. (1997) set. The wavelength coverage of these spectra is $\approx 45$ to $\approx 1.6 \times 10^6$ Å for the low-metallicity Lanz & Hubeny (2003, 2007) TLUSTY grids, from $\approx 12-114$ Å (depending on $T_{\text{eff}}$) to $\approx 3 \times 10^6$ Å for the Pop III TLUSTY grids and $\approx 91$ to $\approx 3 \times 10^6$ Å for the Lejeune et al. (1997) set. While the high-$T_{\text{eff}}$ part of the grid (TLUSTY) does take the effects of non-local thermal thermodynamic equilibrium into account, the low-$T_{\text{eff}}$ part (the Lejeune et al. 1997 set) is based on the assumption of local thermodynamic equilibrium. While this will affect the strengths and shapes of absorption lines (e.g. Bergemann 2014), the impact on the overall shape of the SEDs is likely minor, as the continuum of the stellar atmosphere spectra used in the Lejeune have already been subjected to empirical corrections.

To derive a suitable SED for each extracted data point along the stellar evolutionary tracks, these grids of stellar atmosphere spectra are interpolated in $T_{\text{eff}}$ and log(g). Details on this procedure can be found in Appendix A.

While our models assume an approximately primordial chemical composition for all Pop III stellar atmospheres, one should note that scenarios involving chemically homogeneous evolution and dredge-ups (e.g. Yoon et al. 2012; Song et al. 2020; Liu et al. 2021; Volpato et al. 2023, 2024) are expected to alter the surface composition of the star during their lifetimes (see Section 6.5 for a discussion). While this will affect individual spectral features, the impact on the overall photospheric SED as measured by broad-band filter fluxes is likely to be negligible.

## 3 AN OBSERVATIONAL BIAS IN THE HERTZSPRUNG–RUSSELL DIAGRAM

In this section, we use the Muspelheim models for the JWST broad-band filter fluxes of high-redshift stars, as a function of initial mass and age, to demonstrate the existence of a strong selection bias of where in the HR diagram lensed stars are expected to turn up.

Massive, metal-enriched stars typically spend most of their lifetimes at relatively high $T_{\text{eff}}$ ($\gtrsim 25000$ K; Szécsi et al. 2022), but many of the currently known, lensed $z \gtrsim 1$ stars (e.g. Kelly et al. 2018; Diego et al. 2023a; Meena et al. 2023) are found at considerably lower temperatures ($T_{\text{eff}} \lesssim 15000$). We argue that this is due to a strong selection bias for lensed stars, caused by the combined effects of the detection limits (as a function of wavelength) and the lensing magnification distribution. Because of this bias, JWST surveys for lensed, high-z stars will not be able to populate the Hertzsprung–Russell diagram uniformly, but will instead preferentially pick up stars in certain regions of the diagram (see also Diego et al. 2024b, for a discussion on this). Consequently, stars that evolve into such regions will be easier to detect, and this has direct consequences for the detection probabilities of Pop III stars, since not all proposed evolutionary paths for such stars will venture into the HR diagram regions where the detection prospects are optimized.

In general, the lensed-star selection function (i.e. the bias that determines what stars are lifted above the observational detection threshold through a combination of macro- and microlensing) is

expected to be highly complex. If star formation were unclustered, we argue that this selection function should depend on the wavelength-dependent detection limits, on the shape of the magnification probability distribution, on the shape of the stellar IMF, on the evolution of stars of a given initial mass across the HR diagram, and the dust obscuration and reddening along the line of sight (with potential contributions from the circumstellar medium of the host star, the ISM at its location of birth, the ISM elsewhere in its host galaxy and within the cluster lens). The path that a star of a given initial mass follows across the HR diagram in turn depends on metallicity, rotation rate, the impact of binary evolution, and the potential mergers of stars. The fact that most stars are likely to be born in star clusters brings additional complications, since the difference in magnification between one star and its immediate neighbours dictates whether the observer detects a single lensed star or the blended light from many. IMF sampling effects in small star clusters (e.g. Yan et al. 2023) and the potential ejection of some massive stars from their birth clusters may also play a role. Here, we make no attempt to provide a definitive model for all of these effects, but simply use Muspelheim to demonstrate the likely biases in the highly simplified situation of single-metallicity stars obeying a specific IMF, while neglecting binary evolution, clustering and dust obscuration effects.

As stars age, they typically evolve from a long-lasting high $T_{\text{eff}}$ state on the main sequence, to more short-lived states at lower $T_{\text{eff}}$ (although rapid rotation, binary evolution and mergers can alter this behaviour; e.g. de Mink et al. 2009; Yoon et al. 2012; Szécsi et al. 2022; Wang et al. 2022). This is typically accompanied by an increase in $L_{\text{bol}}$, although this increase can be rather modest for the most massive stars, which are predicted to evolve almost horizontally across the Hertzsprung–Russell diagram (e.g. Yoon et al. 2012; Windhorst et al. 2018; Volpato et al. 2023).

One may perhaps naively expect that most lensed stars would be detected in their more-long lasting main-sequence states, but there are observational biases that favour the detection of more evolved stars. At fixed $L_{\text{bol}}$, a low-$T_{\text{eff}}$ star attains a higher $f_v$ flux at the peak of its SED compared to that of a high-$T_{\text{eff}}$ star<fnref n="3" />. Since the JWST/NIRCam detection limits are fairly constant in $f_v$ units across the $\approx 1-5$ µm range<fnref n="4" />, this means that a high-z, low-$T_{\text{eff}}$ star requires a lower gravitational magnification $\mu$ to be rendered detectable compared to a high-$T_{\text{eff}}$ star of the same $L_{\text{bol}}$ and at the same redshift. Since the probability $P$ to attain a certain magnification approximately goes as $P(> \mu) \propto \mu^{-2}$ in the source plane (up to some limit; see Section 3.1), short-lived, low-$T_{\text{eff}}$ evolutionary states that require low $\mu$ to be detected can in lensed samples trump long-lived, high-$T_{\text{eff}}$

<footnote n="3">³This is because the luminosity in a fixed wavelength interval $\Delta\lambda$ goes as $L(\Delta\lambda) \propto f_v \lambda^{-2} \Delta\lambda$, so that $f_v$ fluxes at long wavelengths contribute less to $L_{\text{bol}}$ than a similar $f_v$ at short wavelengths. A low-$T_{\text{eff}}$ star, for which the SED peaks at long wavelengths, will therefore need a higher peak $f_v$ flux to reach the same $L_{\text{bol}}$ as a high-$T_{\text{eff}}$ star, for which the SED peaks at short wavelengths. The increase in peak brightness with decreasing $T_{\text{eff}}$ at fixed $L_{\text{bol}}$ holds for both stellar atmosphere spectra and blackbody spectra if flux is expressed in $f_v$ units, but not if it's expressed in $f_\lambda$ units.</footnote>

<footnote n="4">⁴Across the wide NIRCam filters F090W, F115W, F150W, F200W, F277W, F356W, and F444W filters used in this paper, the $f_v$ detection limits are within a factor of 3 of each other, at fixed signal-to-noise ratio and exposure time (Rieke et al. 2023). When the peak $f_v$ flux of a stellar SED is redshifted into the wavelength range of this filter set, the difference between the maximum continuum $f_v$ seen in the SED and that captured photometrically is moreover typically minor, due to the shape of stellar SEDs.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The reason why lensed stars form a biased sample of the intrinsic distribution in the HR diagram. (a) A selection of evolutionary tracks at different ZAMS masses from the BoOST SMC set at $Z/Z_\odot \sim 0.1$. Coloured markers along each track indicate time-steps of $\approx 10^5$ yr. Please note how, in the case of the 15 and 30 $M_\odot$ tracks, there are no age markers between log $T_{\text{eff}} \approx 4.4$ and $\approx 3.6$, as the evolution between these states is very rapid. (b) Redshift $z=6$ SEDs corresponding to the evolutionary stages labeled 1–4 along the 55 $M_\odot$ track in panel (a). The vertical dashed line indicates the Lyman $\alpha$ limit at this redshift, shortward of which the flux is expected to be absorbed by the neutral IGM. As seen, the SEDs become considerably brighter with decreasing $T_{\text{eff}}$ throughout the $\approx 1–5$ µm range of JWST/NIRCam, even though the evolution in $L_{\text{bol}}$ is modest. (c) The brightest NIRCam wide-band magnitude attained as a function of $T_{\text{eff}}$ along the 55 $M_\odot$ track in panel a, at redshifts $z=1, 3, 6$, and 10. While the brightness increases with decreasing $T_{\text{eff}}$ at $z=1–3$, this trend starts to turn over at $z=6$, and at $z=10$, the lowest-$T_{\text{eff}}$ point at the end is no longer the brightest, since the peak of the $\approx 4000$ K (log $T_{\text{eff}} \approx 3.6$) SED no longer falls within the NIRCam bands at $z=10$. (d) The relative probability of detection as a function of $T_{\text{eff}}$ along the 55 $M_\odot$ track for a detection threshold of 28 AB mag across the JWST/NIRCam broad-band filters. In the case where a 55 $M_\odot$ star is detected through lensing (and a $P(> \mu) \propto \mu^{-2}$ magnification probability distribution is assumed in the source plane), the coloured orange, green, and blue lines represent the observed probability distribution for different $T_{\text{eff}}$ at $z=3, 6$, and 10. The black line instead describes the intrinsic distribution of $T_{\text{eff}}$ across the lifespan of the star, without any observational bias. While most of the lifetime of this star is spent at log $T_{\text{eff}} \gtrsim 4.6$ (black line; this means that there are more stars with this initial mass at log $T_{\text{eff}} \gtrsim 4.6$ than at lower $T_{\text{eff}}$ in any source plane with active star formation), such high temperatures are very unlikely to be detected. Instead, it is the red supergiant state at log $T_{\text{eff}} \approx 3.6$ that is the most likely to be detected at $z=3$ and $z=6$.</figcaption></figure>

states that require high $\mu$ for detection (see Appendix B for further discussion on this).

We demonstrate this effect for a 55 $M_\odot$ star at $z=1–10$ in Fig. 1. In general, this bias will favour detection of the lowest-$T_{\text{eff}}$ states along the stellar evolutionary tracks, as long as these are not too short-lived, and as long as the peak of the SEDs of these low-$T_{\text{eff}}$ stars falls within the wavelength range of JWST/NIRCam at the redshift of the star. As the redshift is increased, the most favoured $T_{\text{eff}}$ gradually shifts towards higher $T_{\text{eff}}$, since the $f_v$ peak of lower-$T_{\text{eff}}$ stars redshifts out of NIRCam range. However, the selection function still looks very different from one based simply on how long stars intrinsically occupy different $T_{\text{eff}}$ states.

## 3.1 Modelling the selection function of lensed stars

Here, we use Muspelheim to simulate the expected distribution of low-metallicity ($Z \sim 0.1 Z_\odot$) lensed stars at $z=1–10$ in the HR diagram for a JWST/NIRCam survey of cluster-lens fields. To do this, we assume a fixed detection limit (in the 27–30 AB magnitude range) in the NIRCam filters F090W, F115W, F150W, F200W, F277W, F356W, and F444W.

For plotting purposes, we at each redshift interpolate the Muspelheim models for the AB magnitudes of stars along the Szécsi et al. (2022) evolutionary tracks to a finer mass resolution, and derive the brightest AB magnitudes in the JWST/NIRCam bands as a function of redshift at each data point along these interpolated tracks.
To model the intrinsic distribution of stellar masses, we adopt a stellar initial mass function (IMF) with the same high-mass slope as the Kroupa (2001) universal IMF ($dN/dM \propto M^{-2.3}$) throughout the 10–150 M$_{\odot}$ mass range.

We then run a Monte Carlo simulation that draws a sample of stars that obey the adopted IMF, randomly assigns ages up to 30 Myr and rejects stars which, based on their initial mass, have expired at their selected age. This corresponds to a situation where $M_{ZAMS} \ge 10M_{\odot}$ stars of any age may drift across the source-plane caustic of the clusters surveyed. This is a fair assumption in a case where the sample of lensed stars comes from many highly magnified galaxies with different star formation histories. However, this may not accurately reflect a situation in which the sample of lensed stars comes from a single, highly magnified galaxy, since the actual age distribution of stars would then be dictated by the detailed star formation history of that particular object.

For every Monte Carlo generated star, a random magnification $\mu$ is assigned under the assumption of a source-plane probability distribution function $P(\mu) \propto \mu^{-3}$ (or equivalently $P(> \mu) \propto \mu^{-2}$; see Schneider, Ehlers & Falco 1992 for an analytical derivation of this distribution), and the brightest NIRCam wide-band flux is shifted accordingly to predict the brightest apparent magnitude across the NIRCam bands of the star. Stars that do not reach the adopted JWST detection threshold in at least one of the NIRCam filters assumed are rejected from the simulated sample.

We note that in the resulting sample of the lensed stars that make it above the JWST detection threshold, the $P(> \mu) \propto \mu^{-2}$ distribution will no longer be obeyed, even though stars were initially assigned magnifications from this distribution. This happens because some stars can be rendered detectable at lower magnifications than others. The reason for this bias can be understood through the following thought experiment: Let's assume that there are only two types of stars, A and B, which both retain fixed positions in the HR diagram until they expire. Due to having different intrinsic brightness, star A may be detected when $\mu > 1000$, whereas star B requires $\mu > 10000$ to be lifted above the detection threshold. If these stars were equally common, then one would expect to detect B one hundred times less often than A ($P_B(\mu > 10000)/P_A(\mu > 1000) = 10^{-2}$). However, if B is ~100 times more common than A (due to the IMF or a longer lifetime), then the probability to detect A and B become comparable, and hence magnifications of $\mu > 1000$ and $\mu > 10000$ turn out about equally likely in the observed sample of lensed stars, which violates the $P(> \mu) \propto \mu^{-2}$ distribution. In practice, however, the magnification distributions of our simulated samples are never flat, and the highest magnifications are still the rarest, but the mean or median magnifications may vary substantially with the adopted detection limit and choice of stellar evolutionary tracks.

Throughout this paper, we only consider stars with magnifications in the $\mu = 1 \times 10^2-5 \times 10^4$ range. The lower limit has been set in order not to exclude cases which may resemble some of the faintest lensed stars detected at low redshifts with JWST through microlensing-induced variability (e.g Fudamoto et al. 2024). The upper limit is set by the fact that higher magnifications are highly unlikely in realistic settings, as argued by Diego et al. (2018). In some figures, we also present results exclusively for the high-magnification regime of $\mu = 3 \times 10^3-5 \times 10^4$, which corresponds to a range of magnifications that $\gtrsim 1$ pc star clusters cannot reach [see equation (9) and Fig. C1 and the associated discussion]. This higher-magnification interval may be more appropriate in situations where lensed-star candidates do not show significant variability, and there may be some ambiguity whether these objects are high-magnification stars or lower-magnification star clusters. In such cases, the distinction

can be made on the likely macro-magnification at the position of the object. We stress, however, that many of the variability-selected lensed stars detected at low redshifts by JWST may also be in this high-magnification range (Broadhurst et al. 2024).

While the absolute scaling of the $P(> \mu) \propto \mu^{-2}$ magnification distribution is important for the absolute number of lensed stars expected to be detected in a given survey, we are in this section only concerned with the relative distribution across the HR diagram of the lensed stars that are detected (i.e. where in the HR diagram lensed stars are the most likely to turn up, not the number of stars that turn up in those regions), for which the scaling factor does not matter.

Stars that attain their high magnifications due to microlensing by stars in a foreground galaxy cluster, for instance because they were selected based on variability criteria in multi-epoch observations, may show a slightly different HR diagram bias compared to stars selected due to their close proximity of a cluster macrolensing caustic. In general, the size of a star sets an upper limit on its microlensing-induced magnification, which would make microlensing disfavour detection of stars with the largest radii.

While the exact size–magnification limit depends on the gradient of the lensing potential at the position of the event, and the surface mass density in stars at this location in the cluster lens (e.g. Venumadhav, Dai & Miralda-Escudé 2017; Meena et al. 2023), we here adopt an approximate size–magnification limit of the form

$$ \max(\mu_{\text{micro}}) \approx A_0(D_S/1\,\text{Gpc})^{1/2}(R/10\,\text{R}_\odot)^{-1/2}, \quad (1) $$

with $A_0$ either $1 \times 10^4$ or $5 \times 10^4$, where the higher value corresponds to a lower contribution from microlenses to the surface mass density. Here, $D_S$ is the angular size distance to the source at redshift $z$, and $R$ is the radius of the star. The $A_0 = 1 \times 10^4$ value is similar to the lowest value considered by Meena et al. (2023), while the $A_0 = 5 \times 10^4$ is the highest we can consider given our overall cut-off in overall the magnification distribution. This range of $A_0$ corresponds to a change in the surface mass density of microlenses by a factor of $\approx 8$ (Venumadhav et al. 2017).

## 3.2 Selection bias for metal-enriched stars obeying a standard IMF

To illustrate the effects of the resulting selection biases, we in Fig. 2 compare the intrinsic distribution of stars with initial masses 10–150 M$_{\odot}$ in the HR diagram with the most likely regions for lensed stars to appear at $z = 2$ in a JWST survey reaching 27.0 AB magnitudes, given our $Z/Z_{\odot} \sim 0.1$ models. The intrinsic distribution is dominated by $M_{ZAMS} \lesssim 55M_{\odot}$ stars on the main sequence at $\log T_{\text{eff}} > 4.3$, and by $M_{ZAMS} \lesssim 20M_{\odot}$ stars that have evolved to red supergiants at $\log T_{\text{eff}} \approx 3.6$.

In a lensed sample with $\mu = 1 \times 10^2 - 5 \times 10^4$ and a microlensing size–magnification limit with $A_0 = 5 \times 10^4$ (Fig. 2a), the observed distribution instead gets dominated by cooler stars in the $T_{\text{eff}} < 20000$ K ($\log T_{\text{eff}} < 4.3$) range, and with initial masses from $M_{ZAMS} \gtrsim 10 M_{\odot}$ all the way up to $\approx 150 M_{\odot}$. Adopting a more restrictive microlensing size–magnification limit (equation 1) with $A_0 = 1 \times 10^4$ (Fig. 2b) removes some of the lower-mass red supergiants with initial masses $< 30 M_{\odot}$, which have large radii ($\approx 300-1000 R_{\odot}$) and require $\mu > 1000$ to be rendered detectable.

In Figs 2 (c) and (d), we show the corresponding distributions for the highest-magnification stars with $\mu \ge 3000$. In this subset, one finds an increase in the fraction of $T_{\text{eff}} \gtrsim 20000$ K ($\log T_{\text{eff}} \gtrsim 4.3$ stars, which require these very high magnifications to get boosted above the detection threshold. In the $A_0 = 5 \times 10^4$ (Fig. 2c) case,
<figure><img src="image_4.png" /></figure>

**Figure 2.** The selection bias of lensed, low-metallicity ($Z \approx 0.1 Z_\odot$), $M_{ZAMS} = 10–150 M_\odot$ stars in the Hertzsprung-Russell diagram at $z=2$, for two different microlensing size-magnitude relations, and two different lower limits on the magnification. The coloured, filled contours indicate the regions in which ≈ 95 per cent of the stars will fall, under different assumptions on the selection function. In all panels, the dark orange regions mark the intrinsic distribution of the $M_{ZAMS} = 10–150 M_\odot$ stars (i.e. only considering the stellar IMF and the time spend in different evolutionary stages, without considering any brightness or lensing bias). The light blue region represents the regions occupied by the stars most likely to end up in JWST/NIRCam surveys reaching ≤ 27.0 AB mag in lensed fields in the case of a source-plane $P(> \mu) \propto \mu^{-2}$ magnification distribution throughout the range $\mu = 1 \times 10^{2–5} \times 10^4$ (upper panels) or $\mu = 3 \times 10^3–5 \times 10^4$ (lower panels), with an additional microlensing size-magnification limit (equation 1) with $A_0 = 1 \times 10^4$ (left panels) and $A_0 = 5 \times 10^4$ (right panels). In all panels, lensed samples produce a strong selection bias, as seen in the difference between the dark orange and light blue regions. This bias typically favours stars that have evolved significantly away from their ZAMS effective temperatures. Only the most massive stars are likely to be detected at $\log T_{\text{eff}} > 4.0$ whereas less massive stars may be detected on the red supergiant branch $\log T_{\text{eff}} \approx 3.6$. The prospects of detecting red supergiants in the lensed sample are, however, sensitive to the adopted microlensing size-magnification constraint, since a lower $A_0$ (as in panel b) makes it harder for stars with large radii to reach the magnification required for detection. The subset of high-magnification stars ($\mu \ge 3000$; panels c and d) displays a wider range in effective temperature that extends to higher $\log T_{\text{eff}}$, since very high magnifcaitons are required to detect the hottest stars. In the case of the most restrictive microlensing size-magnification relation considered ($A_0 = 1 \times 10^4$), the $\mu \ge 3000$ sample (panel d) rejects low-$T_{\text{eff}}$ stars due to their large radii, leaving only the $\log T_{\text{eff}} > 4.2$ stars detectable.

stars with initial masses ≈ 100–150 $M_\odot$ can in fact be detected at their ZAMS temperatures of $T_{\text{eff}} \gtrsim 50000$ K ($\log T_{\text{eff}} \gtrsim 4.7$). In the case of $A_0 = 1 \times 10^4$ (Fig. 2d), the $T_{\text{eff}} < 12000$ K stars fall out of the sample because of their large radii and the size-magnification limit set by equation (1). Some of the high-mass, highest-log $T_{\text{eff}}$ stars also get removed by this criterion, since these stars have radii ~ 10 $R_\odot$ already on the ZAMS and require $\mu > 10^4$ to be detectable because of their high $T_{\text{eff}}$.

Using $A_0 = 5 \times 10^4$, and $\mu = 1 \times 10^2–5 \times 10^4$, we in Fig. 3 indicate the most likely regions in the HR diagram for $Z/Z_\odot \sim 0.1$ stars to appear at $z = 1, 3, 6$, and 10, in the case of two deeper JWST photometric surveys that reach 28.0 and 30.0 AB magnitudes.

At $z=1$ and $z=3$, the resulting distributions are similar to those at $z=2$ in Fig. 2(a), but extend to higher $T_{\text{eff}}$ at these fainter detection thresholds. At $z=6$, the chance of detecting red supergiants from stars with initial masses $M < 55 M_\odot$ (28 AB mag survey limit) or $M < 30 M_\odot$ (30 AB mag survey limit) diminishes, and stars with $T_{\text{eff}}$ up to their ZAMS values enter the sample for initial masses $\gtrsim 100 M_\odot$ (28 AB mag) or $\gtrsim 40 M_\odot$ (30 AB mag). At $z=10$ the chance of detecting red supergiants at $\log T_{\text{eff}} < 3.6$ is gone, since the SED peaks of these stars have redshifted far outside the JWST/NIRCam wavelength range.

In the redshift range $z=1–6$, it is interesting to note that for the brighter detection limits of 27–28 AB mag [Figs 2 (a–c) and 3
<figure><img src="image_5.png" /></figure>

**Figure 3.** The selection bias of lensed, low-metallicity ($Z \approx 0.1 Z_\odot$) $M_{ZAMS} = 10–150 M_\odot$ stars in the Hertzsprung–Russell diagram at (a) $z = 1$, (b) $z = 3$, (c) $z = 6$, and (d) $z = 10$ under the assumption of a source-plane $P(> \mu) \propto \mu^{-2}$ magnification distribution with $\mu = 1 \times 10^2–5 \times 10^4$ and an additional size–magnification limit (equation (1)) with $A_0 = 5 \times 10^4$. The coloured, filled contours indicate the regions in which $\approx 95$ per cent of the lensed stars will appear, in the case of a 28.0 AB mag (blue patch, thin contour) and 30.0 AB mag (red patch, thick contour) survey limit, respectively. The prospects of detecting red supergiants drops drops at the highest redshifts, since the SED peak of these stars moves out of JWST/NIRCam range, and at $z = 10$ (panel d), no such stars are expected to be seen.

(a–c)], the relative probability of lensed stars to turn up in the 4000–15 000 K range ($\log(T_{\text{eff}}) \approx 3.6–4.2$) becomes substantial, in agreement with recent detections of lensed stars (Kelly et al. 2018; Chen et al. 2022; Welch et al. 2022b; Furtak et al. 2024; Meena et al. 2023; Diego et al. 2023a, b) at these redshifts, even though stars in this temperature range are exceedingly rare.

We stress that we in Fig. 3 only show the expected HR diagram distribution of stars with ZAMS mass $\ge 10 M_\odot$. However, at the lowest redshifts (e.g. $z = 1–3$), it is possible that stars with ZAMS mass $< 10 M_\odot$ could also enter surveys reaching 28–30 AB mag, since the contours extend down to the supergiant branch of the lowest-mass tracks plotted.

### 3.3 Selection bias for Pop III stars with top-heavy IMFs

In Fig. 4, we depict the same selection bias (here with $A_0 = 5 \times 10^4$, 28.0 and 30.0 AB mag survey limits) as in Fig. 3, but now in the case of Pop III tracks from Murphy et al. (2021) and Yoon et al. (2012) with no initial rotation, under the assumption of a top-heavy, log-normal

IMF (Tumlinson 2006, with characteristic mass $M_c = 10 M_\odot$ and $\sigma = 1$), for $M = 10–120 M_\odot$ stars at $z = 6$.

A fundamental difference between the Murphy et al. (2021) and Yoon et al. (2012) tracks is that the latter predict non-rotating Pop III stars to reach a $\approx 4000$ K stage for $M_{ZAMS} \ge 30 M_\odot$ whereas the former stay above 15 000 K until the end of the tracks. As we have previously shown that low-$T_{\text{eff}}$ stages, even if very short-lived, can have a significant impact on the detectability of lensed stars, this difference leads to very different predictions on where in the HR diagram lensed Pop III stars are the most likely to be detected.

As a result, the Yoon et al. (2012) tracks favour the appearance of red supergiants in samples with detection limits at both 28.0 and 30.0 AB magnitudes, whereas the Murphy et al. (2021) tracks produce distributions that are limited to the $T_{\text{eff}} \gtrsim 15000$ K range.

It is worth noting that lensed stars at the peak temperatures of Pop III stars during early stages of evolution ($\sim 10^5$ K) only make up a significant fraction of the simulated sample in the case of the Murphy et al. (2021) tracks and a survey limit of 30.0 AB magnitudes. This extension towards very high $T_{\text{eff}}$ represents an important difference
<figure><img src="image_6.png" /></figure>

**Figure 4.** The most likely positions in the HR diagram of lensed $z = 6, 10–120 M_{\odot}$ Pop III stars following the (a) Murphy et al. (2021) and (b) Yoon et al. (2012) stellar evolutionary tracks for non-rotating stars, under the assumption of a source-plane $P(> \mu) \propto \mu^{-2}$ magnification distribution with $\mu = 1 \times 10^{2}-5 \times 10^{4}$ and a size–magnification limit (equation 1) with $A_{0} = 5 \times 10^{4}$. The coloured, filled contours indicate the regions in which $\approx 95$ per cent of the lensed stars will appear, in the case of a 28.0 AB mag (blue patch, thin contour) and 30.0 AB mag (red patch, thick contour) survey limit, respectively. A top-heavy, log-normal IMF (Tumlinson 2006, with characteristic mass $M_{c} = 10 M_{\odot}$ and $\sigma = 1$) is assumed for all scenarios depicted. In both panels, stars captured by a 30 AB mag sample extend to somewhat higher $T_{\text{eff}}$ than those captured in the 28 AB mag sample. It is only in the 30 AB mag case, and for the Murphy et al. (2021) tracks (panel a), that lensed stars at the characteristically high main-sequence temperatures ($T_{\text{eff}} \sim 10^{5}$ K) of $\gtrsim 100 M_{\odot}$ Pop III stars make up any significant fraction of those rendered detectable. Because the Murphy et al. (2021) tracks do not extend to $T_{\text{eff}} < 15000$ K for Pop III masses in the range plotted, whereas the Yoon et al. (2012) ones do, the contours in panel a appear more compact along the $T_{\text{eff}}$ directions than in panel (b).

between the predicted HR–diagram distribution of lensed 10–120 $M_{\odot}$ Pop III stars (Fig. 4) and the corresponding distribution for lensed metal-enriched stars at $z = 6$ (Fig. 3c), as the latter do not appear at $T_{\text{eff}} > 60000$ K (log $T_{\text{eff}} > 4.8$) given our chosen set of tracks. However, even if $\sim 10^{5}$ K objects were to appear in a sample of lensed stars, measuring their temperatures sufficiently well to constrain their location to this part of the HR diagram may be extremely challenging at faint magnitudes, as discussed in Section 5. Moreover, Wolf–Rayet stars can also attain temperatures of $\sim 10^{5}$ K, and alternative sets of tracks for metal-enriched stars that include this phase (e.g. Eldridge & Vink 2006) may shift the predictions to higher $T_{\text{eff}}$ than shown in Fig. 3(c).

The predicted distribution of lensed Pop III stars across the HR diagram also depends on the shape and upper mass limit of the IMF, which we demonstrate in Fig. 5 in the case of Yoon et al. (2012) tracks and a survey limit of 30 AB mag. Here, we assume that Pop III stars can form up to a mass of 500 $M_{\odot}$ and show the results for two different top-heavy IMFs: a log-normal IMF (Tumlinson 2006, with characteristic mass $M_{c} = 10 M_{\odot}$ and $\sigma = 1$) and a $dN/\log M = $ constant IMF (e.g. Schauer et al. 2022).

The higher upper mass limit adopted for Pop III stars compared to Fig. 4(b) does not have any dramatic effects on the $T_{\text{eff}}$ distribution of lensed stars across the HR diagram for either the log-normal $M_{c} = 10 M_{\odot}$ IMF, or the $dN/\log M = $ constant IMF. However, the probability density function does extend to higher $L_{\text{bol}}$ and $M_{\text{ZAMS}}$, reaching $\approx 500 M_{\odot}$ and $\log(L/L_{\odot}) \approx 7$. If one assumes that the birth of such very massive stars is effectively limited to the Pop III mode of star formation, then the detection of lensed stars in this part of the HR diagram would constitute an indication that Pop III objects have been detected even without detailed metallicity measurements (Section 5). However, ensuring that a given lensed star belongs to this part of the diagram requires that $T_{\text{eff}}$ is well constrained and that strong upper limits on the magnification can be set to infer a robust lower limit on the luminosity.

<figure><img src="image_7.png" /></figure>

**Figure 5.** IMF effects on the HR diagram distribution of lensed Pop III stars at $z = 6$. The coloured, filled contours indicate the regions in which $\approx 95$ per cent of the lensed stars are predicted to appear, in the case of a $dN/\log M = $ constant IMF (indigo patch, thin contour) and a lognormal $M_{c} = 10 M_{\odot}$ and $\sigma = 1$ IMF (blue patch, thick contour), both extending from 10–500 $M_{\odot}$. We here adopt Yoon et al. (2012) stellar evolutionary tracks for non-rotating Pop III stars, a source-plane $P(> \mu) \propto \mu^{-2}$ magnification distribution with $\mu = 1 \times 10^{2}-5 \times 10^{4}$, a size–magnification limit (equation 1) with $A_{0} = 5 \times 10^{4}$ and a JWST/NIRCam survey with detection limit 30.0 AB mag. The higher upper mass limit adopted here compared to Fig. 4 (b) does not have any dramatic impact on the $T_{\text{eff}}$ distribution, but favours the detection of more massive and luminous stars.

In summary, the identification of candidates for lensed Pop III stars based solely on their position in the HR diagram, without any additional metallicity constraints (Section 5), may be very difficult.
## 4 PREDICTED DETECTION RATES FOR HIGHLY MAGNIFIED POP III STARS

While Section 3.1 focused on the relative probability distribution of stars in the HR diagram of lensed stars<fnref n="5" />, we here attempt to estimate the absolute probabilities of detecting Pop III stars at $z \gtrsim 6$ in a JWST survey covering a large number of cluster-lens fields.

When imaging cluster-lens fields, highly magnified stars usually manifest themselves as point sources within lensed arcs, but could in principle also appear as point sources along the macrolensing critical curve of the cluster lens but without any obvious association to an arc. The latter case may arise when the magnification of the star is sufficient to render it detectable, but the host galaxy (which would be subject to a much smaller magnification) remains below the detection threshold. This situation has so far not been seen, but could potentially become increasingly important at the very highest redshifts, where a significant fraction of the cosmic star formation activity may take place in galaxies too faint to be seen through macrolensing even in very deep JWST images (e.g. Wu & Kravtsov 2024). This scenario could also be relevant for Pop III stars forming early on in minihalos, which are not expected to have any association with large galaxies.

In what follows, we will focus on the momentary brightness of lensed Pop III stars at high redshifts, and attempt to estimate the probability that a lensed Pop III star appears in a single-epoch imaging campaign covering several cluster-lens fields. The time-scale over which stars detected this way will remain observable once discovered – which depends on the time-scale of the crossing of macrolensing caustics and the impact of microlensing – is outside the scope of the current paper.

Following, Windhorst et al. (2018), we adopt the relation for fold caustics for the macrolensing magnification $\mu_{\text{macro}}$ as a function of source–plane distance $d_{\text{arcsec}}$ from the caustic in arcseconds:

$$ \mu_{\text{macro}} = \frac{B_0}{(d_{\text{arcsec}})^{1/2}}, \qquad (2) $$

where $B_0$ is a scaling factor that varies from cluster to cluster, in the $B_0 \approx 10-20$ arcsec$^{1/2}$ range for some of the most well-known cluster lenses (Windhorst et al. 2018). Strictly speaking, $B_0$ also varies with source redshift, but the dependence is expected to be very weak in the $z \gtrsim 6$ regime studied in this paper, and we therefore treat it as a constant throughout this paper.

Given an imaging survey with a set of iJWST/NIRCam filters with detection limits $m_{ABlimit, i}$, one can—for a star with initial (ZAMS) mass $M$, age $t$ and intrinsic apparent magnitude $m_{M, t, i}(z)$ – derive the minimum macrolensing magnification $\mu_{\text{detection}}$ required to bring the star with mass $M$ and age $t$ above the detection limit in at least one filter as:

$$ \mu_{\text{detection}}(z) = 10^{\frac{\min(m_{M, t, i}(z)-m_{ABlimit, i})}{2.5}}. \qquad (3) $$

Neglecting microlensing, the maximum distance from the caustic at which such a star may be detected in this survey is then:

$$ \max(d_{\text{arcsec}(z)}) = \left( \frac{B_0}{\mu_{\text{detection}(z)}} \right)^2. \qquad (4) $$

Given this relation and an assumed source–plane length $L_{\text{arcsec}}$ of the caustic in arcseconds, we may estimate the angular surface area $\theta^2$ (in arcsec$^2$) over which a given star will attain a magnification equal to

or higher than $\mu_{\text{detection}}(z)$ (at separations $0 < d_{\text{arcsec}} < \max(d_{\text{arcsec}})$) as

$$ \theta^2 (\mu \ge \mu_{\text{detection}}(z)) = \max(d_{\text{arcsec}}(z))L_{\text{arcsec}}. \qquad (5) $$

By adopting a certain redshift survey depth $\Delta z$, this angular surface area $\theta^2$ can be converted into a redshift-dependent comoving volume $V_{z, \Delta(z)}$ inside which objects will attain $\mu \ge \mu_{\text{detection}}(z)$:

$$ V_{z, \Delta(z)}(\mu \ge \mu_{\text{detection}}(z)) = V(z, \Delta z)_{\text{arcsec}^{-2}} \theta^2(\mu \ge \mu_{\text{detection}}(z)), \qquad (6) $$

where $V(z, \Delta z)_{\text{arcsec}^{-2}}$ is the comoving volume at redshift $z$ that corresponds to depth $\Delta(z)$ per arcsec$^2$. We stress that this treatment only takes macrolensing into account, and neglects the way microlensing is likely to alter the magnification distribution (e.g. Diego 2019; Palencia et al. 2024).

The number of expected detections of lensed Pop III stars in that $\Delta z$ slice can then be derived by estimating the number of luminous Pop III stars within the corresponding volume. Since the magnification required for detection changes dramatically as a function of age of Pop III stars, we sum over $j$ time steps that represent different evolutionary states of the star, each with a duration $\Delta t_j$, and over all redshift slices $z$ for a single cluster:

$$ N_{\text{detection}}(M) \approx \sum_z \sum_j \frac{\text{SFRD}(z) f(M) V_{z, \Delta(z)} (\mu \ge \mu_{\text{detection}}(z)) \Delta t_j}{M}. \qquad (7) $$

Here, $N_{\text{detection}}(M)$ is the expected number of lensed stars of initial mass $M$, SFRD($z$) is the comoving cosmic star formation rate density (in units of $M_{\odot} \text{yr}^{-1} \text{cMpc}^{-3}$) of Pop III stars at redshift $z$, and $f(M)$ is the IMF-dependent mass fraction locked up in Pop III stars with initial mass $M$.

Throughout this paper, we adopt $B_0 = 15$ and $L_{\text{arcsec}} = 100$ arcsec as fiducial values for strong-lensing clusters. As long as $N_{\text{detection}} \ll 1$, $N_{\text{detection}}$ can also be directly interpreted as the approximate probability of detection ($P_{\text{detection}} \approx N_{\text{detection}}$) per cluster field surveyed.

Detecting Pop III stars in JWST surveys of cluster lens fields requires both a relatively high cosmic star formation rate density (SFRD) of Pop III stars at the relevant redshifts, and a top-heavy Pop III IMF. Here, we limit ourselves to $z \gtrsim 6$, since this is the regime that simulations of Pop III stars typically focus on, and for which numerous Pop III SFRD predictions have been published.

Based on the computational machinery outlined, we find that a SFRD of $\sim 10^{-4} M_{\odot} \text{cMpc}^{-3} \text{yr}^{-1}$ somewhere in the $z \approx 6-10$ range is required to make detections realistic in the case of very top-heavy IMFs, and that even higher SFRDs are required for the less extreme IMFs. Adopting Liu & Bromm (2020) as our reference SFRD (which reaches a peak of $\approx 10^{-4} M_{\odot} \text{cMpc}^{-3} \text{yr}^{-1}$ at $z \approx 10$ and drops by a factor of $\approx 3$ by $z = 6$) for Pop III stars, we in Fig. 6(a) show the expected detection rate in a 28.5 AB mag JWST/NIRCam survey for five different IMFs covering the $1-500 M_{\odot}$ mass range, under the assumption of the Yoon et al. (2012) stellar evolutionary tracks for non-rotating Pop III stars. These IMFs include a $dN/dM \propto M^{-2.35}$ Salpeter slope IMF, three Tumlinson (2006) log-normal IMFs with $\sigma = 1$ and characteristic masses $M_c = 10 M_{\odot}$, $30 M_{\odot}$, and $60 M_{\odot}$, plus the Magg et al. (2016) $dN/d\log(M) = \text{constant}$ IMF. Since we are working with discrete stellar evolutionary tracks, each track is in equation (7) assigned a mass fraction $f(M)$ derived from the continuous $1-500 M_{\odot}$ IMF. The lower and upper limits of the mass interval attached to each track are set by the midpoints in initial mass

<footnote n="5">⁵Which provides answers to questions of the form ‘if a lensed star is detected at a given redshift, what properties is it most likely to have’?</footnote>
<figure><img src="image_8.png" /></figure>

**Figure 6.** (a) Expected number of detections of lensed Pop III stars at $z = 6–10$ in a single cluster lens field observed with JWST/NIRCam to a depth of 28.5 AB mag. The five differently coloured lines represent different options for the shape of the Pop III IMF throughout the 1–500 $M_{\odot}$ range. All models adopt the Liu & Bromm (2020) Pop III SFRD throughout the $z = 6–10$ interval, and Yoon et al. (2012) stellar evolutionary tracks for non-rotating Pop III stars. Detections of lensed Pop III stars in JWST surveys covering a few tens of cluster-lens fields are only realistic for $log_{10} N(m \le M) \gtrsim -2$, which is only possible for the most massive Pop III stars ($M \gtrsim 300 M_{\odot}$) and the most top-heavy IMFs (purple, blue, and green lines). (b) The redshift probability distribution for the stars that would make it above the detection threshold in the case of the $dN/\log M = \text{constant}$, 1–500 $M_{\odot}$ IMF (but the distribution is similar for all IMFs used in panel a). The distribution is dominated by detections at $z \approx 6–7$, since many of the lowest-$T_{\text{eff}}$ stars have been redshifted out of the JWST/NIRCam wavelength range at higher redshifts.

between adjacent tracks, except for the lowest-mass track in the Yoon et al. (2012) set ($M = 10 M_{\odot}$), for which the assigned interval is 7.5–12.5 $M_{\odot}$, and the highest-mass track considered ($M = 500 M_{\odot}$), for which the assigned interval is 400–500 $M_{\odot}$.

In the following, we consider the $\mu_{\text{detection}}(z) = 100–5 \times 10^4$ magnification range and adopt a microlensing size–magnification limit with $A_0 = 5 \times 10^4$. The latter limit prevents stars in evolutionary stages where $\mu_{\text{detection}}(z)$ from equation (3) exceeds $\max(\mu_{\text{micro}})$ from equation (1) to enter the summation when computing $N_{\text{detection}}(M)$ using equation (7). However, we note that this only affects the detectability of Yoon et al. (2012) Pop III stars with masses $< 50 M_{\odot}$, since the more massive stars can reach the 28.5 AB mag detection limit with much lower magnifications than constrained by this limit. In fact, a $500 M_{\odot}$ star at $z = 6$ can be rendered detectable at a magnification as low as $\mu \approx 100$ in its brightest phase. At magnifications that low, there may admittedly be observational ambiguities between individual lensed stars and lensed star clusters, as many stars within the same star cluster may experience similar macro-magnifications. Hence, our resulting estimate may include cases where the light from several lensed Pop III stars in practice would blend together.

If we take $log_{10} N_{\text{detection}} \gtrsim -2$ as the threshold above which the detection of lensed Pop III stars can be assumed to be realistic (giving a $\gtrsim 10$ per cent detection probability for a survey of $\gtrsim 10$ clusters fields), Fig. 6(a) suggests that it is only the most massive stars ($\gtrsim 300 M_{\odot}$) that have any realistic chance of being detected, and then only for the most top-heavy IMFs considered (Tumlinson 2006 $M_c = 30 M_{\odot}$ and 60 $M_{\odot}$, plus the Magg et al. 2016 $dN/\log M = \text{constant}$ IMF), where the mass fraction in $\gtrsim 100 M_{\odot}$ stars is $\gtrsim 10$ per cent.

The most top-heavy IMF considered here ($dN/\log M = \text{constant}$) produces a total probability for the detection of a Pop III star of $P_{\text{detect}} \approx N_{\text{detection}}(m \le 500 M_{\odot}) \approx 10^{-1.4}$. In this scenario, a survey covering $N_{\text{clusters}}$ cluster-lensing fields would result in a combined detection probability $P_{\text{detect},N_{\text{clusters}}} = 1 - (1 - P_{\text{detect}})^{N_{\text{clusters}}}$ across

the whole survey. This implies $P_{\text{detect},N_{\text{clusters}}} \approx 33$ per cent for a survey covering 10 cluster fields, $\approx 70$ per cent for 30 cluster fields and $\approx 87$ per cent for 50 clusters.

$P_{\text{detect}}$ is lower by $\approx 0.3, 0.6$, and $1.3$ dex, respectively, for the Tumlinson (2006) $M_c = 60 M_{\odot}$, 30, and $10 M_{\odot}$, $\sigma = 1$ IMFs and $\approx 1.1$ dex lower for the Salpeter-slope 1–500 $M_{\odot}$ IMF (please note that the Tumlinson 2006 IMF with $M_c = 10 M_{\odot}$ and $\sigma = 1$ actually predicts slightly fewer $\gtrsim 300 M_{\odot}$ stars compared to the Salpeter-slope IMF extending to 500 $M_{\odot}$). Even in the most pessimistic case ($M_c = 10 M_{\odot}$ and $\sigma = 1$), $P_{\text{detect},N_{\text{clusters}}}$ becomes $\approx 10$ per cent for an $N_{\text{clusters}} = 50$ survey.

Although the Liu & Bromm (2020) pop III SFRD peaks at $z \approx 10$, we in Fig. 6(b) show that the predicted redshift distribution of lensed stars still peaks at $z \approx 6$ (the lowest redshift considered here) and decreases rapidly towards $z=10$ (the highest redshift considered), primarily due to the progressive loss of the SED peak of the lowest-$T_{\text{eff}}$ stars (which, we have argued, are often the easiest to detect) from the NIRCam wavelength window as the redshift is increased. For the Yoon et al. (2012) set of tracks, this redshift evolution in detection probability only show variations at the level of a few percent among the IMFs considered.

## 4.1 Impact of the Pop III star formation rate density evolution

The redshift evolution of the SFRD of Pop III stars in not well-constrained, and different simulations give rise to radically different SFRD($z$) scenarios – some attain SFRDs in the $z \approx 6–10$ range lower than that of the Liu & Bromm (2020) model ($\sim 10^{-4} M_{\odot} \text{cMpc}^{-3} \text{yr}^{-1}$) adopted in Fig. 6, and some attain higher values (see Klessen & Glover 2023 and Venditti et al. 2023 for recent compilations of simulation results).

To approximately estimate the Pop III IMF and SFRD($z$) limits that would render lensed Pop III stars detectable, we for simplicity here treat the Pop III IMF and the Pop III SFRD($z$) as independent.
Naively, one would expect the Pop III IMF to have some impact on the SFRD evolution via feedback and chemical enrichment, but simulations differ significantly in their predictions on this connection. For example, Maio et al. (2010) and Sarmento, Scannapieco & Côté (2019) find that a high characteristic mass of Pop III masses leads to lower Pop III SFRD, whereas Mebane et al. (2018) find the opposite trend, and Pallottini et al. (2014) find only very small variations in SFRD with different Pop III IMF prescriptions.

In general, keeping the Pop III IMF fixed but varying the SFRD results in a situation in which the $P_{\text{detect}}$ estimates roughly scale with Pop III SFRD, so that a version of the Liu & Bromm (2020) SFRD scaled up or down by some factor would result in an updated probability $P_{\text{detect}}$ times this factor. Adopting a constant SFRD($z$) at $3 \times 10^{-4} \text{ M}_\odot \text{ cMpc}^{-3} \text{ yr}^{-1}$ throughout the $z=6-10$ interval, similar to the Pop III SFRD of Venditti et al. (2023), leads to a redshift distribution similar to that of Fig. 6(b), but a $P_{\text{detect}}$ a factor of $\approx 6$ higher than in the Liu & Bromm (2020) case, thereby making even the detection of stars obeying the $M_c = 10 \text{ M}_\odot$ and $\sigma = 1$ IMF possible ($P_{\text{detect,Nclusters}} \approx 40$ per cent) in a 50-cluster survey.

On the other hand, if the $z=6-10$ Pop III SFRD is much lower ($10^{-5} \text{ M}_\odot \text{ cMpc}^{-3} \text{ yr}^{-1}$), i.e. similar to that derived from the Renaissance void simulations of Xu et al. (2016), then even the most-top heavy IMF we consider here ($dN/d \log M = \text{constant}$) would give only $P_{\text{detect,clusters}} \approx 30$ per cent for a 50-cluster survey.

We note that, if stronger constraints on the Pop III IMF were available, it may become possible to use the number of JWST detections of any lensed stars (regardless of their Pop III nature) at $z \gtrsim 6$ to set upper limits on the Pop III SFRD at these redshifts, as some SFRDs could predict a higher rate that that observed. Such Pop III IMF constraints may for instance come from the study of second-generation stars in the local Universe (e.g. Koutsouridou, Salvadori & Skúladóttir 2024).

## 4.2 Impact of late stages of stellar evolution

As we have argued in Section 3, the post-main-sequence evolution of stars has a pronounced impact on the prospects of detecting lensed stars. Since not all Pop III stellar evolutionary tracks show the low-$T_{\text{eff}}$ extension of the Yoon et al. (2012) set, we explore the impact of removing $T_{\text{eff}} < 15000$ K stars from the calculation. This has the effect of lowering the detection prospects of Fig. 6(a) by more than order of magnitude. Hence, the chances of detecting lensed III stars hinges crucially not just on the IMF, the SFRD but also on the evolutionary paths of these stars in the HR diagram. Because the Volpato et al. (2023) stars also evolve to $T_{\text{eff}} < 15000$ K, we find that the detection prospects at $100-500 \text{ M}_\odot$ are very similar to those of the Yoon et al. (2012) tracks without rotation. For Pop III evolutionary tracks like those of Windhorst et al. (2018), which at $10-1000 \text{ M}_\odot$, never evolve to $T_{\text{eff}} < 30000$ K, and those of Murphy et al. (2021), where non-rotating stars in the $10-120 \text{ M}_\odot$ range never reach $T_{\text{eff}} < 15000$ K, the detection prospects are considerably worse.

In Fig. 7, we illustrate this by the exploring the brightest JWST/NIRCam broad-band magnitudes attained at $z=6$ as a function of age for these four sets of tracks at $100-120 \text{ M}_\odot$, and also show the effects of rotation on the Yoon et al. (2012) and Murphy et al. (2021) versions of these stars. The Windhorst et al. (2018) and the Murphy et al. (2021) Pop III stars attain their peak brightness at levels $\approx 1-1.5$ magnitudes fainter than those of the Volpato et al. (2023) and non-rotating Yoon et al. (2012) stars and hence require magnifications $\approx 2.5-4$ times higher to be detected (with corresponding lensing probabilities a factor of $\approx 6-16$ lower). The effect of rotation on these results are small for the Murphy et al. (2021) stars at this mass.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Peak brightness in JWST/NIRCam wide-band filters at $z=6$ versus age for $100-120 \text{ M}_\odot$ stars based on the stellar evolutionary tracks of Yoon et al. (2012), Windhorst et al. (2018), Murphy et al. (2021), and Volpato et al. (2023). Circles mark the end points along each track. Because the Volpato et al. (2023, red line) and non-rotating Yoon et al. (2012, orange line) models for these stars attain much lower temperatures at late stages of evolution ($\approx 4000-5000 \text{ K}$ at their end points) than the other models, they become significantly brighter in the NIRCam filters and therefore more easily detectable. The rotating $v_K = 0.6$ (Yoon et al. 2012, black line) star evolves in the opposite direction, towards increasing $T_{\text{eff}}$, as it ages and therefore attains a considerably lower peak flux in the JWST/NIRCam bands at an age above $\approx 1$ Myr.</figcaption></figure>

However, the Yoon et al. (2012) $100 \text{ M}_\odot$ models with $v_K \ge 0.3$ behave entirely differently, as these have lifetimes extended by up to $\approx 1$ Myr and evolve towards higher $T_{\text{eff}}$ when they age, eventually reaching $\approx 250000$ K. This makes these model stars significantly fainter than the others in the NIRCam filters at $z=6$ after the first million years of evolution.

## 4.3 Impact of the survey detection threshold

The use of deeper JWST exposures than assumed in Fig. 6 (28.5 AB mag) increases the probability for detecting lensed Pop III stars per cluster field, thereby allowing for deeper surveys covering a smaller number of clusters to also reach substantial probabilities for detections of Pop III stars. For example, pushing the survey depth to 30.0 AB mag allows the two top-heaviest IMFs to reach $P_{\text{detect}} \gtrsim 10$ per cent per cluster field for $\gtrsim 300 \text{ M}_\odot$ Pop III stars. In this case, the calculation does however become significantly affected by the requirement that $\mu \ge 100$, since Pop III stars at $\gtrsim 150 \text{ M}_\odot$ may in fact reach the detection limit at even lower magnification. It should be noted, however, that lensed Pop III stars that appear at brightness levels as faint as $\approx 29-30$ AB mag would be extremely difficult to characterize in detail using the photometric SED, and perhaps impossible to study spectroscopically throughout the foreseeable future (Section 5). For reference, the PEARLS program reached detection limits of $\approx 28.5-29.0$ AB mag across several cluster-lens fields using the 7 JWST/NIRCam filters considered here with a total exposure time of 3–4 h (Windhorst et al. 2023). A corresponding programme that reaches 30.0–30.5 AB mag would require $\approx 60$ h.

While surveys covering several filters are required for a photometric characterization of the SED, some filters are more efficient in catching lensed, high-$z$ Pop III stars than others. Given stellar evolutionary tracks like those of Yoon et al. (2012) and Volpato et al. (2023), which predict evolution to $T_{\text{eff}} \lesssim 6000$ K, the F448W filter
<figure><img src="image_10.png" /></figure>

**Figure 8.** NIRCam colour differences between Pop III stars at different temperatures in JWST wide-band filters at $z = 6$ (left) and $z = 10$ (right). The photometric fluxes from Muspelheim SEDs for Pop III stars at $T_{\text{eff}} = 4 \times 10^{3}-2 \times 10^{5}$ K have here been scaled to the same F444W flux. At each $T_{\text{eff}}$, the SEDs are based on the lowest log(g) allowed by the stellar atmosphere grid. The vertical dashed line indicates the position of the Ly $\alpha$ break, to the left of which no detections are possible due to absorption in the neutral intergalactic medium. As seen, colours differences are small at high $T_{\text{eff}}$ but become progressively larger at low $T_{\text{eff}}$. We stress that in the case of the reddest (lowest $T_{\text{eff}}$) SEDs, detections in more than a couple of NIRCam filter are unlikely, since this would require the peak flux to appear many magnitudes above the survey detection limit.

provides the largest fraction of detections, since this is where these low-$T_{\text{eff}}$ stars are the brightest in the $z = 6-10$ range. However, the filter that provides the second-largest fraction of detections gradually shifts across the NIRCam wavelength range from F356W at $z = 6$ to F150W at $z = 10$, as an increasing fraction of the lowest-$T_{\text{eff}}$ stars shift out of NIRCam detection range at the higher redshifts, and higher-$T_{\text{eff}}$ stars with bluer SEDs start to contribute more to the total detection rate, similar to what is seen in Figs 3 (c) and (d) for metal-enriched stars.

# 5 THE OBSERVATIONAL SIGNATURES OF LENSED POPULATION III STARS

## 5.1 Photospheric signatures

In the case where the light from a lensed, high-redshift star is dominated by its photosphere, its JWST broad-band-filter SEDs is expected to contain relatively little information on the nature of the object. While the rest-frame shape of the Muspelheim SEDs do show some dependence on $T_{\text{eff}}$, log(g) and chemical composition$^6$, the effects of $T_{\text{eff}}$ are the dominant ones. Once $T_{\text{eff}}$ and redshift has been estimated from the observed shape of a lensed-star SED, the observed flux levels of the SED can be converted into a constraint on $\mu L_{\text{bol}}$, which will allow the object to be placed into the HR diagram, although with substantial error bars coming from uncertainties on $\mu$, and potentially also from $T_{\text{eff}}$ (see below).

The level at which $T_{\text{eff}}$ can be determined from the SED of a lensed star at high redshifts depends on its intrinsic position in the HR diagram. In Fig. 8, we demonstrate how the JWST/NIRCam colours depend on $T_{\text{eff}}$ in the case of lensed stars at $z = 6$ and $z = 10$. In general, the prospects of setting meaningful constraints are optimized when $T_{\text{eff}}$ is in the range $\approx 8000-20000$ K (F-type

to B-type), since this is the range where the SED develops a prominent Balmer/4000 Å break – a feature that cannot easily be mimicked by dust attenuation along the line of sight. At higher temperatures, the coarse spectrum traced by the JWST broad-band fluxes appears smooth and blackbody-like. The overall slope of this spectrum becomes increasingly steep in the blueward direction as $T_{\text{eff}}$ is increased, but this slope evolution saturates at $T_{\text{eff}} \gtrsim 50000$ K which makes it very challenging to determine higher temperature accurately. Any dust attenuation that remains uncorrected for will also act to make the SED seem to have a lower $T_{\text{eff}}$. In the absence of other dust constraints, the $T_{\text{eff}}$ derived from a featureless blue slope of a lensed star therefore merely serves as a lower limit on the actual $T_{\text{eff}}$ of the object. At $T_{\text{eff}} < 8000$ K, the SED once again turns almost featureless, with a pronounced red peak towards the longest NIRCam wavelengths. In this case, it becomes unrealistic to determine the redshift from the Ly $\alpha$ limit of the star itself, since the flux directly longward of this edge will be much too low to be detected. Provided that the redshift can be obtained from the spectrum or photometric SED of a gravitational arc to which the star belongs, or from the location of the critical curve in the cluster, the temperature may still be fairly accurately constrained at $\approx 6$. At $z \approx 10$, this would require a situation where the star is magnified to a brightness several magnitudes above the detection level in the F444W filter, since this will otherwise be the only NIRCam filter (among those considered here) in which a $T_{\text{eff}} < 8000$ K star can be detected.

The analysis of a $T_{\text{eff}} < 8000$ K star is less sensitive to interstellar dust attenuation and reddening along the line of sight than for a high-$T_{\text{eff}}$ star, since the peak brightness falls in the rest-frame optical and not the rest-frame ultraviolet. However, circumstellar dust produced by the stars themselves may still complicate the analysis.

While metallicity does have a slight impact on the shape of photometric SEDs of lensed, high-redshift stars in the rest-frame UV where metal absorption bands can cause a slight reduction of flux longward of the Ly $\alpha$, these changes are typically small and difficult to distinguish from effects of dust attenuation (further discussed in

<footnote>$^6$Please note that in the present version of Muspelheim, the effects of rotation and mass-loss on the model spectra, beyond those imprinted in the $L_{\text{bol}}$, $T_{\text{eff}}$, log(g) evolution, are ignored.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulated JWST/NIRSpec spectrum around the CII 1335 feature of a $z = 8$, $T_{\text{eff}} = 15000$ K, log(g) = 1.75 star gravitationally magnified to 26.5 AB mag in the JWST/NIRCam F150W filter. The differently lines correspond to stellar atmosphere spectra at different metallicities from our metal-free (Pop III) TLUSTY grid (black line) and the TLUSTY B-star grid by Lanz & Hubeny (2007) at $Z = 0.1 Z_{\odot}$ (blue) and $Z = Z_{\odot}$ (red). (a) Noise-free, intrinsic spectrum at the resolution of the JWST/NIRSpec G140M disperser and the F100LP filter. (b) Mock spectrum with noise level corresponding to 50 h of single-slit JWST/NIRSpec G140M/F100LP observations based on v.2.2 of the JWST exposure time calculator, after rebinning over 2 adjacent wavelength bins to reduce the noise. Even these very deep spectroscopic observations of an exceptionally bright $z = 8$ lensed star only allow for a low-confidence detection of the C II 1335 line (one of the strongest UV lines at 15 000 K; Leitherer et al. 2010) and weak metallicity constraints (at best $Z \lesssim 0.1 Z_{\odot}$ in the case of a non-detection of the C II 1335 feature).</figcaption></figure>

Section 6). The alternative is to attempt to constrain the metallicity of lensed stars through spectroscopy.

Stars with $T_{\text{eff}}$ around 15 000 K are particularly interesting in this context, both because candidate Pop III stars may potentially be detected at these temperatures (Fig. 5; Section 4), because such star remain relatively bright in the part of the rest-frame UV shortward of the Balmer break and because metals imprint a number of absorption features at rest-frame wavelengths $\approx 1216-2000$ Å that can be probed both by JWST/NIRSpec (in principle up to $z \approx 25$) and by future near-IR spectrographs on ground-based ELTs (up to $z \approx 11$, but potentially with higher sensitivity than JWST).

However, absorption line studies of lensed, high-redshift stars are very challenging in the JWST era, due to relatively low continuum signal-to-noise ratios that one can achieve even for exceptionally bright lensed stars.

Some of the most interesting UV absorption features, in the sense of being relatively strong and hence potentially detectable even at lower continuum signal-to-noise ratios ($S/N$) than usually required for absorption-line studies ($S/N \gtrsim 30$), include the Si II 1260 Å line, the C II 1335 Å line, the Si IV 1393, 1403 Å lines, and the Fe III 1893 Å for a supergiant at $T_{\text{eff}} \approx 15000$ K (Leitherer et al. 2010). In Fig. 9, we show mock JWST/NIRSpec, $R \sim 1000$ spectra corresponding to 50 h observations centered on the C II 1335 Å feature in the case of a $z = 8$, $T_{\text{eff}} = 15000$ K star lensed to 26.5 mag in the JWST/NIRCam F150W filter (this gives a peak brightness in the NIRCam bands of $\approx 26$ AB mag, which is exceptionally bright; $\approx 1$ mag brighter than the $z \approx 6$ lensed star Earendel, Welch et al. 2022b). Even for these optimistic conditions (super-deep observations of an unusually bright $z = 8$ star), metal lines can only be measured at low confidence. Hence, demonstrating the absence of such lines, which would indicate a primordial chemical composition, will be very difficult. Based on the mock observations presented in Fig. 9, it would seem possible to constrain the metallicity to $Z < 0.1 Z_{\odot}$ based on the non-detection of the C II 1335 feature, but the distinction of a Pop III star from a metal-poor star (e.g. $Z \sim 0.01 Z_{\odot}$)

would remain difficult in the JWST era. These mock observations are based on TLUSTY spectra which neglect the effect of winds on the UV lines. However, a similar calculation for the CII 1335 features using the PoWR models with various levels of mass-loss (Hainich et al. 2019) shows that while winds may alter the profile of the line and present a more complicated metallicity dependence with line strength, this does not make the line more easily detectable in any radical way. To improve the constraints, data on several metal absorption features would need to be combined. Alternatively, more sensitive spectroscopic data obtained with a future ground-based ELT in the near-IR may also improve the situation.

## 5.2 Nebular signatures

Pop III stars with $T_{\text{eff}} \gtrsim 30000$ K produce enough ionizing photons to produce an emission nebula in the surrounding gas, which in principle may be intrinsically brighter than the star itself at the rest wavelengths probed by JWST at high redshifts (e.g. Rydberg et al. 2013). This nebula is formed as the surrounding gas absorbs the Lyman continuum radiation at rest wavelengths $< 912$ Å and transforms it into nebular continuum and emission lines at $> 912$ Å, with a brightness that strongly depends on $T_{\text{eff}}$ – for a $T_{\text{eff}} = 30000$ K Pop III star, the Lyman continuum makes up $\lesssim 10$ per cent of the bolometric luminosity of the star, but this fraction rises to $> 90$ per cent at $T_{\text{eff}} = 100000$ K.

Under the assumption that an ionization-bounded nebula is formed close to the star, we in Fig. 10 plot the brightest intrinsic AB magnitude attained in the JWST/NIRCam wide-band filters for 60–1000 $M_{\odot}$, non-rotating Pop III stars from the Yoon et al. (2012) set at an age of $\sim 10^5$ yr, along with the brightest JWST/NIRCam AB magnitudes of their associated nebulae, at $z = 6, 10, \text{and } 15$. Here, we have adopted the Pop III nebular spectrum of a 120 $M_{\odot}$ Pop III ZAMS star derived by Rydberg et al. (2013) for the ZAMS parameters from Schaerer (2002), and scaled it to the LyC luminosities of Pop III stars from our sets. As seen in Fig. 10, Pop III nebulae may intrinsically
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The brightest intrinsic (unlensed) JWST/NIRCam wide-band filter fluxes attained by the Yoon et al. (2012) non-rotating Pop III stars (blue lines) and their surrounding H II-regions (red lines) as a function of ZAMS mass, at an age of $10^5$ yr. Solid, dashed, and dotted lines represent redshifts $z = 6, 10, \text{and } 15$, respectively. As seen, unlensed Pop III nebulae intrinsically reach peak JWST/NIRCam fluxes up to $\approx 2-3$ mag brighter than those of their Pop III stars at these redshifts. Please note that the stellar and nebular SEDs typically do not attain their brightest magnitudes in the same NIRCam filters, since the shapes of these two SEDs differ substantially.</figcaption></figure>

be up to $\approx 3$ mag brighter than the stars that power them in the NIRCam broad-band filters considered here. Even so, any potential impact of this nebular emission on observations of lensed high-z stars critically depends on the differential magnification across the star and its nebula, as well as the time-scale over which the nebula can be expected to remain compact – issues that Rydberg et al. (2013) did not address in detail.

### 5.2.1 Differential magnification of Pop III stars and their HII-regions

Because H II-regions are much larger (pc-sized) than stars ($\lesssim 100 R_\odot$, i.e. $\lesssim 5 \times 10^{-4}$ pc, for $\lesssim 1000 M_\odot$ Pop III stars at $T_{\text{eff}} \gtrsim 30000$ K; Windhorst et al. 2018), they cannot attain as large macrolensing magnifications and can furthermore not be significantly magnified by stellar microlensing. Under the simplifying assumption of a constant–density H II-region in which ionization is balanced by recombination, the size of the nebula may be estimated from the Strömgren radius ($R_{\text{str}}$):

$$R_{\text{str}} = \left( \frac{3Q_{\text{tot}}}{4\pi\alpha_{\text{B}}n_{\text{H}}^2} \right)^{1/3}, \qquad (8)$$

where $Q_{\text{tot}}$ is the production rate of LyC photons per second, $\alpha_{\text{B}}$ is the case B recombination rate and $n_{\text{H}}$ is the hydrogen number density.

If one adopts an ISM density close to the Pop III star of $n_{\text{(H)}} = 1000 \text{ cm}^{-3}$ (often considered a standard value for O star nebulae), and the initial $Q_{\text{tot}}$ values for non-rotating $60-1000 M_\odot$ from Yoon et al. (2012), the radii of the Strömgren sphere around these stars become $\approx 1-4$ pc.

Here, we approximate the maximum macrolensing magnification as function of source size as:

$$\max(\mu_{\text{macro}}) \approx 90 B_0 R_{S,\text{pc}}^{-1/2} (D_S/1\text{Gpc})^{1/2}, \qquad (9)$$

derived from equation (2) at the distance from the caustic at which a circular source attains its maximum magnification ($0.65 R_\odot$; Miralda-Encade 1991).

This limits the maximum magnification for nebulae with radii $\approx 1-4$ pc at $z = 6-20$ to $\max(\mu_{\text{macro}}) \approx 500-1500$ in the case of $B_0 = 15$, i.e. more than an order of magnitude less than what lensed stars may attain. For a JWST/NIRCam photometric detection limit of $\approx 29$ AB mag, the faintest Pop III nebula detectable on its own would (without the additional boost of the star) therefore need to reach an intrinsic brightness $\lesssim 35.7-36.9$ AB mag.

From Fig. 10, we see that this limit is, at $z \approx 6$, only met by the nebula of the most massive Pop III stars ($\gtrsim 150 M_\odot$). At $z > 6$, the Pop III nebulae quickly become fainter because the brightest spectral feature, the H $\alpha$ emission line, redshifts out of JWST/NIRCam range, and at $z = 10$, no lensed nebulae of $M \lesssim 500 M_\odot$ Pop III stars would be detectable on their own in a survey reaching 29 AB mag.

In Fig. 11(a), we show the intrinsic SEDs of a 0.1 Myr old, 500 $M_\odot$ Pop III star from Yoon et al. (2012), and its associated nebula, at $z = 6$. In this case, the nebula is intrinsically already $\approx 2$ mag brighter than the star at $\approx 1$ µm, and progressively becomes more dominant at longer wavelengths due to the different slopes of the SEDs.

In Fig. 11(b), we show the differentially magnified SED of a $z = 6, 0.1$ Myr old, 500 $M_\odot$ star at $\mu = 6000$ and its nebula at $\mu = 300$. In this case, parts of the JWST/NIRCam SED of the compound object can be detected above an $\approx 29$ AB mag detection threshold at both $\approx 1-2$ µm and at $\approx 4.4$ µm, with the overall shape of the SED having significant contributions from both the star and its nebula (including a notable upturn in F444W due to the nebular H $\alpha$ emission line). If such situations were common, this would have interesting consequences for the prospects of constraining the metallicity of a lensed star, since strong emission lines from hydrogen and helium could be imprinted on the spectrum, whereas metal lines may be absent.

To illustrate this, we in Fig. 12, show the simulated JWST/NIRSpec R=1000 spectrum that corresponds to the differentially magnified Pop III star + nebula featured in Fig. 11(b). Both the H $\alpha$ and H $\beta$ emission lines are clearly visible, whereas the typically prominent [O III]4959,5007 lines are absent. Upper limits on these oxygen lines should make it possible to set considerably stronger upper limits on the metallicity compared to the case for absorption line spectroscopy of the Pop III photosphere. We note that it would be considerably more challenging to detect the He II 1640 emission line, which is often highlighted as a strong Pop III indicator, due to its relatively low intrinsic flux.

### 5.2.2 Temporal evolution of the H II-regions around Pop III stars

While the analysis presented in Section 5.2.1 suggests that nebular features in principle could be imprinted in the spectra of lensed $M \gtrsim 150 M_\odot$ Pop III stars, the probability of detecting such features critically depends on the time-scales over which the H II-region around Pop III stars can remain sufficiently compact (pc-scale) to allow magnification factors of several hundreds or more.

Even though the gas densities required to keep the Strömgren sphere contained at parsec scales ($n_H \gtrsim 1000 \text{ cm}^{-3}$) are likely to be retained in the vicinity of the star shortly after formation, the hot photoionized H II region that forms around the star will be overpressured compared to the neutral, cooler surroundings and expand. Models for the density evolution of the ionized gas in the H II region around both Pop III stars in $\sim 10^6 M_\odot$ minihaloes (e.g.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Spectral energy distributions of Pop III nebulae. (a) Unlensed spectral energy distributions of a 0.1 Myr old, 500 $M_{\odot}$ non-rotating Yoon et al. (2012) Pop III star (red line) and its surrounding H II-region (blue line) in JWST/NIRCam and MIRI wide filters out to 10 µm at $z = 6$, based on the Rydberg et al. (2013) models for the nebulae of Pop III stars. The black solid line indicates the underlying spectra on which the photometric SEDs of the star and its H II-region are based. The dashed vertical line marks the position of the Ly α limit. (b) The gravitationally lensed versions of the SEDs shown in (a), assuming $\mu = 6000$ for the star (blue line) and $\mu = 300$ for its H II-region (red line). The green line represents the SED of the compound object (differentially magnified star + nebula). At these magnifications, either the lensed star, the lensed nebula, or the differentially magnified compound object could potentially be detected by JWST in deep exposures.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Simulated JWST/NIRSpec G395M/F290L fixed-slit spectrum at resolution $R \approx 1000$ after 20 h of observations of the $z = 6, 500 M_{\odot}, 0.1$ Myr old differentially magnified star ($\mu = 6000$) plus nebula ($\mu = 300$) from Fig. 11(b). The H α and H β lines from the nebula are clearly detectable (at peak $S/N \approx 18$ and 7, respectively), whereas the absence of the emission lines [O III]4959,5007, indicates a very low metallicity.</figcaption></figure>

Whalen, Abel & Norman 2004; Alvarez, Bromm & Shapiro 2006; Sibony et al. 2022) and high-mass, metal-enriched stars in more normal H II regions (e.g. Mellema et al. 2006; Geen et al. 2021) all suggest that the H II region may expand beyond $> 1$ pc already at $\lesssim 10^5$ yr, which makes the time-scale for detection perilously short in the lensing scenario considered here.

In Appendix C, we present a simple toy calculation that suggests that if the H II-region around Pop III stars indeed only remains compact (pc-sized) for $\approx 10^5$ yr, then the probability for such H II-regions to be magnified to $\approx 29$ AB mag (as in Figs 11b and 12) is $N_{\text{detect}} < 10^{-3}$ per cluster-lens field, which is too low for such effects to realistically be detected with JWST. However, the situation could

potentially improve if the expansion of the H II-region is hampered by the gravitational field of a more massive halo (Kitayama et al. 2004), as could be the case if Pop III star formation largely shifts to $\gtrsim 10^7-10^8 M_{\odot}$ H I-cooling haloes at $z \lesssim 15$, or by the gravitational field of the star itself (Jaura et al. 2022). Moreover, in the case of significant mass loss from a Pop III star (e.g. due to rotation), the wind/outflow may also imprint emission lines on the spectrum of a lensed star.

# 6 DISCUSSION

## 6.1 Comparison with other studies

Rydberg et al. (2013) explored the prospects of detecting highly magnified high-redshift Pop III stars with very deep JWST imaging of a single cluster-lens field and came to the conclusion that detections would be highly unlikely, unless essentially all Pop III stars formed at very high masses $M \gtrsim 300 M_{\odot}$ (and only when considering a maximal flux boost due to nebular emission). However, Rydberg et al. (2013) did not consider Pop III stars beyond the main sequence and therefore did not include the substantial boost in detection rates provided by potential low-$T_{\text{eff}}$ states of evolution. Also, the range of Pop III SFRDs explored at $z = 6-10$ was also significantly lower than that predicted by some recent simulations and models (e.g. Jaacks, Finkelstein & Bromm 2019; Liu & Bromm 2020; Chantavat, Chongchitnan & Silk 2023; Venditti et al. 2023).

The prospects of detecting lensed Pop III stars was revisited by Windhorst et al. (2018), who did consider post-main sequence of Pop III stars, although not evolution to $T_{\text{eff}}$ states as low as those explored in the current paper. Based on the assumption that the cosmic infrared background receives a non-negligible contribution from high-redshift Pop III stars, Windhorst et al. (2018) argue that there is significant chance of detecting highly magnified Pop III stars or, alternatively, the accretion discs around the black holes they leave behind, through caustic-crossing events with JWST if a handful of strong-lensing galaxy clusters are monitored over a few years.
However, in scenarios where the contribution from these objects lie far below the current upper limits on the diffuse infrared background, the *JWST* detection prospects for lensed Pop III stars or their remnant black holes become progressively worse.

Recently, Larkin, Gerasimov & Burgasser (2023) presented an extensive study of the magnifications required to detect 1–1000 $M_{\odot}$ Pop III stars at $z = 3–17$ with *JWST*, albeit only for the ZAMS and without assessing the probability for such detections in realistic surveys of lensed fields.

## 6.2 The effect of dust attenuation

The results presented in this paper assume that lensed stars are not significantly affected by dust attenuation. However, since the stars detected through gravitational lensing are all likely to have high-masses and short lifetimes ($\approx 2–30$ Myr) they are likely to be sitting in young, star-forming regions, where the interstellar dust attenuation may be substantial, especially in the case of high-metallicity environments. Moreover, red supergiants may be surrounded by circumstellar dust, as seen in local samples (Massey et al. 2005; Drout, Massey & Meynet 2012). This attenuation of the intrinsic flux of lensed stars has the effect of producing an additional, and likely metallicity-dependent bias affecting where in the HR diagram such stars are likely to be found. To first order, the dust attenuation can potentially be assessed from photometric and/or spectroscopic studies of gravitational arcs in which lensed stars have so far been identified, but it must be kept in mind that the specific star-forming region in which a lensed star is sitting may not be the same as that derived for its host galaxy. In general, the effects of dust attenuation on the observed SED of a lensed star candidate is to render it fainter and redder, thereby – at fixed estimate of the magnification – shifting it to what (in the absence of corrections for dust) would appear to be a lower $T_{\text{eff}}$ and lower $L_{\text{bol}}$ compared to its intrinsic position.

Since Pop III stars are unlikely to form in dusty surroundings, dust reddening, and attenuation is generally expected to be a more severe problem for metal-enriched lensed stars. Even so, Pop III supernovae are expected to produce dust (e.g. Marassi et al. 2015), and dust generated by the most massive and short-lived members of a small cluster of Pop III stars could potentially affect the line of sight to more long-lived Pop III members in the same region.

## 6.3 Rotation of Pop III stars

A significant fraction of Pop III stars may experience rapid rotation (Hirano & Bromm 2018), which is expected to give rise to a number of complicated phenomena such as non-spherical stars, mechanical mass loss and possibly the formation of decretion discs (e.g. Marigo, Chiosi & Kudritzki 2003; Ekström et al. 2008; Yoon et al. 2012; Murphy et al. 2021). Plane-parallel stellar atmospheres in hydrostatic and radiative equilibrium of the type used here are ill-equipped to model such phenomena, and while such stellar atmosphere models are often used to model the non-ionizing continuum of rotating high-mass stars, the problem becomes more complex when modelling ionizing fluxes, since the outflow of gas from the surface could potentially have a significant impact on the emergent ionizing spectrum (e.g. Schaerer 2002).

## 6.4 Binary Pop III stars

Almost all massive stars are in binary systems (see De Marco & Izzard 2017, for a review), and simulations suggest that Pop III stars may be as well (e.g. Stacy & Bromm 2013; Sugimura et al. 2020).

Aside from having a significant impact on the rotation and evolution of Pop III stars, and hence the overall probability of detecting highly magnified Pop III stars, this also opens up the possibility that what is observed in the case of a Pop III binary located close to a caustic, is the blended contributions from two separate Pop III stars with different properties and potentially slightly different magnifications. This becomes another complication in the analysis of the SEDs of lensed-star candidates, although the two components could potentially be disentangled in the case where $T_{\text{eff}}$ of the two components are substantially different (as has been suggested to be the case for the $z \approx 6$ object Earendel; Welch et al. 2022b) or where the radial velocity between the two stars allows for the detection of two sets of absorption lines.

## 6.5 Surface enrichment of Pop III stars

Even though Pop III stars form at primordial chemical composition, heavy elements are produced from fusion within these stars. Dredge-ups during late stages of Pop III stellar evolution (e.g. Ekström et al. 2008; Volpato et al. 2023, 2024), or mixing with fast-rotating Pop III stars (Yoon et al. 2012; Liu et al. 2021) would allow nucleosynthetic products to reach the surface. Since late, low-$T_{\text{eff}}$ stages of evolution are favourable for detection, this leads to an additional problem when attempting to verify the Pop III nature stars based on spectroscopy. Absorption lines from the stellar surface, or a combination of absorption and emission lines from a wind/outflow may in such cases reveal non-zero CNO levels, which would complicate the identification of the lensed star as Pop III object.

## 6.6 Unsolved problems in the evolution of high-mass stars

We have throughout this paper stressed that existing evolutionary models differ greatly in their predictions for how Pop III stars move across the HR diagram, and that this has pronounced effects on the detectability of such stars. It must be pointed out, however, that current models for metal-enriched, high-mass stars also suffer from a number of problems, including the lack of observed O-stars close to the predicted ZAMS (e.g. Holgado et al. 2020) and the observed excess of blue supergiants compared to models (e.g. Bellinger et al. 2023). Whatever the resolution to these issues, this is likely to affect the time high-mass stars spend at different $T_{\text{eff}}$ and hence forecasts of where in the HR diagram lensed stars are the most likely to be found (Figs 2 and 3). Until such issues have been resolved, caution must be exercised when attempting to interpret observed distributions of lensed stars across the HR diagram in terms of quantities like the star formation history, metallicity, or the IMF of the lensed stellar population.

# 7 CONCLUSIONS

Our results can be summarized as follows:

(i) Surveys for highly magnified, high-redshift stars with *JWST* are subject to strong selection biases related to the *JWST* wavelength range and detection limits, the lensing probability distribution and size–magnification limits (Section 3). For stars that require the highest magnifications to be rendered detectable, these biases tend to favour the detection of objects in late, lower-$T_{\text{eff}}$ stages of evolution even though stars only spend a small fraction of their lifetime in such states. For example, based on the $Z \sim 0.1 Z_{\odot}$ single-star stellar evolutionary tracks used in this paper, $z \approx 3$ stars with initial masses $\lesssim 55 M_{\odot}$ found by *JWST* at $\lesssim 28$ AB magnitudes are only likely to be detected as red supergiants (Fig. 3b).
(ii) We argue that highly magnified Pop III stars at $z \gtrsim 6$ may be detectable in small numbers with JWST surveys covering a large number ($\approx 10-50$) of cluster-lens fields. However, to make such detections likely, Pop III stars need to obey a number of conditions: The Pop III SFRD at $z \approx 6-10$ needs to be sufficiently high ($\gtrsim 10^{-4}$ M$_\odot$ cMpc$^{-3}$ yr$^{-1}$), the Pop III IMF needs to be very top-heavy and allow a significant mass fraction of the gas converted into Pop III stars to be in $\gtrsim 100$ M$_\odot$ objects (mass fraction $\gtrsim 10$ per cent for SFRD $\sim 10^{-4}$ M$_\odot$ cMpc$^{-3}$ yr$^{-1}$), and the evolution of Pop III stars must be such that the $M \gtrsim 100$ M$_\odot$ stars evolve to temperatures $< 15000$ K during their lifetimes (Section 4; Fig. 6). Several simulations predict Pop III SFRDs up to an order of magnitude higher than our limit (see e.g. Venditti et al. 2023), and the recent stellar evolutionary models by Volpato et al. (2023) do predict evolution into low-$T_{\text{eff}}$ states for $\gtrsim 100$ M$_\odot$ Pop III stars. Whether the required Pop III IMF is likely to be sufficiently top-heavy to meet our constraints is potentially more controversial, since studies of the abundance patterns of very metal-poor stars in the local Universe do not seem to favour the enrichment signatures that would result from the most top-heavy Pop III IMFs (e.g. Tumlinson 2006; Koutsouridou et al. 2023).

(iii) Even if the most massive, lensed Pop III stars at $z \gtrsim 6$ were to be lifted above the detection threshold of JWST while close to their ZAMS states at $T_{\text{eff}} \approx 10^5$ K (Fig. 4a) – pinning down $T_{\text{eff}}$ from JWST photometry will be difficult at $T_{\text{eff}} \gtrsim 50000$ K due to modest change in spectral slope in this $T_{\text{eff}}$ regime (Fig. 8).

(iv) Verifying the Pop III nature of highly magnified stars at high redshifts using metallicity constraints will be very challenging. In situations where the lensed image is dominated by the stellar photosphere, very deep JWST and/or ELT spectroscopy can be used to target strong metal absorption lines to set an upper limit on the metallicity of the candidate. However, the constraint is unlikely to be better than $Z/Z_\odot < 0.1$ given the capabilities of JWST, which makes it difficult to separate a very low-metallicity star from a genuine Pop III star (Section 5.1; Fig. 9). Moreover, Pop III stars may also be self-polluted by metals in the late evolutionary states where many of these stars are likely to be found, thereby complicating the analysis.

(v) During early stages of evolution, where a Pop III star may be surrounded by a compact H II region, the differentially magnified nebula could make the object brighter and imprint strong hydrogen emission lines in the spectrum, without any associated metal emission lines (Fig. 12). While this opens up the possibility to detect Pop III stars in their earliest states and probe the metallicity more easily than what is possible using absorption line spectroscopy, the probability of such detections is deemed very low, unless the expanding H II region around the star is confined to $\sim \text{pc}$ scales for longer than $\sim 10^5$ yr (Section 5.2).

(vi) The fact that individual Pop III stars may be rendered detectable by JWST through gravitational lensing in some models has important consequences even in the absence of metallicity constraints on any lensed stars that are found at $z \gtrsim 6$. As JWST continues to observe an increasing number of cluster-lens fields, it should be possible to use the accumulated photometric detections (or even non-detections) of lensed stars at these redshifts to set combined upper limits on parameters that describe the SFRD, IMF, and evolution of Pop III stars, thereby ruling out the most extreme scenarios.

## ACKNOWLEDGEMENTS

The referee is thanked for helpful comments on the manuscript. EZ acknowledges project grant 2022-03804 from the Swedish Research Council (Vetenskapsrådet) and has also benefited from a sabbatical at the Swedish Collegium for Advanced Study. EZ and

AV acknowledge funding from the Swedish National Space Agency. AN acknowledge funding from the Olle Engkvist Foundation. GV acknowledges support from Padova University, through the research project PRD 2021. GC acknowledges support from the Agence Nationale de la Recherche grant POPSYCLE number ANR-19-CE31-0022. JMD acknowledges the support of project PID2022-138896NB-C51 (MCIU/AEI/MINECO/FEDER, UE) Ministerio de Ciencia, Investigación y Universidades. AZ acknowledges support by grant No. 864/23 from the Israel Science Foundation; grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and grant No. 2109066 from the United States National Science Foundation (NSF); and by the Ministry of Science & Technology, Israel.

## DATA AVAILABILITY

Publicly available grids of Muspelheim models are available from: https://www.astro.uu.se/~ez/muspelheim/muspelheim.html.

## REFERENCES

Alvarez M. A., Bromm V., Shapiro P. R., 2006, ApJ, 639, 621

Bellinger E. P., de Mink S. E., van Rossem W. E., Justham S., 2024, ApJ, 967, L39

Bergemann M., 2014, in Niemczura E., Smalley B., Pych W., eds, Determination of Atmospheric Parameters of B-, A-, F- and G-Type Stars. Springer, Cham, p. 187

Broadhurst T. et al., 2024, preprint (arXiv:2405.19422)

Chantavat T., Chongchitnan S., Silk J., 2023, MNRAS, 522, 3256

Chen W. et al., 2022, ApJ, 940, L54

de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott I., Yoon S. C., 2009, A&A, 497, 243

De Marco O., Izzard R. G., 2017, PASA, 34, e001

Diego J. M., 2019, A&A, 625, A84

Diego J. M. et al., 2018, ApJ, 857, 25

Diego J. M. et al., 2023a, A&A, 672, A3

Diego J. M. et al., 2023b, A&A, 679, A31

Diego J. M. et al., 2024a, preprint (arXiv:2404.08033)

Diego J. M. et al., 2024b, A&A, 681, A124

Drout M. R., Massey P., Meynet G., 2012, ApJ, 750, 97

Ekström S., Meynet G., Chiappini C., Hirschi R., Maeder A., 2008, A&A, 489, 685

Eldridge J. J., Vink J. S., 2006, A&A, 452, 295

Fudamoto Y. et al., 2024, preprint (arXiv:2404.08045)

Furtak L. J. et al., 2024, MNRAS, 527, L7

Geen S., Bieri R., Rosdahl J., de Koter A., 2021, MNRAS, 501, 1352

Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS, 474, 2757

Hainich R., Ramachandran V., Shenar T., Sander A. A. C., Todt H., Gruner D., Oskinova L. M., Hamann W. R., 2019, A&A, 621, A85

Hirano S., Bromm V., 2018, MNRAS, 476, 3964

Holgado G. et al., 2020, A&A, 638, A157

Hubeny I., Lanz T., 1995, ApJ, 439, 875

Inoue A. K., 2011, MNRAS, 415, 2920

Jaacks J., Finkelstein S. L., Bromm V., 2019, MNRAS, 488, 2202

Jaura O., Glover S. C. O., Wollenberg K. M. J., 2004, ApJ, 613, 631

Klessen R. S., Glover S. C. O., 2023, ARA&A, 61, 65

Koutsouridou I., Salvadori S., Skúladóttir Á., 2024, ApJ, 962, L26

Koutsouridou I., Salvadori S., Skúladóttir Á., Rossi M., Vanni I., Pagnini G., 2023, MNRAS, 525, 190

Kroupa P., 2001, MNRAS, 322, 231

Lanz T., Hubeny I., 2003, ApJS, 146, 417
Lanz T., Hubeny I., 2007, ApJS, 169, 83

Larkin M. M., Gerasimov R., Burgasser A. J., 2023, AJ, 165, 2

Leitherer C., Ortiz Otálvaro P. A., Bresolin F., Kudritzki R.-P., Lo Faro B.,
Pauldrach A. W. A., Pettini M., Rix S. A., 2010, ApJS, 189, 309

Lejeune T., Cuisinier F., Buser R., 1997, A&AS, 125, 229

Liu B., Bromm V., 2020, MNRAS, 497, 2839

Liu B., Sibony Y., Meynet G., Bromm V., 2021, MNRAS, 506, 5247

Magg M., Hartwig T., Glover S. C. O., Klessen R. S., Whalen D. J., 2016, MNRAS, 462, 3591

Magg M., Hartwig T., Agarwal B., Frebel A., Glover S. C. O., Griffen B. F., Klessen R. S., 2018, MNRAS, 473, 5308

Maio U., Ciardi B., Dolag K., Tornatore L., Khochfar S., 2010, MNRAS, 407, 1003

Marassi S., Schneider R., Limongi M., Chieffi A., Bocchio M., Bianchi S., 2015, MNRAS, 454, 4250

Marigo P., Chiosi C., Kudritzki R. P., 2003, A&A, 399, 617

Mas-Ribas L., Dijkstra M., Forero-Romero J. E., 2016, ApJ, 833, 65

Massey P., Plez B., Levesque E. M., Olsen K. A. G., Clayton G. C., Josselin
E., 2005, ApJ, 634, 1286

Mebane R. H., Mirocha J., Furlanetto S. R., 2018, MNRAS, 479, 4544

Meena A. K. et al., 2023, ApJ, 944, L6

Mellema G., Arthur S. J., Henney W. J., Iliev I. T., Shapiro P. R., 2006, ApJ, 647, 397

Miralda-Escude J., 1991, ApJ, 379, 94

Murphy L. J. et al., 2021, MNRAS, 501, 2745

Nakajima K., Maiolino R., 2022, MNRAS, 513, 5134

Palencia J. M., Diego J. M., Kavanagh B. J., Martinez J., 2024, A&A, 687, A81

Pallottini A., Ferrara A., Gallerani S., Salvadori S., D'Odorico V., 2014, MNRAS, 440, 2498

Rieke M. J. et al., 2023, PASP, 135, 028001

Rydberg C.-E., Zackrisson E., Lundqvist P., Scott P., 2013, MNRAS, 429, 3658

Sanyal D., Grassitelli L., Langer N., Bestenlehner J. M., 2015, A&A, 580, A20

Sarmento R., Scannapieco E., Côté B., 2019, ApJ, 871, 206

Schaerer D., 2002, A&A, 382, 28

Schaerer D., 2003, A&A, 397, 527

Schauer A. T. P., Drory N., Bromm V., 2020, ApJ, 904, 145

Schauer A. T. P., Bromm V., Drory N., Boylan-Kolchin M., 2022, ApJ, 934, L6

Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses. Springer-
Verlag, Berlin, Heidelberg, New York, p. 319

Sibony Y., Liu B., Simmonds C., Meynet G., Bromm V., 2022, A&A, 666, A199

Song H., Meynet G., Li Z., Peng W., Zhang R., Zhan Q., 2020, ApJ, 892, 41

Stacy A., Bromm V., 2013, MNRAS, 433, 1094

Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ, 892, L14

Surace M. et al., 2018, ApJ, 869, L39

Surace M., Zackrisson E., Whalen D. J., Hartwig T., Glover S. C. O., Woods
T. E., Heger A., Glover S. C. O., 2019, MNRAS, 488, 3995

Szécsi D., Agrawal P., Wünsch R., Langer N., 2022, A&A, 658, A125

Trussler J. A. A. et al., 2023, MNRAS, 525, 5328

Tumlinson J., 2006, ApJ, 641, 1

Tumlinson J., Shull J. M., 2000, ApJ, 528, L65

Tumlinson J., Giroux M. L., Shull J. M., 2001, ApJ, 550, L1

Tumlinson J., Shull J. M., Venkatesan A., 2003, ApJ, 584, 608

Venditti A., Graziani L., Schneider R., Pentericci L., Di Cesare C., Maio U., Omukai K., 2023, MNRAS, 522, 3809

Venumadhav T., Dai L., Miralda-Escudé J., 2017, ApJ, 850, 49

Volpato G., Marigo P., Costa G., Bressan A., Trabucchi M., Girardi L., 2023, ApJ, 944, 40

Volpato G., Marigo P., Costa G., Bressan A., Trabucchi M., Girardi L., Addari
F., 2024, ApJ, 961, 89

Wang C. et al., 2022, Nat. Astron., 6, 480

Welch B. et al., 2022a, Nature, 603, 815

Welch B. et al., 2022b, ApJ, 940, L1

Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14

Windhorst R. A. et al., 2018, ApJS, 234, 41

Windhorst R. A. et al., 2023, AJ, 165, 13

Wu Z., Kravtsov A., 2024, Open J. Astrophys., 7, 56

Xu H., Ahn K., Norman M. L., Wise J. H., O'Shea B. W., 2016, ApJ, 832, L5

Yan Z., Jerabkova T., Kroupa P., 2023, A&A, 670, A151

Yoon S. C., Dierks A., Langer N., 2012, A&A, 542, A113

Zackrisson E., Rydberg C.-E., Schaerer D., Östlin G., Tuli M., 2011, ApJ, 740, 13

Zitrin A. et al., 2013, ApJ, 762, L30

APPENDIX A: STELLAR ATMOSPHERE GRIDS

Here, we provide further information on the technical procedure for
generating SEDs of stars in the Muspelheim models. In Fig. A1,
we show selected stellar evolutionary tracks overlaid on the T<sub>eff</sub>-
log(g) coverage of the stellar atmosphere grids used at Z = 0.1 Z<sub>☉</sub>
(a) and Z = 0 (b). This illustrates a common problem, in that the
most massive stars often venture into a region of parameter space at
very low log(g) where stellar atmosphere models are absent. For the
lowest mass tracks, this happens during a very small fraction of their
lifetimes, but at M ≳ 50 M<sub>☉</sub>, this fraction grows to ≳ 10 per cent,
and for the highest mass tracks, a majority of their evolution occurs
at lower log(g) than formally covered by the grid. The low-log(g)
boundary coarsely traces the Eddington limit for these stars, but since
the grids are largely based on a fixed Δ log(g) step size, the actual
limit occurs at log(g) slightly outside these grids. In the case of the
300 M<sub>☉</sub> track featured in panel A1b, ≳ 90 per cent of the lifetime is
spent within one Δ log(g) grid step of the current grid boundary, and it
is only at the late, short-lived, lower-T<sub>eff</sub> stages where the discrepancy
becomes more substantial. On the other hand, in the case of the 250
M<sub>☉</sub> track in panel a, the majority of the lifetime is spent farther away
than this from the grid. As argued by Sanyal et al. (2015), many
massive stars are likely to straddle the Eddington limit at their surface.
The effects of forcing 1D, hydrostatic atmosphere models on to such
stars, which may also be non-spherical due to rotation, will inevitably
lead to inaccuracies when it comes to individual spectral features.
There is, however, no set of published stellar atmosphere models that
can deal with such effects over the large regions of parameter space
required by sets of stellar evolutionary tracks of the type used here.

In cases where data points along a stellar evolutionary tracks fall
outside the log(g) limits of the grid of stellar atmosphere spectra, an
interpolation between the two closest log(g) grid models is performed
in T<sub>eff</sub>, instead of an attempted extrapolation. In the procedure of
interpolating the stellar atmosphere grids, wavelength interpolations
are also occasionally required due to differences in wavelength
coverage for stellar atmosphere spectra in different parts of these
grids. While this is expected not to have any significant impact
on broad- and medium-band magnitudes, individual absorption or
emission features could get distorted. For this reason, an alternative
mode is also implemented, in which the stellar atmosphere spectra
that lies the closest in T<sub>eff</sub> and log(g) is instead selected, without any
interpolation. While this ensures that the native spectral sampling is
always used and that no spurious spectral features are introduced by
the interpolation procedure, this implies that quantities like colours,
that depend solely on the shape of the SEDs, evolve in a less smooth,
stepwise fashion as as these stars evolve through the T<sub>eff</sub>-log(g) grid.
By comparing the results from the two strategies, we have checked
that the results presented in this paper are not strongly affected by
the wavelength interpolation procedure.
<figure><img src="image_15.png" /></figure>

**Figure A1.** $T_{\text{eff}}$–log(g) diagram showing selected stellar evolutionary tracks for massive stars overlaid on the corresponding grid of stellar evolutionary spectra. (a) $Z \approx 0.002$ (SMC) metallicity tracks from Szécsi et al. (2022) at 9 $M_{\odot}$ (orange line), 30 $M_{\odot}$ (green), and 250 $M_{\odot}$ (blue) overlaid on a $Z \approx 0.1$ $Z_{\odot}$ stellar atmosphere grid. (b) Pop III tracks without rotation from Yoon et al. (2012) at 10 $M_{\odot}$ (orange line), 30 $M_{\odot}$ (green), and 300 $M_{\odot}$ (blue), along with one 30 $M_{\odot}$ track with rotation ($v_{\text{init}}/v_{\text{k}} = 0.4$; red line) overlaid on the $Z \approx 0$ stellar atmosphere grid. Along each evolutionary track in both panels, stars mark the start (youngest age) and squares the end. Grey markers indicate the locations of TLUSTY grid points and black points the positions of Lejeune et al. (1997) grid points. As seen in both panels, the $M = 30 M_{\odot}$ tracks venture outside the boundaries of these grids during certain points during their evolution, and the $M \approx 250–300 M_{\odot}$ tracks spend a large fraction of their evolution below the lowest log(g) boundary of the grids.

## APPENDIX B: RELATION BETWEEN BRIGHTNESS AND EFFECTIVE TEMPERATURE

Here, we provide further explanations for why observations of lensed stars will favour the detection of stars in relatively short-lived low-$T_{\text{eff}}$ states over the more long-lived high-$T_{\text{eff}}$ states associated with the stellar main sequence. In Fig. B1, we show the rest-frame evolution in monochromatic luminosity $L_{\nu}$ and observed-frame evolution of $f_{\nu}$ (expressed in AB magnitudes) of the SED of an artificial star at $z=6$ that is assumed to evolve from $T_{\text{eff}} = 200000$ to 4000 K at constant $L_{\text{bol}} = 10^6 L_{\odot}$. While realistic stars admittedly do not evolve at fixed $L_{\text{bol}}$, this is not far from the truth of $M \gtrsim 100 M_{\odot}$ Pop III stars, which are often predicted to evolve almost horizontally in the HR diagram (Windhorst et al. 2018; Volpato et al. 2023).

Fig. B1(a) indicates a factor of ~100 difference in peak $L_{\nu}$ when $T_{\text{eff}}$ goes from high to low. However, the full wavelength range is not directly observable (at any redshift) due to foreground gas absorption at < 912 Å (neutral ISM) or < 1216 Å (neutral IGM at $z \gtrsim 6$) and the wavelength coverage of the most sensitive JWST instruments also limits the relevant range. Fig. B1(b) shows that there is an ≈ 5.5 magnitude difference at $z = 6$ in peak brightness in the ≈ 1–5 µm window available for JWST NIRCam or NIRSpec observations, between the highest and lowest $T_{\text{eff}}$ stars plotted.

As a toy example of how these effects convert into detection probabilities, consider a scenario where this artificial star only has two significant $T_{\text{eff}}$ stages during its evolution across the HR diagram – one long-lived $10^5$ K state and one short-lived $10^4$ K state. In Fig. B1(b), we see that the latter state could be ≈ 4 mag

<figure><img src="image_16.png" /></figure>

**Figure B1.** The (a) rest-frame monochromatic luminosity $L_{\nu}$ and (b) observed-frame brightness (AB magnitude) at $z = 6$ for a star at different $T_{\text{eff}}$ but fixed bolometric luminosity $L_{\text{bol}} = 10^6 L_{\odot}$. As $T_{\text{eff}}$ decreases, the peak $L_{\nu}$ and observed brightness increase substantially. In panel b, the vertical dashed line marks the $z = 6$ Ly α limit, to the left of which the radiation is expected to be largely absorbed by the neutral IGM. Hence, the part accessible to JWST lies to the right of this line.
brighter in the JWST window at $z = 6$. Since both of these stars would intrinsically lie far below the detection threshold of JWST, a substantial magnification boost is required to bring them into the detectable range. Assuming that the magnification probability distribution follows $P(> \mu) \propto \mu^{-2}$, the magnification required for the $10^4$ K state would be $\approx 40$ times lower than for the $10^5$ K state, which would make it $40^2 = 1600$ times more likely to occur. Hence, even if the star only spends 1 per cent of its lifetime in the short-lived $10^4$ K state, it would still be a factor of $\approx 16$ times more likely to be found in this state rather than in the more long-lived $10^5$ K state.

## APPENDIX C: PROBABILITY OF DETECTING A LENSED POP III NEBULA

Here, we present a rough estimate of the probability to detect a Pop III nebula lensed above the detection limit of a 29 AB mag JWST/NIRCam detection limit, arguing that this probability is no more than $\sim 10^{-3}$ per cluster field if the H II-region remains confined to parsec scales for time-scales of $\sim 10^5$ yr.

In Fig. C1, we use the inverse ray-shooting code of Meena et al. (2023) to plot the $z = 6$ macrolensing magnification for different source sizes as a function of distance from the caustic in the MACS J0416.1-2403 galaxy cluster based on the mass model from Zitrin et al. (2013). The behaviour of circular sources with two different sizes relevant for Pop III nebulae ($R_S = 0.3$ pc and 1 pc), are compared to the point-source approximation (equation (2), with $B_0 \approx 10$ at the chosen position along the critical curve) which, within the range of magnifications plotted, illustrates the magnification curve of a stellar-sized source.

While the magnification continues to increase to $\mu \gtrsim 10^4$ as the stellar-sized source approaches the caustic, the two nebular-sized sources reach their maximum magnifications at a distance of $0.65 R_S$ (Miralda-Escude 1991). Moreover, while the stellar-sized source becomes unobservable at the point of crossing the caustic, the larger sources can retain a relatively high magnification even when the midpoint has passed slightly beyond the caustic.

At a distance $\approx 1.5$ pc from the caustic, where the macrolensing magnification in the situation plotted is $\mu \approx 600$, a young Pop III star would – owing to its lower intrinsic flux (Fig. 10) – typically not be rendered detectable by JWST from macrolensing alone, whereas an ionization-bounded nebula around a $\gtrsim 200$ M$_\odot$ Pop III star at $z = 6$ would, by being boosted to $\gtrsim 29.5$ AB mag (Fig. 10). In the absence of an additional microlensing boost to the flux of the star, this would correspond to a case where the observed SED would be nebular dominated. If the star, on the other hand, is highly microlensed at this position, both components could contribute significantly to the SED.

To provide a rough estimate of the probability to detect the lensed Pop III nebula of a single star, we here consider the very optimistic case where all Pop III stars form with a mass of 500M$_\odot$. For the Yoon et al. (2012) non-rotating tracks, we find that a 500 M$_\odot$ star at age $10^5$ yr would produce an ionizing flux $Q_{tot} \approx 7.4 \times 10^{50}$ s$^{-1}$ and a nebula that reaches a peak JWST/NIRCam flux corresponding to $\approx 35.3$ AB mag at $z = 6$. To reach a NIRCam detection limit of

29 AB mag, we therefore need a magnification of $\mu \approx 330$, which according to equation (9) can be achieved for source sizes with radii up to $\approx 20$ pc at this redshift. To contain the Strömgren radius of the H

<figure><img src="image_17.png" /><figcaption>**Figure C1.** Macrolensing magnification $\mu$ in the absence of microlenses as a function of distance from the cluster caustic (position marked by the dashed vertical line) at $z = 6$ for a point source (blue line, also suitable for the photosphere of a Pop III star within the range of magnifications plotted), and for sizes more characteristic of H II-regions around Pop III regions (red and yellow lines, for 0.3 and 1.0 pc, respectively). The situation shown corresponds to a position along the cluster caustic at which $B_0 \approx 10$ (equation 2).</figcaption></figure>

II region within this radius would require a gas density $n_H \ge 50$ cm$^{-3}$ (equation 8), which we assume can be retained for $t_{HII} \sim 10^5$ yr.

Using equation (4) with $B_0 = 15$ implies that we can reach $\mu = 330$ at a maximum separation from the caustic of $\max(d_{arcsec}) \approx 2 \times 10^{-3}$ arcsec. If we assume that the magnification is retained at this level or higher from this location and all the way up to the caustic (see Fig. C1), and moreover adopt a length for the cluster caustic of $L_{arcsec} = 100''$, we get a source-plane area (equation 5) $\theta^2 \approx 0.2 \text{ arcsec}^2$ over which the nebula may be rendered detectable in a single cluster-lens field. For $\Delta(z) = 1$ at $z = 6$, this corresponds to a volume of $V_{z,\Delta(z)} \approx 0.13 \text{ cMpc}^{-3}$ (equation 6). Adopting the Liu & Bromm (2020) Pop III SFRD of $\approx 4 \times 10^{-5} \text{ M}_\odot \text{ yr}^{-1} \text{ cMpc}^{-3}$ at $z = 6$ and using equation (7) (with $f(M) = 1$ and $\Delta t_j = t_{HII} = 10^5$ yr), we find that only $N_{\text{detect}} \approx 10^{-3}$ Pop III nebulae are expected to be detected per cluster field.

Hence, unless the H II-region can remain sufficiently compact to attain magnifications in the hundreds for significantly longer than $\sim 10^5$ yr, it is unlikely to be detectable. This result is strictly valid only for the volume spanning $\Delta(z) = 1$ at $z = 6$, but higher redshifts contribute relatively little to the $N_{\text{detect}}$ estimate since $H \alpha$ gets redshifted out of the NIRCam window at $z \ge 6.9$, thereby rendering the Pop III nebula intrinsically much fainter (see Fig. 10).

This paper has been typeset from a TeX/LATEX file prepared by the author.