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

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 inter-orbital 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.


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, e.g., 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 et al. (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 semi-major axes and ★ Contact e-mail: singhal.myank@matfyz.cuni.czare initially co-planar (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 synchronised, i.e., the initial co-planar 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 co-planar configuration.We refer to the temporal evolution of the specific four-body setup introduced by Haas et al. (2011b) as the VHS mechanism 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 Kozai-Lidov oscillations, which reduces its applicability in observed astrophysical systems.However, Haas et al. (2011a,b) introduced a physically realistic setup in which the VHS mechanism is applicable.They studied a system in which the super-massive 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 circum-nuclear 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 co-planar orbits to be the observed stars from the young stellar disc that is within the distances of 0.04 -0.4 pc from the central super-massive black hole (Levin & Beloborodov 2003;Paumard et al. 2006; Bartko et al. 2009Bartko et al. , 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.Firstly, 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, i.e., 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 generalised VHS dynamics in Section 5.

SETUP
We study a hierarchical four-body system with a dominant central body, characterized only by its mass,  • .The system further consists of a distant perturber of mass  p on a circular orbit with radius  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  and  ′ where ,  ′ ≪  p on orbits around the central body with a semi-major axes  and  ′ , which are much smaller than  p and  ′ < .These light bodies are in inclined orbits, having inclinations  and  ′ 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  =  ′ 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  • ,  p and  p .We set the system that may correspond to the situation in the vicinity of the SgrA * black hole, i.e.,  • = 4 × 10 6 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  p = 10 4 M ⊙ and a semi-major axis of  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.

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 Sec. 2. For this, the mean interaction potential of the system, 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.

Potential of the distant / outer perturber
The averaged interaction potential between a perturbing body on a circular orbit with radius  p and a particle on an orbit with semimajor axis , eccentricity  and inclination  with respect to the orbital plane of the perturber reads (Kozai 1962): where  is the argument of periapses of the orbit.Suppose R 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,  ≡ √ 1 −  2 cos , which, together with the semi-major axis () and R 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 R in the - space, which for fixed values of  and  give sets of possible evolutionary tracks (see Figure 1).These sets form two qualitatively different topologies: For  > √︁ 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 Figure 1).If  ≤ √︁ 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 panel of Figure 1).The lower the value of , 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  > √︁ 3/5, while it undergoes periodic variations when  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 semi-major 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.The potential isolines when the potential due to post newtonian approximation is added to the perturbing potential are displayed in three panels, representing two cases with identical parameters except for the semi-major axis of the test particle.In these examples the value of  = 0.34 which is smaller than the limiting value for K-L dynamics, and the separatrix is shown in red.The panel on the left depicts the instability at  = 0 when  = 0.2  , while the panel on the right illustrates the stability at  = 0 when  = 0.03  .The central panel shows the boundary area when  = 0.145  and the K-L oscillations are damped.

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  • .Although these sources differ, they have very similar effects and impact the evolution of the two bodies in a similar manner.Haas et al. (2011b) and Haas et al. (2011a) 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,

Extended Mass
where   stands for the integral of the extended mass density within the orbit of the perturber ( p ) and The coefficients   are given by with  1 = (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 R p , the potential of the extended mass generally leads to damping of the Kozai-Lidov 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 et al. (2009) showed that for damped K-L oscillations, the rate of change in longitude of the ascending node is given by dΩ d This equation shows that when the K-L oscillations are damped, dΩ d depends on the semi-major axis through  K (see Equation 2) and will result in differential precession for different orbits.

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 et al. 1997;Blaes et al. 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 pericenter due to the relativistic effect of the central body using an additional spherically symmetric potential, where ℎ ≡ √︁   • (1 −  2 ) is the specific angular momentum of the test particle and  stands for the speed of light.Formally, this potential is equivalent to spherical mass distribution with density profile  ∝  −5 , 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 ℎ, and it has one less parameter ( • vs.   and   ).Also note that Equation 8 diverges as  approaches unity.
We visualise the damping effect of the post-Newtonian corrections using the isocontours of the perturbing potential in the - space in Figure 2. In particular, we show three examples of R (, ) with R = R p + R GR for a randomly selected value of  = 0.34 and the properties of the central body and perturber are same as SgrA * ( • = 4 × 10 6 M ⊙ ) and the CWD ( p = 10 4 M ⊙ ,  p = 0.1pc) as described in Section 2. We change the value of the semi-major axis of the inner body, i.e., with variable strength of R GR with respect to R p .
In the left panel of Figure 2, the topology is very similar to that of the middle panel of Figure 1, which means that R p dominates over R GR in absolute value for most of the parameter space.The middle panel of Figure 2 shows a setup with a smaller value of semi-major axis, leading to a decrease in the absolute value of R p while, at the same time, it leads to a growth in the absolute value of R GR , which means that it contributes considerably to R. The topology of the isocontours of 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 semi-major axis, as shown in the right panel of Figure 2, leads to R GR fully dominating over R p , and hence the isocontours of R form nearly circular shapes as a consequence of the independence of R 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 R GR and R p on system parameters, this condition must be evaluated from case to case.
In Figure 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 semi-major axis of the inner body is sampled within the range 0.01 − 0.5 p , which falls into the region of the S-stars for the example setup described in Sec 2. We quantify the damping of K-L oscillations by evaluating the maximum value of eccentricity  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  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  for  ≲ 0.14 p .At the same time, for  ≳ 0.3 p the K-L dynamics is less affected; that is, the zero eccentric orbit is stable only for  ≳ √︁ 3/5, shown by the white dashed line in Figure 3.

Inter-particle potential
In order to describe the four body setup, Haas et al. (2011b) evaluated the averaged inter-particle potential for circular orbits, where  ≡  ′ /. and  ′ are the unit vectors that are normal to the mean orbital plane for the two stars, which can be parameterized as  = [sin  sin Ω, − sin  cos Ω, cos ]  and  ′ = [sin  ′ sin Ω ′ , − sin  ′ cos Ω ′ , cos  ′ ]  .We can define the function Ψ as where   () are the Legendre polynomials.
We can express the potential energy due to the interaction of inner circular orbits and the outer perturber,

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 et al. 2003): here  and  ′ are the mean motion frequncies 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, Equation 14) & 15.
The set of Equation 14& 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,  and  ′ , of the orbital planes (Haas et al. 2011b, Equations 21-26). where and The frequencies  I ,  ′ I and  p ,  ′ 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 semi-major 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  ′ I to  p and  ′ p .In the strong mode,  I ,  ′ I >>  p ,  ′ p which means that evolution of the orbital planes described by Equation 16is governed by the mutual interaction of the inner orbits.On the other hand, if  I ,  ′ I <<  p ,  ′ 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 ,  ′ I ,  p and  ′ p are constant over time.Hence, determination of which mode realises cannot be reliably determined from their initial values.

Weak mode of the inter-particle 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 coplanar 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 Figure 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 et al. 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,  char , of the weak mode of the VHS mechanism comes from that (i) the precession of Ω is dominated by the distant perturber, i.e., 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  char for which Ω( char ) − Ω ′ ( char ) = 2, we can apply Equation 6 independently to both the inner and outer orbits, approximating inclinations and eccentricities with their initial values ( =  ′ =  0 and  =  ′ =  0 ) which yields: For the case of  0 = 0, this simplifies to the formula given by Haas et al. (2011b, Equation 34).Plugging the inital conditions of M1 in Equation 20we get  char ≈ 192 Myr which is same order of magnitude of  M1 = 123.42Myr.The lighter version of the lines is the result of integration of secular equations while the darker versions are the result of integration of equations of motion.The black dashed lines highlight that the mutual inclination becomes 0 when the value of ΔΩ is a multiple of 2  in the integration of equations of motion.Similarly the grey line is for the integration of secular equations.

Strong mode of the inter-particle 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, inter-particle interaction surpasses the differential precession of Ω and Ω ′ induced by the distant perturber, and the orbits co-rotate.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.
Figure 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 timescale.Instead, the fact that  I and  ′ I dominate in Equation 16 can be used to estimate period of the orbital evolution as  char = 2/( I + ′ I ).For M3 we get  char ≈ 0.98 Myr, while the numerical integrations give a period of 0.97 Myr.

NUMERICAL SOLUTIONS
The secular theory discussed above relies on the assumption of constant zero eccentricity of the two nearby orbits.This section aims to The stronger mutual interaction between the two stars results in strong mode of VHS mechanism, resulting in a constant ΔΩ = 0 • , which is a complete overlap for integration of both secular and equations of motion.
There are tiny inclination oscillations for both integration of equations of motion (darker) and secular equations (lighter), which have slightly different amplitude and period.
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.

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 Figure 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 Sec.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 20for characteristic time-scale on eccentricity.For model M2 it gives  char ≈ 75 Myr, while the period determined directly from the numerical integrations is  M2 ≈ 47 Myr.

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. Figure 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 semi-major 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 semi-major 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 Δ ≡  ′ −  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, i.e., 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 Myr in the setup of model M4.However, in the bottom panel of Figure 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, i.e., small amplitude mirrored oscillations of inclinations and suppressed differential precession in terms of ΔΩ, we state that the system described in model M4 undergoes a generalised mode of the strong mode of VHS mechanism with non-zero eccentricity.

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 / 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 Figure 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 Ω and Ω ′ which indicates that we see the two particles moving in the regime where strong mode of VHS mechanism is present, i.e. 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 Figure 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  K of the two bodies so that the new  K is between the  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.

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 Figure 9 show the temporal evolution of the setup M6.
From the beginning, until  ≈ 56 Myr, it shows an evolutionary pattern similar to that of model M4, i.e., inclinations of the two inner particles undergo mirrored oscillations with Δ periodically changing its sign.At the same time, ΔΩ oscillates around zero value, meaning the two orbits co-rotate and are almost co-planar, i.e., the orbits undergo the strong mode of VHS mechanism.Also similar to model M4 is the stochastic (though rather subtle) evolution of semi-major axes and eccentricities.At  ≈ 56 Myr, another close encounter of the two inner particles leads to a more substantial perturbation of their orbits in semi-major 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, i.e., they are co-planar 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  = 55.5 Myr in the right panels of Figure 9.These lighter bodies then have the following orbital parameters:  = 0.0183,  ′ = 0.0146,  = 0.21,  ′ = 0.11, and Δ = 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.

Disc like structures
Let us demonstrate the VHS mechanism in the evolution of a 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 disk of 50 stars orbiting around SgrA * , a supermassive black hole of mass  • = 4 × 10 6 M ⊙ .The disk is perturbed by a massive perturber of mass  p = 1 × 10 4 M ⊙ orbiting SgrA * on a circular orbit at  p = 0.1 pc.The masses of the stars in the disk are sampled from a Salpeter distribution function,  () ∝  −2.35 , in the mass range 1 − 15M ⊙ .
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  ∈ [0.0035, 0.02) pc,  ∈ [0.0, 1.0),  ∈ [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.
Figure 10 illustrates the temporal evolution of the orbital elements of all 50 stellar orbits.Figure 11 shows the projection of the normal vectors of the orbital planes of the same orbits at  = 0, 2.5 and 20 Myrs.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 disclike structure as the orbits are co-rotating with each other.The stars  Evolutionary tracks depicted with red solid lines correspond to orbits that are part of the coherent structure the whole integration time.Blue dotted lines correspond to orbits that get more separated from the coherent structure for at least some period of time.
depicted in blue are objects whose orbits rotate independently and visually occupy a more spread out region in Figure 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  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.

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 co-planar 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.Nonetheless, this does not change the general statement that the orbits co-rotate (ΔΩ ≈ 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. 1he 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) and Haas et al. (2011b), that the VHS mechanism applies to more complex systems with a larger set of initially co-planar 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 co-rotating 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 co-planar and co-rotating 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 disk-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.

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 co-planar 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 Equation 14 and 15 through the dependence of R i on . ′ .
We don't have secular equations for the VHS mechanism with eccentric orbits in hands, still, we may assume that dependence of d/d and d ′ /d on ΔΩ ≡ Ω ′ − Ω is similar to the circular case.Figure 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  =  ′ .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,  0 =  ′ 0 = 0.08, in Figure A2.The oscillatory pattern of ΔΩ is preserved, but now with (i) several orders of magnitude smaller amplitude and (ii) near zero value at the instances of  ≈  ′ and (iii) positive derivative at those instances.Evolution of  and  ′ 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 nonzero 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.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Isocontours of R for different fixed values of the Kozai integral (, indicated above individual panels) in a three body system showing pure Kozai-Lidov dynamics.The shape of these isocontours is independent on .The leftmost panel shows isocontours in  −  space for  < √︁ 3/5, which shows large changes in eccentricity over the orbit.The middle panel is for  still smaller than √︁ 3/5 but with separatrix (displayed with red line) reaching smaller values of eccentricity.In the rightmost panel we have  > √︁ 3/5 and we see oval evolutionary tracks, meaning small changes in eccentricity.

Figure 3 .
Figure 3. Heatmap of largest variation in eccentricity for a 1 ⊙ star initially on an orbit of eccentricity  = 10 −4 in the combined potential of a central dominant body ( • = 4 × 10 6 M ⊙ ) and a distant perturber ( p = 10 4 M ⊙ ,  p = 0.1 pc) on a circular orbit.The lower value of maximum eccentricity in a significant range of the parameter space shows that the relativistic corrections due to the SMBH partially or entirely dampen the Kozai-Lidov oscillations.

Figure 4 .
Figure 4. Evolution of model M1 showing the weak mode of evolution in VHS mechanism.Weak mode of VHS mechanism results in different rate of change in evolution of ΔΩ along with oscillations in inclination.The lighter version of the lines is the result of integration of secular equations while the darker versions are the result of integration of equations of motion.The black dashed lines highlight that the mutual inclination becomes 0 when the value of ΔΩ is a multiple of 2  in the integration of equations of motion.Similarly the grey line is for the integration of secular equations.

Figure 5 .
Figure5.These graphs show the evolution of the orbital parameter of model M3.The stronger mutual interaction between the two stars results in strong mode of VHS mechanism, resulting in a constant ΔΩ = 0 • , which is a complete overlap for integration of both secular and equations of motion.There are tiny inclination oscillations for both integration of equations of motion (darker) and secular equations (lighter), which have slightly different amplitude and period.

Figure 6 .
Figure 6.Example of weak mode of VHS mechanism with non-zero eccentricity.It shows the evolution of the model M2 which is similar to M1 but with non-zero eccentricity.The evolution is similar to Figure 4 and we again see Δ = 0 when ΔΩ = 2  as shown by the black dashed lines.

Figure 7 .
Figure 7.This figure exemplifies how the strong mode of VHS mechanism behaves with non-zero eccentricity.It show the evolution of model M4 which is similar to M3 but with orbits with eccentricity of  = 0.77.

Figure 8 .Figure 9 .
Figure8.The left panel shows an example run of strong mode of VHS mechanism with dynamically evolving eccentricity due to K-L oscillations in model M5.The interaction between the two stars results in a combination of the strong mode of the VHS mechanism and K-L oscillations in both bodies.The strong modes constant ΔΩ is still present and the characteristic oscillations in the inclination overlap with the K-L oscillations.However, the right panel shows the evolution of orbital elements when we remove the effects of VHS mechanism by decreasing masses of the two inner particles.The constant zero ΔΩ changes to a systematic growth while the two particles have independent K-L oscillations in inclination and eccentricity.

Figure 10 .
Figure10.Evolution of orbital elements of individual stars within the model described in Section 4.5.Evolutionary tracks depicted with red solid lines correspond to orbits that are part of the coherent structure the whole integration time.Blue dotted lines correspond to orbits that get more separated from the coherent structure for at least some period of time.

∆ΩFigure A1 .Figure A2 .
Figure A1.Evolution of the orbital elements within model M4 for a short period of time.The two points where Δ = 0 are marked with green (dotted) and red (dash-dot) vertical lines, with appropriate markings in ΔΩ.