Simulation of the loss-cone instability in spherical systems. II. Dominating Keplerian potential

In the pioneer paper, Polyachenko (1991) studied a simple analytical model of a low-mass stellar disc in a dominating point-mass potential. He argued that if the radial-orbit instability (for review, see Polyachenko & Shukhman, 2015) is suppressed, a distant relative of the loss-cone instability in plasma (Rosenbluth & Post, 1965) still may occur.

Tremaine (2005) examined spherical and disc systems and found that stability properties are determined by the dependence of the distribution function (DF) $F$ on angular momentum $L$. In particular, flattened, nonrotating systems were found to be secularly stable if $\partial F/\partial L<0$ at constant energy, while for $\partial F/\partial L>0$ they are generally unstable. Based on the analysis of dipole and quadrupole distortions, the spherical systems in which $F=0$ at $L=0$ (an empty loss cone) were claimed to be stable for any sign of $\partial F/\partial L$.

Using a newly elaborated matrix method for spherical systems analogous to the matrix method for discs (Polyachenko, 2005), Polyachenko et al. (2007) found that the instability in spherical systems considered by Tremaine is nevertheless possible, but for the higher-order distortions starting from octupole. In the subsequent paper (Polyachenko et al., 2008, hereafter Paper I), we refined our claim that it is only relevant to non-monotonic dependence of the DF on $L$. In analogy with plasma, we called it gravitational loss-cone instability (gLCI).

These analytical studies assume that a cluster of mass $M_{*}$ is embedded in the dominating potential of a central point mass $M$, so that $\varepsilon\equiv M_{*}/M$ is a small parameter. The gravitational potential

consists of the dominant Keplerian part, and $\Phi_{*}$ due to selfgravity of the cluster. If the latter is neglected, a particle travels along the fixed ellipses with semiaxes $a=GM/(2|E|)$, $b=L/(2|E|)^{1/2}$ and eccentricity $e=(1-\alpha^{2})^{1/2}$, where $\alpha\equiv L/L_{\rm circ}$, $L_{\rm circ}=GM/(2|E|)^{1/2}$ is the angular momentum of a particle with energy $E$ on a circular orbit. The radial and orbital frequencies are equal and independent of $L$: $\Omega_{1}=\Omega_{2}=(2|E|)^{3/2}/(GM)$.

Accounting for the selfgravity leads to slow orbit precession, which is always retrograde (Polyachenko et al., 2007). Generally, there are four characteristic timescales in such a system (Tremaine, 2005): dynamical time $t_{\rm dyn}\sim\Omega_{1}^{-1}$, precession time $t_{\rm pr}\sim\varepsilon^{-1}t_{\rm dyn}$, resonant relaxation time $t_{\rm res}\sim N\varepsilon^{-1}t_{\rm dyn}$ (Rauch & Tremaine, 1996), where $N$ is the number of particles, and two-body relaxation time $t_{\rm relax}\sim N\varepsilon^{-2}t_{\rm dyn}$. For $N\gg 1$ these timescales are well separated:

If the system is dynamically stable, its evolution is driven by relaxation on timescales $t_{\rm res}$ or $t_{\rm relax}$. However, if there is an instability due to collective mechanism in the distribution of orbits, its timescale $t_{\rm ins}$ is expected to be close to the orbital movement time, i.e. $t_{\rm pr}$. As found for several models in Paper I, the instability time is in fact only a few times larger than the precession time. Thus instability becomes the main change driving process overtaking both kinds of relaxation. However, the instability is slow in the dynamical timescale and requires by a factor of $\varepsilon^{-1}$ longer simulations than standard simulations of systems in which the dynamical time is largely determined by selfgravity.

This paper is aimed to simulate numerically one of the models studied earlier theoretically using linear perturbation theory (Paper I). It serves two purposes: to confirm our theoretical finding and to go beyond the linear approximation and study the consequences of the instability development. In a companion paper (Polyachenko et al., 2019), we simulated the instability in a dominating harmonic potential and showed that unlike, e.g., the well-known bar instability in discs, gLCI did not change the shape of the cluster, but led to the enhanced diffusion resulting in the repopulation of the loss cone and cease of the instability.

In addition to the peculiarities of the harmonic case, the dominating central point mass requires additional care for particles close to the centre. It is worthwhile to emphasize that popular N-body codes incorporating the black hole dynamics or the external central mass potential (e.g., Nbody6++, , Wang et al.2015) are proved to be effective in the opposite limit $M\ll M_{*}$, but not in the case of interest here. Thus additional investigations were needed to improve the available numerical schemes for the current task.

The paper is organized as follows. The model and the set-up are described in Section 2, along with details of N-body simulations. Results containing a comparison of the instability growth rate with the linear perturbation theory, and analysis of peculiarities of the cluster evolution are given in Section 3. The final section 4 contains conclusions, final remarks, and describes future perspectives.

In order to compare results of our numerical simulations with theoretical predictions, we evaluate the spherical harmonics of the density distribution

where $Y_{l}^{m}(\theta,\phi)$ are fully normalised spherical functions, for each snapshot. The coefficients $C_{l}^{m}$ depend on radius and the orientation of the frame. In order to make a comparison simpler, we evaluate amplitudes

invariant to the frame rotation (see Appendix A and Polyachenko et al., 2019).

From Sec. 1 and Sec. 2 we infer that the dynamical time scales with $M$ as $t_{\rm dyn}\sim M^{-1/2}$. Therefore, the instability time scales as $t_{\rm ins}\sim t_{\rm pr}\sim M^{1/2}$. Fig. 3 compares the amplitudes of the octupole perturbations ($l=3$) for two runs in the instability time scale. We see the exponential growth of the amplitudes with the predicted slope, however the amplitudes are modulated by oscillations. These oscillations are due to i) the oscillation character of the eigenfrequency $\omega_{r}\neq 0$, and ii) the existense of another eigenfrequency with the same growth rate. Indeed, the integral equation (2.18) of Paper I has the form $\omega^{2}\psi_{l}={\hat{L}}_{l}\,\psi_{l},$ where ${\hat{L}}_{l}$ is the linear integral real operator. Hence, if $\omega^{2}$ is an eigen-value, then $(\omega^{2})^{*}$ is also tan eigen-value of the problem. This means that if $\omega_{r}+i\gamma$ is a solution, then $-\omega_{r}+i\gamma$ also is a solution. A general growing solution then reads

and

where $\phi$ is due to the difference of the complex phases of the amplitudes $C_{1,2}$. These oscillations we see in Fig. 3.

We also observe early saturation of the instability so that amplitudes remain below a level of 0.03, similar to the gLCI in a dominating harmonic potential (Polyachenko et al., 2019). There, it was explained by rapid repopulation of the loss cone, i.e. the region of low angular momentum, due to diffusion enhanced by the instability.

For near-Keplerian models, we estimate the relaxation time using the standard formula for two-body relaxation

so that in dimensionless form the two-body relaxation time in units of the slow time $M^{1/2}$ is

where the value range covers the parameters shown Fig. 3.

For the resonant relaxation (Rauch & Tremaine, 1996; Tremaine, 2005)

so that in dimensionless form the resonance relaxation time in units of slow time $M^{1/2}$ is

One can conclude that these two mechanisms alone cannot be responsible for a particle repopulation during the time of simulations.

To quantify the speed of the diffusion in $E$ and $L$, we evaluated the spread $\sigma_{t}[E]$ as the difference between the third and the first quartiles of the distribution of energy of individual particles $E_{i}(t)$. This method is more robust to outliers than the calculation of the standard deviation $\sigma$, while it gives $1.35\,\sigma$ for a normal distribution. Similarly, $\sigma_{t}[L]$ was obtained for the angular momentum shifts $L_{i}(t)-L_{i}(0)$. In the diffusion driven by two-body relaxation, one would expect

where $\zeta$, $\eta$ are constants of order unity, $W$ is some energy characteristic of the system, e.g. total kinetic energy, or virial of Clausius $W_{\rm C}$ (e.g., Sellwood, 2015). If the resonant relaxation takes place for angular momentum, one would expect $\eta\sim(\varepsilon\ln N)^{-1}$.

The character of the diffusion can be understood by analyzing the plots summarised in Fig. 4. The upper left panel shows the normalized energy diffusion coefficient

By construction, it should give the constant value $\zeta$ for two-body relaxation. In fact, we observe an almost constant value of order unity, which indicates that spreading in energy is due to two-body relaxation. In the upper right panel, we present the energy relaxation time calculated as

for nine available runs in which the central mass $M$ varied from 3 to 100, and the number of particles $N$ varied from $25\cdot 10^{3}$ to $10^{6}$. Dashed lines show fits

for various $M$, demonstrating that the scaling perfectly corresponds to the two-body relaxation time.

The profiles for angular momentum diffusion consist of two distinctive parts: nearly constant in the beginning and modulated exponential growth and fall later. Using the first constant part, we evaluated the angular momentum relaxation time

which is presented in the lower right panel of Fig. 4. Color dashed lines show fits

which is close to the predicted scailing of resonant relaxation $t_{\rm res}\propto M^{1/2}N$ (Rauch & Tremaine, 1996), shown by the black dashed line.

Using the dependence of Eq. (3.15), we construct the normalized angular momentum diffusion coefficient

The profiles are shown in the lower left panel of Fig. 4. Due to normalization, the constant parts are close to unity. The growing parts are obviously due to instabilty, as follows from a comparison with the theoretical slope for the octupole harmonic given in Fig. 3.

The gravitational loss-cone instability (gLCI) in a dominating Keplerian (near-K) potential predicted earlier theoretically (Polyachenko et al., 2007, 2008), is now revealed for the first time in numerical simulations. In the limit of slow precessing orbital motion, it proved possible to connect the sign of the orbital precession rate with the derivative of the DF with respect to the angular momentum. It was found that gLCI is possible if $\partial F/\partial L\cdot\Omega_{\rm pr}<0$. For the dominating Keplerian potential, the precession is always retrograde ($\Omega_{\rm pr}<0$), i.e., independently of the cluster DF, so gLCI requires $\partial F/\partial L>0$ at small $L$, or a deficit of particles at small $L$ (the loss cone) must be present. This instability has a well-known counterpart in plasma physics called loss-cone instability (Rosenbluth & Post, 1965).

Our N-body runs typically consist of $N\sim 10^{5}$ particles evaluated during $\lesssim 10^{4}$ dynamical times. The main issue is the dominating central mass which requires special treatment of the orbital integration near the center. The simulations reproduce well the linear regime for which instability of the octupole harmonic was predicted. In particular, we obtained the growth rate of the harmonic amplitude and periodic modulations of the amplitudes due to the presence of the complementary unstable eigen-mode (i.e. the mode with the eigen-frequency $-\omega^{*}$). However, the shape of the cluster does not change due to early instability saturation. The latter occurs due to diffusion enhanced by the instability.

A comparison to the near-harmonic (near-h) case (Polyachenko et al., 2019) shows the following distinctive features:

In the future, we plan to explore a diversity of more realistic DFs to detect the gLCI. We also plan to modify the numerical code so that the stars diffused to the loss cone will be excluded from the simulations. This may help to sustain the instability.

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject A06), by the Volkswagen Foundation under the Trilateral Partnerships grants No. 90411, 97778, and the Basic Research program II.16 (Ilia Shukhman). Peter Berczik acknowledges support by the Chinese Academy of Sciences through the Silk Road Project at NAOC, through the “Qianren” special foreign experts program, and the President’s International Fellowship for Visiting Scientists program of CAS, the National Science Foundation of China under grant No. 11673032 and also the Strategic Priority Research Program (Pilot B) “Multi-wavelength gravitational wave universe” of the Chinese Academy of Sciences (No. XDB23040100). The special GPU accelerated supercomputer Laohu at NAOC has been used and we thank the Center of Information and Computing of NAOC for support. Peter Berczik also acknowledges the special support by the NASU under the Main Astronomical Observatory GRID/GPU computing cluster project. This work benefited from support by the International Space Science Institute, Bern, Switzerland, through its International Team program ref. no. 393 “The Evolution of Rich Stellar Populations & BH Binaries” (2017-18).