Secular orbital dynamics of the innermost exoplanet of the $\upsilon$-Andromed{\ae} system

We introduce a quasi-periodic restricted Hamiltonian to describe the secular motion of a small-mass planet in a multi-planetary system. In particular, we refer to the motion of $\upsilon$-And $b$ which is the innermost planet among those discovered in the extrasolar system orbiting around the $\upsilon$-Andromedae A star. We preassign the orbits of the Super-Jupiter exoplanets $\upsilon$-And $c$ and $\upsilon$-And $d$ in a stable configuration. The Fourier decompositions of their secular motions are reconstructed by using the Frequency Analysis and are injected in the equations describing the orbital dynamics of $\upsilon$-And $b$ under the gravitational effects exerted by those two external exoplanets (expected to be major ones in such an extrasolar system). We end up with a $2+3/2$ degrees of freedom Hamiltonian model; its validity is confirmed by the comparison with several numerical integrations of the complete $4$-body problem. Furthermore, the model is enriched by taking into account also the relativistic effects on the secular motion of the innermost exoplanet. We focus on the problem of the stability of $\upsilon$-And $b$ as a function of the parameters that mostly impact on its orbit, i.e. the initial values of its inclination and the longitude of its node. We study the evolution of its eccentricity, crucial to exclude orbital configurations with high probability of (quasi)collision with the central star in the long-time evolution of the system. Moreover, we also introduce a normal form approach, that further reduces our Hamiltonian model to a system with $2$ degrees of freedom, which is integrable because it admits a constant of motion related to the total angular momentum. This allows us to quickly preselect the domains of stability for $\upsilon$-And $b$, with respect to the set of the initial orbital configurations that are compatible with the observations.


Introduction
The υ-Andromedae system was the first ever to be discovered among the ones that host at least two exoplanets. In fact, a few years after the discovery of the first exoplanet, the evidence for multiple companions of the υ-Andromedae system was announced (see [21] and [1], respectively). In particular, the observations made with the detection tecnique of the Radial Velocity (hereafter, RV) revealed the existence of orbital objects with three different periods, 4.6, 241 and 1267 days, which revolve around υ-Andromedae A, that is the brightest star of a binary hosting also the red dwarf υ-Andromedae B. Such exoplanets were named υ-And b, c and d in increasing order with respect to their distance from the main star. Since υ-Andromedae B is very far with respect to these other bodies (i.e., ∼ 750 AU), then it is usual to not consider its negligible gravitational effects when the dynamical behavior of the planetary system orbiting around υ-Andromedae A is studied.
None of the present detection methods allows us to know all the orbital parameters of an extrasolar planet. For instance, the RV technique does not provide any information about both the inclination and the longitude of node. In the υ-Andromedae case, these two orbital elements were measured (although with rather remarkable error bars) for both υ-And c and υ-And d thanks to observational data taken from the Hubble Space Telescope (see [22]). The information provided by such an application of astrometry significantly complemented the knowledge about this extrasolar system; in fact, it has led to the evaluation of the masses of υ-And c and υ-And d (ranging in 13.98 +2.3 −5.3 M J and 10.25 +0.7 −3.3 M J , respectively) and of their mutual inclination (29.9 • ±1 • ). It is well known that only minimum limits for the masses can be deduced by observations made with the RV method (due to the intrinsic limitations of such a technique). Moreover, the mutual inclination between planetary orbits plays a crucial role for what concerns the stability of extrasolar systems (see, e.g., [28]). Thus, the orbital configuration of υ-Andromedae is probably one of the most accurately known among the extrasolar multi-planetary systems which have been discovered so far.
The question of the orbital stability of υ-Andromedae is quite challenging. Numerical integrations revealed that unstable orbits are frequent. Moreover, these extensive explorations allowed to locate four different regions of initial values of the orbital parameters (consistent with all observational constraints) yielding dynamically stable orbital configurations for the three planets of the system (see [5]). All these sets of parameters correspond to values of the mass of υ-And c that are relatively small, in the sense that they are much closer to the lower bound of the range 13.98 +2.3 −5.3 M J than to the upper one. On the other hand, according to a numerical criterion inspired from normal form theory and introduced in [18], the most robust orbital configurations correspond to the largest possible value of the mass of υ-And c in the above range. In the vicinity of the initial conditions giving rise to the most robust orbital configuration, the existence of KAM tori for the dynamics of a secular three-body problem including υ-And c and υ-And d was proved in a rigorous computer-assisted way (see [3]).
In the present paper we aim to extend the study of the stability to the orbital dynamics of υ-And b, still adopting a hierarchical approach. In the case of the particular extrasolar system under consideration, this means that we assume the mass of υ-And b so small 3 with respect to the ones of υ-And c and υ-And d, that the motion of the innermost planet b can be modeled with a good approximation via a restricted four-body problem. More precisely, in order to study the dynamical behavior, we preassign the secular motions of the Super-Jupiter exoplanets c and d in correspondence to the quasi-periodic orbit that is expected to be the most robust. The Fourier decompositions of the secular motions of c and d are reconstructed by using the well known technique of the (so called) Frequency Analysis (see e.g. [14]) and are injected in the equations describing the orbital dynamics of υ-And b, under the gravitational effects exerted by those two external exoplanets. This way to introduce a quasi-periodic restricted model has been recently used to study the long-term dynamics of our Solar System (see [24] and [12]).
The advantage of introducing a secular quasi-periodic restricted Hamiltonian looks evident. In our present case (referring to the υ-Andromedae system with planets b , c , d ), we start with a Hamiltonian model having 9 degrees of freedom, ending up with a simpler one with 2 + 3/2 degrees of freedom, where the short periods are dropped. This explains why numerical explorations of the restricted Hamiltonian model are much faster. Our main purpose is further promoting this procedure, in such a way to introduce a new simplified model with just two degrees of freedom. Indeed, in the present paper, we show that this can be done by adopting a suitable normal form approach, which is so accurate to produce an integrable Hamiltonian that can be used to efficiently characterize the stability domain with respect to the unknown orbital parameters of υ-And b, i.e., the inclination and the longitude of the node.
The present work is organized as follows. In Section 2, Frequency Analysis is used to reconstruct the Fourier decompositions of the secular motions of the outer exoplanets υ-And c and υ-And d. In Section 3, the secular quasi-periodic restricted Hamiltonian model (with 2 + 3/2 degrees of freedom) is introduced and validated through the comparison with several numerical integrations of the complete four-body problem, hosting planets b , c , d of the υ-Andromedae system. The double normalization procedure allowing to perform a sort of averaging which further simplifies the model is described with a rather general approach in Section 4. In Section 5 this normal form procedure is applied to the quasi-periodic restricted Hamiltonian, in such a way to derive an integrable model with 2 degrees of freedom describing the secular orbital dynamics of υ-And b . Such a simplified model is used to study υ-And b stability domain in the parameters space of the initial values of the inclination and the longitude of node. All this computational procedure is repeated in Section 6 starting from a version of the secular quasi-periodic restricted Hamiltonian model which includes also relativistic corrections; this allows us to appreciate the effects on the orbital dynamics due to General Relativity.

Determination of the outer planets motion
To prescribe the orbits of the giant planets υ-And c and υ-And d, we start from the Hamiltonian of the three-body problem (hereafter, 3BP) in Poincaré heliocentric canonical variables, using the formulation based on the reduced masses β 2 , β 3 , that is where m 0 is the mass of the star, m j , r j , p j , j = 2, 3, are the masses, astrocentric position vectors and conjugated momenta of the planets, respectively, G is the gravitational constant and β j = m 0 m j /(m 0 + m j ), j = 2, 3, are the reduced masses. Let us remark that, in the following, we use the indexes 2 and 3 respectively, for the inner (υ-And c) and outer (υ-And d) planets between the giant ones, while the index 1 is used to refer to υ-And b In order to set up a quasi-periodic restricted model for the description of the motion of υ-And b, we need to characterize the motion of the giant planets; this can be done through the Frequency Analysis method, starting from the numerical integration of the complete 3BP Hamiltonian, reported in equation (1) (i.e., before any expansions and averaging 4 ). Thus, we numerically integrate the complete Hamiltonian (1) using a symplectic method of type SBAB 3 , which is described in [16]. As initial orbital parameters for the outer planets, we adopt those reported in Table 1, corresponding to the most robust planetary orbit compatible with the observed data available for υ-And c and υ-And d (see [22]), according to the criterion of "minimal area" explained in [18]. 229.325 7.374 Table 1: Chosen values of the masses and of the initial orbital parameters for υ-And c and υ-And d, compatible with the observed data available, as reported in [22].
Having fixed as initial orbital parameters the ones described in Table 1, it is possible to compute their correspondent values in the Laplace reference frame (i.e., the invariant reference frame orthogonal to the total angular momentum vector r 2 × p 2 + r 3 × p 3 ) and to perform the numerical integration of the full 3BP corrisponding to these initial values. Then, it is possible to express the discrete results produced by the numerical integrations in the canonical Poincaré variables (ξ j , η j ), (P j , Q j ) (momenta-coordinates) given by 5 where Λ j = β j √ µ j a j , β j = m 0 m j /(m 0 + m j ) , µ j = G (m 0 + m j ) , and e j , i j , ω j , Ω j , j = ω j + Ω j refer, respectively, to the eccentricity, inclination, argument of the periastron, longitudes of the node and of the periastron of the j-th planet. However, the numerical integrations do not allow to obtain a complete knowledge of the motion laws t → (ξ j (t), η j (t)), t → (P j (t), Q j (t)) (j = 2, 3), producing only discrete time series made by sets of finite points computed on a regular grid in the interval [0, T ] . The computational method of Frequency Analysis (hereafter, FA) allows however to reconstruct with a good accuracy the motion laws by using suitable continuous in the time variable t functions. This has been done recently in [24] and [12], in order express the motion of the Jovian planets of our Solar System as a Fourier decomposition including just a few of the main terms. In the present Section, we basically follow that approach; through a canonical transformation, corresponding to a second order (in the mass ratios) averaging (see [19]). However, we have observed that the Fourier decomposition given by the second order in masses numerical integration is not enough accurate for such a kind of model and that a more detailed approximation of the orbits of the outer giant planets υ-And c, υ-And d is required. 5 The definition of the Poincaré variables (ξ1, η1), (P1, Q1) will be useful in the following Sections.
therefore, here we limit ourselves to report some definitions which are essential in order to make our computational procedure well definite (see, e.g., [14] for an introduction and a complete exposition about FA). We consider analytic quasi-periodic motion laws t → z(t). This means that the function z : R → C admits the following Fourier series decomposition: where ω ∈ R n is the so called fundamental angular velocity vector, while a k ∈ R + ∪ {0} and ϑ k ∈ T , ∀ k ∈ Z n ; moreover, the sequence {a k } k∈Z n is assumed to satisfy an exponential decay law, i.e., |a k | ≤ c e −|k|σ ∀ k ∈ Z n with c and σ real positive parameters. The FA allows us to find an approximation of z(t) of the following form: where N C is the number of components we want to consider. The equation (4) is an approximation of the motion z(t) in the sense that if N C → +∞ and T → +∞ , the right hand side of (4) converges to (3). Moreover, a s;T ∈ R + and ϑ s;T ∈ [0, 2π) are called respectively the amplitude and the initial phase of the s-th component, while ν where W is a suitable weight function such that T 0 W(t) dt = T . Following [14], we use the so called "Hanning filter", defined (in [0, T ]) as W(t) = 1 − cos[π(2t/T − 1)] .
The numerical integration of the 3BP (equation (1)) produces a discretization of the signals 6 t → ξ j (t) + iη j (t) and t → P j (t) + iQ j (t) , that are {(ξ j (s∆), η j (s∆))} N P s=0 and {(P j (s∆), Q j (s∆))} N P s=0 (j = 2, 3), where s = 0, . . . , N P and the sampling time is ∆ = T /N P . These discretizations allow to compute (by numerical quadrature) the integral in (5) and, consequently, a few of local maximum points of the function (5) considering, as z(t), ξ 2 (t) + iη 2 (t) , ξ 3 (t) + iη 3 (t) , P 2 (t) + iQ 2 (t) and P 3 (t) + iQ 3 (t) . Then, we use the FA to compute a quasi-periodic approximation of the secular dynamics of the giant planets υ-And c and υ-And d, i.e., j = 2, 3 , where the angular vector 6 Actually, the numerical integration of the complete problem allows to determine also a discretization of the fast variables √ 2Λ2 cos(λ2) + i √ 2Λ2 sin(λ2) and √ 2Λ3 cos(λ3) + i √ 2Λ3 sin(λ3) ; however we are not interested in these variables. Thus, we do not report their decompositions as they are provided by the FA. depends linearly on time and ω ∈ R 3 is the fundamental angular velocity vector whose components are listed in the following: Hereafter, the secular motion of the outer planets t → (ξ j (t), η j (t), P j (t), Q j (t)), j = 2, 3, is approximated as it is written in both the r.h.s. of formula (6). The numerical values of the coefficients which appear in the quasi-periodic decompositions 7 of the motions laws are reported in Tables 2, 3, 4 and 5.      In order to verify that the numerical solutions are well approximated by the quasi-periodic decompositions computed above, we compare the time evolution of the variables ξ 2 , ξ 3 , η 2 , η 3 , P 2 , P 3 , Q 2 , Q 3 as obtained by the numerical integration and by the FA. Figure 1 shows that the quasi-periodic approximations nearly perfectly superpose to the plots of the numerical solutions.

The secular quasi-periodic restricted (SQPR) Hamiltonian
Having preassigned the motion of the two outer planets υ-And c and υ-And d , it is now possible to properly define the secular model for a quasi-periodic restricted four-body problem (hereafter, 4BP). 7 Of course, the exact quasi-periodic decompositions include infinite terms in the Fourier series. In order to reduce the computational effort, we limit ourselves to consider just a few components, which are the main and most reliable ones, according to the following criteria. We take into account those terms corresponding to low order Fourier armonics, i.e., 5 j=3 |kj| ≤ 5 or 5 j=3 |kj| ≤ 5 , and such that there are small uncertainties on the determination of the frequencies as linear combinations of the fundamental ones, i.e., |ν (their definition is reported in (2)) as it is computed by numerical integration of the complete (non secular) three-body problem and by the quasi-periodic approximation provided by the FA; the corresponding plots are in blue and red, respectively.
We start from the Hamiltonian of the 4BP, given by We recall that the so called secular model of order one in the masses is given by averaging with respect to the mean motion angles, i.e., In the l.h.s. of the equation above, we disregard the dependence on the actions Λ, because in the secular approximation of order one in the masses their values Λ j = β j √ µ j a j , j = 1, 2, 3, are constant.
Due to the d'Alembert rules (see, e.g., [26] and [25]), it is well known that the secular Hamiltonian can be expanded in the following way: where N is the order of truncation in powers of eccentricity and inclination. We fix N = 8 in all our computations.
Since we aim at describing the dynamical secular evolution of the innermost planet υ-And b , it is sufficient to consider the interactions between the two pairs υ-And b, υ-And c and υ-And b, υ-And d.
In more details, let H i−j sec be the secular Hamiltonian derived from the three-body problem for the planets i and j, averaging with respect to the mean longitudes λ i , λ j . Its expansion writes as Therefore, a restricted non-autonomous model which approximates the secular dynamics of υ-And b can be defined by considering the terms where ξ 2 (t), η 2 (t), . . . P 3 (t), Q 3 (t) are replaced with the corresponding quasi-periodic approximations written in both the r.h.s. appearing in formula (6). Let us stress that, having prescribed the motion of the two outermost planets υ-And c and υ-And d , at this stage the Hamiltonian H 2−3 sec does not need to be reconsidered; indeed, it would introduce additional terms that disappear in the equations of motion (see formula (15) which is written below).
We can finally introduce the quasi-periodic restricted Hamiltonian model for the secular dynamics of υ-And b ; it is given by the following 2 + 3/2 degrees of freedom Hamiltonian: where the pairs of canonical coordinates referring to the planets υ-And c and υ-And d (that are ξ 2 , η 2 , . . . P 3 , Q 3 ) are replaced by the corresponding finite Fourier decomposition written in formula (6) as a function of the angles θ, renamed 8 as q, i.e., q = (q 3 , q 4 , q 5 ) := (θ 3 , θ 4 , θ 5 ) = θ .
Let us focus on the summands appearing in the first row of (13), i.e., the Hamiltonian term ω · p , where ω is the fundamental angular velocity vector (defined in formula (8)) and p = (p 3 , p 4 , p 5 ) is made by three so called "dummy variables", which are conjugated to the angles q. The role they play is made clear by the equations of motion for the innermost planet, which write in the following way in the framework of the restricted quasi-periodic secular approximation: Due to the occurrence of the term ω · p in the Hamiltonian H sec, 2+ 3 2 , the first three equations admit q(t) = ωt as a solution, in agreement with formulae (7) and (14). This allows to reinject the wanted quasi-periodic time-dependence in the Fourier approximations ξ 2 (q), η 2 (q), . . . P 3 (q), Q 3 (q). As a matter of fact, we do not need to compute the evolution of (p 3 (t), p 4 (t), p 5 (t)) because they do not exert any influence on the motion of υ-And b ; they are needed just if one is interested in checking that the energy is preserved, because it is given by the evaluation of H sec, 2+ 3 2 .
We also recall that, in order to produce a restricted quasi-periodic secular model, it is possible to apply the closed-form averaging, which is compared in [20] with the computational method that is adopted here and is based on the expansions in Laplace coefficients. Finally, we emphasize what is discussed below. Thus, it admits a constant of motion that could be reduced, so to decrease 9 by one the number of degrees of freedom of the model.
In order to clarify the statement above, it is convenient to introduce a complete set of action-angle variables, defining two new pairs of canonical coordinates referring to a pair of orbital angles of υ-And b, i.e., 1 and Ω 1 , that are the longitudes of the pericenter and of the node, respectively (see definition (2) of the Poincaré canonical variables). Thus, it is possible to verify the following invariance law: Therefore, p 3 + p 4 + p 5 + Γ 1 + Θ 1 is preserved; such a quantity, apart from an inessential additional constant, is equivalent to the total angular momentum. The above invariance law is better understood, recalling that q 3 and q 4 correspond to the longitudes of the pericenters of υ-And c and υ-And d, respectively, while q 5 refers to the longitude of the nodes of υ-And c and υ-And d (that in the Laplace frame, determined by taking into account only these two exoplanets, are opposite one to the other). This identification is due to the way we have determined (q 3 , q 4 , q 5 ) by decomposing some specific signals of the secular dynamics of the outer exoplanets (this is made by using the Frequency Analysis, as it is explained in Section 2). Thus, the aforementioned invariance law is due to the fact that the dynamics of the model we are studying does depend just on the pericenters arguments of the three exoplanets and on the difference between the longitude of the nodes of υ-And b and υ-And c, i.e., Ω 1 − Ω 2 = Ω 1 − Ω 3 − π. Therefore, the Hamiltonian is invariant with respect to any rotation of the same angle that is applied to all the longitudes of the nodes; as it is well known, by Noether theorem, this is equivalent to the preservation of the total angular momentum.

Numerical validation of the SQPR model
In order to validate our secular quasi-periodic restricted (hereafter SQPR) model describing the dynamics of υ-And b, we want to compare the numerical integrations of the complete 4BP with the ones of the equations of motion (15). Let us recall that the chosen values of parameters and initial conditions for the two outer planets are given in Table 1. For what concerns the orbital elements of the innermost planet υ-And b, both the inclination i 1 and the longitude of the ascending node Ω 1 are unknown (see, e.g., [22]). The available data for υ-And b are reported in Table 6. Among the possible values of the initial orbital parameters of υ-And b, we have chosen a 1 , e 1 , M 1 and ω 1 as in the stable prograde trial PRO2 described in [5]. They are reported in Table 7 and are compatible with the available ranges of values appearing in Table 6. Let us recall that the dynamical evolution of the SQPR model does not depend on the mass of υ-And b, therefore, the choice about its value is not reported in Table 7. For what concerns the unknown initial values of the inclination and of the longitude of nodes, we have decided to vary them so as to cover a 2D regular grid of values , dividing I i and I Ω , respectively, in 20 and 60 sub-intervals; this means that the widths of the grid-steps are equal to 1.35675 • and 6 • in inclination and longitude of nodes, respectively. Let us recall that the lowest possible value of the interval I i , i.e. i 2 (0) = 6.865 • , corresponds to the inclination of υ-And c. Considering values smaller than i 2 (0) could be incoherent with the assumptions leading to the SPQR model we have just introduced; indeed, the factor 1/ sin(i 1 (0)) increases the mass of the exoplanet by one order of magnitude with respect to the minimal one. Therefore, for small values of i 1 (0) the mass of υ-And b could become so large that the effects exerted by its gravitational attraction on the outer exoplanets could not be neglected anymore. On the other hand, it will be shown in the sequel that the stability region for the orbital motion of υ-And b nearly completely disappears for values of i 1 (0) larger than 34 • . These are the reasons behind our choice about the lower and upper limits of I i .
44.106 ± 25.561 Table 6: Available data for the orbital parameters of the exoplanet υ-And b. The values above are reported from Table 13 of [22].
51.14 Table 7: Values of the initial orbital parameters for υ-And b as they have been selected in the stable prograde trial PRO2 of [5] (Table 3 ).
We emphasize that the study of the stability domain, as it is deduced by the numerical integrations, can help us to obtain information about the possible ranges of the unknown values (i 1 (0), Ω 1 (0)). Moreover, the comparisons between the numerical integrations of the complete 4BP and the ones of the SQPR model aim to demonstrate that the agreement is sufficiently good so that it becomes possible to directly work with the latter Hamiltonian model, that has to be considered easier than the former, because the degrees of freedom are 2 + 3/2 instead of 9.

Numerical integration of the complete 4-body problem
For each pair of values (i 1 (0) , Ω 1 (0)) ∈ I i × I Ω ranging in the 20 × 60 regular grid we have previously prescribed, we first compute the corresponding initial orbital elements of the three exoplanets in the Laplace-reference frame, then we perform the numerical integration of the complete 4BP Hamiltonian (9) by using the symplectic method of type SBAB 3 . Contrary to the SPQR model, the numerical integrations of the 4BP are affected by the mass of υ-And b; its value is simply fixed in such a way that m 1 = 0.674/ sin(i 1 (0)) . The largest value reached by the eccentricity can be considered as a very simple numerical indicator about the stability of the orbital configurations. The maximum eccentricity obtained along each of the 21 × 60 numerical integrations is reported in the left panel of Figure 2. In particular, since we are interested in initial conditions leading to regular behavior, i.e., avoiding quasi-collisions, every time that the eccentricity e 1 exceeds a threshold value (fixed to be equal to 0.85), in the color-grid plots its maximal value is arbitrarily set equal to one. Moreover, since we expect that υ-And c is the most massive planet in that extrasolar system and being it the closest one to υ-And b, it is natural to focus the attention also on the mutual inclination between υ-And b and υ-And c . Let us recall that it is defined in such a way that cos(i mut bc ) = cos(i 1 ) cos(i 2 ) + sin(i 1 ) sin(i 2 ) cos(Ω 1 − Ω 2 ) ; (16) therefore, for each numerical integration it is also possible to compute the maximal value reached by i mut bc . The results are reported in the right panel of Figure 2. In both those panels the color-grid plots are provided as functions of the initial values of the longitude of nodes Ω 1 and the inclination i 1 , which are reported on the x and y axes, respectively. By comparing the two plots in the Figure panels 2a and 2b, one can easily appreciate that the regions which have to be considered as dynamically unstable, because the eccentricity of e 1 can grow to large values, correspond also to large mutual inclinations of the planetary orbits of υ-And b and υ-And c.
We remark that the value of the initial mean anomaly M 1 (0) is missing among the available observational data reported in Table 6. As a matter of fact, mean anomalies of exoplanets are in general so poorly known that usually their values are not reported in the public databases. 10   Moreover, the same conclusion applies also to the increasing factor 1/ sin(i 1 (0)) (with i 1 (0) ∈ I i ) which multiplies the minimal mass of υ-And b in such a way to determine the value of m 1 . In fact, substantial differences are not observed between Figures 2 and 4.

Numerical integration of the secular quasi-periodic restricted model
We want now to compare the previous results with those found in the SQPR approximation of the 4-body problem, performing numerical integrations of the system of equations (15). In order to make  Figure 4: Color-grid plots of the maximal value reached by the eccentricity of υ-And b (on the left) and by the mutual inclination between υ-And b and υ-And c (on the right). The maxima are computed during the symplectic numerical integrations of the 4BP which cover a timespan of 10 5 yr. As the only difference with respect to the numerical integration whose results are reported in Figure 2, here the mass of υ-And b is always kept fixed so as to be equal to its minimal value m 1 = 0.674 . these comparisons coherent, also here we consider the data listed in Table 7 as initial conditions for the orbital elements of υ-And b which are completed with the values of (i 1 (0) , Ω 1 (0)) ranging in the 20×60 regular grid that covers At the beginning of the computational procedure, the initial values of the orbital elements are determined in the Laplace reference frame, which is fixed by taking into account only the two outermost planets (i.e., the total angular momentum of the system is given only by the sum of the angular momentum of υ-And c and υ-And d ). Of course, this is made in agreement with our choice to consider a restricted framework, because we are assuming that the mass of υ-And b is so small that can be neglected.
For each numerical integration we compute the maximal value reached by the eccentricity e 1 and the mutual inclination i mut bc . The results are reported in the color-grid plots of the left and right panels of Figure 5, respectively. Once again, they are provided as functions of Ω 1 (0) and i 1 (0) , whose values appear on the x and y axes, respectively.
Comparing Figures 2a with 5a and 2b with 5b, respectively, we can immediately conclude the striking similarity of the color-grid plots, implying the same dependence of the dynamics on the initial values of the orbital elements i 1 (0) and Ω 1 (0) in either model. In particular, the regions of stability located at the two lateral sides of the plots, where the orbit of υ-And b does not become very eccentric, are identical. This occurs also for what concerns the plots of the maximal mutual inclination. However, some discrepancies are evident in the central parts of the panels, i.e. for values of Ω 1 (0) ranging between 90 • and 270 • . We stress that this lack of agreement between the results provided by the two models is expected in these central regions of the panels. Indeed, let us recall that the SQPR model has been introduced starting from some classical expansions in powers of eccentricities and inclinations. Therefore, it is reasonable to expect a deterioration of the accuracy of the SQPR model in the orbital dynamics depicted in the central regions of the plots where large values of the eccentricity e 1 and the mutual inclination are attained. We emphasize that similar remarks about the very strong impact of the initial value Ω 1 (0) on the orbital stability of υ-And b can be found in Section 4.2 of [27].
A further exploration of the stable and chaotic regions of Figure 5a can be done by applying the so called Frequency Map Analysis method (see, e.g., [13]), in order to study the signal ξ 1 (t) + iη 1 (t) produced by the numerical integration of the system (15), i.e., in the SQPR approximation. We perform the numerical integrations as prescribed at the beginning of the present Section, taking into account only a few values in I i for the initial inclinations, i.e., i 1 (0) = 6.865 • , 8.22175 • , 9.5785 • , 10.9353 • and Ω 1 (0) ∈ I Ω . In Figure 6 we report the behaviour of the angular velocity corresponding to the first component of the approximation of ξ 1 (t) + iη 1 (t), as obtained by applying the FA computational algorithm; therefore, this quantity is related to the precession rate of 1 . As initial value for the inclination i 1 (0) we fix 6.865 • for Figure 6a and 10.9353 • for Figure 6b. We do not report the cases (i 1 (0), Ω 1 (0)) ∈ {8.22175 • , 9.5785 • } × I Ω , since the behaviour of those plots is similar to the ones in Figure 6. The situation is well described in Figure 6b and analogous considerations can be done for Figure 6a. For what concerns the values of Ω 1 (0) in the range [0, ∼ 50 • ] and [∼ 325 • , 360 • ] we can observe a regular behaviour of the angular velocity ν which is also monotone with the only exception of the local minimum. According to the interpretation of the Frequency Map Analysis (in the light of KAM theory), such a regular regime is due to the presence of many invariant tori which fill the stability region located at the two lateral sides of the plot 5a. Instead for values of we observe a strongly irregular behaviour, which corresponds to the lateral green stripes and the internal green region of Figure 5a. Thus, they represent chaotic regions in proximity of a secular resonance. Indeed, in Figure 6b Figure 5a). More precisely, the value of ν is equal to −1.04 · 10 −3 , that is ω 4 , i.e., one of the fundamental angular velocities which characterize the quasi-periodic motion of the outer planets  Figure 6: Behaviour of the fundamental angular velocity ν as obtained by applying the Frequency Map Analysis method to the signal ξ 1 (t) + iη 1 (t) , computed through the RK4 numerical integration of the SQPR model (15), covering a timespan of 1.31072 · 10 5 yr. We take, as initial conditions, (i 1 (0), Ω 1 (0)) ∈ {6.865 • } × I Ω for the left panel and (i 1 (0), Ω 1 (0)) ∈ {10.9353 • } × I Ω for the right one.
(see Eq. (8)). This allows us to conclude that they represent the stable central part of a resonant region.

Introduction of a secular model by a normal form approach
This Section aims at manipulating the Hamiltonian with normal form algorithms in order to define a new model that is more compact; this allows us to simulate the secular dynamics of υ-And b with much faster numerical integrations. In fact, we describe a reduction of the number of degrees of freedom (DOF) of our Hamiltonian model. For such a purpose, we apply two normal form methods: first, we perform the construction of an elliptic torus, hence, we proceed removing the angles (q 3 , q 4 , q 5 ) whose evolution is linearly depending on time. The latter elimination is made by applying a normalization methodà la Birkhoff in such a way to introduce a so called resonant normal form 11 that includes, at least partially, the long-term effects due to the outer planets motion.

Normal form algorithm constructing elliptic tori
In [9] the existence of invariant elliptic tori in 3D planetary problems with n bodies has been proved by using a normal form method which is explicitly constructive. However, such an approach does not look suitable to be directly applied to Hamiltonian secular models, because in this latter case the separation between fast and slow dynamics is lost. Therefore, we follow the explanatory notes [17], where the algorithm constructing the normal form for elliptic tori is compared with the classical onè a la Kolmogorov, which is at the basis of the original proof scheme of the KAM theorem. We first summarize this general procedure leading to the construction of elliptic tori. We then add some comments explaining how this general method can be suitably adapted to our problem.
We start considering a Hamiltonian H (0) written as follows: where E (0) is a constant term, representing the energy, (p, q) ∈ R n 1 × T n 1 , (I, α) ∈ R n 2 ≥0 × T n 2 are action-angle variables and (ω (0) , Ω (0) ) ∈ R n 1 × R n 2 is the angular velocity vector. The symbol f (r,s) l is used to denote a function of the variables (p, q, I, α) , such that l is the total degree in the square root of the actions (p , I), s is the index such that the maximum trigonometric degree, in the angles (q, α) , is sK (for a fixed positive integer K) and r refers to the normalization step. In more details, we can say that f (0,s) l ∈ P l,sK , which is a class of functions that we introduce as follows. A few remarks about the above definition are in order. First, since we deal with real functions, the complex coefficients must be such that c m,l,−k,− l =c m,l,k, l . Moreover, the rules about the integer coefficients vector l are such that, ∀ j = 1, . . . , n 2 , the j-th component of the Fourier harmonic l j (that refers to the angle α j ) must have the same parity with respect to the corresponding degree l j of I j and must satisfy the inequality 12 | l j | ≤ l j . Let us here emphasize that our SQPR model of the secular dynamics of υ-And b can be reformulated in such a way to be described by a Hamiltonian of the type (17); this will be explained in detail at the beginning of Section 5.
The following statement plays a substantial role, since it ensures that the structure of the functions f (r,s) l is preserved while the normalization algorithm is iterated.
Lemma 4.1. Let us consider two generic functions g ∈ P l,sK and h ∈ P m,rK , where K is a fixed positive integer number. Then where we trivially extend the definition 4.1 in such a way that P −2, sK = P −1, sK = {0} ∀ s ∈ N.
The algorithm constructing the normal form for elliptic tori is applied to Hamiltonians of the type (17), where the terms appearing in the second row (namely, s≥1 2 l=0 f (0,s) l (p, q, I, α) ) are considered as the perturbation to remove. Therefore, one can easily realize that such a perturbation must be sufficiently small so that the procedure behaves well as regards convergence. There are general situations where this essential smallness condition is satisfied. For instance, this occurs for Hamiltonian systems in the neighborhood of a stable equilibrium point; in fact, it is possible to prove that, f (0,s) l = O(ε s ) , where ε is a small parameter which denotes the first approximation of the distance (expressed in terms of the actions) between the elliptic torus and the stable equilibrium point. The elimination of the small perturbing terms can be done through a sequence of canonical transformations, leading the Hamiltonian in the following final form: with f (∞,s) l ∈ P l,sK . Therefore, for any initial conditions of type (0,q 0 , 0,α) (whereq 0 ∈ T n 1 and the value ofα ∈ T n 2 does not play any role 13 ), the motion law (p(t),q(t),Ĩ(t),α(t)) = (0,q 0 +ω (∞) t, 0,α) is a solution of the Hamilton's equations related to H (∞) . This quasi-periodic solution (having ω (∞) as constant angular velocity vector) lies on the n 1 -dimensional invariant torus such that the values of the action coordinates arep = 0,Ĩ = 0 .
The generic r-th step of the algorithm is defined as follows. Let us assume that after r − 1 normalization steps the expansion of the Hamiltonian can be written as with f (r−1,s) l ∈ P l,sK . By comparing formula (17) with (19), one immediately realizes that the assumption above is satisfied in the case with r = 1 for what concerns the expansion of the initial Hamiltonian H (0) .
The r-th normalization step consists of three substeps, each of them involving a canonical transformation which is expressed in terms of the Lie series having χ 2 as generating function, respectively. Therefore, the new Hamiltonian that is introduced at the end of the r-th normalization step is defined as follows: where exp (L χ ) · = j≥0 (L j χ ·)/j! is the Lie series operator, L χ · = {·, χ} is the Lie derivative with respect to the dynamical function χ , and {·, ·} represents the Poisson bracket. 13 Indeed, whenĨ = 0 ∀α ∈ T n 2 , the canonical coordinates ( 2Ĩj cos(αj) , 2Ĩj sin(αj)) are mapped into the origin of the j-th subspace that is transversal to the elliptic torus. This fictitious singularity of the action-angle variables (I, α) is completely harmless just because all the normalization algorithm can be performed working on Hamiltonians whose expansions are made by terms belonging to sets of functions of type P l,sK . We stress that all the algorithm could be reformulated using polynomial canonical coordinates to describe the dynamics in the subspaces transversal to the elliptic torus; in particular, this is done with complex pairs of canonical coordinates in [2]. In the sequel, we adopt an exposition entirely based on the use of the action-angle coordinates, which makes the algorithm easier to understand.

First substep (of the r-th normalization step)
The first substep aims to remove the term depending only on the angles 14 q up to trigonometric degree rK , i.e., f (r−1,r) 0 (included in the first sum of the second row of (19)), which has to be considered as O(ε r ) . The first generating function χ (r) 0 (q) is determined by solving the following homological equation: Since f (r−1,r) 0 ∈ P 0,rK , its Fourier expansion can be written f Because of the homological equation (21), we find the above solution is well defined if the non-resonance condition is satisfied. We can now apply the Lie series operator exp L χ (r) 0 to H (r−1) . This allows us to write the expansion of the new intermediate Hamiltonian as follows: where (by abuse of notation) for the new canonical coordinates we adopt the same symbols as the old ones. From a practical point of view, the new Hamiltonian terms can be conveniently defined in such a way to mimic what is usually done in any programming language. First, we introduce the new summands as the old ones, so that f (I; r,s) l = f (r−1,s) l ∀ l ≥ 0 , s ≥ 0 . Hence, each term generated by Lie derivatives with respect to χ (r) 0 is added to the corresponding class of functions. By a further abuse of notation, this is made by the following sequence 15 of redefinitions: where with the notation a ← b we mean that the quantity a is redefined so as to be equal a + b . In fact, since χ (r) 0 ∈ P 0,rK , Lemma 4.1 ensures that each application of the Lie derivative operator 14 This first substep of the algorithm is basically useless when the explicit construction of the normal form related to an elliptic torus is started from the Hamiltonian H sec, 2+3/2 described in (51). Indeed, in the case under study just the so called dummy actions are affected by this kind of canonical change of variables, which is defined by a Lie series with a generating function depending on the angles q only. Aiming to make a rather general discussion of the computational procedure, we keep in the algorithm the description of this first normalization substep. 15 From a practical point of view, since we have to deal with finite series, that are truncated in such a way that the index s goes up to a fixed order called NS, we have to require also that 1 ≤ j ≤ min l/2 , (NS − s)/r . L χ (r) 0 decreases by 1 the degree in p (that is obviously equivalent to 2 in the square root of the actions), while the trigonometrical degree in the angles q is increased by rK . By using repeatedly such a simple rule, one can easily verify that f (I; r,s) l ∈ P l,sK ∀ l ≥ 0, s ≥ 0 . Moreover, due to the homological equation (21), we set f (I; r,r) 0 = 0 and update the energy value in such a way that , where · q is used to denote the angular average with respect to q.
Second substep (of the r-th normalization step) The second substep aims to remove the term that is linear in √ I and independent on p , i.e. f (I;r,r) 1 , which is included in the second sum appearing in the second row of (23). The second generating function χ (r) 1 (q, I, α) is determined solving the following homological 16 equation: Since f (I; r,r) 1 ∈ P 1,rK , we can write its expansion as due to the homological equation (25), we find (26) The above expression is well defined provided that the first Melnikov non-resonance condition is satisfied, i.e., min 0<|k|≤rK−1 |l|=1 k · ω (r−1) + l · Ω (r−1) ≥ γ (rK) τ and min |l|=1 l · Ω (r−1) ≥ γ , for a pair of fixed values of γ > 0 and τ > n 1 − 1 (see [17] and reference therein).
We can now apply the transformation exp L χ (r) 1 to the Hamiltonian H (I; r) . By the usual abuse of notation (i.e., the new canonical coordinates are denoted with the same symbols of the old ones), the expansion of the new Hamiltonian can be written as  16 In the r.h.s. of (25) we do not need to put any term produced by an angular average (similar to that appearing, for instance, in the r.h.s. of the homological equation (21)), because f all the terms include the dependence on e ±iα j with j = 1, . . . , n2 , leading to a null mean over the angles.
where in the last row of the previous formula, it is possible to start the first sum from r + 1 instead of r , being f (II; r,r) 0 = f (I; r,r) 0 = 0 . In an analogous way as in the first substep, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., f (II; r,s) l = f (I; r,s) l ∀ l ≥ 0 , s ≥ 0 . Hence, each term generated by the Lie derivatives with respect to χ (r) 1 is added to the corresponding class of functions. This is made by the following sequence 17 of redefinitions: In fact, since χ

Third substep (of the r-th normalization step)
The last substep aims to remove the term f (II; r,r) 2 which is quadratic in the square root of the actions (i.e., either quadratic in √ I or linear in p ) and included in the third sum appearing in the second row of (28). The third generating function χ i k · ω (r−1) + l · Ω (r−1) , provided that both the non-resonance condition and the Melnikov one of second kind are satisfied, i.e., with the same values of the constant parameters γ > 0 and τ > n 1 − 1 appearing in (27).
We can now apply the transformation exp L χ (r) 2 to the Hamiltonian H (II; r) . By the usual abuse of notation (i.e., the new canonical coordinates are denoted with the same symbols as the old ones), the expansion 18 of the new Hamiltonian can be written as Once again, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., f (r,s) l = f (II; r,s) l ∀ l ≥ 0 , s ≥ 0 . Hence, each term generated by the Lie derivatives with respect to χ (r) 2 is added to the corresponding class of functions. This is made by the following sequence 19 of redefinitions: In fact, since χ (r) 2 ∈ P 2,rK is either quadratic in √ I or linear in p, each application of the Lie derivative operator L χ (r) 2 does not modify the degree in the square root of the actions, while the trigonometric degree in the angles is increased by rK . By applying Lemma 4.1 one can verify also that f (r,s) l ∈ P l,sK , ∀l ≥ 0, s ≥ 0 .
Because of the homological equation (30), it immediately follows that the term that cannot be removed, that is f (r,r) 2 = f (II; r,r) 2 q,α ∈ P 2,0 , is exactly of the same type with respect to the main term that is linear in the actions, i.e., ω (r−1) · p + Ω (r−1) · I. It looks then natural to update the angular velocity vectors so that where, as usual, the symbols ∇ p and ∇ I denote the gradient with respect to the action variables p and I, respectively, and to set f (r,r) 2 = 0 . Therefore, the expansion of the Hamiltonian H (r) can be 18 In the third row of (33), it is possible to start the second sum from r + 1 instead of r , being f where f (r,s) l ∈ P l,sK and E (r) ∈ P 0,0 is a constant. It is now possible to iterate the algorithm, by performing the (next) (r + 1)-th normalization step. The convergence of this normal form algorithm is proved in [2] under suitable conditions.
In order to implement such a kind of normalization algorithm with the aid of a computer, we have to deal with Hamiltonians including a finite number of summands in their expansions in Taylor-Fourier series. To fix the ideas, let us suppose that we set a truncation rule in such a way as to neglect all the terms with a trigonometric degree greater than N S K, for a fixed positive integer value of the parameter N S . After iteratively performing N S steps of the constructive algorithm, we end up with an approximation of the Hamiltonian which is in the normal form corresponding to an elliptic torus, i.e., The Hamiltonian H (N S ) represents the natural starting point for the application of a second (Birkhofflike) algorithm, which aims to produce a new normal form in such a way to remove the dependence on the angles q, as explained in the next Subsection.

Construction of the resonant normal form in such a way to average with respect to the angles q
where E B is a constant term, representing the energy, (p, q) ∈ R n 1 ×T n 1 , (I, α) ∈ R n 2 ≥0 ×T n 2 are actionangle variables, (ω B , Ω B ) ∈ R n 1 × R n 2 are the frequencies, N S is a fixed positive integer (ruling the truncations in the Fourier series) and g (0,s) l ∈ P l,sK , ∀ l ≥ 0, 0 ≤ s ≤ N S . In practice, we are starting from the normalized Hamiltonian of the previous Subsection H (N S ) , given by equation (37), where we have defined H ∈ P l,sK ∀ l ≥ 0, 0 ≤ s ≤ N S ; this is done also in order to simplify the notation. By comparison with equation (37), it is easy to remark that g (0,s) l The aim of the present algorithm is to delete the dependence of H (0) B on the angles q , reducing by n 1 the number of degrees of freedom. In order to do this, we have to act on the terms g (0,s) l (p, q, I, α) such that s ≥ 1 and l ≥ 3, removing their dependence on q ; indeed, for s = 0 , the sum l≥3 g (0,0) l (p, I) does not depend on the angles, thus it is already in normal form. This elimination can be done by a 20 We use the symbol H sequence of canonical transformations. If convergent, this would lead the Hamiltonian to the following final normal form: where (p,Ĩ,α) denote the new variables; it is evident that, having removed the dependence onq , the conjugate momenta vectorp is constant. However, as typical of the computational proceduresà la Birkhoff, the constructive algorithm produces divergent series if the normalization is iterated infinitely many times. For this reason, it is convenient to look for an optimal order of normalization to which the algorithm is stopped. In our approach, we have not to consider such a problem, because we are dealing with truncated series; this is done in order to keep our discussion as close as possible to the practical implementations where the maximal degree in actions of the expansions is usually rather low. The generic r-th step of this new normalization algorithm is defined as follows. After r − 1 steps, the Hamiltonian (38) takes the form with g (r−1,s) l ∈ P l,sK . The r-th normalization step consists of a sequence of N S substeps, each of them involving a canonical transformation which is expressed in terms of the Lie series having χ (j; r) B as generating function, with j = 1, . . . , N S . Therefore, the new Hamiltonian introduced at the end of the r-th normalization step of this algorithm is defined as follows: (41) The generating functions χ (j; r) B are defined so as to remove the dependence on q from the perturbing term of lowest order in the square root of the actions, i.e., N S s=1 g (r−1,s) r+2 (p, q, I, α) .

j-th substep of the r-th step of the algorithm constructing the resonant normal form
After j − 1 substeps, the Hamiltonian can be written as follows: where, for j = 1 , we set H In a similar way to what has been done previously, it is convenient to first define the new Hamiltonian terms as the old ones, i.e., g (j; r,s) l = g (j−1; r,s) l ∀ l ≥ 0 , 0 ≤ s ≤ N S ; hence, each term generated by the Lie derivatives with respect to χ (j; r) B is added to the corresponding class of functions. This is made by the following sequence 21 of redefinitions g (j; r,s+mj) l+mr In fact, since χ (j; r) B ∈ P r+2,jK , each application of the Lie derivative operator L χ (j; r) B increases the degree in square root of the actions and the trigonometrical degree in the angles by r and jK , respectively. Moreover, thanks to the homological equation (43) and the second rule included in formula (47) (in the case with m = 1), one can easily remark that g (j; r,j) r+2 = g (j−1; r,j) r+2 q . By applying Lemma 4.1 one can verify also that g (j; r,s) l ∈ P l,sK , ∀l ≥ 0, s ≥ 0 . The r-th step of the algorithm constructing the resonant normal form is completed at the end of the iterative repetition of the j-th substep for j = 1, . . . , N S . Therefore, the expansion of the Hamiltonian can be written in the following form: where g (r,s) l := g (N S ; r,s) l , ∀ l ≥ 0 , 0 ≤ s ≤ N S . Then, the normalization algorithm can be iteratively repeated. Since we are interested in the computer implementation, we consider finite sequences of Hamiltonians whose expansion is truncated up to a finite degree, say, N L in the square root of the actions. Therefore, the iteration of N L − 2 normalization steps of the algorithm constructing the resonant normal form are sufficient to obtain The Hamiltonian (49) does not depend on the angles q. Therefore, the corresponding actions p are constant and they can be considered as parameters whose values are fixed by the initial conditions; this allows us to decrease the number of degrees of freedom by n 1 , passing from n 1 + n 2 to n 2 .
5 Application of the normalization algorithms to the secular quasiperiodic restricted model of the dynamics of υ-And b The SQPR model can be reformulated in such a way as to resume the form of a Hamiltonian of the type (17), to which we can sequentially apply both normalization procedures described in the two previous Subsections. In fact, the canonical change of variables allows to rewrite the expansion of the SQPR Hamiltonian (13) as follows: H sec, 2+3/2 (p, q, I, α) = ω 3 p 3 + ω 4 p 4 + ω 5 p 5 c l,k ( I 1 ) l 1 ( I 2 ) l 2 e i(k 1 α 1 +k 2 α 2 +k 3 q 3 +k 4 q 4 +k 5 q 5 ) , (51) where k = (k 1 , . . . , k 5 ) ∈ Z 5 . The r.h.s. of the above equation can be expressed in the general and more compact form described in equation (17), by setting n 1 = 3 , n 2 = 2 , ω (0) = ω = (ω 3 , ω 4 , ω 5 ) ∈ R 3 , that are the fundamental frequencies of the two outer planets (described in equation (8)), while Ω (0) ∈ R 2 can be easily determined by performing the so called diagonalization of the Hamiltonian part that is quadratic in the square root of the actions I and not depending on the angles q (see, e.g., [8]). In the equation above, the parameters N L and N S define the truncation order of the expansions in Taylor and Fourier series, respectively, in such a way to represent on the computer just a finite number of terms that are not too many to handle with; in our computations we fix N L = 6 as maximal power degree in square root of the actions and we include Fourier terms up to a maximal trigonometric degree of 8, putting N S = 4 , K = 2 . We recall that setting K = 2 is quite natural for Hamiltonian systems close to stable equilibria as it is for models describing the secular planetary dynamics, see, e.g., [10]. Let us also remark that a simple reordering of the summands according to the total trigonometric degree |k| in the angles (q, α) allows us to represent the second row of formula (51) as a sum of Hamiltonian terms each of them is belonging to a functions class of type P l,sK , which is unique for any positive integer K if we ask for the minimality of the index s. These comments can be used all together in order to formally verify that the new expansion of H sec, 2+3/2 in (51) can be finally reexpressed in the same form as H (0) in (17).
Furthermore, in the case of our SQPR model of the secular dynamics of υ-And b, the only term depending on the action variables p (that are the so called dummy variables) is ω (0) · p ; thus, none of the Hamiltonian term f (0,s) l depends on p . This fact would allow to introduce some simplification in the computational algorithm. For instance, the value of the angular velocity vector ω (0) is not modified during the first normalization procedure (i.e. the algorithmic construction of the elliptic tori) and it remains equal to its initial value, given by the fundamental frequencies described in (8). Therefore, the expansion of the starting Hamiltonian in the special case of our SQPR model can be rewritten as however, in our opinion, for what concerns the general description of the previous Subsections it has been worth to consider also an eventual dependence of f (0,s) l on p in order to keep the discussion of the constructive procedure as general as possible.
The first algorithm to be applied aims to construct the normal form corresponding to an invariant elliptic torus. It starts from the Hamiltonian H sec, 2+3/2 rewritten in the same form as H (0) in (17) (more precisely as in (52)) and its computational procedure is fully detailed in Subsection 4.1. Therefore, we perform N S normalization steps of this first normalization algorithm. This allows us to bring the Hamiltonian in the following (intermediate) normal form: where f (N S ,s) l ∈ P l,sK ∀ l = 3, . . . , N L , s = 0 , . . . , N S and the angular velocity vector related to the angles q is such that ω (N S ) = ω (0) = (ω 3 , ω 4 , ω 5 ), whose components are given in (8).
It is now possible to apply the second algorithm aiming to construct a resonant normal form where the dependence on the angles q = (q 3 , q 4 , q 5 ) is completely removed. Such a normalization starts from the Hamiltonian H (N S ) obtained after the first normalization procedure. Therefore, we perform N L − 2 normalization steps of the above algorithm, each of them involving N S substeps as described in Subsection 4.2; this allows us to bring the Hamiltonian in the following (final) normal form: where g (N L −2,s) l ∈ P l,sK ∀ l = 3, . . . , N L , s = 0 , . . . , N S and it still holds true that ω B = ω (0) . All the algebraic manipulations that are prescribed by the normal form algorithms have been performed by using the symbolic manipulator Mathematica as a programming framework.
We emphasize that H 2DOF is an integrable Hamiltonian. In fact, due to the preservation of the total angular momentum, discussed in Remark 3.1, the following invariance law 22 is satisfied: thus, from the Hamilton's equations for H 2DOF , we can immediately deduce that I 1 + I 2 is a constant of motion. Therefore, H 2DOF is integrable because of the Liouville theorem (see, e.g., [7]), since it admits a complete system of constants of motion in involution, that are the dummy variables p (which could be disregarded in (53), reducing the model to 2 DOF), I 1 + I 2 and the Hamiltonian itself. In view of the numerical explorations of the dynamical evolution of our new model described by the integrable Hamiltonian H 2DOF (p, I, α), it is convenient to introduce the canonical transformations related to the so called semianalytic method of integration for the equations of motion (see, e.g., [10]). In order to fix the ideas, let us focus on the second algorithm, designed to construct a resonant normal form. This normalization procedure can be summarized by the transformation that is obtained by 22 Equation (54) can be easily checked, by explicitly performing the derivatives on the expansions (53) which are computed using Mathematica. However, it is worth to sketch also a more conceptual justification. Indeed, it would not be difficult to verify that all the Lie series introduced in Subsections 4.1-4.2 preserve the invariance law described in Remark 3.1. By comparing the definitions of the canonical coordinates in (50) and (2), one can immediately realize that the angles −α1 and −α2 are nothing but the longitudes of the pericenter and of the node, respectively, of υ-And b. Therefore, taking into account that H2DOF does not depend on the angles q, the invariance law discussed in Remark 3.1 can be rewritten in the form (54).
iteratively applying all the Lie series to the canonical variables. This is done as follows: for i = 1, 2 . We introduce the symbol C B to denote the change of coordinates 23 defined by the above expressions, i.e., (I, α) = C B (q,Ĩ,α) . We can proceed in the same way for what concerns the algorithm constructing the normal form corresponding to an invariant elliptic torus. In fact, we first introduce the application of all the Lie series to the canonical variables in such a way to write, ∀ i = 1, 2 , finally, we use the symbol C to summarize the whole change of coordinates that is defined by the whole expression above, i.e., (I, α) = C (q,Î,α) . Let us now introduce the function F : where we have omitted to put the˜symbol on top of q in order to stress that the angles q are not affected by the change of coordinates. Moreover, let also introduce the symbol A to denote the usual canonical transformation defining the action-angle variables for the harmonic oscillator, i.e., by formula (50), in our case this means that By applying the Exchange Theorem (see [11] and [6]), the solutions of the equations of motions related to H 2DOF can be mapped to those for H sec, 2+3/2 . Indeed, assume that t → p(t),q(t),Ĩ(t),α(t) 23 Since none of the generating functions χ (j; r) B depends on p , the way that these dummy variables are modified by the application of the Lie series does not really matter, because they do not enter in Hamilton's equations of motion (15), under the Hamiltonian H sec, 2+3/2 . Since, however, the generating functions do depend on q (but not on their conjugate actions p, as we have remarked just above) in the arguments of CB we have included also the angles q that are not affected by any modification due to the application of the Lie series.
is an orbit corresponding to the integrable flow induced by H 2DOF ; in particular, in our model we have thatq(t) = ω B t = ωt, where the components of the angular velocity vector ω are given in equation (8). Therefore, the orbit is an approximate 24 solution of the Hamilton's equations (15). For our purposes, it is also useful to construct the inverse of the function F , which maps from the original canonical coordinates to the ones referring to the resonant normal form. Therefore, it is convenient to replace all the compositions of Lie series appearing in the r.h.s. of (56) with the following expressions, ∀ i = 1, 2 : gathering all the corresponding changes of coordinates allows us to define 25 C −1 (q, I, α). Proceeding in an analogous way, we can introduce the inverse function of C B ; in more detail, we can start from formula (55), by reversing the order of all the Lie series and by changing the sign to all the generating functions, then we can define C −1 B (q, I, α). Therefore, we can introduce also We are now ready to exploit the (cheap) numerical solutions of the 2 DOF integrable Hamiltonian, which is described in (53), in order to retrieve information about the secular dynamics of υ-And b through our SQPR model. This can be done thanks to the knowledge of the approximate solution (59). The initial conditions are selected in the same way as in Subsection 3.1.2: we consider the initial 24 There are at least two substantial reasons for which this motion law, produced by a (so called) semi-analytic integration scheme, is not an exact solution of the equations (15). Let us recall that Lie series define near-to-the-identity canonical transformations that are well defined on suitable restrictions of the phase space. However, we are always working with finite truncated series; therefore, the corresponding changes of variables cannot preserve exactly the solutions because infinite tails of summands are neglected. Moreover, in the resonant normal form H2DOF described in (53) we do not include the remainder terms; let us recall that they become dominant if the Birkhoff algorithm is iterated infinitely many times, making the series expansion of the normal form to be divergent. Therefore, the semi-analytic solutions are prevented to be exact also because of this second source of truncations acting on the series expansion of the Hamiltonians (instead of the Lie series defining the canonical transformations). As a final remark, let us also recall that, in order to be canonical, the change of coordinates should include also the dummy actions p, in which we are not intested at all because they do not exert any role in the equations of motion (15). 25 Of course, since also C −1 : , C and C −1 share the same domains and codomains, which are different between them), then C −1 cannot be considered as the inverse function in a strict sense. However, if we extend trivially both these functions, in such a way to introduceĈ (q, I, α) = q, C (q, I, α) andĈ −1 (q, I, α) = q, C −1 (q, I, α) , thenĈ −1 would really be the inverse function ofĈ (where elementary properties of the Lie series described in Chap. 4 of [6] are also used and the small effects due to the truncations are neglected). Therefore, it is by a harmless abuse of notation that we are adopting the symbol C −1 . The same abuse will be made for what concerns the symbols C −1 B and F −1 .
orbital elements reported in Table 7 and the minimal possible value of the mass of υ-And b , i.e., m 1 = 0.674 M J . These data are completed with the values of (i 1 (0) , Ω 1 (0)) ranging in the 20×60 regular grid that covers I i ×I Ω = [6.865 • , 34 • ]×[0 • , 360 • ]; moreover, all these initial values of the orbital elements are translated in the Laplace reference frame, which refers only to the two outermost exoplanets. Hence, we can compute a set of 21×60 initial conditions of type I(0), α(0) = A −1 ξ 1 (0), η 1 (0), P 1 (0), Q 1 (0) , by using formula (2) with j = 1 , (50) and the definition (58). Finally, we can translate the initial conditions to initial values of the canonical coordinates found after the resonant normal form, by computing Ĩ (0),α(0) = F −1 0, I(0), α(0) . As shown below, an important information is obtained by a criterion allowing to identify those domains of initial conditions in which the series are either divergent or slowly converging. We introduce such a criterion to preselect initial conditions that are admissible. From a mathematical point of view, the identity (I(0), α(0)) = ( holds in a domain where the normalization procedure is convergent, provided that no truncations are applied to the series F and F −1 and that the computation of the series is not affected by any round-off errors. Due to the errors and truncations introduced in the computation, however, in general we obtain that (I(0), α(0)) = (I O (0), α O (0)) . In the domain where the series expansions are rapidly converging the difference (I(0), α(0)) − (I O (0), α O (0)) is small. When, instead, we obtain a large difference, this is an indicator that we are outside the domain of convergence of the series. The situation is represented graphically below.
(I(0), α(0)) In view of the above, we define the following preselection criterion of admissible initial conditions. For any initial condition (I(0), α(0)) we compute the quantities The use of the quantities √ I 1 and √ I 2 is motivated by the fact that they are of the same order of magnitude as the eccentricity and the inclination of υ-And b , respectively. Moreover, it is useful to define also the following ratios where t → ωt,ξ 1 (t),η 1 (t),P 1 (t),Q 1 (t) := ωt , A F (ωt,Ĩ(t),α(t)) is the approximate solution of Hamilton's equations (15), as produced by the semi-analytic integration scheme summarized in formula (59). Comparing formula (63) with the definition of the Poincaré canonical variables in (2), it is easy to realize that R 1 and R 2 are functions of the time that describe the behavior of the orbital excursions with respect to the eccentricity and the inclination of υ-And b , respectively. We then investigate the behavior of the following function: Note thatẽ 1 would be equal to the eccentricity of υ-And b ifĨ 1 was replaced by ξ 2 1 + η 2 1 /2, with (ξ 1 , η 1 ) defined in (2). However, the new actionĨ 1 is conjugated to ξ 2 1 +η 2 1 /2 wich is only nearly equal to ξ 2 1 + η 2 1 /2 , since the composition of the transformations C and C B is near-to-identity. Therefore, we can considerẽ 1 as an approximate evaluation of the eccentricity under the resonant normal form model.
For each pair i 1 (0), Ω 1 (0) of the 21 × 60 points definining the grid which covers I i × I Ω = [6.865 • , 34 • ] × [0 • , 360 • ] we determine the corresponding initial conditions of type (I(0), α(0)), as explained above, and we proceed as follows: • if max{r 1 , r 2 } > 1, then the corresponding initial condition is considered as "non-admissible", i.e. outside the domain of applicability of the series. Then, we skip the step below and pass directly to consider the next initial conditions of the grid; • If the initial conditions is admissible, we numerically 26 solve the equations of motion for the integrable Hamiltonian model with 2 DOF described in (53), using a RK4 method and starting from 0,Ĩ(0),α(0) ; during such a numerical integration, we compute the maximal values attained by the three previously defined quantitative indicators, that are The results about the maxima of the functions defined in (63)-(64) are reported in Figures 8-9. The white central regions of those pictures correspond to those pairs i 1 (0), Ω 1 (0) for which we obtain failure of the preliminary test, i.e. max{r 1 , r 2 } > 1. We immediately recognize that the missing part of the plots (where the determination of the initial conditions is considered so unreliable that the corresponding numerical integrations are not performed at all) nearly coincides with the central region of Figure 2a, where the orbital eccentricity of υ-And b reaches critical values. We conclude that the stability domain in the space of the initial values of i 1 (0) and Ω 1 (0) (which are unknown observational data) can be reconstructed in a reliable way through the application of the above criterion, which only involves the series transformations, as well as through the numerical solutions of our integrable secular model with 2 DOF. We emphasize that this allows to reduce significantly the computational cost with respect to the long-term symplectic integrations of the complete 4-body problem, which is a 9 DOF Hamiltonian system. Comparing the regions of the stability domain at the border near the (white) central ones, we see that all three numerical indicators plotted in Figures 8-9 increase their values when the unstable zone  Since the dependence of the Hamiltonian on the angles q (which describes the dynamics of the outer exoplanets) is removed from the 2DOF model, it seems reasonable that some of the resonances are not present in the normal form generated by the algorithmà la Birkhoff, although they play a remarkable role in the dynamics of more complex models.
6 Secular orbital evolution of υ-And b taking also into account relequations for the orbital motion of the innermost planet can be written as (68)

Numerical integration of the SQPR-GR model
Similarly as in Subsection 3.1.2, we numerically integrate the equations of motion for the secular quasi-periodic restricted Hamiltonian with general relativistic corrections, defined in formula (68).
As initial values of the orbital parameters a 1 (0), e 1 (0), M 1 (0) and ω 1 (0) we take those reported in Table 7; moreover, we set m 1 = 0.674 as value for the mass of υ-And b and (i 1 (0) , Ω 1 (0)) ranging in the 20 × 60 regular grid that covers Hence, it is possible to compute the corresponding initial values of the orbital elements in the Laplace reference frame (which is determined taking into account υ-And c and υ-And d only) and to perform 21 × 60 numerical integrations starting from all these initial data. Once again, for each numerical integration, we are interested in determining the maximal values reached by the eccentricity of υ-And b and by the maximal mutual inclination between υ-And b and υ-And c . The results are reported in the color-grid plots of Figure 10. By comparing Figure 10a with Figure 5a, one can immediately realize that the regions colored in blue are much wider in the former than in the latter one. Indeed, the darker regions refer to initial conditions which generate motions with maximal values of the eccentricity of υ-And b that are relatively low, while the zones colored in red or yellow correspond to such large values of the eccentricity implying that those orbits have to be considered unstable. Therefore, our numerical explorations highlight that the effects due to general relativity play a stabilizing role on the orbital dynamics of the innermost planet. This conclusion is in agreement with was already remarked about the past evolution of our Solar System, in particular for what concerns the orbital eccentricity of Mercury (see [15]).
Moreover, as already done in Section 3.1.2, in order to further explore the stable and chaotic regions of Figure 10a, we apply the Frequency Map Analysis method to the signal ξ 1 (t) + iη 1 (t) as produced by the numerical integration of the system (68), i.e., in the SQPR-GR approximation. We perform the numerical integrations as described at the beginning of the present Section, taking into account only a few values in I i for the initial inclinations, i.e. i 1 (0) = 6.865 • , 8.22175 • , 9.5785 • , 10.9353 • and Ω 1 (0) ∈ I Ω . In Figure 11 we report the behaviour of the angular velocity ν corresponding to the first component of the approximation of ξ 1 (t) + iη 1 (t), as obtained by applying the FA computational algorithm; we recall that this quantity is related to the precession rate of 1 . As initial value for the inclination i 1 (0) we fix 6.865 • for Figure 11a and 10.9353 • for Figure 11b. Also here, we do not report the cases (i 1 (0), Ω 1 (0)) ∈ {8.22175 • , 9.5785 • } × I Ω , since the behaviour of these plots is similar to the ones in Figure 11.  Figure 11: Behaviour of the fundamental angular velocity ν as obtained by applying the Frequency Map Analysis method to the signal ξ 1 (t) + iη 1 (t) , computed through the RK4 numerical integration of the SQPR-GR model (68), covering a timespan of 1.31072 · 10 5 yr. We take, as initial conditions, (i 1 (0), Ω 1 (0)) ∈ {6.865 • } × I Ω for the left panel and (i 1 (0), Ω 1 (0)) ∈ {10.9353 • } × I Ω for the right one.
The situation is well described by Figure 11a and analogous considerations hold for Figure 11b apart a few main differences which will be highlighted in the following discussion. When the values of Ω 1 (0) are ranging in [0, ∼ 120 • ] and [∼ 260 • , 360 • ] we can observe a regular behaviour of the angular velocity ν which is also nearly monotone, with the only exception around a local minimum. According to the Frequency Map Analysis method, such a regular regime is due to the presence of many invariant tori which fill the stability region located at the two lateral sides of the plot 10a. In the case of Figure 11a, this also applies when Ω 1 (0) is ranging in [∼ 150 • , ∼ 220 • ], which corresponds to the stable blue internal area of Figure 10a. On the other hand, in the case of Figure 11b, for the same range of initial values of the node longitude of υ-And b, the behaviour is not so regular; this is in agreement with the fact that in correspondence with i 1 (0) ∼ 11 • the plot of the maximal values of e 1 in the central region highlights the occurrence of chaotical phenomena. Moreover, for what concerns values of Ω 1 (0) in [∼ 120 • , ∼ 150 • ] and [∼ 220 • , ∼ 260 • ] (corresponding to the green stripes of Figure 10a), Figure 11a shows a behaviour typical of the crossing of a resonance in the chaotic region surrounding a separatrix. The value of the angular velocity for which this phenomenon takes place is, again, related to ω 4 −1.04 · 10 −3 (as it can be easily appreciated looking to the small plateau appearing in Figure 11b).
Comparing Figure 11 with Figure 6 the enlargement of the stable region is evident. Moreover, we can also see how much this phenomenon is influenced by the modification of the pericenter precession rate of the inner planet due to relativistic effects. Indeed, looking at the values reported on the y-axis of Figures 11 and 6, one can appreciate that the fundamental angular velocity, in the case of the SQPR model, takes values remarkably closer to zero with respect to those assumed in the case of the SQPR-GR model.
where E B;GR ∈ R, Ω B;GR ∈ R 2 and h (N L −2,s) l ∈ P l,2s ∀ l = 3, . . . , N L , s = 0 , . . . , N S . We emphasize that also H (GR) 2DOF is integrable because of the same reasons already discussed in Section 5; indeed, after having checked that ∂H (GR) 2DOF /∂α 1 + ∂H (GR) 2DOF /∂α 2 = 0 , we can apply the Liouville theorem, because there are two independent constants of motion, i.e., I 1 + I 2 and H (GR) 2DOF itself. Moreover, also for this new model we can reproduce the same kind of numerical exploration described in Section 5. In particular, we can compute the values of the numerical indicators R 1 MAX , R 2 MAX andẽ 1 MAX corresponding to each pair i 1 (0), Ω 1 (0) of the 21 × 60 points definining the regular grid which covers I i × I Ω = [6.865 • , 34 • ] × [0 • , 360 • ]. In the following, we analyze the color-grid plots for a few different values of the parameter ruling the truncation in the trigonometric degree, namely N S , and in the square root of the action, i.e., N L . The color-grid plots for the maximal value reached by R 1 andẽ 1 are reported in Figures 12-14.
Let us recall that R 2 MAX is an evaluation of the maximal excursion of the inclination of υ-And b. For the sake of brevity, its plots are omitted and they are not included in Figures 12-14, because in our numerical explorations the ranges of values experienced by R 2 MAX are so narrow that their analysis does not look so significant. Therefore, it is better to focus on the plots of R 1 MAX andẽ 1 MAX ; let us recall that both of them refer to the eccentricity of υ-And b. By comparing Figures 12-14, one can appreciate a well known phenomenon concerning the constructive algorithmsà la Birkhoff: the greater the number of normalization steps (i.e., N L − 2), the smaller the domain of applicability (see, e.g., [6] for the discussion about the determination of the optimal step).   Figure 10a, we observe that in the cases with N L = 4, 5 our computational procedure is able to reconstruct with a good accuracy the U -shaped border of the stability domain. Note that the horizontal strip at the bottom of the plots 28 corresponds to orbital motions which look stable, since the eccentricity of υ-And b does not reach large values (with the eventual exception of the narrow green areas that in Figure 10a are expected to correspond to resonant regions). This highlights a main difference with the non-relativistic model discussed in Section 5, because in that case there is an interval of values of Ω 1 (0) centered about 180 • for which none of the initial inclinations i 1 (0) ∈ [6.865 • , 34 • ] corresponds to a stable orbital configuration (see Figure 5a). The reliability of our simplified 2 DOF Hamiltonian model (which is defined in formula (69) and takes into account also GR corrections) is enforced also by the fact that it is able to capture also this phenomenon.