Long-term perturbations in four-body systems with mutual highly inclined orbits

Astrophysics and Space Science (2023) 368:42 https://doi.org/10.1007/s10509-023-04200-7 RESEARCH Long-term perturbations in four-body systems with mutual highly inclined orbits C.V. Monzón1,2 ·J.A. Docobo1,2,3,4 Received: 14 November 2022 / Accepted: 21 May 2023 / Published online: 29 May 2023 © The Author(s) 2023, corrected publication 2023 Abstract Hierarchical Three Body systems with mutual highly inclined orbits have been well-known since the 1960s. Lidov-Kozai cycles arise within them where the inner orbit eccentricity acquires extreme values. In particular, we focus our research on the motion of exoplanets and exomoons on different Three Body stellar scenarios. Our goal is to study how the LK cycles are perturbed by a fourth body (which we called perturbed LK). We analyze the evolution of the eccentricity and inclination of the inner orbit in two cases: the first involves an exoplanet and the second involves an exomoon. Due to the possible stable configurations of a four-body system, we treat two subcases as well: the Totally Hierarchical Configuration and the 2+2 configuration. According to that derived from the particular scenarios of study discussed in the present research, the LK perturbed in exomoon orbits around exoplanets seem to exhibit, in general, way less alterations than the exoplanet orbit around its star. Keywords Celestial mechanics 1 Introduction The Three Body Problem (TBP) begins with the Lunar problem (Gutzwiller 1998). Their objective was to predict the position of the Moon with the highest possible degree of precision. Undoubtedly, the Lunar problem is the most studied case of the TBP and future, present and past of Celestial Mechanics, in general. Lately, as a result of the Lunar problem, the Stellar TBP emerged (Docobo et al. 2021) which became very relevant, especially in the first half of the twentyfirst century. In this scenario, a distant third star perturbs the orbit of a binary star. With these conditions and, unlike in the Moon Problem, the masses of the three bodies are generally comparable and the mutual inclination between the inner and outer orbit is not necessarily small. We know that the Hamiltonian related to the hierarchical TBP can be expanded as a power series depending on a small parameter, which is equal either to the quotient between the semi-axis of the inner and outer orbits, or the quotient between the inner orbit semi-axis and the pericenter distance of the outer orbit. With that in mind, we can study the problem with different perturbation orders and, at the same time, take into account periodic (short and long period) and secular terms. In that regard, some averaging methods for periodic variables were developed (Lagrange 1782; Laplace 1784). Some of the most used averaging methods are: Von Zeipel (1916), Krylov and Bogoliubov (1947), and Deprit (1969). In particular, Harrington (1968,1969,1970) applied the von Zeipel method to average the short period variables. This author also verified that the perturbations are much larger when the mutual inclination is between 39◦and 141◦, as Lidov (1961) and Kozai (1962) had stated some years before, hence the name of these perturbations: the Lidov-Kozai cycles (LK hereafter). These cycles have been observed, for C.V. Monzón carlos.v[email protected] J.A. Docobo [email protected] 1Observatorio Astronómico R. M. Aller, Universidade de Santiago de Compostela, Avda. das Ciencias s/n, Campus Vida, Santiago de Compostela, Galicia, Spain 2Departamento de Matemática Aplicada, Facultade de Matemáticas, Universidade de Santiago de Compostela (USC), Santiago de Compostela, Galicia, Spain 3CITMAga, Campus Vida, E-15782 Santiago de Compostela, Galicia, Spain 4Facultad de Ciencias, Real Academia de Ciencias de Zaragoza, C/ Pedro Cerbuna 12, E-50009 Zaragoza, Spain 42 Page 2 of 16 C.V. Monzón, J.A. Docobo example, in the irregular satellites of Uranus. Later on, Docobo (1977) applied the Deprit method (Deprit 1969)to transform canonical systems into systems that are easier to integrate. In the 1990s, the detection of the first exoplanets orbiting pulsars (Wolszczan and Frail 1992) and main sequence stars (Mayor and Queloz 1995) promoted, once again, the study of the hierarchical TBP, especially those in which the third body is an exoplanet. In these configurations, there are three types of stable1orbits (Dvorak 1984): the S-type (the exoplanet orbits one component of the binary star), the P-type (the exoplanet orbits the two components of the binary star), and the L-type (the exoplanet orbits the L4 or L5 Lagrange point). The stability of an exoplanet in binary star systems has been extensively studied (Harrington 1977; Szebehely 1980; Rabl and Dvorak 1988;Dvoraketal.1989; Pendleton and Black 1983) and several stability criteria have been proposed, for example, Domingos et al. (2006). Ford et al. (2000), Takeda and Rasio (2005), Naoz et al. (2013), and Docobo et al. (2021) studied the evolution of the elements of the inner orbit of a hierarchical Three Body system, taking into account the long term effects due to the Lidov-Kozai cycles. In Docobo et al. (2021), the authors utilized the TIDES numerical integrator to study the solution of the equations in the hierarchical TBP. Takeda and Rasio (2005) and Naoz et al. (2013) averaged the Hamiltonian related to that problem, truncating it at the second term (the so-called quadrupole approximation) and at the third term (octupole approximation), respectively. The main goal of the present work is to analyze how the Lidov-Kozai cycles of the TBP are perturbed when a fourth body is added to the system. We call this as perturbed LK. Although the dynamics of 3 +1 and 2 +2 configurations do appear in the literature (Hamers et al. 2015; Hamers and Lai 2017), their authors tackle the problem in a general manner, and no distinctions are made between exoplanets and theoretical exomoons, for example. In this report, we compare the differences in the long-term evolution of the eccentricity and inclination of the inner orbit between the three and fourbody systems, which is something that seems to be absent in the literature. 2Methods In order to run all of the necessary simulations, we needed to integrate the full n-body problem equations. For that, we used the Mathematica package TIDES, which is specialized software that can solve differential equation systems whose computation time is very long. It was developed in 2011 by 1By stable, we refer to the Lagrange sense, meaning that the orbit is bounded. Alberto Abad, Roberto Barrio, Fernando Blesa, and Marcos Rodríguez in the Grupo de Mecánica Espacial of the Universidad de Zaragoza (Abad et al. 2011a,2012,2015). The code is available at: https://sourceforge.net/projects/tidesodes/. 2.1 TIDES algorithm This package is based on the Taylor Series Method (TSM) to solve the Initial Problem Value (IPV): dy dt =f(t,y(t);p),y(t0)=y0,t∈R,y∈Rn,p∈Rm, (1) y0being the initial conditions, pis the parameters and t,the independent variable. Assuming that fis infinitely differentiable, we can find the approximate solution in a lapse of time ti+1=ti+hi, expanding the Taylor series in y(t) around the point tiand evaluating it in ti+1: y(t0)=y0 y(ti+1)≈y(ti)+dy(ti) dt hi+1+1 2! d2y(ti) dt2h2 i+1+··· +1 n! dny(ti) dtnhn i+1 (1) ≈y(ti)+f(t i,y(ti))hi+1 +1 2! df (ti,y(ti)) dt h2 i+1+···+ 1 n! dn−1f(t i,y(ti)) dtn−1hn i+1. (2) TIDES calculates the n−1 successive derivatives of f to solve the IVP via the so-called automatic differentiation or AD, which permits a very efficient calculation. Since TIDES has an adaptive step, the method is sufficiently robust with no need for a huge amount of integration points. This is a big advantage because, in many cases, we will have to integrate over millions or even tens of millions of years. Hence, the variable step permits a smaller amount of integration nodes without losing precision. This allows TIDES to be very stable because it takes into account the sensitivity of the differential equation with respect to initial conditions and/or parameters. One advantage of TIDES relative to other typical integration algorithms (such as the Runge-Kutta method) is its efficiency, even when we account for high-precision integrations. It has been applied in several astronomical problems and, when comparing TIDES to other packages such as odex, based on the interpolation method, and dop853, based on the RK method (Hairer et al. 1993), it performs considerably better. Technically, the RK method is found to be faster than TIDES when dealing with relatively small tolerances. It is important to note that the precision error is Perturbations in four-body systems Page 3 of 16 42 around 10−16 (16 digits) in our simulations. For that precision level both odex and dop853 required much more CPU time (Abad et al. 2011a,b,2012,2015; Barrio et al. 2011;Simó2001). In summary, TIDES is capable of dealing with very small precision errors without sacrificing CPU time efficiency. 2.2 TIDES structure TIDES consists basically of MathTIDES, a Mathematica package that writes the differential equation system, the initial conditions, and the parameters in C code (or Fortran) and the library LibTIDES, which we link and compile to the C (Fortran) code generated by MathTIDES. Moreover, TIDES includes the GMP and MPFR libraries. These libraries allow double and multiple precision, respectively. 2.3 Other features of TIDES Other interesting features are: – Search of the zeros of a function of the IVP solution that we choose. – Search of the local minima and local maxima. – Integration of the partial derivatives over the variables or the parameters of the IVP that we choose. 3 Long-term perturbations in Three Body systems Long-term perturbations in hierarchical, Three Body systems with mutual highly inclined orbits have been well known for several decades (Lidov 1961 and Kozai 1962). The latter observed that some Asteroid Belt objects with orbits inclined with respect to the ecliptic had a remarkable eccentric orbit as well, due to the action of Jupiter. For example, the orbital eccentricity and inclination of (1373) Cincinnati (a=3.42 au) oscillates between 0.25 −0.6 and 25◦−42◦, respectively. This effect has also been noted in the irregular satellites of Uranus, since none of these objects (except one) have an orbital inclination larger than 50◦or smaller than 140◦. This coincides significantly with the oscillation range where the LK cycles are maximum: 39.23◦−140.76◦(Kozai 1962). In general, when the inclination of the inner orbit with respect to the outer orbit falls between these range of inclinations, LK cycles can excite initially low eccentricities into almost 1 in the long-term (Lidov 1961; Naoz et al. 2013; Takeda and Rasio 2005; Docobo et al. 2021). In this Section, we will examine the differences between the exoplanet and exomoon orbits undergoing LK cycles. Since we are interested in studying the exoplanet and exomoon dynamics, we will consider two cases: one where the inner orbit is that of an exoplanet orbiting a star (hereafter, the exoplanet case) and the other, where the inner orbit is that of an exomoon orbiting an exoplanet (hereafter, the exomoon case). The mass choices for the exoplanet and exomoon, which we will maintain through this and the following Section, are 1 MJand 1 M⊕, respectively. In both cases, we will assume the outer orbit is that of a star orbiting the inner pair. We will consider all of the masses to be point-like and we will not take into account tidal or relativistic effects throughout this work. We have made a series of assumptions, both in the exoplanet and exomoon cases. For the former, a) eout = 0.6; since, in the observed wide binaries distribution, we noted some considerable eccentric orbits, the average being around 0.59 (Tokovinin and Kiyaeva 2015), b) aout < 100 au, just to set an upper limit to the outer semi-axis, c) the eccentricities of the inner orbit are low, d) the initial inclination between the outer and inner orbits is between 60◦and 70◦,2and e) the mass of the third body is that of a brown dwarf. With these assumptions in mind, we have carried out several simulations using TIDES for certain initial values of aout ∈{20,30,50,75,90}au, ein ∈{0.001,0.1,0.2}and the mutual inclination between the inner and outer orbits ∈{60◦,65◦,70◦}. The initial values of T,and ωare the same as in Table 1. We noticed two things. On one hand, orbital flips seem to occur in all cases. On the other hand, there are periods of time where the orbital eccentricity of the exoplanet approaches values extremely close to 1. We carried out an analogous study in the exomoon case but, this time, we also used three stellar mass values (brown dwarf, red dwarf, 1 M). We assumed that a) eccentricities of the inner orbit are low, b) eccentricities of the outer orbit can be moderate (0 ≤eout ≤0.6), c) the initial inclination between the outer and inner orbits is between 60◦and 70◦, and d) the stellar mass is not larger than of the Sun. The different values we used are: aout ∈{3,5,10}au, eout ∈{0,0.1,0.2,0.4,0.6},ein ∈ {0.05,0.1,0.2,0.3}, and the inclination between the outer and inner orbits ∈{60◦,65◦,70◦}. The initial values of T, and ωare the same as in Table 2. Unlike the exoplanet case, we have found no orbital flips in the exomoon orbit. Consequently, the eccentricity, although it is highly excited to values close to 0.8 in most cases, hardly ever reached close to 1 values, except when we placed the exomoon very close to the stability limit.3 2This choice is made to ensure the presence of the LK cycles. 3Throughout this work, the stability limit will always be obtained empirically, by finding the biggest semi-axis in which the exoplanet/exomoon orbit is Lagrange-stable for 30 (exoplanet case) or 0.5 Myr (exomoon case). 42 Page 4 of 16 C.V. Monzón, J.A. Docobo Table 1 Initial orbital elements of the baseline example for the exoplanet case. In each case, the orbital elements are referenced to the invariable plane of the system T(years) ea(au) i(◦)(◦)ω(◦) Inner orbit 0 0.001 6 64.7 240 45 Outer orbit 0 0.6 100 0.3 60 0 Table 2 Initial orbital elements of the baseline example for the exomoon case. In each case, the orbital elements are referenced to the invariable plane of the system T(years) ea(au) i(◦)(◦)ω(◦) Inner orbit 0 0.05 0.1 60 240 45 Outer orbit 0 0 10 0.5 60 0 Fig. 1 Exoplanet initial orbital inclination: 64.7◦. Outer orbit initial inclination: 0.3◦. Inner initial orbital semi-major axis: 6 au. Outer initial orbital semi-major axis: 100 au. Exoplanet initial orbital eccentricity: 0.001. Outer orbit initial eccentricity: 0.6. Exoplanet initial orbit argument of the pericenter: 45◦. Outer orbit argument of the pericenter: 0◦. Exoplanet initial orbit ascending node longitude: 240◦. Outer orbit ascending node longitude: 60◦. All periastron epochs are equal to 0. Example taken from Naoz et al. (2013) to check if the results match and our integration is correct As baseline examples upon which we will build the simulations in the next Section, we chose that of Tables 1and 2. Figure 1shows the eccentricity and inclination4evolution that undergoes the orbit of a Jupiter-type exoplanet that revolves around a Solar-type star and, at the same time, of a brown dwarf (with mass of 40 MJ) that orbits the barycenter of both. The eccentricity, initially 0, is excited close to 1 in 5 Myr. The inclination, on the other hand, oscillates between 40◦and 140◦, which produces orientation changes in the exoplanet orbit, periodically changing from prograde to retrograde motion, and vice versa. In addition, the argument of the pericenter circularizes. Figure 2shows the eccentricity and inclination evolution of the orbit of an Earth-type exomoon revolving around a Jupiter-type exoplanet whose barycenter orbits a brown dwarf. In this example, although the eccentricity of the exomoon orbit does increase significantly from 0 to 0.8, we can see that this increase is less extreme than in the exoplanet case. The same goes for the inclination, where we can see that it never changes its orientation, oscillating between 35◦ and 65◦. Moreover, the argument of the pericenter, as in the previous case, circularizes in the given examples (Fig. 3). However, it can be proven that is not something characteristic of the LK cycles because, in our exoplanet case, if we assume that both stars are Solar-type with the same initial orbit inclination, the argument of the pericenter librates around 90◦, which indicates that the behavior of the argument of the pericenter will depend on the orbital and physical parameters of the bodies in the system.5 In summary, in all of the examples we have considered throughout this Section, the LK cycles undergoing the exomoon orbit seem to be less extreme than those of the exoplanet orbit. It is worth remarking that this is not a general conclusion but a result obtained from a particular set of initial conditions associated with each case of study. 4In this report, the orbital inclination iis referenced to the invariable plane of the system. Moreover, the ascending node longitude  is measured from the positive x-axis of the barycentric reference system adopted for the present study, while the argument of pericenter ω is defined on the orbital plane from the nodal line. 5It is well-known that, in a restricted TBP, the argument of the pericenter librates for such range of orbital inclinations, but our scenarios of work represent a general TBP, since none of the masses are negligible. Perturbations in four-body systems Page 5 of 16 42 Fig. 2 Initial exomoon orbital inclination: 60◦. Initial exomoon orbital eccentricity: 0.05. Initial exomoon orbital semi-axis: 0.1au.Initial outer orbit semi-axis: 10 au. Initial outer orbit eccentricity: 0. Initial outer orbit inclination: 0.5◦. Initial exomoon orbit argument of the pericenter: 45◦. Initial outer orbit argument of the pericenter: 0◦.Initial exomoon orbit ascending node longitude: 240◦. Initial outer orbit ascending node longitude: 60◦. All periastron epochs are equal to 0 Now, why do we see such differences in the LK cycles between the exoplanet and the exomoon cases? We consider that it has to do with the difference in the stability regions. Due to the masses of the bodies in the systems, in order for the bodies to be in a stable orbit, the quotient between the inner orbit and outer orbit semi-axis (the small parameter we mentioned in the Introduction, ) has to be much smaller in the exomoon case. Therefore, if we expand the Hamiltonian of the hierarchical TBP as a function of this parameter, we have the next expression (Harrington 1968; Docobo 1977; Ford et al. 2000; Naoz et al. 2013): H=Gm1m2 2a1+Gm3(m1+m2) 2a2(3) +G a2 ∞  j=2 jMjr1 a1ja2 r2j+1 Pj(cos), Fig. 3 Evolution of the argument of the pericenter for the example in Fig. 1(exoplanet case, up) and for the example of Fig. 2(exomoon case, down). In both cases, the argument of the pericenter circularizes Fig. 4 Diagram of the Totally Hierarchical Configuration where r1and r2are the radius vectors of the inner and outer orbit, respectively, =a1/a2,Pjis the Legendre polyno- 42 Page 6 of 16 C.V. Monzón, J.A. Docobo Table 3 Exoplanet case. Initial orbital elements of the fourth body. In each case, the orbital elements are referenced to the invariable plane of the system Mass (MJ)T(years) ea(au) i(◦)(◦)ω(◦) Configuration Figure Example 1 40 0 0.2 2500 0.5 30 0 THC 6 Example 2 40 0 0.3 2500 0.5 30 0 THC 7 Example 3 40 0 0.2 2500 20 30 0 THC 8 Example 4 209.4 0 0.2 2500 10 30 0 THC 9 Example 5 40 0 0.1 3 0.5 30 0 2 +213 Example 6 40 0 0.1 5 40 30 0 2 +214 Fig. 5 Hierarchical four-body system with an exoplanet orbiting one component of the inner binary. It corresponds to the case 12 of Campo and Docobo (2014) mial of degree j,is the angle between r1and r2,a1and a2are the inner and outer orbit semi-axis, respectively, and Mj=m1m2m3mj−1 1−(−m2)j−1 (m1+m2)j. As usual, in order to handle the Hamiltonian, we need to truncate the sumatory terms. Truncating them at j=2, we obtain the quadrupole approximation and, at j=3, the octupole approximation. Obviously, the smaller , the smaller will be the difference between the quadrupole or octupole approximation and the entire Hamiltonian. So, in the exomoon case, the contribution of the summands such j≥3 will be negligible and, likewise, the octupole approximation will be very similar to the quadrupole approximation. On the other hand, if is not very small, orbital orientation changes may be absent in the quadrupole approximation but not in the octupole approximation (Naoz et al. 2013). Hence, we deduce that the reason why these changes are non-existing in the exomoon case is because is very small and the contribution of the summands such j>2is not enough to change the orbit orientation. 4 Long-term perturbations in four-body systems Now, we analyze how Lidov-Kozai cycles are perturbed by a fourth body. In order to do that, we add it in both the exoplanet (Fig. 1) and the exomoon case (Fig. 2). In a fourbody system, there are at least two possible stable configurations: the Totally Hierarchical Configuration, where the bodies are organized in nested pairs (Fig. 4), and the 2 + 2 configuration, in which the bodies are organized in two separated pairs, where both revolve around their center of mass. Next, we illustrate each system using the decomposition in subsystems depending on their hierarchy, following the methodology described in Abad and Ribera (1984) and Abad and Docobo (1987). This will help us visualize each configuration. 4.1 Exoplanet case We will assume, for the first three bodies, the same orbital elements as in the example of the Fig. 1and we will study different examples depending on the mass, eccentricity, and orbital inclination of the fourth body. We will assume the mass of the fourth body is equal to that of a low-mass star (a red or brown dwarf, depending of the example). The initial eccentricity will be low to avoid orbital instabilities, between 0.1 and 0.3, as well as the initial inclination, which will be between 0.5◦and 40◦.6The initial semi-axis, in the THC (Totally Hierarchical Configuration), will be chosen to be as close as the stability limit as possible. The initial and ωwill be, in both exoplanet and exomoon cases, always equal to 30◦and 0◦, respectively.7 4.1.1 Totally hierarchical configuration We associate each case with another already covered by Campo and Docobo (2014) and Campo (2019) in which they took into account each and every four-body system where one of the bodies is an exoplanet and/or an exomoon (Fig. 5). From all the cases of study (Table 3), we can state some facts: 6The addition of a massive body such as a red or brown dwarf will change significantly the invariable plane of the examples in Sect. 3. However, we maintain the same initial inclinations with respect to that new plane, so the initial mutual inclinations of the three bodies orbits will be the same as in the Three Body case. 7We took those values arbitrarily as we have verified they have negligible effect in the long-term stability of the bodies orbit. Perturbations in four-body systems Page 7 of 16 42 Fig. 6 Temporal evolution of the eccentricity (top panel) and the orbital inclination (bottom panel) for the exoplanet associated with Table 1 in a totally hierarchical four-body system. The mass and the initial semi-major axis, eccentricity, and inclination of the fourth body are 40 MJ, 2500 au, 0.2, and 0.5 degrees, respectively 1. The examples 1 and 2, in which the initial orbital inclination of the fourth body is 0.5◦(Figs. 6and 7, respectively), still preserve the LK mechanism, as in the TBP: the inclination oscillations cause the orbit to change from prograde to retrograde motion, and vice versa. We notice one difference though: the periods in which these orbital flips happen become dramatically larger or shorter in ∼17 Myr (see Figs. 6 and 7). Likewise, its orbital eccentricity undergoes huge oscillations, periodically approaching 1, and the duration of these periods can increase or decrease over time. We think that the presence of the fourth body is the cause of such behavior but a detailed explanation about this will be explored in future studies. 2. The example 3, in which the initial orbital inclination of the fourth body is 20◦(Fig. 8), still preserves the LK behavior at first, with the exception that, in this case, its inclination undergoes higher amplitude oscillations (30◦to 180◦). However, after 10 Myr, the orientation changes beFig. 7 Temporal evolution of the eccentricity (top panel) and the orbital inclination (bottom panel) for the exoplanet associated with Table 1 in a totally hierarchical four-body system. The mass and the initial semi-major axis, eccentricity, and inclination of the fourth body are 40 MJ, 2500 au, 0.3, and 0.5 degrees, respectively come more chaotic and rapid as its eccentricity constantly oscillates between 0 and 0.9999 in less time. 3. When we modified the initial orbital inclination of the fourth body in both examples 3 and 4 to be greater than 25◦, we verified that the exoplanet orbit become unstable. Therefore, in these examples, the orbital inclination of the fourth body cannot be very large. 4. In example 4, when the mass of the fourth body is 209.4 MJ(red dwarf), the exoplanet orbit becomes unstable when we modify the initial orbital inclination of the fourth body to be 0.5◦. However, when its orbit is initially slightly inclined (between 10◦and 20◦), the exoplanet orbit becomes stable. So, if the initial orbital inclination of the fourth body is 10◦(Fig. 9), the LK mechanism can be noticed in the initial orbital flip (15 Myr). Later on, these flips happen chaotically, without the pattern present in the LK cycles (symmetrical flips with respect to 90◦). 42 Page 8 of 16 C.V. Monzón, J.A. Docobo Fig. 8 Temporal evolution of the eccentricity (top panel) and the orbital inclination (bottom panel) for the exoplanet associated with Table 1 in a totally hierarchical four-body system. The mass and the initial semi-major axis, eccentricity, and inclination of the fourth body are 40 MJ, 2500 au, 0.2, and 20 degrees, respectively 5. With regard to the argument of the pericenter, in all cases, periods of libration and circularization alternate over time. Circularization seems to be more frequent (Figs. 10(a), 10(b), 10(c) and 10(d)). 4.1.2 2 +2configuration In this case, which is illustrated in Figs. 11 and 12, we will assume that the initial exoplanet orbital inclination is 70◦. The rest of the initial orbital and physical parameters will be the same as those of the Three Body case for all three bodies (see Fig. 1).8 8Regarding Figs. 6,7,8and 9, the results must be carefully interpreted since this is a consequence of assuming point-like masses. In fact, if the central star was assumed with a real radius of 0.0046 au, a collision would be very likely. Moreover, the tides and general relativity effects can play a key role in the dynamical behavior of bodies. Fig. 9 Temporal evolution of the eccentricity (top panel) and the orbital inclination (bottom panel) for the exoplanet associated with Table 1 in a totally hierarchical four-body system. The mass and the initial semi-major axis, eccentricity, and inclination of the fourth body are 209.4 MJ, 2500 au, 0.2, and 10 degrees, respectively In examples 5 and 6, where the mass of the fourth body is 40 MJ, something curious happens: the orbital flips observed in the Three Body case are suppressed. This occurs when a=3au,i=0.5◦(Fig. 13), and a=5au,i=40◦ (Fig. 14). Therefore, the presence of a fourth body in these examples severely restricts the oscillation of the exoplanet orbit inclination with respect to the Three Body case. The argument of the pericenter exhibits behavior analogous to the THC configuration in Example 5. However, in Example 6, the argument of the pericenter librates around 270 degrees (Figs. 10(e) and 10(f)). 4.2 Exomoon case For the first three bodies, we will assume the same orbital elements as in Fig. 2. In the exoplanet case, we study different examples depending on the mass, eccentricity and orbital inclination of the fourth body. For the exomoon case, the examples are represented in Table 4. Perturbations in four-body systems Page 9 of 16 42 Fig. 10 Evolution of the argument of the pericenter for the examples of the exoplanet case 42 Page 16 of 16 C.V. Monzón, J.A. Docobo Abad, A., Barrio, R., Blesa, F., et al.: Algorithm 924: TIDES, a Taylor series integrator for differential equations. ACM Trans. Math. Softw. 39(1), 1–28 (2012) Abad, A., Barrio, R., Marco-Buzunariz, M., et al.: Automatic implementation of the numerical Taylor series method a Mathematica and Sage approach. Appl. Math. Comput. 268, 227–245 (2015) Barrio, R., Rodríguez, M., Abad, A., et al.: Breaking the limits: the Taylor series method. Appl. Math. Comput. 217(20), 7940–7954 (2011) Campo, P.P.: Dynamics of exoplanets and exosatellites in binaries. Ph.D. thesis, Universidade de Santiago de Compostela (2019) Campo, P.P., Docobo, J.A.: Analytical study of a four-body configuration in exoplanet scenarios. Astron. Lett. 40(11), 737–748 (2014) Deprit, A.: Canonical transformations depending on a small parameter. Celest. Mech. 1(1), 12–30 (1969) Docobo, J.A.: Aplicación de la teoría de Perturbaciones al estudio de sistemas estelares triples. Ph.D. thesis, Universidad de Zaragoza (1977) Docobo, J.A., Piccotti, L., Abad, A., et al.: A study about the secular evolution of the hierarchical three-body problem using the numerical integrator TIDES. Astron. J. 161(1), 43 (2021) Domingos, R.C., Winter, O.C., Yokoyama, T.: Stable satellites around extrasolar giant planets. Mon. Not. R. Astron. Soc. 373(3), 1227–1234 (2006) Dvorak, R.: Numerical experiments on planetary orbits in double stars. In: The Stability of Planetary Systems, pp. 369–378. Springer, Berlin (1984) Dvorak, R., Froeschlé, C., Froeschle, C.: Stability of outer planetary orbits (P-types) in binaries. Astron. Astrophys. 226, 335–342 (1989) Ford, E.B., Kozinsky, B., Rasio, F.A.: Secular evolution of hierarchical triple star systems. Astrophys. J. 535(1), 385 (2000) Gutzwiller, M.C.: Moon-Earth-Sun: the oldest three-body problem. Rev. Mod. Phys. 70(2), 589 (1998) Hairer, E., Nørsett, S.P., Wanner, G.: Solving Ordinary Differential Equations. 1, Nonstiff Problems, 2nd edn. Springer Series in Comput. Math., vol. 8 (1993) Hamers, A.S., Lai, D.: Secular chaotic dynamics in hierarchical quadruple systems, with applications to hot jupiters in stellar binaries and triples. Mon. Not. R. Astron. Soc. 470(2), 1657–1672 (2017) Hamers, A.S., Perets, H.B., Antonini, F., et al.: Secular dynamics of hierarchical quadruple systems: the case of a triple system orbited by a fourth body. Mon. Not. R. Astron. Soc. 449(4), 4221–4245 (2015) Harrington, R.S.: Dynamical evolution of triple stars. Astron. J. 73, 190–194 (1968) Harrington, R.S.: The stellar three-body problem. Celest. Mech. 1(2), 200–209 (1969) Harrington, R.S.: Encounter phenomena in triple stars. Astron. J. 75, 1140 (1970) Harrington, R.S.: Planetary orbits in binary stars. Astron. J. 82, 753–756 (1977) Kozai, Y.: Secular perturbations of asteroids with high inclination and eccentricity. Astron. J. 67, 591 (1962) Krylov, N.M., Bogoliubov, N.N.: Introduction to Non-linear Mechanics. Princeton University Press, Princeton (1947) Lagrange, J.L.: Théorie des variations séculaires des élémens des planètes. Seconde partie. Contenant la détermination de ces variations pour chacune des planètes principales (1782) Laplace, P.S.: Mémoire sur les inégalités séculaires des planètes et des satellites. Académie royale des sciences (1784) Lidov, M.: Evolution of artificial planetary satellites under the action of gravitational perturbations due to external bodies. Iskusstviennye Sputniki Zemli 8, 5–45 (1961) Mayor, M., Queloz, D.: A Jupiter-mass companion to a solar-type star. Nature 378(6555), 355–359 (1995) Naoz, S., Farr, W.M., Lithwick, Y., et al.: Secular dynamics in hierarchical three-body systems. Mon. Not. R. Astron. Soc. 431(3), 2155–2171 (2013) Pendleton, Y.J., Black, D.C.: Further studies on criteria for the onset of dynamical instability in general three-body systems. Astron. J. 88, 1415–1419 (1983) Rabl, G., Dvorak, R.: Satellite-type planetary orbits in double stars—a numerical approach. Astron. Astrophys. 191, 385–391 (1988) Simó, C.: Global dynamics and fast indicators. In: Global Analysis of Dynamical Systems, pp. 373–389 (2001) Szebehely, V.: Stability of planetary orbits in binary systems. Celest. Mech. 22(1), 7–12 (1980) Takeda, G., Rasio, F.A.: High orbital eccentricities of extrasolar planets induced by the Kozai mechanism. Astrophys. J. 627(2), 1001 (2005) Tokovinin, A., Kiyaeva, O.: Eccentricity distribution of wide binaries. Mon. Not. R. Astron. Soc. 456(2), 2070–2079 (2015) Von Zeipel, H.: Ark. Mat. Astron. Fys. 11(1) (1916) Wolszczan, A., Frail, D.A.: A planetary system around the millisecond pulsar PSR1257 +12. Nature 355(6356), 145–147 (1992) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.