Evolution of eccentric stellar disks around supermassive black holes: the complex disk disruption dynamics and the milliparsec stars

We study the 10 Myr evolution of parsec-scale stellar disks with initial masses of $M_{\mathrm{disk}} = 1.0$ - $7.5 \times 10^4 M_\odot$ and eccentricities $e_\mathrm{init}=0.1$-$0.9$ around supermassive black holes (SMBHs). Our disk models are embedded in a spherical background potential and have top-heavy single and binary star initial mass functions (IMF) with slopes of $0.25$-$1.7$. The systems are evolved with the N-body code $\texttt{BIFROST}$ including post-Newtonian (PN) equations of motion and simplified stellar evolution. All disks are unstable and evolve on Myr timescales towards similar eccentricity distributions peaking at $e_\star \sim 0.3$-$0.4$. Models with high $e_\mathrm{init}$ also develop a very eccentric $(e_\star\gtrsim0.9)$ stellar population. For higher disk masses $M_\mathrm{disk} \gtrsim3 \times10^4\;\mathrm{M_\odot}$, the disk disruption dynamics is more complex than the standard secular eccentric disk instability with opposite precession directions at different disk radii - a precession direction instability. We present an analytical model describing this behavior. A milliparsec population of $N\sim10$-$100$ stars forms around the SMBH in all models. For low $e_\mathrm{init}$ stars migrate inward while for $e_\mathrm{init}\gtrsim0.6$ stars are captured by the Hills mechanism. Without PN, after $6$ Myr the captured stars have a sub-thermal eccentricity distribution. We show that including PN effects prevents this thermalization by suppressing resonant relaxation effects and cannot be ignored. The number of tidally disrupted stars is similar or larger than the number of milliparsec stars. None of the simulated models can simultaneously reproduce the kinematic and stellar population properties of the Milky Way center clockwise disk and the S-cluster.


INTRODUCTION
Stellar disks with non-zero eccentricities around supermassive black holes (SMBHs) appear to be a common feature of galactic nuclei in the local Universe.In the Local Group, both the nuclei of the Milky Way and M31 galaxies host stellar disks within the gravitational influence radii of their central SMBHs.For M31, the observed double nucleus (Light et al. 1974;Lauer et al. 1993;Kormendy & Bender 1999;Statler et al. 1999;Bacon et al. 2001;Sambhus & Sridhar 2002) is well explained by a apsidially aligned (i.e.lopsided) eccentric disk model (Tremaine 1995;Salow & Statler 2001, 2004).A schematic illustration of such a disk is provided in Fig. 1.The central SMBH of M31 appears to be surrounded by blue, bright stars (Lauer et al. 1998;Brown et al. 1998;Lauer et al. 2012) with an estimated age of 200 Myr (Bender et al. 2005).In addition, a number of complex observed features in the nuclei of local non-obscured, both dwarf and massive, early-type galaxies, such as offset or double nuclei and nuclei with central minima (Binggeli et al. 2000;Lauer et al. 2002Lauer et al. , 2005;;Houghton et al. 2006), can be explained by aligned ★ E-mail: anttiran@mpa-garching.mpg.deeccentric disks and viewing angle effects (e.g.Sridhar & Touma 1999;Madigan et al. 2018).
At the distance of the Milky Way center, ∼ 8.2 kpc (GRAVITY Collaboration et al. 2019), stellar kinematics can be studied in unprecedented level of detail compared to other galaxies.The central parsec of the Milky Way around the SMBH Sgr A* (e.g.Morris & Serabyn 1996;Alexander 2005;Genzel et al. 2010;Alexander 2011;Mapelli & Gualandris 2016) contains ∼ 100-200 young, bright and massive stars (e.g.Krabbe et al. 1991;Genzel et al. 1994;Blum et al. 1996;Lu et al. 2006) of mainly types O and WR.The ages of these stars must be very young, ≲ 10 Myr from life-times O-type stars.Paumard et al. (2006) gives an age estimate of 6 ± 2 Myr while Lu et al. (2013) prefers a somewhat younger age of 2.5-6 Myr.The total mass of the young stars can be as high as 1.4 × 10 4 M ⊙ -3.7 × 10 4 M ⊙ above 1 M ⊙ (Lu et al. 2013).The stellar initial mass function (IMF) in the Milky Way center is very top-heavy, the IMF ( () ∝  −  ) slope estimates ranging from  = 0.45±0.3(Bartko et al. 2010) to  = 1.7 ± 0.2 (Lu et al. 2013) in various structures in the central 0.5 pc.Resolved stellar kinematics (Genzel et al. 2000;Ghez et al. 2000;Levin & Beloborodov 2003;Beloborodov et al. 2006;Lu et al. 2006;Gillessen et al. 2009; Bartko et al. 2010) has revealed multiple disky and filamentary structures around the central SMBH Sgr A*, and most of the bright stars (up to ∼ 75%, von Fellenberg et al. 2022) reside in these disk structures.The main structures are the two disks at very large angle with respect to each other, the clock-wise (CW) and the counter-clockwise disk (CCW) (e.g.Genzel et al. 2003;Tanner et al. 2006;Paumard et al. 2006;von Fellenberg et al. 2022).The CW disk is significantly warped (Bartko et al. 2009, see also Yelda et al. 2014) and is moderately eccentric with a median eccentricity of 0.4-0.5.The CCW disk has similar median eccentricity as the CW disk but also contains higheccentricity stars, though not beyond  > 0.9 (von Fellenberg et al. 2022).Additionally, two outer disky or filamentary structures have been recently reported (von Fellenberg et al. 2022).The median eccentricity of these structures is higher than in the CW and CCW disks,  ∼ 0.7.
For the formation and evolution of such nuclear disks, hydrodynamical simulations have been used to study fragmentation and star formation both in circular (Nayakshin et al. 2007) and eccentric (Alexander et al. 2008) gas disks and rings around SMBHs.The molecular cloud infall and disruption simulations (Sanders 1998;Bonnell & Rice 2008;Mapelli et al. 2008Mapelli et al. , 2012;;Lucas et al. 2013;Trani et al. 2018) have shown to result in clumpy, filamentary, eccentric and rapidly star-forming disks with a top-heavy IMF.Various other dissipative scenarios have been investigated, including gas clouds engulfing SMBHs and forming accretion disks (Wardle & Yusef-Zadeh 2008;Alig et al. 2011), collisions of clouds with circumnuclear gas rings (Alig et al. 2013) or clouds colliding with other clouds (Hobbs & Nayakshin 2009).
In this work we specifically focus on the dynamics of lowmass asymmetric eccentric disk systems for which  disk ≪  • .The stability of such disks around SMBHs depends on how massive the disk is compared to the mass of the background stellar cusp at the spatial scale of the disk.Low-mass disks with  disk ≪  cusp are unstable due to the so-called secular eccentric disk instability (Madigan et al. 2009(Madigan et al. , 2011, see also Haas & Šubr 2016;Šubr & Haas 2016).The present-day individual Milky Way center disks around the  • = 4 × 10 6 M ⊙ SMBH (GRAVITY Collaboration et al. 2019, 2020;Event Horizon Telescope Collaboration et al. 2022) are in this range as  disk ∼ 10 4 M ⊙ for the CW disk and  disk ∼ 5 × 10 3 M ⊙ for the CCW disk (Bartko et al. 2010) and the background cusp mass within  0 = 1 pc is  cusp ( 0 ) ∼ 5 × 10 5 M ⊙ (e.g.Genzel et al. 2003;Schödel et al. 2007).These Milky Way values are collected into Table 1.On the other hand, massive eccentric disks ( disk ≫  cusp ) such as the asymmetric M31 center structure can be stable over gigayear timescales (Tremaine 1995;Jacobs & Sellwood 2001;Madigan et al. 2018).
A long-time motivation for studying the stellar dynamics in the central parsec of the Milky Way has been the goal of understanding its very central arc-second (∼ 0.04 pc).Going towards the central SMBH from the inner edge of the CW disk, the properties of the orbiting stars dramatically change.The innermost O/WR stars have a semi-major axis close to  ★ ∼ 0.04 pc inside which most of the observed stars are of type B. This compact cluster around the Milky Way SMBH is the famous S-cluster (Eckart & Genzel 1996;Ghez et al. 1998).Rather than residing in a disky structure, the orbits of the S-stars are isotropic.The eccentricities of the S-stars are predominantly high compared to the outer disk stars, following a thermal or even superthermal distribution (e.g.Schödel et al. 2002;Ghez et al. 2003;Eisenhauer et al. 2005;Ghez et al. 2005;Gillessen et al. 2008Gillessen et al. , 2009Gillessen et al. , 2017)).Just as the O/WR stars in the disks, the B-stars of the S-cluster must be young, having an age below ≲ 400 Myr from their main sequence lifetimes.Whether the S-cluster stars and the disk originate from a common formation scenario still remains an open question.The S-cluster lacks O/WR stars, which are prominently present in the disks, and being isotropic and very eccentric the orbits of S-cluster stars considerably differ from the orbits of the disk stars (Gillessen et al. 2017;von Fellenberg et al. 2022).
One of the scenarios for the formation of the S-cluster relies on the Hills mechanism.Given a small enough pericenter distance  p , SMBHs can disrupt binary stars (Hills 1988).One binary component is captured by the SMBH on a very eccentric orbit whereas the other component becomes unbound is ejected as a high-velocity star (Bromley et al. 2012;Kobayashi et al. 2012;Zhang et al. 2013; Chen & Amaro-Seoane 2014; Madigan et al. 2014;Zhang et al. 2013;Generozov & Madigan 2020;Generozov 2021).An A-type high-velocity star with a velocity of 1800 km/s and trajectory consistent with ejection by the Milky Way SMBH ∼ 5 Myr ago has recently been discovered (Koposov et al. 2020), which might indicate that the Hills mechanism indeed operates in the Milky Way center.However, a single high-velocity star still does not prove that the Hills mechanism is the dominant formation pathway for the S-stars.Various competing scenarios to the Hills picture involve a low-eccentricity origin for the S-stars (Genzel et al. 2010).The low-eccentricity scenarios include disk migration in the initial (gaseous) disk (Levin 2007) and migration aided by a massive perturber (Perets et al. 2009), such as an intermediate-mass black hole (e.g.Mastrobuono-Battisti et al. 2014).It should be noted that both the high and the low-eccentricity origin channels require an efficient relaxation mechanism to change their initial eccentricity distribution to the observed one, either from above or below (Genzel et al. 2010).
Various N-body simulation studies have been devoted to studying the origin and dynamics of the S-cluster stars (e.g.Berukoff & Hansen 2006;Kenyon et al. 2008;Perets et al. 2009;Madigan et al. 2009;Fujii et al. 2010;Perets & Gualandris 2010;Antonini & Merritt 2013;Hamers et al. 2014;Šubr & Haas 2016;Generozov & Madigan 2020).The simulations can be classified in two broad categories: the S-star origin simulations and S-cluster relaxation simulations.There is substantial overlap between the two categories, but only few simulation studies have explored the full dynamical picture using a single simulation code in a self-consistent manner.For the disk binary origin scenario we highlight the simulations of Šubr & Haas (2016) as the one of the most self-consistent setups to date.Their simulations begin from an eccentric disk with an IMF and a binary star population, follow the disk disruption and Hills mechanism of binary disruption, and finally the relaxation of the captured stars using the NBODY6 code version of Aarseth (2003).However, their full setup did not include a stellar background population, post-Newtonian effects or stellar evolution.Thus, the eccentric disk instability and Hills mechanism scenario for the S-cluster formation has not yet been fully explored in the literature in a self-consistent manner.
In this study we examine the evolution of initially asymmetric eccentric disks and the dynamical consequences of the eccentric disk instability using direct N-body simulations.Our simulations have a number of important differences compared to the previous studies regarding the simulation setup and the numerical code used.Our initial conditions include an asymmetric eccentric disk model (with single and binary stars) around a SMBH embedded in a smooth external background potential representing the background stellar cusp.The gravitational dynamics of the system from disk disruption to Hills mechanism and the evolution of the population of captured stars around the SMBH is integrated using a single simulation code BIFROST (Rantala et al. 2021(Rantala et al. , 2023)).The code combines a hierarchical version of the fourth-order forward symplectic integrator (Chin 1997;Chin & Chen 2005;Chin 2007;Dehnen & Hernandez 2017) with specialized secular and regularized solvers for binaries, close fly-bys, multiple systems and small clusters around massive black holes (Rantala et al. 2020).As described above, most simulation studies have concentrated in separate aspects such as the eccentric disk instability, Hills mechanism or the relaxation of the S-cluster, or included several but have switched from one simulation code to another during the study.We also employ single stellar evolution, tidal disruption event (TDE) prescription and post-Newtonian equations of motion for the stars in our simulation setups.
Next, we study the dynamics of the asymmetric eccentric disk instability in the intermediate-mass regime, i.e. the disk is still less massive than the background cusp, but its mass is not negligible (with still assuming  disk ≪  • ).Thus, the inequality  disk <  cusp still holds, but the relation  disk ≪  cusp does not.Because of this, the disk disruption process is more complex than the standard secular eccentric disk instability scenario of Madigan et al. (2009), as the assumption of the negligible contribution of the disk mass itself to the precession of the disk star orbits does not hold anymore.We argue that the Milky Way center CW and CCW disks might have been in this mass regime at their formation time ∼ 6 Myr ago, and being top-heavy they have since lost mass due to stellar evolution, reaching  disk ∼ 10 4 M ⊙ at present day.Simulation studies in the literature have focused on studying the disruption of asymmetric eccentric disks with a mass of the present-day CW disk (e.g.Madigan et al. 2009;Šubr & Haas 2016) and not considered the fact that the disks were more massive at their formation time in the past.
The article is structured as follows.After the introduction we briefly review the relevant stellar-dynamical processes around SMBHs for this study in Section 2. The updated simulation code and the initial conditions for numerical simulations are described in Section 3 and Section 4. Next, we discuss the simulation results.The evolution of the intermediate-mass asymmetric eccentric disks, their disruption and later evolution are examined in Section 5 while Section 6 presents the properties of high-velocity stars and stars accreted by the SMBH in the simulations.The origin and the properties of milliparsec-scale stellar population are investigated in Section 7. We discuss our results in the context of Milky Way center observations and previous simulation studies in 8. Finally, we summarize and conclude in Section 9. component present-day mass notes cusp ∼ 5 × 10 5 M ⊙ old stars within 1 pc Table 1.The observed present-day masses of the Milky Way SMBH (GRAV-ITY Collaboration et al. 2019), the clockwise and the counterclockwise disk (Bartko et al. 2010) and the old spherical cusp population (Genzel et al. 2003;Schödel et al. 2007).
physical effect timescale at 0.01 pc at 0.  The orbital eccentricity is assumed to be  ★ = 0.5.The Lense-Thirring timescale is calculated for a rapidly spinning black hole with  = 1.The relaxation timescales  NR ,  SRR and  VRR are adopted from the recent estimates in the literature (Panamarev & Kocsis 2022).Lower estimates for the scalar resonant relaxation timescale have also been proposed, down to 10 Myr near 0.01 pc (Antonini & Merritt 2013).

Relativistic precession
We very briefly review the most important stellar dynamics processes around supermassive black holes relevant for this study.If the SMBH is considerably more massive than the background stellar cusp ( • ≫  cusp ), the motion of stars around the SMBH is nearly Keplerian.The Keplerian orbital period of a star with a mass  ★ and a semi-major axis  ★ around a SMBH with mass  • is  2) while the reduced angular momentum vector ì ℎ is defined as ì ℎ = ì  × ì .At the apocenter the torque ì  due to a small non-Keplerian perturbation ì If the perturbing background is spherically symmetric (left panel), the torque and thus the precession direction is always retrograde.In the asymmetric case (right panel) both retrograde and prograde torques may occur at different points of the orbit, and the net torque integrated over the entire orbit determines the precession direction.
in which  =  • +  ★ .If the stellar cusp is spherically symmetric, the reduced angular momentum vector ì ℎ = ì  × ì  of an example star is constant and the orbit of the star is confined into a plane.The eccentricity vector (also known as the Laplace-Runge-Lenz vector) ì  of the star pointing towards the pericenter defined as is another constant of motion in the pure Keplerian case, and the direction of the periapsis of the orbit does not change over time.
The eccentricity of the stellar orbit is the norm of the eccentricity vector,  ★ = ∥ ì ∥.Both Newtonian (non-Keplerian) and relativistic effects can cause the orbit of the star to precess, i.e. the direction of the eccentricity vector will change with time.In the following we assume  • ≫  ★ .In the lowest-order post-Newtonian order PN1.0 the precession time-scale  GR of the orbit is in which  g is the gravitational radius of the SMBH.The PN1.0 precession is commonly referred to as the geodetic, mass, or de Sitter precession.The precession direction is famously prograde (Einstein 1915).The spin ì  • of the SMBH is the source of the frame dragging or the Lense-Thirring nodal precession (Lense & Thirring 1918).The precession time-scale (Poisson & Will 2014) is On the milliparsec scales we are interested in this study the Lense-Thirring precession is weak compared to the PN1.0 precession for moderately eccentric stellar orbits and would become important only on spatial scales close to the gravitational radius  g (Merritt 2013).

Mass precession -the spherical case
Precession effects are present in the Newtonian case as well as long as the potential deviates from the Keplerian point mass potential.Following e.g.Madigan et al. (2018), the rates of change of the reduced angular momentum vector ì ℎ and the eccentricity vector ì  of a star orbiting the SMBH are in which ì  is a small non-Keplerian contribution to the total acceleration from the extended mass distribution around the SMBH.This is schematically illustrated in Fig. 2. In the simple case of a spherically symmetric but extended mass distribution, the non-Keplerian perturbing acceleration will always point towards the SMBH due to Newton's shell theorem, so the change of the eccentricity vector will always be in the retrograde direction, ì  ′ • ì  p < 0.Here ì  p is the periapsis velocity.If the stellar orbit is embedded in a stellar cusp with mass  cusp ≪  • , the precession time-scale (Merritt et al. 2011;Merritt 2013) is The direction of this so-called mass precession is retrograde.For a spherical power-law cusp with a slope of  the mass precession time-scale scales as so for  = 3/2 the precession rate is independent of  ★ .It is very important to note that the mass precession timescale  M monotonically increases with increasing orbital eccentricity as  ()/ > 0 in Eq. ( 6).Thus, stars on more eccentric orbits precess at a lower rate than stars on more less eccentric orbits, and almost radial stellar orbits in the limit  ★ → 1 do not precess at all.Estimates for both the relativistic and Newtonian precession timescales for the Milky Way center are collected in Table 2.

Mass precession -the asymmetric disk case
Another relevant setup for our study is an asymmetric eccentric disk around a SMBH.Such a setup is illustrated in the right panel of Fig. 1.The mass precession process of a stellar orbit in an asymmetric eccentric disk is more complicated than in the simple spherical case.Naively, one might expect that the contribution of the disk mass within the orbit of the star would cause retrograde precession just as in the spherical case.However, the eccentric disk models are asymmetric and thus intuition from Newton's shell theorem does not hold.The perturbing non-Keplerian acceleration ì  as in Fig. 2 can point to different directions along the orbit.The net torque on a single stellar orbit depends on the semi-major axis and eccentricity distribution of the disk (i.e. its surface density profile).The change of the eccentricity vector direction over one orbit and  M can be in principle can be calculated using Eq. ( 5).
In the inner parts of the disk the resulting precession is prograde ( ì  ′ • ì  p > 0) as most of the mass lies outside the orbits in an asymmetric configuration, and retrograde in the outermost parts of the disk as most disk mass is located within the orbits.The exact details where the prograde rotation turns into retrograde depend on the disk properties.For example, Madigan et al. (2018) reports a disk model with up to 10% of the disk mass in the outer parts precessing retrograde while the majority of the disk stars precess in the prograde direction.

Asymmetric eccentric disks
Circular and axisymmetric eccentric disks can be stable over long periods of time, but this is necessarily not the case for asymmetric eccentric disks for which the eccentricity vectors are apsidially aligned.In the simple case of a SMBH, a spherically symmetric stellar background cusp and an asymmetric eccentric disk (with a narrow eccentricity distribution), the stability of the asymmetric disk depends on how it precesses.Disks massive compared to the background stellar population ( disk ≫  cusp ) precess prograde and are stable (Tremaine 1995;Jacobs & Sellwood 2001;Madigan et al. 2018) while low-mass disks ( disk ≪  cusp ) precess in the retrograde direction and rapidly disrupt (Madigan et al. 2009;Šubr & Haas 2016).

Stable massive disks
We show a schematic illustration explaining the stability of massive ( disk ≫  cusp ) asymmetric eccentric disks in the two bottom panels of Fig. 3.The prograde precession of the disk orbits is caused by the mass of the disk itself.The prograde precession originates from the fact that in most points of the orbit of the star the perturbation from the disk to the star points away from the SMBH, causing a prograde torque.A disk star with an eccentricity higher than the mean eccentricity of the majority of disk stars will precess at a slower prograde rate than the disk stars, and thus begins to lag behind the bulk of the disk.Now the disk perturbs the lagging star and torques its orbit, increasing the magnitude ∥ ì ℎ∥ of the angular momentum vector ì ℎ = ì  × ì  of the star and thus decreasing its eccentricity.Lower eccentricity leads to a faster prograde precession rate, and the star catches up the bulk of the disk and is reabsorbed into its main body.
The opposite process occurs for stars on initial orbits with lower eccentricity.They precess in the prograde direction faster than the bulk of the disk and thus begin to begin to lead the main disk.The torque from the main disk now works to decrease ∥ ì ℎ∥ and increase  ★ , lowering the prograde precession rate.The main body of the disk can thus catch up the leading stars, again leading to the re-absorption of the star to the main disk.Thus, massive eccentric disks are stable against this type of eccentricity oscillations.

Unstable low-mass disks
The unstable case of the low-mass ( disk ≪  cusp ) asymmetric eccentric disk is presented in the top panels of Fig. 3. Now the massive background population causes the orbits the disk stars to precess in the retrograde direction.Compared to the prograde case, the sense of stars leading or lagging behind the disk is reversed as leading and lagging behind are defined with respect to the precession direction.Thus, the retrograde lagging behind case corresponds to the prograde leading case, and the retrograde leading case the prograde lagging behind case.
A star with a slightly higher eccentricity than the most disk stars will precess at a slower rate than the disk in the retrograde direction, and thus will begin to lag behind the main part of the disk.Now, the torque from the bulk of the disk causes ∥ ì ℎ∥ to decrease, increasing the eccentricity of the star even further.This leads to an increasingly slow precession rate, and the star further recedes from the main part of the disk.
On the other hand, stars with lower eccentricity than in the bulk of the disk will precess at a faster retrograde rate and begin to lead the disk.The torque from the disk acts to increase ∥ ì ℎ∥ and decrease  ★ .This further increases the retrograde precession rate, again leading into a runway process.Thus, asymmetric eccentric disks in the regime ( disk ≪  cusp ) are unstable against eccentricity oscillations.The initially narrow eccentricity distribution is broadened while the disk is disrupted.The stars lagging behind the disk will reach high eccentricities while the stars leading the disk become more circular.The entire process is termed the secular eccentric disk instability following Madigan et al. (2009).

Binary stars, the Hills mechanism and tidal disruption of stars
SMBHs can unbind binary stars resulting in the ejection of the one binary component and the capture of the another component around the SMBH (Hills 1988;Yu & Tremaine 2003).This is known as the Hills mechanism.Only binary stars close enough to the SMBH can become unbound, the approximate limit being the tidal radius  t defined as in which  b and  b are the mass and semi-major axis of the binary.Assuming a parabolic orbit for the binary center-of-mass, the semimajor axis  s of the captured star around the SMBH is in which  ∼ 0.3-1 depends on the eccentricity and the phase of the binary (Generozov & Madigan 2020).The orbit of the captured star is very eccentric: with  k ∼ 1-2.For example, a  b = 20 M ⊙ binary with a semimajor axis of  b = 1 AU disrupted by the Milky Way SMBH with  • = 4 × 10 6 M ⊙ will result in a captured star with orbital elements of  s ∼ 17-55 mpc and  s ∼ 0.97-0.98.The tidal forces of SMBHs affect single stars as well.A single star can be fully or partially disrupted at the separations smaller than the order-of-magnitude tidal (or Roche) radius  t defined as in which  ★ is the radius of the star.Part of the disrupted stellar material can remain bound to the SMBH, and falls back into the SMBH in an accretion flow powering a luminous observable transient with a characteristic decay time-scale, a tidal disruption event (Hills 1975;Rees 1988;Gezari 2021).The modeling of tidal disruption events (TDEs) is in general challenging and numerically expensive, requiring advanced (magneto)hydrodynamical solvers, general relativity and detailed models stellar structure.In addition, the parameter space of the problem is large (including the stellar models and orbits A schematic illustration to understand the instability and and stability of asymmetric eccentric disks.Here all the precession processes have a Newtonian origin.The main body of the asymmetric eccentric disk is depicted here as a gray background with a point cloud perturbing a target star (yellow symbol), which is either leading the disk or lagging behind it.It is very important that the terms leading and lagging behind are defined with respect to the precession direction (curved arrow).The precession direction (with respect to the stars apocenter velocity) entirely defines whether the asymmetric eccentric disk is stable or unstable.If the precession is retrograde (top row) as in the case  cusp ≫  disk , the torques  = ì  × ì  from the disk will result in the star receding from its main body.Stars lagging behind the disk will increase their eccentricities while leading stars end up on more circular orbits.In the prograde case (bottom row) when  cusp ≪  disk the situation is exactly the opposite.Leading stars will circularize and stars lagging behind the disk main body will reach higher eccentricities.In both prograde cases the star ends up being reabsorbed to the main body of the disk, so the disk is stable despite the eccentricity and precession rate oscillations.

Relaxation processes
While the Hills mechanism produces stars on very eccentric orbits around the SMBH, the cumulative eccentricity distribution of the Sstars in the Milky Way center is  () =   with  ∼ 2.6 (Gillessen et al. 2009).Thus, an efficient mechanism would be required to relax the orbital eccentricities of the Hills mechanism origin Scluster stars withing the age of the S-cluster, ∼ 10 Myr (Habibi et al. 2017).
The common non-resonant two-body relaxation changes both the energy and angular momentum of the stars (and hence  ★ and  ★ ) in a random-walk like manner.The two-body relaxation timescale  NR of  stars around a SMBH can be estimated as in which m = ⟨ 2 ⟩/⟨⟩ (Binney & Tremaine 2008).Here  2 = ln ( • / ★ ) is a factor of the order of unity.The two-body relaxation time-scale  NR for the S-stars is considerably longer than the ages of the stars (≲ 10 Myr), so the effect of two-body relaxation on the S-stars is weak.
In near-Keplerian potentials the orbits of the stars are constrained by the symmetries of the potential, and torques on the stellar orbits are correlated and build up in a coherent manner (Rauch & Tremaine 1996).While the energy (semi-major axis) is relatively unaffected by resonant relaxation, both the magnitude (scalar resonant relaxation, SRR) or the direction (vector resonant relaxation, VRR) of the angular momentum vector can change.Thus, scalar resonant relaxation can relax the eccentricities of the S-cluster stars while vector resonant relaxation can change the orientation of their orbital planes.
The scalar resonant relaxation timescale  SRR is approximately when mass precession (not the relativistic precession) is the dominant orbit precession mechanism.It has been argued that the scalar resonant relaxation timescale can be as low as  SRR ∼ 10 Myr, comparable to the age of the S-cluster, especially when the effect of compact remnants from the old spherical stellar population is included in the calculations (e.g.Antonini & Merritt 2013;Bar-Or & Fouvry 2018;Generozov & Madigan 2020).
The timescale for vector resonant relaxation (VRR) in a spherical system around a SMBH (Eilon et al. 2009;Kocsis & Tremaine 2015;Fouvry et al. 2019) can be estimated as in which  () is the mass enclosed by the orbit of the star, and  v is again a constant of the order of unity.The  VRR in the Milky Way center is less than a Myr the VRR being the fastest relaxation process in the S-cluster region (Panamarev & Kocsis 2022).Estimates for the relaxation timescales in the Milky Way center are presented in Table 2.

Gravitational dynamics
For the numerical simulations of this study we use the novel N-body simulation code BIFROST (Rantala et al. 2021(Rantala et al. , 2023)).BIFROST is an accurate GPU-accelerated direct-summation N-body code based on the hierarchical (Rantala et al. 2021) fourth-order forward symplectic integrator technique (Chin 1997;Chin & Chen 2005;Chin 2007;Dehnen & Hernandez 2017).The main integrator of the code is supplemented with specialized secular and regularized few-body integrators for binary systems, close fly-bys, multiple systems and small clusters around supermassive black holes (Rantala et al. 2020).
The code includes optional post-Newtonian (PN) equations of motion for the simulation particles, the option for external background potential, simple prescriptions for stellar evolution, stellar mergers and tidal disruption events (TDEs), and supports arbitrary fraction of binary stars due to its effective MPI parallelization.If PN equations of motion are enabled, the interactions between the stars and the supermassive black hole are post-Newtonian everywhere in the simulation domain, not just in the subsystem regions as it is commonly done in the literature.We include only the lowest-order PN1.0 term and no spin PN terms as their dynamical effect on the simulation results would be small in the setups we are interested in.
The user-given accuracy parameters of BIFROST in this study are displayed in Table 3.In addition to the code options in Rantala et al. (2023), we introduce a number of additional code functionalities and parameters.Massive particles (such as SMBHs) can now optionally have larger subsystem radii  rgb,• than the subsystem radii of stars  rgb,★ .We also introduce a time-step limiter for small particles (i.e.stars) in the vicinity of a massive particle such that all particles within  ngb,• <  <  Δt,• acquire the minimum time-step of the particles within the region.Both these small updates increase the integration accuracy of stellar orbits around SMBHs.

Simple stellar evolution
The current version of the BIFROST code is not coupled to a binary stellar evolution module.In order to include a simple stellar evolution prescription into this study we use the SEVN rapid population synthesis code (Spera et al. 2015;Spera & Mapelli 2017;Spera et al. 2019;Mapelli et al. 2020;Iorio et al. 2023) to calculate the evolution of mass  ★ () and radius  ★ () of stars in our simulation.
We use the metallicity  =  ⊙ throughout the study.The SEVN code itself is based on interpolating pre-tabulated sets of latest PARSEC (Bressan et al. 2012;Chen et al. 2015;Nguyen et al. 2022) stellar tracks.
The main effect of stellar evolution in our study is the rapid mass loss and death of massive stars.First stars die at the age of ∼ 2.5 Myr.At the inferred age of the CW disk O-type stars in the Milky Way center,  = 6 Myr, stars more massive than ∼ 29 M ⊙ have already ended their lives.At the end of our simulations,  = 10 Myr, no stars more massive than  ★ ∼ 18 M ⊙ are alive anymore.Obviously, including or not including stellar evolution has a large effect on the number of O-type stars.Thus, incorporating even a very simple stellar evolution prescription into our simulations allows a better comparison with observations.Another motivation for including simple stellar evolution in the runs is the study of single stars in the S-cluster.Massive stars rapidly lose mass in a few Myr before the end of their lives, and entering the giant phase makes the stars extremely susceptible for tidal disruption by the SMBH due to the significantly increased stellar radius.The effect of binary stellar evolution on the results of this work will be addressed in a future study.

Tidal disruption events
We briefly describe here the updated the TDE prescription of BIFROST (Rizzuto et al. 2022;Rantala et al. 2023).First, we do not use the tidal capture and drag force routines used in the intermediate mass black hole growth study of (Rizzuto et al. 2022) in this study.The stars within the tidal disruption radius from the SMBH are accreted by the SMBH and removed from the simulation.For this study we update the definition of the tidal disruption radius R t using the fitting formulas from detailed hydrodynamical tidal disruption simulations of Ryu et al. (2020).The tidal disruption radius is defined as in which  ★ is the radius of the star,  t is the order-of-magnitude tidal disruption radius (Kochanek 1992) and the correction factors are given in Ryu et al. (2020).The SMBH term is of the order of unity (Ψ • ∼ 1.32) for the Milky Way SMBH mass of  • = 4 × 10 6 M ⊙ , and the star term Ψ ★ ( ★ ) ∼ 1/2 for stars considerably more massive than 1 M ⊙ .Thus, massive stars such as B-type and O-type stars we are most interested in this study disrupt closer to the SMBH compared to the standard order-of-magnitude tidal disruption radius estimate.In all TDEs the star is fully accreted by the SMBH, i.e. we do not include a partial tidal disruption prescription in our code for this study.Possible caveats of the adopted simple tidal disruption model are elaborated in Section 8.4.2.

Motivation
The dynamics and instability of asymmetric eccentric stellar disks of masses around  disk ∼ 10 4 M ⊙ have been widely studied in the literature (e.g.Madigan et al. 2009;Šubr & Haas 2016), as this disk mass is close the inferred mass of the CW and CCW disk stellar disks in the Milky Way center (Bartko et al. 2010).However, there are caveats in using simulations of  disk ∼ 10 4 M ⊙ stellar disks to study the structure and long-term evolution of the Milky Way center.The inferred age of the young, bright and massive stars (6 Myr) is enough that stars more massive than ∼ 29 M ⊙ have already reached the end of their life-times.Thus, a top-heavy  disk ∼ 10 4 M ⊙ stellar disk today was more massive at its formation.For extremely top-heavy power-law IMF slopes the mass loss of the stellar population in 6 Myr is substantial.With  = 0.25 and  max = 120 M ⊙ , a stellar population can lose 90% of its mass in only 6 Myr.For more moderately top heavy IMFs of  = 1.3 and  = 1.7 the mass loss in the same time period is ∼ 70% and ∼ 50%.This encourages studying the dynamics of initially more massive eccentric disk systems in the context of the Milky Way center.However, the more massive initial disks do not necessarily fulfill the  disk ≪  cusp criterion anymore for the secular eccentric disk instability (Madigan et al. 2009).As the eccentric disks of this mass scale do not fulfill the high-mass stability criterion  disk ≫  cusp either, we expect that the disks in this intermediate mass range are unstable, but their disruption mechanism might differ from the model of Madigan et al. (2009).Understanding the instability of these intermediate-mass asymmetric eccentric disks is the key motivation for this study.

Disk properties
We construct asymmetric eccentric disk models of stars orbiting a SMBH embedded in a static spherical background cusp potential.In this study the mass of the SMBH  • is set to the mass of the SgrA* SMBH,  • = 4 × 10 6 M ⊙ (e.g.Do et al. 2019;GRAVITY Collaboration et al. 2019, 2020;Event Horizon Telescope Collaboration et al. 2022) for Milky Way center comparison.
We sample the semi-major axes  ★ for the disk stars and binary star center-of-masses from a distribution of  () in such a way that the surface density Σ() of the disk is The power-law slope of −2 is chosen to enable more straightforward comparison to previous simulation studies and Milky Way center observations (Paumard et al. 2006;Madigan et al. 2009) The effect of choosing a different disk surface density profile is discussed later in Section 8.4.3.The scale radius of the disk  0 and the surface density at the scale radius Σ 0 = Σ( 0 ) are fixed by setting the inner  in and the outer  out edge of the disk as well as its total mass  disk .We simulate in total three different disk masses, 1.0×10 4 M ⊙ (sample L for low-mass and literature setup), 3.75×10 4 M ⊙ (sample M for mid-mass) and 7.5 × 10 4 M ⊙ (sample H for high-mass).We fix the inner and outer edge of the disk models to  in = 0.05 pc and  out = 0.5 pc.In each disk model all the stellar orbits share a common eccentricity  init .We model in total nine different disk eccentricities from  init = 0.1 to  init = 0.9 with constant increments of 0.1.The initial inclinations of the orbits are small but non-zero, || < 0.01.The other orbital elements, the longitude of the ascending node Ω, the argument of periapsis  and the mean anomaly The initial asymmetry of the disk models is enforced by aligning the eccentricity (Laplace-Runge-Lenz) vectors of Eq. ( 2) of the disk stars, i.e. the pericenters of the stellar orbits lie all in the same direction.This is illustrated in the right panel of Fig. 1.

Stellar population properties
The masses of the stars in the disk model are follow a standard power-law initial mass function

Binary star population properties
Star formation in the vicinity of SMBHs is poorly understood, and binary star formation even less so, and there is a very limited amount of observations available to constrain the properties of the binary population in the disks around SMBHs.In order to proceed, we assume that some very general binary population properties from star clusters and field stars remain valid in extreme environment of the galactic nuclei.
First, we assume that massive stars prefer having (massive) companion stars with mass ratios of  =  2 / 1 > 0.1.The binary fraction as a function is chosen to reflect the observations in less extreme star formation environments (Moe & Di Stefano 2017;Winters et al. 2019).The binary fraction of B-type primary stars is > 90% while practically every O-star in the initial conditions is in a binary.Low-mass stars are predominantly single stars.
Next, the eccentricity distribution  ( b ) of binaries is close to the thermal eccentricity distribution, i.e.

𝑓 (𝑒
and that close binaries cannot have high eccentricities, Finally, the distribution of the initial semi-major The minimum  min and maximum  max semi-major axes are selected based on survivability of binary systems in the disk.Wide binaries are evaporated by the repeated perturbations from encountered stars, and very hard, especially eccentric binaries are destroyed by stellar mergers (e.g.Generozov & Madigan 2020).For this study we use two options, the standard population with  min = 10 −6 pc and  max = 10 −4 pc and and the hard population for which each binary shares a common semi-major axis of  ★ = 3 × 10 −6 pc.

Background potential
We model the effect of the old stellar population of the nuclear star cluster of the Milky Way (e.g.Schödel et al. 2014) on the orbits of the stars by including an external potential  ext in the simulation.The external potential corresponds to a spherical power-low background density distribution of We use the slope  = 1.4 throughout this study motivated by the observations of Genzel et al. (2003).The background scale radius  0 and density at the scale radius  0 = ( 0 ) are chosen so that the enclosed mass  () within  0 = 1 pc is  ( 0 ) = 5 × 10 5 M ⊙ (Genzel et al. 2003;Schödel et al. 2007).
The main effect of the spherical external potential is to cause the orbits of the disk stars to precess and thus contribute to the instability of the initially asymmetric eccentric disks.The fact that the spherical background cusp population is represented by a smooth external potential instead of live particles significantly lowers the computational costs, but also weakens the scalar and vector resonant relaxation effects in the simulations.The consequences of this are discussed later in this article.

List of simulations
The number of simulations performed for this study is in total 36.The simulations are divided into two sets.The first sample consists of the three sub-samples L, M and H (nine runs each) exploring the effect of initial disk mass  disk and eccentricity  init on the disk disruption mechanism and milliparsec star cluster formation around the SMBH.The main properties of these simulation runs are collected in Table 4.
The smaller simulation sample consists of seven runs further investigating the effect of the IMF, binary population properties, TDE prescription and PN equations of motion on the milliparsec stellar population around the SMBH, its formation and physical properties.The main parameters of the runs in this simulation sample are presented in Table 5.

The rapid instability of asymmetric eccentric disks
As expected from the fact that our disk models do not lie on the asymmetric eccentric disk stability regime  disk ≫  cusp , they are unstable and rapidly evolve into axisymmetric eccentric disks.That is, the directions of the disk star eccentricity (Laplace-Runge-Lenz) vectors lose their initial alignment and are rapidly randomized.The process lasts less than 1 Myr for all the models in simulation setups L, M and H.
In order to quantify the disruption process, we study the distribution of eccentricities of the stars in the disk, and the distribution of directions of the eccentricity (Laplace-Runge-Lenz) vectors.We define the pericenter directions of the stellar orbits as in which  x and  y are the planar components of the eccentricity vector ì .Here  0 is the reference direction of the initial eccentricity vectors.We show how the eccentricity vector directions  lrl of stars with different semi-major axes  ★ evolve in Fig. 4. The figure presents the initial evolution of the models L7, M7 and H7 during the first 0.4 Myr of the simulations with intervals of 0.1 Myr.
The lowest-mass model L7 evolves as similar setups in the literature.The asymmetric eccentric precesses in the retrograde direction with the precession rate being almost independent of the semi-major axes of the stars, as expected from Eq. ( 6) and the used background stellar profile power-law slope of  = 1.4.Stars with large  ★ precess somewhat faster than stars with smaller  ★ .The secular eccentric disk instability as in Madigan et al. (2009) proceeds in the simulation as the distribution of the eccentricity vector directions  lrl widens around its mean value as time goes on.After 0.4 Myr of evolution the mean pericenter direction of the model L7 is ⟨ lrl ⟩ = −2.7 ± 0.6 radians.
The models M7 and H7 behave qualitatively differently.While the outer parts of the disks above  ∼ 3 in precess in the retrograde direction just as the model L7, the inner parts of the eccentric disks either show little precession (M7) or actually precess in the prograde direction (H7).This is in contrast to the behavior of the setup L7 in which the stars precess retrograde at an almost the same rate, and the secular eccentric disk instability gradually widens the eccentricity vector direction distribution.In the model M7 ⟨ lrl ⟩ = −2.6 ± 1.4 at  = 0.4 Myr.The disk disruption especially in the model H7 is very different from the model L7, and is clearly not caused by the standard secular eccentric disk instability.While in L7 the pericenter direction distribution only gradually broadens, in H7 the initially narrow distribution rapidly flattens as the different parts of the disk are precessing at opposite directions.After 0.4 Myr in the model H7 the mean eccentricity vector direction is ⟨ lrl ⟩ = −3.1 ± 1.6.
As explained in Sections 2.1.2and 2.1.3,the spherical background always causes retrograde precession, and the disk selfcontribution to the precession can be prograde in the inner parts of the disk.Thus, we argue that the assumption of negligible disk selfprecession of Madigan et al. (2009) i.e. ( disk ≪  cusp ) is broken in the setups M and H.The disruption mechanism of the highermass disks differ from low-mass disks as the disks disrupt when its inner and outer parts precess at different directions.The prograde precession direction is not enough to stabilize the intermediate-mass disk models as opposed to the high-mass models (Madigan et al. 2018).We attribute this to the facts that a sizeable fraction of disk stars still precess retrograde, and that the prograde precession rate is different at different  ★ from the SMBH.Each panels shows the eccentricity vector direction  lrl of the disk particles at different semi-major axes from  ★ from the SMBH.The low-mass disk model L7 precesses prograde at a rate almost independent of  ★ while the standard secular eccentric disk instability of (Madigan et al. 2009) proceeds as the distribution of  lrl around its main value broadens.For the model M7 and especially the model H7 the beginning of the asymmetric disk disruption process is more rapid and qualitatively different than in the low-mass model.The reason for this is that in the models 7 and 7 the distribution of the disk to the orbit precession becomes non-negligible, breaking one of the assumptions of the eccentric disk instability model of (Madigan et al. 2009).While the outer parts of the disk precess retrograde, the inner parts of the disk either precess at a slow rate (M7) or precess prograde (H7).The initial alignment of the eccentricity vectors in setups M7 and H7 is rapidly lost as the different parts of the disk precess to opposite directions at different rates depending on  ★ .For the three models the mean values for  lrl after Precessing elliptical orbits in thin, axisymmetric systems can give rise to spiral patterns, (e.g.Kalnajs 1973;Binney & Tremaine 2008;Harsoula et al. 2021) the most famous example being grand design spiral arms of disk galaxies.Spiral patterns indeed develop in the runs of our more massive setups  and , but they are shortlived as the disk disruption proceeds.One such spiral pattern from the run H7 is shown in Fig. 5 compared to the corresponding spiralfree disk in the low-mass model L7 at the same time.The spiral patterns are always one-sided and right-handed as the outer parts of the disks precess retrograde and the inner parts in the prograde direction.

Eccentric disk instability in the intermediate-mass regime
In the previous section we demonstrated that the nature of the eccentric disk instability is different in the low-mass disk model L compared to the mid-and high-mass models M and H.This is because the self-contribution of the disk to the orbit precession becomes important in models M and H, breaking one of the as- sumptions of the secular eccentric disk instability of Madigan et al. (2009).
We develop a simple semi-numerical torque model to explain the initial orbit precession directions and rates in our simulations.The model is mostly analytical containing a single numerical step.The instantaneous time derivatives of the reduced angular momentum vector ì ℎ and the eccentricity vector are shown in Eq. ( 5).Assuming that the non-Keplerian perturbations ì  are small, the total change of ì  per orbital period is obtained by integration the expressions over the unperturbed orbit.After a single orbital period the eccentricity vector is ì If the perturbation is planar then the reduced angular momentum vector ì ℎ remains constant.In order to proceed, one should obtain the expression for the perturbation ì ) at each point of the orbit, and integrate over a full orbital period to obtain Δ ì .For the spherically symmetric background the non-Keplerian perturbation ì  cusp is simply and can be evaluated from the cumulative mass profile  cusp () of the cusp in a straightforward manner.
As opposed to the spherical case, for perturbing potentials without spherical symmetry calculating the perturbing acceleration ì  disk is challenging and we instead proceed numerically.This is the numerical step of the simple torque model.In the asymmetric eccentric disk case the change of the eccentricity vector of a star depends on the disk model and the semi-major axis of the star  ★ .We construct disk representations following the recipe described in Section 4.2 consisting of  m mass elements with the mass of each element being  =  disk / m .Next, we discretize the unperturbed stellar orbit ì  (), ì () into  p discreet points ì  k , ì  k , equidistant in mean anomaly M. The perturbing force ì k is obtained by direct summing the Newtonian gravitational acceleration with a small Plummer softening of  = 10 −3 pc from each mass element .Thus, the change of the eccentricity vector Δ ì  over a single orbit is We also include the effect of the post-Newtonian PN1.0 perturbation from Eq. ( 3) into our analysis.The orbits of stars within approximately  ★ ≲ 10 −2 pc of the SMBH will always precess in the prograde direction, depending on the eccentricity of the star.In The simple numerical precession showing the change of the direction of the eccentricity vector per orbit Δ lrl / for disk masses 10 4 M ⊙ ≤  disk ≤ 9 × 10 4 M ⊙ and initial eccentricities 0.1 ≤  ≤ 0.9.The precession rate is measured at four different semi-major axes around the SMBH at  ★ = 0.05 pc, 0.1 pc, 0.25 pc and 0.5 pc.For  disk = 10 4 M ⊙ stars at all radii precess retrograde and the precession rate decreases by a factor of ∼ 3 from low to high eccentricities.As Eq. ( 6) yields only a decrease by ∼ 40% from the spherical background, the rest is attributed to the disk mass contribution itself.The inner edge of the disk ceases to precess retrograde around  disk ∼ 4 × 10 4 M ⊙ at all disk eccentricities, after which it always precesses prograde with an increasing magnitude.At high disk masses, the precession rate at the inner edge of the disk is approximately constant.This is further elaborated in the main text.When orbiting outside the inner edge, the precession rate is always more prograde and less retrograde at high disk eccentricities at all disk masses.The outer parts of the disk almost always rotate in the retrograde direction, except with high eccentricities ( init = 0.9) and high disk masses ( disk ≳ 7 × 10 4 M ⊙ ).
practise, we use  m = 10 5 mass elements per disk and  p = 10 4 sample points per stellar orbit.The final results are averaged by repeating the procedure 10 times with different random realization of the disk in each sample.
We present the precession direction analysis described above for a system consisting of a SMBH, spherical stellar background cusp and an asymmetric eccentric disk model in Fig. 6.The mass of the SMBH and the properties of the cusp are identical to the ones in the simulations of these study.The disk masses range from  disk = 10 4 M ⊙ to  disk = 9 × 10 4 M ⊙ with intervals of 10 4 M ⊙ .The initial eccentricity of the disk is  = 0.7.With low disk masses  disk ≤ 3×10 4 particle orbits the precess retrograde, except in near the SMBH with semi-major axes  ★ < 10 −2 pc where relativistic precession is important.With increasing disk mass, its contribution of the disk to orbit precession becomes increasingly large.The prograde precession rate decreases, especially in the vicinity of the disk inner edge, and becomes dependent on  ★ .At disk mass of  disk ∼ 4 × 10 4 the prograde picture qualitatively changes as a region of weak prograde precession appears around 0.05 pc ≲  ★ ≲ 0.1 pc.Now the disk has two prograde regions and two retrograde regions as the original outer retrograde region is split into two.The initially narrow prograde region gradually widens and the magnitude of precession increases, until at  disk ∼ 9 × 10 4 M ⊙ the inner retrograde region vanishes, and precession is prograde within 0.2 pc of the SMBH.
The simple precession model presented explains the results of the N-body eccentric disk simulations of Fig. 4. The low-mass literature setup L7 precesses retrograde at all radii, as expected, and the rate of precession is slightly larger at the outer parts of the disk.The outer parts of the model M precess retrograde just as in model L7, but the very inner edge of the disk slowly precesses in the prograde direction.This is exactly what is also seen in the panel  disk = 4 × 10 4 M ⊙ of Fig. 6: slow inner prograde precession and stronger outer retrograde precession.The behaviour of the highestmass model H7 is also well explained by the simple precession models with masses  disk ∼ 7 × 10 4 M ⊙ - disk = 8 × 10 4 M ⊙ : prograde precession within ∼ 3 in ∼ 0.15 pc of the SMBH and retrograde precession of equal magnitude outside this radius.
Whereas Fig. 6 shows the results of the simple precession model with initial eccentricity of  = 0.7, Fig. 7 presents the results for nine initial eccentricities between  = 0.1 and  = 0.9 with intervals of 0.1.The precession rates and directions of the stellar orbits are measured at four different semi-major axes from the SMBH, at  ★ = 0.05 pc (inner edge of the disk),  ★ = 0.10 pc,  ★ = 0.25 pc and  ★ = 0.50 pc (outer edge).
With low disk mass  disk = 10 4 M ⊙ , the precession is retrograde at every radii, and the precession rate is larger with lower eccentricities, as expected from Eq. ( 6).Going towards higher disk masses, the retrograde precession rates decrease and eventually turn prograde, beginning with disks of higher eccentricities  ≳ 0.7 as the self-contribution of the disk to the precession becomes gradually more important.The inner edge of the disk precesses prograde when  disk > 4 × 10 4 M ⊙ at all eccentricities while at  ★ = 0.1 pc this occurs after  disk > 8 × 10 4 M ⊙ .At larger semi-major axes from the SMBH, at  ★ = 0.25 pc and 0.50 pc, the precession direction is almost always retrograde, except at high disk masses ( disk ≳ 7 × 10 4 M ⊙ ) and high eccentricities ( ≳ 0.8) where the precession rate is small or weakly prograde.
When not orbiting at the inner edge of the disk, the precession rate always shows a dependence on the disk eccentricity.However, when the disk mass is high ( disk ≳ 7 × 10 4 M ⊙ ), the precession rate at the inner edge of the disk is almost independent of the disk eccentricity.As the contribution of the spherical cusp is eccentricitydependent (see e.g. the panel  disk = 10 4 M ⊙ in Fig. 7), the disk contribution at the inner edge must be eccentricity-dependent as well to result in a constant precession rate.The prograde precession rate being higher at lower eccentricities can be understood in the following terms.When orbiting at the inner edge of the disk, there is no disk mass inside the pericenter distance from the SMBH.An orbit with low eccentricity also has a low apocenter and most of the disk mass is located outsize of the orbit, causing prograde torque (see Fig. 2).With higher eccentricities and thus higher pericenters, there is increasingly more mass inside the orbit near the apocenter, making the precession less prograde by retrograde torques.We emphasize that the approximately constant precession rate at the inner edge of the disk is not a generic feature, but depends on the chosen masses and (surface) density profiles of the disk and the cusp, and the SMBH mass.The initially narrow eccentricity distributions rapidly broaden in all the models.The broadening phase is rapid, just as the disruption of the alignment of the eccentricity vectors.This reflects the nature of the asymmetric eccentric disk instability in the disks with different masses just as seen in Fig. 4 for the spread of the eccentricity vector directions.Model L disrupts due to the secular eccentric disk instability while especially the model H disrupts more rapidly as the different regions to the disk precess to different directions.In the models L the rapid phase is essentially over in ∼ 0.6 Myr while in the models M and H the corresponding times are ∼ 0.4 Myr and 0.2 Myr.With the higher initial eccentricity of  init = 0.7, stars can reach very high eccentricities up to unity while for the initial eccentricity of  init = 0.3, the majority of stars have eccentricities always less than  ★ ∼ 0.8.

The long-term evolution of the disk eccentricity distribution
The histograms presenting the eccentricity distributions of the disks in model M, in total nine different initial eccentricities, are shown in Fig. 9.For the analysis the binary systems in the disk are again treated as their center-of-masses.The eccentricity distributions are presented at two times, at  = 0 Myr and later at  = 5 Myr long after the unstable asymmetric eccentric disk has disrupted.At this point the eccentricity structure of the disks is in approximate equilibrium, as seen in Fig. 8 for the initial disk eccentricities of  init = 0.3 and  init = 0.7.After  = 5 Myr of evolution all the nine models produce a broad disk eccentricity distribution from the initially narrow distributions.The broad distribution always peaks at 0.2 ≲  ≲ 0.5 regardless of the initial disk eccentricity.All the models have stars on almost circular orbits with  ★ < 0.1.However, only eccentric disks models with moderate to high initial eccentricities,  init ≳ 0.4 produce stars with extremely high eccentricities of  ★ > 0.9.A secondary peak of high-eccentricity stars develops in the models with initial eccentricity higher than  init ≳ 0.7 with the distribution being weakly bimodal in the model with initial eccentricity of 0.7.In the initially very eccentric model with  init = 0.9 the high-eccentricity peak of the disk eccentricity distribution becomes very prominent.

The long-term evolution of the vertical disk structure
Next, we examine the evolution of the shape, and especially the vertical structure of the disk models.We first characterize the shape of the disk system by the ratios of its principal axis ,  and , with  being the longest and  the shortest axis.The axis ratios are calculated from the inertia tensor of the disk stars and binary centerof-masses within the 3D half-mass radii of the models ( h ∼ 0.2 pc) as in Zemp et al. (2007).Fig. 10 presents the evolution of principal axis ratios / and / in each of the nine models of the simulation sets L, M and H.
The evolution of the azimuthal direction, /, reflects the instability of the initially asymmetric eccentric disks.The intermediate axis ratio / rapidly reaches unity in all the models on a time-scale corresponding to the broadening time-scale of the initially narrow eccentricity distributions of 0.6 Myr, 0.4 Myr and 0.2 Myr for the models L, M and H, respectively.The minor axis ratio / first increases in from an initial near-zero value to moderate values of ∼ 0.2-0.5 within 0.2 ≲  ≲ 1 Myr in the models.After 1 Myr, the vertical evolution proceeds on a different, slower time-scale as.The minor axis ratio of the disk models keeps increasing in all the models, albeit more slowly than in the early disk disruption phase.The initially most eccentric disk models produce the most vertically extended disks with the final values of / being 0.22 ≲ / ≲ 0.65, 0.25 ≲ / ≲ 0.51 and 0.33 ≲ / ≲ 0.55 in the models L, M and H after 10 Myr of evolution.The highest values of / are obtained perhaps surprisingly in model L.There are two reasons for this.First, the in low-mass disk models of set L the vertical accelerations due to the disk gravity of the disk stars are smaller.Second, in the model L we do not include stellar evolution, so massive stars in the run do not die, so they can keep perturbing towards higher inclinations stars for a longer time.
Fig. 11 shows the root-mean-square (rms) inclinations  rms of the disk stars for the simulation setups L, M and H.The rms inclinations of the disk models evolve as expected from the evolution principal axis ratios in Fig. 10 with an initial rapid phase before ∼ 1 Myr, and more gradual increase afterwards.The initial rapid phase approximately follows the theoretical expectation of the  1/4 evolution of the rms inclination evolution (Lissauer 1993).The final rms inclinations after 10 Myr of evolution are 12°<  rms < 32°, 14°<  rms < 31°and 17°<  rms < 31°for the three models, with the initially most eccentric disk models resulting in higher rms inclinations of disk stars.

Counter-rotating stars in the disks
Fig. 12 shows the number fraction of counter-rotating stars  c =  counter / in the disk as a function of time in the disk models L, M and H. Initially, all the stars are rotating in the same direction.The disk models with initial eccentricities less than ∼ 0.4 end up having a very small (≲ 1%) fraction of counter-rotating stars.Disks initially more eccentric than this rapidly acquire a fraction of counter-rotating stars in the brief eccentric disk instability phase during the first Myr of the simulations, the time-scale being shorter for more massive disk models.This is consistent with the time-scales of the broadening of the initially narrow eccentricity distributions and the principal axis ratio / in Fig. 8 and Fig. 10.After the initial rapid phase the  c only mildly evolves.
The final  c strongly depends on the initial disk mass and eccentricity.The maximum counter-rotation fractions are  c ∼ 0.15,  c ∼ 0.25 and  c ∼ 0.34 for the three disk models L, M and H, respectively.The stars than flip their rotation direction have preferentially high eccentricities when they do so.For example in the run H7, approximately 90% of the stars have  ★ > 0.9 at the moment of they turn counter-rotating.A majority of stars flipping their rotation direction are located in the inner half of the disk,  ★ ≲ 0.3 pc.For the analysis the orbits of the binary stars around the SMBH are treated as the orbits of their center-of-masses.The originally narrow disk eccentricity distributions (in orange) considerably broaden in 5 Myr, with the initial rapid evolution phase lasting only ∼ 0.5 Myr.Interestingly, after 5 Myr, the eccentricity distribution is broad for all initial eccentricities peaking between 0.2 ≲  ★ ≲ 0.5.While all setups produce stars on almost circular orbits, setups with initial eccentricities less than ≲ 0.4 do not have extremely eccentric stars beyond  ★ > 0.9.For more eccentric initial setups with  ★ ≳ 0.7, a secondary high-eccentricity peak begins to develop, becoming very prominent at initial eccentricity of  init = 0.9.

Disk stellar content at later times
in the models M7b and H7b, the B-star disks are somewhat more vertically extended than the disks of O-stars in the less top-heavy models.The late-type stars have even a more vertically extended distribution, especially in the disk setup M 1.7 .As all the disk models were originally very thin, we attribute the differences of the vertical disk heights of stars of different masses to the standard disk heating process.Interestingly, the O-star disk of the run M 1.7 shows a weak warp structure.No warp structures are present in the other disks.
Finally, the number of O-type and B-type stars in binaries is still large in the disks at  = 6 Myr,  O+B = 87 (M7b),  O+B = 178 (H7b),  O+B = 567 (M 1.3 ) and  O+B = 798 (M 1.7 ).Initially, most of the B-type and O-type stars are located in binary systems.When comparing these numbers to the total numbers of B-type and O-type stars, we can see that the overall binary faction is still roughly ∼ 0.5 in the models at  = 6 Myr.

Surviving binary stars in the disks at later times
By the end of the simulations, the disk models still retain binary stars.We present the radial dependence of the binary fraction of massive stars ( ★ > 2 M ⊙ ) at  = 6 Myr in our simulations in Fig. 15.The radial separation of the binary star center-of-mass ranges from 0.04 pc to 0.5 pc in the analysis.In our models the binary fraction defined as  b =  b / N ranges from 0.4 ≲  b ≲ 0.6 near the outer edge of the disk to 0.05 ≲  b ≲ 0.25 at its inner edge.When comparing the orbital eccentricities of single stars and binary center-of-masses, we find that their distributions correspond to each other, except at high eccentricities which lack binary stars.This is well in line with the expectation that SMBHs drive nearby binaries to mergers, or unbind them.
Even though the binary fraction decreases towards the inner edge of the disk, there is still a small number of binaries orbiting at these small separations from the SMBH.The innermost binary center-of-mass in our simulations has a semi-major axis of  8.However, the vertical structure of the disks evolves more slowly than their azimuthal structure.We attribute this process to the dynamical heating of the disks in addition to the disk instability itself.The most eccentric initial models yield the most puffed-up, vertically extended thick disks.The final axis ratios / after 10 Myr of evolution range from ∼ 0.25 to ∼ 0.5 depending on the initial disk eccentricity while in the low-mass model L the highest / can reach even somewhat higher values of ∼ 0.6.0.04 ≲  ★ ≲ 0.065 pc with respect to the SMBH while having an eccentricity of 0.2 ≲  ★ ≲ 0.5, consistent with orbital elements of typical inner-disk single stars.

Escapers and high-velocity stars
In the disk simulations the binary star center-of-masses can reach high eccentricities and thus low periapsis distances from the SMBH.The binary components can then become unbound to each other due to the Hills mechanism.The component unbound to the SMBH can reach high velocities, as explained in Section 2.3.Other mechanisms to accelerate stars into radial outward orbits with high speeds include disruption of binaries when one component reaches the end of its life-time, and complex and chaotic strong single-binary and binarybinary interactions.In our simulations, unbound stars are removed from the simulation at the separation of 10 pc from the SMBH, corresponding to a minimum speed of  esc ∼ 59 km/s at this radial separation.When an escape event occurs, we record the ID number, the mass, the escape direction and the velocity of the escaping star.
In the simulations of this study the escapers are all single stars.Components of a previously disrupted binary may escape at very similar times, but the two stars are always well separated from each other at the escaper removal radius.
We show the number of escapers  esc as a function of the initial disk eccentricity in the models L, M and H in Fig. 16.We divide the simulation time into two parts, the early period before  < 2.5 Myr and the later period between 2.5 Myr <  < 10 Myr.The division time approximately corresponds to the life-time of most massive stars in our simulation, after which they are replaced with compact remnants (sets M and H).
Before 2.5 Myr, the primary reason for stars reaching the escape velocity is the binary break-up due to the Hills mechanism.We explicitly check the ID numbers of the escaping stars and stars on tightly bound orbits around the SMBH to confirm the Hills origin of the escaping stars in our simulations.Almost all the Hills binary break-ups occur at early times after the disruption of the asymmetric eccentric disks.In the low-mass setup L there are no escapers at early times in the simulations with the initial disk eccentricity less than  init < 0.7, after which the number of escapers steadily increases with increasing initial eccentricity to  esc = 13 at  init = 0.9.In the two more massive simulation setups M and H the trend is similar.Only a few escapers ( esc < 8) originate from strong few-body interactions not involving the SMBH at low initial disk eccentricities, after which  esc begins to increase.For the mid-mass set M this occurs around  init ∼ 0.5-0.6 for the set M and for the most massive set H after  init ∼ 0.3.The number of escapers increases rapidly, reaching  esc = 44 and  esc = 111 for the two sets.
After 2.5 Myr the dominating mechanism of producing escaping stars in the disk simulations is the binary break-up due to death and removal of massive stars.In the low-mass simulation setup L, in which the removal is not performed, there are only few escaping stars after  > 2.5 Myr, all due to strong interactions involving  ≥ 3 bodies.For the simulation setups M and H, there are ∼ 10-20 escaping stars in each simulation, with no dependence on the initial disk eccentricity, as expected.
We present the distribution of speeds of the escaping stars at  = 10 pc from the SMBH in Fig. 17.As before, we divide the escape events into early ( < 2.5 Myr) and late ( > 2.5 Myr) ones.As the number of escapers in many simulations is small, we study only the distribution of combined escape speeds from all escapes in our 27 simulations.The both early and late escape speed distributions are well characterized by a log-normal distribution.The early-time mostly Hills-driven escape speed distribution has a mean logarith-Figure 11.The rms inclinations of the disk stars and binary center-of-masses within 0.04 <  ★ < 0.5 pc in the simulated disk models.Already at  = 1 Myr the rms inclinations range between 6°<  rms < 25°, 8°<  rms < 27°and 10°<  rms < 24°in the models L, M and H, respectively.As with the minor axis ratio / in Fig. 10, the evolution of the rms inclination in the disk models slows down after 1 Myr.The final rms inclinations after 10 Myr of evolution are 12°<  rms < 32°, 14°<  rms < 31°and 17°<  rms < 31°for the three simulation sets.The early inclination feature around  ∼ 0.1 Myr in models M and H coincides with the appearance and dissolution of the spiral feature in the disk illustrated in Fig. 5.Besides the initial eccentricity, the fraction of counter-rotating stars strongly depends on the initial disk mass.The build-up phase of  c is rapid, consistent with the eccentric disk instability timescale in each model.After this, the evolution of  c is slow.mic escape speed of ⟨log 10 ( esc )⟩ = 2.70±0.38 corresponding to a speed of ∼ 500 km/s.The distribution is broad, and approximately 22% of the early escapers have a high escape speed greater than 1000 km/s.The two most extreme high-velocity stars move at the speed of 3500 km/s and 5000 km/s when exiting the simulation domain.After 2.5 Myr, the escape speeds are on average lower, but the early and late distributions overlap.The log-normal fit yields a mean logarithmic escape speed of ⟨log 10 ( esc )⟩ = 2.34 ± 0.31.This corresponds to a speed of 217 km/s.Only few late particles (< 4) reach velocities higher than 1000 km/s, and at two of these are late-time Hills ejections.
As pointed out by Šubr & Haas (2016), the anisotropic direction distribution of the escaping stars is characteristic for the Hills mechanism with the binary supply from an asymmetric eccentric disk.We confirm this result in our simulations, presented in Fig. 18.Most of the stars ejected by the Hills mechanism early in the simulations escape to the directions close to the disk mid-plane  = 0.At the azimuthal direction, the distribution shows a clear hemispherical asymmetry.The scatter of one standard deviation around the in the in-plane escape direction  is   = 57°,   = 45°and   = 32°f or the disk models H7, H8 and H9, respectively.
The anisotropic escaper direction distribution reflects the na-

Hills other escapes
Figure 18.The anisotropic distribution of early Hills escapers (in black) versus the more isotropic distribution of the late non-Hills escapers (in red).
The data shown combines the simulations H7, H8 and H9.In the Mollweide projection, the initial apoapsis direction of the disk stars lies on the left of the image while the initially aligned periapsis direction is at the center.The stars ejected by the Hills mechanism preferentially escape near the disk midplane  = 0 while the azimuthal direction  shows a clear hemispherical asymmetry.
ture of the supply of binary stars for the Hills mechanism.The direction of the escaping component of the disrupted binary is typically roughly opposite to the direction of the eccentricity vector of the just disrupted binary center-of-mass with respect to the SMBH.Thus, if the disk star eccentricity vectors are not yet completely misaligned even though the asymmetric eccentric disk is disrupting, the escaper direction distribution can still retain information of the orientation of the eccentric disk from the time the stars were ejected.
The late-time escapers not originating from the Hills mecha- nism show a more isotropic distribution, both having escapers all around the azimuthal angle  and at high inclinations far above and below the disk plane  = 0.

Accreted stars and TDEs
We to study the total number of stars accreted by the SMBH and star-star mergers in the simulations of the samples L, M and H as a function of disk mass and initial disk eccentricity.These results are shown in Fig. 19.There are more than ∼ 5 star-SMBH mergers in the disk models with initial eccentricity higher than 0.7 (L), 0.5 (M) and 0.1 (H).The number of total accreted stars strongly depends on the initial disk eccentricity.The vast majority of the accreted stars originate from almost parabolic orbits.The maximum number of total accreted stars is  acc = 162 in the simulation H9.Most accretion events occur early, immediately after the instability of the eccentric disks and the occurrence of stars on highly eccentric orbits.As most of the stars are accreted during the first Myr of the simulations, the merger rate in the simulation samples L, M and H roughly corresponds to a maximum TDE rate of 60-160 Myr −1 .
Star-star mergers occur in each of the simulation setups as presented in the bottom panel of Fig. 19.In contrary to the star-SMBH merger rate, the star-star merger rate does not show any dependence of the initial disk eccentricity.The actual number of star-star mergers ranges from ∼ 60 mergers in the models L to ∼ 180 in the setup H.
The number of accreted stars depend on the chosen merger prescription, the IMF and the assumed binary star population properties.We present the number of stars accreted by the SMBH ( acc ) as a function of the IMF slope and the initial disk mass of the simulation models in Fig. 20.As already demonstrated in Fig. 19, more massive eccentric disks result in a larger number of accreted stars,  acc = 7 (L7),  acc = 37 (M7) and  acc = 95 (H7).First, the number of accreted stars increases with the increasing disk mass as the number of stars available for disruption increases, as already shown in Fig. 23 and Fig. 19.Next, with a constant initial disk mass, the number of accreted stars with increasing (less top-heavy) IMF slope.The number of accreted stars ranges from less than  acc < 10 events in the setup L7 to almost  acc ∼ 10 3 in the model M 1.7 .The binary population properties have a relatively small effect on  acc when comparing runs M7, M7b, H7 and H7b.Runs without PN effects have less accreted stars as the mean orbital eccentricities are lower and pericenter distances higher compared to runs with PN effects included.This comparison can be only performed for the sample M 1.3 as in the setups with other IMF slopes the runs without PN didn't include a TDE prescription either.
Introducing the hard binary population has only a small to moderate effect (up to ∼ 30%) on the number of accreted stars as the setup M7b has  acc = 34 accretion events and H7b has  acc = 126.Changing the IMF of the mid-mass disk model into a less extremely top-heavy one (M 1.3 ) results in  acc = 115 stars accreted by the SMBH, a number comparable to the high-mass setup H7b.Turning off the PN equations of motion in the model M 1.3 noPN reduces the number of accreted stars by over a factor of 2.5 as the stars end up on average on less eccentric orbits due to the missing suppression of the resonant relaxation effects.Finally, the least top-heavy model M 1.3 undergoes  acc = 847 star-SMBH mergers.The comparison of  acc in the other models with and without PN than M 1.3 is not possible as the runs without PN also had TDEs not enabled.
The vast majority of the stars are again accreted early.If the number of accreted stars within the first Myr is translated into a TDE rate, the models M 1.3 PNTDE and M 1.7 PNTDE have TDE rates of  ∼ 90 Myr −1 and  ∼ 770 Myr −1 within 1 Myr of the simulation start.The model 7 has a TDE rate comparable to the setup M 1.3 PNTDE while the other extremely top-heavy models have rates of  ∼ 30 Myr −1 (M7) and  < 10 Myr −1 (L).Note that as before, the given TDE rates are not instantaneous but instead averaged over 1 Myr, and that instantaneous TDE rates can be considerably higher.At later times the TDE rates are by several orders of magnitude lower.Overall, the first-Myr TDE rates considerably differ in the simulations with different disk masses and IMF slopes.
Finally, we estimate whether the TDEs in our simulations would be in principle observable.If the tidal disruption radius R t is smaller than the direct capture radius  cap = 4  • / 2 (assuming a zero-spin SMBH), the star disrupts within the horizon and the event is not observable (e.g.Coughlin & Nixon 2022).By equating the order-of-magnitude tidal radius  t and the capture radius we estimate that in the main sequence all our accreted stars disrupt outside the tidal radius.This is because the relatively low (Milky Way) mass of the SMBH in our simulations.For more massive SMBHs the tidal forces acting on a star would be weaker, and direct captures become more common.In the giant phase the disruption is further facilitated by the increased stellar radii, so direct captures are not expected in the latter phases of the simulations either.

Radial surface density and eccentricity profiles
During the simulation time a number of stars end up having a smaller semi-major axis around the SMBH than the original disk inner edge radius,  in = 0.05 pc.Hereafter, we term these stars the milliparsec stars.We present the radial surface density profile Σ() of two of our models H1 ( init = 0.1) and H7b ( init = 0.7) in Fig. 21 after  = 6 Myr of evolution.The radial bin size used is Δ = 0.01 pc.Outside  in , the two surface density profiles are similar.However, inside the initial disk inner edge the profiles differ.In the model H1, the surface density profile flattens inside  ★ ∼ 0.05 pc, and finally drops inside  ★ ∼ 0.02 pc.The innermost stars in this model orbit at  ★ = 0.014 pc.However, the model H7b behaves qualitatively differently as the Hills mechanism has deposited stars at separations of few tens of milliparsecs from the SMBH.The model has a local minimum in the surface density profile near  ★ ∼ 0.04 pc, inside which the surface density steadily increases as a power-law cusp.
The central concentration of stars in the run H7b within  ★ ≲ 0.04 pc has a very steep three-dimensional density profile ( ★ ) ∝  −  ★ with a slope of  ∼ 2.4.This is roughly consistent with the results of (Fragione & Sari 2018), who predict a power-law density profile () ∝  −  with  = 9/4 = 2.25 when Hills mechanism deposits stars around the SMBH.This central profile is steeper than the steady-state Bahcall & Wolf (1976) cusp which has a density profile slope of 7/4.
Next, we study the radial mean eccentricity profile of the two models in the same Δ = 0.01 pc bins.The expectation is that in the models with a higher initial disk eccentricity the mean eccentricity of the milliparsec stars should be higher than in the initially more circular setups.This is because the disk migration origin milliparsec stars have similar mean eccentricity as the disk, and for the Hills mechanism origin milliparsec stars should be initially have  ★ ∼ 0.98.The mean eccentricity profiles are presented in Fig. 22.As expected, the model H1 has an almost constant mean eccentricity of  ★ ∼ 0.4 at all radii, including within  in .The mean eccentricity of the milliparsec stars in the model H1 is consistent with the main disk.However, in the model H7b the mean eccentricity strongly increases within  ★ ≲ 0.03-0.04pc, reaching  ★ ≳ 0.9 at the center.Again, this is consistent with the preferentially Hills mechanism origin of milliparsec stars in the model H7b.We note that the central mean eccentricity strongly depends whether PN effects are included in the simulations (H7b) or not (H7).This will be elaborated later in the study.

The number of milliparsec stars
In addition to producing high-velocity escaper stars, the Hills mechanism also deposits stars on tightly bound, very eccentric orbits around the SMBH.We trace the orbits of each star in the simulations from BIFROST snapshots (snapshot interval Δ snap = 10 −3 Myr) in order to find the stars which had semi-major axes smaller than  ★ < 0.04 pc.The choice of the limiting semi-major axis of  ★ = 0.04 pc is motivated by the local minimum of the surface density profile of the model with  init = 0.7 in Fig. 21.We also require that for the analysis the star has to reside within the milliparsec region for at least  = 100Δ snap .This procedure helps to filter out candidate stars which had an artificially small osculating  ★ in the snapshot.We also trace the origin of the milliparsec to find out which fraction on the stars originated from a Hills mechanism binary break-up.For a star to be classified as a Hills origin milliparsec star we require that it was originally part of a binary system, and that its initial companion has escaped the simulated system with a high velocity ( >  esc ).
The results of the milliparsec star tracing analysis are presented in Fig. 23.In the low-mass sample L, the total number of milliparsec never exceeds  tot = 16 stars.There are no stars of Hills origin in the models with a low initial eccentricity of  init ≲ 0.5, and the maximum number of the Hills mechanism origin milliparsec stars in the sample L is  Hills = 7 in the initially very eccentric models.It is also worthwhile to note that models with a low initial eccentricity produce more milliparsec stars of disk inner edge origin than more eccentric models.
The more massive simulation setups M and H have more milliparsec stars, including those with Hills mechanism origin, compared to the low-mass setup L, as expected.Disk models with initial eccentricities lower than  init ≲ 0.4 have no Hills origin milliparsec stars while the maximum number of disk inner edge origin stars is  disk = 17 and  disk = 66 in the models M9 and H9, respectively.With higher initial disk eccentricities the number of Hills mechanism origin milliparsec stars grows rapidly, and is in fact exponential as their numbers follow well the relation  Hills ∝ exp ( init ).The maximum number of  Hills reaches  Hills = 31 and  Hills = 71 in the setups M9 and H9.
Next, we study the effect of the IMF, PN equations of motion and including a TDE prescription on the number of milliparsec stars.The number of milliparsec stars in the simulation setups M7b, H7b, M 1.3 and M 1.7 is presented in Fig. 24.As opposed to Fig. 23, here we only show the total number of milliparsec (regardless of their Hills or disk edge origins) but now as a function of time.We categorize the stars based on their mass into types <B ( The number of B-type stars in the milliparsec population ranges from  B ∼ 10 in the setup M7b to  B ∼ 65-80 in the setups M 1.3 PNTDE and M 1.3 noPNnoTDE.In the runs with stellar accretion onto the SMBH enabled (label TDE in Fig. 24), the number of stars of each type is somewhat smaller than without the TDE prescription.The PN equations of motion (label PN) result in a slightly smaller overall number of milliparsec stars.This is because PN precession suppresses relaxation effects and thus the Hills origin milliparsec stars retain their high eccentricities (and hence low pericenter distances) and are more susceptible to tidal disruption.This is further discussed in the Section 7. an extremely top-heavy IMF (M7b, H7b) and setups with only a moderately top-heavy IMF (M 1.3 , M 1.7 ) is the number of late-type stars in the milliparsec population.While the extremely top-heavy models less than  <B ≲ 5 late-type stars, the moderately top-heavy models have  <B ∼ 130 and  <B ∼ 300, respectively.

Milliparsec population eccentricity distribution: the importance of the PN equations of motion
We characterize the evolution of the eccentricity distribution of the models M7b, H7b, M 1.3 and M 1.7 in Fig. 25.As the number of milliparsec stars is initially ( ≲ 0.2 Myr) small and the distribution does not always necessarily follow an exponential distribution  () ∝   , we instead use a non-parametric proxy for the eccentricity distribution.Namely, we investigate the fraction of eccentric milliparsec stars which we define as  0.5 =  ( ★ > 0.5)/, in which  is the total number of stars.For the common thermal eccentricity distribution  0.5 = 3/4.The fraction of high-eccentricity stars  0.5 is high at the time of the milliparsec population formation ( ∼ 0.5 Myr), reflecting the predominant Hills origin of milliparsec stars in the simulations.At later times whether the PN equations of motion are enabled in the simulations determines how the fraction of high-eccentricity stars evolves.Initial disk mass, number of milliparsec stars, TDE prescription and the IMF have a smaller effect than the PN.Without PN, the fraction of high-eccentricity stars ranges becomes smaller than the thermal fraction of 3/4 already at  ∼ 1 Myr, and ranges afterwards from  0.5 = 0.46 to  0.5 = 0.75.While the semi-major axis and hence the orbital energy of the milliparsec stars change very slowly, the norm of the angular momentum and thus eccentricity evolve rapidly.This is characteristic for scalar resonant relaxation (e.g.Merritt et al. 2011).With PN, the fraction of high-eccentricity stars always remains high,  0.5 > 0.8.This is because PN precession effects efficiently suppresses SRR effects.
We present the cumulative eccentricity distribution  () of the milliparsec stars at  = 6 Myr of the runs H1-H9, M7b, H7b and the in total five variants of the setups M 1.3 (three variants) and M 1.7 (two variants) in Fig. 26 and in Fig. 27.In the set of nine runs H1-H9 without PN, the cumulative eccentricity distribution is always more circular that the  () =  2 thermal distribution.The models with a large number of initially very eccentric milliparsec stars H7 and H9 have relaxed to have a mildly subthermal eccentricity distribution.The initially more circular models H1 and H3 are strongly subthermal, and especially the model H1 includes no milliparsec stars with eccentricities higher than  ★ ≳ 0.75.In Fig. 27, in the runs with PN equations of motion enabled, the eccentricity distribution remains strongly superthermal, i.e. there is little evolution after the formation of the milliparsec population as the scalar resonant relaxation effects are suppressed.With PN, simulation setups with more top-heavy IMFs result in more superthermal eccentricity distributions at  = 6 Myr.Without PN, the eccentricity distributions of the milliparsec stars are close to thermal (M 1.3 noPNnoTDE) or subthermal.The less top-heavy IMFs result in more subthermal eccentricity distributions.The effect of including or not including the TDE prescription appears to have only a small effect on the eccentricity distributions.
To summarize, the inclusion of the PN equations of motion in the simulations has a profound effect on the evolution of the eccentricity distribution on the milliparsec stars.With PN, the Hills mechanism origin milliparsec stars can retain their high eccentricity as the vector resonant relaxation effects are properly suppressed.Not including PN equations of motion in the simulations can result in close to thermal eccentricity distributions, such as in the simulation run M 1.3 noPNnoTDE.We emphasize that PN effects at least of the order of PN1.0 should always be included in the simulations of stars around SMBHs at milliparsec separations to properly capture the suppression of the resonant relaxation effects.The fraction of milliparsec stars  0.5 =  ( ★ > 0.5)/ with eccentricities higher than  ★ > 0.5 as a function of time in the various simulation setups.The thermal eccentricity distribution  () = 2 (dashed black line) has a fraction of  0.5 = 3/4.Regardless of the disk mass, the IMF and the TDE prescription, the simulation setups with PN equations of motion enabled (solid green, blue and purple lines) have a higher fraction of high-eccentricity stars than the thermal distribution.On the other hand, in the runs without PN (solid orange, yellow and red lines), the fraction of high-eccentricity stars is smaller than for the thermal distribution.This is due to the initially very eccentric orbits of the Hills origin milliparsec stars and the suppression of relaxation effects due to the PN precession, as discussed in the text.With the PN effects enabled (orange lines), the VRR is suppressed and the eccentricity distributions remain superthermal.The eccentricity distributions are the most concentrated to very high eccentricities ( ★ > 0.9) with the setups including extremely top-heavy IMFs (M7b, H7b).With less extremely top-heavy IMFs (M 1.3 PNTDE and M 1.7 PNTDE), the eccentricity distribution becomes gradually less and less superthermal.Without PN (blue and yellow lines), relaxation proceeds and the resulting eccentricity distributions are subthermal.Whether TDEs are included or not appear to have only a small effect on the overall results.

The (an)isotropy of the milliparsec stellar population: the importance of PN and background modeling
We show the distribution of the angular momentum vector directions of the milliparsec stars at  = 6 Myr in Fig. 28.The two extremely top-heavy models ( = 0.25) yield either disky (M7b) or lopsided (M7b) angular momentum vector direction distributions.
The milliparsec stars become more isotropic with less extremely top-heavy IMF, i.e. with increasing .Comparing the simulations without PN equations of motion to their PN counterparts, we see that without PN, the milliparsec star populations are more isotropic.
In this study, the spherical background cusp potential is represented as a smooth external potential instead of live particles in order to limit the computational costs.Including a live population of particles representing the old stellar population around the SMBH would enhance the relaxation effects, especially VRR (Panamarev & Kocsis 2022), and make the milliparsec stellar populations more isotropic than in the runs with a smooth external background.It has been shown in idealized simulation setups that relaxation effects from the background stellar population can rapidly render the stars at milliparsec sepatations from the SMBH isotropic, even when PN effects are properly modeled (Hamers et al. 2014).

Background and motivation
Even though stellar disks within the gravitational influence radii of SMBHs are presumed to be common in the galactic nuclei in the local Universe, resolved stellar disk dynamics at the level of single stars can only be observed from the nucleus of the Milky Way galaxy.
The second-closest galactic nucleus, the center of M31, contains a lopsided eccentric disk structure, but it is very massive compared to the surrounding spherical stellar cusp ( disk ≫  cusp ).Thus, comparison of the M31 eccentric disk dynamics to the eccentric disk simulations of this study for which  disk ≲  cusp is not very meaningful.Hence, the Milky Way center is the only available comparison for our simulations.
Recent observations indicate that the central parsec contains two almost perpendicular eccentric disk structures (CW and CCW disks, e.g.Bartko et al. 2009) and possibly two other outer filamentary structures (von Fellenberg et al. 2022).Obviously, our simulation results starting from initial conditions with a single asymmetric eccentric disk are far too simplistic to explain these complex disk structures of young, massive stars.However, we can still compare a number of observable properties of the individual disks to our simulation results.We focus on the mainly on the CW disk and to a lesser extent to the CCW disk and the outer filamentary structures.
The main disk structural parameters varied in our asymmetric eccentric initial conditions are the initial disk mass, the initial disk eccentricity and the IMF of the stellar population.Our full simulation sample probes the initial disk mass range from  disk = 10 4 M ⊙ to 7.5 ×  disk = 10 4 M ⊙ .Due to stellar evolution, extremely topheavy disk models can lose up to 50%-90% of their mass in  = 6 Myr, as discussed in Section 4.1.Thus, the initially massive disk models have masses corresponding to the present-day CW and CCW disk masses at later times.The initial eccentricities cover the almost entire parameter space from almost circular to very eccentric disks, i.e. 0.1 ≤  init ≤ 0.9.The IMF determines the number of stars in the disk, as well as number ratios of stars of different masses and spectral types.We tested both extremely top-heavy (power-law slope of  = 0.25) and less top-heavy IMFs ( = 1.3,  = 1.7).For most of the comparisons we choose the simulation state at  = 6 Myr which is consistent with the age estimates of the young, massive stars in the Milky Way center.Unless otherwise stated, the observed disk structural parameters are from von Fellenberg et al. (2022).

Disk star eccentricity distribution
We have shown in Fig. 9 that the late-time eccentricity distribution of the disk stars is to a large extent almost independent of the initial eccentricity of the disks.This is despite of the fact that the initial eccentricity distributions are very narrow and distinct from each other.Initial disk setups with eccentricities of 0.1 ≤  init ≤ 0.6 result in qualitatively very similar eccentricity distributions with a broad peak around 0.2 ≲  ★ ≲ 0.5.In the setups with a higher initial disk eccentricity than  init ≳ 0.6, a secondary peak of higheccentricity stars appears.This feature becomes more prominent with increasing initial disk eccentricity.The simulated eccentric disks at  = 6 Myr are consistent with the observed eccentricity distributions of the CW and CCW disks and the outer filaments.The CW and CCW disks have a median eccentricity of 0.4-0.5 while the outer filaments have a somewhat higher median eccentricity, 0.7.The CCW disk also contains high-eccentricity stars, but not beyond ≳ 0.9.However, the highest initial eccentricity of  init = 0.9 produces a large number of extremely eccentric stars ( ★ > 0.9) which appear to be missing from the Milky Way center disks, and it is thus disfavoured.Based on these results, it is difficult to infer the initial disk eccentricity 6 Myr ago with certainty from the CW disk eccentricity distribution observed today.

The missing disk warps
The CW disk appears to be significantly warped in both angles on the sky (Bartko et al. 2009;Genzel et al. 2010).The angle between the angular momentum vectors at the inner edge of the CW disk compared to its outer edge is ∼ 60°.Our simulated disks do not show any strong warp features.The disks of simulated O-type stars show warp-like features in Fig. 14 in the run M 1.7 and possibly in H7b, but both of these features are weak.Hence, we conclude that the simulation setups of this study do not reproduce the kind of warped disk features observed in the Milky Way center CW disk.
There are two plausible explanations for the missing disk warps in our simulations, both connected to the chosen initial conditions.First, our initial setups only include a single stellar disk embedded in a smooth spherical external background.If the disk stars were all test particles, their orbital planes would be conserved, also preventing any warp structure formation.In addition, the Milky Way center has at least two sub-parsec scale disks, the almost perpendicular and coeval CW and the CCW disks.It has been shown using Nbody simulations that the interaction of such two disks can lead to significant warping of the disks (Löckmann & Baumgardt 2009).
Second, Kocsis & Tremaine (2011) have shown using analytical arguments and Monte Carlo simulations that an initially circular, razor-thin disk embedded in a spherical background stellar population can rapidly develop a significant warp due to vector resonant relaxation.Spiral structure and warp formation in disks is also possible in N-body simulations, e.g.Perets et al. (2018) with self-consistent field methods.This brings back the question of the modeling of the spherical background in our simulations.In this study, we modeled the background using a smooth static external potential instead of a live cusp of stellar particles for computational cost reasons.Thus, our simulations lack the VRR effects on the disk stars originating from the spherical background.We already discussed the need of including a live background population in the context of the simulated milliparsec star eccentricity distribution in Section 7.3 and angular momentum vector direction distributions in Section 7.4.The missing disk warp issue adds a third motivation to pursue these significantly computationally more expensive simulations in the future.

Disk stellar population and the IMF
The observed numbers of different type of massive and evolved stars can be used to constrain the IMF and age of the Milky Way center young stellar population.Both the ratio of O/WR stars to B stars (  OWR−B =  O/WR / B ) and the ratio of O and B stars to WR stars  OB−WR =  OB / WR are sensitive to the IMF and the age of the stellar population (e.g.Lu et al. 2013).
In our simulations, the disk stellar type content at the observed age of the CW disk is largely determined by the IMF and the assumed solar metallicity.It should be noted that the definition of a B-type or an O-type star in this study is extremely simplified and only based on the initial mass of the star.The initial conditions of the simulation sets M and H with an extremely top-heavy IMF ( = 0.25) are constructed to have an equal number of type B and O-stars after  ∼ 6 Myr of evolution.Consequently, all our other models with a less topheavy IMF ( = 1.3 in M 1.3 and  = 1.3 in M 1.7 ) will have an excess of B-stars compared to the number of O-stars which is not observed in the Milky Way center (Perets & Gualandris 2010).While the models M and H have more O-type stars than B-type stars in the disks at  ≲ 6 Myr and thus qualitatively match the observed O/WR-B ratio in the disk, a more profound question is whether these models can simultaneously match the other observational constraints from the Milky Way center.
Finally, a crucial assumption in setting up our initial conditions was that the IMF is constant everywhere in the disk from the inner edge of  in = 0.05 pc to the outer edge of  out = 0.5 pc.Recent observations show that in the Milky Way center, the inner disk structures have a higher O-B star ratio than the outer filaments (von Fellenberg et al. 2022) while the IMF slope of all the stars within the central ∼ 0.5 pc is certainly higher than  = 0.25 (Lu et al. 2013).The possibility of an IMF gradient is also consistent with a number of predictions that the IMF could have been more top-heavy closer to the SMBH than in the outer parts of the star-forming gas disk (Hobbs & Nayakshin 2009).
Including an IMF gradient in our disk setup will open new possibilities to better satisfy the various observational constraints from the Milky Way center.Hence, the IMF gradient in the initial disk is an important step towards better-motivated numerical simulations of the Milky Way center disks in the future.The key question here is whether a simulation with an initial IMF gradient simultaneously reproduce the O-B number ratios in the disks and more other observational constraints than a simulation with a fixed IMF slope.

Binary stars in the disks
Little is known about the binarity of stars in the Milky Way center.Still, a few stars have been undoubtedly identified as binary stars, most of which are eclipsing binaries with O/WR components within the central 0.2 pc (Ott et al. 1999;Martins et al. 2006;Rafelski et al. 2007;Pfuhl et al. 2011).No binary stars have been detected in the S-cluster region (e.g.Chu et al. 2018).There are indications of decreasing binary fraction below  ≲ 0.02 pc from the central SMBH compared to the field binary fractions of massive stars (Chu et al. 2023), which is consistent with the simple theoretical expectation that SMBHs unbind binary stars or make them merge into single stars.
Our simulations in this study indicate that the binary fraction of massive stars the disks is decreased from further away, at  ≲ 0.1-0.2pc.Further away from the SMBH, the binary fraction remains unchanged and is  b ∼ 0.5 at the disk age of  = 6 Myr.
However, these results strongly depend on the assumptions of the binary population properties in the disk initial conditions.Still, the results show a possibility that the binary fraction of the massive stars in the CV and CCV disks and outer filaments in the Milky Way center at the present day could be moderate or high.Most of the O-stars in our simulations that end up milliparsec population via the Hills mechanism have tidally disrupted either in the main sequence or in the giant phase and, or are massive enough and simply have died of old age by  = 6-7 Myr.This is illustrated in Fig. 24.One of our simulations, M7b, has all its milliparsec O-type stars tidally disrupted already early in its evolution.A simulation sample far larger than included in this study would be needed to assess how common this early disruption scenario actually is.

The number of S-stars and the O/WR-B ratio of the disks
The S-cluster is observed to contain approximately 30-40 B-type stars (e.g.Gillessen et al. 2017) at the inferred age of  = 6 Myr of the disks (Paumard et al. 2006).Our simulation that produce milliparsec stars typically have from a few stars up to  ∼ 80 milliparsec stars.The number of milliparsec stars in our simulations depends on multiple parameters.In the following the notation  ↑ means that the increasing the value of the parameter  or enabling physical process  increases the number of milliparsec stars.The notation  ↓ indicates the opposite.The most important parameters affecting the number of milliparsec stars are the initial disk eccentricity ( init ↑, exponentially), initial disk mass ( disk ↑), IMF ( ↑) and binary star population properties (e.g. b ↓), the inclusion of the PN equations of motion (PN ↓) and the TDE prescription (TDE ↓).Thus, we conclude that our models can produce milliparsec stellar populations with a comparable number of stars as in the S-cluster in the Milky Way center.
We emphasize the most important question about the number of milliparsec stars in our simulations is not whether a model produces an exact match to the observed number of S-stars.Rather we argue that a more relevant question is whether a model can produce a correct number of S-stars while simultaneously also matching the numbers of O and B-type stars and their ratio at  = 6 Myr.The same applies for other observational constraints from the Milky Way center disk structures.The work in the literature has not extensively focused in this important point while studying the Hills channel of S-cluster formation via direct numerical simulations (Perets & Gualandris 2010).It seems that there is no easy, evident solution for the problem.
The initial conditions of the simulation sets M and H with an extremely top-heavy IMF ( = 0.25) are constructed to have approximately an equal number of type B and O-stars (100-200) in the disk after  = 6 Myr of evolution.The specific models M7b and H7b yield  B ∼ 10 and  B ∼ 30 B-type stars in the milliparsec region, and as demonstrated above their numbers can be adjusted up or down by changing the various parameters in new runs.We argue that these runs are our best match in both the number of S-stars and the numbers of B and O stars and the O-B ratio in the eccentric disk.

IMF gradients in the initial disk?
It is probable that the IMF slope  in the central parsec of the Milky Way is not as top-heavy as  = 0.25 used in our most topheavy models.As the models have an equal number of O-stars and B-stars approximately at  = 6 Myr, a population consisting predominantly of O/WR stars at that time as observed in the CW disk (von Fellenberg et al. 2022) would require a close-to-flat IMF.The most top-heavy IMF in the literature for the central disks is  = 0.45 ± 0.3 (Bartko et al. 2010).Less top-heavy IMF slopes have been proposed, e.g. = 0.85 (Paumard et al. 2006) and  = 1.7±0.2(Lu et al. 2013).We find no apparent radial mass segregation in the eccentric disks during the simulation runs i.e. more massive Ostars do not migrate towards the inner edge of the disk as shown in Fig. 13.This suggests that observed O/WR-B ratio in the CV disk might be primordial and not caused by stellar migration.An IMF gradient suggested by hydrodynamical star formation simulations (e.g.Hobbs & Nayakshin 2009) is an attractive solution to the O/WR-B ration and the global IMF slope issue.Near the inner edge of the disk the IMF is almost flat ( ≈ 0 with a cutoff  max ) and the stars there are predominantly O-type stars.Moving outwards from the inner the disk  increases, resulting in a less extreme top-heavy global IMF in the central parsec.Performing simulations using disk initial conditions including an IMF gradient is an interesting prospect and we will perform such simulations in future work.

Late-type stars in the S-cluster
Main sequence late-type stars (in the sense of later type than B) i.e. stars below the mass of a few M ⊙ cannot be yet resolved as individual stars in observations of the Galactic center.Our simulations show that the number of late-type stars is strongly dependent on the IMF of the initial disk.As demonstrated in Fig. 24, there are less late-type stars than B-type stars if the IMF is extremely top-heavy ( = 0.25) by a factor of ∼ 3-5.On the contrary, in the simulations with a less extremely top-heavy IMF ( = 1.3 and  = 1.7) there are more late-type stars than B-stars in the milliparsec population.The ratio of the number of late-type stars and the B stars  <B / B decreases from  <B / B ∼ 5 to  <B / B ∼ 3 when the IMF slope increases from  = 1.3 to  = 1.7.Thus, the number of late-type stars in the S-cluster, if possible to directly observe in the future, can provide useful information from the IMF of the Milky Way center disks.

S-stars and the original disk eccentricity
We noted in Section 8.1.2that based on the observed present-day CW and CCW disk eccentricity distributions it is difficult to infer the original eccentricity distribution of the disks at the time of their births 6 Myr ago.This is because almost all of our initial disk setups yield eccentricity distributions comparable to the observed distribution after 6 Myr of evolution despite the very large differences of the initial models.Only the most eccentric models with  init = 0.9 are disfavoured as they result in a very large number of stars with  ★ > 0.9, which is not observed in the disks today.On the other hand, in order to Hills mechanism to operate in the simulations, binary star center-of-masses on highly eccentric orbits are required.In our simulations, Hills origin milliparsec population stars begin to be produced when the initial eccentricity of the disk exceeds  init ≳ 0.5-0.6.Based on these facts, an approximate simulation-based constraint can be placed on the initial disk eccentricity, 0.6 ≲  init ≲ 0.9.It should be noted that this simulation constraint is relies on the assumption that Hills mechanism is the source of a majority of young stars in the S-cluster.

The eccentricity and angular momentum vector direction distributions: the need for a live background stellar cusp
The fact that the eccentricity distribution of our simulated milliparsec stars is too superthermal and their angular momentum direction distribution anisotropic too compared to observations of the S-cluster is not necessarily an outstanding problem.This is because of the missing live stellar background population in this study and the use of a smooth static spherical external potential instead.It has been established in the literature in simplified models that the spherical background can render the inclination (angular momentum direction) distribution isotropic (Šubr & Haas 2016) and relax the initially very eccentric stellar orbits close to the observed eccentricity distribution (Hamers et al. 2014).As already discussed in the Section 8.1.3,the vector resonant relaxation from the background could also naturally alleviate the missing disk warp issue in our simulations.Including a live spherical background population is thus an utmost priority for future simulation studies.

Asymmetric eccentric disk instability mechanism
Numerical eccentric disk studies in the mass range  disk ≲  cusp have often focused in studying systems resembling one of the disky structures and the central S-star cluster in the Milky Way center.However, only relatively few self-consistent simulations including the evolution and the disruption of the eccentric disks and following the formation of the S-cluster via the Hills mechanism have been reported in the literature.
Representative simulation studies such as Madigan et al. (2009) or more recently Generozov & Madigan (2020) and Generozov (2021) start from a single eccentric disk with a mass comparable to the present-day Milky Way center CW disk ( disk ∼ 10 4 M ⊙ ).After the eccentric disk instability the most eccentric orbits are identified and populated with binary stars.The binary star encounters with the SMBH are then integrated in isolation using a few-body code.Importantly the chosen initial disk mass is the inferred presentday Milky Way center CW disk mass, even though the disk mass at its formation 4 Myr -8 Myr ago was higher due to subsequent stellar mass loss and deaths of massive stars.The mass loss in 6 Myr can be substantial, especially if the IMF of the original stellar population was extremely top-heavy.We have shown in Section 5.1 and Section 5.2 that if the initial mass of the asymmetric eccentric disk is higher than few times ∼ 10 4 M ⊙ , its exact disruption mechanism differs from the standard secular eccentric disk instability of Madigan et al. (2009).Most importantly, the assumption of the negligible contribution of the disk on the precession rates of the stars breaks down.The disk self-contribution has important non-trivial consequences.The asymmetric eccentric disk is still unstable in this intermediate-mass regime  disk ∼  cusp within a certain radius, as opposed to the stable case  disk ≫  cusp as observed in the M31 center.The nature of the instability in the intermediate-mass disk regime is a precession-direction instability: the asymmetric eccentric disk is very rapidly ( ≲ 0.6 Myr) disrupted as different parts of the disk precess at different directions.While the outer parts of the asymmetric eccentric disk precess to the retrograde direction just as in Madigan's secular eccentric disk instability, the inner parts of the disk precess in the prograde direction due to torques from the disk itself.A curious consequence of the precession direction instability is the appearance of a short-lived right-handed one-sided spiral pattern in the disrupting disk in Fig. 5 due to precessing elliptical orbits of the disk stars.

TDE rates in eccentric disk setups
We note that TDE rates from our simulations can appear somewhat lower than in the eccentric instability studies in the literature, namely Madigan et al. (2018), who report a very high TDE rate of 0.3-3 yr −1 per galaxy.In our simulation models, the TDE rates range roughly from 10 Myr −1 up to 1000 Myr −1 depending on the initial setup.
There are several reasons for the difference.Briefly, Madigan et al. ( 2018) assumes a significantly more massive stellar disk ( disk = 10 6  ⊙ ) than we have adopted for this study, and the given TDE rates are averaged over somewhat smaller time windows (0.03 Myr-0.3Myr) than in our study.In more detail, first, their TDE rate estimate does not originate directly from a simulation, but is based on a number of assumptions, namely the M31-inspired SMBH-disk mass ratio yielding  disk = 10 6  ⊙ .Next, based on numerical simulations, it is estimated that ∼ 10% of the stars (roughly  ★ ∼ 10 5 ) could disrupt during 10 2 -10 3 orbital periods.As the orbital period at the initial inner edge of the disk is  ★ ∼ 300 yr, the estimate of 0.3-3 yr −1 per galaxy is reached.
In our simulations, the ratio of disrupted stars to the total number of stars in the simulation typically ranges from ∼ 5%-10%, a value only slightly lower than adopted by Madigan et al. (2018), due to our more detailed TDE procedure.Thus, the differences of the TDE rates can be fully attributed to the different assumed IMFs and chosen initial disk masses, as both less top-heavy IMFs and more massive disks provide more stars to be disrupted by the SMBH.

The importance of PN effects
The work closest resembling the approach of our study is Šubr & Haas (2016).The authors start with a disk setup with a topheavy IMF and a high binary fraction, but include no external background potential, stellar evolution or PN equations of motion.The asymmetric eccentric disk unstable due to axisymmetric planar von Zeipel-Lidov-Kozai effect and produces eccentric binary centerof-mass orbits and results in a formation of a star cluster within the inner edge of the disk.However, we note that in simulations studying the Milky Way center S-cluster formation at least a smooth background potential should be included as one otherwise ignores important mass precession effects around the SMBH.Šubr & Haas (2016) find that a thermal S-cluster eccentricity distribution can be reached within a reasonable time, and the observed isotropic angular momentum vector distribution can be reproduced as well if a small spherical particle population ( = 500 stars with  ★ = 1 M ⊙ ) is included.
However, we caution that these results may be overly optimistic because of the pure Newtonian nature of the simulations.We have shown in Sections 7.3 and 7.4 that the milliparsec star eccentricity and angular momentum direction distributions at  = 6 Myr in our models crucially depend on whether PN equations of motion are included in the simulation or not.Without PN, the relaxation effects are artificially strong, rendering the simulated central star clusters on average less eccentric and more isotropic than they should in reality be.We emphasize that PN accurate equations of motion should always be used in simulation codes when studying the stellar dynamics in the vicinity of SMBHs.

The importance of spherical background modeling
It has been questioned in the literature whether the observed CW disk thickness supports the idea of an initially dynamically cold CW disk.Cuadra et al. (2008) argued against the cold-disk origin of Milky Way center stellar disks because of too low rms inclinations and the lack of counter-rotating stars in their simulations.The main difference between the simulations of this study and Cuadra et al. (2008) is that their study uses equal-mass stars in the disks, and does not include external potential.We note that the maximum rms inclinations in our studies are approximately two times higher, and our simulation setup includes runs with moderate counter-rotation fractions up to ∼ 30%.We attribute this difference to the eccentric disk instability, which is partially driven by the spherical background which the runs of Cuadra et al. (2008) do not include.
The simulations of this study use a smooth spherical external background potential instead of a live stellar background population, and no N-body studies using a setup similar to ours and including the live background have been reported in the literature.As the background is smooth, relaxation effects of the disk and the milliparsec stars caused by the background are neglected.A number of issues in our simulation results can be related to the missing relaxation from the background.First, our disks do not show any significant warping.It has been shown using Monte Carlo simulations that vector resonant relaxation induced by the background stellar population can cause strong warping of the disk (Kocsis & Tremaine 2011).The second issue is the S-cluster relaxation problem: which physical relaxation processes can thermalize and isotropize the S-cluster in its age?Hamers et al. (2014) argues that including an observationally motivated (e.g.Hopman & Alexander 2006) background live stellar cusp is enough to relax the eccentricity distribution of the to the observed superthermal distribution ( () =  2.6 , Gillessen et al. 2009), even when PN terms are included in the equations of motion of the stars.Their setup involves integrating  = 19 S-stars with  ★ ≲ 0.032 pc on initially highly eccentric orbits 0.93 ≤  ★ ≤ 0.99 as test particles in a background of  = 4.8 × 10 3 stellar remnants with a density profile slope of  = 2.Each of the remnant has a mass of  ★ = 10 M ⊙ .Already at  = 6 Myr, a superthermal distribution with  () =  2.6 is reached, consistent with observations of the orbits of the S-cluster stars.
Testing the realistic live background relaxation scenario would be extremely desirable for our eccentric disk and milliparsec popula-tion models as well in order to test whether the additional relaxation effects are strong enough to make our milliparsec star eccentricity distributions less superthermal and more isotropic.Unfortunately, including a stellar component extending well beyond the outer edge of the eccentric disk would significantly increase the particle number and thus computational costs.Assuming the smooth external density profile used in this study extending to 4 out = 2 pc and a background particle (stellar remnant) mass of 10 M ⊙ the background particle number is  bkg ∼ 1.5 × 10 5 , two-three orders of magnitude more than in the pure disk simulations of the current study.While simulating a dense star cluster of a few times 10 5 -10 6 particles for 10 Myr is nowadays not a problem per se (Wang et al. 2016(Wang et al. , 2020)), the fact that the star cluster is fully embedded in the influence radius of the SMBH in which the particle velocities and accelerations are high makes the problem far more computationally challenging than the traditional N-body star cluster simulations (e.g.Aarseth 2003).Even though the simulations of million-body systems around SMBHs are computationally expensive, they are not prohibitively so and have recently become feasible.We highlight the simulations direct N-body simulations of Panamarev et al. (2019) as an example of such studies.Motivated by this, we will perform eccentric disk simulations with a live nuclear star cluster background in a future study.

Top-heavy IMF cutoff and very massive stars
The IMF in the Milky Way center disks is top-heavy (e.g.Paumard et al. 2006;Bartko et al. 2010;Lu et al. 2013), and at the inferred age of 6 Myr (Paumard et al. 2006) stars with initial masses above ∼ 30 M ⊙ have already died.If the IMF is extremely top-heavy ( ≲ 1.0), in addition to the power-law slope of the IMF, another important question is whether the IMF has a high-mass cut-off at  max .For these extemely top-heavy IMFs the cutoff determines the number of very massive stars (VMSs).This study adopted a cut-off of  max = 120 M ⊙ but this choice is of course somewhat arbitrary.For the most top-heavy IMF we studied,  = 0.25, the fraction of very massive stars (> 100 M ⊙ ) is ∼ 13%) for  max = 120 M ⊙ .If the cutoff mass is increased to  max = 250 M ⊙ , the fraction of stars with > 100 M ⊙ is ∼ 50% and > 200 M ⊙ is ∼ 16%.With an even larger cutoff mass of  max = 500 M ⊙ the fractions of stars over 100 M ⊙ , 200 M ⊙ and 400 M ⊙ are 70%, 50% and 16%, respectively.Thus, the number of VMSs is extremely dependent on the IMF cutoff mass for extremely top-heavy IMFs.From the Nbody code point of view, higher stellar masses pose no problem in the mass range in which stellar evolution tracks from fast population synthesis codes are available.In the literature single stellar evolution tracks up to  max = 500 M ⊙ exist (e.g.Szécsi et al. 2022).
Adopting a high value for the cut-off mass  max for an extremely top-heavy IMF would result in VMSs both in the disk and the milliparsec region early in their evolution.The maximum residence time of such massive stars near the SMBH is only a few Myr due to their short life-times, smaller than the typical estimates of the current age of the S-cluster stars.The very massive stars naturally act as massive perturbers in the S-cluster, enhancing ZLK oscillations and resonant relaxation effects, and are conveniently already dead by the present day.The effect of a massive perturber (typically in intermediate-mass black hole) on the relaxation of the S-cluster has been studied in the literature (e.g.Generozov & Madigan 2020).The remnants of such stars (e.g.Gravity Collaboration et al. 2023) are an attractive candidate for the detection of gravitational waves from the galactic center.

Prescriptions for tidal disruption events
The massive stars in the S-cluster can prematurely die if they are tidally disrupted by the SMBH, ending their perturbing effect on the S-stars.TDEs and the dynamics of S-stars are inherently coupled.This is because after full tidal disruption, the disrupted star does not contribute to resonant relaxation effects or ZLK oscillations of the S-star orbital elements.As TDEs occur in reality, also simulations of S-cluster stars should always have a recipe for tidal disruption of stars.The chosen TDE prescription and its details affect the longterm evolution of the simulated S-clusters.We show in the middle panel of Fig. 24 that whether a TDE prescription is included in the simulation can affect the number of S-stars by up to 20%-30%.
In this study, we adopted a somewhat simple prescription based on the tidal radius argument (Kochanek 1992) with a correction factor based on detailed hydrodynamical simulations of TDEs (Ryu et al. 2020).We note that the tidal disruption simulations of Ryu et al. ( 2020) assume a main-sequence stellar structure.It has been shown that the tidal disruption radius strongly depends on the stellar age due to the evolved structure of the star, and stars close to the end of their main sequence lifetime may not be disrupted even at pericenter separations of  p ∼  t /3 (Golightly et al. 2019).In addition, Coughlin & Nixon (2022) suggest that the tidal disruption radius of likely <  t /5 for massive and evolved stars.Furthermore, in our simple model the stars can only fully disrupt, e.g.there are no partial nor repeated TDEs.Based on these considerations, our simulations most probably overestimate the TDE rate while somewhat underestimating the number of S-star like milliparsec stars.However, we emphasize that most of the tidal disruption events in our simulations occur during or just after the eccentric disk instability when  < 1.5 Myr.At this point even the most massive stars have not had time to evolve away from the main sequence, and the assumption of the main sequence tidal disruption is likely still valid.

The initial surface density profile of the asymmetric eccentric disk
For this work we chose a power-law eccentric disk surface density profile with a slope of −2 in order to facilitate comparison with previous simulation studies (e.g.Madigan et al. 2009) and Milky Way center observations (e.g.Paumard et al. 2006), and to restrict the already fairly large parameter space of the simulation models.
In the following we briefly discuss the implications of choosing a different initial asymmetric eccentric disk surface density profile.Consider a single star in the asymmetric eccentric disk with a certain semi-major axis .If the star lags or leads the main bulk of the disk, the star is pulled towards the main disk.This perturbing acceleration can be divided into a radial and a tangential component.The radial component affects the precession rate of the orbit while the tangential component dictates how rapidly the star recedes from the main disk, or is re-absorbed back into it, as illustrated in Fig. 3.
First, if the disk self-contribution to the precession rate of the stars is negligible (as in the case of low-mass disks), the different disk surface density profile does not affect the precession rate of the star.However, it has an effect on the time-scale of the secular eccentric disk instability as it affects the tangential component of the perturbing acceleration from the bulk of the disk.Second, if the disk itself is massive enough to affect the precession rate of the star (the case of intermediate-mass disks), steeper disk surface density profiles lead to increased retrograde precession, as more disk mass is located within the orbit of the star with a fixed semimajor axis .The opposite is true for shallower disk surface density profiles.Thus, we expect that for intermediate-mass asymmetric eccentric disks with a steeper surface density profile than the slope of −2 the instability of the disk follows more the characteristics of the secular instability of Madigan et al. (2009) than the more complex precession direction instability presented in this work.The opposite is again true for eccentric disks with a shallower surface density profiles: with a given stellar orbit there is less disk mass within the orbit, and prograde orbit precession is more prevalent.This facilitates the onset of the more complex precession direction instability.
We note that the considerations presented here are only approximate and simplified, and additional detailed N-body simulations would be required to assess the effect of the varying eccentric disk surface density profile slope on the results of this study.As the parameter space of the presented study is already fairly large, we do not perform such simulations for the purposes of this work.

SUMMARY AND CONCLUSIONS
We study the 10 Myr evolution of asymmetric eccentric disks with top-heavy single and binary stellar populations around SMBHs in galactic nuclei using direct N-body simulations.First we would like to better understand the eccentric disk instability in the regime in which the self-contribution of the disk to the orbit precession cannot be ignored and study the properties of the resulting disks.Second, we want to self-consistently study the formation of milliparsec populations of stars within the inner edge of the disk around the SMBH.Previous N-body studies have focused on more simplified initial conditions or individual aspects of the formation process (e.g.Madigan et al. 2009;Hamers et al. 2014;Šubr & Haas 2016;Generozov & Madigan 2020).Finally, we compare our simulated disks and milliparsec stars to observations of the Milky Way central parsec with the conclusion that no model studied simultaneously matches the various observational kinematic and stellar population constraints from the S-cluster and the central stellar disks.
We run a simulation sample of 36 asymmetric eccentric disk models with varying initial disk masses, the initial disk eccentricities, the stellar IMFs and the binary star population properties.The initial disk masses range from  disk = 1.0 × 10 4 M ⊙ to  disk = 7.5 × 10 4 M ⊙ .In the literature, most simulation studies focused on the low mass regime to match the present day observed Milky Way center CW disk mass with no stellar evolution.However, the disk is  ∼ 6 Myr old, and being top-heavy may already have lost a substantial amount of mass due to stellar mass loss and massive stars dying.Thus, the original initial masses of the Milky Way center disks most probably lie in the upper mass range.We simulate a wide range of initial disk eccentricities from  init = 0.1 to  init = 0.9 and three different power-law IMF slopes, the extremely top-heavy  = 0.25, and the moderately top heavy IMFs with  = 1.3 and  = 1.7.We also perform simulations with and without post-Newtonian equations of motion, and with and without a tidal disruption event prescription.The fixed parameters include the SMBH mass and the power-law static external potential representing the spherical background stellar cusp, both matched to Milky Way measurements.
The lowest-mass asymmetric eccentric disks undergo the sec-ular eccentric disk instability described in Madigan et al. (2009) and produce a thick, axisymmetric disk in less than 1.5 Myr.Our more massive disk models are unstable as well, but the physical mechanism for the asymmetric eccentric disk disruption is different from Madigan et al. (2009).For disk masses higher than a few times 10 4 M ⊙ the contribution of the disk itself to the precession of disk star orbits becomes non-negligible, breaking one of the assumptions of the disk instability model of Madigan et al. (2009).The negligible contribution is usually stated as  disk ≪  cusp .Based on our results, the disk contribution to the orbit precession already becomes important when  disk / cusp ∼ 0.07.The reasons for this surprisingly low disk-cusp ratio are that the mass distribution of the disk is planar as opposed to the more massive spherical background, and that the disk is asymmetric.While in the standard secular eccentric disk instability the (retrograde) orbit precession is only caused the spherical background cusp, in the more massive setups the selfcontribution of the disk to the orbit precession becomes important.
As the eccentric disks are asymmetric, the net torque over an orbit can result in prograde precession in the disk.The prograde precession proceeds from the inner towards the outer parts of the disk with increasing  disk .Thus, the alignment of the eccentricity vectors of the disk stars is rapidly lost, in less than 1 Myr, as the inner and outer parts of the disk precess in different directions.We term this form of the asymmetric eccentric disk instability the precession direction instability.The elliptical orbits precessing at different rates and directions also generate a short-lived one-sided right-handed spiral pattern in the disrupting disk.This structure torques the eccentricities of the disk stars, and the initially narrow eccentricity distributions broaden on shorter time-scales than in the standard eccentric disk instability.
We construct an analytical model for estimating the precession rates and directions of the stars in the disks as a function of  disk ,  init , the semi-major axis of the star  ★ and the density profile of the background cusp.The only numerical component in the model is the calculation of the acceleration of the star by the other disk stars at different points of the orbit.The simple model agrees well with the simulation results.The model can be used to predict the precession rates and directions of disk stars and thus the nature of the eccentric disk instability of arbitrary disk setups without running simulations.
The late-time eccentricity distribution of the disk stars always peaks between 0.2 ≲  ★ ≲ 0.5, regardless of the initial eccentricity of the disk model.Overall, the disk eccentricity distributions are very similar.Only when  init ≳ 0.7, a secondary high-eccentricity peak begins to appear in the distribution.This peak reaches the height of the moderate-eccentricity peak at  init ≳ 0.9.If the initial disk eccentricity is very high, the eccentricity distribution is bimodal with a secondary peak at  ★ ∼ 0.9.The vertical height of the disk rapidly increases from its initial thin configuration, and keeps increasing at a slower rate after the eccentric disk instability is over.The rms disk star inclinations reach values typically between 10°≲  rms ≲ 30°in the models.Isolated asymmetric eccentric stellar disks can attain a sizeable fraction of counter-rotating stars, but only if their initial eccentricity is  init > 0.4.Initially more eccentric and more massive disk models end up having up to 30% of stars on counter-rotating orbits.The majority of stars which flip their rotation direction have very high eccentricities ( ★ > 0.9) at the moment of the flip, consistent with the results of Madigan et al. (2018).
Simulation models with initial eccentricity higher than  init ≳ 0.5-0.6 produce up to ∼ 100 escapers and high(hyper)-velocity stars.The escaper and high-velocity star direction distribution from the Hills mechanism binary break-ups in our simulations is strongly anisotropic, consistent with the simulation results of Šubr & Haas (2016).While the escape direction distribution of non-Hills escapers is close to isotropic, stars ejected by the Hills mechanism predominantly escaper close to the disk mid-plane  = 0.In the azimuthal direction  the escaper direction distribution is a moderately narrow cone reflecting direction of the eccentricity vectors of the stars near the inner edge of the disrupting disk.The standard deviation of the in-plane escaper direction distribution gets narrower with the increasing initial disk eccentricity and can be as narrow as   ∼ 30°.Around ∼ 20% of the escaping stars have speeds exceeding 1000 km/s at  = 10 pc from the SMBH.
If a tidal disruption event prescription is included in the runs, a large number of stars can be disrupted by the SMBH, most during the first Myr of the simulation.Typically, the number of disrupted stars is at least comparable to the number of stars in the simulated milliparsec cluster.The TDE rate in our simulations during the first Myr strongly depends on the IMF, and varies from  < 10 Myr −1 to  < 770 Myr −1 .At later times the TDE rate drops by several orders of magnitude.
Populations of stars closely orbiting the SMBH form in our simulations inside the original inner edge of the disk.We term these simulated stars the milliparsec stars.For initial disk eccentricities  init ≲ 0.5 the milliparsec population consists predominantly of stars originating from the inner edge of the disk.With higher  init there are enough binary stars on a high-eccentricity low-pericenter orbits around the SMBH so that the Hills mechanism can operate.The number of Hills mechanism origin milliparsec stars depends on the initial disk mass, initial disk eccentricity (exponential dependence) the IMF and binary star population properties.The maximum number of B-type milliparsec stars is ∼ 80 in our simulations.Most our simulated milliparsec populations contain a few O-type stars until  = 6-7 Myr.Compared to the observed Milky Way center S-cluster with ∼ 30-40 members (Gillessen et al. 2017) we conclude that our simulation models can produce central clusters with a comparable number of B-stars in the S-cluster.However, an IMF slope of  = 0.25 (or lower) is required to simultaneously match the observed number of O/WR stars and the high O/WR-B ratio in the CW disk.The global IMF slope of the young stars within the central parsec of the Milky Way is not that low, for example Lu et al. (2013) provide an estimate of  = 1.7 ± 0.2.When comparing the eccentricity and angular momentum distributions of our simulated milliparsec stars to the S-cluster observations (e.g.Gillessen et al. 2009Gillessen et al. , 2017;;Burkert et al. 2023), we find that our clusters are too superthermal in eccentricity and too anisotropic.This fact points to a missing source of relaxation in our simulations.The most obvious candidate for the missing relaxation effect is the resonant relaxation from a realistic stellar cusp background population.In this study as we model the background with a static smooth external potential.
We emphasize the importance of including PN equations of motion in the simulations of stars around SMBHs, which is not always the case in the literature.Without PN effects, the resonant relaxation effects are not properly suppressed and relaxation is too strong compared to reality.In this unrealistic case, the simulated milliparsec can have a close to thermal eccentricity distribution, just as in the simulations of Šubr & Haas (2016) who did not include PN terms in their equations of motion.With the proper PN treatment, we find that the resulting milliparsec stellar population is strongly superthermal compared to the S-cluster observations.We note that simulation setups with a single, asymmetric eccentric disk embedded in a spherical background as our models are too simplistic to reproduce all the properties of the Milky Way center with all its complexity, such as the S-cluster, the warped CW disk, the almost perpendicular CCW disk and the outer filamentary structures.While our simulations can reproduce a number of individual properties of Milky Way center, such as the number of S-stars and the eccentricity distribution of the central disks, it is challenging to simultaneously fulfill the various observed kinematic and stellar population constraints in a single simulation.While the disk binary Hills mechanism origin of the S-stars remains plausible, more refined setups are needed to better understand the Milky Way center using numerical simulations.
We identify three key improvements for future models.First, resonant relaxation effects in our simulations might be underestimated as we used a fixed spherically symmetric external potential instead of including a live stellar background in the study.Due to the too weak scalar and vector resonant relaxation, the simulated Sclusters might remain superthermal in eccentricity and anisotropic in angular momentum vector direction distribution at  = 6 Myr.In has been shown in the literature that the scalar resonant relaxation originating from the background stars can indeed relax the initially very eccentric S-cluster star orbits within the age of the stars (Antonini & Merritt 2013;Hamers et al. 2014;Generozov & Madigan 2020).In addition, it has been argued that vector resonant relaxation from the background stars can induce strong disk warps (Kocsis & Tremaine 2011), which are observed in the Milky Way center CW disk but not appear in our simulation.Thus, including a live stellar background seems a promising approach to solve multiple different issues at once.At present, representing the full background population with individual stars is computationally too challenging.
Second, including an IMF gradient in our initial disk setups.Matching the observed number of O/WR stars and the O/WR-B star ratio in the CW disk in our simulations with a global IMF would require an almost flat IMF slope ( < 0.25).This substantially differs from the overall observed IMF slope in the central parsec, which is less top-heavy with  = 1.7 ± 0.2 (e.g.Lu et al. 2013).An IMF gradient would allow having the high O/WR-B ratio in the inner parts of the disks while globally having a less extreme total IMF.Such a variation of the IMF within the disks has been suggested by both recent observations (von Fellenberg et al. 2022) and star formation simulations (Hobbs & Nayakshin 2009).
Finally, a single and binary stellar evolution code properly coupled to our simulation code would allow a more straightforward comparison of observed Milky Way center stars of various spectral types and evolutionary stages to simulation results.We plan to couple the very recent SEVN binary stellar evolution package (Iorio et al. 2023) to our BIFROST code in a future study.

DATA AVAILABILITY STATEMENT
The data relevant to this article will be shared on reasonable request to the corresponding author.

Figure 1 .
Figure 1.Three types of disks around a SMBH.The most simple type of a disk is a Solar system like circular disk (left panel) in which all the disk stars have almost circular orbits.The orbits of the stars can also be eccentric in an axisymmetric disk (middle panel) with non-aligned the pericenter directions (eccentricity vectors) ì  i .If the pericenter directions are apsidially aligned (right panel), i.e. ì  i = ì  aligned , the eccentric disk becomes asymmetric.

Figure 2 .
Figure 2. A schematic illustration of Newtonian retrograde (left panel) and prograde (right panel) torques on a star (yellow symbol) orbiting the central SMBH (solid black circle).The eccentricity vector ì pointing towards the pericenter is defined in Eq. (2) while the reduced angular momentum vector ì ℎ is defined as ì ℎ = ì  × ì .At the apocenter the torque ì  due to a small non-Keplerian perturbation ì is ì  = ì  × ì ℎ.If the perturbing background is spherically symmetric (left panel), the torque and thus the precession direction is always retrograde.In the asymmetric case (right panel) both retrograde and prograde torques may occur at different points of the orbit, and the net torque integrated over the entire orbit determines the precession direction.
from a minimum stellar mass of  min = 0.08 M ⊙ to the maximum mass of  max = 120 M ⊙ .We study three different IMF slopes, the extremely top-heavy  = 0.25 and two moderately top heavy options  = 1.3 and  = 1.7.These choices are motivated by the observations of the Milky Way center which indicate a top-heavy IMF in the range from  = 0.45 ± 0.3 (Bartko et al. 2010) (central disks) to  = 1.7±0.2(overall innermost 0.5 pc).The high O/WR-B ratios observed in the CW and CCW disks (e.g.von Fellenberg et al. 2022) point towards even more extreme IMF in the innermost parts of the Milky Way center disk structures.

Figure 4 .
Figure 4.The beginning of the disruption of the disk models L7 (top row), M7 (middle row) and H7 (bottom row) at five times from  = 0.0 Myr to  = 0.4 Myr.Each panels shows the eccentricity vector direction  lrl of the disk particles at different semi-major axes from  ★ from the SMBH.The low-mass disk model L7 precesses prograde at a rate almost independent of  ★ while the standard secular eccentric disk instability of(Madigan et al. 2009) proceeds as the distribution of  lrl around its main value broadens.For the model M7 and especially the model H7 the beginning of the asymmetric disk disruption process is more rapid and qualitatively different than in the low-mass model.The reason for this is that in the models 7 and 7 the distribution of the disk to the orbit precession becomes non-negligible, breaking one of the assumptions of the eccentric disk instability model of(Madigan et al. 2009).While the outer parts of the disk precess retrograde, the inner parts of the disk either precess at a slow rate (M7) or precess prograde (H7).The initial alignment of the eccentricity vectors in setups M7 and H7 is rapidly lost as the different parts of the disk precess to opposite directions at different rates depending on  ★ .For the three models the mean values for  lrl after  = 0.4 are ⟨  lrl ⟩ = −2.7 ± 0.6 (L7), ⟨  lrl ⟩ = −2.6 ± 1.4 (M7) and ⟨  lrl ⟩ = −3.1 ± 1.6 (H7).

Figure 5 .
Figure 5.The number surface density of stars in the models L7 (top row) and H7 (bottom row) at  = 0 Myr and at  = 0.1 Myr.Each panel is stacked from five BIFROST snapshots.The SMBH position is marked with a cross symbol in each panel.A short-lived right-handed spiral structure (below the SMBH in the bottom right panel) forms due to the outer parts of the disk precessing in the retrograde direction, the mid-disk not precessing and the inner parts of the disk precessing into prograde direction.The low-mass disk does not show such a feature.

Figure 6 .
Figure6.The numerically estimated change of the eccentricity vector direction per orbital period Δ lrl of a star as a function of its semi-major axis  ★ around a SMBH in nine asymmetric eccentric disk models.The eccentricity of the disk is  = 0.7 and the masses of the disks range from 10 4 M ⊙ to 9 × 10 4 M ⊙ .The mass of the SMBH and the properties of the stellar background cusp are the same in each panel.The two white tics at the bottom of each panel mark the inner and outer edges of the disks.In the panels of the top row ( disk ≤ 3 × 10 4 M ⊙ ), near the SMBH ( ★ ≤ 10 −2 pc) the precession is prograde (blue color) due to relativistic periapsis advance while further away from the SMBH the precession is retrograde (in red) due to mass precession from the spherical background.With higher disk masses (middle and bottom rows) the picture qualitatively changes as the disk mass itself becomes an important source of precession.First, a new region of weak prograde precession appears between 0.1 pc ≲  ★ ≲ 0.2 pc when  disk = 4 × 10 4 M ⊙ .This prograde region broadens and the magnitude of Δ lrl further increases in higher-mass disks until the inner retrograde region completely vanishes when  disk = 9 × 10 4 M ⊙ .This simple model for initial disk precession directions explains well the features of the N-body simulations of Fig.4.

Fig. 8
Fig.8presents the evolution of the disk star eccentricities until the end of the simulations at  = 10 Myr.Binary systems are treated as their center-of-masses.Stars with semi-major axes smaller than  ★ = 0.04 pc or larger than 0.5 pc are excluded from the analysis.Two initial eccentricities of  init = 0.3 and  init = 0.7 are shown from the three models L, M and H of different disk masses each.

Figure 8 .
Figure 8.The broadening of the initially narrow eccentricity distributions of the disks with initial eccentricities of  init = 0.3 (red lines) and  init = 0.7 (in black).The lines show the Lagrangian eccentricities (in the sense of the Lagrangian radii of a cumulative mass profile)  5 ,  20 ,  40 ,  60  80 and  95 .For example, 95% of the particles have a lower eccentricity than  95 .Binary stars are treated as their center-of-masses in the analysis.With the initial eccentricity of 0.7, the disk stars can reach almost circular or radial orbits in less than ∼ 0.6 Myr.Disk with initially lower eccentricity also has stars with moderately eccentric orbits, but only after several Myrs of evolution.The main difference between the models L, M and H is how rapidly the eccentricity distribution spreads.For the model L7 the initial phase of rapid evolution is over after ∼ 0.6 Myr, while the time-scale is ∼ 0.4 Myr and ∼ 0.2 Myr for the models M and H, respectively.

ForFigure 9 .
Figure9.The eccentricity distributions of the nine simulations of the disk model M at 0.04 pc <  ★ < 0.5 pc at two times,  = 0 Myr and  = 5 Myr.For the analysis the orbits of the binary stars around the SMBH are treated as the orbits of their center-of-masses.The originally narrow disk eccentricity distributions (in orange) considerably broaden in 5 Myr, with the initial rapid evolution phase lasting only ∼ 0.5 Myr.Interestingly, after 5 Myr, the eccentricity distribution is broad for all initial eccentricities peaking between 0.2 ≲  ★ ≲ 0.5.While all setups produce stars on almost circular orbits, setups with initial eccentricities less than ≲ 0.4 do not have extremely eccentric stars beyond  ★ > 0.9.For more eccentric initial setups with  ★ ≳ 0.7, a secondary high-eccentricity peak begins to develop, becoming very prominent at initial eccentricity of  init = 0.9.

Figure 10 .
Figure10.The evolution of the shape of the disk models L, M, and H as characterized by the principal axis ratios / and / determined from the inertia tensors of the systems.The intermediate axis ratio / rapidly reaches unity, as expected, as the initial alignment of the disk star eccentricity vectors is lost.The time-scale of / reaching unity is very similar than the time-scale for the broadening of the initially narrow eccentricity distributions, from ∼ 0.6 Myr in model L to ∼ 0.2 Myr in model H, shown in Fig.8.However, the vertical structure of the disks evolves more slowly than their azimuthal structure.We attribute this process to the dynamical heating of the disks in addition to the disk instability itself.The most eccentric initial models yield the most puffed-up, vertically extended thick disks.The final axis ratios / after 10 Myr of evolution range from ∼ 0.25 to ∼ 0.5 depending on the initial disk eccentricity while in the low-mass model

Figure 12 .
Figure 12.The fraction of counter-rotating stars  c =  counter / in the simulated disks.The models with initial eccentricity  init ≲ 0.4 have less than a percent of counter-rotating stars.Above this  init , the  c rapidly increases to peak values of  c ∼ 0.15,  c ∼ 0.25 and  c ∼ 0.34 in the models L, M and H, respectively.Besides the initial eccentricity, the fraction of counter-rotating stars strongly depends on the initial disk mass.The build-up phase of  c is rapid, consistent with the eccentric disk instability timescale in each model.After this, the evolution of  c is slow.

Figure 13 .Figure 14 .Figure 15 .Figure 16 .Figure 17 .
Figure 13.The face-on projection of the disk showing the positions of Btype (black dots) and O-type (small red circles) stars after  = 6 Myr of evolution.

Figure 19 .
Figure 19.Number of stars merging with the SMBH during the simulations (top panel) and with other stars (bottom panel) in the simulation setups L, M and H.While the number of stars accreted by the SMBH strongly depends on of the initial disk eccentricity, the number of star-star mergers has no such dependence.

Figure 20 .
Figure 20.The number of stars accreted by the SMBH ( acc ) in selected simulation runs as a function of the IMF slope of the stellar population if the setup.First, the number of accreted stars increases with the increasing disk mass as the number of stars available for disruption increases, as already shown in Fig.23and Fig.19.Next, with a constant initial disk mass, the number of accreted stars with increasing (less top-heavy) IMF slope.The number of accreted stars ranges from less than  acc < 10 events in the setup L7 to almost  acc ∼ 10 3 in the model M 1.7 .The binary population properties have a relatively small effect on  acc when comparing runs M7, M7b, H7 and H7b.Runs without PN effects have less accreted stars as the mean orbital eccentricities are lower and pericenter distances higher compared to runs with PN effects included.This comparison can be only performed for the sample M 1.3 as in the setups with other IMF slopes the runs without PN didn't include a TDE prescription either.

Figure 21 .
Figure21.The radial face-on surface mass density Σ ( ★ ) of the stars between the semi-major axis  ★ = 0.01 pc and  ★ = 0.5 pc in the runs H1 ( init = 0.1, in blue) and H7b ( init = 0.7, in orange) at  = 6 Myr.The initial inner edge of the eccentric disk  in = 0.05 pc is indicated with a vertical dashed line.The two surface density profiles are similar outside  in .Inside this radius, the surface density profile of the initially more circular model H1 first flattens and then begins to decrease within  ★ ∼ 0.02 pc.The initially more eccentric setup H7 has a local minimum in the surface density around  ★ ∼ 0.04 pc, inside which the surface density again increases.These milliparsec stars within this region originate from the Hills mechanism.

Figure 22 .
Figure22.The radial mean eccentricity profile of the simulations H1 ( init = 0.1, in blue) and H7b ( init = 0.1, in blue) at  = 6 Myr averaged in bins with a width of Δ = 0.1 pc.The initial inner edge of the disk is indicated with a vertical dashed line as in Fig.21.The initially more circular model H1 has a relatively constant mean eccentricity ∼ 0.4 at every part of the disk.In the model H7b the mean eccentricity has a similar value outside  in , but rapidly increases within  ★ ≲ 0.04 pc.This is again due to very eccentric stars deposited from binary systems into this region by the Hills mechanism.
and O ( ★ > 20 M ⊙ ) as before.The milliparsec population always has considerably more Btype stars than O-type stars.The number of of O-stars is always  O ≲ 12 and decreasing rapidly with time as O-stars die due to single stellar evolution or get tidally disrupted by the SMBH.The last O-type stars in the milliparsec population are dead by  ∼ 6 Myr -7 Myr.We specifically highlight the simulation setup M7b in which all the few initial O-stars ( O = 4) were disrupted by the SMBH before  = 0.5 Myr, and the number of B-type milliparsec stars was  B ∼ 10 until the end of the simulation at  = 10 Myr.A larger simulation sample beyond a single realization of the each initial model would be required to assess how common the early O-star disruption scenario actually is.

Figure 23 .
Figure 23.The number of stars which resided at least 0.1 Myr as a milliparsec star in the simulation setups L, M and H.The total number of stars is show in orange while the number of stars which originate from a Hills mechanism binary break-up is presented with the blue line.With initial disk eccentricities smaller than  init ≲ 0.4 there are no stars originating from the Hills mechanism.With higher initial eccentricities the number of Hills origin stars around the SMBH grows exponentially (dashed black line in the right plot) and is  Hills ∼ 70 at maximum.

Figure 24 .
Figure 24.The number of simulated O-type (O, in yellow), B-type (B, in orange) and late-type (<B, in blue) milliparsec stars as a function of time.The top panel shows the results of the setups M7b (solid lines) and H7b (dashed lines) The number of O-stars never exceeds  O = 12 and reaches zero at  = 6.5 Myr.In the setup M7b all the O-stars are early disrupted by the SMBH.The number of B-stars is close to constant  B ∼ 10 (M7b) and  B ∼ 30 (H7b) while the number of late-type stars is small,  <B ≲ 4. The middle panel shows the runs M 1.3 with a less extreme IMF ( = 1.3).The variant without both PN terms and TDEs (dashed line) has somewhat more stars in the milliparsec population than the runs which include TDEs (dotted line) and both PN and TDEs (solid line).The number of O-stars behaves as in the run H7b while the number of B-type stars is  B ∼ 10-20 in the run variant M 1.3 PNTDE.The number of late-type stars is now larger than the number of B-stars,  <B ∼ 100.The bottom panel shows the results of the setup M 1.7 which are qualitatively similar as in setup M 1.3 .The main difference of the two setups is the somewhat larger number of B-type and late-type stars in the setup M 1.7 due to its less top-heavy IMF.
Figure25.The fraction of milliparsec stars  0.5 =  ( ★ > 0.5)/ with eccentricities higher than  ★ > 0.5 as a function of time in the various simulation setups.The thermal eccentricity distribution  () = 2 (dashed black line) has a fraction of  0.5 = 3/4.Regardless of the disk mass, the IMF and the TDE prescription, the simulation setups with PN equations of motion enabled (solid green, blue and purple lines) have a higher fraction of high-eccentricity stars than the thermal distribution.On the other hand, in the runs without PN (solid orange, yellow and red lines), the fraction of high-eccentricity stars is smaller than for the thermal distribution.This is due to the initially very eccentric orbits of the Hills origin milliparsec stars and the suppression of relaxation effects due to the PN precession, as discussed in the text.

Figure 26 .
Figure 26.The cumulative eccentricity distributions  () of the milliparsec stars in the simulation sample H at  = 6 Myr.Even though most of the stars originate from the Hills mechanism in runs H7 and H9, the milliparsec stars thermalize in  = 6 Myr.Initially more circular setups H1 and H3 have a clearly sub-thermal distribution, and the most circular model H1 does not have stars with  ★ > 0.75.These runs do not include PN equations of motion and relaxation effects are overestimated as discussed in the main text.Simulations including PN effects are presented in a later figure.

Figure 27 .
Figure27.The cumulative eccentricity distribution of  () of the milliparsec stars at  = 6 Myr.The thermal distribution  () =  2 is indicated using the dashed black line.With the PN effects enabled (orange lines), the VRR is suppressed and the eccentricity distributions remain superthermal.The eccentricity distributions are the most concentrated to very high eccentricities ( ★ > 0.9) with the setups including extremely top-heavy IMFs (M7b, H7b).With less extremely top-heavy IMFs (M 1.3 PNTDE and M 1.7 PNTDE), the eccentricity distribution becomes gradually less and less superthermal.Without PN (blue and yellow lines), relaxation proceeds and the resulting eccentricity distributions are subthermal.Whether TDEs are included or not appear to have only a small effect on the overall results.

Figure 28 .
Figure 28.The angular momentum vector directions of the milliparsec stars at  = 6 Myr.The original angular momentum vector directions lie at the top of each panel.The extremely top-heavy models ( = 0.25) in the top row (especially M7b) yield disky, anisotropic milliparsec populations.Models without PN equations of motion are more isotropic than their PN counterparts because of the missing suppression of the relaxation effects.There milliparsec populations become increasingly is isotropic with less top-heavy IMFs (increasing ) of the setups.As further discussed in the text, modeling the spherical background as live particles instead of a smooth external potential would enhance relaxation effects, rendering the milliparsec stars more isotropic.

8. 2
Comparison to observations -the S-stars in the Milky Way center8.2.1 O-type stars and the S-clusterThe innermost massive O/WR stars in the Milky Way center are located outside the S-cluster, orbiting at the sharp inner edge of the present-day CW disk.In our simulations, the milliparsec stars in this region originate both from the Hills mechanism (initially high eccentricity) and from the inner edge of the disk (initially low to moderate eccentricity).Both mechanism bring O-type stars in the young, central milliparsec population predominantly before  = 1 Myr.A number of solutions to the missing O/WR stars in the S-cluster have been proposed.First, the S-cluster could be old enough that the O-stars are already dead.However, this does not explain why the CW disk prominently still contains O-type stars if the S-cluster and the CW disk are indeed coeval.A refined solution is that the O-stars in the S-cluster have already entered the giant phase of their evolution, and have been tidally disrupted by the SMBH.As another related and more complex solution, Chen & Amaro-Seoane (2014) proposed a formation of a so-called rapidly evolving region (RER) near the inner edge of the disk which can drive massive stars migrate close to the SMBH due to ZLK-like resonances and are tidally disrupted.Finally, it has been noted that the Hills mechanism deposits more massive (and thus more Otype) captured companions at larger separations from the SMBH(Generozov & Madigan 2020).

Table 2 .
The approximate estimates for a number of precession and relaxation timescales in the Milky Way center.

Table 3 .
The bifrost user-given parameters relevant for simulations of this study.

Table 4 .
The main properties of the simulation samples L (low-mass, literature), M (mid-mass) and H (high-mass) containing nine simulations each.

Table 5 .
The simulations exploring the effect of the IMF, the binary population properties and PN and TDE prescriptions compared to the setup M from Table4.