Dynamical coupling of Keplerian orbits in a hierarchical four-body system: from the Galactic Centre to compact planetary systems

ABSTRACT

This study focuses on the long-term evolution of two bodies in nearby initially coplanar orbits around a central dominant body perturbed by a fourth body on a distant Keplerian orbit. Our previous works that considered this setup enforced circular orbits by adding a spherical potential of extended mass, which dampens Kozai–Lidov oscillations; it led to two qualitatively different modes of the evolution of the nearby orbits. In one scenario, their mutual interaction exceeds the effect of differential precession caused by a perturbing body. This results in a long-term coherent evolution, with nearly coplanar orbits experiencing only small oscillations of inclination. We extend the previous work by (i) considering post-Newtonian corrections to the gravity of the central body, either instead of or in addition to the potential of extended mass, (ii) relaxing the requirement of strictly circular orbits, and (iii) removing the strict requirement of complete Kozai–Lidov damping. Thus, we identify the modes of interorbital interaction described for the zero eccentricity case in the more general situation, which allows for its applicability to a much broader range of astrophysical systems than considered initially. In this work, we scale the systems to the orbits of S-stars; we consider the clockwise disc to represent the perturbing body, with post-Newtonian corrections to the gravity of Sagittarius A* playing the role of damping potential. Considering post-Newtonian corrections, even stellar-mass central bodies in compact planetary systems can allow for the coupled evolution of Keplerian orbits.

1 INTRODUCTION

The study of dynamics in Keplerian potentials is an old yet very progressive area of research. The secular orbital evolution of light (test) particles in the dominating central potential accompanied by a distant perturber is one of the classical problems in celestial mechanics. According to the pioneering works of Kozai (1962) and Lidov (1962), the orbital solution within this hierarchical three-body setup is often called Kozai–Lidov (K-L) dynamics. Various works have extended its original formulation, which supposed a non-evolving circular orbit of the perturber, for example, eccentric perturber (Naoz et al. 2011; Lithwick & Naoz 2011), relativistic effects (Naoz et al. 2013; Lim & Rodriguez 2020), mass loss, and transfer (Michaely & Perets 2014).

Considering the four-body setup brings new degrees of freedom and also more variants of the general setup (Huang 1960; Simó 1978; Scheeres 1998; Baltagiannis & Papadakis 2011). A possible configuration has recently been investigated by Haas, Subr & Vokrouhlicky (2011b). Similarly to K-L dynamics, their setup consists of a dominating central body and a massive perturber on a circular orbit. Contrary to K-L dynamics, they considered the orbital evolution of two light, mutually gravitationally interacting bodies inner to the orbit of the massive perturber. In their work, Haas et al. (2011b) focused on the case when the two inner orbits are close to each other in terms of semimajor axes and are initially coplanar (with arbitrary inclination with respect to the orbit of the perturber). An additional assumption, primarily imposed due to limitations of the used calculus, was the non-evolving zero eccentricity of the two inner orbits. Haas et al. (2011b) argue that this assumption is relevant if another non-Keplerian spherically symmetric potential is present within the system, being strong enough to damp the K-L oscillations of the inner bodies enforced by the massive perturber. Within this setup, Haas et al. (2011b) developed a secular theory showing that the two inner orbits periodically exchange their angular momentum such that their inclinations oscillate. If their mutual interaction is strong enough (which depends on their mass and separation), the precession of their orbits is synchronized, that is, the initial coplanar structure is nearly preserved. In the other case, orbital planes of the inner bodies precess differentially due to the perturbing force of the outer body, which leads to disruption of the coplanar configuration. We refer to the temporal evolution of the specific four-body setup introduced by Haas et al. (2011b) as the VHS mechanism (based on the authors of the Haas et al. (2011b) paper, D. Vokrouhlický, J. Haas and L. Šubr) throughout this paper.

A shortcoming of the secular theory of VHS dynamics is the requirement of the spherically symmetric external potential needed to dampen the K-L oscillations, which reduces its applicability in observed astrophysical systems. However, Haas, Šubr & Kroupa (2011a), Haas et al. (2011b) introduced a physically realistic setup in which the VHS mechanism is applicable. They studied a system in which the supermassive black hole (SMBH) in the Galaxy’s centre, Saggitarius A* (SgrA*, Ghez et al. 2003; Eisenhauer et al. 2005; Gillessen et al. 2009a, b; Yelda et al. 2010), represents the dominant body, and the additional spherical potential is due to the surrounding nuclear star cluster. They considered the perturbing body to be the circumnuclear gaseous disc (Martín et al. 2012; Liu et al. 2012; Hsieh et al. 2017; Tsuboi et al. 2018; Goicoechea et al. 2018; Hsieh et al. 2021) and the bodies on inner nearby coplanar orbits to be the observed stars from the young stellar disc that is within the distances of 0.04–|$0.4\, \mathrm{pc}$| from the central SMBH (Levin & Beloborodov 2003; Paumard et al. 2006; Bartko et al. 2009, 2010). Haas et al. (2011a, b) suggested that the four-body dynamics in the spherically symmetric external potential can explain the specific, near-perpendicular orientation of the stellar disc with respect to the distant perturber.

Our study aims to expand the scope of the VHS dynamics described in Haas et al. (2011b) and to explore its applicability in a broader range of astrophysical systems by relaxing some of the assumptions of the underlying secular theory. First, we develop the idea, suggested in the original work, that the non-Keplerian spherical potential can be omitted if we consider post-Newtonian terms in the gravity of the central body while still working within the secular approach. Secondly, we investigate the evolution of systems with a non-zero eccentricity of the two inner orbits by directly integrating the equations of motion. Finally, we consider a scenario in which the eccentricity of the inner orbits evolves over time, that is, when the K-L oscillations are not entirely damped.

The paper is structured as follows: in Section 2, we provide a detailed description of the four-body setup we are studying. Section 3 provides a summary of the secular theory developed in Haas et al. (2011b), along with a discussion of the damping of K-L oscillations due to the effects of general relativity. Section 4 describes several examples of systems with non-zero eccentricity that were integrated. Finally, we present our conclusions on the generalized VHS dynamics in Section 5.

2 SETUP

We study a hierarchical four-body system with a dominant central body, characterized only by its mass, M_•. The system further consists of a distant perturber of mass M_p on a circular orbit with radius R_p around the central body. The orbit of the perturber defines the reference plane. We can choose any line within this plane to define our reference axis to calculate the longitude of the ascending node, Ω. Finally, we consider two light particles of masses m and m′ where m, m′ ≪ M_p on orbits around the central body with a semimajor axes a and a′, which are much smaller than R_p and a′ < a. These light bodies are in inclined orbits, having inclinations i and i′ with respect to the reference plane. The last important parameters we consider are the longitudes of the ascending nodes of the two bodies, Ω and Ω′. Initial conditions are set up such that i = i′ and Ω = Ω′.

As an example, in this work we use the objects observed in the Galactic Centre as an astrophysical system to provide us realistic values for M_•, M_p, and R_p. We set the system that may correspond to the situation in the vicinity of the SgrA* black hole, that is, M_• = 4 × 10⁶ M_⊙ (Ghez et al. 2003; Eisenhauer et al. 2005; Gillessen et al. 2009a, b; Yelda et al. 2010). We consider the distant perturber of mass of M_p = 10⁴ M_⊙ and a semimajor axis of R_p = 0.1 pc, aiming to mimic the overall gravitational influence of the observed clockwise young stellar disc (CWD, Lu et al. 2013; von Fellenberg et al. 2022). The two light particles could be representatives of the S-stars that are observed in the Galactic Centre.

3 SECULAR THEORY

In this section, we follow the mathematical approach used in Haas et al. (2011b) and briefly sketch the main ideas. In particular, we consider a secular approach to describe the long-term evolution of the system described in Section 2. For this, the mean interaction potential of the system, |$\overline{\mathcal {R}}$|⁠, averaged over fast changing variables needs to be specified. As it can be given as a direct sum of the individual terms describing different components of the system, we discuss these separately in the following sections.

3.1 Potential of the distant/outer perturber

The averaged interaction potential between a perturbing body on a circular orbit with radius R_p and a particle on an orbit with semimajor axis a, eccentricity e, and inclination i with respect to the orbital plane of the perturber reads (Kozai 1962):

where ω is the argument of periapses of the orbit. Suppose |$\overline{\mathcal {R}}_\mathrm{p}$| is the only component of the total perturbing potential (i.e. the system is reduced to a three-body setup). In that case, the body on the inner orbit is subject to quadrupole K-L dynamics (Kozai 1962; Lidov 1962). Depending on the initial conditions, its eccentricity and inclination may undergo large periodic variations that are mutually coupled through the so-called Kozai integral, |$C \equiv \sqrt{1 - e^2}\cos {i}$|⁠, which, together with the semimajor axis (a) and |$\overline{\mathcal {R}}_\mathrm{p}$|⁠, is a conserved quantity along the orbit evolution.

The number of known integrals of motion allows for an effective insight into the K-L dynamics through plots of isocontours of |$\overline{\mathcal {R}}$| in the e–ω space, which for fixed values of a and C give sets of possible evolutionary tracks (see Fig. 1). These sets form two qualitatively different topologies: for |$C \gt \sqrt{3/5}$|⁠, they consist of concentric ovals, which means that the eccentricity oscillates slightly along the evolutionary path and ω rotates within the whole range (0, π) (see right panel of Fig. 1). If |$C \le \sqrt{3/5}$|⁠, the topology qualitatively changes: a separatrix crosses the central point; It divides the diagram into zones with ω librating around the value of π/2 or 3π/2 and the outer rotation zone (left and middle panels of Fig. 1). The lower the value of C, the larger the eccentricity oscillations. The characteristic time scale for these oscillations is given by (Kozai 1962; Lidov 1962):

An important result from the isocontour plots is that the zero eccentric orbit is stable for |$C \gt \sqrt{3/5}$|⁠, while it undergoes periodic variations when C below the limiting value.

Finally, let us note that the longitude of the ascending node, Ω, rotates monotonically in the full range of (0, 2π) for arbitrary initial conditions. The rate of precession depends on the other orbital elements, as well as on the mass and semimajor axis of the perturber. However, the value of Ω does not affect the evolution of the other orbital elements, which is a natural consequence of the axial symmetry of the problem.

3.2 Spherical potential

In our study of four-body systems, we examine two distinct sources of an external spherical potential. The first source is the presence of an extended mass around the central body, while the second source is an approximation of the first-order post-Newtonian corrections to the gravity of M_•. Although these sources differ, they have very similar effects and impact the evolution of the two bodies in a similar manner.

3.2.1 Extended mass

Haas et al. (2011a, b) considered such an astrophysical context involving an extended mass around the central body, influencing the secular dynamics of the inner orbit(s). In particular, the authors provide an analytic form for the mean potential corresponding to the mass density distribution with power-law profile, ρ_c ∝ r^{β − 2},

where M_c stands for the integral of the extended mass density within the orbit of the perturber (R_p) and

The coefficients a_n are given by

with a₁ = β(1 + β)/4.

From the spherical symmetry of this perturbing potential we get that its only manifestation on the orbit evolution is a monotonous (retrograde) rotation of the argument of pericentre, ω. When combined with the potential of the distant perturber |$\overline{\mathcal {R}}_\mathrm{p}$|⁠, the potential of the extended mass generally leads to damping of the K-L oscillations (see Haas & Šubr 2021, for a detailed discussion). This damping stabilizes the zero eccentricity orbit for arbitrary inclination for a suitable choice of system parameters. Note also that in such a situation, monotonous rotation of the longitude of the ascending node remains the primary manifestation of the influence of the distant perturber. Šubr, Schovancová & Kroupa (2009) showed that for damped K-L oscillations, the rate of change in longitude of the ascending node is given by

This equation shows that when the K-L oscillations are damped, |$\frac{\text{d}\Omega }{\text{d} t}$| depends on the semimajor axis through T_K (see equation 2) and will result in differential precession for different orbits.

3.3 Post-Newtonian corrections

It has already been discussed in the literature that relativistic corrections to the gravity of the central body can play a role similar to the spherical potential of the extended mass in secular dynamics (Holman, Touma & Tremaine 1997; Blaes, Lee & Socrates 2002; Karas & Šubr 2007), enforcing a (prograde) rotation of the argument of the pericentre, ω. A straightforward way to implement this relativistic effect within the framework presented above is to use the approximation given by Rubincam (1977). This approximation mimics the rotation of the argument of pericentre due to the relativistic effect of the central body using an additional spherically symmetric potential,

where |$h \equiv \sqrt{GM_\bullet a(1-e^2)}$| is the specific angular momentum of the test particle and c stands for the speed of light. Formally, this potential is equivalent to spherical mass distribution with density profile ρ ∝ r⁻⁵, that is, the form of the averaged potential given by equation (3) with β = −3 may be directly used, giving us the mean potential of the first-order post-Newtonian correction,

Note that, in comparison to the general mean potential for extended mass distribution, equation (8) contains additional dependence on eccentricity through h, and it has one less parameter (M_• versus M_c and a_P). Also note that equation (8) diverges as e approaches unity.

We visualize the damping effect of the post-Newtonian corrections using the isocontours of the perturbing potential in the e–ω space in Fig. 2. In particular, we show three examples of |$\overline{\mathcal {R}}(e,\omega)$| with |$\overline{\mathcal {R}}= \overline{\mathcal {R}}_\mathrm{p} + \overline{\mathcal {R}}_\mathrm{GR}$| for a randomly selected value of C = 0.34 and the properties of the central body and perturber are same as SgrA*(M_• = 4 × 10⁶ M_⊙) and the CWD (M_p = 10⁴ M_⊙ and R_p = 0.1 pc) as described in Section 2. We change the value of the semimajor axis of the inner body, that is, with variable strength of |$\overline{\mathcal {R}}_\mathrm{GR}$| with respect to |$\overline{\mathcal {R}}_\mathrm{p}$|⁠.

In the left panel of Fig. 2, the topology is very similar to that of the middle panel of Fig. 1, which means that |$\overline{\mathcal {R}}_\mathrm{p}$| dominates over |$\overline{\mathcal {R}}_\mathrm{GR}$| in absolute value for most of the parameter space. The middle panel of Fig. 2 shows a setup with a smaller value of semimajor axis, leading to a decrease in the absolute value of |$\overline{\mathcal {R}}_\mathrm{p}$| while, at the same time, it leads to a growth in the absolute value of |$\overline{\mathcal {R}}_\mathrm{GR}$|⁠, which means that it contributes considerably to |$\overline{\mathcal {R}}$|⁠. The topology of the isocontours of |$\overline{\mathcal {R}}$| remains the same as in the previous case, but the overall structure changes so that the separatrix does not reach smaller eccentricity values. Further reduction of the semimajor axis, as shown in the right panel of Fig. 2, leads to |$\overline{\mathcal {R}}_\mathrm{GR}$| fully dominating over |$\overline{\mathcal {R}}_\mathrm{p}$|⁠, and hence the isocontours of |$\overline{\mathcal {R}}$| form nearly circular shapes as a consequence of the independence of |$\overline{\mathcal {R}}_\mathrm{GR}$| on ω. In this case, the K-L oscillations are strongly damped, and the zero eccentricity orbit becomes stable and does not evolve.

For the sake of the analytic treatment of the four-body dynamics described in the following sections, the system configuration must be such that the zero eccentricity orbit is stable. However, due to the non-trivial dependence of |$\overline{\mathcal {R}}_\mathrm{GR}$| and |$\overline{\mathcal {R}}_\mathrm{p}$| on system parameters, this condition must be evaluated from case to case.

In Fig. 3, we evaluate it for parameters of the system that may correspond to the situation in the vicinity of the SgrA* black hole and the semimajor axis of the inner body is sampled within the range 0.01–0.5R_p, which falls into the region of the S-stars for the example setup described in Section 2. We quantify the damping of K-L oscillations by evaluating the maximum value of eccentricity e_max reached by the system during its evolution when starting from near-zero eccentricity. The K-L oscillations are successfully damped when we obtain smaller values of e_max, as the only source of change in eccentricity in these systems is the K-L dynamics. We see that in this setup, the GR effects damp the K-L oscillations for the entire range of C for a ≲ 0.14R_p. At the same time, for a ≳ 0.3R_p, the K-L dynamics is less affected; that is, the zero eccentric orbit is stable only for |$C\gtrsim \sqrt{3/5}$|⁠, shown by the white dashed line in Fig. 3.

3.4 Interparticle potential

In order to describe the four-body setup, Haas et al. (2011b) evaluated the averaged interparticle potential for circular orbits,

where α ≡ a′/a. |$\boldsymbol{n}$| and |$\boldsymbol{n^\prime }$| are the unit vectors that are normal to the mean orbital plane for the two stars, which can be parametrized as |$\boldsymbol{n} = [\sin {i}\sin {\Omega }, -\sin {i}\cos {\Omega }, \cos {i}]^T$| and |$\boldsymbol{n^\prime }=[\sin {i^\prime }\sin {\Omega ^\prime }, -\sin {i^\prime }\cos {\Omega ^\prime }, \cos {i^\prime }]^T$|⁠. We can define the function Ψ as

where P_l(x) are the Legendre polynomials.

We can express the potential energy due to the interaction of inner circular orbits and the outer perturber,

3.5 VHS mechanism

The total averaged potential of the four-body setup described in Section 2 is:

and the classical orbital elements are assumed to evolve according to the Lagrange equations (see e.g. Bertotti, Farinella & Vokrouhlický 2003):

here η and η′ are the mean motion frequencies of the two bodies. Although the average potential due to either the extended mass or relativistic corrections plays an essential role in damping the K-L oscillations of the circular orbits, we may omit it here as it does not contribute to the target subset of Lagrange equations (14) and (15).

The set of equations (14) and (15) with mean perturbing Hamiltonian (equation 13) has been first studied by Haas et al. (2011b) and we refer to their solution in general as the VHS mechanism. These equations may be translated to equations for normal vectors, |$\boldsymbol{n}$| and |$\boldsymbol{n^\prime }$|⁠, of the orbital planes (equations 21–26 of Haas et al. 2011b).

where

and

The frequencies ω_I, |$\omega ^\prime _\mathrm{I}$| and ω_p, |$\omega ^\prime _\mathrm{p}$| correspond to the frequencies caused by the mutual interaction of the two bodies and the perturber, respectively.

Haas et al. (2011b) have shown both by means of analysis of the averaged Hamiltonian as well as direct integration of the Lagrange equations that there exist two qualitatively distinct classes of solutions. On a qualitative level, if the masses of the inner orbits are small enough, or their separation (in terms of semimajor axes) is sufficiently large or a combination of both, we call the regime of interaction weak. In the opposite case, we call the interaction strong. Explicit formula defining boundary between the two modes is not known, nevertheless, an estimate for particular setup can be obtained comparing the frequencies ω_I and |$\omega ^\prime _\mathrm{I}$| to ω_p and |$\omega ^\prime _\mathrm{p}$|⁠. In the strong mode, |$\omega _\mathrm{I},\ \omega ^\prime _\mathrm{I} \gt \gt \omega _\mathrm{p},\ \omega ^\prime _\mathrm{p}$| which means that evolution of the orbital planes described by equation (16) is governed by the mutual interaction of the inner orbits. On the other hand, if |$\omega _\mathrm{I},\, \omega ^\prime _\mathrm{I} \lt \lt \omega _\mathrm{p} , \ \omega ^\prime _\mathrm{p}$|⁠, the weak mode of the VHS mechanism takes place in which the two planes precess differentially due to the gravitational influence of the outer orbit. Note that none of the frequencies ω_I, |$\omega ^\prime _\mathrm{I}$|⁠, ω_p, and |$\omega ^\prime _\mathrm{p}$| are constant over time. Hence, determination of which mode realizes cannot be reliably determined from their initial values.

3.5.1 Weak mode of the interparticle interaction

In the weak mode of the VHS mechanism, the two orbits periodically interchange their angular momenta such that their magnitudes stay constant, but mutual inclination changes. The longitudes of their ascending nodes rotate at different rates while still being mutually influenced. This independent rotation of Ω disrupts the original co-planar configuration. At the moments when ΔΩ ≡ Ω′ − Ω reaches the value of a multiple of 2π, the relative inclination of the two orbits drops to zero, and the planar structure is re-established for the moment.

An example of this solution is shown in Fig. 4 with parameters of the system given in Table 1 under the label M1. Besides showing the orbital evolution according to the secular approximation, we also plot the evolution of the orbital elements coming from direct integration of the equations of motion. For the latter case, we utilize the arwv integrator (Chassonnery, Capuzzo-Dolcetta & Mikkola 2019) since it allows for integrations of a few-body system with up to 2.5 order post-Newtonian approximation. Both solutions share the same qualitative properties with slight differences in the amplitude and period of oscillations of the inclinations which indicate the quality of the secular approximation in this particular configuration.

Estimate of characteristic time-scale, T_char, of the weak mode of the VHS mechanism comes from that (i) the precession of Ω is dominated by the distant perturber, that is, it is nearly constant but different for the two inner bodies and (ii) the period of oscillation of inclinations is determined by the time instances when Ω − Ω′ = 2π. To find T_char for which Ω(T_char) − Ω′(T_char) = 2π, we can apply equation (6) independently to both the inner and outer orbits, approximating inclinations and eccentricities with their initial values (I = I′ = I₀ and e = e′ = e₀) which yields:

For the case of e₀ = 0, this simplifies to the formula given by Haas et al. (2011b, equation 34). Plugging the initial conditions of M1 in equation (20) we get T_char ≈ 192 Myr which is same order of magnitude of |$T_ {\sf {}M1}=123.42$| Myr.

3.5.2 Strong mode of the interparticle interaction

The strong mode of the VHS mechanism occurs when the masses of the inner orbits are large and/or their orbits are closer to each other. In this case, interparticle interaction surpasses the differential precession of Ω and Ω′ induced by the distant perturber, and the orbits corotate. Similarly to the weak case, the two inner orbits keep the magnitude of total angular momenta constant, yet exchange angular momentum so that their inclinations undergo mirrored oscillations. The amplitude of these oscillations is typically smaller than in the weak mode and, therefore, the two orbits stay nearly coplanar during the whole course of the secular orbit evolution.

Fig. 5 shows a typical example of strong mode which occurs in the setup with initial conditions labelled as model M3 in Table 1. Just like the weak case, we show results of the orbital evolution according to the secular approximation as well as direct integration of the equations of motion. The solutions are qualitatively the same which proves that the secular theory is suitable for understanding the nature of the VHS mechanism.

As the differential precession of Ω is suppressed in this mode of the VHS mechanism, it cannot be used to define any characteristic time-scale. Instead, the fact that ω_I and |$\omega ^\prime _\mathrm{I}$| dominate in equation (16) can be used to estimate period of the orbital evolution as |$T_\text{char} = 2 / (\omega _\mathrm{I} + \omega ^\prime _\mathrm{I})$|⁠. For M3, we get T_char ≈ 0.98 Myr, while the numerical integrations give a period of 0.97 Myr.

4 NUMERICAL SOLUTIONS

The secular theory discussed above relies on the assumption of constant zero eccentricity of the two nearby orbits. This section aims to investigate the four-body dynamics that relax this strict constraint. Since no analytic theory is formulated for the non-zero eccentricity case (and is expected to be non-trivial as the dynamics of nearby eccentric orbits is susceptible to chaotic behaviour induced by close encounters), we study our desired setups with sufficient accuracy using direct numerical integrations.

The lack of analytic theory makes it difficult to define distinct classes of possible evolution. Our strategy is then to perform a set of integratons with different initial conditions and compare the results with the ideal cases (zero eccentricity) for which we have an analytic insight. Therefore, the set of examples presented below is likely to be incomplete in terms of all the possible outcomes but it still shows that the two basic modes of the VHS mechanism have identifiable effects in more general setups.

Table 1 lists the initial conditions of the four setups we discuss in this section, along with the zero eccentricity cases discussed in the previous section. A large number of direct integrations with relativistic corrections using arwv were conducted; We selected a subset of the runs to clearly demonstrate the strong and weak modes of the VHS mechanism when we relax certain requirements for secular theory. These individual cases are discussed separately in the following sections.

4.1 Weak mode with e > 0

Let us start by relaxing the condition of zero eccentricity of the two nearby orbits. The model M2 is then straightforwardly derived from M1 simply by changing the initial eccentricity values from zero to 0.721. Since the K-L oscillations are damped, the eccentricity of the orbits does not evolve. We can see this in the temporal evolution of selected orbital elements of the two particles for this setup in Fig. 6. When Ω − Ω′ reaches a multiple of 2π, the relative inclination of the two particles drops to zero. This directly agrees with the angular momentum exchange in the weak mode of VHS mechanism between the two bodies as described in Section 3.5.1.

Period of the secular evolution within the weak mode of the VHS mechanism for model M2 is clearly shorter with respect to the circular case (M1). This is in accord with the dependence of equation (20) for characteristic time-scale on eccentricity. For model M2, it gives T_char ≈ 75 Myr, while the period determined directly from the numerical integrations is |$T_ {\sf {}M2} \approx 47$| Myr.

4.2 Strong mode with e > 0

In another example, we consider a system based on M3, but with an initial eccentricity of 0.77 and we refer to this model as M4. Fig. 7 shows the evolution of the orbital elements of this model. We observe that the inclinations of the two bodies exhibit mirrored oscillations, while the value of ΔΩ oscillates around zero. These two signatures suggest that the system is influenced by the strong mode of VHS mechanism, albeit with some qualitative differences compared to the zero eccentricity case.

Contrary to the previous cases (M1–M3), the orbits undergo non-periodic changes of their semimajor axes, which means that there is a stochastic energy exchange occurring between the two particles. We attribute this to particles on two nearby eccentric orbits occasionally getting so close to each other that the instantaneous two-body scattering noticeably affects their semimajor axes and eccentricities. These scattering events mean that we cannot treat the orbital evolution as secular.

A clear distinction between M3 and M4 is the evolution of the inclination of the two particles. In M3, the orbits evolve in accordance with the secular theory of Haas et al. (2011b), which implies that the inner of the two coplanar orbits is pushed to higher values of inclination, while the inclination of the outer orbit decreases. The evolution is more complex in M4 compared to M3. In M4, the value of Δi ≡ i′ − i periodically changes its sign (see Appendix  A for further discussion). On the other hand, the (quasi)periodic mirrored oscillations of inclinations of the two orbits suggest that the angular momentum transfer between them is secular. The magnitude of the change in inclination is also higher in M4 compared to M3, but still smaller compared to the inclination oscillations present in weak mode of VHS mechanism (models M1 and M2).

Finally, let us focus on the evolution of the longitudes of the ascending nodes Ω and Ω′ of the two particles. If these were test particles, that is, not interacting with each other, Ω and Ω′ would evolve at different constant rates according to equation (6), which means that ΔΩ would grow monotonically in time, reaching a value of ≈28° on the time-scale of |$20\, \mathrm{Myr}$| in the setup of model M4. However, in the bottom panel of Fig. 7, we see limited oscillations of ΔΩ around zero with maximum amplitude ≈10°. Small ΔΩ means that differential precession is suppressed, although not as ideal as in model M3 with zero eccentricity.

Considering the two necessary signatures in the evolution of the orbital elements, that is, small amplitude mirrored oscillations of inclinations and suppressed differential precession in terms of ΔΩ, we state that the system described in model M4 undergoes a generalized mode of the strong mode of VHS mechanism with non-zero eccentricity.

4.3 Strong mode on the top of Kozai–Lidov cycles

Now that we have seen examples of systems with nonzero eccentricity showing either the weak or strong mode of VHS mechanism, we now try to relax the requirement of having constant eccentricity by reducing the damping of K-L dynamics. We can do this by increasing the ratio of a/R_p, which strengthens the perturbing potential due to the outer body with respect to the damping potential due to the post-Newtonian corrections. We study model M5 (Table 1) to explore the VHS mechanism with variable eccentricity.

The left panels of Fig. 8 show the evolution of orbital elements for this system, with the eccentricity oscillations of the two particles now sharing a common period and amplitude. Their inclinations have a more complex evolution, but it is straightforward to identify short-term mirrored oscillations around the mean value. The mean value of the inclination oscillates due to K-L dynamics, which is on a much longer time scale than the strong mode of VHS mechanism. In this case, the inclinations evolve according to the secular theory of Haas et al. (2011b) in that the inclination of the inner body is always greater than that of the outer one. Finally, it is the suppressed differential precession of |$\Omega $| and Ω′ which indicates that we see the two particles moving in the regime where strong mode of VHS mechanism is present, that is, with a mutually locked orientation of their orbital planes while undergoing typical long-term K-L cycles.

Since the particles undergo two independent types of secular evolution at once, we find it beneficial to demonstrate how the orbits will evolve without the VHS mechanism. We can achieve this in the test-particle regime, that is, when mutual interaction between the two inner bodies is suppressed, as shown in the right panels of Fig. 8. Both particles undergo independent K-L oscillations in the test-particle regime with different periods and amplitudes. Difference of the longitudes of the ascending nodes, ΔΩ, systematically (though not monotonically) grows over time. We can also see the period of the K-L oscillations are different between the left and right figures. This means that the VHS mechanism changes T_K of the two bodies so that the new T_K is between the T_K of the two bodies if they were evolving independently.

Let us also point out the apparent regular nature of this setup contrary to the above-discussed model M4 in Section 4.2. This property, however, is not generic as the system is chaotic; slightly modified initial parameters of the system may lead to dramatically different evolution of orbits.

4.4 Transition from the strong to weak mode

It has been demonstrated already in Section 4.2 (model M4) that the systems with non-zero eccentricity may be subject to slightly chaotic evolution due to stochastic close encounters between the two inner particles. Model M6 in Table 1 is another example of a system where such encounters play an essential role. One notable difference from M4 is that the initial eccentricity in the current setup is close to zero but not precisely zero. The left panels of Fig. 9 show the temporal evolution of the setup M6.

From the beginning, until |$T \approx 56\, \mathrm{Myr}$|⁠, it shows an evolutionary pattern similar to that of model M4, that is, inclinations of the two inner particles undergo mirrored oscillations with Δi periodically changing its sign. At the same time, ΔΩ oscillates around zero value, meaning the two orbits corotate and are almost coplanar, that is, the orbits undergo the strong mode of VHS mechanism. Also similar to model M4 is the stochastic (though rather subtle) evolution of semimajor axes and eccentricities.

At |$T\approx 56\, \mathrm{Myr}$|⁠, another close encounter of the two inner particles leads to a more substantial perturbation of their orbits in semimajor axes and eccentricities. Subsequent evolution shows that this event led to the transition from the weak to the strong mode: the inclinations of the two particles exhibit larger amplitude mirrored oscillations. At the same time, longitudes of the ascending nodes precess differentially. At the moments when ΔΩ reaches a natural multiple of 2π, both orbits share the same value of inclination, that is, they are coplanar for that short period.

Another remarkable feature during the phase of weak mode is short-periodic oscillations of eccentricity and inclination of the outer particle. These are K-L oscillations induced by the outer perturbing body that now become less damped because of a suitable angular momentum and energy change. To confirm the nature of these oscillations, we show the evolution of a system of two test particles in the external potential with initial conditions taken from the state of M6 shortly after the two-body scattering event at T = 55.5 Myr in the right panels of Fig. 9. These lighter bodies then have the following orbital parameters: a = 0.0183 pc, a′ = 0.0146 pc, e = 0.21, e′ = 0.11, and Δi = 0.202°. The outer particle, which is more influenced by the distant perturber, undergoes coupled regular oscillations of eccentricity and inclination. In contrast, the oscillations of the inner particle are strongly damped due to the stronger effect of the relativistic precession.

4.5 Disc-like structures

Let us demonstrate the VHS mechanism in the evolution of an N-body system. We study a setup inspired by Haas et al. (2011a) but with two significant differences. First, the initial eccentricities of the orbits are uniformly distributed within the range [0,1), while Haas et al. (2011a) initially considered circular orbits. Second, the post-Newtonian corrections to the central body’s gravity dampen the K-L oscillations instead of the extended mass distribution included in Haas et al. (2011a).

We consider a hypothetical disc of 50 stars orbiting around SgrA*, an SMBH of mass M_• = 4 × 10⁶ M_⊙. The disc is perturbed by a massive perturber of mass M_p = 1 × 10⁴ M_⊙ orbiting SgrA*on a circular orbit at R_p = 0.1 pc. The masses of the stars in the disc are sampled from a Salpeter distribution function, ξ(m) ∝ m^−2.35, in the mass range 1–15 M_⊙.

For all orbits, the initial values of the argument of pericentre ω and the longitude of the ascending node Ω are set to zero. At the same time, other orbital elements are sampled uniformly with |$a \in [0.0035,0.02)\, \mathrm{pc}$|⁠, e ∈ [0.0, 1.0), i ∈ [65°, 75°), and the true anomaly ν ∈ [0, 2π). We integrate this setup with the same integration code, arwv, which we used in the previous sections.

Fig. 10 illustrates the temporal evolution of the orbital elements of all 50 stellar orbits. Fig. 11 shows the projection of the normal vectors of the orbital planes of the same orbits at T = 0, 2.5, and 20 Myr. We currently separate the stars whose Ω stays within 20° of the median of the whole sample throughout the course of evolution and mark them in red in both figures. We refer to this as the disc-like structure as the orbits are corotating with each other. The stars depicted in blue are objects whose orbits rotate independently and visually occupy a more spread out region in Fig. 11.

Although the current configuration differs from the model presented in Haas et al. (2011a), the main dynamical effects are qualitatively similar: approximately 2/3 of the orbits, predominantly from the inner region of the disc, maintain the disc-like configuration, characterized by similar values of both i and Ω, throughout the entire course of evolution. The remaining outer orbits precess differentially in terms of Ω, resulting in a scattered structure. However, this structure still exhibits a specific feature, as the inclinations of these orbits are typically smaller than their initial values.

In contrast, the inclination of the coherent structure grows with respect to the initial value, becoming nearly perpendicular to the orbital plane of the outer perturbing body. We interpret this evolution similarly as was done in Haas et al. (2011b). Specifically, we suggest that the inner orbits mutually interact in the strong VHS mode. Furthermore, the inner and outer parts of the disc initially act as two bodies that mutually interact in the weak VHS mode. After some time, the outer body loses initial coherency due to the differential precession of the orbits of its individual members, which suppresses the weak mode of VHS mechanism between the inner and outer bodies of the disc.

It is worth noting that the model presented here is scaled so that the coherent structure has spatial dimensions similar to those of the system of S-stars observed in the Galactic Centre. While this paper cannot provide any insight into the role of the VHS mechanism in the dynamical evolution of stars in the Galactic Centre, recent research by Ali et al. (2020) suggests that coherent disc-like structures can be identified within the S-star cluster. This presents an opportunity to observe the potential effects of the VHS mechanism on the stars in the Galactic Centre.

5 CONCLUSIONS

In this work, we built on the previous study conducted by Haas et al. (2011b) that explored the dynamical evolution of two nearby, Keplerian, and initially co-planar orbits under the influence of a massive, distant perturber. The secular theory proposed in Haas et al. (2011b) assumes constant zero eccentricity of orbits of all bodies (the two inner objects close to the dominant body and the distant perturber). This assumption is only applicable to systems where an additional non-Keplerian spherically symmetric potential is not only present, but is also strong enough to damp K–L oscillations of the two inner bodies caused by the gravity of the distant perturber. The secular theory provides two qualitatively different solutions of orbital evolution of the inner bodies, which we refer to as the weak and strong modes of the VHS mechanism.

Generally, the weak mode applies when the masses of the bodies on the inner orbits are small, and/or their separation in semi-major axes is large. This mode results in independent rotations of the longitudes of the ascending nodes of the two orbits due to the influence of the distant perturber. Additionally, the two orbits periodically exchange their angular momentum, leading to periodic coupled oscillations of their inclinations. However, when ΔΩ is an integer multiple of 2π, both inner orbits become coplanar again.

For systems with more massive bodies and/or minor separations between the two inner orbits, the strong mode applies. In this mode, the inner orbits have a common rotation rate of Ω, accompanied by oscillations of small-amplitude inclinations.

This paper demonstrates that the qualitative features of the two modes of VHS mechanism are identifiable in systems where some of the critical assumptions of the secular theory are relaxed. Instead of the external potential of some extended mass distribution, post-Newtonian corrections to the dominant body’s gravity can dampen the K-L oscillations. This damping is well understood within the original secular theory of Haas et al. (2011b) with the first-order post-Newtonian approximation given by Rubincam (1977). By relaxing the need for the extended mass to dampen the K-L oscillations, VHS mechanism applies to a broader range of astrophysical systems, such as compact planetary systems or the innermost regions of galactic nuclei.

We have further studied systems with non-zero eccentricity of the inner orbits. We cannot use the secular theory of Haas et al. (2011b) to study such a setup. Nevertheless, by directly integrating the equations of motion, we have identified key features of both the weak and strong modes of the VHS mechanism. The main difference we found in these setups compared to the zero eccentricity case is within the strong mode. In this mode, the orbital inclinations of the inner particles may swap, meaning that in some setups, they oscillate around the common starting value. None the less, this does not change the general statement that the orbits corotate (ΔΩ ≈ 0) within this evolutionary mode.

In order to achieve a more general setup, we have partially relaxed the assumption of constant eccentricity, which assumes complete damping of K-L oscillations of the inner orbits due to the gravity of the outer perturber. We have presented examples of systems where we observe only partially damped K-L oscillations of the inner orbits.¹ The typical features of the VHS mechanism’s weak or strong modes are identifiable in these systems.

Finally, we have demonstrated, similarly to Haas et al. (2011a, b), that the VHS mechanism applies to more complex systems with a larger set of initially coplanar bodies in a relativistic potential. Recent research by Ali et al. (2020) suggests that coherent disc-like structures can be identified within the S-star cluster. This opens avenues for observing the possible effects of VHS mechanism in the stars in the Galactic Centre.

In summary, the analytical expression of VHS mechanism described in Haas et al. (2011b) appears to be a robust phenomenon that can even govern the evolution of systems that do not meet the assumptions of the analytic theory. We have demonstrated through several examples that the VHS mechanism patterns can be found even in systems where instantaneous close encounters significantly affect the orbital evolution. Specifically, the persistent near corotating configuration within the strong mode may have straightforward, observationally detectable consequences for a broad range of astrophysical systems, such as compact planetary systems or stellar structures in the innermost regions of galactic nuclei. However, it is essential to note that the strong mode of the VHS mechanism does not create coplanar and corotating structures within our current understanding; instead, it allows for the survival of existing such structures for extended periods. The weak mode may lead to a specific evolution of its orientation, as shown in Section 4.5, which was discussed for a particular setup in Haas et al. (2011a).

The result of a more general understanding of the VHS mechanism is a potential application in the Galactic Centre to orbits of the S-star cluster. A consequence of evolving eccentric orbits is the introduction of chaos in these systems, which needs to be understood better. Studying this in more detail can facilitate a deeper understanding of the evolution of disc-like structures with the VHS mechanism. These studies will lead to significant insights into the behaviour of astrophysical systems and contribute to a better understanding of the underlying mechanisms that govern their evolution.

FUNDING

MS is supported by the Grant Agency of Charles University under the grant no. 179123. LŠ and JH acknowledge support from the Grant Agency of the Czech Republic under the grant 20-21855S.

Acknowledgement

We thank Sai Sasank Chava and Yugantar Prakash for feedback on the manuscript. We thank David Vokrouhlický for his input on using the Rubincam approximation.

DATA AVAILABILITY STATEMENT

The data and tools used to produce the plots in this paper will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: INCLINATION CROSSING IN ECCENTRIC STRONG MODE

In Section 4.2, we describe qualitative difference of the strong mode of the VHS mechanism with eccentric orbits in comparison to the circular case. It has been argued in Haas et al. (2011b) that, starting from coplanar configuration, the inclination of the inner orbit always grows, while that of the outer one decreases. An important piece of their argument is that precession of the outer orbit due to the distant perturber is always faster which leads to positive value of sin (Ω′ − Ω) which implicitly occurs in equations (14) and (15) through the dependence of |$\overline{\mathcal {R}}_\mathrm{i}$| on |$\boldsymbol{n}.\boldsymbol{n}^\prime$|⁠.

We do not have secular equations for the VHS mechanism with eccentric orbits in hands, still, we may assume that dependence of di/dt and di′/dt on ΔΩ ≡ Ω′ − Ω is similar to the circular case. Fig. A1 shows zoomed-in evolution of model M4 for a short period of time. Indeed, we see that, in contrary to the circular case, ΔΩ reaches non-zero (both positive and negative) vales at the instances of i = i′. Depending on the sign of ΔΩ, inclination of the inner orbits either grows similarly to the circular case (ΔΩ > 0) or decreases. For comparison, we also show detailed view of evolution of orbital elements for setup similar to model M4, but now with small eccentricities of the two inner orbits, |$e_0 = e_0^\prime = 0.08$|⁠, in Fig. A2. The oscillatory pattern of ΔΩ is preserved, but now with (i) several orders of magnitude smaller amplitude, (ii) near zero value at the instances of i ≈ i′, and (iii) positive derivative at those instances. Evolution of i and i′ is then in accord with the analytic argumentation for the zero eccentricity case.

Let us note that due to lack of analytic secular theory for the non-zero eccentricity case, it is hard to discriminate, whether evolution of ΔΩ in the strong mode of the VHS mechanism is primarily due to non-uniform precession of the orbits in the field of the distant perturber, or whether it is mainly governed by their mutual torques.