The generalized non-conservative model of a 1-planet system revisited

We study the long-term dynamics of a planetary system composed of a star and a planet. Both bodies are considered as extended, non-spherical, rotating objects. There are no assumptions made on the relative angles between the orbital angular momentum and the spin vectors of the bodies. Thus, we analyze full, spatial model of the planetary system. Both objects are assumed to be deformed due to their own rotations, as well as due to the mutual tidal interactions. The general relativity corrections are considered in terms of the post-Newtonian approximation. Besides the conservative contributions to the perturbing forces, there are also taken into account non-conservative effects, i.e., the dissipation of the mechanical energy. This dissipation is a result of the tidal perturbation on the velocity field in the internal zones with non-zero turbulent viscosity (convective zones). Our main goal is to derive the equations of the orbital motion as well as the equations governing time-evolution of the spin vectors (angular velocities). We derive the Lagrangian equations of the second kind for systems which do not conserve the mechanical energy. Next, the equations of motion are averaged out over all fast angles with respect to time-scales characteristic for conservative perturbations. The final equations of motion are then used to study the dynamics of the non-conservative model over time scales of the order of the age of the star. We analyze the final state of the system as a function of the initial conditions. Equilibria states of the averaged system are finally discussed.

effects. These perturbations cause the periapses rotation and the precession of the orbital nodes as well as spins of the rotating bodies. Considering the longest time scale which is comparable to the age of the parent star, one has to take into account also a dissipation of the mechanical energy. It is believed, that the most important physical mechanism dissipating the energy is the tidal perturbation of the velocity field in parts of bodies possessing non-zero turbulent viscosity (e.g., Zahn 1977). This takes place in the convective zones of stars and planets. The mechanisms of the energy dissipation, particularly active in stars and Jupiterlike planets, were studied by many authors (e.g., Goldreich and Soter 1966;Ogilvie and Lin 2004;Wu 2005a,b;Ogilvie and Lin 2007;Miller et al. 2009;Gu and Ogilvie 2009;Arras and Socrates 2010). An open problem is to estimate values of physical parameters characterising the strength and time-scales of these processes, and in turn, the time-scale of the dissipative evolution of planetary systems.
It is well known, that the energy loss leads to a variation (a decrease in general) of both the semi-major axis and eccentricity, as well as to evolution of spins, their directions and magnitudes. The planetary dynamics of systems with energy loss were considered both in cases with one and two planets (e.g., Mardling and Lin 2002;Witte and Savonije 2002;Dobbs-Dixon et al. 2004;Mardling 2007;Barker and Ogilvie 2009;Leconte et al. 2010;Rodríguez and Ferraz-Mello 2010;Michtchenko and Rodríguez 2011;Rodríguez et al. 2011;Correia et al. 2011;Laskar et al. 2012). The equations of motion of such systems are usually derived through introducing a dissipative force. For instance, the well known model of Hut (1981), assumes that this force emerges due to a time delay in forming the tidal bulge and the orbital motion of a perturber. That implies a non-zero angle between the radius vector of the deformable body with respect to its companion, and the axis of symmetry of the tidally deformed object. In more elegant way, the tidal force was derived from the energy loss functionĖ defined by Eggleton et al. (1998).
In this paper, we found a more straightforward derivation of the dissipative model that relies on the Lagrangian equations of the second kind. As we will show, this approach makes it possible to obtain quite simply the dissipative forces acting in the N -body system; however, we limit here the derivation to N = 2. Moreover, to derive the equations governing the evolution of angular velocities, both in conservative, as well as in non-conservative models, we should not apply the Euler equations, which hold only for a specific form of the potential energy V that has to be then a function of Euler angles φ, θ, ψ, and should not depend on their time derivativesφ,θ,ψ. This is only true in the case of the rigid body. For deformable objects, V is a function of the angular velocity Ω Ω Ω = Ω Ω Ω (φ, θ, ψ,φ,θ,ψ) and thus it does not fulfill these assumptions. Hence, the Euler equations stating that the time derivative of the rotational angular momentum equals to the torque acting on the rigid body do not hold in general. However, we will show here that it is still possible to obtain the equations of the evolution of the angular velocities in vectorial form, which is reminiscent of the classic Eulerian equations.
In Sect. 3, we find expressions for these functions in a particular model considered in this paper. We use the polytropic model of Chandrasekhar (1933a,b) to calculate the internal structure and a deformation of figures of both objects. Relativistic correction to the Lagrangian is taken from Brumberg (2007). To determine the energy loss function, we use a simple model by Eggleton et al. (1998).
In Sect. 4, the derived equations of motion are averaged out over all angles that vary in time-scales related to the conservative evolution. These angles, ordered from the fastest to the slowest one are the following: the mean anomaly, the horizontal angle of the precessing planetary spin in the orbital frame, and the argument of pericenter. As the result we obtain the equations of motion describing the dissipative dynamics only, which are then studied in Sect. 5. We analyze the final state of the planetary system in terms of the initial conditions. As a particular solution, we discuss the equilibrium permitted in the system.

The equations of motion
We shall consider a planetary system in terms of the mechanical system with the Lagrangian L = L(q 1 , . . . , q n ,q 1 , . . . ,q n ) and energy loss functionĖ =Ė(q 1 , . . . , q n ,q 1 , . . . ,q n ), where q i are the generalized coordinates, andq i are the generalized velocities. Index i spans from 1 to n, where n is the number of the degrees of freedom. The dynamical evolution of the system is governed by solutions to the Lagrangian equations of the second kind (see e.g., Greiner 2003, p. 328): The above equations are correct whenĖ fulfills the following condition (e.g., Goldstein et al. 2002, p. 63) Particularly, it takes place whenĖ is a quadratic form of the generalized velocities. The explicit form ofĖ is given further in the text. Although it is not a quadratic form inq, it fulfills the above condition as well (it is shown in Appendix A).
In the problem considered here, the generalized coordinates are the following: , y, z, s x , s y , s z , p x , p y , p z }. The first three coordinates are components of the position vector of the planet relative to the star, r = [x, y, z] T . The last six coordinates are vectors of three independent components in two unit quaternions, s = [s x , s y , s z ] T , p = [p x , p y , p z ] T : s = s 0 + i s 1 + j s 2 + k s 3 , s k ∈ R, The components of the quaternions as well as their time derivatives are related as follows: s 2 0 + s 2 1 + s 2 2 + s 2 3 = 1, s 0ṡ0 + s 1ṡ1 + s 2ṡ2 + s 3ṡ3 = 0 and similarly, for p. Therefore, only three components of each quantity are independent. We choose these independent components as follows: s x = s 1 , s y = s 2 , s z = s 3 and p x = p 1 , p y = p 2 , p z = p 3 . Thus s 0 ,ṡ 0 , p 0 ,ṗ 0 are functions of s,ṡ, p,ṗ, according to Eq. (3). The Lagrangian of the system, as well asĖ, is assumed to be function of planetary position, velocity and spin vectors of both bodies, i.e., L = L(r,ṙ, Ω Ω Ω 0 , Ω Ω Ω 1 ),Ė =Ė(r,ṙ, Ω Ω Ω 0 , Ω Ω Ω 1 ). It is well known, that the angular velocity may be expressed with the help of quaternions in the following form (see e.g., Heard 2006), p. 49: and similarly for Ω Ω Ω 1 with quaternion p. Thus Ω Ω Ω 0 = Ω Ω Ω 0 (s,ṡ) and Ω Ω Ω 1 = Ω Ω Ω 1 (p,ṗ). There exists also an inverse relation: The angular velocities are then Ω Ω Ω 0 = Ω Ω Ω 0 (s,ṡ) and Ω Ω Ω 1 = Ω Ω Ω 1 (p,ṗ). The full set of Lagrange equations is then the following: It may be shown (see Appendix B), that the second equation in the above set may be transformed into the matrix equation This equation may be solved with respect to X, leading to the following solution: where the star angular velocity vector has components Ω Ω Ω 0 = [Ω 0,x , Ω 0,y , Ω 0,z ] T . Similar expression may be obtained for the angular velocity of the planet, The product in the right-hand side of the equation is simply the vector product Ω Ω Ω 0 × Y. The final equations of the evolution of Ω Ω Ω l have the following form: The above equations are valid for systems, the Lagrangian of which depends on the space orientation of extended objects only through the angular velocities. As we will show, for the case considered here, these equations have similar explicit form as the Euler equations. Nevertheless, they were derived under assumption of a particular form of the potential function V = V (r, Ω Ω Ω 0 , Ω Ω Ω 1 ) suitable to our model.

The Lagrangian and the dissipative function
The Lagrangian of the system is given by the following expression: where are the kinetic and potential energies of two point masses interacting in accord with the Newtonian gravity. The masses of the star and the planet are denoted with m 0 and m 1 , respectively, β ≡ (1/m 0 + 1/m 1 ) −1 is the reduced mass, r ≡ ||r|| and k is the Gauss constant. Terms V rot and V tid are for the perturbing potential energy of two extended, non-spherical objects. We assume, that each object is deformed due to its own rotation, as well as to the mutual tidal interaction with the other body. These two effects are considered separately. Such a simplification is correct to the first order. Moreover, we assume that the planet deforms the shape of the star due to the point mass interaction, and the star is a point-mass perturber deforming the figure of the planet. The rotational kinetic energy of a deformable object is given by T rot , while the relativistic term of the Lagrangian is denoted by L rel .

Potential of the axially symmetric object
Both perturbing terms V rot and V tid have an axial symmetry. The axis of symmetry of V rot coincides with a direction of the angular velocity vector of a particular object, while the axis of symmetry of V tid coincides with a direction of a vector r joining the mass centers of the bodies. It is well known, that the potential energy of a system that consists of a point mass m and extended axially symmetric object of mass M and characteristic radius R may be developed in Taylor series with respect to (R/r ): where J 2l are the Stokes coefficients, and P 2l (x) are the Legendre polynomials [see, e.g., Schaub andJunkins (2003), p. 480 or Murray andDermott (2000), p. 133]. Here, these series are truncated at terms proportional to J 2 . In the case of the rotational deformation, For the tidal perturbations induced by a body at r,ẑ =r ≡ r/r . The zonal Stokes coefficients J 2 are expressed by the following integrals (e.g., Schaub and Junkins 2003, p. 477): where the density ρ (R, ϑ) depends on the magnitude of perturbations due to the rotation and tides. That function may be determined in a coordinate system which has its z-axis alonĝ z. Angle ϑ is measured from this axis, between 0 (the "north pole") to π (the "south pole"), and the radial variable is R ∈ [0, R (ϑ)]. The shape of a body is described in terms of R (ϑ), which may be approximated by R when the perturbation is small.

Calculating the Love numbers
To determine and describe tidal perturbations, we shall need to calculate the Love numbers 1 .
To accomplish this task, we apply remarkable results of Chandrasekhar (1933a), who considered uniformly rotating polytropes 2 and obtained equations for the density function ρ (R, ϑ). He postulated the following form of ρ: where ρ c is the central density, n is the polytropic index and A is a function, which is usually introduced to define dimensionless variable z (see, e.g., Kippenhahn and Weigert 1994, p. 176 for the explicit formulae). He found that for slowly rotating body: Function w(z) may be obtained by solving the Lane-Emden equation (Eq. [12] in the cited paper). Functions U 0 (z) and U 2 (z) may be found by solving equations [37 1 ] and [37 2 ] in the cited paper, respectively (let us note that the notation used here is different from the original one). Coefficients a n depend on the polytropic index n and may be derived as solutions to the problem considered by Chandrasekhar. Using the resulting ρ (R, ϑ), it is relatively easy to show that the Stokes coefficient J 2 defined by Eq. (15) has the following form k L,r = k L,r (n) = 6 a n 5 z 5 where z n is dimensionless size of undisturbed body, i.e., w (z n ) = 0 and k L,r may be attributed to the fluid Love number of the object deformed due to the rotation (Munk and MacDonald 1975, p. 26) which is a function of index n. In the range n ∈ [0, 4], this function is very well approximated by the formulae: log 10 k L,r ≈ f (n) = log 10 3 2 + α 1 n + α 2 n 2 + α 3 n 3 , In the next paper, Chandrasekhar (1933b) considered the deformation of the polytropic body having mass M by the tidal force emerged due to the outer point-mass perturber of mass m. He found the density function, which may be then used to calculate the coefficient J 2 of 1 In stellar astrophysics, the so called apsidal motion constant (instead of the Love number, which is twice as much) is used as a physical parameter giving the rate of the rotation of apsidal line of the binary orbit in space. On the other hand, in planetary astrophysics, the Love number is usual (not only when a planet-moon system is considered, but also for a star-planet system). In this work we use this planetary convention. 2 A polytrope is understood here as a gaseous object being in hydrostatic equilibrium under its own gravity.
The pressure-density relation is given by the formulae P = Kρ (n+1)/n (where pressure is denoted by P, density by ρ, K is a constant factor and n is a polytropic index). A rotating polytrope or/and being attracted by some external force is then called a perturbed polytrope.
tidally distorted object. It may be shown that for small perturbations, the linear approximation is valid, and then: where k L,t is the tidal Love number (Munk and MacDonald 1975, p. 27) and is given by the formulae: Thus the fluid Love numbers of rotationally and tidally deformed object are numerically equal. Let's denote them by k L . It is not surprising, because k L is a physical measure of deformability of an object under the attraction of some force (here, centrifugal and tidal forces are considered). The Love number relates to the quantity Q introduced in Eggleton et al. (1998, Eq. 15c), i.e., k L = Q/(1 − Q). The numerical results obtained above are in agreement with the results stated in the cited work. However, the computation of k L (or the apsidal motion constant) were performed by many authors before (e.g., Brooker and Olle 1955), to make this work self-contained, we present a brief overview of one of the way which leads to numerical values of k L . We choose to make use of Chandrasekhar's results, nevertheless this is not the only approach possible (see, e.g., Eggleton et al. 1998). Equation 15c from this last paper gives a formulae for Q for a more general model of mass distribution. For a polytropic model they find a polynomial expression, Eq. 19. Our numerical results arisen from Eq. (20) are in agreement with the results in the literature (e.g., Brooker and Olle 1955;Eggleton et al. 1998). However, the model of a polytrope being deformed by centrifugal as well as tidal forces is expected to work well for stars and gaseous planets, it is not expected so, when considering rocky planets (see, e.g., Bursa 1984).
To conclude, we write down the following forms of V rot and V tid where k L,0 and k L,1 are the Love numbers of the star and the planet, respectively. Their numeric values are given by Eq. (20). It is believed, that a mass distribution in Sun-like stars is well approximated by a polytropic model of index n ∈ [3, 4]. In a case of Jupiter-like planets, the appropriate value is n ∈ [1, 2].

The kinetic energy of rotating, extended bodies
The kinetic energy of rotations of extended objects is given by a sum: where I 0 and I 1 are the moments of inertia of the star and the planet, respectively. The moment of inertia is defined in a coordinate system with the origin that coincides with the mass center of the body as the following integral (see e.g., Schaub and Junkins 2003, p. 129): wherer points towards the mass element in the body, and [r] is the so called tilde matrix of a vectorr. For a vector a = [a x , a y , a z ] T , it has the form of: The matrix multiplication of [a] and a vector b is equivalent to the the vector product a × b.
A multiplication in Eq. (26) gives: The density function of the l-th object (l = 0, 1) may be written as a sum: where ρ (0) l (r ) is for the density function of undisturbed body,r ≡ ||r||, ρ (rot) l (r) is a correction stemming from the rotational deformation. Thus the moment of inertia may be expressed as the following sum: Each term in that equation represents the moment of inertia for some density function. The first term is the moment of inertia of undistorted body possessing spherically symmetric distribution of its mass. Using a simple polytropic model of that object, we obtain: where E is a unit 3 × 3 matrix, i.e., E = diag(1, 1, 1). For the homogeneous mass distribution (n = 0), coefficient κ = 1. In the range of n ∈ [0, 4], the coefficient is well approximated by the formulae similar to that ones derived for the Love number: log 10 κ ≈ f (n) = α 1 n + α 2 n 2 + α 3 n 3 , (33) α 1 = −0.2061, α 2 = +0.0297, α 3 = −0.0140.
The resulting form of T rot is then:

The relativistic perturbing Lagrangian
The model considered here includes also the relativistic perturbations. We do not include terms containing Ω Ω Ω 0 , Ω Ω Ω 1 , because their contributions to the equations of motion are of the second order. The relativistic contribution to the Lagrangian of the two point masses is given by the formulae (e.g., Brumberg 2007): where the mass parameters have the following form:

The dissipative function
A function describing dissipation of the mechanical energy of the system is the following sum:Ė where the first term comes from the energy dissipation in the star due to tidal interaction with a planet, while the second term expresses the energy dissipation in the planet due to the interaction with the star. The particular term ofĖ is given by a simple model proposed by Eggleton et al. (1998, Eq. 43): for l = 0, 1. Parameters σ 0 , σ 1 are dissipation constants for the star and the planet, respectively 3 . They may be expressed in the following form : where λ 0 and λ 1 are non-dimensional coefficients of the energy loss rates in the star and the planet, respectively. A calculation of these values is complex, and we postpone the derivation to other work. Actually, that is basically beyond the scope of the present paper. Following the literature (e.g., Ogilvie and Lin 2004;Wu 2005a,b;Ogilvie and Lin 2007;Hansen 2010), we may assume these coefficients in the range of λ 0 ∈ [10 −7 , 10 −4 ] and λ 1 ∈ [10 −7 , 10 −4 ]. Particular values of the dissipative coefficients are not well known. In the paper by Hansen (2010), he calibrated the equilibrium tide model comparing the results of calculations with the observed relation of the orbital period and eccentricity of close binaries for which the age is known. He obtainedσ * ≈ 5.3 × 10 −5 , which may be "translated" to λ 0 , through relation λ 0 =σ * (1 + k L,0 ) −2 , which then gives λ 0 ≈ 5 × 10 −5 . The dissipation constants for Jupiter-like planets were derived by analysis of the period-eccentricity relations of known systems with these planets. The resultingσ p ≈ 6.8 × 10 −7 gives λ 1 ≈ 5 × 10 −7 .

The explicit form of the equations of motion
Furthermore, we consider a particular form of the Lagrangian and the dissipative function: The derivative of L overṙ does not depend on Ω Ω Ω l and similarly, the derivative of L over particular angular velocity does not depend on the other Ω Ω Ω, nor on the linear velocity. Hence, the vector equations for r, Ω Ω Ω 0 and Ω Ω Ω 1 are separable. Equations in the top row of (6) describe the orbital motion of the planet. The Eq. (11) with l = 0, 1 are for the evolution of the angular velocities of the star and the planet, respectively.

The equations of translational motion
The top-row vector equation of (6) has to be solved with respect tor to obtain the explicit form of the equations of translational motion. The dependence of L onṙ appears in the T p-p and L rel terms. Only the first one is a uniform quadratic function ofṙ. Nevertheless, we may solve the problem. The solution has the following form in the first order approximation: where the first term is the Keplerian acceleration, the second term has its origin in the perturbing potential of tidally deformed bodies. The term in the second row emerges due to the rotational deformation of both objects. The perturbing acceleration in the third row is for the relativistic contribution, while the last term comes from the dissipation of a mechanical energy. In spite of very different approaches used in Eggleton et al. (1998) and in this work, the final equations of motion, especially their dissipative part, are the same (see their Eq. 34 with f TF given by Eq. 45). According with the above notation

The equations of the evolution of angular velocities
The equations of the evolution of angular velocities of the star and the planet, denoted here by Ω Ω Ω 0 and Ω Ω Ω 1 are obtained through solving Eq. (11) with respect toΩ Ω Ω l (l = 0, 1). To the first order, one derives the following form of these equations: where the first-row term is the same as would be derived directly from the Euler equations. The second-row terms emerge due to dependency of the potential energy on Ω Ω Ω l . We note that the non-rigid-body contribution to the equations are the same order as the rigid-body term. The dissipative term stands in the last row of (50) and it is exactly the same as in Eggleton et al. (1998). To compare the results, see their Eq. 36 with f TF given by Eq. 45.

Effects not included in the model
The equations derived so far do not take into account a few physical and perturbing phenomena, which we would like to list here explicitly. First of all, we assumed that the star does not evolve in time keeping its internal structure constant. However, if we consider the evolutionary time scale, the stellar radius as well as the internal structure change: the time evolution of R 0 and the mass distribution alter the inertia moment, which then may provide an additional term in Eq. (50). Hence, also the Love number is not constant over time. Evolution of the internal structure, for instance, the size of the convective zone, implies that none of coefficients λ 0 remains constant. Also a significant stellar wind would be responsible for yet additional term in Eq. (50), (e.g., Dobbs-Dixon et al. 2004;Barker and Ogilvie 2009). The same discussion applies to the planet. Its radius and internal structure also may vary over time, which is mainly attributed to the tidal and thermal influence by the parent star (e.g., Gu and Ogilvie 2009;Miller et al. 2009;Hansen 2010;Arras and Socrates 2010;Ibgui et al. 2011). However, the tidal evolution of the interiors of hot-Jupiters is not well recognized. Moreover, very close-in planets revolving around their stars with periods of the order of one day, may be close enough to loose their mass through the Lagrangian point L 3 (Li et al. 2010), as they may fill up the Roche lobe.
In the considered model, we assume, that the material is purely viscous, which means that there is no unelasticity nor rigidity in the bodies. The model ignores the so-called -effect, which leads to differential rotation in the stars (see, e.g., Kitchatinov and Rudiger 1993;Kueker et al. 1993). While this effect is well studied for single stars, it is not understood yet how a close planetary companion can tidally affect the differential rotation of the star. On the other hand, the tidal influence of the star on the planet's differential rotation is not studied well enough. Here, we omit this effect. We also ignore the effect of magnetic rigidity (see, e.g., Williams 2006).
As it was noted already, the dissipation parameters σ l (and also λ l ) depend on the tidal frequency, ω t . After Eggleton et al. (1998), we assume, that the energy is dissipated due to turbulent convection in the objects (thus, we consider the so-called equilibrium tide model). Therefore, the parameter σ l depends on the value of the turbulent viscosity coefficient ν t (see Eggleton et al. 1998, Eq. 113). We assumed, that ν t is constant, i.e., does not depend on ω t , which is not true. For small ω t , ν t is nearly constant, while for larger tidal frequency, it is usually believed to be proportional to the inverse of its square, (see, e.g., Goodman and Oh 1997;Terquem et al. 1998). Moreover, ν t is not isotropic, and its value depends on the Coriolis number (Kueker et al. 1994). While for the star of the known internal structure, it is possible to calculate ν t , for the planet it is usually not possible due to very poor knowledge of its interior. In our simplified model, we do not take this effect into account.
The equilibrium tide is not the only contribution to the energy loss. The other mechanism relies on the damping of the internal waves excited by tidal force (the so-called dynamical tide). The waves may be damped by radiative diffusion, convective viscosity as well as nonlinear breaking (see, e.g., Barker and Ogilvie 2010, and references therein). This mechanism produces different dependence of the dissipative parameter, then the equilibrium tide. In the literature, the so-called quality factor Q (or rather the modified quantity Q ) is used to express the energy dissipation rate. If we find formal relation between σ and Q, the equilibrium tide model leads to Q ∝ ω t . Dynamical tide model, on the other hand, gives the inverse relation. The dependence is in fact much more complex and is a subject of many studies in the literature (Zahn 1975;Goodman and Dickson 1998;Goodman and Lackner 2009;Ogilvie 2009;Barker and Ogilvie 2010;Efroimsky 2012).

Averaging over the mean anomaly
The orbital equations of motion have the following general form: where the first term is the Keplerian acceleration and f is a perturbation of small magnitude as compared to the Keplerian term. Therefore, the planet moves on weakly perturbed elliptic orbit. The time scales of variation of its elements are a few orders of magnitude longer than the orbital period. Hence, one might perform the averaging of the equations of motion over fast evolution (i.e., over the mean anomaly) making use of the averaging principle [see, e.g., Arnold et al. (2006), §6.1 or Arnold (1989), §51, §52]. Nevertheless, at first we should verify whether the evolution of Ω Ω Ω 0 and Ω Ω Ω 1 occur in slow-enough time scale. Because the magnitude of the moment of inertia of the planet is much smaller than that one of the star, i.e., I 1 I 0 , the planetary angular velocity vector varies faster than the stellar one. With the help of a relatively simple analysis of (50) 4 , one finds that the relative variation of Ω Ω Ω 1 during one orbital period has its order of: where n is the mean motion. For a one-day orbit and fast rotating Jupiter-like planet (Ω 1 ∼ n), the above quantity is of the order of 6 × 10 −2 . For wider orbits. it obviously decreases. Thus, we have shown that the assumptions of the averaging principle are fulfilled at acceptable level even for rather "extreme" systems. Nevertheless, we should keep in mind that this level of approximation is only acceptable when one would like to compare the results of the mean and exact model. Having the equations of motion in general form of (51), one may find the following equations for the evolution of the orbital angular momentum and Laplace vectors, respectively h and e:ḣ Keplerian orbital elements, such as the semi-major axis, eccentricity, inclination, longitude of ascending node and argument of pericenter, may be easily derived from the components of h and e. We choose this representation instead of the classical Gauss planetary equations, because it has no singularity with respect to the inclination. If the motion is considered to as the Keplerian one, all angles h, e, Ω Ω Ω 0 and Ω Ω Ω 1 are constant and the fast vectors r andṙ may be expressed as the following combinations of e and h × e: where x orb and y orb are Cartesian (x, y)-coordinates in the Kepler frame, i.e, x orb = r cos ν and y orb = r sin ν, where ν is the true anomaly. Thus, the right-hand sides of the equations forḣ,ė,Ω Ω Ω 0 andΩ Ω Ω 1 are some functions of ν, which are parametrised by constant vectors h, e, Ω Ω Ω 0 and Ω Ω Ω 1 . According with the principle of averaging 5 , the equations of motion are averaged out over the mean anomaly M, under assumption that the orbital elements are fixed over the averaging period. Formally, the procedure relies on calculating following integral: where F F F is averaged function. In the above formulae, we changed the integration variable, from the mean anomaly to the true anomaly. Note that we know F F F as a function of ν rather then M [see the discussion in Migaszewski and Goździewski (2008)]. By calculating such integrals over the right-hand sides of the equations of motion, one finds the secular equations that have the following form: Formally, variable names ofḣ,ė, etc., should be encompassed by square brackets . . . , that denote the averaged values. However, it will be clear from the context, that after the averaging all functions should be understood as the secular or the mean quantities. The functions F F F are given by the following formulae: All these functions are equal to 1 for circular orbits and are greater than 1 for elliptic orbits.To verify the correct form of the averaged equations, we compare the time-evolution of the unaveraged system, Eqs. (49) and (50), with the evolution of the secular system, Eq. (55)-(57). . As we can see, the results agree very well, although the precession of the planetary spin is relatively fast, with the period of about 100 days. We will discuss further in this paper, that this time-scale is the fastest one after the orbital period, and the next one, in terms of the magnitude, is the characteristic time-scale of the pericenter advance. The experiment described above concerns the conservative system, with λ 0 = λ 1 = 0. The next Fig. 2 shows the results for the non-conservative models, with λ 0 = 10 −2 and λ 1 = 10 −2 (the parameters are unrealistically large, just to shorten the computations time). The left-hand panel is for the evolution of the ratio Ω 1 /n. As we will show later on, for an eccentric orbit, the rotational frequency does not converge exactly to the mean motion n, but rather to some larger value. The left-hand panel illustrates the time evolution of θ 1 (angle between Ω Ω Ω 1 and h). Clearly, it decreases to zero. Both these quantities evolve in relatively short time-scale; let us recall that λ 0 , λ 1 parameters are by 3-4 orders of magnitude too large in this test as compared to likely values of the real systems. Again, the numerical test shows a perfect agreement between the results derived from both models. Eccentricity e and the semimajor axis a evolve in longer time-scale. Figure 3 shows plots of e(t) (the left-hand panel) and a(t) (the right-hand panel) for the same parameters as taken in the previous experiment. Again, the agreement between the exact and the mean models are excellent.
As we have seen on Fig. 1, the planetary spin evolves in the time-scale that is only two orders of magnitude longer than the orbital period. It is still short, as compared to the full time-scale counted in Gyrs. The integration step would be then still too short to study the long term dynamics of the system CPU-effectively. In the next section, we will show that angle φ 1 may be treated as the second fast angle, and it is possible to perform the second averaging of the secular system.

Time-scale of the evolution and the second averaging
Our next goal is to estimate the time-scale of the evolution of the conservative model. In this section we fix σ 0 = σ 1 = 0. That implies the following properties of the equations of motion: It means that the conservative system possesses at least four first integrals e, h, Ω 0 , Ω 1 . The constant eccentricity and the magnitude of orbital angular momentum imply also the semimajor axis (a) constant. Moreover, it may be easily shown, that the total angular momentum L = βh+ I 0 Ω Ω Ω 0 + I 1 Ω Ω Ω 1 is a constant vector. Another constant quantity is e·h = 0. Therefore, the initial set of 12 scalar equations of motion may be reduced with the help of the 8 first integrals, so it should be possible to transform it to the set of four scalar equations. Nevertheless, we do not attempt to write them down, but rather to make use of the properties of the conservative secular system, in order to simplify the non-conservative system. At first, let us make a simple estimation of the time scale of the evolution of vectors h, e, Ω Ω Ω 0 , Ω Ω Ω 1 in the conservative system. We choose the following (say, representative) values of its parameters: m 0 = 1 m , m 1 = 1 m J , a = 0.02 au, R 0 = 1 R , R 1 = 1 R J , Ω 0 = 2 π/(5 days), Ω 1 = 2 π/(1 day), k L,0 = 0.03, k L,1 = 0.3, κ 0 = 0.2, κ 1 = 0.5. We also take into account the dominant term only. We obtain: As we have shown, Ω Ω Ω 1 becomes the fast vector after the first averaging. Then, we may express this vector in the Keplerian frame, assuming that h, e are constant vectors. We obtain: where θ 1 is an angle between h and Ω Ω Ω 1 , and φ 1 is the azimuthal angle of Ω Ω Ω 1 in the orbital plane, measured from the direction of vector e. Using Eq (57) for l = 1, one may find: The characteristic time-scale for φ 1 may be estimated as follows: The precession rate is constant in the conservative model. Nevertheless, for θ 1 ≈ π/2, φ 1 may not change fast enough to be considered as the fast angle. That introduces a limitation of the analysis. To conclude, we stress that in some cases angle φ 1 may be not slow enough to make the averaging over the mean anomaly valid (that corresponds to too large, too close, too fast rotating planet), or it may be not fast enough to make it possible to perform the second averaging over itself (Uranus-type orientation of the spin vector). As we will show further on, the last limitation is weaker than the first one. After the second averaging, we obtain the following equations of motion: After the second averaging, the fastest angle is the argument of pericenter. Thus the time step-size of the integration increases significantly. To verify the correctness of this step, we perform a new experiments, the results of which are illustrated in Fig. 4. The parameters of the system are the same as in the previous tests. At this time, the plots show the evolution of Ω 0 /n (the left-hand panel) and θ 0 (the right-hand panel). The agreement between the first averaged (black dots) and the second-averaged model (grey curve) are basically perfect.

The third averaging and the dissipative equations of motion
From Eqs. (66)-(68) one can obtain equations governing evolution ofȧ,ė,Ω 0 ,θ 0 and alsȯ Ω 1 ,θ 1 . They read as follows: Ω 1 = 9 4 I 1 σ 1 A 2 1 n a 6 1 − e 2 6 2 E 2 cos θ 1 − E 4 1 − e 2 3 2 Ω 1 n 1 + cos 2 θ 1 , where These equations do not depend on the longitude of the ascending node nor on the azimuthal angle of Ω Ω Ω 0 . A simple estimation of the time-scale of the equations of motion derived in the previous section shows that, after the second averaging, the argument of pericenter ω becomes the fastest angle in the system. It varies typically in a time-scale of 10 4 years. To be more specific, and to find limitations to the choice of ω as the fast angle, we obtain a direct form ofω in the second averaged system (the conservative system). It is a sum of five terms, i.e., ω =ω 3,0 +ω 3,1 +ω 4,0 +ω 4,1 +ω 5 , where particular terms of the sum correspond to functions F F F 3,0 , F F F 3,1 , F F F 4 (these are contributions form the star and the planet), F F F 5 , respectively. We obtain:ω We introduce a new quantity: which is the ratio of the magnitude of the orbital angular momentum to the magnitude of the stellar rotational angular momentum. For typical parameters considered in this work, α is of the order of unity. Contributions form the tidal deformation of the objectsω 4,0 ,ω 4,1 , as

The critical inclinations
We may now introduce the critical inclinations θ (crit) 0 and θ (crit) 1 . The first quantity is a solution to the equationω 3,0 (θ 0 ) = 0, the second one is a solution to the equationω 3,1 (θ 1 ) = 0. One can easily obtain the following expressions: The critical inclination θ (crit) 0 has the following limits: as a function ofᾱ ≡ α/(1 + α). In this paper we consider systems with α ∼ 1 (orᾱ ∼ 0.5), thus magnitudes of the stellar and orbital angular momenta are of the same order. Note that for θ 0 ∼ θ (crit) 0 , a contribution toω 3,0 emerging due to the rotational deformation of the star is very small. If it was the only contribution, ω could not be considered as the fast angle. Similarly, if θ 1 ∼ θ (crit) 1 , the contribution of the rotational deformation of the planet is negligible.

The time-scale of the rotation of pericenter
To justify ω as the fast angle in the second-averaged system, we estimate the time-scale τ 3,l ≡ 2 π/ω 3,l , τ 4,l ≡ 2 π/ω 4,l ( l = 0, 1), τ 5 ≡ 2 π/ω 5 . We obtain: where C l ≡ cos θ l and T rot,0 , T rot,1 denote the rotation periods of the star and the planet, respectively. The typical (characteristic) time-scales were calculated for k L,0 = 0.03 and k L,1 = 0.3. The effect of the tidal deformation of the planet dominates over remaining perturbations in the range of the typical parameters. For a wider orbit, τ 4,1 increases faster then τ 3,0 or τ 3,1 . Nevertheless, for a Sun-like star and a Jupiter-like planet, the effect of the rotational deformation of the bodies may be never the most important contribution toω. For a wider orbit, the dominating perturbation is the relativistic term. We have shown, that in the interesting range of physical and orbital parameters of the system, the argument of pericenter always increases (i.e.,ω > 0), and may be treated as the fast angle in the third averaging. Figure 6 illustrates graphs of characteristic time-scale of the evolution of the conservative system. Thin, black curves are for the time-scale τ 3,l , τ 4,l and τ 5 discussed above. The dashed, red curve is for the joint effect of two dominating perturbations, i.e., the general relativity and the tidal deformation of the planet. The time-scale due to the rotational deformation of the object is always longer than the former one, thusω is always positive. This figure ensures us, that the third averaging (over ω), is valid in the considered range of parameters. Nevertheless, for fast rotating (periods of ∼1 day) and large stars (with R 0 ∼ 2 R ), the time-scale τ 3,0 may be shorter for some orbits (e.g., a = 0.1 au) than τ 4,1 and τ 5 .

The time-scale of the evolution of the system
Because of the small magnitude of the moment of inertia of the planet, it has to be verified whether the dissipative time-scale forΩ 1 andθ 1 may be of the same order, as the time-scale of the rotation of periastron. For the completeness of the analysis, we estimate the time-scales of all variables, i.e., a, e, Ω 0 , θ 0 , Ω 1 , θ 1 . Because the right-hand sides of the equations of motion have rather complex form, we limit this estimation to the case of small e, θ and Ω l /n 1. For the typical parameters of the system, one finds the following characteristic time-scales: where for any quantity denoted as X , the time-scale is defined as τ X ≡ X/|Ẋ |. The equations of motion for a and e are sums of terms representing contributions due to the energy dissipation in the star and in the planet, i.e,ȧ =ȧ ( * ) +ȧ (p) ,ė =ė ( * ) +ė (p) . We calculated these time-scales for contributions coming from each body separately. The dissipative evolution of the planetary spin occurs relatively fast. However, one should keep in mind that these estimates were done under the assumption of Ω 1 /n 1. This means that we neglect terms with Ω 1 /n in the equations of motion. This is not correct in general, and that simplification has been used only to derive the time-scales. Moreover, the magnitude of λ 1 is not well known. A concise discussion of the time-scales present in the system may be summarized graphically in Fig. 7. It illustrates a separation of particular time-scales for Keplerian motion, precession of the planetary spin (or frequency associated with the angle φ 1 ), rotation of the pericenter (frequency of variation of ω), and finally, the time-scale of the dissipative evolution of planetary spin. According with formulae derived for τ Ω 1 , we included also the dependency of the mean motion, see Eq. (76) with e = 0, θ = 0. For extreme systems, with a one-day orbital period, the mean motion may be comparable with the period of the rotation of the pericenter, and Ω 1 , θ 1 may be considered as the fast quantities, similarly to ω.
Nevertheless, the equations forΩ 1 andθ 1 do not depend on ω, and the equations foṙ Ω 0 andθ 0 , which depend on ω, do not depend of Ω 1 , θ 1 . Therefore, we can average out the equations for the evolution of the stellar spin, without considering the evolution of the planetary spin, even if it varies in comparable time-scale.

The third averaging and the quasi-synchronization of the planet's spin
As we have shown, the time derivatives of Ω 0 and θ 0 depend on the fast vector e, i.e., they depend on the fast angle ω. One can easily find that P = E 4 sin θ 0 , thus, after the third averaging, the equations forΩ 0 andθ 0 read as follows: The spin of the planet evolves much faster than the other quantities, i.e., a, e, Ω 0 , θ 0 . Assuming that these parameters are constant, one may find that Ω 1 and θ 1 tend to an equilibrium. There exists only one such an equilibrium, which is linearly stable, i.e., Thus, the rotational velocity of the planet tends to n only for circular orbits, while for e > 0, the equilibrium value is larger than the mean motion. That is why we call this state as quasisynchronous, rather than the synchronous one. For e > 0, the energy is still dissipated in the planetary interior, even if the spin has reached the equilibrium orientation. Time-scales τ a, p , τ e, p are relatively short, ∼ 10 6 years for a one-day orbit. Yet they were obtained under the assumption that the planetary rotational velocity remains separated from the quasi-synchronous state, which is true only for t < τ Ω 1 . When the planet is located at this state, its contribution toȧ andė may be smaller. Putting the equilibrium values for Ω 1 and θ 1 to Eqs. (72) and (73), we find the following expressions: However, the tides in the planet cause a fast orbital decay; it acts in relatively short time-scale of the order of τ e, p . After circularization of the orbit, the tides in the planet do not contribute toȧ. We notice here, that in multi-planet systems the eccentricity as well as inclination of the orbit vary in the conservative time-scale and the spin of the planet may unlikely reach the equilibrium. Therefore, only in single-planet systems, it is possible to use the simplified equations of motion (i.e. those for the equations with fixed θ 1 = 0 and Ω 1 = Ω 1 | eq ).

The dissipative evolution of the system
The final set of the equations of motion, i.e., Eqs. (87), (88), (90), (91) has an integral of the magnitude of the total angular momentum, L ≡ |L|, where To describe the evolution of the system, one has to solve three equations of motion, i.e. those forȧ,ė andΩ 0 , for a fixed value of L. Still, we cannot write down the solution to these equations, and we have to integrate them numerically. Moreover, we may study some particular solutions, like the equilibria, analytically.

The stability of the equilibrium
The final approximation of the evolutionary equations of the system possesses one type of equilibrium. It corresponds to a state, in which the mechanical energy is not lost, i.e., e = 0, Ω 0 = n * = μ/a 3 * , θ 0 = 0, where a * (or n * ) may be treated as a parameter of a family of such equilibria. For some a * , the equilibrium may be stable, while for some other value it may be unstable. To study the stability, we write down the linear variational equations in the vicinity of the stationary solution as follows: where δ a , δ e , δ Ω 0 and δ θ 0 are the variations of a, e, Ω 0 , θ 0 , respectively. Equations for δ e and δ θ 0 are separable, and they admit simple solutions, i.e., both these quantities decrease exponentially to 0: The equations for δ a and δ Ω 0 are conjugate each to other. To solve them, we define a new quantity x ≡ Ω 0 /n and write down the relevant variational equation: This equation also admits a simple solution: δ where λ x may be either positive or negative. The stability condition reads as follows: For typical values of the physical parameters, the equilibrium is stable if The orbital period corresponding to a = 0.074 au is ∼ 7.4 days. Thus, for a Jupiter-like planet orbiting a Sun-like star, with orbital period shorter than ∼ 7.4 days, the planet does not tend to the equilibrium. When the equilibrium is unstable, the initial condition Ω 0 < n forces the planet to fall down into its host star.
Considering the time-scale of the tidal evolution of orbits with semi-major axes larger than this critical value, one may conclude that it is unlikely that a Jovian planet reaches the equilibrium. Moreover, in a more realistic model, this equilibrium does not exist at all (Barker and Ogilvie 2009), because additional effects, like the angular momentum dissipation due to stellar winds, may decrease the rotational velocity of the star. Figure 8 illustrates the evolution of the system initially located in the vicinity of the equilibrium. Black dots are for the solutions to the equations of motion of the second averaged system, the grey curve is for the third-averaged system with the spin of the planet in the synchronous state. This test shows that the model is correct. It also illustrates the behavior of the system near the unstable equilibrium. The system with initially Ω 0 n excites both the semi-major axis and the Ω 0 /n ratio. Still, the rates of variations of a and Ω 0 /n do not suppress dissipation of the total energy. On the other hand, for initial Ω 0 n, the planetary destiny is to fall down into the parent star.

Parametric study of the dissipative evolution
Here, we present the results of the analysis of the evolution of the system for various initial conditions. We adopt typical physical parameters of the system as follows. The masses are m 0 = 1 m and m 1 = 1 m J , for the star and the planet, respectively. Their radii are R 0 = 1 R and R 1 = 1 R J . The mass distribution is described by polytropic models with indices n 0 = 3 and n 1 = 1.5. Parameters of the energy dissipation rates are λ 0 = 5×10 −5 and λ 1 = 5 × 10 −7 in the first case (results are presented on Figs. 9, 10) and λ 0 = λ 1 = 5 × 10 −5 in the second case (Figs. 11, 12). We also consider a few sets of initial T rot,0 and θ 0 . Figure 9 illustrates the results derived for the initial rotational period of the star T rot,0 = 5 days, and for three initial values of the inclination between the stellar equator and the orbit (from the top to the bottom row, it is 0 • , 60 • , 120 • , respectively). Colors encode the final a (the left-hand column) and final T rot,0 (the right-hand column). Maps of the final state of the system as a function of the initial semi-major axis and eccentricity. Panels in the left-hand column show the final a while panels in the right-hand column show the final rotation period of the star. Parameters of the system are the following: m 1 = 1 m , m 1 = 1 m J , R 0 = 1 R , R 1 = 1 R J . Polytropic indices of the star and the planet are equal to 3 and 1.5, respectively. The energy dissipation constants are λ 0 = 5 × 10 −5 , λ 1 = 5 × 10 −7 . The initial rotation period of the star is T rot,0 = 5 days. The initial inclination θ 0 is 0 • , 60 • , 120 • for the top, middle and bottom rows, respectively. The integration time is 5 Gyr At first, let us study panel (a), obtained for initial θ 0 = 0. For the initial (a, e) in the white region, the planet falls down into the star in a time shorter than 5 Gyr. The black thick curve is the border of survival. For the initial conditions located at the (a, e)-plane, on the right side of this curve, the planet does not fall into the star during the whole time of integration 5 Gyr. Other black (thin) curves are for initial conditions for which the planet survives during 10 7 , 10 8 and 10 9 years, respectively. As we can see, for the parameters adopted in this experiment, the border is placed at ∼ 0.046 au for initially circular orbits and at ∼ 0.068 au for initial e = 0.5. Clearly, planets with initially larger e migrate towards the star faster. As we already noticed, the contribution toȧ stemming from the energy dissipation in the interior of the planet is typically larger or similar to that due to the dissipation in the star. Nevertheless,  Fig. 10 The same as shown on Fig. 9, but for initial T rot,0 = 2 days this term is proportional to e 2 and vanishes for circular orbit. The eccentricity is damped in similar time-scale than a for corresponding values of λ 0 , λ 1 , thus the dissipation in the planet is important not only at the early stages of the dissipative evolution.
The evolution of the eccentricity may be also read from the discussed figure. Contours plotted with yellow thin curves are for the levels of constant final e, which are 0.001, 0.1, 0.2, 0.3 and 0.4, respectively. As we can see, for initial conditions located closer than ∼ 0.01 au from the "survival curve", on its right-hand side, the eccentricity is damped significantly, while for larger initial a, its variability remains small. For a 0.074 au, the dissipation timescales τ e are longer than the integration time, and the eccentricity is not altered significantly. Also the decay of the orbit is much slower.
The The same as in Fig. 9, but for λ 0 = 5 × 10 −5 , λ 1 = 5 × 10 −5 . The initial rotational period of the star is T rot,0 = 5 days the left-hand side of the border curve), the rotation of the star becomes faster. For the initial conditions near the border, the rotational period decreases from 5 to ∼ 2.2 days. It may be understood, because the initial orbital angular momentum of the planet increases for larger initial a. For initial a 0.074 au, for which the planet survives, and a variation of a is not large, also a change of Ω 0 is not large.
The next two rows show the results for initial θ 0 = 60 • (the middle row) and 120 • (the bottom row). In the case of panels (c) and (d), the differences between these results and the previous ones are not significant. The border of planetary survival shifts slightly towards larger initial a. The minimal values of the final T rot,0 for θ 0 = 60 • increase also slightly, ∼ 3 days near the border. This is an effect of addition of two vectors, i.e. β h and I 0 Ω Ω Ω 0 , which are mutually inclined. The angular momentum exchange that may manifest itself in the decrease of T rot,0 is smaller here than in a coplanar system. The inclination θ 0 does not  Fig. 12 The same as in Fig. 9, but for λ 0 = 5 × 10 −5 , λ 1 = 5 × 10 −5 . The initial rotational period of the star is T rot,0 = 2 days change significantly during the evolution. The yellow contours shown at panel (d) indicate constant levels of the final θ 0 = 50 • , 57 • , 59 • . The systems near the survival border change their θ 0 by about of 10 degrees, while for systems located closer than ∼ 0.02 au from this border, on its right-hand side, the inclination is almost constant. The observed slow rate of inclination dumping agrees with the determined inclination in hot-Jupiter systems (e.g., Triaud et al. 2010;Pont et al. 2010).
For the retrograde orbits, i.e., with initial θ 0 > π 2 , here θ 0 = 120 • , the final state of the system is different. The survival border shifts again. But most significant differences may be visible in the final rotational period of the star. The picture differs from those obtained for θ 0 = 0 • and θ 0 = 60 • . The larger final T rot,0 are for those initial conditions for which the planet falls down into the star or has survived but stays close to the border. The difference may be explained due to the z -components of initial β h and I 0 Ω Ω Ω 0 have opposite signs. Dissipative parameters are λ 0 = 10 −5 and λ 1 = 10 −4 . The initial system is close to the equilibrium (i.e., e ≈ 0, θ 0 ≈ 0, Ω 0 ∼ 1). Grey curves are for the initial Ω 0 /n < 1, the black curves are for Ω 0 /n > 1 Again, the inclination of configurations which end with a > 0 does not change significantly. The evolution of the eccentricity in all three cases looks like very similar. The main difference between the studied cases concerns the final T rot,0 . Figure 10 illustrates the results derived for parameters used in the previous experiment; only the initial rotational period of the star is now T rot,0 = 2 days. The view of the final state of the system is similar to that one presented in Fig. 9. There is a difference, however, that relies in a shift of the survival border towards smaller initial semi-major axes for θ 0 < π 2 , and towards larger a for θ 0 > π 2 . It may be understood rather easily. For the faster rotation of the star, the border Ω 0 = n shifts at the (a, e)-plane towards larger initial a. A contribution proportional to the ratio Ω 0 /n implies an increase of a, e, θ 0 and a decrease of Ω 0 ; it acts opposite to the remaining terms. Larger Ω 0 /n, for some selected initial condition, means that this term dominates over the other terms. It manifests itself also in the fact that for initially coplanar configurations, there exist regions in (a, e)-plane, in which the eccentricity is excited to larger values. For systems with θ 0 > π/2 the survival border shifts towards smaller a. It is quite clear, because the term Ω 0 /n is multiplied by cos θ 0 which is negative for the retrograde orbit. This term causes a decrease of both a and e, and has larger magnitude for greater Ω 0 (shorter T rot,0 ).
In the next two figures (Figs. 11 and 12), we present a similar study to the previous experiments, but derived for larger values of the energy dissipation constant λ 1 = 5 × 10 −5 . The results are illustrated in the same manner. The only difference relies now in significantly smaller final eccentricity than it was derived in the previous experiments. Again, it may be explained if we recall that the contributionė p increases by two orders of magnitude with respect to the previous tests. For large λ 1 it dominates overė * and forces a fast dumping the eccentricity. On the other hand, the termȧ p is important only at the early stages of the evolution, when e is significantly different from 0. When e becomes very small, the termȧ p does not accelerate the orbit decay. The borders of survival as well as maps of the final T rot,0 are very similar to those ones obtained for smaller λ 1 .
It should be noted here, that in none of the cases illustrated in Figs. 9-12, the system ended in the synchronous state. For a 0.074 au this equilibrium is unstable, while for a 0.074 au the time-scales of the evolution are too long to fix the system in this state. Figure 13 illustrates the evolution of the system which is initially close to the equilibrium for a = 0.08 au (the rotational period of the star T rot,0 ∈ [5.4, 9.4] days). Although the equilibrium is stable, the system tends towards this state very slowly; the time on the x-axis is counted in terayears! Moreover, for initial Ω 0 /n < 1 (the grey curves), the planet likely will fall down onto the star, even if at the beginning the system seems tend to the equilibrium. It may be explained relatively easily. If the initial Ω 0 /n < 1, then for a 0.074 au this ratio increases, tending to unity, and simultaneously a decreases. When a decreases below the limit of stability, the equilibrium becomes unstable and the system evolves outwards from the equilibrium.

Summary and conclusions
In this work we revisited the problem of the generalized model of a single-planet system including the mechanical energy dissipation. We derived the equations of motion from the very basic formulation of the Lagrange equations of the second kind. In that approach, the potential is a function of the angular velocities of the bodies. Our derivation is different from the standard approach used to derive the rotational equations of motion of the rigid body, in which the potential depends only on the Euler angles describing the orientation of the rigid body, and does not depend on their time derivatives. We obtained the equations possessing a different general form, but in the particular case considered here they remain very similar to the Eulerian equations. Moreover, as we show, an additional term appearing in the right-hand side of the equation forΩ Ω Ω l does not introduce any secular contribution.
The derived equations of motion were averaged out over time-scales corresponding to the conservative evolution of the system. We analyzed all characteristic time-scales, and we showed that they may be ordered in terms of a hierarchical set of variables, which then may be treated as fast variables in a recursive averaging process. The fastest component of the evolution is the Keplerian mean motion of the planet. After the averaging over the mean anomaly, we obtained the first-averaged system. This set of equations describes the evolution in which the fastest variable is the precessing angular momentum vector of the planet. Then the second averaging was performed over the azimuthal angle of Ω Ω Ω 1 in the orbital reference frame. In this way we obtained the second averaged system. Finally, the third averaging was performed over the argument of pericenter. Thus, to eliminate secular variability that occurs in the conservative time-scale, and to obtain the dissipative equations of motion, we had to average out the system over three fast angles, i.e, M, φ 1 , ω. The precision of the averaging has been tested at all stages of the procedure, and we demonstrate that this approach leads to correct results.
Next, we have shown that the dissipative evolution of the angular momentum of the planet occurs much faster than the orbital evolution as well as a variability of the stellar angular momentum. The inclination of Ω Ω Ω 1 with respect to h decreases to 0 in a time-scale only slightly longer than the characteristic time-scale of the pericenter rotation, for a close-in planet. Simultaneously, Ω 1 tends to an equilibrium value which is ≥ n, where the frequencies are equal only for the circular orbit. The equilibrium is stable and we can fix θ 1 and Ω 1 at their equilibrium values. In this way we derived the final set of equations governing the long term evolution of the system that admits energy dissipation, Eqs. (87), (88), (90) and (91).
The obtained equations of motion describe a dynamical system with one equilibrium corresponding to the circular, coplanar and synchronized orbit, which is unstable for orbits inside certain critical distance from the star, and which is stable for orbits outside this border. For a Sun-Jupiter system that border is on a ∼ 0.074 au. Studying the evolution of the system for generic values of physical parameters, we have shown however, that this equilibrium is unlikely to be reached. As we demonstrated, the final state of the system is a complex function of the initial conditions and physical parameters of the model. The derived equations make it possible to study this problem very effectively, both analytically and in terms of the CPU overhead, because all the evolution of dynamical variables occurring in the intermediate time-scale related to the conservative corrections have been averaged out.
Acknowledgments I would like to thank Benoît Noyelles, Michael Efroimsky and the anonymous reviewer for the informative reviews that improved the manuscript and Sylvio Ferraz-Mello for comments on the averaging theory. Many thanks to Krzysztof Goździewski for a discussion and corrections of the manuscript. This work was supported by the Polish Ministry of Science and Higher Education grant N/N203/402739. The author is a recipient of the stipend of the Foundation for Polish Science (programme START, editions 2010 and 2011). This research was carried out with the support of the "HPC Infrastructure for Grand Challenges of Science and Engineering" Project, co-financed by the European Regional Development Fund under the Innovative Economy Operational Programme.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Appendices Appendix A: proof thatĖ given by formulae (45) fulfills the condition (2)
In the considered case, the left-hand side of condition (2) reads as follows: It is quite easy to show, that the following equality occurs, whenĖ is given by Eq. (45): Now, to prove a consistency of the Lagrange equations with a dissipative term in the considered problem, we have to show thaṫ s · ∂Ė ∂ṡ = Ω Ω Ω 0 · ∂Ė ∂ Ω Ω Ω 0 (and similarly forṗ and Ω Ω Ω 1 terms). It is again quite elementary and may be checked by using the above formulae for ∂Ė/∂ṡ and the relation between Ω Ω Ω andṡ (Eq. 4) as well as the fact, that the following equality occurs for unit quaternion s (i.e., s 2 0 + s 2 = 1):