© ESO 2020

1 Introduction

Among the main characteristics of the nine well-characterised circumbinary (CB) systems discovered so far by the Kepler mission, probably the most remarkable one is the proximity of the planets to the binaries. This feature implies that most of the circumbinary planets (CBPs) are located very close to the stability limit (Holman & Wiegert 1999) where, in addition, the in situ formation is very unlikely due to perturbations of the binary (Lines et al. 2014).

While it is widely accepted that migration by interaction with a protoplanetary disc took place in these systems (Dunhill & Alexander 2013), it is not yet clear what role tidal forces played in the subsequent dynamical evolution, which have the peculiarity of being exerted by two central objects of comparable masses and may have had an important effect in such compact systems.

In Zoppetti et al. (2019a) we presented a self-consistent weak friction model for tides on a CBP. We used Kepler-38 system as a testing case and found several interesting features: for example the planet synchronised its spin in a sub-synchronous equilibrium state and, on longer timescales, it migrated outward from the binary. We also derived analytical expressions for the variational equations of the orbital and spin evolution based on low-order elliptical expansions in semimajor axis ratio and eccentricities. This analytical approach was reduced to the two-body case when one of the central masses was taken equal to zero and it reproduced the results of our low-eccentric numerical simulations very well. However, the low-order expansion in eccentricities of our analytical model did not allow its direct application to the eccentric Kepler systems such as Kepler-34 (Welsh et al. 2012) and Kepler-413 (Kostov et al. 2014).

In this work, we present and discuss an extended version of our self-consistent weak friction (Constant Time Lag; CTL) tidal model which is valid for any binary and planetary eccentricity. We apply the results to the well-characterised nine Kepler systems in an attempt to reproduce their past tidal evolution, as well as unveiling the role played by tidal forces in determining their current orbital parameters.

This paper is organized as follows. In Sect. 2, we review the construction of the variational equations for the spin and orbital evolution in a self-consistent manner. In Sect. 3, we present the analytical expressions of the variational equations in closed form for the eccentricities and compare their predictions with results obtained from numerical simulations, as well as with our original low-order expanded model. In Sect. 4, we apply our model to the Kepler CB planets by computing the characteristic timescales for their evolution and by reproducing the tidal flow around each system. A simple geometrical interpretation of the main results of our model is offered in Sect. 4.2 to better understand these dynamical characteristics. Finally, Sect. 5 summarises and discusses our main results.

2 Review of the model

We return to the problem of a binary system with stellar masses m₀ and m₁, and a CBP with mass m₂. The binary and the planet share the same orbital plane and their spins are perpendicular to it. We further assume that each body is an extended mass with physical radius ${\mathcal{R}}_i$ and deformable due to the tidal effects of its companions.

We again assume that the internal dissipation due to tides in each of the three bodies is governed by the CTL and is, therefore, subject to the description suffested by Mignard (1979). In our previous work (Zoppetti et al. 2019a) we demonstrated that, in absence of any mean-motion commensurability between m₁ and m₂, the cross tides can be neglected and the only tidal forces that need to be taken into account on each body m_i are those caused by the deformation that the mass m_i itself induces on its companions. For this reason, in a multi-body tidal model like the one we are studying, the tides on each body can be simply computed by summing the forces and torques between deforming-deformed pairs as given in Eqs. (5) and (6) of Mignard (1979).

In an inertial frame with arbitrary origin, where the position of each body m_i is denoted as R_i, the orbital evolution equations may be expressed as the sum of the point-mass and tidal interactions, which is $$\begin{array}{ll}{m}_i{\ddot{R}}_i\hfill & ={\sum }_{j=0,j\ne i}^2\left(\frac{\mathcal{G}{m}_i{m}_j}{|{\Delta }_{ji}{|}^3}{\Delta }_{ji}+(F_{ij}-F_{ji})\right),i=0,1,2\hfill \end{array}$$(1)

where $\mathcal{G}$ is the gravitational constant and Δ_ij ≡R_i −R_j. We denoted by F_ij the tidal force acting on m_i due to the deformation that it induces on m_j and, following Ferraz-Mello et al. (2008), we also took into account the reacting forces when computing the tides on each mass m_i. According to Mignard (1979), these can be expressed by $$F_{ij}=-\frac{{\mathcal{K}}_{ij}}{|{\Delta }_{ij}{|}^{10}}\left[2({\Delta }_{ij}\cdot {\stackrel{.}{\Delta }}_{ij}){\Delta }_{ij}+{\Delta }_{ij}^2({\Delta }_{ij}\times {\Omega }_j+{\stackrel{.}{\Delta }}_{ij})\right]$$(2)

where Ω_j is the spin vector of m_j (assumed parallel to the orbital angular momenta) and the magnitude factor of the tidal force ${\mathcal{K}}_{ij}$ is given by $${\mathcal{K}}_{ij}=3\mathcal{G}{m}_i^2{\mathcal{R}}_j^5{k}_{2j}\Delta t_j,$$(3)

where k_2j and Δt_j are the second degree Love number and the time-lag of the body m_j. Note that we have neglected the contribution of the static tides since, in the absence of any mean-motion resonance, it does not introduce any angular momentum exchange and therefore yields no secular change in the spins nor in the orbital evolution.

While Eq. (1) governs the orbit, the rotational dynamics can be deduced from the conservation of the total angular momentum. Assuming that the variation in the spin angular momenta of the body m_i is only due to the contributions associated to its own deformation, we can uncouple the spin evolution equations of each body to obtain $$\frac{\text{d}{\Omega }_i}{\text{d}t}=\frac{1}{{\mathcal{C}}_i}{\sum }_{j\ne i}\frac{K_{ji}}{|{\Delta }_{ji}{|}^6}[\frac{{\Delta }_{ji}\times {\stackrel{.}{\Delta }}_{ji}}{|{\Delta }_{ji}{|}^2}-{\Omega }_i],i=0,1,2$$(4)

where ${\mathcal{C}}_i$ is the principal moment of inertia of m_i.

3 Analytical secular model for arbitrary eccentricities

To construct the analytical model of the orbital evolution Eq. (1) and the rotational evolution Eq. (4), we first adopt a Jacobi reference frame for the position and velocities of the bodies. In this system, the Jacobi position vectors of the masses ρ_i can be expressed in terms of their inertialcounterparts R_i as $$\begin{array}{ll}{\rho }_0\hfill & =\frac{1}{{\sigma }_2}(m_0{R}_0+m_1{R}_1+m_2{R}_2),\hfill \\ {\rho }_1\hfill & =R_1-R_0,\hfill \\ {\rho }_2\hfill & =R_2-\frac{1}{{\sigma }_1}(m_0{R}_0+m_1{R}_1),\hfill \\ \hfill & \hfill \end{array}$$(5)

where ${\sigma }_i={\sum }_{k=0}^i{m}_k$. The transformation equations for the velocities are analogous.

To construct an analytical secular model valid for any binary and planetary eccentricity, we note that both Eqs. (1) and (2) depend on powers of the inverse of the relative distances, i.e. $|{\Delta }_{ij}{|}^{(-n)}\equiv |R_j-R_i{|}^{(-n)}$, with n a positive integer. Assuming momentarily that t = R_i∕R_j < 1, we can write $$g^{(n)}(t,x)\equiv {\left(\frac{R_j}{\left|{\Delta }_{ij}\right|}\right)}^n={(1+t^2-2tx)}^{-n/2}$$(6)

where x is the cosine of the angle between both position vectors. The expansion of g⁽¹⁾ (t, x) in power series of t is well known $$g^{(1)}(t,x)={\sum }_{k=0}^{\infty }{t}^k{P}_k(x)$$(7)

with P_k(x) the k-degree Legendre polynomial. In a similar manner, for n > 1 we can introduce “generalized” polinomials $Z_k^{(n)}(x)$ such that $$g^{(n)}(t,x)={\sum }_{k=0}^{\infty }{t}^k{Z}_k^{(n)}(x)$$(8)

such that $Z_k^{(0)}(x)=1$ and $Z_k^{(1)}(x)=P_k(x)$ for all values of k ≥ 0. These new polynomials can be obtained for the original Legendre functions through recurrence relations. Writing g⁽ⁿ⁾ (t, x) = g⁽ⁿ⁻¹⁾(t, x) ⋅ g⁽¹⁾(t, x), replacing each function by its series expansion and equating equal powers of t, we find $$Z_k^{(n)}(x)={\sum }_{j=0}^k{P}_{k-j}(x)Z_j^{(n-1)}(x).$$(9)

Introducing these expressions into the equations of motion, transforming the position and velocity vectors to Jacobi variables, and averaging over the mean longitudes λ₁ and λ₂, we find closed analytical expressions in terms of the eccentricities e₁ and e₂. The procedure is analogous to that followed by Mignard (1980) in the two-body case, with two important exceptions. First, we consider expansions in α = a₁∕a₂ up to fourth order (i.e. k = 4); second, we do not average over ϖ₁ and ϖ₂. While this would greatly simplify the resulting equations, in Zoppetti et al. (2019a) we observed that in some simulations, typically associated to low planetary eccentricity, the planet and the secondary star entered in the aligned secular mode Δϖ = ϖ₂ − ϖ₁ ~ 0. Thus, the perturbation terms dependent on Δϖ may have a net contribution in the long-term dynamics, at least for some initial conditions.

The complete procedure was carried out with the aid of an algebraic manipulator. Some intermediate expressions had to be expanded up to fifth order in t = R_i∕R_j to avoid loss of precision in the final results. The process was complex, especially for k > 2 where we found little benefit explicitly working with the $Z_k^{(n)}(x)$ functions. However, we keep their definition and recurrence relations since they could be useful for quadrupole-level theories.

Fig. 1 Time derivative of the planetary spin rate Ω2 in our generic binary system as a function of α = a1∕a2. Each panelassumes different eccentricities for both orbits. Broad gray curves show results from semi-analytical calculationsand predictions from analytical models are depicted in thin lines. Values obtained with our original (Zoppetti et al. 2019a) are presented in blue, while the colors correspond to our current model with different orders in α. See inlaid box in left-hand plot for details.

3.1 Planetary spin evolution

The process discussed previously was applied to both the orbital and spin variational equations. In the latter case, up to fourth order in semimajor axes ratio, we obtain $$<\frac{\text{d}{\Omega }_2}{\text{d}t}>=\frac{3\mathcal{G}{\mathcal{R}}_2^5{k}_{2,2}\Delta t_2}{2{\mathcal{C}}_2{a}_2^6}{\sum }_{i=0}^4{K}_i^{(s)}(A_i^{(s)}{n}_1+B_i^{(s)}{n}_2+C_i^{(s)}{\Omega }_2){\alpha }^i,$$(10)

where n₁ is the mean motion of the binary, n₂ that of the planet, and $$K_i^{(s)}={\gamma }_0^i{m}_0^2+{\gamma }_1^i{m}_1^2,{\gamma }_0=\frac{m_1}{{\sigma }_1},{\gamma }_1=-\frac{m_0}{{\sigma }_1}.$$(11)

The coefficients $A_i^{(s)}$, $B_i^{(s)}$ and $C_i^{(s)}$, functions of the eccentricities and difference of longitudes of pericenter, are listed in Appendix A.

To test the accuracy of our analytical model, we consider a generic circumbinary system with primary mass m₀ = 1 M_⊙ and physical radius ${\mathcal{R}}_0=1{\mathcal{R}}_{\odot }$, plus a secondary star with mass m₁ = m₀∕3 and a physical radius ${\mathcal{R}}_1={(m_1/M_{\odot })}^{0.9}{\mathcal{R}}_{\odot }$ (Demircan & Kahraman 1991). The mass parameter of the binary is thus $\overline{\mu }=m_1/(m_0+m_1)=0.25$ and the semimajor axis was set to a₁ = 0.15 AU. For the planet, we chose m₂ = 10 M_⊕ and ${\mathcal{R}}_2=4{\mathcal{R}}_{\oplus }$. We assumed a planetary principal moment of inertia equal to ${\mathcal{C}}_2=0.25m_2{\mathcal{R}}_2^2$ and a modified tidal dissipation factor of $Q_{2^{\prime }}=3/(2n_2{k}_{2,2}\Delta t_2)=10$.

Figure 1 shows the time derivative of Ω₂ as a function of α = a₁∕a₂ for three different sets of eccentricities. Different curves are used to show results using different models. Values obtained by numerically averaging Eq. (A.1) are depicted in broad grey curves, the analytical model from Zoppetti et al. (2019a) is shown in blue curves, and different variations of our current model are presented in other colors. We observe that our new α⁴ model fits the numerical values very well up to α ~ 0.35. The accuracy of our closed-form expressions in the eccentricities is evident when we analyse cases in which at least one of the orbits is moderately or highly eccentric. These situations were beyond the scope of our previous model. The low-eccentricity regime is well reproduced in all models, although again we obtain more precision with the closed fourth-order approach.

Curiously, due to the complexity of the power series, results obtained with α³ are not always better than those obtained with a quadrupole approximation. Moreover, in some cases (e.g. the central panel) both of these models yield predictions which are less accurate compared to those obtained from Zoppetti et al. (2019a). These results show the importance of using a high order expansion in semimajor axis ratio, even for systems that are not particularly compact.

The value of the stationary spin rate $<{\Omega }_2{>}^{(\text{st})}$ predicted by our model can be calculated equating Eq. (A.1) to zero, which yields: $$<{\Omega }_2{>}^{(\text{st})}=-\frac{{\sum }_{i=0}^4{K}_i^{(s)}{D}_i^{(s)}{\alpha }^i}{{\sum }_{i=0}^4{K}_i^{(s)}{C}_i^{(s)}{\alpha }^i},$$(12)

where we have denoted $D_i^{(s)}=A_i^{(s)}{n}_1+B_i^{(s)}{n}_2$. Expanding this expression in power series of α up to fourth order gives $$<{\Omega }_2{>}^{(\text{st})}={\sum }_{i=0}^4<{\Omega }_2{>}_i^{(\text{st})}{\alpha }^i,$$(13)

where the coefficients $<{\Omega }_2{>}_i^{(\text{st})}$ are again listed in Appendix A.

Figure 2 shows the stationary spins (divided by n₂), as function of α, calculated with the same analytical models considered in Fig. 1. To avoid too much clutter, we did not plot results corresponding to the α³ model. We also considered two different values of the binary mass parameter: results for $\overline{\mu }=0.1$ are shown in the left panels while those obtained for $\overline{\mu }=0.5$ are graphed in the right panels.

Comparisons with numerical estimations (broad gray curves) show that our current α⁴ model predict reliable stationary solution for any binary and planetary eccentricity up to α ~0.35. Larger mass ratios appear to require less sophisticated dynamical models, while at least a α⁴ model is necessary for low-mass secondary components. We also note that the mass of the secondary star (represented by $\overline{\mu }$) competes against the planetary eccentricity in establishing the stationarity solution: higher secondary masses tend to sub-synchronous solutions while, as in the two-body case, high planetary eccentricity leads to super-synchronous solutions. This fact was previously observed in Zoppetti et al. (2019a) for low-eccentric systems.

Fig. 2 Stationary spin of a CB planet as a function of α. Different rows consider different eccentricity of the orbits (assumed equals). Different columns represent different reduced mass: $\overline{\mu }=0.1$ (left panels) and $\overline{\mu }=0.5$ (right panels). The type of the curves represent the same that in Fig. 1.

3.2 Secular evolution for the planetary semimajor axis and eccentricity

The variational equation for the planetary semimajor axis in the Jacobi frame is given by (Beutler 2005) $$\frac{\text{d}{a}_2}{\text{d}t}=\frac{2a_2^2}{\mathcal{G}{\sigma }_2}(\stackrel{.}{{\rho }_2}\cdot \delta f_2)$$(14)

where δf₂ is the total tidal force (per unit mass) affecting the two-body motion of the planet around the barycenter of the binary system. It is explicitly given by $$\delta f_2=\frac{1}{{\beta }_2}[(F_{2,0}-F_{0,2})+(F_{2,1}-F_{1,2})],$$(15)

where β₂ = m₂σ₁∕σ₂ is the Jacobian reduced mass of the planet. We then repeat the procedure employed in Sect. 3.1, now in order to obtain averaged expressions for the orbital evolution which are in closed form with respect to the eccentricities. Up the fourth order in α and averaging over the mean longitudes, we obtain $$<\frac{\text{d}{a}_2}{\text{d}t}>=\frac{n_2}{\mathcal{G}{m}_2{\sigma }_1{a}_2^4}{\sum }_{i=0}^4{K}_i^{(a)}(A_i^{(a)}{n}_1+B_i^{(a)}{n}_2+C_i^{(a)}{\overline{\Omega }}_i^{*}){\alpha }^i,$$(16)

where $$K_i^{(a)}={\mathcal{K}}_0^{(+)}{\gamma }_0^i+{\mathcal{K}}_1^{(+)}{\gamma }_1^i,{\mathcal{K}}_j^{(+)}={\mathcal{K}}_{2,j}+{\mathcal{K}}_{j,2},$$(17)

and the coefficients $A_i^{(a)}$, $B_i^{(a)}$, and $C_i^{(a)}$, are listed in Appendix B. The variational equation for the semimajor axis depends on the spins through a new complex parameter given by $${\overline{\Omega }}_i^{*}=\frac{({\mathcal{K}}_{2,0}{\Omega }_0+{\mathcal{K}}_{0,2}{\Omega }_2){\gamma }_0^i+({\mathcal{K}}_{2,1}{\Omega }_1+{\mathcal{K}}_{1,2}{\Omega }_2){\gamma }_1^i}{{\mathcal{K}}_0^{(+)}{\gamma }_0^i+{\mathcal{K}}_1^{(+)}{\gamma }_1^i}.$$(18)

Finally, the time evolution of the planetary eccentricity may be found from the orbital angular momentum $$L_2=\sqrt{L_2\cdot L_2}={\beta }_2\sqrt{\mathcal{G}{\sigma }_2{a}_2(1-e_2^2)}.$$(19)

Differentiating over time and extracting the eccentricity term, we may write $$\frac{\text{d}}{\text{d}t}(e_2^2)=\frac{1}{a_2}[(1-e_2^2)\frac{\text{d}{a}_2}{\text{d}t}-\frac{2L_2}{\mathcal{G}{\sigma }_2{\beta }_2}(({\rho }_2\times \delta f_2)\cdot {\widehat{z}}_2)],$$(20)

where ${\widehat{z}}_2$ is the unit vector of the orbital angular momentum. Proceeding in the same manner that for Eq. (B.1), the rate of change of the planetary eccentricity can be approximated as: $$<\frac{\text{d}{e}_2^2}{\text{d}t}>=\frac{n_2}{\mathcal{G}{m}_2{\sigma }_1{a}_2^5}{\sum }_{i=0}^4{K}_i^{(a)}(A_i^{(e)}{n}_1+B_i^{(e)}{n}_2+C_i^{(e)}{\overline{\Omega }}_i^{*}){\alpha }^i.$$(21)

We note that both the tidal magnitude coefficient $K_i^{(a)}$ and the spin functions ${\overline{\Omega }}_i^{*}$ are the same as deduced for the semimajor axis. The coefficients $A_i^{(e)}$, $B_i^{(e)}$, and $C_i^{(e)}$, are listed in Appendix B.

Figure 3 shows the derivative of the semimajor axis (left panels) and eccentricity (right panels), as a function of α, of a CBP orbiting our generic binary system. The spin rates of the stars and the planet were set in their stationary values: ${\Omega }_0={\Omega }_1\simeq n_1(1+6e_1^2)$ (Hut 1980) and ${\Omega }_2=<{\Omega }_2{>}^{(\text{st})}$ (Eq. (13)). We observe that our closed fourth-order model is a very precise approximation of the numerical values up to α ~ 0.35. The need of an analytical model closed in eccentricities is again particularly evident in the high-eccentricity regime (bottom plots), where our previous model fails not only to estimate the magnitude but may also predict erroneous signs for the time derivative.

Regarding the dynamical implications, we observe that tides always seem to decrease the planetary eccentricity, with magnitude proportionalto e₂ and also to the proximity to the binary. The effect of tides on the planetary semimajor axis is more diverse: for low-to-moderate planetary eccentricities, the tidal forces leads to an outward migration, while inward migration is expected for the high-eccentricity case. These results confirm our previous findings (Zoppetti et al. 2019a) and a simple geometrical explanation is offered in Sect. 4.2.

Fig. 3 Time derivative of the planetary semimajor axis (left panels) and eccentricity (right panels), as a function of α, for a synchronised CB planet around our generic binary system. As before, results obtained with different models are shown in different colored curves: numerical in broad grey, Zoppetti et al. (2019a) in thin blue, and predictions from our current α4 model in black.

4 Application to Kepler systems

Having developed a model valid for arbitrary binary and planetary eccentricities, we can analyse the tidal evolution of Kepler CB systems.

Table 1 lists the physical parameters (mass and radius) of the binaries and the planets discovered by the Kepler mission, as well as the semimajor axes and eccentricities of the binaries (sub-index 1) and the planets (sub-index 2). The last column specifies the reference works from which the previous values were taken. In several cases, the difficulties in determining the planetary mass has only allowed to establish an upper limit as defined by gravitational perturbation thresholds.

From Table 1, we can infer that most of the Kepler CBPs are in the Neptune to Jupiter range, with expected tidal values ${Q^{\prime }}_2={10}^3-{10}^5$ (Ferraz-Mello 2013; Lainey 2016). For such planets, a CTL-model as the one adopted here should be a good approximation. On the other hand, there are some planets such as Kepler-47b and Kepler-47c that are probably closer to the super-earth range. In such systems, the direct application of our model may not be adequate (Efroimsky 2012, 2015).

Although we are considering stars with very different masses and different possible internal structure, in the following sections we will assume the same value of the tidal parameters: $Q_{0^{\prime }}=Q_{1^{\prime }}=1\times {10}^6$.

4.1 Tidal timescales

In order to evaluate the importance of tidal evolution in these systems, we define characteristic timescales for the planetary spin evolution (τ_s), the semimajor axis evolution (τ_a) and eccentricity evolution (τ_e) according to: $$\frac{1}{{\tau }_{\text{s}}}=\frac{1}{{\Omega }_2}<\frac{\text{d}{\Omega }_2}{\text{d}t}>;\frac{1}{{\tau }_{\text{a}}}=\frac{1}{a_2}<\frac{\text{d}{a}_2}{\text{d}t}>;\frac{1}{{\tau }_{\text{e}}}=\frac{1}{e_2}<\frac{\text{d}{e}_2}{\text{d}t}>.$$(22)

In the case of the planetary spin, the simple form of Eq. (A.1) allow us to accurately estimate the characteristic timescales of a CBP due to tidal effects as $${\tau }_{\text{s}}=\frac{2{\mathcal{C}}_2{a}_2^6}{3\mathcal{G}{\mathcal{R}}_2^5{k}_{2,2}\Delta t_2}{\left({\sum }_{i=0}^4{K}_i^{(s)}{C}_i^{(s)}{\alpha }^i\right)}^{-1}.$$(23)

Results are presented in Table 2, were the first numerical value is the current semimajor axis ratio α of the CBP. The two following columns give the characteristic timescale necessary for the planet to reach a stationary spin (i.e. τ_s) and the corresponding equilibrium value $<{\Omega }_2{>}^{(\text{st})}$, expressed in units of its orbital frequency n₂. The planetary moment of inertia was taken equal to ${\mathcal{C}}_2=0.25m_2{\mathcal{R}}_2^2$.

The values of τ_s were calculated adopting ${Q^{\prime }}_2=1$ so they should be scaled appropriately for other planetary tidal parameters. According to Eq. (23), ${\tau }_{\text{s}}\propto {Q^{\prime }}_2$, making it straightforward to relate both quantities. In agreement with the numerical experiments done for Kepler-38 by Zoppetti et al. (2019a), pseudo-synchronisation appears to be attained rapidly in most Kepler systems even for large value for ${Q^{\prime }}_2$. With the possible exception of Kepler-1647, we expect all Kepler CBP to currently lie in stationary spin-orbit configurations. As can be seen in the following column, all the stationary spins are expected to be sub-synchronous, some of them for a large amount (such as Kepler-35) and others almost in perfect synchronisation with its own mean motion (like Kepler-34).

The four last columns of Table 2 give the average time derivatives and characteristic timescales of the planetary semimajor axis and eccentricity. We assumed stationary values for Ω₂ and current values for a₂ and e₂. Due to the complex functional dependence of the derivative with the initial conditions, the numerical values shown in the table should be considered local and not necessarily indicative of their primordial magnitudes. Even so, they give a qualitative idea of the strength of the tidal effects and, particularly, the sign of the present-day orbital variation.

For the observed orbital and physical parameters, our model predicts that all Kepler circumbinary planets should be migrating outwards as a consequence of their tidal interactions with the binary. However, the timescale necessary for a significantorbital migration is much larger than the age of the host star, even adopting very small values of ${Q^{\prime }}_2$. Conversely, tidal evolution always seems to damp the planetary eccentricity, and the circularization process appears to occur in a slightly smaller timescale than the semimajor axis.

A final word on the case of Kepler-1647. Due to its large distance from the binary, as expressed by its small value of α, the characteristic timescales for tidal evolution are typically three orders of magnitude larger than in any other system. Even so however, it is still possible for the planet to be close to a spin-orbit stationary solution even for moderate-to-large values of ${Q^{\prime }}_2$.

4.2 A geometrical interpretation for outward migration

The prediction that all Kepler CBPs should be experiencing outward migration from the tidal effect raises the question as to why this happens. Although a full explanation is beyond the scope of this paper, we can present a simple geometrical interpretation that may give some insight.

We consider the case in which the binary and the planetary orbit are circular, e₁ = e₂ = 0. In addition, weassume that the planet is sufficiently far from the binary to neglect all terms in Eq. (16) explicitly dependent on α. In this approximation, the stationary spins are trivial and equal to the mean motion: $<{\Omega }_0^{(\text{st})}>=<{\Omega }_1^{(\text{st})}>=n_1$ and $<{\Omega }_2^{(\text{st})}>=n_2$. From expression in Appendix B, we find that the time derivative of the planetary semimajor axis acquires the simple form: $$\frac{1}{a_2}<\frac{\text{d}{a}_2}{\text{d}t}>=\frac{6n_2{m}_2}{m_0+m_1}[k_{2,0}\Delta t_0{(\frac{{\mathcal{R}}_0}{a_2})}^5+k_{2,1}\Delta t_1{(\frac{{\mathcal{R}}_1}{a_2})}^5](n_1-n_2).$$(24)

We can compare this expression with that obtained in the two-body problem for a point-mass planets of mass m₂ in a circular orbit around a star m₀ rotating with arbitrary spin Ω₀. In such a case, the semimajor axis of the planet evolves according to $$\frac{1}{a_2}\frac{\text{d}{a}_2}{\text{d}t}=\frac{6n_2{m}_2}{m_0}[k_{2,0}\Delta t_0{(\frac{{\mathcal{R}}_0}{a_2})}^5]({\Omega }_0-n_2),$$(25)

as can be seen, for example, in Hut (1980). The similarity between these two equations helps explain why CBPs with low-to-moderate eccentricities migrate outwards. In the circumbinary geometry, the tidal effect of the binary star system may be substituted by a single body with a rotational frequency equal to the binary mean motion n₁. Since this quantity is always greater than the mean motion of the planet n₂, the semimajor axis of the planet increases. This situation is analogous to the two-body problem in which the central star m₀ rotates faster than orbital motion of the planet.

4.3 Past tidal evolution

Our next step is to attempt to reconstruct the past tidal evolution of the Kepler CB systems, as well as estimate their future trends. We assume that all bodies (stars and planet) are in stationary spin-orbit configurations and analyse only the changes in semimajor axis and eccentricity. Since the mass or radius of some of the planets are not well constrained,we completed the necessary physical parameters with the empirical mass-radius relation defined by Mills & Mazeh (2017).

Figure 4 shows, for each Kepler CBP, the tidal evolution velocity field around the current location of the planet in the (α, e₂) plane. The velocity field was computed assuming stationary spins for all intervening bodies and disregardingthe tidal evolution of the binary stars. The arrows represent the direction of orbital evolution in α and e₂, and its size was kept constant.

The green curves show initial conditions with $<\text{d}{a}_2/\text{d}t>=0$. All points on the plane above display inward migration (da₂∕dt < 0) while those below the curve lead to outward tidal migration. Kepler-34 shows a second green curve for large values of α that is also observed for some of the other systems for even higher values of α (which lie outside of the range adopted in the panels of Fig. 4). However, above such a curve, our model predicts an increase in the eccentricity of the synchronised CBP. We suspect that this may be a spurious effect consequence of the truncation to fourth-order in the semimajor axis ratio when estimating the stationary CBP spin Eq. (13).

The magenta curves show the capture eccentricity e_2,cap as function of α for each system. Its value is given by the mean-square average between the forced eccentricity (Moriwaki & Nakagawa 2004) and the mean eccentricity calculated for zero-amplitude secular variations (Paardekooper et al. 2012). Explicitly, $$e_{\text{2},\text{cap}}=\sqrt{e_{2,f}^2+e_{2,P}^2},$$(26)

where the classical forced eccentricity is given by $$e_{2,f}=\frac{5}{4}\alpha \frac{(m_0-m_1)}{(m_0+m_1)}{e}_1\frac{(1+\frac{3}{4}{e}_1^2)}{(1+\frac{3}{2}{e}_1^2)},$$(27)

while the mean eccentricity takes the form $$e_{2,P}=\frac{3}{4}{\alpha }^2\frac{m_0{m}_1}{{(m_0+m_1)}^2}\sqrt{1+\frac{34}{3}{e}_1^2}.$$(28)

More details of the calculations leading to these values may be found in Zoppetti et al. (2019b). This so-called capture eccentricity defines the average value expected as the result of non-conservative exterior forces acting on the system.

Finally, initial conditions inside the pink region are dynamically unstable according to the empirical criterion of Holman & Wiegert (1999). Note however, that this criterion is only valid for small-mass planets and circular orbits, so their extent should be considered more qualitative than accurate.

Analysing the velocity vector fields, we first note that the eccentricity usually appears more affected than the semimajor axis, except for quasi-circular orbits where the opposite occurs. When the gravitational perturbations are includedin the model, it is possible to observe that the orbital evolution of the system follows the locus of e_2,cap as function of α, as observed in Zoppetti et al. (2018) for Kepler-38.

In some systems we also observe that the domain associated to outward migration can reach high eccentricities (e.g. Kepler-35, Kepler-413, and Kepler-1647) while in others it is restricted to more circular trajectories (e.g. Kepler-16, Kepler-38, and Kepler-64). An explanation for such dichotomy may be found in Zoppetti et al. (2019a), where we found that in the low-eccentric systems, the size of the outward migrating-CB-region is directly proportional to the reduced mass $\overline{\mu }$ and inversely proportional to e₂. It is straightforward to check from Table 1 that Kepler-35, Kepler-413 and Kepler-1647 satisfy quite these conditions while, for example, in the rest of the systems the value of $\overline{\mu }$ is significantly lower.

The proximity of the observed Kepler circumbinary planets to the stability limit with the binary has been the subject of several recent studies (e.g. Quarles et al. 2018). While it has been proposed that this pile-up is unlikely to be affected by observational bias (Li et al. 2016), the recent discovery of a third and outer planet in Kepler-47 (Orosz et al. 2019) may indicate that we are only detecting the tail of the distribution. Nevertheless, we still need to address the current population and its past orbital evolution.

A key aspect of this question is the tidal evolution of the binary system and, particularly, what were the primordial orbital separation and eccentricity of the binary at the time of the planet’s formation and migration. In Zoppetti et al. (2018) we showed that even for moderate values of the stellar tidal parameters the original binary eccentricity could have ben much larger, thereby pushing the instability barrier closer to the planet. Similar results were also found for Kepler-34. An outward tidal migration of the planet itself would have worked in the same direction leading to a potentially more unstable primordial configuration.

These results could be interpreted as evidence that the orginal location of the CBPs was closely tied to the instability barried, perhaps in the form of sub-dense inner gaps in the protoplanetary disks. If so, then planetary traps such as those proposed by Kley & Haghighipour (2014) and Thun & Kley (2018) would be a more probable stalling mechanism for planetary migration than resonance capture.

Additionally, it is interesting to note that most of the CBPs located tightly packed to the binaries have eccentricities also very close to the capture eccentricity e_2,cap, with the exceptions of Kepler-34 and Kepler-413. Recently, Thun & Kley (2018) suggested that low-mass CBPs are strongly influenced by the protoplanetary disc and their eccentricities may have been very excited during the migration process. If this hypothesis is confirmed, then another non-conservative effect must be invoked to explain their current state. Tidal interactions appears as a possibility, even if they would need to be extremely accentuated. Curiously, Kepler-35 and Kepler-38, two systems with e₂ ~ e_2,cap belong to very old systems with estimated ages of the order of ~10 Gyrs (Welsh et al. 2012; Orosz et al. 2012a). Perhaps accumulated tidal effects over such a long time could be at least partially responsible, and could help in constraining the magnitudes of tidal parameters for both stars and planets.

Fig. 4 Velocity vector fields in the (α, e2) plane depicting the routes of tidal evolution in the vicinity of each Kepler CBP. The arrows show the direction of orbital evolution throughout the plane, although its size was kept constant and is thus not representative of the magnitude of the derivatives. The green curves corresponds to $<\text{d}{a}_2/\text{d}t>=0$) while the magenta curves indicate the location of the capture eccentricity ecap as defined by Eq. (26). The pale pink region correspond to values of the orbital elements identified as unstable according to the criterion by Holman & Wiegert (1999). Current positions of the Kepler CBPs are shown in blue circles.

5 Summary and discussion

We presented an extended version of the analytical tidal model for CB systems introduced in Zoppetti et al. (2019a). Once again, we assumed that all the bodies are extended and viscous in such a way that all interactions between pairs are considered but under the weak-friction regime, in which the tidal forces can be approximated by the classical expressions of Mignard (1979).

Startingfrom the variational equations of the spin and orbital evolution of the planet, we constructed an analytical approach by averaging over the mean longitudes. In all the cases, we expanded up to fourth-order in the semimajor axes ratio α but obtained closed expressions in the binary and planetary eccentricities. The resulting analytical model was compared with the results of numerical simulations and a very good agreement was found for all eccentricities and semimajor axes ratios up to α ~0.35.

Having expressions closed in eccentricities allowed us to apply the model to all well known Kepler CBPs. We investigated their past tidal evolution in two steps. First, we calculated the characteristic tidal timescales and found that the typical time required for a CBP to acquire its stationary spin is typically much lower than the expected age of the host star, even for moderate-to-large values of Q′. Consequently, most of the observed system should lie in spin-orbital stationary solutions with a sub-synchronous rotational frequency. Regarding the orbital evolution we found that the typical tidal timescales in the Kepler systems are much longer than the rotational timescales, and little tidal induced orbital migration is expected to have occurred. The eccentricity damping timescales, however, are one or two orders of magnitude lower than τ_a and tidal interactions could have caused some decrease in the eccentricities during the systems’ age.

In a second part of the work, we obtained some insight on the past tidal histories studying the tidal velocity fields around each planet. We found that all bodies are located in a tidal stream associated to an outward migration, pushing the planet away from the binary. Furthermore, with the exception of Kepler-34, all other systems are located distant from the curve separating inward and outward migrations. This seems to indicate that most or all of the systems lifetime was spent inside the region of outward migration. Consequently, their primordial separation from the binary should have been smaller.

Thun & Kley (2018) recently suggested that low-mass CBPs may have suffered strong excitation from the protoplanetary disk, leading to a highly eccentric final orbit. This contrasts with the current values which, in most cases, are found very close to the capture eccentricity. Perhaps tidal evolution played some role in damping the primordial value, which would help define some constraints on the numerical values of the tidal parameters of the system.

We mention, finally, that most of the CBPs discovered so far by the Kepler mission are in the range of Neptune to Jupiter masses, where CTL-models for tides should be valid. However, we note that not much is known about the internal structure of this bodies and therefore a precise estimation of the ${Q^{\prime }}_2$-values becomes a difficult task.

Acknowledgements

The authors are very grateful to M. Efroimsky for the helpful comments and suggestions. This research was funded by CONICET, SECYT/UNC, and FONCYT.

Appendix A Coefficients for the spin equations

The variational equation for the planetary spin, up to fourth order in α, was given in Eq. (10), which is repeated here for convenience: $$<\frac{\text{d}{\Omega }_2}{\text{d}t}>=\frac{3\mathcal{G}{\mathcal{R}}_2^5{k}_{2,2}\Delta t_2}{2{\mathcal{C}}_2{a}_2^6}{\sum }_{i=0}^4{K}_i^{(s)}(A_i^{(s)}{n}_1+B_i^{(s)}{n}_2+C_i^{(s)}{\Omega }_2){\alpha }^i.$$

The $A_i^{(s)}$ coefficients, multiplying the mean motion n₁ of the binary, acquire the form $$\begin{array}{lll}{A}_0^{(s)}\hfill & =\hfill & 0,\hfill \\ A_1^{(s)}\hfill & =\hfill & 0,\hfill \\ A_2^{(s)}\hfill & =\hfill & -\frac{6\sqrt{1-e_1^2}}{{(1-e_2^2)}^{13/2}}{X}_1,\hfill \\ A_3^{(s)}\hfill & =\hfill & -\frac{168\sqrt{1-e_1^2}}{{(1-e_2^2)}^{15/2}}{X}_3{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ A_4^{(s)}\hfill & =\hfill & -\frac{6\sqrt{1-e_1^2}}{{(1-e_2^2)}^{17/2}}(8Y_1{X}_5+175X_6{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \end{array}$$(A.1)

in terms of eccentricity functions Y_i(e₁) and X_i(e₂) which are explicitly given in Appendix C. The coefficients that multiply the mean motion of the planet are found to be $$\begin{array}{lll}{B}_0^{(s)}\hfill & =\hfill & \frac{2}{{(1-e_2^2)}^6}{X}_1,\hfill \\ B_1^{(s)}\hfill & =\hfill & \frac{72}{{(1-e_2^2)}^7}{X}_3{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ B_2^{(s)}\hfill & =\hfill & \frac{3}{{(1-e_2^2)}^8}(8Y_1{X}_5+175X_6{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ B_3^{(s)}\hfill & =\hfill & \frac{50e_1{e}_2}{{(1-e_2^2)}^9}(24Y_2{X}_8\cos (\Delta \varpi )+49X_9{e}_1^2{e}_2^2\cos (3\Delta \varpi )),\hfill \\ B_4^{(s)}\hfill & =\hfill & \frac{15}{{(1-e_2^2)}^{10}}(8Y_3{X}_{11}+756Y_4{X}_{12}{e}_1^2{e}_2^2\cos (2\Delta \varpi ),\hfill \\ \hfill & \hfill & +\frac{9261}{16}{X}_{13}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(A.2)

Finally, the coefficients accompanying the planetary spin rate are $$\begin{array}{lll}{C}_0^{(s)}\hfill & =\hfill & -\frac{2}{{(1-e_2^2)}^{9/2}}{X}_2,\hfill \\ C_1^{(s)}\hfill & =\hfill & -\frac{45}{{(1-e_2^2)}^{11/2}}{X}_4{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ C_2^{(s)}\hfill & =\hfill & -\frac{9}{{(1-e_2^2)}^{13/2}}(2Y_1{X}_1+25X_7{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ C_3^{(s)}\hfill & =\hfill & -\frac{35e_1{e}_2}{{(1-e_2^2)}^{15/2}}(18Y_2{X}_3\cos (\Delta \varpi )+\frac{175}{8}{X}_{10}{e}_1^2{e}_2^2\cos (3\Delta \varpi )),\hfill \\ C_4^{(s)}\hfill & =\hfill & -\frac{9}{{(1-e_2^2)}^{17/2}}(8Y_3{X}_5+490Y_4{X}_6{e}_1^2{e}_2^2\cos (2\Delta \varpi ),\hfill \\ \hfill & \hfill & +\frac{3675}{16}{X}_{14}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(A.3)

We note that all coefficients are function of the eccentricities, while some are also dependent on the difference in the longitudes of pericenter.

Recalling Eq. (13), the stationary spin rate predicted by our model may be written as $$<{\Omega }_2{>}^{(\text{st})}={\sum }_{i=0}^4<{\Omega }_2{>}_i^{(\text{st})}{\alpha }^i.$$

The different order terms $<{\Omega }_2{>}_i^{(\text{st})}$ are given by $$<{\Omega }_2{>}_0^{(\text{st})}=-\frac{1}{C_0^{(s)}}{D}_0^{(s)},<{\Omega }_2{>}_1^{(\text{st})}=-\frac{K_1^{(s)}}{K_0^{(s)}{C}_0^{(s)}}(D_1^{(s)}+C_1^{(s)}<{\Omega }_2{>}_0^{(\text{st})}),<{\Omega }_2{>}_2^{(\text{st})}=-\frac{K_2^{(s)}}{K_0^{(s)}{C}_0^{(s)}}\text{(}{D}_2^{(s)}+C_2^{(s)}<{\Omega }_2{>}_0^{(\text{st})}+\frac{K_1^{(s)}}{K_2^{(s)}}{C}_2^{(s)}<{\Omega }_2{>}_1^{(\text{st})}\text{)},<{\Omega }_2{>}_3^{(\text{st})}=-\frac{K_3^{(s)}}{K_0^{(s)}{C}_0^{(s)}}\text{(}{D}_3^{(s)}+C_3^{(s)}<{\Omega }_2{>}_0^{(\text{st})}+\frac{K_2^{(s)}}{K_3^{(s)}}{C}_2^{(s)}<{\Omega }_2{>}_1^{(\text{st})}\text{.}+\frac{K_1^{(s)}}{K_3^{(s)}}{C}_1^{(s)}<{\Omega }_2{>}_2^{(\text{st})}\text{)}<{\Omega }_2{>}_4^{(\text{st})}=-\frac{K_4^{(s)}}{K_0^{(s)}{C}_0^{(s)}}\text{(}{D}_4^{(s)}+C_4^{(s)}<{\Omega }_2{>}_0^{(\text{st})}+\frac{K_3^{(s)}}{K_4^{(s)}}{C}_3^{(s)}<{\Omega }_2{>}_1^{(\text{st})}\text{.}+\frac{K_2^{(s)}}{K_4^{(s)}}{C}_2^{(s)}<{\Omega }_2{>}_2^{(\text{st})}+\frac{K_1^{(s)}}{K_4^{(s)}}{C}_1^{(s)}<{\Omega }_2{>}_3^{(\text{st})}\text{)},$$(A.4)

where, for the sake of simplicity, we have defined new auxiliary functions $$D_i^{(s)}=A_i^{(s)}{n}_1+B_i^{(s)}{n}_2.$$

Appendix B Coefficients for the tidal evolution of the semimajor axis and eccentricity

As seen in Eq. (16), the time derivative of the planetary semimajor axis may be writen as $$<\frac{\text{d}{a}_2}{\text{d}t}>=\frac{n_2}{\mathcal{G}{m}_2{\sigma }_1{a}_2^4}{\sum }_{i=0}^4{K}_i^{(a)}(A_i^{(a)}{n}_1+B_i^{(a)}{n}_2+C_i^{(a)}{\overline{\Omega }}_i^{*}){\alpha }^i,$$

where the definitions of $K_i^{(a)}$ and ${\overline{\Omega }}_i^{*}$ were explicitly given in Eqs. (17) and (18). The coefficients are: $$\begin{array}{lll}{A}_0^{(a)}\hfill & =\hfill & 0,\hfill \\ A_1^{(a)}\hfill & =\hfill & 0,\hfill \\ A_2^{(a)}\hfill & =\hfill & \frac{10\sqrt{1-e_1^2}}{{(1-e_2^2)}^8}{X}_5,\hfill \\ A_3^{(a)}\hfill & =\hfill & \frac{400\sqrt{1-e_1^2}}{{(1-e_2^2)}^9}{X}_8{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ A_4^{(a)}\hfill & =\hfill & \frac{25\sqrt{1-e_1^2}}{8{(1-e_2^2)}^{10}}(32Y_1{X}_{11}+9\left(120X_{12}-184e_1^2{X}_{21}+127e_1^4{X}_{22}\right)\hfill \\ \hfill & \hfill & \cdot \frac{e_1^2{e}_2^2\cos (2\Delta \varpi )}{{(1-e_1^2)}^2}).\hfill \end{array}$$(B.1) $$\begin{array}{lll}{B}_0^{(a)}\hfill & =\hfill & -\frac{2}{{(1-e_2^2)}^{15/2}}{X}_{15},\hfill \\ B_1^{(a)}\hfill & =\hfill & -\frac{114}{{(1-e_2^2)}^{17/2}}{X}_{16}{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ B_2^{(a)}\hfill & =\hfill & -\frac{1}{{(1-e_2^2)}^{19/2}}(34Y_1{X}_{17}+1195X_{18}{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ B_3^{(a)}\hfill & =\hfill & -\frac{25}{8{(1-e_2^2)}^{21/2}}(772Y_2{X}_{19}{e}_1{e}_2\cos (\Delta \varpi )\hfill \\ \hfill & \hfill & +2415X_{20}{e}_1^3{e}_2^3\cos (3\Delta \varpi ),\hfill \\ B_4^{(a)}\hfill & \text{}=\text{}\hfill & -\frac{5}{32{(1-e_2^2)}^{23/2}}(\text{}1408Y_3{X}_{23}+195300Y_4{X}_{24}{e}_1^2{e}_2^2\cos (2\Delta \varpi )\hfill \\ \hfill & \hfill & +222075X_{25}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(B.2) $$\begin{array}{lll}{C}_0^{(a)}\hfill & =\hfill & \frac{2}{{(1-e_2^2)}^6}{X}_1,\hfill \\ C_1^{(a)}\hfill & =\hfill & \frac{72}{{(1-e_2^2)}^7}{X}_3{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ C_2^{(a)}\hfill & =\hfill & \frac{3}{{(1-e_2^2)}^8}(8Y_1{X}_5+175X_7{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ C_3^{(a)}\hfill & =\hfill & \frac{50}{{(1-e_2^2)}^9}(24Y_2{X}_8{e}_1{e}_2\cos (\Delta \varpi )+49X_9{e}_1^3{e}_2^3\cos (3\Delta \varpi )),\hfill \\ C_4^{(a)}\hfill & =\hfill & \frac{5}{16{(1-e_2^2)}^{10}}(384Y_3{X}_{11}+36288Y_4{X}_{12}{e}_1^2{e}_2^2\cos (2\Delta \varpi )\hfill \\ \hfill & \hfill & +27783X_{13}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(B.3)

In the same manner, the planetary eccentricity evolution up to fourth order in α is given by $$<\frac{\text{d}{e}_2^2}{\text{d}t}>=\frac{n_2}{\mathcal{G}{m}_2{\sigma }_1{a}_2^5}{\sum }_{i=0}^4{K}_i^{(e)}(A_i^{(e)}{n}_1+B_i^{(e)}{n}_2+C_i^{(e)}{\overline{\Omega }}_i^{*}){\alpha }^i$$

where $K_i^{(e)}=K_i^{(a)}$ and the coefficients are listed below: $$\begin{array}{lll}{A}_0^{(e)}\hfill & =\hfill & 0,\hfill \\ A_1^{(e)}\hfill & =\hfill & 0,\hfill \\ A_2^{(e)}\hfill & =\hfill & \frac{75\sqrt{1-e_1^2}}{{(1-e_2^2)}^7}{X}_3{e}_2^2,\hfill \\ A_3^{(e)}\hfill & =\hfill & \frac{85\sqrt{1-e_1^2}}{{(1-e_2^2)}^8}{X}_{30}{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ A_4^{(e)}\hfill & =\hfill & \frac{25\sqrt{1-e_1^2}}{8{\left(1-e_2^2\right)}^9}[304Y_1{X}_8{e}_2^2\hfill \\ \hfill & \hfill & +\left(380X_{34}-256e_1^2{X}_{35}+443e_1^4{X}_{36}\right)\frac{e_1^2{e}_2^2\cos \left(2\mathrm{\Delta }\varpi \right)}{\left(1-e_1^2\right)2}].\hfill \end{array}$$(B.4) $$\begin{array}{lll}{B}_0^{(e)}\hfill & =\hfill & -\frac{18}{{(1-e_2^2)}^{13/2}}{X}_3{e}_2^2,\hfill \\ B_1^{(e)}\hfill & =\hfill & -\frac{27}{{(1-e_2^2)}^{15/2}}{X}_{26}{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ B_2^{(e)}\hfill & =\hfill & -\frac{5}{4{(1-e_2^2)}^{17/2}}(296Y_1{X}_8{e}_2^2+361X_{28}{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ B_3^{(e)}\hfill & =\hfill & -\frac{25}{8{(1-e_2^2)}^{19/2}}(148Y_2{X}_{31}{e}_1{e}_2\cos (\Delta \varpi )\hfill \\ \hfill & \hfill & +1141X_{32}{e}_1^3{e}_2^3\cos (3\Delta \varpi )),\hfill \\ B_4^{(e)}\hfill & =\hfill & -\frac{15}{8{(1-e_2^2)}^{21/2}}(1504Y_3{X}_{37}{e}_2^2\hfill \\ \hfill & \hfill & +5187Y_4{X}_{38}{e}_1^2{e}_2^2\cos (2\Delta \varpi )+10017X_{39}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(B.5) $$\begin{array}{lll}{C}_0^{(e)}\hfill & =\hfill & \frac{11}{{(1-e_2^2)}^5}{X}_4{e}_2^2,\hfill \\ C_1^{(e)}\hfill & =\hfill & \frac{39}{2{(1-e_2^2)}^6}{X}_{27}{e}_1{e}_2\cos (\Delta \varpi ),\hfill \\ C_2^{(e)}\hfill & =\hfill & \frac{45}{{(1-e_2^2)}^7}(4Y_1{X}_3{e}_2^2+5X_{29}{e}_1^2{e}_2^2\cos (2\Delta \varpi )),\hfill \\ C_3^{(e)}\hfill & =\hfill & \frac{85}{16{(1-e_2^2)}^8}(48Y_2{X}_{30}{e}_1{e}_2\cos (\Delta \varpi )\hfill \\ \hfill & \hfill & +145X_{33}{e}_1^3{e}_2^3\cos (3\Delta \varpi )),\hfill \\ C_4^{(e)}\hfill & =\hfill & \frac{285}{8{(1-e_2^2)}^9}(\frac{2224}{57}{Y}_3{X}_8{e}_2^2\hfill \\ \hfill & \hfill & +112Y_4{X}_{34}{e}_1^2{e}_2^2\cos (2\Delta \varpi )+147X_{40}{e}_1^4{e}_2^4\cos (4\Delta \varpi )).\hfill \end{array}$$(B.6)

Appendix C Eccentricity functions

The eccentricity functions Y_i(e₁) and X_i(e₂) to which we refer in the previous equations are found to be combinations of the Hansen coefficients. Explictly, the first are found to be: $$\begin{array}{lll}{Y}_1\hfill & =\hfill & 1+\frac{3}{2}{e}_1^2;Y_2=1+\frac{3}{4}{e}_1^2\hfill \\ Y_3\hfill & =\hfill & 1+5e_1^2+\frac{15}{8}{e}_1^4;Y_4=1+\frac{1}{2}{e}_1^2,\hfill \end{array}$$(C.1)

while the functions of the planetary eccentricity are given in Table C.1:

References

All Tables

All Figures

Fig. 1 Time derivative of the planetary spin rate Ω2 in our generic binary system as a function of α = a1∕a2. Each panelassumes different eccentricities for both orbits. Broad gray curves show results from semi-analytical calculationsand predictions from analytical models are depicted in thin lines. Values obtained with our original (Zoppetti et al. 2019a) are presented in blue, while the colors correspond to our current model with different orders in α. See inlaid box in left-hand plot for details.

In the text

	Fig. 2 Stationary spin of a CB planet as a function of α. Different rows consider different eccentricity of the orbits (assumed equals). Different columns represent different reduced mass: $\overline{\mu }=0.1$ (left panels) and $\overline{\mu }=0.5$ (right panels). The type of the curves represent the same that in Fig. 1.
In the text

	Fig. 3 Time derivative of the planetary semimajor axis (left panels) and eccentricity (right panels), as a function of α, for a synchronised CB planet around our generic binary system. As before, results obtained with different models are shown in different colored curves: numerical in broad grey, Zoppetti et al. (2019a) in thin blue, and predictions from our current α4 model in black.
In the text

Fig. 4 Velocity vector fields in the (α, e2) plane depicting the routes of tidal evolution in the vicinity of each Kepler CBP. The arrows show the direction of orbital evolution throughout the plane, although its size was kept constant and is thus not representative of the magnitude of the derivatives. The green curves corresponds to $<\text{d}{a}_2/\text{d}t>=0$) while the magenta curves indicate the location of the capture eccentricity ecap as defined by Eq. (26). The pale pink region correspond to values of the orbital elements identified as unstable according to the criterion by Holman & Wiegert (1999). Current positions of the Kepler CBPs are shown in blue circles.

In the text