The Population of Massive Stars in Active Galactic Nuclei Disks

Gravitational instability in the outskirts of active galactic nuclei (AGN) disks leads to disk fragmentation and formation of ∼300 M ⊙ supermassive stars with potentially long lifetimes. Alternatively, stars can be captured ex situ and grow from gas accretion in the AGN disk. However, the number density distribution throughout the disk is limited by thermal feedback as their luminosities provide the dominant heating source. We derive equilibrium stellar surface density profiles under two limiting contexts: in the case where the stellar lifetimes are prolonged, due to the recycling of hydrogen-rich disk gas, only the fraction of gas converted into heat is removed from the disk accretion flow. Alternatively, if stellar composition recycling is inefficient and stars can evolve off the main sequence, the disk accretion rate is quenched toward smaller radii resembling a classical starburst disk, albeit the effective removal rate depends not only on the stellar lifetime, but also the mass of stellar remnants. For AGNs with central supermassive black hole masses of ∼106–108 M ⊙ accreting at ∼0.1 Eddington efficiency, we estimate a total number of 103–105 massive stars and the rate of stellar mergers to be 10−3 to 1 yr−1. We initiate the detailed study of the interaction between a swarm of massive stars through hydro and N-body simulations to provide better prescriptions of dynamical processes in AGN disks, and to constrain more accurate estimates of the stellar population.


INTRODUCTION
Supermassive Black Holes (SMBH) are known to host accretion disks that provide power to active galactic nuclei (AGN) in the centers of massive galaxies (Lynden-Bell 1969).Strong radiative cooling in the outer regions of these disks lead to fragmentation via gravitational instability (Goodman 2003;Goodman & Tan 2004), especially in a radiation pressure dominated environment (Jiang & Goodman 2011;Chen et al. 2023b).Subsequently, the fragments can grow from collision and gas accretion to become massive stars.Alternatively, stars in the nuclear cluster can be captured onto the disk and continue to grow by gas accretion (Artymowicz et al. 1993;MacLeod & Lin 2020;Wang et al. 2023).
The onset of in-situ fragmentation, or star formation, is associated with the fact that heating by MRI and gravito-turbulence are insufficient to maintain energy balance against radiative cooling.Assuming star formation provides the necessary auxiliary power to maintain energy equilibrium at marginal gravitational instability (Toomre parameter Q ∼ 1), Thompson et al. (2005) constructed a star-burst disk model in which the accretion rate Ṁd should decrease towards the SMBH as mass flux is continuously dumped into stars.However, considering that stars embedded in AGN disks can form from massive self-gravitating clumps or grow up to a few hundred M ⊙ from accretion of external gas, their evolution and fates can be significantly different from field stars.In particular, recent works have pointed out that since the accretion disk evolves on a much longer Salpeter timescale (∼ 10 8 years), embedded massive stars may replenish their hydrogen reservoir by recycling fresh gas from the disk within a few parsecs of the disk (Cantiello et al. 2021;Jermyn et al. 2022).Consequently, they undergo little chemical/nuclear evolution and may stay on the hydrogen-burning main sequence for a timescale that is comparable to AGN lifetime τ AGN .These "immortal stars" (IMS) become efficient nuclear power plants that burns whatever amount of hydrogen they accrete from the ambient gas into heat and return all helium byproducts back into the disk.In fact, Jermyn et al. (2022) focused their analysis of massive stars to the inner ∼ 0.03 pc (since that is where disk-capture mechanism is most efficient) and already infer an average count of ∼ 10 2 − 10 4 coexisting immortal stars within an AGN disk, while it should be further realized that star formation is in fact capable of producing massive stars over a wider range of distances (Jiang & Goodman 2011;Chen et al. 2023b).In the IMS sce-Chen & Lin nario, after the formation of the first generation of stars, only a small fraction of gas accreted on to these stars is converted into heat to balance radiative cooling, while the rest is recycled out back into the disk, leading to a negligible reduction in disk accretion rate.
Nevertheless, the IMS track may not be the whole picture.Supporting evidences for ongoing stellar and chemical evolution in AGN disks are provided by the observed red-shift independent super solar metallicity of AGNs' broad line regions and their enhanced heavyelement abundance, relative to the narrow emission line regions (Huang et al. 2023).These properties are inconsistent with the implication of "immortal stars" powered by CNO burning on the main sequence without any net heavy element production.However, if the accreted gas is buffered by a radiative envelope and is not well mixed with the nuclear burning zone deep in the stellar interior, the accumulation of He ashes from CNO burning would lead to the onset of efficient He burning, the transition to post main sequence evolution, and eventually collapse into black holes or supernova explosion with prolific heavy-element yields (Ali-Dib & Lin 2023).In this ongoing "Stellar Evolution and Pollution in AGN Disk" (SEPAD) track, massive stars still accrete and radiate just like the IMS model, but the overall efficiency of stars in extracting gas density from the disk to convert into luminosity is different due to their limited lifespan, and in turn the disk accretion rate is regulated in a different manner.
The coexistence and concentration of a large population of massive stars, along either the IMS or SEPAD tracks, raises the probability of stellar mergers.Stellar coalescence events may continuously increase individual stars' mass and luminosity.If this process proceeds rapidly, an energy equilibrium would not be sustainable at all times.
In this paper, we study the population of massive stars in an AGN based on consistent disk models at thermal equilibrium, in order to provide constraints on their number distribution and merger rates.We formulate the main equations for AGN disk models and stellar heating in §2.We quantify key differences between IMS and SEPAD tracks and their influences to the disk structure.In §3, we estimate the dynamical merger timescales of stellar coalescence as functions of disk and star properties.In §4, we solve specific disk models in the IMS and SEPAD limits as well as the merger timescales associated with their stellar density profiles.In §5, we summarize our findings and discuss future prospects for subsequent numerical work that could refine our model of massive star and their impact on the AGN disk.

General equations
We illustrate our problem based on a parameterized AGN disk model, similar to what is applied in MacLeod & Lin (2020); Davies & Lin (2020).A SMBH of mass M • is fed by a mass flux of where With the conventional α ν prescription (Shakura & Sunyaev 1973) for viscosity (ν = α ν c 2 s /Ω where α ν is a factor representing angular momentum transport efficiency, c s is the sound speed, Ω = GM • /R 3 , P = 2π/Ω are the Keplerian angular velocity and orbital period, at a distance of R, respectively.The disk accretion rate can be written as where Q ≡ c s Ω/πGΣ is the gravitational stability parameter (Safronov 1960;Toomre 1964).The aspect ratio is h = H/R where H is the scale height.The radial velocity is given by If accretion is driven by magneto-rotational instability (MRI) with viscous α ≲ 10 −2 (Balbus & Hawley 1991, 1998;Bai & Stone 2011;Deng et al. 2019) and the radiative cooling is balanced by accretion heating alone, gravitational instability (GI) would occur (with Q ≲ 1) in region with R ≳ R Q=1 ≃ 10 −2 pc around SMBHs with m 8 ∼ 0.01 − 1 (Paczynski 1978).GI leads to gravitoturbulence with α ν ∼ 0.01 − 1 (Lin & Pringle 1987;Lin et al. 1988;Kratter et al. 2008;Zhu et al. 2010b,a;Martin & Lubow 2011;Rafikov 2015;Deng et al. 2020) and prolific rates of in situ star formation (Goodman 2003;Goodman & Tan 2004).
In star-forming disk regions of interest, viscous heating from turbulence is much weaker than radiative cooling, and an energy balance needs to be maintained by radiation from either a permanent population of IMS stars or the persistent formation of SEPAD stars.The latter mode may lead to a strongly (radially) variable Ṁd (see §3).Nevertheless, we still expect the disk to self-regulate and maintain Q ∼ 1 everywhere due to the requirement for a continuous replenishment of heat sources.The viscous stress becomes decoupled from the energy equation and is simply a parametrization of the inflow velocity V R from a maximum level of turbulence that the disk is able to support.Assuming Q = 1 and Ṁd ≃ 2πΣV R R at all R, we determine the aspect ratio and surface density: In outer regions, the sound speed c s = hV K becomes flat for local instability with a given α ν and Ṁd (Sirko & Goodman 2003), where V K = ΩR is the Keplerian velocity.
For a vertically self-similar disk model, the mid-plane density while the midplane characteristic temperature T c satisfies the Eddington quartic when optically thick As we show below, the outer regions of such disks are self-consistently radiation dominated, since radiation to gas pressure fraction is an increasing function of r (the relation is reversed if turbulence provides heating in a constant α ν , gravitationally stable standard disk (Shakura & Sunyaev 1973), which is beyond the scope of our paper).
In dense regions where the disk is optically thick, we can express the radiative cooling rate from two sides of the disk Q − = 2σT 4 eff (to be differentiated with the Toomre Q parameter) as a function of these variables (Hubeny 1990;Jiang & Goodman 2011;Chen et al. 2023b): where κ is a constant opacity set to order unity in cgs unit.This simple prescription approximates the grain opacity within 10 2 − 10 3 K and the electron opacity above 10 4 K (Bell & Lin 1994).We discuss the possible effect of an opacity window close to the grain sublimation temperature in §4.4.Interpolation between optically thin and optically thick expressions yield a more general expression (Sirko & Goodman 2003) The radiation pressure is accordingly generalized as 2.2.Embedded stars as auxiliary heating sources The outer regions of marginally (un)stable (Q = 1) disks are prone to fragmentation and star formation.To provide thermal feedback that re-establishes energy equilibrium, stars convert some portion of the gas mass into heating, through CNO during the main sequence or triple-α reactions during the post-main sequence evolution.In contrast to star forming regions in galaxies, the emerging stars in AGN disks continue to accrete gas on both the IMS and SEPAD tracks.As these embedded stars gain mass, they also release intense winds.In most regions of AGN disks, their growth is suppressed by an accretion-wind balance when they attain equilibrium masses M ⋆ (≃ 10 3 M ⊙ ) with luminosities L ⋆ (≃ 10 7 L ⊙ ) comparable to their own Eddington limit L Edd,⋆ (Cantiello et al. 2021;Jermyn et al. 2022).Assuming one characteristic stellar mass M ⋆ and luminosity L ⋆ , the heating term in an equilibrium state is fundamentally given by where s ⋆ is the surface number density of stars.An energy equilibrium is reached with However, depending on the evolution tracks of stars and details of dynamical interactions, a fixed amount of heating rate corresponds to different reduction rates of gas supply from the disk accretion flow.The disk surface density evolves with a sink term The net removal rate Σnet is linked to The conversion of rest-mass energy into radiation through nuclear burning is given by Σenergy = s ⋆ L ⋆ /c 2 .In an accretion-wind equilibrium (Cantiello et al. 2021), the disk-star mass exchange on the IMS and SEPAD tracks do not contribute to Σnet .But, the formation of Chen & Lin black-hole remnants at the end of stars' lifespan τ ⋆ effectively permanently removes each remnants mass (M rem ) from the disk-gas reservoir at a rate Σrem = s ⋆ M rem /τ ⋆ . (15) The effective channels of gas return to the disk for two stellar evolution scenarios are shown in the schematics of Figure 1.
Figure 1.Differences between the IMS and the SEPAD stellar evolution models.For immortal stars, all gas deposited into stars is eventually recycled out, save for the fraction consumed to generate radiative energy.For stars that evolve off the main sequence within an evolution timescale of τ⋆, there is an additional removal rate of Mrem/τ⋆.
At the end of the SEPAD track, the magnitude of M rem may be a few M ⊙ (≪ M ⋆ ) and the remaining mass M ⋆ − M rem is returned to the disk flow, while a new generation of stars with mass M ⋆ forms to fill the energy deficit on the dynamical timescale ∼ Ω −1 ≪ τ ⋆ , maintaining an equilibrium stellar surface density.As black-hole remnants accumulate in number and accrete mass at the Eddington-limited rate, they deplete mass from the gas at a rate per unit area where τ Sal is the Salpeter timescale and Σ BH (t) is the surface mass density of remnant black-holes.Since τ Sal is generally much longer than other timescales, contributions from ΣBH ramps slowly with t and becomes significant towards the advanced stages (at ∼ τ AGN ) of AGN evolution.On the other hand, in the IMS limit, neither Σrem nor ΣBH exists.

Channels of gas depletion
Neglecting the ΣBH term first in a time-independent calculation, Equation ( 14) reduces to where ϵ ⋆ ∼ 10 −3 is the efficiency factor defined from the stars' characteristic evolution timescale if evolved in isolation τ 0 ≡ ϵ ⋆ M ⋆ c 2 /L ⋆ , which is ≪ τ ⋆ for the IMS stars (due to stellar metabolism) and τ 0 ∼ τ ⋆ for the SEPAD stars.To first order, the influence of feedback on the accretion flow from different stellar evolution models can all be expressed in terms of one gas-recycle efficiency factor the net gas removal rate is small when gas is efficiently recycled and ϵ recyle = 1 as in the IMS case.For stars on the SEPAD track, ϵ recycle ∼ ϵ ⋆ M ⋆ /M rem .Note when M ⋆ /M rem ∼ 1 as the average property of field stars, ϵ recycle ∼ ϵ ⋆ reduces to the Thompson et al. (2005) model, but such a case might not be appropriate for the SEPAD stars of typically more than a hundred solar masses ( §1).With these prescriptions in a thermal equilibrium (s Besides these effects, coexisting stars also undergo close encounters, binary capture, orbital in-spiral and coalescence, which may bring complication in the effective recycle frequency.For massive stars reaching equilibrium mass, their coalescence is followed by shedding of one star's worth of mass to the disk during a short shedding timescale (Jermyn et al. 2022, also see §3.3).If this shedding timescale is much shorter than the typical merger timescale τ shed ≪ τ merge , the stellar properties maintain their equilibrium values.For the IMS, this means no net change in the expression for ϵ recycle since an immortal star remains immortal after merger.For SEPAD track, the mixing of composition between young and old stellar population may introduces subtle changes in the effective τ ⋆ which affects ϵ recycle .
On the other hand, if merger is very frequent and τ shed ≳ τ merge , AGN stars would undergo runaway growth from sequential mergers and result in supermassive star that eventually explode as pair instability supernovae or collapse into massive black holes, but we show below that this scenario is rather unlikely since merger efficiency is limited by the timescale of inspiral after stars are captured into binaries.In the next section, we discuss how to estimate the merger timescale based on given disk and star properties.As a large number of massive stars co-exist in the disk, they encounter each other through dynamical interactions.High stellar surface number density may result in a short merger timescale, τ merger within which most stars can go through one typical merger process.The stellar capture process is analogous to planetesimal coagulation (e.g.Ida & Lin 2004).In particular, without gas damping, the physical collisional frequency is based on the radius of stars R ⋆ enhanced by gravitational focusing (Binney & Tremaine 1987;Palmer et al. 1993) with a corresponding timescale of where Safronov (1972) number, ∆V is star's dispersion velocity amplitude ( §3.1.3).
For gas-poor planetary rings (Goldreich & Tremaine 1978;Stewart et al. 1984) and planetesimals (Palmer et al. 1993), the effective capture impact parameter ∼ (1+Θ) 0.5 R ⋆ is usually much smaller than the Hill radius of these objects.This is because if planetesimals start from a cold configuration, their dispersion (non-circular) velocity ∆V grows while their Θ ⋆ decreases on a diskstar-relaxation timescale (Palmer et al. 1993;Aarseth et al. 1993) On the collisional timescale (τ col ), stars dissipate their kinetic energy at a rate d∆V 2 /dt ∼ ∆V 2 /τ col , i.e. the collisional damping timescale τ damp ∼ τ col .Therefore, a dynamical equilibrium is established with τ col ∼ τ V+ or Θ ∼ 1.
In gas-rich disks, embedded stars are surrounded and fed by mini-disks.Simulations have shown that during their encounters, disk-embedded binaries can become bound after their separation reach of order R H (Wang et al. 2018;Li et al. 2022aLi et al. , 2023)), which means their capture-impact parameter is effectively ∼ R H .This is because the existence of gas torques tend to damp velocity dispersion towards a much lower value corresponding to Θ ≫ 1, such that the effective impact parameter grows towards R H for increasing gas density.However, it must be capped at R H when the velocity amplitude is consistently reduced down to the Hill velocity V H = ΩR H .To connect the gas-poor and gas-rich regimes, it is useful to introduce such that the generalized expression for capture rate and timescale becomes

The migration timescale
Due to their tidal interaction with the AGN disk, isolated low-mass stars are affected by type I migration torque (Goldreich & Tremaine 1980;Ward 1997) with well-established prescriptions for migration timescales with f I ∼ O(1) being a dimensionless efficiency factor determined by Σ and T c distribution as well as the companion mass ratio defined as q = M ⋆ /M  Baruteau & Masset 2013).Provided their R H ≲ H and sufficiently high α, embedded stars undergo random walk (either inward or outward) over a distance ∆R ≃ R(P/τ mig ) per orbital period (Wu et al. 2024), with a diffusion coefficient There is also the possibility of partial gap formation and transition from type I to II migration (Lin & Papaloizou 1986).For more general range of q, the migration timescale in a laminar disk can be expressed as (Kanagawa et al. 2018): The companion follows rapid type-I migration at K ≪ 1, while for large companion mass K ≫ 1 the migration pace is greatly reduced due to the reduction of gas density around the companion vicinity by gap-opening effect the (Lin & Papaloizou 1986).However, this formula no longer applies for very deep gaps and massive companions' migration can be stalled or even reversed (Chen et al. 2020b;Dempsey et al. 2021).
Nevertheless, the above discussions apply to the migration of isolated stars.The net torque and τ mig are modified by the interference between the wake of multiple coexisting stars (Goodman & Rafikov 2001).These processes imply that the population of stars/compact objects in AGNs is unlikely to be uniformly guided towards a migration trap (Bellovary et al. 2016) for merging over a timescale of τ mig .Instead, their merger rates are ultimately governed by their velocity dispersion at each radius.A more precise prescription to account for the effect of migration is to introduce radial advection terms in s ⋆ in a time-dependent calculation.In this paper, we merely provide τ mig to calculate the velocity dispersion ∆V because it is closely connected to the gas damping timescale.

Velocity dispersion
In the AGN context, stars are embedded in their natal disks.Related to the migration torque ( §3.1.2),the stars' mean eccentricity is also damped by their tidal interaction with their host stars on a timescale of (Artymowicz et al. 1993;Kley & Dirksen 2006;Li et al. 2019;Chen et al. 2021) When τ e is significantly shorter than τ damp , the extra effect of gas damping becomes more significant to additionally damp ∆V towards a value ≪ GM ⋆ /R ⋆ and therefore Θ ⋆ ≫ 1.In this regime, we can also write Equation ( 22) as A dynamic equilibrium is established with τ V+ ∼ τ e .The solution to this equation, ∆V e , satisfies This term can be either substituted into R eff or directly compared with V H to ascertain the regime that Equations ( 24) & (23) belong to.If ∆V e > V H , the effective gas damping is not yet sufficient to raise R eff to R H .However, when ∆V e < V H , the capture radius has already reached its upper limit R H and ∆V e in fact decouples from the capture timescale.

Inspiral of binary stars
The gas rich capture process with R eff ≫ R ⋆ is much more efficient than direct collision due to much larger cross section.However, the expense is that even after capture into bound binaries, the stars still need to go through an extra inspiral timescale τ ins towards merger, so the effective timescale of one entire merger process is Due to the technical challenges of demanding numerical resolution over large dynamical range (R H /R ⋆ ), past numerical simulations of embedded binary in disks around SMBH (or central stars) have been only resolved to a fraction of R H (Baruteau & Lin 2010;Dempsey et al. 2022;Li & Lai 2022;Li et al. 2021Li et al. , 2022b)).They show that the magnitude of τ ins is highly uncertain, depending on the circumbinary disk (CBD) properties as well as the eccentricity, inclination of the hierarchical binaries' orbit, softening parameter of the gravitational potential, and the sink-hole radius of the stars or their black hole remnants.
Ignoring uncertainties in the direction of binary movement, a statistically-averaged estimate of the magnitude of the inspiral timescale influenced by the CBD gravitational torque satisfies (Baruteau & Lin 2010;Dempsey et al. 2022;Li & Lai 2022;Li et al. 2021Li et al. , 2022b) where M B ∼ 2M star is the total mass of the binary and ṀB is the accretion flow in the CBD towards the center.To first order, we can assume ṀB is comparable to the sum of accretion rates onto each of the binary component 2 Ṁ⋆ .Moreover, for massive stars reaching wind-accretion equilibrium, the accretion rate is close to the mass loss rate tapped by a near-Eddington luminosity In terms of the Eddington factor λ ⋆ = L ⋆ /L Edd , this translates to an inspiral timescale of where λ ⋆ is the star's luminosity Eddington fraction.Note that this estimate only considers a set of isolated binaries in a laminar disk.It's possible for the inspiral to stall at R ⋆ < a b < R H over a long time (≳ τ cap ), making it favorable for additional-body interaction to either disrupt the binary pair (Hut & Bahcall 1983) or forge hierarchical systems (Portegies Zwart et al. 2004).Secular interaction, including the von Zeipel-Lidov-Kozai, evection, and eviction resonances between the binary and SMBH may also lead to eccentricity excitation to speed up the merger process (Bhaskar et al. 2022;Gautham Bhaskar et al. 2023).When the disk structure is perturbed by turbulence in the limit R H ≲ H, the spin orientation of the gas inside their R H may fluctuate on the turbulent-eddy's turn over time (∼ P ) (Chen & Lin 2023) and slow down the inspiral process.

Post-merger mass shedding
Following the increases in their mass to M⋆ (> M ⋆ ), the merger products either become unstable and collapse, or adjust to a new hydrostatic equilibrium on a dynamical timescale.In the latter case, the elevated pressure and temperature in the stellar nuclei also strongly enhances their nuclear burning rate to a luminosity L⋆ in excess to the Eddington limit associated with their newly acquired M⋆ .Since M ⋆ of the progenitors provides a balance between their accretion and mass loss rates, the merged byproducts with super-Eddington luminosity shed a fraction or all of the excess mass ( M⋆ − M ⋆ ) during the return course to the pre-merger equilibrium state.If this mass readjustment does not occur dynamically during the merger episode, we still expect it to proceed on a timescale of (Jermyn et al. 2022, see their Equation 15): which bears the same expression as τ ins due to the windaccretion equilibrium.The final fate of stars depends on the hierarchy between τ shed and the merger timescale.
We discuss this aspect for specific disk profiles in §4.3.

IMS Track
Neglecting the effect of merger on stellar lifetime and remnant masses, we provide some quantitative analysis on the distribution of AGN stars for a constant ϵ recycle .In particular, we can estimate the magnitude of s ⋆ in IMS or SEPAD models by determining the Q − profile.Along the IMS track with ϵ recycle = 1, hardly any gas density is depleted to generate enough heating for the disk, and the disk Ṁd ≈ Ṁ• .Q − can be determined in terms of M • , λ • , ϵ • , α ν , κ cgs (Eqn 8).In the limit that Q − /ϵ recycle c 2 is negligible, the disk model resembles that of Sirko & Goodman (2003).Specifically, Equations 4 and 5 give, where r pc = R/1pc and Σ = 86m Note that Equations 4 and 5 are derived under the assumption Q = 1.At small radii R ≲ R Q=1 ≃ 0.01pc, so a full disk solution should include a transition to conventional models (with Q − balanced by viscous heating (Shakura & Sunyaev 1973;Frank et al. 2002)).
In the radiation pressure dominated and optically thick limit, the cooling rate is From Equation 11we see this function correspond to a stellar surface density of To compare with these scalings, we also provide the numerical solutions to the disk equations taking into account i) small Σnet due to ϵ recycle = 1; ii) both optically thin and thick cooling and iii) both gas and radiation pressure.The top panels of Figure 2 show h and Σ as functions of r pc , for different SMBH masses m 8 = 0.01, 0.1, 1, which follow closely the power laws given by Eqns 4 and 5.We adopt the parameterization Ṁd = 0.22λ d m 8 ϵ −1 d M ⊙ yr −1 analogous to Eqn 1 to define a radially varying accretion efficiency εd = ϵ d /λ d , and supply an accretion rate corresponding to εd (R out ) = 0.1 from the outer boundary, a value supported by SMBH population models (Yu & Tremaine 2002;Marconi et al. 2004;Shankar et al. 2004Shankar et al. , 2009) ) based on optical and x-ray data (Soltan 1982;Fabian 1999;Elvis et al. 2002).We also adopt κ cgs = α ν = 1 as our fiducial parameters.We solve the disk profile from R out inwards and stop at the radius where viscous heating Q + ν = 3 Ṁd Ω 2 /8π becomes comparable to Q − , since within that radius, the disk should match onto a viscously heated solution with no more star formation.The outer boundary is set to be R out = 10pc.
In the middle panels of Figure 2, we show T c and radiation to gas pressure ratio Π(:= P rad /P gas ) with molecular weight µ = 0.6m p .At r pc > r Π=1 ≈ 0.1m  disk becomes radiation pressure dominated.The lower panels of Figure 2 show s ⋆ and disk accretion rates as functions of r pc , assuming characteristic L ⋆ = 10 7 L ⊙ at their Eddington limits.At small radius the disk may be gas pressure dominated and T c is a constant, meaning At larger radius where radiation pressure is dominant, s * yields to Eqn 38 when optically thick.Finally, the Ṁd profiles in the lower right panel serve as a confirmation that due to large ϵ recycle the reduction in gas flow is negligible and Ṁd = Ṁ• is a valid approximation.

SEPAD Track
In the other limit, disk depletion becomes significant due to ϵ recycle ≪ 1.To estimate typical stellar density profiles under this context, instead of approximating with constant Ṁd , one must close the Equation set with a variable accretion rate (Eqn 19), akin to the procedures of Thompson et al. (2005).For a fixed ϵ recycle one has the following solution for Ṁd : where γ is an exponent associated with detailed opacity and equation of state scalings.The depletion fraction is where f (γ) is an order-unity factor.In the optically thick, radiation pressure dominated scenario, γ = 3/2 and f (3/2) = 8/3.Regardless of the exact value of γ, Eqn 41 represents the fraction of mass that needs to be taken from the disk accretion flow, i.e. subtracted off the accretion rate Ṁd , to support the total surface cooling rate ∼ πR 2 Q − .Another equivalent interpretation of f deplete is the ratio between star forming timescale Σ/ Σ⋆ and the accretion (disk-clearing) timescale πR 2 Σ/ Ṁd (Thompson et al. 2005).The f deplete = 1 condition gives the critical supply, Ṁcrit = f (γ)πQ − (R out )R 2 out /ϵ recycle c 2 , that is need from the outer boundary for non-negligible fraction of Ṁd (R out ) to reach the inner region.
At R ≪ R out , the accretion rate may relax to a constant ∼ Ṁd (R out )(1 − f deplete ) γ for f deplete < 1.Within this radius, the stellar density follows The decrease in mass accretion leads to an asymptotic profile that's reduced from Eqn 38 by a factor of ∼ (1−f deplete ) γ .However, it is worth noting that due to helium enrichment in the stellar core along the SEPAD track, their time-averaged L ⋆ is lowered by a factor of f L⋆ compared to immortal models which would increase s * ,SEPAD accordingly.Figure 3 shows example numerical SEPAD disk models for different SMBH masses as well as two groups of accretion rates, adopting an average L ⋆ = 2.2 × 10 6 L ⊙ informed by stellar evolution models from Ali-Dib & Lin (2023), corresponding to f L⋆ ∼ 4. For disk profiles corresponding to the dotted lines, We still fix εd (R out ) = 0.1 as in Figure 2. We assume M rem ∼ 30M ⊙ (Ali-Dib & Lin 2023) so ϵ recycle = 0.01.The recycle efficiency would be larger and f deplete smaller with a lower residue black hole mass at the end of massive star evolution M rem (Eqn 18).Consistent with Eqn 41, f d decreases with SMBH mass for a given ϵ ⋆ and εd .For the M • = 10 8 M ⊙ model, f deplete ∼ 1 so a small fraction of Ṁd (R out ) can flow into the inner parsec.For lower SMBH mass models, only the outer region is left with non-negligible surface density to contribute to star formation Σ * and radiative cooling term Q − , as well as small amount of AGN stars.The disk is markedly truncated within R trunc , which satisfies 1 − R trunc /R out ∼ 1/f deplete , when f deplete is larger than unity.
Nevertheless, there is no reason for εd at the outer boundary to be subject to rigid constraints since accretion feedback occur in the inner region of AGN disks.With somewhat larger Ṁd (R out ), f deplete would be markedly reduced and a much larger accretion rate fraction would reach the inner regions of the disk with Ṁd (R in ) ∼ Ṁ• .Solid lines in Figure 3 correspond to disk profiles that iteratively select M d (R out ) parameters, such that the inner region asymptotic Ṁd is consistent with εd (R ≪ R out ) = 0.1.For m 8 = 1, Ṁd (R out ) needs to be at least raised above the critical value by a factor of 3 such that f deplete < 1.For m 8 = 0.1, 0.01 the rise in Ṁd (R out ) need to be more significant.Due to the increase in gas supply, these disks have gas and massive star density even higher than the corresponding IMS models in the outer disk, while the profiles converge at small radii save for a factor of f L⋆ .

Quantitative Velocity Dispersions and Merger Timescales
The equilibrium number density of stars greatly influence their velocity dispersion and capture timescales.Informed by discussion in §3.1.1,we plot ∆V e /V H in Figure 4 for the IMS (left panel, disk profiles corresponding to Figure 2) as well as the SEPAD (right panel, disk profiles corresponding to Figure 3) disk models.In all cases, the effect of gas eccentricity damping is important in reducing the velocity dispersion amplitude from the escape velocity on the stellar surface (in the gas-free scenario) towards values close to the Hill velocity, especially so in the limit of very low s * , corresponding to the inner regions of low accretion rate SEPAD models.When ∆V e /V H < 1, the effective capture radius reaches its physical upper limit at R H .
In Figure 5 we plot the capture timescales (solid lines) calculated from the forementioned velocity dispersions, against the estimate of the inspiral/shedding timescales (horizontal red dashed lines).In the IMS models (top panel), since stellar number density is high, the capture timescale τ cap ∝ 1/s ⋆ is short, and for M • ≳ 10 7 M ⊙ , τ ins dominates the merger timescale τ merge = τ cap + τ ins for certain regions where s * is large.In this limit, hierarchical clusters, with a few members, may form.Since τ ins ∼ τ shed , mergers products can usually efficiently shed mass in between mergers and such that the characteristic stellar mass is marginally stable against runaway growth.Nevertheless, dispersion in τ ins and dynamical evolution in the hierarchical clusters may still result in a fraction of stars going through runaway mergers and eventually acquiring sufficient mass to undergo collapse.The byproducts of collapse are either black holes with minimal gas return (ϵ recycle ≪ 1) or pair instabil- Figure 3.The aspect ratio h, gas surface density Σ, temperature T , radiation to gas pressure P rad /Pgas, equilibrium stellar surface number density s * ,SEPAD = Q − /L⋆, as well as the variable disk accretion rate as functions of radius rpc for stars on the SEPAD track.We assume the radiative cooling is balanced by combined luminosity from massive stars with characteristic mass M * and luminosity L * .We assume κcgs = αν = 1.The solid lines correspond to models with inner boundary accretion rate matched to εd = 0.1, while dotted lines correspond to outer boundary accretion rate matched to εd = 0.1.
ity supernovae with nearly all mass returning to the disk (ϵ recycle = 1).
In the SEPAD models (right panel), if εd (R out ) = 0.1 is constrained (dotted lines), the capture timescale becomes much longer due to the lower equilibrium s ⋆ and therefore τ cap would always dominate the merger timescale, being the limiting timestep (τ cap ≫ τ shed , τ inspiral ).In this case, the typical stellar mass is firmly stable against runaway growth.On the other hand, if the outer boundary accretion rate can freely increase to match εd = 0.1 at the inner boundaries (solid lines), the radial dependence for τ cap would be similar to the IMS model.A lower average luminosity results in a longer shedding timescale as well as larger stellar number density, so the inspiral timescale again becomes the limiting timestep in a merger event.

Mass Budgets and Potential Effect of Opacity Window
In Figure 6 we plot the integrated stellar number N ⋆ (> R) for both the IMS (solid lines) and SEPAD (dashed and dotted lines) models, the asymptotic value of which at the inner boundary can reach of order ∼ 10 5 m 8 .For such a high number density of AGN stars, it is also important to examine the mass budget.The lower panel shows N ⋆ (> R)M ⋆ / Ṁd (R out ), providing an estimate of the timescale required for the accretion rate to be converted into this quantity of stars and achieve a quasi-steady state.For the IMS models (solid lines), the formation timescale is of order ≲ 10 7 years, a significant fraction of AGN lifetime (∼ 10 8 years), suggesting a steady state can only be achieved at mature stages of AGN evolution.The formation timescale is reduced for SEPAD models (dashed lines) where the outer accretion rate can increase to match εd = 0.1 at the inner boundary, since the accretion rate at R out is significantly enhanced while the stellar number density is increased by a smaller factor.In fact, it can be seen from Eqn 38 that s ⋆ only scales with ε−1/3 d while the mass accretion rate by definition scales with ε−1 d , such that a high accretion rate shortens the formation timescale ∝ s ⋆ / Ṁd and mitigate the mass budget issue.In fact, if the accretion supply from the outer boundary diminishes over an AGN lifetime, the large εd models may eventually transition to low εd SEPAD models (dotted lines), as the disk gradually hollows out from the inside.With the formation timescale is further prolonged due to the reduced accretion rate, we may expect this depletion to escalate into a runaway process during late stages of AGN evolution.Given that this formation timescale depends on the rate of mass supply, it is worth mentioning that realistic opacities can play a significant role in determining the upper limit of Ṁd (R out ) for a fixed SMBH accretion rate.While we confirm that that applying a more realistic grain opacity treatment κ ∝ T 2 below 100 K (Bell Chen & Lin & Lin 1994) only quantitatively affect our results, s ⋆ becomes sensitive to the opacity profiles around the sublimation temperature around T sub ∼ 1000 − 2000K.In this opacity window, if the midplane temperature T c is used to calculate κ(ρ c , T c ) and in turn the optical depth, radiative cooling becomes very efficient and the star formation rate must be exceptionally high to maintain energy balance.Therefore, even if f deplete < 1 (where f deplete is determined by outer boundary opacities) and a large amount of accretion rate can flow into the inner disk, only a fixed amount of Ṁin ∼ 0.1 − 1M ⊙ /yr, independent of Ṁd (R out ), can penetrate further beyond the opacity window, and the rest would all be converted into an additional population of stars in a small annulus close to the sublimation radius (Thompson et al. 2005, see their Appendix A).If such a steady state can be sustained, the stellar density in the outer disk could scale up with arbitrarily large Ṁd (R out ), while the timescale towards achieving steady state in the outer disk would be shortened due to the aforementioned reason.
However, due to the extreme sensitivity of opacity to temperature in this regime, the vertical structure of the disk plays a crucial role in determining the vertically integrated optical depth, and the average opacity that controls cooling may significantly deviate from κ(ρ c , T c ).For example, the disk photosphere can host an opacity significantly larger than the midplane while having nonnegligible density.Such a thermal stratification would lead to convective instability and an adiabatic vertical structure (Lin & Papaloizou 1980).Moreover, the temperature profile in the vertical direction may be heavily influenced by irradiation (Chiang & Goldreich 1997;Garaud & Lin 2007) and accretion of stars when s ⋆ becomes dense enough to cover the midplane.While the details of this feedback process is subject to investigation through radiation hydrodynamic simulation, we chose not to quantitatively model the opacity window in this paper given this major uncertainty.We only note that the overall stellar number density is already quite large in our estimates without the effect of opacity gap.Should the opacity window mechanism prove effective, our numbers for N ⋆ would only be a conservative estimate for moderate values of Ṁd (R out ) and can further increase for larger outer disk accretion rates.

Integrated Merger Rates
Taking τ merge = τ cap + τ ins , the integrated merger rate beyond a certain radius Γ ⋆ (> R) is given to be In Figure 7 we plot Γ ⋆ (> R) for both the IMS (solid lines) and SEPAD (dashed and dotted lines) models.The asymptotic values for Γ ⋆ (> R) at small radii indicate that the choice of stellar evolution models can have large impact on the distribution of merger rates.For the IMS models, Γ ⋆ (> R) strictly increases with SMBH mass for given α ν and εd , with differential contribution from all disk radii.For the SEPAD models, if εd = 0.1 is constrained at the outer boundaries (dotted lines), all contributions to Γ ⋆ (> R) would be in the outer disk while the asymptotic values are similar to the IMS model.If εd = 0.1 is constrained at the inner boundaries (dashed lines), Γ ⋆ (> R) can be notably larger than the IMS models due to larger s ⋆ at regions close to R out .

SUMMARY AND DISCUSSIONS
In this paper, we have derived equilibrium surface density profiles of massive stars in AGNs under two limiting contexts.In the first scenario the stars are immortal (IMS track) during AGN's active phase and efficiently recycle back all accreted gas if it's not converted into luminosity, such that the disk exhibits nearly smooth disk accretion-rate profiles.Alternatively, in the case where stars evolve off the main sequence into remnants (SEPAD track), their demise effectively removes significant amount of gas from the disk accretion flow and selflimit their number density.We define a physically motivated ϵ recycle factor to differentiate the effective removal of disk gas by stellar formation and evolution in the IMS and SEPAD scenarios (Eqn 18), which also serves as a natural bridge between Sirko & Goodman (2003)like and Thompson et al. (2005)-like disk models.This framework can be incorporated into state-of-the-art 1D disk modeling tools such as those constructed by Gangardt et al. (2024).
The observed populations of young stars near the Galactic center where (mostly B-type) S-stars with semimajor axes a ⋆ ∼ 5 × 10 −3 − 0.04pc (Genzel et al. 2003;Ghez et al. 2003;Gillessen et al. 2017) are surrounded by a population of clockwise-rotating disk (mostly O/WR type) stars with a ⋆ ∼ 0.04 − 0.5 pc, (Levin & Beloborodov 2003;Lu et al. 2006;Yelda et al. 2014;von Fellenberg et al. 2022).A zone of avoidance highlights the absence of stars with pericenter distance and eccentricity smaller than those of the S stars (Burkert et al. 2024).This observed feature is more consistent with the possibility that these young stars formed coevally in SEPAD-type disk models, with Ṁd increasing with R (Zheng et al. in prep).Additionally, we comment that a "hollow" disk due to star formation depletion can also possibly lead to low-luminosity AGNs such as NGC 3167 (Bianchi et al. 2019).
Focusing on constraining number density and merger rate of stars, we apply a simple constant-opacity disk model that captures the essential features of the stellarluminosity-heated outer region of AGNs.Integrating the stellar surface density, we estimate that the total number of AGN stars within a few parsecs to be N ⋆ ∼ O(10 3 ) − O(10 5 ).These objects can encounter each other and merge on a timescale of τ merge , which ranges from O(10 4 )−O(10 6 ) years, while the lower limit is set by the shortest timescale for binaries to harden under the effect of surrounding gas.We also estimate that for M • between 10 6 − 10 8 M ⊙ , the overall merger rates per AGN ranges from Γ ⋆ ∼ 10 −3 to ∼ 1 per year.This inference places some constraints in the rate of massive star mergers in large-scale surveys with ∼ 10 6 AGN samples such as ZTF (Frederick et al. 2021) & LSST (Ivezić et al. 2019) to be ∼ O(10 3 − 10 6 ) per year.Typical timescale is determined by details of the energy dissipation of the collision process (Hu & Loeb 2024), which is currently poorly understood.It may range from as short as the dynamical timescale to as long as the shedding timescale.Our calculation provides motivation for detailed hydrodynamic studies of massive mergers in AGN disks, which could help find traits to identify them from the plethora of transients.When mutual inclination between a massive star and the accretion disk is excited during the inspiral, the star can also collide and shock-heat the disk during each vertical passage and produce Quasi-periodic Eruption(QPE)-like signals (Linial & Metzger 2023;Tagawa & Haiman 2023).
An important characteristic of SEPAD stars is their contribution to the enrichment of the AGN environment.As they transition off the main sequence and evolve into stellar mass black holes (sBHs), each SEPAD star redistributes approximately ∆M Z ∼ a few so-lar masses of α elements back into the disk.Consequently, the overall pollution rate can be estimated as ṀZ ∼ πR 2 ∆M Z s ⋆,SEPAD /τ ⋆ .This production rate of heavy elements results in a significant increase in metallicity, with ∆Z ∼ ṀZ / Ṁd , which is several times higher than the solar value (Ali-Dib & Lin 2023) and consistent with observational data (Huang et al. 2023).On a timescale τ AGN ∼ 10 8 year, the cumulative surface density of produced sBHs grows to s BH ∼ s ⋆,SEPAD τ AGN /τ ⋆ and they eventually becomes the dominating heating source.The Eddington-limited rate of accretion onto these embedded black holes leads to ΣBH ≪ Σrem in Eqn ( 14), significantly reducing the stellar contribution needed to maintain a thermal equilibrium at late stages of AGNs, shutting off star formation (Eqn 12).In contrast, IMS stars neither contribute to the metallicity enrichment nor evolve into sBHs.
It has been suggested that merging black holes detected by LIGO Virgo collaboration may find one another in AGN disks (McKernan et al. 2012(McKernan et al. , 2014;;Tagawa et al. 2020;Samsing et al. 2022), with their merging process possibly associated with electromagnetic (EM) counterparts (Graham et al. 2020;Tagawa et al. 2023).Since disk capture of compact objects from the nuclear star cluster is inefficient due to the small cross section (Wang et al. 2018(Wang et al. , 2022)), the SEPAD track provides a promising channel for the formation of these gravitational wave sources.Chen et al. (2020a) suggested that if sBH mergers occur inside a gaseous medium, the gravitational wave signal can be altered by gas dynamics, which provides another way of directly identifying merger in gas-rich channels apart from EM counterparts.Additionally, if sBHs and massive stars can coexist in large numbers within an AGN disk, massive stars may efficiently capture sBHs successively on a similar timescale as the stars' own τ cap (since the impact parameter is dominated by the star's Hill radius), which then merge inside the core of the massive star after common-envelope evolution.In this scenario, gas dynamics would distort the GW signal even more significantly due high density in the stellar core.This kind of merger scenario is first proposed to explain possible electromagnetic counterpart to event GW150914 (Loeb 2016;Dai et al. 2017), and expected to be rare for field stars due to their short lifetimes.However, such events can be prevalent in AGN disks.
In subsequent works on detailed time-dependent AGN disk models, details of stellar evolution should be taken into account.It is conceivable that the disk may undergo regulated cycles between the two limiting solutions, if the stellar number/heating distribution is allowed to advect radially under effect of migration (Gr- Chen & Lin ishin et al. 2023;Wu et al. 2024).The capture of stars from the surrounding nuclear cluster can also be added as a source term for the stars at small radii (Wang et al. 2023;Generozov & Perets 2023).In SEPAD disk model where stars evolve into sBHs, the effect of accretion from sBHs should be taken into account in the heating and density depletion terms (Gilbaum & Stone 2022;Tagawa et al. 2022;Chen et al. 2023a;Zhou et al. 2024).Moreover, our calculation of binary capture and inspiral are subject to numerous uncertainties.To resolve these un-certainties and provide more accurate estimation of stellar merger rates, detailed investigation of migration and capture processes in a dynamically crowded disk using hydrodynamic and/or N-body simulations is warranted.

Figure 2 .
Figure2.The aspect ratio h, surface density Σ, temperature T , radiation to gas pressure P rad /Pgas, equilibrium stellar surface number density s * ,IMS = Q − /L⋆, as well as the variable disk accretion rate as functions of radius rpc for immortal stars.We assume the radiative cooling is balanced by combined luminosity from massive stars with characteristic mass M * and luminosity L * .We assume κcgs = αν = 1, and fix εd := ϵ d /λ d = 0.1 at the outer boundary.

Figure 4 .Figure 5 .
Figure4.∆Ve/VH on the IMS track for disk profiles presented in §4.1 (top panel) and the SEPAD track §4.2 (lower panel).When ∆Ve/VH < 1, the effective capture radius in Eqn 24 reaches its physical upper limit at RH.

Figure 6 .
Figure6.Upper panel: the integrated number of stars beyond radius R as functions of R or r = R/pc, for the IMS (solid lines) as well as SEPAD models (dashed lines and dotted lines); Lower panel: the integrated number multiplied by the factor M⋆/ Ṁd (Rout), which approximates the timescale towards establishing an equilibrium solution.

Figure 7 .
Figure7.The integrated/cumulative merger rates beyond radius R as functions of R or r = R/pc (Eqn 43), for the IMS (solid lines) as well as SEPAD models (dashed lines and dotted lines).