© The Authors 2026

1. Introduction

Understanding the nature of dark matter (DM) remains one of the most pressing challenges in astrophysics and cosmology. Observational evidence from galactic rotation curves, gravitational lensing, and cosmic microwave background anisotropies strongly suggests the existence of DM, which is believed to constitute about 85% of the matter content of the Universe (Planck Collaboration VI 2020). Despite its elusive nature, DM plays a crucial role in the formation and evolution of cosmic structures.

Among the various DM candidates, weakly interacting massive particles (WIMPs) are particularly compelling due to their natural emergence in extensions of the standard model of particle physics (Jungman et al. 1996; Bertone et al. 2005). WIMPs interact through gravity and potentially via weak-scale interactions, allowing them to scatter off baryonic matter and annihilate with each other. In regions of high DM density, such as the early Universe or the centres of galaxies and stars, WIMP annihilation can release significant amounts of energy (Bergström et al. 1998).

The first generation of stars, known as Population III (Pop III) stars, formed in DM halos at high redshifts (z ∼ 20–30; Bromm & Larson 2004; Maeder & Meynet 2012; Liu et al. 2025). These stars are thought to be massive and metal-free, playing a pivotal role in the re-ionisation of the Universe and the enrichment of the interstellar medium with heavy elements (Liu et al. 2021). The formation and evolution of Pop III stars are influenced by the conditions of the early Universe, including the presence of DM. In particular, the concept of ‘dark stars’ (DSs) has been proposed, where DM annihilation provides an additional energy source that can significantly alter stellar structure and evolution (Spolyar et al. 2008; Freese et al. 2010; Ilie et al. 2023b). These hypothetical stars form when the gravitational collapse of baryonic matter in early DM halos leads to the accumulation of both ordinary matter and DM particles. The annihilation of WIMPs within the star’s core can counteract gravitational contraction, resulting in stars that are cooler, larger, and have extended lifetimes compared to conventional Pop III stars (Freese et al. 2008b; Ilie et al. 2012; Rindler-Daller et al. 2015; Tan et al. 2024; Nandal et al. 2025b). The prolonged lifetimes and unique characteristics of DSs may leave observable imprints, such as distinct spectral features or anomalous elemental abundances, which could be probed by current-generation telescopes like the James Webb Space Telescope (JWST; Zackrisson et al. 2010; Freese et al. 2016; Banik et al. 2019; Ilie et al. 2025).

The interaction between DM and stars has been extensively studied over the past few decades. Gould (1987) developed the theoretical framework for WIMP capture and annihilation within stars, laying the foundation for understanding how DM can influence stellar evolution. This framework was further refined to account for various stellar environments and DM properties (Press & Spergel 1985; Faulkner & Gilliland 1985). Subsequent works have explored the impact of DM on stellar evolution in greater detail, including the potential formation of DSs (Freese et al. 2008a; Ilie et al. 2012; Scott et al. 2009; Ilie et al. 2021). Studies have shown that DM annihilation can alter the temperature gradients within stars, affect nucleosynthesis processes, and modify the stellar evolution tracks on the Hertzsprung-Russell diagram (HRD; Taoso et al. 2008; Casanellas & Lopes 2011b,a). The effects of DM annihilation on Pop III stars, according to both rotating and non-rotating models, have significant implications for our understanding of the first stars and the role of DM in stellar physics (Yoon et al. 2008; Hirano et al. 2015; Stacy et al. 2012).

Tan et al. (2024) sub-classified Pop III stars into Pop III.1 and Pop III.2 types based on their formation environments and the influence of radiation feedback. Pop III.1 stars form in pristine environments without prior star formation, resulting in isolated, massive stars of up to several hundred solar masses (Hirano et al. 2014). In contrast, Pop III.2 stars form later, influenced by radiation and chemical feedback from earlier generations, leading to a broader mass spectrum and different formation pathways (Greif et al. 2012). The presence of DM can further complicate this classification, as DM annihilation may play a more pronounced role in the formation and evolution of Pop III.1 stars due to the higher DM densities in early halos (Natarajan et al. 2009). Recent simulations have investigated the interplay between DM annihilation and stellar rotation, suggesting that rotation can enhance the mixing processes within stars and affect the transport of angular momentum and chemical elements (Haemmerlé et al. 2020; Stacy et al. 2013). Understanding these effects is crucial for predicting the observational signatures of the first stars and their contribution to cosmic re-ionisation and early chemical enrichment (Bromm et al. 2009; Wise et al. 2012; Gondolo et al. 2022).

In this work we investigated the effects of DM annihilation on the evolution of rotating and non-rotating massive Pop III stars. Using updated stellar evolution models, we analysed how DM capture and annihilation influence stellar structure, lifetimes, and nucleosynthesis. By comparing models with and without rotation, we aim to elucidate the interplay between DM physics and stellar evolution in the context of primordial stars.

The paper is organised as follows: In Sect. 2 we describe the methods and numerical setup used in our simulations. In Sect. 3 we present the results of our models and discuss the implications for Pop III star evolution. Finally, in Sect. 5 we summarise our findings and outline prospects for future work.

2. Methods

In this work we modelled the impact of DM annihilation on Pop III stellar evolution by incorporating the capture and annihilation of WIMPs into the Geneva stellar evolution code (GENEC). Our formalism builds on that of Taoso et al. (2008), hereafter T08) while reflecting improvements in the underlying physics and numerical techniques (Ekström et al. 2012; Nandal et al. 2024, 2025a). In the following, we detail the methodology and present the relevant equations, clarifying all assumptions and definitions.

The evolution of the total number of WIMPs in the star, N_χ, is governed by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\chi }}{\mathrm{d}t}=C_c+C_s{N}_{\chi }-A{N}_{\chi }^2-E{N}_{\chi },\end{array}$$(1)

where C_c is the capture rate via WIMP–baryon scattering, C_s is the self-capture rate due to WIMP–WIMP interactions, A is the annihilation coefficient, representing the rate at which pairs of WIMPs annihilate, and E denotes the evaporation rate. For the massive stars studied here, both self-capture and evaporation are negligible; hence, these terms are set to zero in our models. The total number of WIMPs is defined as

$$N_{\chi }={\displaystyle {\int }_0^{R_{\ast }}}{n}_{\chi }(r)\mathrm{d}V,$$

with n_χ (r) being the local WIMP number density and R_∗ the stellar radius.

We adopted a Gaussian profile to describe the spatial distribution of WIMPs within the star:

$$\begin{array}{c}\hfill n_{\chi }(r)=n_0{e}^{-r^2/r_{\chi }^2},\end{array}$$(2)

where the characteristic radius,

$$\begin{array}{c}\hfill r_{\chi }=\sqrt{\frac{3k{T}_c}{2\pi G{\rho }_c{m}_{\chi }}},\end{array}$$(3)

marks the region where WIMP–baryon interactions are most effective. Here T_c and ρ_c are the central temperature and density, m_χ is the WIMP mass, k is Boltzmann’s constant, and G is the gravitational constant.

2.1. Capture rate

Weakly interacting massive particles entering the star lose kinetic energy via scattering off baryons until they become gravitationally bound. The local differential capture rate is expressed as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{C}_c(r)}{\mathrm{d}V}& ={\left(\frac{6}{\pi }\right)}^{1/2}{\sigma }_{\chi ,N}\frac{{\rho }_i(r)}{M_i}\frac{{\rho }_{\chi }}{m_{\chi }}\frac{v^2(r)}{{\overline{v}}^2}\frac{\overline{v}}{2\eta A^2}\hfill \\ \hfill & \times \{(A_{+}{A}_{-}-\frac{1}{2})[\chi (-\eta ,\eta )-\chi (A_{-},A_{+})]\hfill \\ \hfill & +\frac{1}{2}{A}_{+}{e}^{-A_{-}^2}-\frac{1}{2}{A}_{-}{e}^{-A_{+}^2}-\frac{1}{2}\eta e^{-{\eta }^2}\}.\hfill \end{array}$$(4)

The variables are defined as follows:

$$\begin{array}{cccc}\hfill A^2& =\frac{3v^2(r)\mu }{2{\overline{v}}^2{\mu }_{-}^2},\hfill & \hfill A_{\pm }& =A\pm \eta ,\hfill \\ \hfill {\eta }^2& =\frac{3v_{\ast }^2}{2{\overline{v}}^2},\hfill & \hfill \chi (a,b)& =\frac{\sqrt{\pi }}{2}[\text{Erf}(b)-\text{Erf}(a)],\hfill \\ \hfill {\mu }_i& =\frac{m_{\chi }}{M_i},\hfill & \hfill {\mu }_{-}& =\frac{{\mu }_i-1}{2}.\hfill \end{array}$$

In these expressions, ρ_i(r) is the density profile of element i (with atomic mass M_i), σ_χ,N is the WIMP–baryon scattering cross-section (which can be spin-dependent or spin-independent), and ρ_χ and $\overline{v}$ are, respectively, the ambient WIMP density and velocity dispersion. The stellar velocity v_∗ is assumed to be equal to $\overline{v}$. The escape velocity, v(r), is computed from

$$\begin{array}{c}\hfill v^2(r)={\displaystyle {\int }_{\infty }^r}\frac{GM(r\prime )}{r^{\prime 2}}\mathrm{d}r\prime ,\end{array}$$(5)

where M(r′) is the mass enclosed within radius r′. The total capture rate is finally obtained by integrating Eq. (4) over the entire stellar volume:

$$\begin{array}{c}\hfill C_c=4\pi {\displaystyle {\int }_0^{R_{\ast }}}{r}^2\frac{\mathrm{d}{C}_c(r)}{\mathrm{d}V}\mathrm{d}r.\end{array}$$(6)

2.2. Annihilation rate

Assuming that WIMPs are their own antiparticles (Jungman et al. 1996; Aprile et al. 2017), they can annihilate in the dense stellar core. The local energy generation rate per unit mass due to WIMP annihilation is given by

$$\begin{array}{c}\hfill {ϵ}_a(r)=\frac{1}{2}⟨\sigma v⟩m_{\chi }{c}^2\frac{n_{\chi }^2(r)}{\rho (r)},\end{array}$$(7)

where ⟨σv⟩ is the thermally averaged annihilation cross-section, c is the speed of light, and ρ(r) is the local baryon density. The integrated annihilation coefficient is then defined as

$$\begin{array}{c}\hfill A=4\pi {\displaystyle {\int }_0^{R_{\ast }}}{ϵ}_a(r)\rho (r)r^2\mathrm{d}r.\end{array}$$(8)

We assumed that the totality of the annihilation energy is processed by the star. Hence,

$$\begin{array}{c}\hfill L_{\chi }=m_{\chi }{C}_c\end{array}$$(9)

is the net luminosity deposited in the star. Freese et al. (2010), Rindler-Daller et al. (2015), and Ilie et al. (2021, 2025) consider that a third of the energy is lost by neutrinos escaping the star.

2.3. Normalisation and equilibrium

A key unknown in our formulation is the normalisation constant (n₀) of the WIMP density profile. To determine n₀, we assumed that after a characteristic timescale,

$$t_{\chi }=\sqrt{\frac{1}{C_c⟨\sigma v⟩}{\pi }^{-3/2}{r}_{\chi }^{-3}},$$

the system attains equilibrium such that the rate of WIMP annihilation balances out the capture rate. Under the simplifying assumption that self-capture and evaporation are negligible, the equilibrium condition can be approximated by

$$2A{N}_{\chi }^2=C_c.$$

Once C_c is computed via Eq. (6), n₀ is adjusted in GENEC so that the energy injection from WIMP annihilation, as given by Eq. (7), satisfies this equilibrium condition. Although this approach involves approximations regarding the equilibrium timescale t_χ, the resulting uncertainties remain within acceptable limits for our analysis. This methodology self-consistently integrates DM capture and annihilation into the stellar energy budget, allowing us to investigate how WIMP annihilation modifies stellar structure and evolution.

3. Results

We summarise the suite of 20 M_⊙ Pop III star models in Table 1. Columns 1–4 list the model name and key initial parameters: ambient WIMP density ρ_χ, initial rotation (as a percentage of the critical velocity v_crit), and the spin-dependent WIMP scattering cross-section σ_SD. Columns 5–7 give selected results at the end of core hydrogen burning: the helium core mass M_He, the equatorial surface velocity v_eq, and the stellar lifetime. Non-rotating models (v/v_crit = 0%) are listed in the upper part of the table, followed by rotating cases (20%, 40%, and 60% of v_crit). In cases where the main-sequence (MS) evolution did not fully converge (defined here by central hydrogen mass fraction X_c reaching 10⁻³), the table entries are marked with an ‘X’. As seen in Table 1, increasing the ambient WIMP density (for fixed σ_SD) dramatically extends the MS lifetime of a 20 M⊙ star. For example, in the non-rotating sequence at σ_SD = 10⁻³⁸ cm², the lifetime grows from ∼9.6 Myr with no DM up to ∼8.44 × 10³ Myr at the highest DM density (ρ_χ = 3 × 10¹⁰GeVcm⁻³). The helium core mass M_He remains around 5 M_⊙ for those non-rotating models that complete H-burning, indicating that despite the vastly different MS durations, they synthesise similar helium cores. In contrast, the rotating models (which is discussed in Sect. 3.2) develop larger M_He and higher v_eq by MS termination due to rotational mixing. A second set of non-rotating models (Table 1, last rows) explores a lower cross-section σ_SD = 10⁻⁴⁰ cm² with proportionally increased ρ_χ: these yield evolution outcomes very similar to the σ_SD = 10⁻³⁸ cases, reflecting the primary dependence on the product ρ_χσ_SD.

3.1. Non-rotating models

3.1.1. Hertzsprung-Russell evolution

Figure 1 illustrates the evolutionary tracks of the non-rotating models in the HRD. Increasing the ambient WIMP density shifts the zero age main sequence (ZAMS) to markedly lower effective temperatures and slightly lower luminosities. Physically, the additional energy injection from WIMP annihilation in the stellar core inflates the star, increasing internal radiative pressure and yielding a larger radius (cooler T_eff) at the ZAMS. Throughout the early MS, the energy output is dominated by WIMP annihilation rather than nuclear fusion. Consequently, the star undergoes only minimal hydrogen burning during this phase, and its track in the HRD remains near the ZAMS locus for an extended period. In the HRD (Fig. 1), this behaviour is seen as an initial clustering of the tracks at relatively cool T_eff for high ρ_χσ_SD models. As evolution proceeds, however, the DM energy input gradually becomes less dominant. Once a significant amount of hydrogen is converted into helium in the core, the WIMP capture rate drops (since we assume purely spin-dependent capture; helium, with no nuclear spin, is ineffective at capturing WIMPs). Thereafter, nuclear energy production steadily increases and the star resumes a more conventional MS evolution. Accordingly, the HRD tracks for different ρ_χσ_SD eventually bend back towards higher T_eff, converging with the track of a star without DM by the end of the MS. Aside from an initial offset in ZAMS position, the overall morphology of our HRD tracks with DM is similar to that found by Taoso et al. (2008, dashed curves in Fig. 1). Small quantitative differences (e.g. a slight shift in our ZAMS location) can be attributed to updates in the GENEC code since T08. Importantly, once WIMP annihilation becomes subdominant late in the MS, all tracks, with and without DM, converge to the same Hertzsprung–Russell endpoint.

Fig. 1. HRD from this work (solid lines) and T08 (dashed lines) for a static 20 M⊙ Pop III star surrounded by different WIMP densities (see the colour code).

3.1.2. Main-sequence lifetimes

The non-rotating series in Table 1 quantifies the dramatic lifetime extension produced by WIMP heating. The reference model without DM completes core–hydrogen burning in 9.6 Myr, whereas a mild enhancement ρ_χ = 6.3 × 10⁹ GeV cm⁻³ (63d9, ρ_χσ_SD = 6.3 × 10⁻²⁹ GeV cm⁻¹) already doubles the lifetime to 23.5 Myr. Further increases give 50.9 Myr at 1.0 × 10¹⁰ (1d10), 7.2 × 10² Myr at 2.0 × 10¹⁰ (2d10), and 8.4 Gyr at 3.0 × 10¹⁰ (3d10)–a factor of 10³ relative to the standard Pop III case. The alternate sequence with σ_SD = 10⁻⁴⁰ cm² reproduces the same lifetimes when ρ_χ is increased by the same factor, demonstrating that the controlling parameter is the product ρ_χσ_SD.

The reference model without DM completes core–hydrogen burning in 9.6 Myr, whereas a mild enhancement ρ_χσ_SD = 6.3 × 10⁻²⁹ GeV cm⁻¹ (63d29) already doubles the lifetime to 23.5 Myr. Further increases give 50.9 Myr at 1.0 × 10⁻²⁸ (1d28), 7.2 × 10² Myr at 2.0 × 10⁻²⁸ (2d28), and 8.4 Gyr at 3.0 × 10²⁸ (3d28) —a factor of 10³ relative to the standard Pop III case. The alternate sequence with σ_SD = 10⁻⁴⁰ cm² reproduces the same lifetimes when ρ_χ is increased by the same factor, demonstrating that the controlling parameter is the product ρ_χσ_SD.

During the early MS, energy injection from WIMP annihilation maintains a comparatively cool, low-density core; hydrogen fusion therefore proceeds slowly, and the stellar structure evolves on an extended nuclear timescale. As hydrogen is converted to helium, spin-dependent capture weakens, nuclear burning gradually overtakes the DM contribution, and the model ultimately resumes a standard terminal-age MS path. For ρ_χσ_SD ≳ 10⁻²⁸ GeV cm⁻¹ the MS duration exceeds a gigayear, implying that such WIMP-supported Pop III stars could persist to low redshifts if their DM environment remains undisturbed.

3.1.3. Internal structure

Figure 2 presents the response of the stellar core to increasing WIMP capture for seven central hydrogen mass fraction during the MS phase (from X_c = 0.70 to 0.005). The left-hand panel shows the central temperature T_c for a non-rotating star, while the right-hand panel displays a star rotating with initial velocity 20% of v. Both models are plotted against the product ρ_χσ_SD, which parametrises the WIMP capture rate.

Fig. 2. Central temperature as a function of WIMP density for a 20 M⊙ Pop III star without rotation (left) and with a rotation 20% of vcrit (right). Curves corresponding to different evolutionary stages (defined by the central hydrogen abundance) are plotted with spline interpolation; the dashed curves are taken from T08.

At ZAMS (X_c = 0.70) the model without DM occupies (T_c, ρ_c)≃(8.0 × 10⁷ K, 1.2 × 10² g cm⁻³). A tenfold increase in ρ_χσ_SD lowers each quantity by roughly one quarter, and at 3 × 10⁻²⁸ GeV cm⁻¹ both are reduced by nearly a factor of two, indicating that instead of hydrogen fusion, WIMP annihilation provides the dominant luminosity. The mid-MS curves (X_c = 0.40) illustrate how this stabilising influence persists. In the standard model the core has already contracted and heated to (9.5 × 10⁷ K, 1.5 × 10² g cm⁻³), whereas DM–rich cores remain within ∼10% of their ZAMS state, confirming that the additional energy source largely suppresses Kelvin–Helmholtz contraction during the first few 10⁷ yr.

Once hydrogen is exhausted, spin-dependent capture ceases and the tracks for all ρ_χσ_SD converge to a narrow locus around T_c ≃ 1.0–1.1 × 10⁸ K and ρ_c ≃ 1.6–1.8 × 10² g cm⁻³. Thus, despite large differences during most of the extended MS, each model ultimately attains the core conditions required for helium ignition after WIMP heating subsides. In effect, DM energy support delays–but does not prevent–the standard contraction–heating sequence; the star therefore resumes canonical Pop III evolution once its hydrogen reservoir no longer sustains efficient WIMP capture.

3.2. Rotating models

Rotating Pop III models were evolved with GENEC ’s standard meridional-circulation and shear prescription. We first compared stars that begin at 20 % v_crit with the non-rotating sequence, then increased the spin to 40 and 60 % v_crit. Unless stated otherwise, the spin-dependent cross-section was σ_SD = 10⁻³⁸ cm². The companion grid with σ_SD = 10⁻⁴⁰ cm² reproduces the same trends after ρ_χ is raised by the same factor and is noted briefly in Sect. 3.1.2.

3.2.1. Moderate rotation (20 % v_crit) versus no rotation

The right-hand side of Fig. 2 shows the central temperature as a function of the product ρ_χσ_crit for a Pop III rotating star 20% of v_crit. In opposition to the static case the temperature stays lower at every step during the MS. With the mixing the energy produce by WIMPs annihilation is transported outwards, keeping the core cooler.

Figure 3 shows that rotation leaves the ZAMS almost unchanged but shifts the subsequent DS phase upwards in luminosity. wtt_0.2v is ≃0.05 dex brighter and 5.0 × 10² K hotter, from the static model, at hydrogen exhaustion. With 1d28_0.2v the luminosity offset rises to 0.20 dex, while the temperature excess falls to 1.1 × 10⁴ K because centrifugal support partly counteracts DM-driven contraction.

Fig. 3. Similar to Fig. 1 but with an initial rotation velocity of 20 % vcrit. The grey curve represents the non-rotating tracks.

The higher L reflects a larger convective core: rotational mixing carries fresh H inwards and He outwards, raising the nuclear luminosity that supplements WIMP heating. Consequently the envelope expands less than in the non-rotating case, keeping the surface hotter for a given ρ_χ.

Quantitatively, the helium-core mass at MS exit grows from 5.06 M_⊙ (wtt) to 5.29 M_⊙ for wtt_0.2v, and from 5.19 M_⊙ to 7.32 M_⊙ for 1d28_new and 1d28_0.2v, respectively. The factor of 1.4 increase sets the stage for still greater differences at higher spin.

3.2.2. Faster rotation: 40–60 % v_crit

The 40% v_crit (63d29_0.4v) track in Fig. 4 is ≃0.3 dex more luminous and ∼1.6 × 10⁴ K hotter than the 20% (63d29_0.2v) track at the terminal-age MS (X_c ≃ 10⁻³). The luminosity rise follows directly from the larger helium core: M_He = 13.28 M_⊙ versus 7.32 M_⊙, a factor of 1.8 increase that boosts the nuclear energy release once WIMP heating subsides. Centrifugal support lowers the effective gravity at the stellar surface by ≈12% (evaluated at the equator), reducing the photospheric temperature despite the higher L. Rotationally induced meridional circulation also transports entropy outwards, inflating the envelope by ∼20% in radius relative to the 20% model. The net result is a star that is brighter but slightly cooler, sliding upwards and rightwards in the HRD. A higher DM supply (ρ_χσ_SD = 1.0 × 10⁻²⁸ GeV cm⁻¹ (1d28_0.4v)) produces an identical differential shift.

Fig. 4. HRD of a Pop III star of 20 M⊙, from the ZAMS to the end of the MS, for a WIMP density of ρχσSD = 6.3 ⋅ 10−29 GeV ⋅ cm−1, for an initial velocity at 20% and 40% of the critical velocity of the star.

Figure 5 demonstrates that DM longevity enhances the surface spin-up driven by internal angular-momentum transport. At ρ_χ = 6.3 × 10⁹ GeV cm⁻³ the 40 %model (63d29_0.4v) attains Ω/Ω_crit = 0.99 at the end of MS, whereas the same model without DM has reached only 0.9 by that age. The higher DM density extends the Kelvin–Helmholtz spin-up phase to ∼700 Myr, doubling the time available for meridional circulation to transfer core angular momentum outwards. In the 60% sequence the stellar surface meets the break-up limit at ≃90 Myr, long before the core hydrogen fraction falls below 0.6; the computation therefore halts. This behaviour implies that any Pop III star born with v ≳ 0.5 v_crit in a DM halo of ρ_χσ_SD = 1 ⋅ 10⁻²⁸ GeV ⋅ cm⁻¹ will likely reach critical rotation well inside its DM-supported phase.

Fig. 5. Surface Ω/Ωcrit versus time for 20 (bottom), 40 (middle), and 60% (top) vcrit models at several WIMP densities, with a colour scheme similar to that of Fig. 3.

Internal rotation profiles in Fig. 6 show how DM-induced longevity enforces nearly solid-body rotation. All models present a constant rotation profile at the mid-MS (X_c = 0.37), while at the end-MS the contrast between the centre and surface is ≃15% for ρ_χσ_SD < 1.0 × 10⁻²⁸ GeV cm⁻¹ and ≲1% for the all other models, highlighting that DM energy input, by prolonging the MS by several orders of magnitude, gives rotation enough time to homogenise the interior angular momentum.

Fig. 6. Angular-velocity profile, Ω(m), at mid-MS (left) and MS end (right) for different WIMP densities and initial velocities.

3.2.3. Main-sequence lifetimes

Lifetimes in Fig. 7 underline the synergy of rotation and DM. At ρ_χσ_SD = 6.3 × 10⁻²⁹ GeV cm⁻¹ the MS lasts 23.5, 34.1, and 48.4 Myr for 0, 20, and 40% v_crit, respectively; at ρ_χσ_SD = 1.0 × 10⁻²⁸ GeV cm⁻¹ the sequence reads 50.9, 84.2, and 92.0 Myr. The absolute gain from rotation is largest at intermediate DM supply, where nuclear burning still contributes appreciably and can benefit from fresh H mixed inwards. For ρ_χσ_SD ≳ 2 × 10⁻²⁸ GeV cm⁻¹ the non-rotating lifetime already exceeds 0.7 Gyr; rotation adds only a ∼10% increment because WIMP heating so completely suppresses fusion that additional mixing has little leverage. Repeating the 20 and 40% sequences with σ_SD = 10⁻⁴⁰ cm² and proportionally larger ρ_χ reproduces the HRDs, core masses, and lifetimes to within 3% of the values above; see Table 1.

Fig. 7. MS lifetime versus ρχ for 0 (black), 20 (blue), and 40% vcrit (green). The dashed grey line is taken from T08 and describes a static 20 M⊙ star.

Rotation reinforces the effects of DM heating by (i) raising the surface luminosity at fixed ρ_χ, (ii) driving the star towards near-solid-body rotation, (iii) building substantially larger helium cores, and (iv) further extending the already prolonged MS lifetime. At sufficiently high spin and DM density the star approaches critical rotation well before core hydrogen exhaustion, a regime requiring dedicated treatment of centrifugal mass loss that lies beyond the scope of the present work.

3.3. Chemical abundance evolution in DM-supported Pop III stars

Figure 8 shows how DM heating and rotation reshape the chemical abundances in a 20 M_⊙ Pop III star with WIMPs and rotation. Four evolutionary checkpoints are shown: mid–H burning (X_c ≃ 0.37), end of H burning (X_c ≃ 10⁻³), mid–He burning (Y_c ≃ 0.50), and end of He burning (Y_c ≃ 10⁻³). The centre is fully convective so the central abundance, noted A_c with A an arbitrary element, is the same across the core. A_s is taken at the stellar surface.

Fig. 8. Abundances of the main elements inside a 20 M⊙ star with ρχσSD = 1 ⋅ 10−28 GeV ⋅ cm−1 at 20% of vcrit initially.

3.3.1. H and He abundances

For a star without DM and rotation the surface abundances remain constant during evolution with X_s = 0.75 and Y_s = 0.25. Table 2 shows that these abundances do not vary between models without WIMPs, either with rotation (wtt_0.2v) or without rotation (wtt), and the model with WIMPs but no rotation (1d28_new). The only differences appear in the model combining WIMPs and rotation, 1d30_0.2v. In the middle of the MS the gradient from core to surface decreases from ΔX = ΔY = 0.38 to ΔX = ΔY = 0.20, caused by higher He and lower H surface abundances due to homogenisation.

At higher spins (1d28_0.4v, not in Table 2), rotation alone nearly achieves full homogenisation, with ΔX = ΔY = 0.10 in the middle of MS (X_s = 0.47, Y_s = 0.53). At higher WIMP density (3d28_0.2v, not in Table 2) full homogenisation is almost reached, with ΔX = ΔY = 0.02 in the middle of MS (X_s = 0.40, Y_s = 0.60). Adding DM to rotation extends the period over which the star remains mixed. Observable effects include helium-enriched surfaces and weaker core–envelope composition contrasts, both of which alter ionising spectra and late-stage yields.

3.3.2. CNO abundance patterns and primary nitrogen

Figure 8 and Table 2 show the mass fractions of the main CNO isotopes (¹²C, ¹⁴N, and ¹⁶O); neon is not included. The figure shows a sharp CNO spike confined inside the helium core (m/M ≲ 1/3). During the MS the abundances remain below 10⁻⁸ for X_¹²C, 10⁻⁷ for X_¹⁴N, and 10⁻⁹ for X_¹⁶O, which is consistent across all initial conditions. Models without rotation keep metal-free envelopes, while model wtt_0.2v reaches X_¹²C = 2.8 × 10⁻¹⁴, X_¹⁴N = 6.6 × 10⁻¹², and X_¹⁶O = 1.3 × 10⁻¹³ at the end of the MS. Adding DM with rotation, in model 1d28_0.2v, raises the surface abundances by a factor of ≈10² at the end of the MS, giving X_¹²C = 1.33 × 10⁻¹², X_¹⁴N = 2.7 × 10⁻¹⁰, and X_¹⁶O = 4.8 × 10⁻¹².

As with H and He, stronger rotation and higher WIMP density enhance mixing and increase the CNO surface abundances. In model 1d28_0.4v the values are X_¹²C = 4.4 × 10⁻¹¹, X_¹⁴N = 6.9 × 10⁻⁹, and X_¹⁶O = 9.6 × 10⁻¹¹. In model 3d28_0.2v the values are X_¹²C = 2.6 × 10⁻¹¹, X_¹⁴N = 4.3 × 10⁻⁹, and X_¹⁶O = 6.3 × 10⁻¹¹. Higher rotation mixes the heavier elements (¹²C, ¹⁴N, and ¹⁶O) more efficiently than higher WIMP density, which instead enhances the mixing of H and He by producing lower gradients for these elements.

3.3.3. Supernova remnant yields

We did not push the computations until the end of core Si burning to better produce remnant masses, and we are aware that physics of calculating remnant masses is uncertain. However, we can nonetheless provide a first-order estimate on the remnant masses using the two formulations widely available in the literature. The first is Maeder (1992), where they find a remnant mass of 2.14 M_⊙ for an initial non rotating star of 20 M_⊙ (their Table 4). They computed this mass by considering that all the skin layers have been ejected by the supernova explosion. Second one is Patton et al. (2022), they use non rotating models but take binary interactions into account. They find a mass of ≈5.7 M_⊙ (Fig. 1). The discrepancy between the two values comes from different methodology to determine the remnant mass, and highlight the difficulty in estimating these values.

In this work we computed the mass of each element expected in the future supernova remnant (SNR) yields. For this we integrated the mass at the end of He burning (Y_c​ ≃ ​10⁻³) from the edge of the helium core (X_¹⁴C < 10⁻³) to the surface for each element. The results are summarised in Table 3. The integrated masses of H and He are ∼9 M_⊙ and ∼6 M_⊙, except for model 1d28_0.2v where the values are ∼5 M_⊙ and ∼8 M_⊙. This corresponds to an external layer that is 45% poorer in H and 33% richer in He. Nitrogen follows a similar trend, with a mass five times higher for 1d28_0.2v. Oxygen yields are two orders of magnitude larger in rotating stars, wtt_0.2v and 1d28_0.2v (> 10⁻⁶ compared to ∼3 × 10⁻⁸), and in 1d28_0.2v they are eight times higher than in the DM-free rotating model wtt_0.2v.

For the other elements the trends are less clear. The highest carbon yield is in model wtt (7.3 × 10⁻⁵) and the lowest is in wtt_0.2v (1.8 × 10⁻⁵). The trend reverses in DM models, but the difference is negligible. Neon yields are one order of magnitude higher in 1d28_0.2v (4.7 × 10⁻¹⁰) compared to wtt_0.2v (2.7 × 10⁻¹¹), while they remain relatively similar in the static models wtt and 1d28 (9.6 × 10⁻¹⁰ and 2.7 × 10⁻¹⁰).

If we consider that the remnant mass will be the remaining helium core while the surface layers are ejected by the supernova explosion (similarly as in Maeder 1992), we get 4.8 M_⊙ for wtt, 5.1 M_⊙ for 1d28_new, 5.0 M_⊙ for wtt_0.2v, and 7 M_⊙ for 1d28_0.2v. The masses are more similar to Patton et al. (2022) than Maeder (1992), and with a remnant mass ≈40% higher for the star with DM and rotation. However, these results are to be taken carefully as we do not compute the models until pre-supernova stage.

In conclusion, the chemical composition of SNRs may provide a way to probe DM in Pop III stars. The rotating model with WIMPs is the most promising as it shows increased He, N, and O yields together with lower H. Carbon and neon are not good tracers because their SNR abundances differ little from DM-free models. The remnant masses could be a way to detect those objects as well, albeit more work needs to be done regarding this topic.

3.3.4. DM-induced mixing of H and He

The integrated masses of H and He are ∼9 M_⊙ and ∼6 M_⊙, respectively, except for 1d28_0.2v where the integrated masses are ∼5 M_⊙ and ∼8 M_⊙, respectively, i.e. the external part of the star is 45% poorer in H and 33% richer in He. N has a similar trend with a mass 5 times higher for 1d28_0.2v. The oxygen yields are two orders of magnitude more massive for rotating stars, wtt_0.2v and 1d28_0.2v (> 10⁻⁶ compare to ∼3 ⋅ 10⁻⁸), especially for 1d28_0.2v, which is 8 times higher than the DM-less rotating star wtt_0.2v.

For the two other elements, finding a general trend is more difficult. The highest carbon yields is present for wtt (7.3 ⋅ 10⁻⁵) while the lowest is present for wtt_0.2v (1.8 ⋅ 10⁻⁵). The opposite is true when the star is surrounded by DM, but with a negligible difference. The Neon yield is one order of magnitude higher for 1d28_0.2v model (4.7 ⋅ 10⁻¹⁰) compare to wtt_0.2v (2.7 ⋅ 10⁻¹¹), while they are relatively similar for the static stars, wtt and 1d28 (9.6 ⋅ 10⁻¹⁰ and 2.7 ⋅ 10⁻¹⁰).

In conclusion, a way to look for DM in Pop III stars is to look at the chemical composition of the SNR from those stars. The rotating model with WIMPs is the most promising as it exhibits an increase in He, N and O, while decreasing for H. C and Ne do not seem to be good tracers to detect stars surrounded by DM as their proportion in SNR do not differ significantly from the DM-less models.

4. Discussion

The 20 M_⊙ mass prescription chosen in this work is a suitable candidate for type II supernova (core collapse supernova). Additionally, 20 M_⊙ is a classical choice when looking at initial mass functions not only for in local universe but also for high redshift universe, this means if stars powered by DM exist, it would be a much higher chance to observe such a star (ξ(m)Δm ≈ 6 × 10⁻³). Lastly, the lifetime of this star without any DM is long enough for any rotation processes to be significant (≈9 Myr for static and ≈10 Myr for rotating star (Ekström et al. 2012). As an example, if we were to take a 60 M_⊙ star, its lifespan would decrease to a few million years and the chemical mixing would be less significant. A rotating 20 M_⊙ star, whose internal structure is composed of a convective core and a radiative envelope, does not produce a homogeneous evolution. However, adding DM to this star, we indeed find homogeneous evolution due to lifetime extension, and allowing the transport of chemicals even through the radiative zones. In addition, the likelihood of finding a 20 M_⊙ star, despite there being hints of a top heavy initial mass function for Pop III stars, is much higher than, for instance, 60 M_⊙ stars (Chen et al. 2020). This makes the 20 M_⊙ model a viable case for comparison in this work.

The annihilation of WIMPs reshapes the internal composition of our 20 M_⊙ grid. The non-rotating control model builds a steep H–He discontinuity with ΔX ≃ 0.38 at m/M ≃ 0.22 and retains a pristine envelope at X ≈ 0.75. When DM heating is included with v_ini = 0.40 v_crit the star evolves quasi-homogeneously and the gradient falls to ΔX ≃ 0.10 at H exhaustion. The surface nitrogen mass fraction rises from 6.6 × 10⁻¹² to 2.7 × 10⁻¹⁰ at the end of MS. These values exceed the primary nitrogen reported by Yoon et al. (2008) for the same mass. The non-rotating DM case shows weaker enrichment and thus confirms that DM and rotation together drive the extreme CNO yields.

Lifetimes scale with the product ρ_χσ_SD as anticipated by Taoso et al. (2008). The reference model without DM leaves the MS after 9.6 Myr. At ρ_χ = 6.3 × 10⁹ GeV cm⁻³ the lifetime doubles to 23.5 Myr, matching the 20 M_⊙ result of Scott et al. (2009) to within 5%. Raising the density to 1.0 × 10¹⁰ extends the lifetime to 50.9 Myr, i.e. by a factor of 5.3, which agrees with the 4.9 reported by Casanellas & Lopes (2011a). Beyond 2 × 10¹⁰ our models fail to exhaust hydrogen within a Hubble time, reproducing the threshold identified by Taoso et al. (2008). Switching to σ_SD = 10⁻⁴⁰ cm² and scaling ρ_χ by 100 recovers identical lifetimes, confirming that τ_MS depends only on the capture product.

Rotation changes the outcome only when DM prolongs the nuclear timescale. A v_ini = 0.40 v_crit model without DM ends core-H burning with Ω_core/Ω_surf ≃ 3 and M_He = 6.8 M_⊙. Introducing ρ_χ = 1.0 × 10¹⁰ GeV cm⁻³ lowers the ratio to 1.08 (8%) and raises the core mass to 13.3 M_⊙. These values fulfil the quasi-chemically homogeneous condition τ_mix < τ_nuc proposed by Yoon et al. (2008). With DM, the same state appears at 0.20 v_crit. This halves the chemically homogeneous evolution threshold compared with the ≳0.50 v_crit requirement in DM-free Pop III models. We conclude that WIMP energy support, rather than centrifugal mixing alone, dictates the onset of chemical homogeneity in primordial massive stars.

Figure 2 shows the central temperature T_c versus ρ_χσ_SD for a rotating 20 M_⊙ Pop III model with WIMP annihilation, with stages defined by X_c and dashed reference curves from T08. At the ZAMS the WIMP model is 26% cooler and 55% larger with T_eff ≈ 50 kK and R ≈ 2.8 R_⊙ compared to the DM–free model with T_eff ≈ 63 kK and R ≈ 1.8 R_⊙. This inflation follows from annihilation heating, which supplies part of the luminosity and lowers the nuclear rate, and T_c, which in turn increases the core entropy. At later times the annihilation rate drops and static models differ by only 2–8% in radius and by less than 3% in temperature so their structures converge. In the rotating WIMP model the trend reverses after the early MS and the radius becomes smaller by up to 17% while T_eff becomes higher by 14–20%. Rotation maintains strong mixing and flattens the μ gradient, which enriches the envelope in He and lowers the opacity. Therefore, the T_eff rises at fixed luminosity and the star contracts. The core remains cooler throughout, because mixing spreads the annihilation power and throttles nuclear burning.

5. Conclusions

We have computed the first grid of 20 M_⊙ Pop III models that self-consistently include WIMP capture, annihilation, and rotation. The grid spans six ambient DM densities (10⁸–3 × 10¹⁰ GeV cm⁻³) and three initial spins (0, 0.20, and 0.40 v_crit), enabling a controlled comparison of non-rotating, rotating, DM-free, and DM-rich tracks.

• Main-sequence longevity. WIMP annihilation lengthens the lifetime of a 20 M_⊙ Pop III star from ∼10 Myr to gigayear scales once the capture product exceeds ρ_χσ_SD ∼ 2 × 10⁻²⁸ GeV cm⁻¹. The vastly slower nuclear clock lets rotation and mixing act for orders of magnitude longer, reshaping subsequent feedback on the host halo.

• Rotation driven towards solid body and helium-core growth. In DM-rich tracks, a model born at 0.40 v_crit ends core-H burning with only a few-percent shear and a helium core roughly twice the mass of its DM-free counterpart. Near-rigid rotation suppresses shear instabilities yet sustains meridional circulation, enhancing fuel transport and setting the stage for energetic late phases.

• Central structure stabilisation. Continuous DM heating holds T_c and ρ_c nearly constant during most of the hydrogen burning, delaying Kelvin–Helmholtz contraction until spin-dependent capture ceases. This stabilisation maintains a low-density core that favours efficient mixing and prolongs the homogeneous phase.

• Chemically homogeneous evolution at moderate spin rates. DM support lengthens τ_nuc such that Eddington–Sweet circulation overwhelms it even at 0.20 v_crit, while DM-free stars need ≳0.50 v_crit. Quasi-homogeneity flattens the H–He profile, keeps fresh fuel in the core, and produces progenitors suitable for collapsar-type explosions.

• Surface H/He and CNO patterns with primary nitrogen. Rotation and DM mixing lowers the surface hydrogen to X ∼ 0.53, raises helium to Y ∼ 0.47, and boosts ¹²C, ¹⁴N, and ¹⁶O by two orders of magnitude relative to classical Pop III stars.

• Spectral signatures. The best chance we have of detecting DM is to study rotating Pop III stars believed to be surrounded by DM. They live longer than their static counterparts and will produce SNRs rich in He, N, and O and poor in H; this is due to modified outer layer abundances that are in turn due to more efficient mixing. Additionally, these stars begin their evolution as cool and bloated objects at the ZAMS but become hotter and more compact throughout the evolution compared to classical baryonic Pop III stars.

Future work will expand the mass range to 40–120 M_⊙ and probe halo densities above 3 × 10¹⁰ GeV cm⁻³. Spin-independent capture on metals will be added to follow late evolutionary phases. Magnetic torques and mass loss will be coupled to test whether DM-supported fast rotators reach critical spin and launch chemically enriched winds. Synthetic spectra across the UV–IR range will quantify N III and He II diagnostics and calibrate ionising photon budgets for re-ionisation models. Yield calculations will track ²²Ne and heavy-element production, which will be compared with N-rich extremely metal-poor stars. Finally, models of 10³–10⁵ M_⊙ protostars with sustained DM heating will explore the formation of supermassive DSs and the seeds of early quasars, which could be observed by the JWST (see Ilie et al. 2023a for an analysis of recent observations of quasars that could be DS products).

Acknowledgments

The authors would like to thank Dr. Marco Taoso and Prof. Georges Meynet for their valuable comments and feedback on the aspects of DM annihilation and rotation. Anaïs Pauchet would like to thank her supervisor Pr. Stefanie Walch for her supervision and encouragement towards the project. AP was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 500700252 – SFB 1601. DN was supported by the Swiss National Science Fund (SNSF) Postdoctoral Fellowship, grant number: P500-2235464 and by the Virginia Initiative of Theoretical Astronomy (VITA) Fellowship.

References

All Tables

All Figures

	Fig. 1. HRD from this work (solid lines) and T08 (dashed lines) for a static 20 M⊙ Pop III star surrounded by different WIMP densities (see the colour code).
In the text

	Fig. 2. Central temperature as a function of WIMP density for a 20 M⊙ Pop III star without rotation (left) and with a rotation 20% of vcrit (right). Curves corresponding to different evolutionary stages (defined by the central hydrogen abundance) are plotted with spline interpolation; the dashed curves are taken from T08.
In the text

	Fig. 3. Similar to Fig. 1 but with an initial rotation velocity of 20 % vcrit. The grey curve represents the non-rotating tracks.
In the text

	Fig. 4. HRD of a Pop III star of 20 M⊙, from the ZAMS to the end of the MS, for a WIMP density of ρχσSD = 6.3 ⋅ 10−29 GeV ⋅ cm−1, for an initial velocity at 20% and 40% of the critical velocity of the star.
In the text

	Fig. 5. Surface Ω/Ωcrit versus time for 20 (bottom), 40 (middle), and 60% (top) vcrit models at several WIMP densities, with a colour scheme similar to that of Fig. 3.
In the text

	Fig. 6. Angular-velocity profile, Ω(m), at mid-MS (left) and MS end (right) for different WIMP densities and initial velocities.
In the text

	Fig. 7. MS lifetime versus ρχ for 0 (black), 20 (blue), and 40% vcrit (green). The dashed grey line is taken from T08 and describes a static 20 M⊙ star.
In the text

	Fig. 8. Abundances of the main elements inside a 20 M⊙ star with ρχσSD = 1 ⋅ 10−28 GeV ⋅ cm−1 at 20% of vcrit initially.
In the text