A NEW SECULAR INSTABILITY OF ECCENTRIC STELLAR DISKS AROUND SUPERMASSIVE BLACK HOLES, WITH APPLICATION TO THE GALACTIC CENTER

Most well-observed galaxies with stellar bulges show evidence of supermassive black holes (SMBHs) in their centers (Kormendy 2004). An important feature of stellar dynamics in galactic nuclei is that for orbits contained well within the SMBH's radius of influence, the precession time t_prec is many orders of magnitude larger than their period P.¹ Therefore, for a study of long-term (secular) gravitational dynamics, it is useful to think of the stars not as particles but as massive elliptical wires stretched along the stars' eccentric orbits (Rauch & Tremaine 1996). Gravitational torques between the wires can become the dominant mode of orbital eccentricity evolution, as is seen in the case of resonant relaxation (Rauch & Tremaine 1996; Gürkan & Hopman 2007; Touma et al. 2009; Eilon et al. 2008). Secular effects have a significant impact on the stability of stellar cusps (Tremaine 2005; Polyachenko & Shukhman 1981; Polyachenko et al. 2008), on gravitational wave emission from in-spiraling compact objects (Hopman & Alexander 2006a; Hopman & Alexander 2006b), on tidal disruption of stars (Rauch & Ingalls 1998) and, as we show here, on the origin of the S-stars in the Galactic center.

In this Letter, we describe a mechanism for generating high-eccentricity stars via a fast secular instability of an eccentric stellar disk (Section 2). We use N-body simulations to investigate the nonlinear development of the instability (Sections 3 and 4) and to explore its dependence on the initial eccentricity and mass of the disk. We then show that the instability may be relevant for young stars in the Galactic center, and may provide an interesting channel for formation of the S-stars and of high-velocity stars in the Galactic halo (Section 5). We critically discuss our results in Section 6.

Consider a thin stellar disk of mass M_disk around an SMBH of mass M_• which is embedded in a power-law stellar cusp of mass M_cusp. Our two main assumptions are that (1) M_disk ≪ M_cusp, so that precession of individual orbits is driven by the cusp, and (2) initially the stellar eccentricity vectors

$$\begin{equation} {\bf e} = \frac{1}{G M_{\bullet }} {\bf v} \times ({\bf r} \times {\bf v}) - \hat{\bf r} \end{equation} \tag{ 1 }$$

are aligned and similar in magnitude so that the disk precesses coherently as a whole (for a typical power-law cusp, the orbital precession rate depends only weakly on the semimajor axis; see below). Here, |e| corresponds to the magnitude of the eccentricity of the orbit and r and v are the position and velocity vectors, respectively. Both assumptions are realistic for a disk of young stars formed in a galactic center via the gravitational capture of part of an infalling molecular cloud: observations reveal a massive spherical stellar cusp (Genzel et al. 2003; Schödel et al. 2007), while hydrodynamical simulations tentatively show similar initial orbital eccentricities of young stars² (Sanders 1998; Alexander et al. 2008; Bonnell & Rice 2008; Hobbs & Nayakshin 2009).

The orbits precess in the direction opposite to the orbital rotation of the stars. The precession timescale scales as

$$\begin{eqnarray} t_{\rm prec} & \sim \frac{M_{\bullet }}{N({<}a) m} P(a) f(e)\nonumber\\ & \propto a^{\gamma -3/2} f(e), \end{eqnarray} \tag{ 2 }$$

where $P(a) = 2 \pi \sqrt{a^3/GM_{\bullet }}$ is the period of a star with semimajor axis a, N(<a) is the number of stars in the cusp within a, m is the individual mass of the stars, and γ is the power-law index for the space density ρ(r) ∝ r^−γ. For a power-law cusp with γ = 3/2 (Young 1980), t_prec is constant for all a.

Observations of the Galactic center indicate that there is a stellar cusp in the inner parsec, with published measurements of the slope not significantly different from γ = 3/2, e.g., γ = 1.4 ±  0.1 (Genzel et al. 2003) and γ = 1.19 ±  0.05 (Schödel et al. 2007). It is important to stress that these observations show only a small and biased subset of the stellar content in that region. They are restricted to massive young stars, that have not yet relaxed, and to giants, which have relaxed, but are affected by hydrodynamical collisions with other stars (e.g., Alexander 1999). These two stellar types are not expected to dominate the density. Instead, theoretical models accounting for mass segregation show that the density at the regions of interest is mostly determined by stellar black holes and main-sequence stars, with cusp values between γ_MS = 1.4 and γ_BH = 2.0 (Alexander & Hopman 2008). We conclude that both theory and observations yield values of γ close to 3/2 and the dependence of t_prec on a is weak.

A key element of the eccentric instability is that f(e) in Equation (2) is an increasing function of eccentricity (e.g., Binney & Tremaine 2008), so that an orbit that is slightly more eccentric than average will lag behind in precession. Such a star will feel a strong, coherent torque from the other stars in the disk, in the opposite direction of its angular momentum vector. As a result, its angular momentum decreases in magnitude, causing its eccentricity to increase even further. In this way, very high eccentricities can potentially be achieved. Conversely, if a stellar orbit is slightly less eccentric than the average, it will have a higher precession rate, thus experiencing a torque which acts to decrease its eccentricity further.

We thus find that the initial eccentricity vector distribution is unstable, due to a combination of retrograde precession and df/de > 0, and that all eccentricities are driven away from the initial eccentricity of the disk. It can be shown analytically that the initial growth timescale of the eccentric instability scales as

$$\begin{equation} \tau _{\rm growth} \sim t_{\rm prec} \left(M_{\rm cusp}\over M_{\rm disk}\right)^{1/2}\left[e\sqrt{1-e^2} \frac{d}{de} \left(\frac{1}{f(e)}\right)\right]^{-1/2}. \end{equation} \tag{ 3 }$$

For the Galactic center, this amounts to several precession timescales.

We now demonstrate the eccentricity instability using N-body simulations and look at its impact on the young stars in the Galactic center.

We summarize a newly developed code which we have used in these simulations. Our integrator is based on a fourth-order Wisdom–Holman algorithm (Wisdom & Holman 1991; Yoshida 1990), where the forces on a particle are split into a dominant Keplerian force and perturbation forces. We use Kepler's equation to integrate the orbit (Danby 1992), keeping the position of the central object fixed. As we use adaptive time stepping to deal with close encounters, the algorithm loses its symplecticity and is not strictly energy conserving. We implement time symmetry to reduce the energy error and the K₂ kernel from Dehnen (2001) for gravitational softening. The overall fractional energy error for these simulations is typically ∼10⁻⁸ or less. We describe the code in detail in an upcoming paper.

3.1. Initial Conditions

We have three main components in our simulations. (1) An eccentric disk with surface density profile Σ ∝ r⁻² (compare Paumard et al. 2006), consisting of 100 equal-mass stars, with semimajor axes 0.05 ⩽ a ⩽ 0.5 . We vary the masses of the stars for different simulations. We initialize the eccentricity vectors of all stars pointing in roughly the same direction, though they are scattered in their orbital phase. Initial inclinations are distributed between h/r ∼ ±0.01. (2) A SMBH of 4 × 10⁶ M_☉ (Ghez et al. 2008), and (3) a smooth stellar cusp with density profile index γ = 1.5 and a mass of 0.5 × 10⁶ M_☉ within 1 pc (Schödel et al. 2007; Genzel et al. 2003). The stars within the disk precess due to the cusp and experience a prograde general relativistic precession which amounts to apsidal precession of the orbit by an angle

$$\begin{equation} \delta \phi _{\rm prec} = \frac{6 \pi GM_{\bullet }}{ a (1 - e^2) c^2} \end{equation} \tag{ 4 }$$

per orbit. Precession due to self-gravity of the disk is minimal.

To verify our results, we also run simulations with twice the number of half-as-massive stars and various softening parameters. We discuss the use of different values of γ (e.g., Bahcall & Wolf 1976) in Section 6.

With a mass of 8 × 10³ M_☉ and an initial eccentricity of 0.6, we let the disk of stars evolve for several precession times (t_prec∼ 0.6–0.7 Myr or ∼1200 orbits at a = 0.05  ). In Figure 1, we show a selection of stars from this run which best demonstrates the dispersion of the disk's initial eccentricity, while Figure 2 follows the evolution of the eccentricity vectors (this time for M_disk = 4 × 10³ M_☉). The circles in the latter plot indicate the stars experiencing the highest torque (τ) and hence angular momentum changes. Note that the highest eccentricity orbits are precessing behind the bulk of the stars as expected. We find a monotonically decreasing relationship between the rate of change of angular momentum, i.e., torque on each orbit, and semimajor axis a.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effects of the instability on the stellar eccentricity distribution are evident after ∼0.5 − 1t_prec and the strongest torques are experienced over a few precession times before the disk loses coherence.

We explore the dependence of the final range of eccentricities on the mass and initial eccentricity of the disk. To this end, we run a set of simulations with a disk of initial eccentricity 0.6 for various disk masses (4, 6, 8, 10 ×10³ M_☉), and a set with M_disk = 10⁴ M_☉ for various initial eccentricities (0.1, 0.4, 0.5, 0.6, 0.7, 0.9). The results are presented in Table 1. Higher initial eccentricities and masses of the disk are correlated with higher resulting eccentricities. We also see a correlation between the highest eccentricities and inclinations reached, with the largest values occurring at the smallest radii.

In contrast to the secular instability described in Tremaine (2005), relativistic effects in this case do not promote stability once the relativistic prograde precession rate exceeds that of the retrograde precession due to the cusp. The highest eccentricity orbit precesses with prograde motion at a higher rate than the average and again experiences a torque in the opposite direction of its angular momentum.

The region surrounding the SMBH in the Galactic center (Gillessen et al. 2009) is host to several interesting populations of young stars, which include one well-defined clockwise disk of O and Wolf–Rayet (WR) stars, and a group of counter-clockwise moving O and WR stars, possibly a dissolving second disk (Levin & Beloborodov 2003; Nayakshin & Cuadra 2005; Genzel et al. 2003; Paumard et al. 2006; Lu et al. 2009; Bartko et al. 2008). These O and WR stars are located at projected distances of 0.05–0.5 pc, and were most likely formed as a result of gravitational fragmentation of a gaseous disk (Levin & Beloborodov 2003).

More puzzling in terms of origin is the population of B-stars known as the "S-stars" (Schödel et al. 2002; Ghez et al. 2005; Eisenhauer et al. 2005; Gillessen et al. 2009). At distances of 0.003 − 0.03  from the SMBH, in situ formation seems highly improbable due to the tidal field of SgrA* (Levin 2007). In addition, their young age (∼20 − 100 Myr) imposes a tight constraint on formation scenarios in that they cannot have traveled far from their place of birth. One possibility is the formation of the S-stars from the disruption of stellar binaries (Hills 1988; Gould & Quillen 2003). Perets et al. (2007) propose the surrounding stellar bulge as the source of the binaries, while Löckmann et al. (2008, LBK) invoke two stellar disks in the Galactic center to explain the origin. In their simulations, the stars experienced Kozai-type torques which induced strong changes in their inclination and eccentricity. The eccentricity of many of the stars became greater than 0.9, and LBK argued that if these stars were binaries they would be tidally disrupted by the SMBH, thus producing the S-stars and high-velocity stars via Hills' (1988) mechanism.

However, the gravitational influence of a stellar cusp, which was neglected by LBK, strongly suppresses Kozai-type dynamics. This was recently discussed in Chang (2009). Chang's treatment was highly simplified; in particular, his expressions for Kozai evolution were valid only if the stellar orbits were inside the inner edges of both disks. We have performed LBK-type simulations and have confirmed Chang's general conclusion that a realistically massive stellar cusp entirely suppresses the production of high-eccentricity stars via a Kozai-type mechanism.

In contrast, the eccentricity instability is capable of driving several stars to near-radial orbits, as is evident in Figure 1. For capture, the pericenter r_p of the binary's orbit must come within the tidal radius r_t = (2M_•/m_bin)^1/3a_bin. This results in a captured star with semimajor axis a_• that scales as a_• ∼ (M_•/m_bin)^2/3a_bin, where a_bin is the semimajor axis of the binary itself, m_bin is its combined mass and the eccentricity of the captured star can be approximated as 1 − e ∼ (m_bin/M_•)^1/3.

We summarize our results for varying masses and initial eccentricities of the disk in Table 1. All simulations have evolved for ∼ 7 Myr. For a given semimajor axis of a binary system, we calculate the tidal radius at which the binary will be disrupted and finally the percentage of stars that come within this radius. It is clear that for many parameters of our disk, stars have r_p < r_t. For example, a disk with an initial eccentricity of 0.6 and mass of 10⁴ M_☉, 14% of stars in the disk will enter the tidal radius. Thus, if some of these stars are binaries, Hills' mechanism for producing the S- and high-velocity stars (Brown et al. 2006) could be at work. An important element in this scenario is the binary survival rate at the radii of the disk; however, Perets (2009) show that for a = 0.1 AU the evaporation rate is >10⁷ yr.

The scalar resonant relaxation time in the Galactic center is everywhere larger than the age of the disks. It follows that since the postcapture eccentricity of the stars is of order 1 − e ≈ 10⁻², most of the S-stars cannot have formed in the current disk. The lifetime of many S-stars could be of the order of 100 Myr, and thus earlier (∼10–100 Myr) starburst episodes could have contributed to their present population. In that case, the scalar relaxation time is short enough to randomize the eccentricities (Hopman & Alexander 2006a; Perets et al. 2009). We note that the scenario presented by LBK has the same limitation.

5.1. Bimodal Eccentricity Distribution in Disk

In Figure 3, we compare the observed stellar eccentricity distribution of the clockwise disk (Bartko et al. 2008) with the simulated postinstability eccentricity distribution for all stars with an initial disk eccentricity of 0.6 (compare Bonnell & Rice 2008). We draw attention to the qualitative similarities in the two plots, in particular the double-peaked profile. This bimodality is a natural consequence of the eccentric disk instability as all eccentricities are driven away from the mean. Without the instability it is unclear how to generate such a distribution.

Zoom In Zoom Out Reset image size

So far, we have described the eccentricity instability for an idealized initial disk configuration. We now demonstrate the effectiveness of the instability for a broader scope of parameters and extend our disk model to include (1) a significant range in both magnitude and direction of the stellar eccentricity vectors, (2) heavily inclined stellar orbits, and (3) different values of γ, the power-law index of the cusp, which determines the precession rate within the disk. We discuss each of these in turn.

• 1.  
  We simulate disks with initial eccentricity vectors scattered between a range of opening angles θ ∈ [0, 2π]. As anticipated, significant spreading of the vector decreases the effectiveness of the instability. However, we find sufficiently high-eccentricity orbits to generate S-stars up to θ ∼ π rad, which suggests that the instability remains important for substantially less favorable circumstances. Next, we vary the magnitudes of the eccentricities throughout the disk, basing our initial conditions on Runs 1 and 2 of Bonnell & Rice (2008). Run 1 (0.6 < e < 0.7) is in excellent agreement with our previous results. The instability in Run 2 (0 < e < 0.6), however, is not efficient at generating the highest e orbits and hence the S-stars. These conclusions are transparent; the stars which experience the greatest torque are those at small a (see Section 3) and if initially circular, need a greater torque to produce the highest e orbits. We find that the most important condition, for M_disk = 10⁴, is for innermost orbits to have e > 0.4. Future work on gas dynamics of eccentric disks will tell whether this presents a serious problem for our scenario.
• 2.  
  We study disks with stellar inclinations h/r in the range [ − 0.3, 0.3] and find that the instability is robust in all cases. However, above the limit [ − 0.2, 0.2], the instability appears ineffective in generating the highest e (>1–10⁻³) orbits.
• 3.  
  We vary the power-law slope of the cusp. For γ < 3/2, due to the shallowness of the cusp, the innermost orbits lag behind in precession and are more susceptible to being pushed to high e orbits. While this improves the ability of the instability to produce S-stars, the instability is not as effective at generating low e orbits. Conversely, an increase in γ (>3/2) results in the innermost orbits precessing faster than the bulk of the stars and hence suppresses the generation of the highest e orbits and lengthens the timescale on which it occurs. We conclude that the production of S-stars is unlikely for a cusp with γ>7/4. In all cases examined (1 < γ < 2), the instability proves robust and, as emphasized in Section 2, there is both observational and theoretical evidence that the slope of the cusp in the Galactic center is sufficiently close to 3/2 for the instability to be highly effective.

The eccentricity instability described in this work is relevant for an eccentric disk embedded in a stellar cusp. We note that the eccentric disk in M31 (Tremaine 1995) is not subject to our instability, since it is not embedded in any visible cusp. Indeed, we find that such a disk remains coherent over several precession times, although the innermost stars can undergo significant angular momentum changes as described by Cuadra et al. (2008). The situation is different in the Galactic center, where massive stars were likely born from an eccentric gaseous disk, and where there is strong observational evidence for a heavy spherical stellar cusp (Genzel et al. 2003; Schödel et al. 2007). We have shown that the instability may have strongly influenced the development of the eccentricities of the massive stars observed, and in particular, can account for the bimodality in their eccentricity distribution. Furthermore, combined with Hills' (1988) mechanism, the eccentricity instability may have played a role in the formation of the S-stars and of the high-velocity stars in the Galactic halo.

We thank Atakan Gürkan, Richard Alexander, Simon Portegies Zwart, Mher Kazandjian, Jihad Touma, and Sergei Nayakshin for useful discussions, and Mark Hayden for comments on the text. A.M. is supported by a TopTalent fellowship from the Netherlands Organization for Scientific Research (NWO), Y.L. by a Vidi fellowship from NWO, and C.H. by a Veni fellowship from NWO.