A Semi-Analytical Model for Stellar Evolution in AGN Disks
Alexander J. Dittmann Matteo Cantiello



Alexander J. Dittmann
dittmann@ias.edu

0000-0001-6157-6722]Alexander J. Dittmann
Department of Astronomy and Joint Space-Science Institute, University of Maryland, College Park, MD 20742-2421, USA
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

0000-0002-8171-8596]Matteo Cantiello
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

§ ABSTRACT 
Disks of gas accreting onto supermassive black holes may host numerous stellar-mass objects, formed within the disk or captured from a nuclear star cluster. We present a simplified model of stellar evolution applicable to these dense environments; our model exhibits exquisite agreement with full stellar evolution calculations at a minuscule fraction of the cost. Although the model presented here is limited to stars burning hydrogen in their cores, it is sufficient to determine the evolutionary fate of disk-embedded stars: whether they proceed to later stages of nuclear burning and leave behind a compact remnant, reach a quasi-steady state where mass loss and accretion balance one another, or whether accretion proceeds faster than stellar structure can adjust, causing a runaway. We provide numerous examples, highlighting how various disk parameters, and effects such as gap opening, affect stellar evolution outcomes. We also highlight how our model can accommodate time-varying conditions, such as those experienced by a star on an eccentric orbit, and can couple to N-body integrations. This model will enable more detailed studies of stellar populations and their interaction with accretion disks than have previously been possible. 






§ INTRODUCTION
 
The accretion disks that fuel active galactic nuclei (AGN) <cit.> may host numerous stars and stellar mass objects. Stars may form within gravitationally unstable regions of these accretion disk <cit.>, or be captured from nuclear star clusters <cit.>.  

Stellar populations in AGN disks may strongly affect the structure and observable characteristics of the accretion disks around supermassive black holes (SMBHs), particularly in the outer regions of the disk <cit.>. These populations may beget both electromagnetic <cit.> or gravitational wave <cit.> transients, and could be responsible for the supersolar metal content typically inferred in quasar broad line regions <cit.>. Closer to home, an accretion disk may have sculpted chemically and dynamically peculiar massive stellar populations in our own galactic center during the most recent accretion episode onto Sagittarius A^* <cit.>. 

Although numerous studies have considered the influences of embedded objects on accretion disks <cit.>, time-dependent models coupling accretion disks, stellar populations, and stellar evolution have yet to be constructed. Thus far, the cost of stellar evolution simulations has prohibited their self-consistent inclusion in such studies, although some stellar evolution studies have yielded insight into the typical fates of embedded stars as functions of various disk properties that can be applied to simplified models <cit.>. Our hope is that the model presented here or extensions thereof will enable more realistic modeling of the interplay between stars and accretion disks. 

In the following, we present a simple semi-analytical model of stellar evolution within AGN disks, coupling an Eddington standard model to prescriptions for accretion, mass loss, and related processes. We present the resulting set of differential equations in Section <ref>. Section <ref> provides examples of stellar evolution under various conditions, illustrating the close agreement between our model and various stellar evolution calculations. We discuss some potential future applications, model limitations, and potential extensions in Section <ref>, and summarize our work in Section <ref>. An implementation of our model is available at https://github.com/ajdittmann/starsamhttps://github.com/ajdittmann/starsam.




§ SIMPLIFIED STELLAR MODELS

Our model of stars in AGN disks consists of two main components: equations approximating how the structure of a massive star depends on its bulk properties (its mass and composition); and equations describing how the bulk properties of the star evolve over time (through accretion, mass loss, etc.). Concerning the former, we draw on previous calculations of the structure of massive stars presented in <cit.>, which we review in Section <ref>. Concerning the latter, we draw upon the treatments employed in previous full calculations of stellar evolution in AGN disks, such as <cit.>, which we review in Section <ref>. 



 §.§ Stellar Structure

The fundamental quantities in our stellar models are the total masses of hydrogen, helium, and metals constituting the star (M_X, M_Y, and M_Z), which are used to calculate the total stellar mass and mass fractions 

    M_*=M_X+M_Y+M_Z,
    
    X_*=M_X/M_*,   Y_*=M_Y/M_*,   Z_*=M_Z/M_*.

We assume that the star is chemically homogeneous. In nonrotating stars, this may not be satisfied until M_*≳100 M_⊙, at which point the stars become almost fully convective <cit.>. On the other hand, rotating stars are subject to numerous additional instabilities that lead to chemical homogeneity <cit.>. Since stars in AGN disks are likely to spin up as they accrete <cit.>, the assumption of chemical homogeneity is likely to be well-justified even at lower masses. Homogeneity makes these stars particularly well-described by simplified models (as noted in , the structure of stars in AGN disks tend to be well-described by n=3 polytropes). In the following, we describe the Eddington standard model for massive stars <cit.>, drawing heavily from Section 2b of <cit.>. 

The average number of nuclei and electrons per baryon within the star, the inverse of the mean molecular weight, is given by 

    Y_T=μ^-1=1/4(6X_*+Y_*+2).

The photon entropy (s_γ) is then related to the ratio of radiation pressure to gas pressure (σ≡ P_ rad/P_ gas) by

    σ = s_γ/(4Y_T)=(β^-1-1),

where β≡ P_ gas/(P_ gas+P_ rad). In these terms, the familiar form of Eddington's quartic equation 

    1-β=0.003β^4μ^4(M_*/M_⊙)^2

can be rewritten as 

    M_* = 1.141 s_γ^2(1+4Y_T/s_γ)^3/2 M_⊙.

The stellar luminosity L_* is related to s_γ through the Eddington ratio, Γ≡ L_*/L_ Edd, where

    Γ = (4Y_T + 1)^-1
    
    L_ Edd≡ 1.2 × 10^38 Y_e^-1 (M_*/M_⊙)  erg s^-1,

and where Y_e = (1+X)/2 is the number of electrons per baryon. Notably, the Eddington standard model discussed here yields an appropriate mass-luminosity relation L_*∝ M_*^3 at ∼solar masses and L_*∝ M_* at higher masses, and it is our assumption of uniform composition and efficient mixing that breaks down in low-mass stars. 

In order to employ the mass-loss estimates described in Section <ref>, it is necessary to compute an estimate of the stellar radius. In the (n=3) polytropic approximation, the stellar radius is <cit.>

    R_* = 30.4Y_Tσ^1/2(1+σ)^1/2T_c^-1 R_⊙,

where T_c is the central temperature of the star in keV. To estimate the central temperatures, we follow appendix A of <cit.>, assuming that energy generation is solely due to hydrogen burning via the CNO cycle. Specifically, Equations A2 and A3 of <cit.> provide

    T_c/2.67 =(1 - 0.021 logA - 0.021 logX_N^-1 + 0.053logT_c)^-3
    
    A = s_γ[1+(ν_p-3)(1+ς)/6]^3/2[XY_e(1+σ^-1)]^-1,
    ν_p = 22.42T_c^-1/3 + 7/3
    ς = (2s_γ/Y_N+1)^-1.

The above expressions follow from a power-law fit to the CNO nuclear burning energy generation rate, where ν_p is the power-law exponent with respect to temperature in the neighborhood of the central temperature and X_N is the nitrogen mass fraction. The ς term is a small modification to the results of <cit.> accounting for variation of photon entropy within the convective core <cit.>. We approximate the nitrogen mass fraction by scaling the solar <cit.> nitrogen, carbon, and oxygen abundances by Z_*/Z_⊙ and assuming that they are entirely in the form of nitrogen, which is appropriate given the rate-limiting step of the CNO cycle ^14_7N+^1_1H→ ^15_8O+γ.



 §.§ Bulk Evolution

The bulk composition of the star evolves due to three processes: accretion (Ṁ_+), mass-loss (Ṁ_-), and fusion. In the following, we have specialized to a few specific forms of each of these, but it is trivial to extend our model to other treatments or approximations.  


  §.§.§ Accretion

In the following, we have explored four regimes of stellar accretion: spherically symmetric Bondi accretion, tidally-limited accretion, accretion limited by gap formation in an accretion disk, and accretion onto stars with substantial velocity relative to the disk, each of which determine a base accretion rate, Ṁ_0. For a star with zero luminosity, the accretion rate Ṁ_+=Ṁ_0. However, as the star approaches the Eddington luminosity, the accretion rate is reduced by some factor 0≤ f ≤ 1, according to 

    Ṁ_+ = Ṁ_0f(L_*+L_S),

where L_* is the intrinsic stellar luminosity and L_S≈Ṁ_+GM_*/R_* is the luminosity generated by shocks as accreted material strikes the surface (see Section 4.5 of ).[One notational difference between this work and <cit.> is that <cit.> defined L_* to include both the intrinsic stellar luminosity and the shock luminosity, while here we define these separately. ]^,[The assumptions made in <cit.> may break down when the accretion flow is sufficiently optically thick <cit.>.]

When considering spherically symmetric Bondi accretion, the accretion rate is simply given by <cit.>

    Ṁ_0 = πρ c_s R_B^2
    
    R_B ≡ 2GM_*/c_s^2,

where ρ and c_s are the density and sound speed of the AGN disk around the star, and the Bondi radius (R_B) characterizes the distance from the embedded star within which the gravity of the star can overcome the thermal response of the fluid. If modeling the evolution of a star moving with significant velocity relative to the background fluid (v≳ c_s), then the Bondi-Hoyle-Lyttleton accretion rate may be more appropriate, leading to <cit.>

    Ṁ_0 = 4πρG^2M_*^2/(c_s^2 + v^2)^3/2.


Within an accretion disk, in the presence of an SMBH, gas simply being within the Bondi radius is not a sufficient condition for it to accrete onto the star, as that gas may still be gravitationally bound to the SMBH. The Hill radius <cit.> defines the region within which gas is bound to the star rather than the SMBH, and can be defined in terms of the angular velocity of the stellar orbit about the SMBH as

    R_H = (GM_*/3Ω^2)^1/3.

For a given c_s and Ω, R_B < R_H at low stellar masses, and R_H < R_B at higher masses. When we consider “tidally-limited” accretion, we simply take whichever of these is smaller as the limiting scale, using <cit.>

    Ṁ_0 = πρ c_s min(R_B^2,R_H^2).

It is also possible to generalize Equation (<ref>) to account for the stellar Hill radius and the vertical extent of the disk, as in <cit.>, prescribing accretion according to 

    Ṁ_0 = 4πρ R_ accmin(R_ acc,H)
    
    R_ acc=GM_*/v_ eff^2
    
    v_ eff = (c_s^2 + R_H^2Ω^2 + v^2)^1/2.


If mass ratio between the star and the SMBH (q≡ M_*/M_∙) is sufficiently large (for a given disk scale height and effective viscosity) then the presence of the star can begin to alter the large-scale structure of the accretion disk. We consider a simple scenario, where this interaction decreases the average surface density of the disk by some factor ξ≡Σ'/Σ_0, the ratio of the perturbed to unperturbed surface densities, which in turn decreases the accretion rate by a factor of ∼ξ <cit.>. We simply approximate that <cit.>

    ξ = 1/1+0.04K
    
    K≡ q^2/(h^5α),

although stars may open even deeper (more depleted) gaps in very low-viscosity disks <cit.>. Thus, if we consider “gap-limited” accretion, we take the base accretion rate to be

    Ṁ_0 = πρ c_s min(R_B^2,R_H^2)ξ(α,M_∙,h).


Thus, after choosing some approximation for the unmodified (by radiation) accretion rate, we are left to make some approximation for how radiative feedback further limits accretion. Here, we limit ourselves to an approximation introduced in <cit.>, inspired by multidimensional flows where accretion and outflows may coexist. In this case, 

    f(L_*+L_S) = 1-tanh[(L_*+L_S)/L_ Edd].

Thus, the accretion rate onto the star can be calculated given its mass, luminosity, and some number of extrinsic parameters dependent on calculation one aims to conduct. 



  §.§.§ Mass Loss

As stars approach the Eddington limit (L_*∼ L_ Edd), they are thought to lose appreciable mass via continuum-driven winds <cit.> as radiation pressure overcomes the pull of gravity. In this work, we follow a common approach used in stellar evolution simulations <cit.>, assuming an outflow at the escape velocity (v_ esc=(2GM_*/R_*)^1/2) only when the stellar luminosity is very near the Eddington limit, specifically

    Ṁ_- = L/v_ esc[1 + tanh(L_*+L_S-L_ Edd/0.1L_ Edd)].

Although we have only investigated the above prescription in this work, our formalism can accommodate any mass loss prescription that only depends on bulk properties of the star. 



  §.§.§ Nuclear Burning

We treat nuclear burning in a greatly simplified manner, approximating energy generation due to hydrogen fusion alone. Unlike full stellar evolution calculations, we cannot extend our simplified model beyond core hydrogen burning. A sufficient approximation is that for each helium atom that is produced, four protons are consumed and 27 MeVs are liberated, contributing to the intrinsic stellar luminosity. Thus, helium and hydrogen are created and depleted at a rate

    Ṁ_He,L=-Ṁ_H,L=L_*/27 MeV/4 m_p,

where m_p is the proton mass.



 §.§ Piecing the Model Together

To model the evolution of a star within an AGN disk or other dense environment, one must select initial properties of the star, specifically its initial mass, hydrogen abundance, helium abundance, and metal abundance. One must also select at least an ambient density ρ, an ambient sound speed c_s, and the mass fractions of the ambient gas (X, Y, Z) — in principle these could be treated as constants, or functions of time. One must then decide which approximation to make of the stellar environment. The above parameters are sufficient to model spherically symmetric  accretion, but if one would like to study stars within accretion disks, one must also chose some angular velocity for the stellar orbit Ω; and possibly a disk aspect ratio h, viscosity parameter α, and SMBH mass M_∙ if one aims to account for gap opening — in our models, we allow Ω and h to be constants or functions of time, but assume α and M_∙ to be constant.

Once these choices have been made, simulating stellar evolution in AGN disks boils down to solving a trio of ordinary differential equations (ODEs) describing how the hydrogen, helium, and metal content of the star evolve over time, and at each step solving a series of nonlinear equations to determine the rate of change of each component. Ultimately, we arrive at the system

    dM_H/dt = XṀ_+ - X_*Ṁ_- - L_*/27 MeV/4 m_p
    dM_He/dt = YṀ_+ - Y_*Ṁ_- + L_*/27 MeV/4 m_p
    dM_Z/dt = ZṀ_+ - Z_*Ṁ_-,

where L_*, Ṁ_-, and Ṁ_+ can be calculated following section <ref> and <ref>. These equations can be solved using standard methods for solving ODEs and a nonlinear one-dimensional root-finding algorithm. 





In addition to a set of ODEs and their initial conditions, we must decide when to terminate each integration. We consider three potential outcomes of stellar evolution. First, stars may simply run out of hydrogen fuel, at which point they would go on to later stages of nuclear burning on a comparatively rapid timescale, likely leaving behind some form of compact remnant <cit.>; however, the models presented here are not appropriate to stages of stellar evolution beyond core hydrogen burning, so we terminate each integration if M_H≤0. Additionally, it is possible that stellar accretion and mass loss might balance each other, leading to a quasi-steady “immortal” state <cit.>; the integration of these models typically ends upon reaching a user-imposed time limit. Additionally, it is possible that stars might accrete more quickly than they can thermally adjust; we do not attempt to model the evolution of these stars subsequent to reaching such a condition, although it seems reasonable that these stars might continue to increase their mass without appreciable radiative feedback until they accrete the entirety of the local gas supply <cit.>.
In the following calculations, we have estimated the thermal relaxation timescale of massive stars to be 

    τ_ KH=3/2GM_*^2/R_*L_*,

appropriate for a star with n=3 polytropic structure. We assume that the star undergoes runaway accretion when τ_ acc≡M_*/Ṁ_+<τ_ KH. However, it would be trivial to incorporate a different runaway condition if desired. 
The aforementioned stellar evolution outcomes are illustrated in Figure <ref>, where the higher-density model shown in yellow-green undergoes runaway accretion, the medium-density model shown in blue reaches an immortal state, and the lower-density model shown in blue runs out of hydrogen fuel. 

As a word of caution, in some situations it may be unreasonable to terminate calculations when τ_ acc<τ_ KH without further consideration. For example, for near-solar masses L_*∝ M^3 so that τ_ KH decreases with mass, but if accretion is limited by tidal effects or gap opening, τ_ acc increases with mass. This interplay at low masses might lead to the early termination of a calculation of a star than undergoes a brief period of rapid accretion prior to the accretion timescale increasing. If accretion is limited by both gap opening and tidal effects, the accretion rate (Ṁ_0) itself decreases with stellar mass. 



§ EXAMPLES

We present here a number of example calculations that can be compared to existing stellar evolution calculations. Each of these was computed using a publicly available implementation of our model, which can be found at https://github.com/ajdittmann/starsamhttps://github.com/ajdittmann/starsam.[Scripts to reproduce each of our calculations in this paper will be made public upon publication.] 
Unless otherwise noted, we use 
M_*(t=0)=10M_⊙, X_*(t=0)=0.74, Y_*(t=0)=0.24, X=0.72, Y=0.27, and c_s=10^6  cm s^-1.
In most of the following, we have used the basic RK4(5) integrator <cit.> as implemented in scipy <cit.>; in Section <ref>, when coupling the above differential equations to a short term orbital integration, we use the standard 4th-order method due to <cit.>. 



 §.§ Tidally Limited Accretion

We begin by comparing our models to existing stellar evolution calculations. In particular, <cit.> investigated a number of approximations for how deviation from spherical symmetry in AGN disks could affect stellar evolution, and found that the tidal limiting described by Equation (<ref>) was often significant and the most widely applicable. 

We present in Figure <ref> a number of M_*(t) tracks, paralleling the calculations in <cit.>, although assuming very slightly different mass fractions. We find very good agreement between the results of our simplified model and full stellar evolution calculations. For example, at both Ω=10^-10  s^-1 and Ω=10^-12  s^-1 we find that models transition from exhausting their helium supply to reaching a quasi-steady state at the same ambient density, although the exact shape of the M_*(t) curve deviates slightly in the ρ=10^-16  g  cm^-3, Ω=10^-10  s^-1 case.

Overall, given the enormous simplifications of the present model, the level of agreement they achieve with full stellar evolution calculations is extremely encouraging. They seem to reproduce the typical evolutionary outcomes of a given star, as well as their masses, to good precision. For example, the masses of immortal stars in Figure <ref> usually deviate from those in the relevant models of <cit.> by factors of at most a factor of a few. 

One perhaps counter-intuitive result displayed in Figure <ref>, particularly the simulations at Ω=10^-12 s^-1, is that immortal stars embedded in high-density media reach lower masses than stars embedded in comparatively lower-density media. In these cases, the accretion shock luminosity can reach a substantial fraction of the total stellar luminosity, resulting in substantially higher mass loss rates at a given mass. This trend has been observed in previous full in stellar evolution simulations <cit.>. 



 §.§ Gap Opening

The possibility of stars opening gaps in AGN disks and subsequently decreasing the rate of accretion onto embedded stars was discussed in <cit.>, although they did not carry out any calculations since Equation (<ref>) illustrates that gap opening will only be relevant in disks that are either very thin h≲ 10^-3, nearly inviscid α≲ 10^-4, or those orbiting very low-mass SMBHs M_∙≲ 10^6 M_⊙. Although gap opening is more important for intermediate-mass black holes and extremely massive (M_*≳10^4) stars <cit.>, it can also be relevant for stars and stellar-mass black holes around less-massive SMBHs <cit.>.



Figure <ref> illustrates the effects of gap opening explicitly for disks with α=0.1, c_s=10^6  cm s^-1, and Ω = 10^-10  s^-1. In the top panel, showing calculations assuming an SMBH mass of 10^6 M_⊙, the disk aspect ration h=H/r strongly affects the evolution of embedded stars. In the h=10^-3 case, accretion is entirely inhibited, since initially K=10^6. In the h=10^-2 case, K was initially ∼10, so accretion was not suppressed until the stars reached somewhat higher masses. One other hand, in h=0.1 disks, stellar evolution was unaffected. 

The bottom panel of Figure <ref> illustrates the same calculations for a more massive SMBH. In this case, the h=0.1 and h=10^-2 disks lead to nearly indistinguishable stellar evolution, while accretion onto stars in the h=10^-3 case is heavily suppressed. Intriguingly, in the ρ=10^-17  g cm^-3 case, gap opening actually leads to a higher-mass immortal star by decreasing the accretion rate as discussed in Section <ref>. It is conceivable that gap opening might lead to disparate stellar populations around low- and high-mass SMBHs. However, we leave a detailed investigation of this and its observational consequences for future studies. 





 §.§ The Critical Density to Sustain Immortal Stars

One of the questions we hope to answer with our models is what environmental conditions lead to which stellar evolutionary outcomes — what sort of disk parameters should lead to stars reaching quasi-steady states, which should lead to runaway accretion, and which should should produce stars poised to undergo supernovae? Here, we examine the boundary between stars that reach a quasi-steady state and those that run out of hydrogen fuel. This depends on whether or not the stellar nuclear burning timescale can become shorter than the accretion timescale prior to accretion and mass loss achieving a quasi-steady state (see Figure 1 of ). Figure <ref> collects calculations determining this critical density in a number of scenarios using 25 bisection iterations at each point.

The top panel of Figure <ref> examines  the effect of the background helium abundance of the AGN disk on the critical density for immortal stars undergoing spherically symmetric (Bondi) accretion. Here, because stars with higher mean molecular weights tend to be more luminous at constant mass than those with lower mean molecular weights, stars that accrete material with a higher helium abundance will have shorter nuclear burning timescales at a given mass, requiring higher densities to achieve a short enough accretion timescale to sustain a star in a quasi-steady state. The results of this calculation agree very closely with the analogous stellar evolution calculations presented as the orange curve in Figure 5 of <cit.>. 

The middle panel of Figure <ref> plots the critical density as a function of angular velocity and sound speed for stars undergoing tidally-limited accretion. At sufficiently low angular velocities, accretion occurs in a spherically symmetric manner, requiring much higher densities at high sound speeds, since R_B∝ c_s^-2. However, at higher angular velocities, the accretion rate is proportional to the sound speed, so lower densities are required to sustain immortal stars at fixed Ω. The orange curve, for c_s=10^6  cm s^-1, agrees very well with the analogous result based on stellar evolution calculations (Figure 8 and Equation 22 of ) in the Ω≳10^-11  s^-1 regime where those simulations were conducted. 

The bottom panel of Figure <ref> shows the critical density as a function of disk aspect ratio and viscosity for a star orbiting a 10^8 M_∙ SMBH at Ω=10^-10  s^-1 within a c_s=10^6  cm s^-1 disk. The values of α chosen here roughly bracket those used to model compact object accretion disks <cit.> and the effective values of α derived from global simulations of the magnetorotational instability <cit.>. At this SMBH mass, gap opening makes it more difficult to sustain stars in a quasi-steady state in very thin disks (h<10^-2), although it is still quite possible.  



 §.§ A Star on an Eccentric Orbit

The model presented here should facilitate coupled N-body and stellar evolution calculations. To illustrate this capability, we study here the evolution of a 2-body star-SMBH system and a toy disk model. 
Here, we assume an initially 10 M_⊙ star on a prograde eccentric (e=0.4) orbit around a 10^7 M_⊙ SMBH with a period of ∼ 1  Myr. The star moves though an isothermal Keplerian accretion disk with c_s=10^5 cm s^-1 and a density profile given by
ρ(r) = ρ_0(r/r_0)^-1, 
normalized to 3×10^-13g cm^-3 at 10 AU. 



Figure <ref> illustrates the results of an example calculation over four orbital periods under the assumptions of Equations (<ref>) and (<ref>), where we have neglected the change in stellar velocity due to accretion or drag forces, which over this dynamically short integration have little effect. In this case, the accretion timescale is shorter than the orbital timescale and the star essentially moves between two quasi-steady states of differing mass. Notably, in a multi-body system, scattering events could easily drive stars into regions where their evolution is irrevocably altered. 





 §.§ Disk Models

Our model will enable more precise feedback between models of AGN disks and the stellar populations they harbor. To illustrate this possibility, albeit without incorporating the effects of stellar feedback on the disk, we illustrate in Figure <ref> the stellar evolutionary outcomes of a set of 8×10^5 calculations; each simulation assumed the star to reside at a fixed location in a set of <cit.> disk models calculated in <cit.>, allowed a maximum duration 10^7 years, and assumed tidally-limited accretion (without considering gap opening). We classify stars as `runaway' (τ_ acc<τ_ KH), `immortal' (reaching a quasi-steady state), `massive' (achieving a mass greater than 8 M_⊙ and accreting at too slow a rate to achieve a quasi-steady state), and `ordinary' otherwise. A shorter calculation, effectively a smaller upper limit on the AGN disk lifetime, would primarily limit the number of higher-mass stars. The calculations in this subsection assumed an initial stellar mass of M_*=M_⊙, which greatly increased the initial accretion timescale compared to the calculations in prior sections.

The same general trends occur at each SMBH mass, although the extent of each region varies with mass. Gas densities are typically very low at large radii, usually resulting in very low accretion rates that do not affect stellar evolution significantly. Moving inwards, gas densities increase and the accretion timescale can become comparable or shorter than the initial stellar nuclear burning timescale and stars grow massive. This trend continues, and stars eventually accrete quickly enough to reach quasi-steady states or undergo runaway accretion. However, at smaller radii, the gas density varies non-monotonically, while the disk temperature increases and tidal forces from the SMBH become stronger, limiting accretion. 

Our model predicts slightly fewer massive and immortal stars than the analytical calculations of <cit.> but a similar number of stars undergoing runaway accretion. This is most likely because of the slightly different criteria used to differentiate between `ordinary' and `massive stars,' <cit.> comparing the initial accretion timescale to the expected disk lifetime in contrast to our cut based on the final stellar mass. Notably, more massive stars have shorter accretion timescales, so these results depend on our choice of initial stellar masses. 



§ DISCUSSION

The models presented in this work, thanks to their close agreement with full stellar evolution simulations at a minuscule fraction of the cost, will enable new avenues of research detailed below, which we hope will facilitate rigorous comparisons between theory and observation. Specifically, these models should make it possible to carry out time-dependent studies of stellar populations in AGN disks with minimal sacrifices compared to full stellar evolution calculations. Furthermore, the model presented here may be extended in a number of ways, to broaden its applicability to other stages of stellar evolution, or to incorporate other physical effects. 
However, when applying the present models, one should keep in mind some of the assumptions used to derive them and the resulting model limitations, which we collect below.



 §.§ Applications

The model developed here should make it possible to conduct detailed time-dependent studies of stellar populations in AGN disks, circumventing the cost of full stellar evolution simulations. For example, numerous studies have coupled N-body or smoothed-particle-hydrodynamics codes to either full <cit.> or approximate <cit.> stellar evolution. Because the model developed here is a relatively straightforward ordinary differential equation, it should be straightforward to couple it to N-body calculations. Presumably, an accretion prescription along the lines of Equation (<ref>) or Equation (<ref>) would be appropriate for such studies. 

Because the model presented here self-consistently predicts stellar luminosities and mass loss rates given local disk properties, it should improve models of stellar feedback on AGN disks <cit.>. It may also be possible to incorporate these stellar models into time-dependent disk model calculations <cit.>, or one-dimensional Monte Carlo calculations similar to those used to predict stellar-mass black hole merger rates in AGN disks <cit.>. 



 §.§ Caveats

The models presented here rely on numerous assumptions that may cause their predictions to deviate from reality. Many of these, such as the highly approximate prescriptions for how stars accrete and lose mass, and how the accreting matter interacts with a near-Eddington radiation field, follow from the necessity of describing intrinsically three dimensional and turbulent effects within a simplified ODE model. Such ad hoc assumptions are also shared by full stellar evolution calculations; although such assumptions cannot be circumvented, their effects on predictions can be gauged through parameter studies, which are far more efficient using the model presented here. 

Compared to full stellar evolution calculations, the simplified model presented in Section <ref> has primarily sacrificed spatial and chemical fidelity. The assumption of CNO-dominated burning limits its applicability to stars more massive than the Sun and with appreciable metallicity (making them ill-suited to studying stars in the very early Universe). The assumption of an n=3 polytropic structure is quite suitable for stars embedded in AGN disks <cit.>, but is nevertheless a limitation. The assumption of chemical homogeneity is not reasonable for stars of all masses: below M≲100M_⊙, nonrotating stars still possess a radiative envelope, which may preclude freshly accreted material from mixing into the stellar core <cit.>. However, stars accreting from AGN disks are likely to reach substantial rotation rates <cit.>, at which point rotation-induced mixing should lead to chemical homogeneity at virtually all stellar masses <cit.>.



 §.§ Extensions

Many of the model components presented here, such as the recipes for the accretion rate onto stars given external accretion disk conditions, the interplay between near-Eddington radiation fields and accretion, and continuum-driven mass loss, are quite flexible. Similarly one may choose to implement different criteria for runaway accretion. For example, in the theoretical study of early-Universe supermassive black hole formation, some have considered the balance between the accretion and Kelvin-Helmholtz timescales of the stellar envelope specifically, as opposed to the global timescales employed here <cit.>. 

Beyond adjusting the components already present within the model, it is possible to incorporate other considerations entirely. For example, <cit.> has suggested modifications to the stellar accretion rate applicable when the accretion flow onto the star is very optically thick. Incorporating an accretion rate along the lines suggested in Section 3.2 of <cit.>, suitably modified when R_H<R_B as discussed in Section <ref> <cit.>, would be a straightforward way to extend the applicability of our model to other regimes. Judging by Figure <ref>, limiting the stellar accretion rate has the potential to substantially increase the number of massive and immortal stars within the disk

One natural extension of the models presented here is to explicitly account for stellar angular momentum. Thanks to the background shear flow of AGN disks, accretion may spin stars up to considerable angular velocities. Prescriptions for these effects have been incorporated into previous stellar evolution calculations <cit.>, which would be straightforward to incorporate into the semi-analytical model presented here by introducing an additional differential equation for the total stellar spin. One crucial ingredient for such considerations is the stellar moment of inertia, which for the Eddington standard model discussed here is approximately 

    I_*≈ 7.529×10^-2 M_*R_*^2.


One of the limitations of the model presented in this work is its restriction to burning via the CNO cycle. At least under the assumption of chemically homogeneous evolution <cit.> which seems fairly applicable to the evolution of stars in AGN disks <cit.>, it may be reasonable to assume energy generation due to triple-α burning after hydrogen exhaustion. Such considerations might improve the estimates made by these models of the amount of metals ejected into the ambient medium. However, modeling later stages of stellar evolution would require evolving the abundances of more elements, at which points the costs would likely outweigh potential benefits. 
Moreover, the time a star spends burning heavier elements is significantly shorter than the time it spends burning hydrogen, which usually accounts for 90% of the stellar lifetime. Later stages of burning make up a smaller fraction of the the lives of AGN stars, where accretion and mixing result in relatively longer hydrogen burning phases.



§ CONCLUSIONS

In this work, we have combined a simple model of stellar structure, an Eddington standard model <cit.>, with approximations for the interplay between stellar radiation, accretion, and mass loss and AGN disks employed in previous stellar evolution calculations <cit.> to yield a simple model of stellar evolution in AGN disks or other dense environments. As shown in Section <ref>, our model very closely reproduces the results of full stellar evolution simulations, typically at less than a thousandth of the cost, although the reduced model is only applicable to core hydrogen burning via the CNO cycle in its present form. We hope that our distillation of the equations of stellar evolution in AGN disks to a trio of ordinary differential equations will enable future studies to investigate the interplay between stellar populations and their environments in much more detail than had been possible before now. 



§ SOFTWARE

 <cit.>,  <cit.>, scipy <cit.>



§ ACKNOWLEDGMENTS

AJD is grateful for the hospitality of the Center for Computational Astrophysics throughout a visit during which this work was initiated.
The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.
AJD gratefully acknowledges support from LANL/LDRD under project number 20220087DR, and NASA ADAP grants 80NSSC21K0649 and 80NSSC20K0288. The LA-UR number is LA-UR-24-28455.

aasjournal


