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

We study the 10 Myr evolution of parsec-scale stellar discs with initial masses of M_disc = 1.0-$7.5 \times 10^4\, \mathrm{M}_\odot$ and eccentricities e_init = 0.1-0.9 around supermassive black holes (SMBHs). Our disc 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 BIFROST, including post-Newtonian (PN) equations of motion and simplified stellar evolution. All discs are unstable and evolve on Myr time-scales towards similar eccentricity distributions peaking at e_⋆ ~ 0.3-0.4. Models with high e_init also develop a very eccentric (e_⋆ ≳ 0.9) stellar population. For higher disc masses M_disc ≳ 3 × 10⁴ M_⊙, the disc disruption dynamics is more complex than the standard secular eccentric disc instability with opposite precession directions at different disc radii - a precession direction instability. We present an analytical model describing this behaviour. A milliparsec population of N ~ 10-100 stars forms around the SMBH in all models. For low e_init, stars migrate inward while for e_init ≳ 0.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 centre clockwise disc and the S-cluster.