The impermanent fate of massive stars in AGN discs

ABSTRACT

Stars are likely to form or to be captured in active galactic nucleus (AGN) discs. Their mass reaches an equilibrium when their rate of accretion is balanced by that of wind. If the exchanged gas is well mixed with the stellar core, this metabolic process would indefinitely sustain an ‘immortal’ state on the main sequence (MS) and pollute the disc with He byproducts. This theoretical extrapolation is inconsistent with the super-solar α element and Fe abundances inferred from the broad emission lines in AGNs with modest He concentration. We show this paradox can be resolved with a highly efficient retention of the He ashes or the suppression of chemical blending. The latter mechanism is robust in the geometrically thin dense sub-pc regions of the disc where the embedded-stars’ mass is limited by the gap-formation condition. These stars contain a radiative zone between their mass-exchange stellar surface and the nuclear-burning core. Insulation of the core lead to the gradual decrease of its H fuel and the stars’ equilibrium masses. These stars transition to their post-MS (PostMS) tracks on a chemical evolution time-scale of a few Myr. Subsequently, the triple-α and α-chain reactions generate α and Fe byproducts which are released into their natal discs. These PostMS stars also undergo core collapse, set off type II supernova, and leave behind a few solar-mass residual black holes or neutron stars.

1 INTRODUCTION

1.1 Background and motivations

Studying active galactic nuclei (AGNs) is fundamental in understanding the formation and evolution of galaxies and their central supermassive black holes (SMBHs; Kormendy & Ho 2013). These very bright objects are powered by the release of gravitational energy through viscous dissipation in gaseous accretion discs around SMBHs with masses |$M_\bullet \sim 10^{4-9} \ \mathrm{M}_\odot$| (Lynden-Bell 1969; Rees 1984).

Based on optical and x-ray data (Soltan 1982; Fabian & Iwasawa 1999; Elvis, Risaliti & Zamorani 2002), population synthesis models (Yu & Tremaine 2002; Marconi et al. 2004; Shankar et al. 2004; Shankar, Weinberg & Miralda-Escudé 2009) infer average values for the Eddington factor [λ_• ≡ L_AGN/L_E• ≃ 0.6, although see Raimundo & Fabian (2009) who proposed the existence of three distinct AGN populations with λ_• ranging from 0.05 to 0.6 as a function of their obscuration and ϵ_•] and efficiency factor (⁠|$\epsilon _\bullet \equiv L_{\rm AGN} /{\dot{M}}_\bullet c^2 \simeq 0.06$|⁠) where L_AGN, L_⊙, L_E• = M_•L_E⊙/M_⊙, and L_E⊙ = 3.2 × 10⁴ L_⊙ are the AGN, solar, Eddington luminosity for SMBH and the Sun, respectively. The accretion rate |${\dot{M}}_\bullet \simeq 2 f_{\bullet } m_8\ {\ \mathrm{M}_\odot \, {\rm yr}^{-1}}$| where |$m_8 \equiv M_\bullet /10^8 \ \mathrm{M}_\odot$| and f_• ≡ (λ_•/0.6)(0.06/ϵ_•).

In a steady state, SMBHs’ |${\dot{M}}_\bullet$| is fed by their surrounding discs with same accretion rate |${\dot{M}}_{\rm d}$|⁠. The structure of these discs is determined by the efficiency of angular momentum transfer (Lynden-Bell & Pringle 1974; Pringle 1981). With an α_ν prescription (Shakura & Sunyaev 1973) for viscosity (ν = α_νc_sH where α_ν is an efficiency factor, |$\Omega =\sqrt{G M_\bullet /R^3}$| is the Keplerian frequency), the radial distribution of the gas surface density Σ_g = 2ρ_cH, midplane density ρ_c, temperature T_c, sound speed c_s, disc thickness H ≃ c_s/Ω, aspect radio h = H/R, surface temperature T_e, and radiative flux |$Q^- = 2 \sigma T_{\rm e}^4$| can be determined (Frank, King & Raine 2002; for either dust κ_dust or electron scattering κ_es) opacity under the assumption of thermal equilibrium in which Q⁻ is balanced by the rate of viscous dissipation |$Q^+ _\nu = 9 \Sigma _{\rm g} \nu \Omega ^2/4$|⁠.

With the conventional α_ν prescription and the inferred |${\dot{M}}_{\rm d} (M_\bullet)$|⁠, the outer regions (at radii R ≡ r_pc pc ≳ 10⁻² pc) of discs around SMBHs with m₈ ∼ 0.01–1 are prone to gravitational instability GI (Paczynski 1978) which leads to gravito-turbulence with α_ν ∼ 0.1–1 (Lin & Pringle 1987; Lin, Pringle & Rees 1988; Kratter, Matzner & Krumholz 2008; Zhu et al. 2010a; Zhu, Hartmann & Gammie 2010b; Martin & Lubow 2011; Rafikov 2015). GI also leads to prolific rates of in situ star formation (Goodman 2003; Goodman & Tan 2004; Thompson, Quataert & Murray 2005) and to capture from the omnipresent nuclear clusters (Syer, Clarke & Rees 1991; Artymowicz, Lin & Wampler 1993; Davies & Lin 2020). Here we introduce and adopt a generic stellar evolution and pollution in AGN discs (SEPAD) model under the assumption that these processes lead to a self-regulated marginally stable state with the gravitational |$Q (\simeq {c_{\rm s} \Omega / \pi G \Sigma _{\rm g}} = {M_{\bullet } / {\sqrt{8}} \pi \rho _{\rm c} R^3})$| value ∼1 (Safronov 1960; Toomre 1964).

For the SEPAD model, we assume the embedded stars provide an auxiliary power |$Q^+ _\star (\gt Q^+ _\nu)$| generated by the conversion of hydrogen H into helium He during their main sequence (MS) evolution and that of He into α element along their post MS tracks. The deduced h(R) (equation 2) and Q⁻(R) (Goodman 2003; Thompson et al. 2005) for the modified thermal equilibrium (with |$Q^+ _\star \simeq Q^-$| in Section 4.2) are consistent with those inferred from the reverberation mapping data (Starkey et al. 2022), the microlensing light curve (Pooley et al. 2007; Morgan et al. 2018; Cornachione et al. 2020), and the observed spectral energy distribution, especially over the infrared wavelength range(Sanders et al. 1989). A direct implication of the SEPAD scenario is that embedded stars also return He, α, and Fe byproducts to the AGN discs and these pollutants are mixed with gas in the main-disc flow towards the SMBH. Consequently, the contaminants do not accumulate with the AGNs’ evolution and the discs’ metallicity is independent of z_γ (Artymowicz et al. 1993).

In a recent work, Cantiello, Jermyn & Lin (2021) and Jermyn et al. (2022) studied the evolution of stars embedded in an AGN disc with the mesa code (Paxton et al. 2011, 2013, 2015, 2018, 2019). Their simulations include both mass accretion and stellar wind loss with rates |${\dot{M}}_{\rm Bondi}$| and |${\dot{M}}_{\rm wind}$| respectively. They found that, within the inner few pc’s, where the density in the disc midplane is relatively high (ρ_c ≳ 6 × 10⁻¹⁸g cm⁻³), the Bondi accretion time-scale (⁠|$\tau _{\rm B} \equiv M_\star / {\dot{M}}_{\rm Bondi}$|⁠) is much shorter than the nuclear burning time-scale. Consequently, newly formed or trapped stars undergo fast mass growth. As their mass |$M_\star (\equiv m_\star \ \mathrm{M}_\odot)$| increases to |$\sim 10^3 \ \mathrm{M}_\odot$|⁠, their nuclear fusion (through the CNO cycle) generates a luminosity L_⋆ ∼ L_E⋆ ≡ m_⋆L_E⊙ (their Eddington limit), under which radiation pressure suppresses the accretion rate and intensifies the stellar wind. Cantiello et al. (2021) further assume that the accretion and wind lead to inflow of pristine H-rich disc gas and outflow of He-rich byproducts, albeit the stars’ mass growth is stalled when these processes reach an equilibrium. In addition, they assume that the fresh H replenishment and He byproduct are fully mixed throughout the stellar interior, including their radiative zones (Section 3.1). Under these two conditions, massive stars retain their mass M_eq and composition, and perpetually remain on the MS as ‘immortal’ stars until the AGN phase is terminated with a severe depletion of the disc gas. In the stellar core, H is converted into He through the CNO cycle at a rate |${\dot{M}}_{\rm He} \simeq L_\star /\epsilon _{\rm He} c^2$| where ϵ_He ≃ 0.007 is the H-to-He fusion efficiency. Multiple stars coexist to maintain a thermal equilibrium (⁠|$Q^- \simeq Q^+ _\star \gt \gt Q^+ _\nu$|⁠) with a surface density s_⋆ ≃ Q⁻/L_⋆ ≃ Q⁻/m_eqL_E⊙. The ‘immortal’ stars’ wind return their He ashes to the disc at the same rate (⁠|${\dot{M}}_{\rm He}$|⁠). This feedback process increases the He surface density in the disc at an average rate up to |${\dot{\Sigma }}_{\rm He} = s_\star {\dot{M}}_{\rm He} = Q^- / \epsilon _{\rm H} c^2$| and the He mass fraction in the disc gas to |${\it Y}_{\rm d} \simeq 2 \pi \int {\dot{\Sigma }}_{\rm He} R \mathrm{ d}r/ {\dot{M}}_{\rm d}$|⁠. Although the CNO burning process markedly increases the abundance ratio of N/(C+O) from that produced by the triple-α and α-chain reactions during the PostMS evolution of some previous generation stars, it does not lead to changes in the star’s α (mostly C+N + O) and Fe mass fraction (Z_α and Z_Fe) from their pre-existing values, as long as the high local and global disc density is maintained. Consequently, the perpetual massive MS stars do not lead to increases in the ongoing-AGN discs’s α and Fe mass fraction (⁠|${\mathcal {Z}}_{\alpha }$| and |${\mathcal {Z}}_{\rm Fe}$|⁠).

These predictions can be tested by the observational inferred metallicity of AGN discs. With the aid of the cloudy models for emission-line physics, the magnitude of Y_d, |${\mathcal {Z}}_{\alpha }$|⁠, and |${\mathcal {Z}}_{\rm Fe}$| of the disc gas can be inferred from the ratios of the measured strength of various broad emission lines associated with AGNs (Osterbrock & Ferland 2006) at a range cosmological red-shifts, z_γ. An expository survey of spectroscopic data associated with AGNs’ broad-line regions (BLRs; Huang, Lin & Shields 2023) indicates that (i) Y_d is |$\sim 30{{\ \rm per \, cent}}$| higher than its solar value (Osterbrock & Shuder 1982; Shang et al. 2007), (ii) α (C+N + O) element abundance |${\mathcal {Z}}_{\alpha }$| is ∼3–5 that of the solar value (Nagao, Marconi & Maiolino 2006a; Temple et al. 2021; Lai et al. 2022), (iii) N/(C + O) is larger than the solar value and increases with |${\mathcal {Z}}_{\alpha }$| (Hamann et al. 2002), and (iv) Fe mass fraction |${\mathcal {Z}}_{\rm Fe}$| is comparable to its solar value (Dietrich et al. 2003; Maiolino et al. 2003; Wang et al. 2022).

The modest Y_d elevation and the substantial |${\mathcal {Z}}_{\alpha }$| and |${\mathcal {Z}}_{\rm Fe}$| enhancements contradict the expected predominantly He yield from a population of hypothetically ‘immortal’ massive MS stars which do not evolve onto their post-MS (PostMS) phase until either the disc is nearly hydrogen free with |${\it Y}_{\rm d} \sim {\mathcal {O}} (1)$| (Jermyn et al. 2022) or the stars’ accretion rate is extensively suppressed by the severe depletion of the disc gas (Cantiello et al. 2021). The former scenario is inconsistent with the observed modest value of Y_d whereas the latter possibility is inconsistent with the observationally inferred super-solar |${\mathcal {Z}}_{\alpha }$| and nearly solar |${\mathcal {Z}}_{\rm Fe}$| for ongoing AGNs, independent of the cosmological redshift.

1.2 Main goal: ending the MS

In order to resolve the theoretical and observational tensions, we re-visit stellar evolution in the context of the SEPAD scenario. Although stellar-evolution models themselves do not depend on the SEPAD assumptions, they can provide useful constraints on the SEPAD models in terms of both the disc-environment boundary conditions and the implications on the He, α and Fe byproducts. The primary motivation in this paper is to investigate some potential channels for the massive evolving stars to transition from the MS to PostMS phase. This transition requires that the stellar hydrogen mass fraction X_⋆ be severely depleted in the hydrogen burning zone. Since the mass content of these embedded stars is recurrently reprocessed by accretion-wind metabolism, in contrast to the conventional standalone stars, it is essential to suppress the resupply of fresh hydrogen fuel into the nuclear furnace.

In light that the modest Y_d inferred from observational data, we explore two potential interruption mechanisms for the hydrogen supply chain:

1. an effective retention of the helium byproducts carried by the stars’ wind between their Bondi R_B and Hills R_H radii, or

2. the isolation of the CNO burning core from the outer mass-exchange surface by a mixing-free radiative layer in the stellar interior.

The second scenario is attainable provided a mass-growth limit (⁠|$\lesssim 6-700 \ \mathrm{M}_\odot$|⁠) is imposed by the suppression of stellar accretion rate due to the formation of tidally induced gaps in the disc.

We carry out a series of numerical simulation with the mesa code (Section 2). Results of these simulations are presented and analyzed in Section 3. We examine the feedback of the mass loss process on the chemical composition of the accreted material, evaluate its effectiveness in enabling the MS-to-PostMS transition, and determine the duration of the MS evolution τ_MS.

1.3 Secondary goals

1.3.1 Quantifying evolving stars’ PostMS evolution

The evolving stars’ transition from MS to PostMS evolution (Section 4) not only activates the production of α elements through the triple-α and α-chain reactions but also provides a robust channel for the in situ production of solar or supersolar |${\mathcal {Z}}_{\rm Fe}$| with subsolar Fe/α through (core-collapse) type II supernovae (SN II). In the hot cores of very metal-deficient (with Log Z_α/Z_⊙ ≲ −1.75 and presumably tenuous mass-loss rates), relatively massive (⁠|$\gtrsim 10^2 \ \mathrm{M}_\odot$|⁠) PostMS stars, prolific pair production excites instabilities which trigger super novae and eject a large amount of α, and Fe remnants without residual black holes (Spera & Mapelli 2017; Woosley 2017). These instabilities are suppressed by substantial mass losses in stand-alone stars, although the stability metallicity threshold is mass and model dependent, with estimates ranging from log Z_α/Z_⊙ ≳ −1.75 to ≳ −0.47 (Heger & Woosley 2002; Langer et al. 2007; Chatzopoulos & Wheeler 2012a, b).

As He is converted to α and Fe through triple-α and α-chain reaction on their PostMS track, the magnitude of He mass fraction Y_⋆ inside the stars reduces while both Z_α and Z_Fe increase. These byproducts are returned to the disc through PostMS winds and they contribute to the BLR’s observed Y_d, |${\mathcal {Z}}_\alpha$|⁠, and |${\mathcal {Z}}_{\rm Fe}$|⁠. The intensity of the stellar wind also determines the mass of the pre-collapse cores, the possibility of SN II, and the mass of remnant black hole remnants rBHs. Along with SN II, accretion onto rBHs also contributes to |$Q^+ _\star$|⁠. The rBHs also provide seeds for intense VIRGO/LIGO gravitational wave events through mergers.

A second goal of this investigation, with the aid of the mesa code (Section 4.1), is to estimate embedded PostMS stars’ mass and luminosity evolution, the amount of He, α, and Fe yield returned to the disc, pre-collapse core mass M_pc, their duration Δτ_PostMS, and stars’ life span τ_⋆. Under the assumption of uninterrupted state of thermal equilibrium (⁠|$Q^+ _\star \simeq Q^-$|⁠), we estimate the rate of stars formation per unit surface area to be ∼s_⋆/τ_⋆ and that of the heavy element pollution rate for the disc. We also set upper limits on the mass range of emerging rBHs after the SN II to be M_pc. Since they are likely to be retained by the disc, these rBHs may gain substantial mass through subsequent gas accretion and multiple merger events.

1.3.2 Multiple generations of evolving stars

If embedded stars evolve through the MS-PostMS-SN II or collapse impermanent life cycle, accumulative pollution from multiple generation of stars is also needed to elevate Y_d, |${\mathcal {Z}}_{\alpha }$|⁠, and |${\mathcal {Z}}_{\rm Fe}$| to their observationally inferred values. Moreover, these abundances are larger in the BLR than those in the narrow-line region (Nagao, Maiolino & Marconi 2006b; Terao et al. 2022). Their independence of z_γ(∼ 0–7) (Nagao et al. 2006a; De Rosa et al. 2011; Xu et al. 2018) suggests that the rich α and Fe contents produced by the embedded stars not only pollute the discs but are also carried by the gas flow and consumed by the central SMBH. These inferences support the SEPAD’s basic assumption of in situ chemical evolution (Section 1.1).

Our final objective is hence to gather evidences for the in situ formation of multiple generation of evolving stars. Based on our numerical models and analysis, we show, in Section 4.4, that the enhanced N/(C+O) ratio implies secondary processes, i.e. the α (C+N + O) production from the triple-α burning during the PostMS phase of an initial-generation stars must be reprocessed through the CNO cycle during the MS stage of next-generation stars. Finally, we discuss some model uncertainties in Section 5.1 and summarize our results and discuss some observable tests in Section 5. We note that the symbols and variables used throughout this paper are all defined in Table E1.

2 METHODS

2.1 AGN disc structure

The SEPAD model is constructed with the widely adopted α_ν-disc prescription (Shakura & Sunyaev 1973) for effective viscosity ν(= α_νc_sH). In a steady state, the mass flux in disc is given by

(Goodman 2003). With Q ∼ 1 and α_ν ∼ 0.1–1 for gravito-turbulence being marginally maintained, we find

The magnitude of h can be observational measured from the reverberation mappings of fluctuating lampposts within a few gravitational radii R_• ≡ GM_•/c². In exposed outer regions of the disc, the aspect ratio of the photo-surface for the reprocessed radiation is typically a few h (Garaud & Lin 2007) and that measured (∼10⁻²) for NGC 5548 (Starkey et al. 2022) at r_pc ∼ 10⁻² (a few light days) is consistent with equation (2).

2.2 mesa numerical model

For all of our simulations we use the same modified mesa code used and provided by Cantiello et al. (2021). We hence define the classic Bondi accretion rate as

This estimate does not include the effect of stars’ tidal interaction with the disc gas which leads to gap formation and reduction of the accretion rate (see further discussion in Sections 3.3 and  C).

For trans-Eddington luminosities, the wind suppresses accretion resulting in a reduced Bondi mass-accretion rate that can be written parametrically as

where L_⋆ and L_E⋆(= m_⋆L_E⊙) are the stellar luminosity and its Eddington limit. We adopt the same ‘smoothing’ prescriptions as in Cantiello et al. (2021) for comparison purpose and it is an approximation of the modified |${\dot{M}}_{\rm Bondi}$| of equation (A3) in Appendix  A. The mass loss rate is also given parametrically by

with an escape velocity V_⋆ = (2GM_⋆/R_⋆)^1/2 from the stellar surface at radius R_⋆ [this prescription is equivalent to the analytic approximation for |${\dot{M}}_{\rm wind}$| in equation (B1) in Appendix  B]. Using these prescriptions and by setting the AGN density to ρ_c = 6 × 10⁻¹⁷ g cm⁻³ and sound speed to c_s = 10⁶ cm s⁻¹, the star initially grows from |$1 \ \mathrm{M}_\odot$| to an accretion-wind equilibrium state with an asymptotic equilibrium mass |$M_{\rm eq}~\sim 630 \ \mathrm{M}_\odot$|⁠.

We find that the AGN density parameter does not affect our results qualitatively as long as it is in the ‘runaway accretion’ regime, reported by Cantiello et al. (2021) to be ρ_c ≥ 6 × 10⁻¹⁸ g cm⁻³. In a disc with marginal stability (i.e. Q ∼ 1) around |$10^8 \ \mathrm{M}_\odot$| SMBHs, this density includes all the region |$\lesssim 2.3 m_8 ^{1/3}$| pc (equation 3). In contrast, typical radial extent of the BLR inferred from reverberation mapping data gives a BLR radius ∼10–40 light days (Bentz et al. 2013; Horne et al. 2021). However, the formation of gaps in the geometrically thin, dense (equations 2 and 3) sub-pc disc region near embedded stars’ orbits may severely quench the gas supply from the disc and impose growth limits on their masses (see further discussion in Sections 3.3 and Appendix  C).

Throughout this paper, we follow Cantiello et al. (2021) in assuming that the Eddington luminosity is given by L_E⋆ = m_⋆L_E⊙, hence ignoring its dependency on the hydrogen abundance through the electron-scattering opacity: |${L_{\rm E}}_\star = 4\pi GM_\star c/(0.2\, (1 + X_\star))$|⁠. This approximation does not introduce any significant qualitative effects on our results when we reran the simulations using the more self-consistent opacity definition.

Finally, we emphasize that all of our simulations were done with solar metallicity and composition for the accreted disc gas, in contrast with Cantiello et al. (2021) who used metal-free gas. The importance of this assumption is discussed in Section 4.1.

3 TRANSITION FROM MS TO POSTMS TRACKS FOR THE EVOLVING STARS

3.1 Isolating the stellar core

In Cantiello et al. (2021), mixing throughout the radiative zone (if it is present) is imposed under the assumption that its efficiency equals to that in the convective region. This prescription was justified by noting that with the density and temperature distributions being nearly those of a γ = 4/3 polytrope, the stars’ internal structure is prone to rotational instability and meridional circulation. This extra mixing in principle can allow the freshly accreted material to cross the radiative zone and to reach the core. It also enables the nuclear byproducts in the core to be dredged up to their outer envelope. Without this extra blending, the radiative zone becomes a buffer which insulates the nuclear furnace in the stellar core from the (accretion-wind) mass-exchange zone near the stellar surface. Here we relax this assumption and investigate the role of this parameterized mixing, in interrupting the supply of hydrogen to the core and enabling the MS-to-PostMS transition.

In Fig. 1 we show the evolution tracks for models starting from 1 M_⊙ star with a solar composition, for different AGN disc densities. In all but one case, the radiative zone mixing has been turned off, allowing only for mixing in the convective zone. For the nominal (ρ_c = 6 × 10⁻¹⁷ g cm⁻³) model, the star’s evolution is nearly identical to the runaway cases in Cantiello et al. (2021). This star quickly gains mass on a time-scale ≲ 1 Myrs and reaches a steady state with M_⋆ ∼ 630 M_⊙ and the radius R_⋆ ∼ 21.5 R_⊙. The stellar luminosity L_⋆ increases from 1 L_⊙ to 10^7.5 L_⊙. This asymptotic L_⋆ is the Eddington luminosity of a star with this mass (L_E⋆ = m_⋆L_E⊙), and hence this state represents the star’s equilibrium point.

A slightly larger steady-state value of |$M_{\rm eq}~(\sim 810 \ \mathrm{M}_\odot$|⁠) is obtained in the ρ_c = 2 × 10⁻¹⁷ g cm⁻³ model. This difference is mostly due to the difference in the boundary condition near the surface of the accreting stars. Fractional M_eq-deviation between these and previous (Cantiello et al. 2021) models is due to the solar versus metal-free composition adopted for the disc gas.

In a lower-ρ_c (1 × 10⁻¹⁷ g cm⁻³) model, the star grows to |$\sim 900 \ \mathrm{M}_\odot$|⁠, before the emergence of a radiative zone. The suppression of mixing through this insulating layer leads to He enrichment in the core and net mass loss over a time-scale of few Myrs. Eventually, the stars’ mass is reduced to |$\lesssim 30 \ \mathrm{M}_\odot$| as they undergo a MS-to-PostMS transition. This evolutionary pathway is physically distinct from the one found by Cantiello et al. (2021) for ρ_c = 5 × 10⁻¹⁸ g cm⁻³. Transition in our model is caused by the state of mixing in the radiative zone instead of the competition between H-to-He conversion and H replenishment rates in the core. We furthermore find that, for solar-metallicity gas, the area where extra mixing is relevant is limited to ρ_c = 6 × 10⁻¹⁸–1 × 10⁻¹⁷ g cm⁻³ or over a radial extent |$\sim 2 m_8 ^{1/3}$| pc in gap-free AGN discs. This domain is much further away from the SMBH than the region at 10-40 light days where the high α-elements abundance is inferred from the broad-line ratios. For densities lower than 6 × 10⁻¹⁸ g cm⁻³, even with extra mixing, stars evolve into PostMS over ∼10⁸ yr as found by Cantiello et al. (2021).

To understand the origins of this mixing-driven evolutionary divergence between the models of immortal versus evolving stars, we inspect the ratio of the star’s convective to total masses in Fig. 1, where we notice that while for the high density case the star is fully convective when it reaches the equilibrium state (⁠|${\dot{M}}_{\rm Bondi, \Gamma }={\dot{M}}_{\rm wind, \Gamma }$| with λ_⋆ ≃ 1), the lower density case maintains a permanent (outer) radiative zone. Without extra mixing, this radiative buffer prevents the freshly accreted gas from being transported to the core that burns its residual hydrogen before entering the PostMS stage.

Comparing this case against one with identical disc density (ρ_c = 1 × 10⁻¹⁷ g cm⁻³) but including radiative zone extra mixing (also Fig. 1), we see that the star reaches and maintain a high-mass ‘immortal’ state, similar to the higher disc density without extra mixing case. This dichotomy confirms that mixing in the radiative zone is controlling the fate of the star in the low density case. Finally, Fig. 1 also shows that embedded stars are fully convective when they reach an equilibrium state in disc regions with ρ_c ≥ 2 × 10⁻¹⁷ g cm⁻³.

In summary, with the possible absence of gas mixing mechanisms other than convection, stars born in regions of the AGN disc where ρ_c ≤ 1 × 10⁻¹⁷ g cm⁻³ [at |$\gtrsim 3.3 m_8^{1/3}$| pc in unperturbed regions of the disc, equation (3) or in a gap where the disc gas is locally depleted, equation (C2) in Appendix  C] evolve beyond the MS track. Since additional mixing mechanisms may be present in massive, rapidly rotating stars, we need to quantitatively access the likelihood of rapid rotation for the massive stars and determine the efficiency of compositional blending under such circumstances.

3.2 Re-accretion of He-byproducts in the wind

3.2.1 mesa implementation

In this section we investigate the re-accretion of the He-byproducts ejected through the stars’ wind as a possible mechanism for the MS-to-PostMS transition (Appendix  D).

For illustration purpose, we adopt the nominal-ρ_c (= 6 × 10⁻¹⁷ g cm⁻³) model with an initial mass m_⋆ = 1. Once the star’s mass has reached an aymptotic equilibrium value, the model is saved using mesa’s save_model_when_terminate inlist flag. This stellar model is then loaded and resumed using a slightly modified code that incorporates an evolving chemical composition for the accreted material.

Assuming that the accreted material is a mixture of wind-ejected and pristine-disc gas, the H mass fraction of the accreted material (controlled directly as a free parameter in mesa) becomes

where X_⋆ is the H mass fraction in the star (assumed to be the same as that carried by the wind from of the star’s surface), and X_d is that of the AGN disc (assumed to be its solar value). The retention efficiency f_↺ is implemented as a free model parameter in mesa and f_↺ = 0 or 1 for a totally open or closed system, respectively. Analogue prescriptions are also applied for helium and metals. Finally, we let the star to evolve till either to a new equilibrium state is reached or to finish its evolution (we use mesa’s default stopping condition when the core temperature reaches 10^9.5 K).

Here we consider the case where significant fraction of the star’s accreted material originates from locally retained winds (inside the star’s Roche radius, Appendix  C). With f_↺ = 0.95 (equation 7), we restart the simulation from the nominal equilibrium model with |$M_\star = M_0 = 630 \ \mathrm{M}_\odot$| and R_⋆ = 21.5 R_⊙. Although the extra mixing through the radiative zone is switched off (Section 3.1), the star with this boundary condition (ρ_c = 6 × 10⁻¹⁷ g cm⁻³) and initial mass is fully convective and its internal chemical composition is fully homogenized between the outer mass-exchange zone and the nuclear-burning core. In the absence of efficient retention (i.e. with f_↺ < <1), fresh H replenishment would maintain this state indefinitely.

However, efficient retention and recycling of He increases the star’s He mass fraction Y_⋆ and molecular weight μ_⋆ over a time-scale of ∼ 2–3Myr. As the star adjusts to new hydrostatic quasi mass-exchange equilibria (with increasing Y_⋆), it losses mass, at a fraction of |$\dot{M}_{\mathrm{wind, \Gamma }}$| (equation 6 and Appendix  E). Its radius also shrinks while its luminosity is maintained near its Eddington limit with L_⋆ ≃ L_E⋆, (solid blue curves in the left-hand, central, and right-hand panels of Fig. 2). Over ∼3 Myr, the star loses a substantial fraction of its initial mass on its MS track. When the star’s interior eventually runs out of hydrogen, its evolution makes a transition from the MS to the PostMS phase.

3.2.2 Preserve a quasi accretion-wind equilibrium in chemically evolving stars through net mass loss

Evolving MS stars with monotonically declining H fraction loss mass despite having established a quasi accretion-wind equilibrium with |${\dot{M}}_{\rm Bondi, \Gamma } \simeq {\dot{M}}_{\rm wind, \Gamma }$| and λ_⋆ ≃ 1. The rate of net mass loss is a small fraction of |${\dot{M}}_{\rm wind, \Gamma }$|⁠. To better understand the physical origins of this net mass loss along this MS evolutionary track, we plot (right-hand panel of Fig. 3) the evolution of the stellar luminosity L_⋆ during the growth of embedded star due to the accretion of disc gas with a range of Y_d (and thus Y_⋆). These stellar models are constructed with the same code and parameters (with no extra radiative zone mixing) as the growth phase of our nominal model. Contributions from stellar wind is also included with the f_↺ = 0 open box prescription. These growth tracks approaches the Eddington limit (dashed line) where the stars either attain an equilibrium mass and their mass growth is stalled in a ‘immortal state’ for Y_d ≲ 0.5, or start to loss mass and eventually evolve into PostMS for Y_d ≳ 0.5. Note that Jermyn et al. (2022) found a transitionary Y_d ∼ 0.9, and this difference is due to the exclusion of extra mixing in our simulations.

Fig. 3 shows clearly that higher helium abundance in the star (Y_⋆) lead to lower Eddington-luminosity mass, and in consequence lower equilibrium or gain–loss transition mass (M_eq). The time evolution of our nominal simulation with wind-accretion (mass loss curve) is hence simply tracking this Y_⋆–M_eq trend, as also shown on the left-hand panel of the same plot. Physically, this is caused by the dependency of the luminosity-mass (L_⋆–M_⋆) relation on the mean molecular weight μ_⋆ of the star, and thus its helium abundance. For a constant mass M_⋆, stars with increasing μ_⋆ have higher L_⋆. Therefore, stars accrete from disc gas with different Y_d reach the Eddington luminosity L_E⋆ with different masses M_⋆, leading to the Y_⋆–M_eq relation we found in the left-hand panel of the same plot. Analytic approximation for the quantitative M_eq dependence on the molecular weight μ_⋆ is presented in Appendix  E.

We also show (right-hand panel, Fig 2) the λ_⋆ ≡ L_⋆/L_E⋆ controlling whether the star is accreting or losing mass. Overall, throughout the star’s evolution, λ_⋆ is close to unity as the star is simultaneously gaining and losing mass, although the trend is downwards as the star’s luminosity decreases slightly steeper with mass than the Eddington luminosity. It is noticeable that while the luminosity is trans-Eddington, λ_⋆ is actually less than 1 for most of the simulation, even though the star is net losing mass. This result is caused by the tanh smoothing terms in equations (5) and (6) that allow the star to accrete mass for λ_⋆ larger but close to 1, and to lose mass for λ_⋆ smaller but close to 1. For λ_⋆ = 1, |$\dot{M}_{\rm wind,\Gamma }$|/|$\dot{M}_{\rm Bondi,\Gamma }$| = 0.84× |$\dot{M}_{\rm wind}$|/|$\dot{M}_{\rm Bondi}$| instead of 0. On the other hand, for λ_⋆ = 1.1, |$\dot{M}_{\rm wind,\Gamma }$|/|$\dot{M}_{\rm Bondi,\Gamma }$| = 8.8× |$\dot{M}_{\rm wind}$|/|$\dot{M}_{\rm Bondi}$|⁠. This smoothing prescription has no qualitative influence on our results, and is discussed further in Section 5.1.

3.2.3 Compositional-mixing buffer by an outer radiative zone around a nuclear burning core

An unexpected feature of the mass-loss curve is the change in slope around m_⋆ ≃ 485. To understand this change better we plot in Fig. 4 the ratio of the star’s total and convective core masses. The change in mass loss curve slope coincides with the expansion of an initially small outer radiative zone due to the gradual stellar-mass decrease and change in the opacity. The presence of this radiative zone influences the star’s M_⋆–R_⋆ and M_⋆–L_⋆ relations (as functions of Y_⋆ and μ_⋆), changing the mass loss curve slope. Moreover, it introduces a barrier for the freshly accreted material from the outer mass-exchange region to reach the H-burning cores. In the absence of extra mixing through the radiative zone, this insulating layer leads to a close-box with f_↺ = 1 for stellar core (Section 3.1). Thereafter, H in the nuclear furnace at the stellar core is exhausted on a time-scale of Δt_MS ∼ 2.5 Myr (Fig. 4 and Appendix  E) regulated by the chemical evolution (Appendix  E) and the star undergoes a MS-to-PostMS transition when the H fuel in its core is exhausted.

3.2.4 Dependence on the retention efficiency

In Fig. 2 we also show the mass evolution curve for an identical simulation but with f_↺ = 0.90. In this case, while the star initially loses ∼ 100 M_⊙ of mass, it quickly re-establishes a new equilibrium around m_⋆ = m_eq ∼ 500, X_⋆ ∼ 0.51, and Y_⋆ ∼ 0.48. The star remains fully convective such that the fractional amount of fresh supply from the disc gas can fully mix with the nuclear burning core. This modest amount of elemental metabolism enables the star to become ‘immortal’ and to stay on the MS. This numerical result is consistent with the inference from analytic approximation (equation D6 and Appendix  D).

The transitional stellar mass (below which there is a significant radiative zone) may vary slightly depends on μ_⋆, κ, and ρ_c. For discussion purpose below, we adopt the value of the nominal model. It is imperative to note that this critical mass is only comparable to M_eq at tenuous gap-free outer (r_pc ≳ 5 Section 2.2 and Fig. 1) region or M_gap due to gap formation in the inner dense, geometrically thin regions of the disc (below).

3.3 PostMS transition under growth limit due to disc gaps

In the previous subsection, we show that local retention and re-accretion of wind-ejected helium are viable mechanisms to enable stars to undergo MS-to-PostMS transition in the tenuous (with sufficiently low ρ_c) outer (R ≳ a few pc) or in the dense inner disc region provided the retention factors f_↺ ≳ 0.95. Although such a high retention efficiency may be challenging to accomplish, we note that it is only required for the fully convective massive stars (Section 3.2.3). Less massive stars have a radiative zone which can prevent the accreted disc gas from reaching the nuclear burning core and mixing with the He-byproducts. For stars embedded in regions with ρ_c = 6 × 10⁻¹⁷g cm⁻³, radiative zone occurs in stars with |$M_{\rm eq}~\lesssim 485 \ \mathrm{M}_\odot$| (Fig. 4). In regions with lower ρ_c, the onset of radiative zone occurs in stars with larger M_eq (Fig. 1).

A natural physical process to stall accretion is the formation of gaps in the disc around the embedded stars (Appendix  C). In the inner (r_pc < 1) dense, thin-disc region, gap-formation sets an upper growth limit to M_gap ≲ a few |$10^2 \ \mathrm{M}_\odot$| [equation (C2) with λ_⋆ < <1]. Note that the gap-formation and growth-stalling conditions are determined by the discs’ aspect ratio rather than by the unperturbed ρ_c (Appendix  C) such that they can set an upper limit for M_⋆ over a wide radial range in the disc’s sub-pc inner region. Nevertheless, gap formation decreases the ambient disc density which is also favourable for the preservation of the embedded stars’ radiative zone (cf Figs 1 and 4). We also note that as shown in equation (C2), M_gap scales as h³. As a nominal case we assumed h ≃ 10⁻² as discussed in Section 2.1 based on the reverberation mapping of the metal-rich BLR in NGC 5548. In the outer regions of AGN discs (beyond a few pc), h may be a few times higher, leading to the possibility of very high M_gap’s. However, regions with such high h’s and marginal gravitational instabilities are likely to fragment and to break up into systems of clumpy clouds (similar to AGNs’ dusty tori). In such porous medium, the stellar accretion rate is substantially smaller than |${\dot{M}}_{\rm Bondi}$| (equation 4). (Moreover, in these regions, ρ_c is sufficiently small for radiative zone to be preserved in very massive stars; see Section 3.1 and Fig. 1.) Based on the results of Section 3.2.3, we now consider the possibility that the accretion is quenched with |$M_\star \lesssim M_1 \simeq 485 \ \mathrm{M}_\odot$| so that the presence of radiative zone would lead mixing insulation and MS-to-PostMS transition regardless the magnitude of f_↺ for the outer mass-exchange region.

For illustration, we present in Fig. 5, three models with the same parameters (ρ_c = 6 × 10⁻¹⁷ g cm⁻³, sound speed to c_s = 10⁶ cm s⁻¹, f_↺ = 0) and accretion and wind rates (equations 5 and 6) as the nominal model. We introduce a cutoff in the accretion rate for stars with m_⋆ ≥ 300. For the first and second models, we start with the nominal (with m_⋆ = 630 and enhanced Y_⋆) and a lower-initial-mass (m_⋆ = 300 and solar composition) model respectively. We do not include extra mixing in the radiative layer. In both models, stars evolve onto PostMS phase with similar evolutionary tracks, transitional and pre-core-collapse masses as those of the nominal model with f_↺ = 0.95.

We simulated another mesa model with an initial mass m_⋆ = 300. In this model, we include the extra mixing prescription for the radiative zone following the prescription by Cantiello et al. (2021). In this case, the star evolves into a ‘immortal’ equilibrium state with m_⋆ ≃ 300, and X_⋆ ≃ 0.5. Despite the enhancement in Y_⋆, μ_⋆, L_⋆ < <L_E⋆ (i.e. λ_⋆ < <1). This result is consistent with our analytic approximation (Appendices  D and  E). Transition from MS to PostMS occurs provided the freshly accreted disc gas does not replenish and refresh the nuclear burning core. These results clearly indicate the importance of the radiative buffer zone which prevents the replenishment of fresh supply of disc gas from reaching the nuclear burning core.

3.4 Stellar heating rate, population, and mass function

In Section 1.1, we indicated that the main motivation for the SEPAD model is to provide adequate auxiliary power with a stellar surface number density s_⋆ ≃ Q⁻/L_⋆. For ‘immortal’ stars, L_⋆ ≃ m_eq L_E⊙. But since the luminosity of evolving stars decreases with time, s_⋆ is determined by its time averaged value (from the mesa model) over their life span τ_⋆ ≃ 5.1 Myr,

which is equivalent to the Eddington luminosity of a time-average |${\bar{M}}_{\rm eq}~\simeq 143 \ \mathrm{M}_\odot \simeq 0.23 M_0$| stars. In order to provide adequate auxiliary power to main marginal self gravity for the disc, the magnitude of |$s_\star ({\bar{M}}_{\rm eq}~)= Q^-/ {\bar{L}}_\star$| for this evolving population needs to be ∼4.4 times larger than that, s_⋆(M₀) = Q⁻/L_⋆(M₀), of the ‘immortal’ stars with |$M_0=630 \ \mathrm{M}_\odot$|⁠, X₀ = 0.69, Y₀ = 0.29, and Z₀ = 0.02. Since the α and Fe released per individual embedded PostMS stars is approximately the same and the total amount of He yield is determined by the required auxiliary power to maintain a state of marginal gravitational stability, this enlargement of the stellar population (s_⋆) enhances the disc-pollution rates of α and Fe relative to that of He (Section 4.3).

Within a disc radius R_d (where constant-|${\dot{M}}_{\rm d}$| flow can be approximately maintained), the total population of stars is |$N_\star \simeq 2 \pi \int _0 ^{R_{\rm d}} s_\star R \mathrm{ d}r$| in this model. In a steady state with a constant N_⋆ or s_⋆, the formation rate |${\dot{N}}_\star \simeq N_\star /\tau _\star$| is a constant. The stellar mass function

varies with m_⋆ (rather than a single power law) which can be evaluated numerically [from the first part of equation (9) and the left-hand panel of Fig. 2] to be |$\sim {m_\star / \vert {\dot{m}}_\star } \vert \tau _\star \sim 0.2$| for the mass ranges m₁ ≥ m_⋆ ≥ m₂ and m₂ ≥ m_⋆ ≥ m₃ on the MS and PostPM tracks, respectively (where m₁, m₂, and m₃ are defined in Section 4.2). The mass function can also be analytically approximated with a similar value [from the second equation in equation (9) where the first and second terms on the denominator corresponds to MS and PostMS stars and A_y, A_z, |${\dot{Y}}_\star (= -{\dot{X}}_\star)$|⁠, and |${\dot{Z}}_\star$| are estimated in Appendices  D and  E]. This top-heavy mass function is similar to the IMF inferred by Miller & Scalo (1979).

Based on this mass function, we can compute the average luminosity

for a population of quasi-equilibrium stars with N₀ as its normalization factor. This total luminosity is equivalent to the Eddington luminosity of a time-average |${\bar{M}}_{\rm eq}~\simeq 143 \ \mathrm{M}_\odot \simeq 0.23 M_0$| stars (equation 8) for the fiducial stellar model. The equivalent mass |${\bar{M}}_{\rm eq}$| would be modified by M_gap if the stellar mass growth is quenched by gap formation (a smaller |${\bar{M}}_{\rm eq}$| would lead to a larger N₀ and s_⋆). This mass function also has implications on the merger probability of coexisting stars (with |$s_\star ({\bar{M}}_{\rm eq}) \sim 4.4 Q^-/m_0 L_{\rm E \odot }$| and total population N_⋆), star formation rate per unit area (⁠|${\dot{s}}_\star = s_\star / \tau _\star$|⁠), and over the entire disc |${\dot{N}}_\star$| which will be considered in future investigations.

4 EFFECTS ON THE COMPOSITION OF THE DISC

4.1 mesa model of PostMS evolution

In this section we quantify the star’s PostMS evolution and its effects on the chemical composition of the AGN disc. We hence investigate the star’s chemical yield and the composition of the remnant. We focus on the local wind retention scenario, as we find it to be the more realistic mechanism to take the ‘immortal’ stars off the MS. For our nominal model with f_↺ = 0.95, the MS-PostMS transition occurs and He burning starts after 4.7 Myr. This corresponds to the first sharp peak in the λ_⋆ ≡ L_⋆/L_E⋆ curves of Fig. 2. The second peak, at 5.1 Myr, corresponds to the onset of Ne burning and the end of the simulation as the star’s centre reach a temperature of 10^9.5 K. Above this temperature, our 21 elements nuclear network becomes insufficient. At this point the star is ∼ 13 M_⊙ and it is expected that silicon burning will follow till exhaustion, leading to a possible core collapse and SN II explosion.

4.2 Disc’s He pollution rate by evolving and immortal stars

In the limit of high retention efficiency (f_↺ ≳ 0.95), negligible He, α, and Fe-byproducts are released to the disc and reaccreted onto the stars through stellar metabolism. Nevertheless, there would be a net loss of He from the stars to the disc if M_eq declines due to the accumulation of He and the increases in both Y_⋆ and μ_⋆.

In the fiducial model with f_↺ = 0.95, stars’ evolution can be characterized by three stages (Appendix  E). Stage 1: the radiative region of high masses from M_⋆ ∼ M₀ (with Y₀ = 0.29) to |$M_1 = 485 \ \mathrm{M}_\odot$| (with Y_⋆ = Y₁ ≃ 0.5) and their M_⋆–Y_⋆ dependence can be approximated by A_μ ≡ dlnM_eq /∂lnμ_⋆ ≃ −2 (equation D3). Stage 2: MS stars loss mass between |$M_1 \le M_\star \le M_2 \simeq 28 \ \mathrm{M}_\odot$| (with Y₂ = 1 and A_μ ∼ −4.7) and an expanding outer radiative zone which insulate stars’ nuclear burning core from their stellar surface and impose f_↺ = 1. When the residual H in the core is exhausted, they evolve onto the PostMS track. Stage 3: PostMS stars continue to reduce their mass from M₂ to |$M_3 \simeq 13 \ \mathrm{M}_\odot$| as they convert He into mainly α elements with an efficiency ϵ_α ≃ 10⁻³. Once again, the increase in μ leads to a decrease in the accretion-wind equilibrium mass. Soon after C burning is initiated. the remaining cores collapse and undergo SN II with α, Fe ejecta and rBHs.

As He is produced during the MS (stages 1 and 2), their release rate into the disc is |${\dot{M}}_{\rm He}=- Y_\star {\dot{M}}_{\rm eq}$| (equation D3) where |${\dot{M}}_{\rm eq}~= A_{\rm y} M_{\rm eq}~{\dot{Y}}_\star / Y_\star$| is the changing rate of the equilibrium mass |$m_{\rm eq}~(Y_\star) \propto Y_\star ^{A_{\rm y}}$| (equation E5). Neglecting variations in A_y, the total He return to the disc is

which is comparable to the He yield from the mesa model (⁠|$355 \ \mathrm{M}_\odot$|⁠, left-hand panel, Fig. 6).

During MS stages 1 and 2, |$\Delta M_{\rm H} = 275.4 \ \mathrm{M}_\odot$| of hydrogen was not burned and was returned to the disc. This amount is slightly larger than the difference between the star’s initial mass (⁠|$\sim 630 \ \mathrm{M}_\odot$|⁠, Section 3.1) and ΔM_He + M₂(≃ 383). This small discrepancy can be attributed to the much reduced but still finite accretion of disc gas during the MS evolution. Subtracting the initial He mass (with initial Y₀ = 0.29, Z_⋆ = Z₀ = 0.02, and |$M_\star =630 \ \mathrm{M}_\odot$|⁠), the total change ΔY ≃ ΔM_He(1 − Z_⋆)/(ΔM_H + ΔM_He) − Y₀ ≃ 0.26.

Averaged over the stellar lifespan τ_⋆ ∼ 4–5 Myr (Fig. 2 and Appendix  E), the He pollution rate in the disc is |${\dot{M}}_{\rm He} \sim ({\Delta M}_{\rm He}- M_0 \times Y_0) /\tau _\star \sim 4 \times 10^{-5} \ \mathrm{M}_\odot \, {\rm yr}^{-1}$| per star. In constrast, the maintenance of a ‘immortal’ equilibrium state (with λ_⋆ ∼ 1, f_↺ = 0, Section 3.1, and |$M_\star \simeq M_0 =630 \ \mathrm{M}_\odot$|⁠, Fig. 1) requires a H-to-He conversion rate

This difference is due to |${\bar{M}}_{\rm eq}~\simeq 0.22 M_0$| so that |${\bar{L}}_\star \simeq {L_{\rm E}}_\star ({\bar{M}}_{\rm eq}) \simeq 0.22 {L_{\rm E}}_\star (M_0) \simeq 0.22 L_\star (M_0)$| during stage 2 of MS evolution (equation 8 and Section 3.4). In order for the evolving stars to provide an adequate auxiliary power to maintain a state of marginal stability (Q ∼ 1 with |$Q^+ _\star \simeq Q^-$|⁠), stars which undergo the MS-to-PostMS transition need to be more numerous than the ‘immortal’ stars (i.e. with a constant |$s_\star ({\bar{M}}_{\rm eq}~) {\bar{L}}_\star$| determined by Q⁻ in Section 3.4). This compensation implies that the difference in the total He pollution rate for the disc by the evolving stars may be indistinguishable from that by the ‘immortal’ stars. Although it is not possible to eliminate the possibility of ‘immortal’ stars, solely based on the modest observed value of Y_d, it does eliminate predominant He pollution as a potential mechanism to enable the MS-to PostMS transition (Jermyn et al. 2022).

4.3 Pollution rate of α elements

MS stars’ luminosity is powered by the CNO process which does not increase the abundance of C+N + O from that of the disc gas (⁠|$\sim 2 {{\ \rm per \, cent}}$|⁠). Nevertheless, it leads to the conversion of C + O to N as secondary elements during the initial mass rump up to M₀ and subsequent MS stages 1 and 2 mass loss to M_a and M₂, respectively (left-hand and middle panels of Fig. 7) with substantial increase in N/(C + O) to super solar values (≳ 20). Along with He, 5.1 M_⊙ of N and negligible amounts of C and O are returned to the disc during the MS stages 1 and 2 of evolving stars (left-hand panel of Fig. 6 and Section 4.4). The N/(C+O) ratio may increase further for the immortal stars if they re-accrete gas already laden with secondary byproducts (i.e. super solar N/(C + O)).

The main difference between the ‘immortal’ and the evolving stars is in the production of α elements and Fe during the PostMS evolution of the latter. As they make a transition between stages 2 and 3, the evolving stars are primarily composed of He. During PostMS stage 3, He is converted into α at a rate |${\dot{M}}_\alpha \simeq L_\star / \epsilon _{\rm He} c^2 = \lambda _\star M_\star / \epsilon _\alpha \tau _{\rm Sal}$|⁠. On the time-scale |$\Delta t_{\rm PostMS} \sim M_\star /{\dot{M}}_\alpha \simeq \epsilon _\alpha \tau _{\rm Sal} \sim 0.5$| Myr, all the He contents would be exhausted (Appendix  E) in agreement with mesa model (Fig. 2). Although this time-scale is an order of magnitude longer than the super-Eddington |$M_\star /{\dot{M}}_{\rm wind}$| (with λ ≳ 1 and A_λ ∼ 1 in equation B1), the duration of stage 3 is prolonged by the quasi accretion-wind equilibrium (with λ_⋆ ≲ 1; Fig. 2). Despite the re-accretion of the disc gas, the presence of an outer radiative region (Fig. 4) prevent the replenishment of fresh fuel to the He-burning core. The net reduction of the PostMS stars M_eq is due to the ramp up of α, Fe, and μ_⋆, analogous to the quasi-equilibrium mass loss during the MS evolution (Section 3.2.2).

As calculated in Appendix  E, during PostMS stage 3, A_μ ≃ −1.5, X_⋆ = 0, μ_⋆ ≃ 1.3 → 2, A_z ≡ ∂lnM_eq/∂lnZ_⋆ = A_μ∂lnμ_⋆/∂lnZ_⋆ ≃ A_μZ_⋆/(3 − Z_⋆) ≃ −0.01 → −0.75, the equilibrium mass,

as Y_⋆ ≃ 1 → 0 and Z_⋆ ≃ 0.02 → 1. The rate of α and Fe elements return to the disc is |${\dot{M}}_{\rm Z}=- Z_\star {\dot{M}}_{\rm eq}$| where |${\dot{M}}_{\rm eq}~= A_{\rm z} M_{\rm eq}~{\dot{Z}}_\star / Z_\star$| is the changing rate of the equilibrium mass M_eq (Z_⋆) (equation 13). Neglecting variations in A_μ, the total Z_⋆ return to the disc is

which is comparable to the heavy element return from the mesa model (∼ 4.3 M_⊙ of C and 3 M_⊙ of O, along with small amounts of Ne and virtually no N; left-hand panel of Fig. 6). A fraction of the stellar mass |$M_2-\Delta M_{\rm Z}-M_3 \simeq 9 \ \mathrm{M}_\odot$| is released as He. Most of the α (⁠|$9.6 \ \mathrm{M}_\odot$| O, |$0.6 \ \mathrm{M}_\odot$| Ne, |$1.6 \ \mathrm{M}_\odot$| Mg, |$0.7\ \mathrm{M}_\odot$| Si) and Fe are contained in the pre-collapsing core (with a total mass |$M_3\simeq 13 \ \mathrm{M}_\odot$|⁠) at the end of stage 3 (right-hand panel of Fig.  6). A fraction of M₃ is released as α and Fe ejecta (with sub-solar [Fe/α]) from the subsequent SN II and the remainder is contained in the collapsed rBHs (Sukhbold et al. 2016). Note that the mass of the rBH is assumed to be the mass of the star at the end of the mesa simulation, when core collapse is starting as seen in Fig. 2 (central panel). This mass is of course not necessarily the final mass of the remnant, so it should be taken as an upper limit. The total mass of extra α return (⁠|$\sim \Delta M_{\rm Z}+M_3 \simeq 20 \ \mathrm{M}_\odot$|⁠, Fig. 6) is |$\sim 3{{\ \rm per \, cent}}$| that of H + He, i.e. ΔZ ∼ 1.5Z₀, leading to the α enrichment of the disc gas.

The ratio ΔY/ΔZ ∼ 9 which is substantially above that (∼1.8) inferred from galactic chemical evolution (Carigi & Peimbert 2008; Peimbert, Peimbert & Luridiana 2016; Valerdi, Peimbert & Peimbert 2021). However, the magnitude ΔM_He (and therefore ΔY) of individual evolving star would be reduced if the initial equilibrium mass (M₀) was limited by gap formation with a mass M_gap < M_eq (X₀, Y₀, Z₀) (Appendix  C). Moreover, the total α and Fe production rate is proportional to |${\dot{N}}_\star \Delta M_{\rm Z} \simeq N_\star \Delta M_{\rm Z} / \tau _\star$| and the population (s_⋆ and N_⋆) of coexisting evolving stars have much larger (by a factor |$\sim M_0/{\bar{M}}_{\rm eq}~\sim 4.4$| in the fiducial model) than those of the ‘immortal stars’ population (Sections 3.4 and 4.2). But the H-to-He conversion rate is regulated by the auxiliary power criterion which does not depend on the stars’ M_⋆ and N_⋆ (Section 4.2). Consequently ΔY/ΔZ is much reduced from that inferred from the nominal model, and more in line with the observed values (Huang et al. 2023).

4.4 Multiple generations of evolving stars

The pre-collapse cores of evolving stars contain the byproducts of triple-α reaction during their PostMS stage 3 which are mostly O with negligible amount of N. However, the byproducts of CNO burning during the MS evolution of subsequent generation of stars efficiently converts O + C into N as secondary products. This conversion is evident in the ramp up to m_eq = 630. Although the disc gas is assumed to have a solar composition (with N/(C+O) ≲ 0.08), CNO burning substantially ramp up the stellar N abundance (with N/(C + O) ∼20 in M₀, Fig. 6). During MS stages 1 and 2, the N return is an order of magnitude larger than the O+C return to the disc. With additional α byproducts during the PostMS stage 3, the total N/(C + O) in the returned gas is ∼ 0.7, which is an order of magnitude larger than its solar value. This secondary production process is in good agreement with that observed (Edmunds & Pagel 1978; Barbuy 1983; Tomkin & Lambert 1984; Thuan, Izotov & Lipovetsky 1995; Henry, Edmunds & Köppen 2000) and expected from the galactic chemical evolution (Talbot & Arnett 1974; Arnett 1996; Meynet & Maeder 2002).

5 SUMMARY AND DISCUSSIONS

In the high-density environment of AGN discs, embedded stars readily accrete gas and gain mass. The growth is halted either by the formation of gap through tidal truncation or the onset of mass loss through intense wind when their luminosity due to nuclear burning approaches the Eddington limit. The latter stellar metabolism generally leads to a high-equilibrium mass. If the re-accreted gas is well mixed with the nuclear burning core, these ‘immortal’ stars remain on the MS and do not contribute to the super solar α and Fe abundance commonly observed in AGNs.

In an attempt to resolve this discrepancy between theory and observation, we investigate the effect of re-accretion of material which is pre-ejected by the star. The main outstanding issues are whether embedded stars can be preserved as ‘immortal stars’ and sufficient conditions for them to undergo transition from MS to PostMS evolution.

We carry out numerical simulations with the mesa code. Adopting the same prescriptions on accretion rate, wind flux, boundary conditions, and internal structure as (Cantiello et al. 2021), we reproduce previous models for ‘immortal stars’ in the dense inner region and for evolving stars in the outer tenuous regions (beyond a few AU) of AGN discs. We then modify our mesa model with: (1) the adaptation of the solar composition rather than prestine composition for the AGN discs, (2) the suppression of mixing in the radiative region, (3) the introduction of a retention factor, and (4) the consideration of a upper-limit for the stellar mass due to gap formation. We found effects (2), (3), and (4) can lead to the transition from ‘immortal’ to evolving stars.

High retention efficiency recyles most He byproducts on the MS and it is equivalent to a ‘closed-box’ model. After creating an ‘immortal’ star with mesa and then turning on 95 per cent retention efficiency, we find that this set of initial and boundary conditions causes the stars to lose mass on a chemical evolution time-scale of few Myr. This net mass loss is due to the gradual He accumulation and molecular weight increase which change the L_⋆–M_⋆ relationship and the value of M_eq when the Eddington limit is approached. The star runs out of H as its M_eq is reduced to |$\sim 30\ \mathrm{M}_\odot$| and it transitions to PostMS evolution with the onset of triple-α reaction. Subsequently, nuclear synthesis generates α elements through α chain reaction. This model shows that with high retention efficiency ‘immortal’ stars are impermanent and have a short MS lifetime of few Myr.

In relatively thin AGN discs (Starkey et al. 2022), the embedded stars’ tidal torque leads to gap formation and their mass reaches an upper-growth-limit. Provided this mass is less than |$\lesssim 600-700 \ \mathrm{M}_\odot$|⁠, the evolving stars have a radiative zones in their envelopes which can prevent compositional mixing between the nuclear burning cores and the surface region where gas is being exchanged between the stars and the disc. This insulation also enables the accumulation of He ashes, quenches H replenishment, and lead to the MS-to-PostMS transition on the time-scale of a few Myr.

During the PostMS stage, the evolving stars continue to loss mass as their μ_⋆ increases further until silicon is ignited when their mass is reduced to |$\sim 13-15 \ \mathrm{M}_\odot$|⁠. Due to technical limitations, we terminate our calculation. By the end of their PostMS track, individual evolving stars would have ejected hundreds of M_⊙ of He, several M_⊙ of α elements, and a fraction of |$\ \mathrm{M}_\odot$| of Fe back into the AGN discs, significantly affecting their chemical composition. When the polluted disc gas is being recycled, reprocessed, and returned by the next generation stars, their N/(C + O) ratio is reset by the CNO reaction on the MS track. These byproducts are mixed with the disc gas and together they diffuse toward the SMBH at the centre of AGN discs. In combination, these processes determine the abundance distribution of AGN discs.

5.1 Uncertainties of the numerical model

In our numerical model multiple simplifications were used to enhance numerical stability. Here we discuss the effects of relaxing these simplifications. In Fig. 2 we show curves (in orange) where we removed the tanh smoothing terms from the accretion and mass loss equations, and replace them with a manual condition where |$\dot{M}_{\rm Bondi}$| is 0 for L_⋆ > L_E⋆, and |$\dot{M}_{\rm wind}$| is 0 for L_⋆ < L_E⋆. We find that this approach, coupled with a manually set maximum limit on the time-step, is enough to keep the numerical calculations stable. The qualitative behaviour of the star is unchanged, where it loses mass due to helium buildup till it reaches PostMS. Comparing the mass loss curves for the nominal and ‘numerically unsmoothed’ models, we find that they reach the same final mass around 13 M_⊙, although the ‘numerically unsmoothed’ model reaches it ∼ 1.5 Myr earlier. The Eddington factor for the unsmoothed model oscillates more closely around unity as can be seen in the middle panel of Fig.  2, since the star is not allowed to lose mass for L_⋆/L_E⋆ < 1, in contrast with our nominal model. In this case, L_⋆/L_E⋆ keeps oscillating between two states where it gains mass during a step increasing its L_⋆/L_E⋆ to unity, then loses mass during the next step decreasing its L_⋆/L_E⋆ to below unity. We emphasize that, in reality, no star can exceed the surface Eddington limit. Stars at the Eddington luminosity are not in hydrostatic equilibrium. The reasons our models do briefly exceed L_E⋆, are almost entirely numerical as explained above and not due to the physical assumptions of the model or the treatment of L_E⋆. Fundamentally, any similar 1D numerical model will suffer from transient Super-Eddington luminosities unless the time-step used is extremely small (hindering long term model evolution). In a hypothetical global fully hydrodynamical 3D simulation, the situation is partly rectified by the 3D flow where the star, at each time-step, can simultaneously accrete mass at the cool equator, and loose mass at the hotter poles. However, even this model is still dependent on time-stepping that might lead to transient Super-Eddington luminosities.

5.2 Effects of suppressed winds and the μ_⋆ gradient

Here we start by investigating the role the wind efficiency factor A_λ in equations (6) and (B1). So far we assumed A_λ ∼ 1, implying maximally efficient winds. In reality, effects such as the ‘porosity’ of the star’s atmosphere and photon tiring can suppress the winds efficiency (Owocki, Gayley & Shaviv 2004). Manually decreasing A_λ by a factor 10 in our nominal simulation, however, does not lead to qualitative changes. We find that winds efficiency need to be suppressed by a factor 30–100 for the star to evolve PostMS. A suppression of this magnitude is unlikely, and thus reduced winds efficiency probably does not play a central role in evolving the stars to the PostMS phase. However, it is possible that moderately suppressed winds, coupled with some other effects such as local wind-retention discussed above, might be enough to force the star to evolve through the MS to PostMS phase. Finally, we also investigated the effects of using the Ledoux convection criteria instead of the default Schwarzschild criteria and allowing for semiconvection, and found that this does not affect our nominal model either.

5.3 Retention efficiency, mixing barrier, mass growth limit

The analyses and simulations presented here highlight two sufficient criteria for MS-to-PostMS transition: (i) a high retention efficiency or (ii) a mixing barrier. We postulate the possibility that in the limit of λ_⋆ ∼ 1, there is a region between the Bondi and Roche radius where the momentum and mass fluxes of the outgoing-wind and incoming-accretion flux are well mixed and subsonic. High retention efficiency is based on the assumed stability of this two-flow pattern. This hypothesis needs follow-up verification and examination with some 3D radiative hydrodynamic simulations.

Under our assumption that there is no compositional mixing in the radiative zone. This structure provides a buffer between the mass-exchanging outer layers and the nuclear-burning core of the stars. In contrast, Cantiello et al. (2021) suggest the possibility of extra mixing due to rotational instability. Rapid spin (a necessary condition for rotational instability and hence compositional mixing throughout the radiative zone) is commonly observed in stand-alone massive stars. Whether similar conditions may be common for the massive stars which grow through the accretion of turbulent disc gas is another area of natural follow-up studies with 3D hydrodynamic simulations.

We introduce a scenario that gas accretion onto embedded stars is quenched with an upper limit on their masses M_gap. For illustration purposes, we choose |$M_{\rm gap}=300 \ \mathrm{M}_\odot$| (Fig. 5). Provided this mass is sufficiently modest to preserve a significant radiative zone, we can bypass the stringent retention criterion and ensure MS-to-PostMS transition for stars embedded over a large radial range in AGN discs. The gap formation criterion in Appendix  C is derived and simulated for isolated low-mass companions (such as planets) embedded in their natal discs on nearly circular orbits (Lin & Papaloizou 1993; Kley & Nelson 2012). Tidal interaction between embedded stars and AGN discs endures perturbation from neighboring coexisting stars. The effectiveness of mass-growth limit needs to be further studies and simulated.

ACKNOWLEDGEMENTS

We thank Andrew Cumming for suggestions, guidance, and major contributions throughout this investigation. We thank Jiamu Huang, Greg Shields, Stan Woosley, Yixian Chen, and Jamie Law-Smith, and Chris Fryer for useful conversation. We thank the anonymous referee for their insightful comments that significantly improved the manuscript. This work is supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CAP³.

DATA AVAILABILITY

The data underlying this article (numerical simulations output files) will be shared on reasonable request to the corresponding author.

References

APPENDIX A: BONDI ACCRETION RATE NEAR THE EDDINGTON LIMIT

We first consider the conventional analysis for Bondi accretion in the limit that the stars’ luminosity is near their Eddington limit. The influence of radiation pressure from stellar luminosity on the accretion rate and wind flux are separately approximated by |${\dot{M}}_{\rm Bondi, \Gamma }$| (equation 5) and |${\dot{M}}_{\rm wind, \Gamma }$| (equation 6), in the absence of each other. Under the assumption that these two opposing flows are fully mixed and combined into a single spherically symmetric radial (r) flow, its radial velocity U_⋆ relative to the star (positive/negative values correspond to outflow/inflow, respectively) is determined by the momentum equation

where the radiation pressure term is incorporated by the λ_⋆GM_⋆/r² term. The magnitude of λ_⋆ ≡ L_⋆/L_E⋆ is an increasing function of m_⋆. In a steady state (when the mass flux |${\dot{M}}_{\rm Bondi} =4 \pi \rho U_\star r^2$| is independent of r; equation 4), the sonic point (where U_⋆ = c_s) occurs at the Bondi radius

Stellar wind is launched from MS stars’ R_⋆ ≪ R_B when λ_⋆ < <1. But, with λ_⋆ → 1, R_B → R_⋆, the entire accretion flow becomes subsonic (with U_⋆ ≲ c_s), and the effective accretion rate is greatly quenched

where c_s and ρ_c are obtained from equations (2) and (3). This derivation provides the basis for the prescription for equation (5).

APPENDIX B: ACCRETION-WIND EQUILIBRIUM

A simplifies prescription for the wind

where |$R_\star \simeq m_\star ^{0.6} \mathrm{ R}_\odot$| and R_⊙ are the radius of star and the Sun, |$r_\bullet = G \ \mathrm{M}_\odot /c^2$| and |$\ \mathrm{M}_\odot$| are the gravitational radius and mass of the Sun, and |$\tau _{\rm Sal} (= {L_{\rm E}}_\odot / \ \mathrm{M}_\odot c^2= 4.5 \times 10^8$| yr) is the Salpeter time-scale. In comparison with equation (6), the factor A_λ = [1 + tanh (10λ_⋆ − 10)]/2 which = 0, 0.5, and 1 for λ_⋆ < <1, 1, and >>1, respectively.

An accretion-wind equilibrium (⁠|${\dot{M}}_{\rm Bondi} = {\dot{M}}_{\rm wind}$|⁠) is attained with

For m_⋆ ∼ 600–1000 at r_pc ∼ 1, λ_⋆ ∼ 0.9–1 and these results are consistent with the asymptotic equilibrium mass M_eq and Eddington factors computed with the mesa models without recycling f_↺ = 0 (Fig. 1). In this equilibrium (with M_⋆ = M_eq), the Bondi accretion time-scale |$\tau _{\rm Bondi} \equiv {M_\star / {\dot{M}}_{\rm Bondi}}$| (equations A3 and B1) becomes

APPENDIX C: ROCHE-LOBE FLOW MODIFIED BY RADIATION PRESSURE

Analogous to the emergence of proto-gas-giant planets, the embedded stars also tidally interact with the disc near their Roche/Hills radius

This equation takes into account the reduction of the star’s gravity by the radiation pressure imposed by its luminosity. At relative small radii, it is possible for embedded stars’ M_⋆ towards an M_eq  ≳ 3^1/3h³M_•/(1 − λ_⋆) with λ_⋆ < 1. But their growth would be quenched (D’Angelo, Kley & Henning 2003; Dobbs-Dixon, Li & Lin 2007; Bodenheimer et al. 2013; Tanigawa & Tanaka 2016; Rosenthal et al. 2020; Li et al. 2021) by the formation of a gap in the proximity of their orbit (Lin & Papaloizou 1993) when their mass M_⋆ → M_gap ≃ h³M_•. Substitute h (equation 2) gaps form when their normalized m_⋆ exceeds

This gap-opening process also reduces the local Σ_g, ρ_c (Duffell & MacFadyen 2013; Kanagawa et al. 2015; Chen et al. 2020) and limits stars’ asymptotic R_B ∼ R_H ∼ H and impede their accretion with M_⋆ ∼ M_gap < M_eq and λ_⋆ <<1.

At relatively large R, embedded stars continue to gain mass until m_⋆ → m_eq with λ_⋆ ≃ 1. In this limit, m_gap ≳ m_eq, gap formation does not occur, and in the absence of radiative feedback, |$G M_{\rm eq}/c_s^2 \lt (M_{\rm eq}/3 M_\bullet)^{1/3} R \lt H$|⁠. Although R_H also retreats as the stars’s luminosity reaches its Eddington limit (i.e. λ_⋆ → 1), it does not decrease as steeply as R_B (equations A2 and C1). When |${\dot{M}}_{\rm Bondi}$| is limited by |${\dot{M}}_{\rm wind}$| with λ_⋆ ∼ 1 and R_B reduced to ∼R_⋆ while R_H remains >R_B. We adopt the approximation f_↺ ≃ 1 (Section 3.2) under the assumption that these embedded stars would recycle the He polluted gas in the surrounding nearly static pressure-supported envelope between R_B and R_H that is tidally bound to them. The fully mixed wind-accretion assumption needs to be verified by future investigations.

APPENDIX D: CRITICAL RETENTION EFFICIENCY

An accretion-wind equilibrium is maintained with an influx of fresh disc gas at a rate |$(1-f_\circlearrowleft) {\dot{M}}_{\rm Bondi}$|⁠. The rate of change in the total hydrogen mass due to accretion and wind is

To power L_⋆ on the MS, H is also converted to He with an efficiency ϵ_He(≃ 0.007). The net H mass change rate is

In addition, changes in composition (X_⋆ and μ_⋆) modify the equilibrium mass (Sections 3.2.2 and  E)

is a positive factor of order unity (Fig. 3 and Appendix  E). The net compositional change

From equations (D1), (D3), and (D4), it follows

Inadequate replenishment of H (first term inside the bracket) to compensate the conversion of H to He (second term inside the bracket) would lead to a decline in X_⋆ (i.e. |${\dot{X}}_\star \lt 0$|⁠) and stars’ equilibrium mass (i.e. |${\dot{M}}_\star \lt 0$|⁠).

With a substitution for λ_⋆τ_Bondi/τ_Sal (equation B3), |${\dot{X}}_\star \lt 0$| (X_⋆ declines well below X_d = 0.7 over time) provided

With A_λ ≃ 0.5 (for λ_⋆ ≃ 1) and m_⋆ ∼ 600 (from mesa models in Section 3.1), we find the condition for total H depletion is f_↺ ≃ 0.99.

If we set f_↺ = 0.95 as shown Section 3.2, a minimum X_⋆ ∼ 0.5 would be reachable before |${\dot{X}}_\star$| vanishes (equation D5). Since |${\dot{Y}}_\star = - {\dot{X}}_\star$|⁠, the magnitude of μ_⋆ increases and the equilibrium mass M_eq decreases ≲ 485 with the He enhancement (equation B2). Below this mass, further decline in the X_⋆ leads to a rapidly expanding radiative region in the outer envelope of the star (Figs 4 and 3.2.2). This zone segregates the replenished gas and the stellar interior which is equivalent to resetting insulation with f_↺ = 1 for the nuclear burning core. This internal structure modification enables the star to undergo MS-to-PostMS transition. But with f_↺ ≲ 0.9, the X_⋆ decline is quenched with X_⋆ only slightly less than X_d while m_⋆ ≳ 485 with negligible radiative barrier to prevent the mixing of the accreted material and the burning core. In this case, stars may preserve their longevity.

APPENDIX E: EQUILIBRIUM-MASS AND MOLECULAR-WEIGHT EVOLUTION

In Section 3.2.2, we attribute the evolution of equilibrium mass obtained from the mesa model to increases in the molecular weight μ. For an analytic approximation, we start with the definition of the radiative pressure gradient,

which taken with dP/dr = −GMρ_⋆/r² gives

With P = P_rad + P_gas and β = P_gas/P = 1 − P_rad/P,

where κ_γ is the opacity inside the star. For an ideal gas P_gas = R_gρ_⋆T_⋆/μ_⋆ where R_g is the gas constant and the molecular weight for fully ionized gas μ_⋆ ≃ 4/(3 + 5X_⋆ − Z_⋆) = 4/(8 − 5Y_⋆ − 6Z_⋆) reduces to

during the initial MS evolution when Z_⋆ < <1. During the PostMS with X_⋆ < <1,

which is negative |$\sim {\mathcal {O}} (-1)$|⁠.

With radiation pressure (⁠|$P_{\rm rad}= a T_\star ^4/3$|⁠) and radiation transfer (⁠|$F_{\rm rad} =(ac / 3 \kappa _\gamma \rho _\star) d T_\star ^4/ \mathrm{ d}r$|⁠), Owocki & Shaviv (2012) adopt the Eddington approximation with scaling relations |$\rho _\star (r) \sim 3 M_\star /4 \pi {\rm R}_\star ^3$|⁠, T_⋆ ∼ (1 − λ_⋆)μ_⋆GM_⋆/R_gR_⋆ for a n = 3 polytrope, |$L_\star \sim 4 \pi F_{\rm rad} {\rm R}_\star ^2$|⁠, and all scale heights are comparable to R_⋆ and derive

where the μ_⋆-dependent opacity κ_⋆ drops out in the above expression. With the accretion-wind equilibrium (⁠|${\dot{M}}_{\rm Bondi} \simeq {\dot{M}}_{\rm wind}$| and M_⋆ = M_eq) being maintained at a nearly constant λ_⋆ (slightly less than unity), the stellar mass |$M_{\rm eq}~\propto \mu _\star ^{-2}$| (i.e. A_μ ≃ −2).

There are several stages of stellar evolution (Section 4.2). The n = 1.5 polytrope model is a good approximation for stage 1 (MS with m_⋆ ≳ 485) and stage 3 (PostMS with m_⋆ ≲ 28) as in the case of fully convective stars and cores of red giants (Chandrasekhar 1939; Hansen, Kawaler & Trimble 2004). The equilibrium masses obtained from the mesa model (Fig. 3) indicate A_μ ≡ dlnm_eq/dlnμ_⋆ ≃ −1.5 for m_⋆ = m_eq  ≳ 500 (stage 1) and A_μ ≃ −2.3 for m_⋆ ≲ 30 (stage 3) which are in general agreement with the above analytic approximation. They also imply A_x ∼ 0.7 for the H-to-He conversion (MS stage 1) and A_y ∼ 0.8 → 0 for the He-to-α conversion (PostMS, stage 3) lead to the decline of equilibrium masses [equations (D3), (E4), (E5), (E6) and Fig. 2].

But for stage 2 (MS with m_⋆ ≃ 28–482), the n = 3 polytrope model is a better approximation. For any polytropic star (with P = Kρ^{(1 + n)/n}), the mass can be written as

where ρ_core is the star’s core density and |$\xi _1^2 {\mathrm{ d} \theta }/{\mathrm{ d} \xi })_{\xi _1}$| is a constant set by the polytropic index n and

is the adiabatic constant with the Avogadro and Boltzmann constants N_A and k_β, and the molecular weight for a fully ionized gas is μ_⋆ ≃ 4/(3 + 5X − Z). Taking these two equations for n = 1.5 polytrope to describe a fully convective star with |$\rho _{\rm core} \propto M_\star / R_\star ^3$|⁠, we find

Equations (E3) and (E9) lead to |$L_\star \propto M_\star ^{1.25} R_\star ^{0.75} \mu _\star ^4 \beta ^4$|⁠. Assuming |$R_\star \propto M_\star ^{0.6}$| as in our case, one finally gets

such that |$\lambda _\star \sim {\mathcal {O}} (1)$| is maintained with |$M_{\rm eq}~\propto \mu _\star ^{-5.7}$| (i.e. A_μ ≃ −5.7). This analytic approximation agrees well with the mesa result for equilibrium mass in the range of ∼30–485 (A_μ ≡ dlnm_eq/dlnμ_⋆ ∼ −4.7 during stage 2 in Fig. 3). It also implies A_x ≃ 2.1 → 0 and A_y ≃ −2.1 → −7.8 (as X_⋆ ≃ 0.5 → 0) and H-to-He conversion (MS stage 2) decreases the equilibrium masses [equations (D3), (E4), and (E10) and Fig. 2].

Based on the above analytic approximation, we can estimate the time-scale for changing X_⋆ by an amount ΔX_⋆ during MS stage 2 (as X_⋆ ≃ 0.5 → 0, μ_⋆ ≃ 0.72 → 1.3, A_x ≃ −A_y ≃ 2.6 → 0, and m_eq  ≃ 485 → 28) to be Δt_MS ∼ ϵ_Heτ_Sal ∼ 3 Myr (equation D5). This characteristic time-scale agrees well with that obtained with the mesa model. Minor discrepancies may be due to the nuclear burning core has only a fraction of the total stellar mass, dependence of A_mux on X_⋆ and A_muy on Y_⋆, and the negligence of the energy needed to power the loss of M_eq in the analytic approximation (equivalent to a slightly larger λ_⋆ in equation D5). In conclusion, the quasi-equilibrium mass loss time-scale is determine by the chemical-evolution time-scale in the core.