Magnetic fields limit the mass of Population III stars even before the onset of protostellar radiation feedback Open Access

ABSTRACT

The masses of Population III stars are largely unconstrained since no simulations exist that take all relevant primordial star formation physics into account. We perform the first suite of radiation magnetohydrodynamics (RMHD) simulations of Population III star formation, with the POPSICLE project. Compared to control simulations that only include magnetic fields (MHD), protostellar ionizing and dissociating feedback, or neither, the RMHD simulation best resembles the MHD simulation during the earliest stages of collapse and star formation. In 5000 yr, the mass of the most massive star is |$65\, {\rm M}_{\odot }$| in the RMHD simulation, compared to |$120\, {\rm M}_{\odot }$| in simulations without magnetic fields. This difference arises because magnetic fields act against gravity, suppress mass transport, and reduce compressional heating. The maximum stellar mass of Population III stars is thus already limited by magnetic fields, even before accretion rates drop to allow significant protostellar radiative feedback. Following classical main sequence stellar evolution with MESA reveals that it is difficult to create Population III stars with masses larger than |$600\, {\rm M}_{\odot }$| in typical dark matter mini-haloes at |$z \gtrsim 20$|⁠, with maximum stellar masses |$\sim 100\, {\rm M}_{\odot }$| more likely due to expected negative feedback from both magnetic fields and stellar radiation. This work lays the first step in building a full physics-informed mass function of Population III stars.

1 INTRODUCTION

Understanding the initial mass function (IMF) of Population III (Pop III) stars is of paramount importance, as evidenced by numerous theoretical works that examine the formation of these stars (Bromm 2013; Klessen & Glover 2023, and references therein), as well as indirect observational evidence from metal-poor stars (e.g. Frebel & Norris 2015; Nordlander et al. 2019; Skúladóttir et al. 2021) and |$z>10$| galaxies (e.g. Harikane et al. 2023; Yajima et al. 2023; Maiolino et al. 2024). The IMF is also crucial for determining whether Pop III stars can be observed with JWST, or if their luminosity function can be differentiated from Pop II stars in integrated light measurements (Bromm, Kudritzki & Loeb 2001; Schaerer 2002; Zackrisson et al. 2011; Trussler et al. 2023; Zackrisson et al. 2024; Fujimoto et al. 2025). Equally important is to understand what sets the maximum possible mass (or, the upper mass cutoff of the IMF) of Pop III stars (e.g. Chantavat, Chongchitnan & Silk 2023; Bovill et al. 2024; Liu et al. 2024), which is an essential input for black hole seeding and supermassive star formation (e.g. Haemmerlé et al. 2018; Latif et al. 2022). The upper mass cut-off is also crucial to assess whether Pop III stars of masses between 140 and |$270\, {\rm M}_{\odot }$| existed, and could have exploded as pair-instability supernovae (PISNe; Heger & Woosley 2010; De Bennassuti et al. 2017; Koutsouridou, Salvadori & Skúladóttir 2024). Simulating the collapse of gas and formation of metal-free stars in dark matter mini-haloes provides a robust way to constrain the Population III IMF and its upper mass cut-off.

However, radiation-magnetohydrodynamics (RMHD) simulations in the era of Pop III star formation remain largely absent. Both magnetic fields and protostellar radiation feedback are critical ingredients that influence (massive) star formation across all metallicities (e.g. Tanaka et al. 2018; Chon, Omukai & Schneider 2021; Sharda & Krumholz 2022; Chon et al. 2024). Therefore, conclusions from prior numerical work aimed at deriving the masses of Pop III stars are likely subject to major uncertainties since they either exclude magnetic fields (e.g. Hosokawa et al. 2011, 2016; Sugimura et al. 2020; Jaura et al. 2022) or radiation feedback (e.g. Turk et al. 2012; Sharda, Federrath & Krumholz 2020; Prole et al. 2022; Saad, Bromm & El Eid 2022; Sadanari et al. 2023). In this work, we use the POPSICLE project (Sharda et al., in preparation) to extend the recent RMHD simulations of Pop III star formation by Sharda & Menon (2024) to include far-UV (FUV) molecule-dissociating radiation feedback in addition to extreme-UV (EUV) ionizing feedback. Our aim in this Letter is to show that magnetic fields significantly limit the mass growth of massive Pop III stars, even before radiative feedback becomes dominant. We arrange the remainder of the paper as follows: Section 2 summarizes the setup we use to develop and run the simulations, and Section 3 describes the results. In Section 4, we look at the main-sequence evolution of simulated stars using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019). Finally, we summarize in Section 5.

2 POPSICLE SIMULATIONS

We briefly describe the POPSICLE¹ project setup employed in this work, and point the reader to Sharda & Menon (2024) for details of the numerical implementation. We use a custom version of the adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000; Dubey et al. 2008) which employs an approximate Riemann solver for MHD (Waagan 2009; Waagan, Federrath & Klingenberg 2011). We include non-equilibrium primordial chemistry (with deuterium) from the KROME astrochemistry package (Grassi et al. 2014). We use the VETTAM radiation hydrodynamics scheme, which uses a non-local variable Eddington tensor (VET) closure obtained with a ray-trace solve to close the radiation moment equations (Menon et al. 2022). We use the GENEVA Pop III protostellar model grid to evolve radiative properties of the protostars as a function of their mass and accretion rates (Haemmerlé et al. 2016, 2018). The protostars are represented by sink particles in the simulation, following the criteria described in Federrath et al. (2010b, 2011).

Population III protostars contract and produce significant radiation once accretion rates drop below |$0.01\, {\rm M}_{\odot }\, \mathrm{yr}^{-1}$| (Omukai & Palla 2001, 2003). We consider the ionization of H and H|$_2$| due to extreme-UV (EUV) photons released from the protostar(s) with energies upwards of |$13.6\, \mathrm{eV}$| (upwards of |$15.2\, \mathrm{eV}$| for H|$_2$|⁠). The key physics we add to the simulations by Sharda & Menon (2024) is the dissociation of H|$_2$| by FUV photons in the Lyman–Werner (LW) band, between 11.2 and 13.6 eV, which can significantly affect the thermodynamic state of dense gas in the vicinity of the stars. Importantly, we consider both self-shielding of H|$_2$| as well as cross shielding by H (Wolcott-Green & Haiman 2011); the latter has been ignored in previous works. We adopt the fitting functions for self- and cross-shielding of H|$_2$| from Wolcott-Green, Haiman & Bryan (2011). We do not invoke an external LW background in addition to in-situ LW radiation from the protostars.²

We keep the initial conditions identical to Sharda & Menon (2024): the initial cloud mass and radius are |$1000\, {\rm M}_{\odot }$| and |$1\, \mathrm{pc}$|⁠, respectively. We also include initially trans-sonic turbulence within the cloud, and impose solid body rotation with rotational energy equal to 3 per cent of the gravitational energy. We refine using 64 cells per Jeans length, significantly higher than published simulations that include radiation feedback even at Solar metallicity. With 10 levels of grid refinement based on the Jeans length, the maximum spatial resolution we achieve is |$\Delta x = 7.5\, \mathrm{au}$|⁠. Our sink particle density threshold is |$10^{13}\, \mathrm{cm^{-3}}$|⁠. Based on earlier results that show how even initially weak magnetic fields are quickly amplified to saturation due to the small-scale dynamo (e.g. Schober et al. 2012; Turk et al. 2012; Schober et al. 2015; Sharda et al. 2020; Sharda et al. 2021), we set our initial magnetic field strength to be |$28\, \mu \mathrm{G}$|⁠, equivalent to 10 per cent of the initial turbulent kinetic energy, appropriate for parsec scales in the presence of trans-sonic turbulence and low magnetic Prandtl numbers (for details, see Hirano & Bromm 2018; Sharda et al. 2020, their fig. 4). The power spectrum of the trans-sonic turbulence we initially drive follows |$P_{\mathrm{v}} \propto k^{-1.8}$|⁠, in between the Kolmogorov (⁠|$k^{-5/3}$|⁠; Kolmogorov 1941) and Burgers (⁠|$k^{-2}$|⁠; Burgers 1948) scalings. The turbulence consists of a mixture of solenoidal and compressive modes (Federrath et al. 2010a). The initial magnetic field is completely random with no preferred orientation, as expected for a small-scale turbulent dynamo.

We compare to control simulations without magnetic fields or radiation hydrodynamics (HD), and with only magnetic fields (MHD) from Sharda & Menon (2024). The turbulent realization we select from Sharda & Menon (2024) to control for stochasticity is the one that produces only one star in the HD case, since this case is straightforward to analyse. In addition to this, we run a control simulation with only radiation feedback for the same realization (including both EUV and FUV feedback, termed RHD). These control simulations are instructive as they can help distinguish between the impact of magnetic fields and radiation feedback on Pop III star formation.

3 RESULTS

We evolve the simulations until 5000 yr post the formation of the first star. We find that both runs including magnetic fields show fragmentation, leading to the formation of Pop III star clusters. This means that the evolution of the most massive Pop III star in the MHD and RMHD runs is influenced by companion stars. However, with only one turbulent realization, we lack the statistics to quantify the impact of magnetic fields on fragmentation in primordial clouds. In fact, other turbulent realizations presented in Sharda & Menon (2024) fragment even in the HD case (see also, discussions in Wollenberg et al. 2020; Sharda et al. 2020 on stochastic fragmentation due to turbulence).

Fig. 1 plots the density-weighted projections of the number density of the gas at the end of the simulations. The projection window is |$0.01\, \mathrm{pc}$| wide. It is centred on the isolated star in the HD and RHD runs, and uses the centre of mass of all stars in the MHD and RMHD runs. We plot the density-weighted projections of the gas temperature for the four simulations in Fig. 2, and phase diagrams for the entire cloud in Fig. 3. We see from Figs 1 and 2 that the central star in the HD and RHD runs is fed by a well-defined, thermal pressure-supported accretion disc that remains stable against fragmentation for the entirety of the simulation. The high temperatures in these runs is a result of dissociation of H|$_2$| by shocks, a process that is only resolved in simulations with sufficiently high Jeans resolution (e.g. Turk et al. 2012; Sharda et al. 2021; Sharda & Menon 2024). The key difference between the HD and RHD runs is the presence of cooler, H|$_2$|-dominated gas near the star in the latter, where H|$_2$| dissociation is prevented due to radiation pressure slowing down accretion shocks (see also, Sharda & Menon 2024). Given the relatively high accretion rates, the star is unable to contract and produce significant ionizing (EUV) photons that can ionize H or H|$_2$| (e.g. Omukai & Palla 2001; Hosokawa & Omukai 2009; Haemmerlé et al. 2018). In contrast to HD and RHD simulations, the gas is cooler in MHD and RMHD runs across the entire central envelope. Magnetic fields suppress gravitational collapse, thereby also reducing the rate of compressional heating. Reduced heating allows the gas to remain molecular, and H|$_2$| cooling further ensures gas temperatures remain low.

These figures show, for the first time, the simultaneous effects of magnetic fields, and ionizing as well as dissociating radiation feedback during the evolution of primordial clouds. Note that we do not include the effects of radiative heating caused by accretion luminosity, which could lead to higher gas temperatures close to the protostars, thereby reducing the accretion rates and limiting protostellar mass growth (e.g. Smith et al. 2011; Wollenberg et al. 2020). However, it cannot hinder accretion for long time periods since low accretion rates in turn lead to lower accretion heating. Further, accretion luminosity heating only becomes significant when the opacity of (primordial) gas is large, which occurs for |$T > 5000\, \mathrm{K}$| at the densities we resolve (Mayer & Duschl 2005, table E3). Given that magnetic fields lead to lower average gas temperatures, we expect accretion luminosity to be less important in the MHD and RMHD simulations as compared to HD and RHD simulations.

We now turn to look at the time evolution of gas properties that dictate the mass growth of the protostars. The two panels of Fig. 4 plot the evolution of the mass enclosed within an envelope surrounding the accretion disc–star system, and the accretion disc itself in the four simulations. Following Sharda et al. (2021) and Sharda & Menon (2024), we define the envelope to be a |$0.01\, \mathrm{pc}$| spherical region centred on the most massive star. Similarly, we define the accretion disc as a cylindrical region of radius 500 and |$50\, \mathrm{au}$| in height (from the midplane). The dashed curves in the MHD and RMHD runs mark the onset of fragmentation, emphasizing that the most massive star does not evolve in isolation thereafter.

We see from Fig. 4 that the amount of mass reaching the envelope initially increases in the MHD and RMHD runs, but then turns over and starts to decline. Despite the envelope and the disc containing larger masses early on in the MHD and RMHD runs, the corresponding accretion rates onto the protostars are lower (see Fig. 5). Close to the protostar, the magnetic field is sub-Alfvénic (plasma |$\beta < 1$|⁠), significantly inhibiting mass transport. In contrast, the rate of mass transfer from the envelope to the accretion disc is initially quite fast in the HD and RHD runs, resulting in a decline of mass in the envelope and buildup of mass in the disc. This mass is consequently accreted by the star at a high rate. Fragmentation events correlate with a large buildup of mass in the accretion disc such that the ratio of mass in the disc to the protostellar mass increases far beyond unity. The presence of multiple stars reduces the mass available within the disc of the primary star in the MHD and RMHD runs, whereas the disc mass in the HD and RHD runs continues to build up.

The top panel of Fig. 5 shows that the accretion rates (of the most massive star) are of the order of |$0.03\, {\rm M}_{\odot }\, \mathrm{yr}^{-1}$| in the first |$1000\, \mathrm{yr}$|⁠, and decline as time progresses. The accretion rates in the MHD and RMHD runs remain systematically lower than the HD and RHD runs for the most part, initially due to strong magnetic fields inhibiting accretion, and later due to fragmentation. The key impact of lower accretion rates in the MHD and RMHD runs is that the most massive stars only build up half as much mass as the HD and RHD runs within the same time period. We show this in the bottom panel of Fig. 5. The mass of the (isolated) star at the end of the simulations in the HD and RHD runs is |$127\, {\rm M}_{\odot }$| and |$120\, {\rm M}_{\odot }$|⁠, respectively. On the other hand, the mass of the most massive star in the MHD and RMHD runs is |$48\, {\rm M}_{\odot }$| and |$67\, {\rm M}_{\odot }$|⁠, respectively. However, the integrated star formation efficiency (SFE; defined as the cumulative mass of all stars normalized by the initial cloud mass; solid orange and purple curves in the bottom panel of Fig. 5) is lower by 30 per cent in the RMHD run as compared to all the control runs.

There is a subtle but important difference between magnetic fields suppressing mass transport, and fragmentation-induced starvation. The former turns on earlier, and affects both the mass of the primary star and the integrated SFE (dashed and solid lines in the bottom panel of Fig. 5), whereas the latter emerges later, and affects the mass of the primary star mass without necessarily affecting the integrated SFE. Thus, the ‘mass-limiting’ effects of the runs including magnetic fields are not simply due to fragmentation-induced starvation. The slower mass growth and consequently lower star formation efficiency in the presence of magnetic fields has also been observed in simulations of Population I massive star formation, although the differences arise much later on, and are smaller in magnitude (e.g. Rosen & Krumholz 2020; Kim, Ostriker & Filippova 2021). Our findings are also in qualitative agreement with Sharda & Menon (2024), where the authors simulated Pop III star formation with magnetic fields and only ionizing feedback.

4 EVOLUTION ON THE MAIN SEQUENCE

In this section, we explore the implications for the long-term mass growth of Pop III stars and the upper mass cut-off of the Pop III IMF. Fully simulating this process is beyond the scope of this work. Nonetheless, we can use the information available from our simulations to predict the time-scales to reach the zero-age main sequence (ZAMS), as well as the ultimate fate as these stars reach the terminal age main sequence (TAMS). For this purpose, we run stellar structure evolution calculations using MESA, version 23.05.1 (Paxton et al. 2011, 2013, 2015, 2018, 2019). Appendix  A lists the details of the modelling.

We assume no accretion in our MESA models for the first |$\sim 25$| yr, corresponding to the amount of time the most massive star in the RMHD simulation takes to reach the adopted initial stellar mass in MESA (⁠|$1\, \mathrm{M}_\odot$|⁠). After that, we use the RMHD accretion history of said most massive star, and consider three increasingly realistic scenarios for the evolution beyond the final time of our simulations (⁠|$5000\, \mathrm{yr}$|⁠). In Case A (Fig. 6), we assume the star continues to accrete at a constant rate of |$0.005\, {\rm M}_{\odot }\, \mathrm{yr}^{-1}$|⁠, corresponding to the average accretion rate over the final |$1000\, \mathrm{yr}$| of the simulation. Such an accretion history for the entirety of the stellar lifetime is rather unrealistic, but we include it to provide a baseline for comparison. In Case B, we assume the average accretion rate beyond |$5000\, \mathrm{yr}$| follows the trend observed in the simulations, decreasing with time as |$t^{-0.65}$|⁠. In Case C, we empirically include the effects of radiation feedback that can halt accretion at late times. To do so, we define an accretion rate that goes as |$e^{-t^2}$| between |$5000 \le t \le 2\times 10^4\, \mathrm{yr}$| and linearly declines as |$t^{-3}$| beyond |$2\times 10^4\, \mathrm{yr}$|⁠, mimicking the trend seen in radiation hydrodynamics simulations of Hosokawa et al. (2011, fig. 3).

We plot the resulting mass growth from the three cases as blue, orange, and green curves in Fig. 6, respectively. Despite different accretion histories, the star in all the three cases reaches ZAMS around |$13\, 000\, \mathrm{yr}$| (marked by diamonds in Fig. 6) with stellar mass between 90 and |$100\, {\rm M}_{\odot }$|⁠. However, the main-sequence evolution is significantly different, and the TAMS mass spans a large range. The stars are expected to turn into black holes soon after TAMS (denoted by circles in Fig. 6). The star in Case A oscillates between two branches, leading to wiggles in the HR diagram, because its evolution is sensitive to atmospheric parameters (Herrington, Whalen & Woods 2023). This phenomenon only occurs for Case A because the accretion rate is close to the critical value of |$0.01\, {\rm M}_{\odot }\, \mathrm{yr}^{-1}$| which separates the blue and red branches of supermassive stellar evolution (Hosokawa & Omukai 2009; Haemmerlé et al. 2018). The period we simulate with different physics, although short, matters for the final fate of the star. To show this, we use Case C to extrapolate the growth of the star in the RHD simulation (dashed green curve in the left panel of Fig. 6). Due to differences in accretion histories in the first |$5000\, \mathrm{yr}$|⁠, this star ends into the regime of PISNe, rather than ending its life as a black hole. This is unlikely to occur in reality since magnetic fields are excluded in the RHD simulation, but it shows that omitting key physical processes even during the early stages of evolution can profoundly impact the ultimate fate of Pop III stars.

Given the high stellar effective temperatures and luminosities on the main sequence (right panel of Fig. 6), radiation feedback combined with competitive accretion due to fragmentation will potentially limit the final mass to |$80\, {\rm M}_{\odot } \le M_{\star } \lesssim 600\, {\rm M}_{\odot }$|⁠, making Case A highly unlikely. The magnetic field strength is expected to remain significant beyond the period we simulate here due to a mixture of turbulent and mean field dynamos (Liao, Turk & Schive 2021; Sharda et al. 2021), and it can provide additional negative feedback by suppressing accretion onto the star, or generating outflows (Machida & Doi 2013). Together with the effects of heating due to accretion luminosity we discuss in Section 3, this will enhance feedback effects and further limit the maximum stellar mass, such that the actual stellar mass is closer to |$80\, {\rm M}_{\odot }$| than |$600\, {\rm M}_{\odot }$|⁠.

5 SUMMARY

We use the POPSICLE simulations suite to perform RMHD simulations of Population III star formation. We simultaneously include non-equilibrium primordial chemistry, turbulent magnetic fields, ionizing and dissociating stellar feedback, following the evolution |$5000\, \mathrm{yr}$| post the formation of the first star. We also carry out control simulations where we use identical initial conditions but only include magnetic fields (MHD), protostellar radiation feedback (RHD), or exclude both (HD).

We find that, during the earliest stages we simulate, magnetic fields suppress gravitational collapse, leading to less compressional heating and inefficient mass transport from the cloud to the protostar. As a result, the gas temperature in the MHD and RMHD runs are lower than the HD and RHD runs, which makes the gas more susceptible to fragmentation. Both the runs including magnetic fields (MHD and RMHD) fragment, however, with only one turbulent realization, we lack the statistics to make quantitative conclusions on how fragmentation occurs in the presence of both magnetic fields and radiation feedback (e.g. Sharda et al. 2020; Wollenberg et al. 2020). The combined effect of suppression of mass transport to the star and fragmentation-induced starvation is that the mass of the most massive star in the RMHD run is |$65\, {\rm M}_{\odot }$|⁠, factor of two lower than that in the HD and RHD runs.

Radiation feedback has long been proposed as the primary mechanism that halts the growth of Pop III stars and sets the upper mass cutoff of the Pop III IMF (Hosokawa et al. 2011; Stacy, Greif & Bromm 2012; Hosokawa et al. 2016; Stacy, Bromm & Lee 2016). Here, we show that magnetic fields can also limit mass growth, even before accretion rates drop to facilitate significant ionizing or dissociating feedback (Haemmerlé et al. 2018). Because magnetic fields also slow down accretion onto protostars, they can induce strong radiation feedback earlier than expected. If a mean field dynamo operates and is sustained for long periods of time, magnetic fields will also enable launching protostellar outflows, further reducing the maximum possible mass of Pop III stars (Machida & Doi 2013; Liao et al. 2021; Sharda et al. 2021).

Stellar structure modelling with MESA using the accretion history from the RMHD simulation, extrapolated considering a range of scenarios for subsequent mass growth, shows that the |$65\, {\rm M}_{\odot }$| star reaches ZAMS early on at an age of |$13\, 000\, \mathrm{\, yr}$|⁠, and likely continues to accrete on the main sequence. The star then evolves off the main sequence in about |$2.5\, \mathrm{Myr}$|⁠. The initial evolutionary phase has implications both for the final mass and fate (supernova versus black hole) of the star. The final mass spans a large range based on possible accretion histories, from |$80\, {\rm M}_{\odot }$| to |$600\, {\rm M}_{\odot }$|⁠, although combined effects of radiation feedback, magnetic fields and fragmentation will render the actual mass closer to the lower bound. This work lays the foundation for constructing a Pop III IMF with all relevant star formation physics. Future work will involve simulating multiple realizations to build a mass distribution, following the accretion history for longer time periods to quantify feedback at late times, and exploring black hole seeding from Pop III stars.

ACKNOWLEDGEMENTS

We thank Zoltan Haiman, Simon Glover, Mark Krumholz, Mariska Kriek, and Matthieu Schaller for useful discussions, and the anonymous referee for their feedback that helped improve this work. We also thank SURF for the support in using the Dutch National Supercomputer Snellius. PS is supported by the Leiden University Oort Fellowship and the International Astronomical Union–Gruber Foundation (TGF) Fellowship. SHM and BB acknowledge support through NASA grant No. 80NSSC20K0500 and National Science Foundation (NSF) grant AST-2009679, and the CCA at the Flatiron Institute. RG is supported by the Society of Science Fellowship provided by the University of Notre Dame. BB is grateful for generous support by the David and Lucile Packard Foundation and Alfred P. Sloan Foundation. BDW acknowledges support from NSF grant AAG-2106575. The simulations and data analyses presented in this work used high-performance computing resources provided by the Dutch National Supercomputing Facility SURF via project grant EINF-8292 on Snellius, as well as high-performance computing resources provided by the Centre for Computational Astrophysics (CCA) at the Flatiron Institute, and the Australian National Computational Infrastructure (NCI) through project jh2 in the framework of the National Computational Merit Allocation Scheme and the Australian National University (ANU) Allocation Scheme. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. We acknowledge using the following softwares: flash (Fryxell et al. 2000; Dubey et al. 2013), VETTAM (Menon et al. 2022), KROME (Grassi et al. 2014), petsc (Balay et al. 1997; Zhang et al. 2022), Astropy (Astropy Collaboration 2013, 2018, 2022), numpy (Oliphant 2006; Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), yt (Turk et al. 2011), cmasher (Van der Velden 2020). This research has made extensive use of NASA’s Astrophysics Data System Bibliographic Services (ADS).

DATA AVAILABILITY

Movies associated with this article are available as supplementary (online only) material. To access the suite of simulations or properties of the stars, please contact the authors.

Footnotes

REFERENCES

APPENDIX A: DETAILS OF MESA MODELLING

In general, MESA simulations cannot account for the early phases of stellar evolution, as the simulated object is required to be in approximate hydrostatic equilibrium (up to a set of implemented dynamical corrections) and obey the equations of stellar structure. For this reason, it is necessary to choose a non-zero initial mass and the initial structure of the star at the zero age of the simulation. We initialize all MESA simulations with a metal-free composition (⁠|$X=0.75$|⁠, |$Y=0.25$|⁠, |$Z=0$|⁠), initial mass |$1\ \mathrm{M}_\odot$|⁠, and initial central temperature |$T_c=61\,500\ \mathrm{K}$|⁠. For the accretion rate predicted by the RMHD simulation, the chosen values approximately correspond to the lowest initial mass and |$T_c$|⁠, for which the model can converge. In order to capture as much of the early stellar evolution as possible, it is desired to choose the smallest possible values of initial mass and |$T_c$|⁠, as they likely represent the earliest phases of the pre-main sequence contraction. MESA uses an adaptive time-step that depends on the rate of evolution; however, we found the default implementation of the adaptive algorithm too generous to capture the subtleties of the pre-main sequence evolution with rapid accretion. In all three cases, we average the accretion history provided to MESA in |$10\ \mathrm{yr}$| steps to suppress discontinuities in the evolutionary tracks.

We define the ZAMS point as the earliest age, at which 0.1 per cent of the hydrogen content at the centre of the star has been used up by nuclear fusion. The TAMS is typically defined by the complete exhaustion of hydrogen in the core, which is subsequently followed by the onset of the red giant branch, characterized by ignition of nuclear fusion away from the centre of the star (Gamow & Keller 1945). However, in a massive Pop III star, the latter process may begin before the core hydrogen reservoir is fully exhausted, making the TAMS point less distinct. The carbon-nitrogen-oxygen (CNO) cycle serves as the primary energy source in the cores of massive Population III stars (Larkin, Gerasimov & Burgasser 2023) even for initially metal-free compositions. On the other hand, a rapid increase in the energy production rate by the proton–proton chain may be used as an early signature of nuclear burning outside the core, where the conditions are less favourable for the CNO cycle. We therefore define TAMS as the earliest age, at which (1) the energy production rate due to the proton–proton chain is increasing, and (2) the central hydrogen fraction has dropped below 5 per cent. The first condition is necessary as the energy production of the proton–proton chain is also expected to increase during the pre-main sequence evolution.