© ESO 2019

1 Introduction

As of 2019, the Kepler mission has discovered approximately ten circumbinary (CB) planetary systems. All binary components define compact systems with orbital periods less than ~40 days and a wide range of eccentricities and mass ratios. The planets surrounding them also have a diversity of masses (between super-Earths to Jupiter masses) but they are all almost coplanar with the binary. With the exception of Kepler 34b and Kepler 413b, all CB planets seem characterized by small semimajor axis low eccentricities.

While the low inclinations suggest that these planets formed in a CB disc aligned with the orbital plane of the central binary, it is well accepted that in situ formation so close to the binary is unlikely due to the strong eccentricity excitation induced by the secondary star (e.g. Lines et al. 2014; Meschiari 2012). However, as we move away from the binary, the gravitational potential approaches that of a single star and planetary formation appears to be easier, following usual core-accretion models. This suggest that CB planets could have formed farther out, later migrated inwards due to interaction with a primordial disc and finally stalled near their current orbits by some mechanism (Dunhill & Alexander 2013).

In a previous work (Zoppetti et al. 2018), we tested the possibility that the circumbinary planets may have halted its inward migration due to a capture in a high-order mean-motion resonance (MMR) with the binary and, once the disc is dissipated, slowly escaped from the commensurability due to tidal forces. We applied this hypothesis to Kepler-38, a very old system in which captures in the 5/1 MMR had been reported with hydro-simulation (Kley & Haghighipour 2014). Tidal interactions were modelled following Mignard (1979) and incorporated to a N-body integrator following the prescription detailed in Rodriguez et al. (2011). We observed that while the binary orbit shrinks due to tidal interactions, the planet seemed to increase its semimajor axis, even after the system reached stationary solutions in the spin rates. We were unable to explain these findings, which in principle could have been caused by a non-consistent treatment of the tidal interactions between the different bodies of the dynamical system.

In this article we present and discuss a self-consistent tidal model for a multi-body system, in which all tidal forces between pairs are computed adopting a weak-friction (Mignard-type) model. While the model is general, we will focus primarily on the spin and orbital evolution of the CB planet. To allow for a simpler comparison with our previous results, we will once again employ Kepler-38 as a reference system (Orosz et al. 2012). However, we will also explore a wider range of system parameters as well as different initial orbital elements and spin rates.

This paper is organized as follows. In Sect. 2 we present the model in two steps: in Sect. 2.1 we first discuss which tidal forces have a net effect onto the dynamical evolution of a three-extended-body system, while in Sect. 2.2 we show how these forces are incorporated self-consistently into our tidal model. Section 3 presents a series of numerical integrations of the full spin and orbital equations of motion. We concentrated on two different time-scales: the early dynamical evolution of the system before the spins reached stationary solutions, and the subsequent long-term orbital evolution of the CB planet in spin stationarity. In Sect. 4, we constructed analytical expressions for the orbital and spin evolution, averaged over the orbital periods but retaining secular terms, including those containing the difference between longitudes of pericenter. These allow us to estimate the stationary spin rate of CB planets, as well as the direction and magnitude of the orbital migration. We compared these predictions with full N-body simulations. Finally, Sect. 5 summarises our main results and discusses their implications.

2 The model

We first considered a binary system in which m₀ and m₁ are the masses of the stellar components and m₂ is a circumbinary planet. We supposed that all the bodies lie in the same orbital plane and their spins are perpendicular to it. We also assumed that all the bodies are extended masses with physical radii ${\mathcal{R}}_i$ and are deformable due to tidal effects between them.

For the gravitational interactions between each pair of bodies we adopted the classical weak-friction tidal model (Mignard 1979). However, since the tidal deformation of each body had to incorporate the gravitational potential generated by both of its companions, we first needed to address two issues: (i) identifying which tidal deformations have a net effect on the long-term dynamical evolution of the system and, (ii) understanding how the different forces should be incorporated into a self-consistent physical model. These issues are addressed in the next two subsections.

2.1 The Mignard forces revisited

We began by considering our three-body system as having two simplifications. First, we neglected the gravitational perturbations generated by m₂ on the other two bodies, as well as the effects of m₁ on m₂. Second, only m₀ was assumed to be an extended mass while m₁ and m₂ were taken as point masses. As a consequence of these approximations, the dynamics of both m₁ and m₂ around m₀ were defined by the point-mass approximation plus the tidal deformation of m₀ generated solely by m₁. The role of m₂ was thus reduced to serve as a tracker of the dynamical effect of the tidal bulge on any generic orbit in the configuration plane.

A schematics of this scenario is presented in Fig. 1, where r_i are the m₀ -centric position vectors of the other masses. Following Mignard (1979), the tidal bulge of m₀ considered is displaced with respect to the instantaneous position of m₁ by a constant time-lag Δt₀. We assume the lag is sufficiently small to expand the gravitational potential U generated by m₀ in anywhere in the space up to first-order in Δt₀, such that $$U(r,r_1)=U^{(0)}(r,r_1)+U^{(1)}(r,r_1)+\mathcal{O}(\Delta t_0^2),$$(1)

where U⁽⁰⁾ and U⁽¹⁾ are the expanded tidal potentials oforder $\mathcal{O}(0)$ and $\mathcal{O}(\Delta t_0)$, respectively.In particular, if we evaluate Eq. (1) on the position of m₂ (i.e. r = r₂), we obtain

Fig. 1 Tidal-lagged bulge generated on m0 due to m1 and its effecton test body m2, from m0 -centric coordinate frame.

$$\begin{array}{lll}{U}^{(0)}(r_2,r_1)\hfill & =\hfill & \frac{\mathcal{G}{m}_1{\mathcal{R}}_0^5}{2r_1^5{r}_2^5}{k}_{2,0}\left[3{(r_2\cdot r_1)}^2-r_2^2{r}_1^2\right],\hfill \\ U^{(1)}(r_2,r_1)\hfill & =\hfill & \frac{3\mathcal{G}{m}_1{\mathcal{R}}_0^5}{r_1^5{r}_2^5}{k}_{2,0}\Delta t_0[\frac{(r_1\cdot \stackrel{.}{r_1})}{2r_1^2}\left[5{(r_2\cdot r_1)}^2-r_2^2{r}_1^2\right]\hfill \\ \hfill & \hfill & -(r_2\cdot r_1)[r_1\cdot ({\Omega }_0\times r_2)+r_2\cdot \stackrel{.}{r_1}]],\hfill \end{array}$$(2)

where $\mathcal{G}$ is the gravitational constant, Ω₀ is the spin vector of m₀ and k_2,0 is the second degree Love number.

The tidal force per unit mass $f$ generated by m₀ at a generic position vector r can be obtained as $$f={\nabla }_{\text{r}}(U^{(0)}+U^{(1)})=f^{(0)}+f^{(1)},$$(3)

where explicit expressions evaluated on m₂ are given by $$\begin{array}{lll}{f}^{(0)}\hfill & =\hfill & \frac{3\mathcal{G}{m}_1{\mathcal{R}}_0^5}{2r_2^5{r}_1^5}{k}_{2,0}[2({\text{r}}_2\cdot {\text{r}}_1){\text{r}}_1+(r_1{}^2-\frac{5}{r_2^2}{({\text{r}}_2\cdot {\text{r}}_1)}^2){\text{r}}_2],\hfill \\ f^{(1)}\hfill & =\hfill & \frac{3\mathcal{G}{m}_1{\mathcal{R}}_0^5}{r_2^5{r}_1^5}{k}_{2,0}\Delta t_0[\frac{({\text{r}}_1\cdot \stackrel{.}{{\text{r}}_1})}{r_1^2}[5{\text{r}}_1({\text{r}}_2\cdot {\text{r}}_1)-{\text{r}}_2{r}_1{}^2]\hfill \\ \hfill & \hfill & -[{\text{r}}_1\cdot ({\Omega }_0\times {\text{r}}_2)+{\text{r}}_2\cdot \stackrel{.}{{\text{r}}_1}]{\text{r}}_1-({\text{r}}_1\times {\Omega }_0+\stackrel{.}{{\text{r}}_1})({\text{r}}_2\cdot {\text{r}}_1)\hfill \\ \hfill & \hfill & +\frac{5{\text{r}}_2}{r_2^2}[({\text{r}}_2\cdot {\text{r}}_1)[{\text{r}}_1\cdot ({\Omega }_0\times {\text{r}}_2)+{\text{r}}_2\cdot \stackrel{.}{{\text{r}}_1}]\hfill \\ \hfill & \hfill & -\frac{({\text{r}}_1\cdot \stackrel{.}{{\text{r}}_1})}{2r_1^2}\left[5{({\text{r}}_2\cdot {\text{r}}_1)}^2-r_2^2{r}_1{}^2\right]]].\hfill \end{array}$$(4)

Finally, the torques per unit mass can be calculated as T(r, r₁) ≃ r × (f⁽⁰⁾ + f⁽¹⁾) = T⁽⁰⁾(r, r₁) + T⁽¹⁾(r, r₁). As before, evaluating on the position of m₂ yields $$\begin{array}{lll}{T}^{(0)}({\text{r}}_2,{\text{r}}_1)\hfill & =\hfill & \frac{3\mathcal{G}{m}_1{k}_{2,0}{\mathcal{R}}_0^5}{r_2^5{r}_1^5}({\text{r}}_2\cdot {\text{r}}_1)({\text{r}}_2\times {\text{r}}_1)\hfill \\ T^{(1)}({\text{r}}_2,{\text{r}}_1)\hfill & =\hfill & \frac{3\mathcal{G}{m}_1{k}_{2,0}{\mathcal{R}}_0^5}{r_2^5{r}_1^5}\Delta t_0[5\frac{({\text{r}}_1\cdot \stackrel{.}{{\text{r}}_1})}{r_1^2}({\text{r}}_2\cdot {\text{r}}_1)({\text{r}}_2\times {\text{r}}_1)\hfill \\ \hfill & \hfill & -[{\text{r}}_1\cdot ({\Omega }_0\times {\text{r}}_2)+{\text{r}}_2\cdot \stackrel{.}{{\text{r}}_1}]({\text{r}}_2\times {\text{r}}_1)\hfill \\ \hfill & \hfill & -({\text{r}}_2\cdot {\text{r}}_1)[({\text{r}}_2\cdot {\Omega }_0){\text{r}}_1-({\text{r}}_2\cdot {\text{r}}_1){\Omega }_0+{\text{r}}_2\times \stackrel{.}{{\text{r}}_1}]].\hfill \end{array}$$(5)

In the classical two-body tidal problem, the acceleration f and the torque T are computed on the position of the deforming body m₁. It is easy to see that in such a case, the zero-order torque reduces to zero (i.e. T⁽⁰⁾ (r = r₁, r₁) = 0) and the only net contribution to the orbital and spin evolution (notwithstanding a precession term) stems from the first-order expressions f⁽¹⁾ and T⁽¹⁾ (see Eqs. (5) and (6) of Mignard 1979). However, it is not immediately clear what occurs if r ≠ r₁. In other words, we wish to analyse what are the (long-term) dynamical effects of a tidal bulge generated on m₀, due to the perturbing potential of m₁, on the orbit of another body m₂.

To address this question, let the m₀-centric orbits of m_i be characterized by semimajor axis a_i, eccentricity e_i, mean longitude λ_i and longitude of pericenter ϖ_i. Furthermore, let n_i denote the mean-motion (orbital frequency) of each body. We then compute the net secular torques (⟨T⁽⁰⁾(r₂, r₁)⟩ and ⟨T⁽¹⁾ (r₂, r₁)⟩) for different values of a₂, assuming fixed values for (a₁, e₁, e₂, ϖ₁, ϖ₂). The secular torques are calculated averaging over the short-period terms associated with λ₁ and λ₂. Since we will not restrict our analysis to non-resonant configurations between m₁ and m₂, we cannot assume that both mean longitudes are necessarily independent. We thus substitute the classical double averaging over λ_i with a time average over time, such that $$\langle T^{(i)}({\text{r}}_2,{\text{r}}_1)\rangle =\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\displaystyle {\int }_0^{\tau }{T}^{(i)}}({\text{r}}_2(t),{\text{r}}_1(t))\text{\hspace{0.17em}d}t,$$(6)

with i = 0, 1. This technique allows us to evaluate the net secular contribution in both resonant and secular configurations of both bodies. In particular, the classical Mignard expressions should be obtained assuming equal orbits (and orbital frequencies) for m₁ and m₂.

Results are shown in Fig. 2 for m₀ = 1, a₁ = 1, and ϖ₁ = ϖ₂ and e₁ = e₂. The values of the averaged torques are normalized with respect to m₁ and Δt₀, and plotted as function of the mean-motion ratio n₁ ∕ n₂. The left-hand panels correspond to the zero-oder torque ⟨T⁽⁰⁾⟩, while the first-order contributions ⟨T⁽¹⁾⟩ are depicted in the right-hand graphs.

In the upper panels we analyse the circular case (e₁ = e₂ = 0) and in the lower panel eccentric orbits (e₁ = e₂ = 0.1). Eccentricities are assumed to be fixed throughout the numerical averaging. All plots show distinct peaks, where the net torque is different from zero, overlaid with respect to background values that decrease smoothly as n₁ ∕ n₂ →∞. This background trend is a consequence of the numerical approximation employed to evaluate the time integral (Eq. (6)), which basically consisted of a discrete sum over a finite time interval equal to 500 orbital periods of the outermost body.

For circular orbits we observe that the only value of a₂ for which m₂ receives a non-zero net torque corresponds to a 1 ∕ 1 MMR, that is in a coorbital position with the deforming body m₁. In particular, the case in which the position of both masses coincide (i.e. λ₁ = λ₂) yields results analogous to those obtained from the two-body tidal model. Different initial values of the mean longitudes would, in principle, allow us to estimate both torques in other coorbital configurations, such as that occurring for m₂ located in a Trojan-like orbit with the other masses. Finally, although |⟨T⁰(r₂)⟩| is different from zero in the 1 ∕ 1 MMR, its dynamical effect on the orbit reduces to a tidal precession term (e.g. Correia et al. 2011) and does not contribute to any secular changes in the orbits or spin rates.

The eccentric case, depicted in the lower frames, exhibits a richer diversity. Non-zero torques are found in several MMRs and not only in the coorbital region. This seems to imply that the tidal deformation generated by m₁ on m₀ should affect the dynamical evolution of m₂ whenever there exists a commensurability relation between the orbital frequencies. This is an important finding, indicating that tidal models for resonant bodies could require the full tidal deformation on each body as generated by all the other bodies of the system.

The numerical results described in Fig. 2 were confirmed introducing elliptic expansions for the position and velocity vectors in the Eqs. (4) and (5), truncated up to fourth order in semimajor axis ratio and eccentricities, and integrating the resulting expressions analytically.

In conclusion, in the absence of any mean-motion relation between m₁ and m₂, the only tidal forces that need to be considered on m_i are those stemming from the deformation that m_i generates onm_j and m_k (with i≠ j≠k). Since the tidal deformation generated on m_j by m_k may be neglected, a multi-body tidal model may be constructed simply by adding the forces and torques between deformed-deforming pairs as given in Eqs. (5) and (6) of Mignard (1979). The effect of the tidal torques on resonant orbits will be investigated in a forthcoming work.

Fig. 2 Secular normalized torques of zero-order |⟨T0(r2)⟩| (left column) and first-order |⟨T1(r2)⟩| (right column), computed on m2 due to the tidal deformation on m0 induced by m1, plotted as a function of the mean-motion ratio n1 ∕ n2. We considered m0 = 1, a1 = 1 and varied a2 to include orbits both interior and exterior to m1. Upper panels: circular orbits (e1 = e2 = 0). Lower panels: eccentric orbits with e1 = e2 = 0.1. Light brown vertical lines highlight the location of some important MMRs. We note that the first-order torques in the right panels are also normalized respect to the time-lag Δt0.

2.2 The equations of motion

Having identified the tidal forces affecting the long-term and secular dynamical evolution of the system, we now discuss how they should be incorporated into the equations of motion of the circumbinary system in a self-consistent manner.

We return to our full circumbinary system where now all bodies are considered extended and gravitationally interacting. As shown in Fig. 3, the equilibrium deformation of body m_i is the sum of two ellipsoids, each generated by the gravitational potential of the other two masses. As shown in Folonier & Ferraz-Mello (2017), the sum of two ellipsoidal bulges can be approximated by a single ellipsoidal bulge with its own flattening and orientation. However, as discussed in Sect. 2.1, we only need to consider the direct distortion between pairs.

Let us denote by R_i the position vector of m_i in an inertial reference frame. Then, the complete equations of motion, including both tidal and point-mass terms, may be expressed as: $$\begin{array}{lll}{m}_0\stackrel{\mathrm{..}}{{\text{R}}_0}\hfill & =\hfill & \frac{\mathcal{G}{m}_0{m}_1}{|{\Delta }_{10}{|}^3}{\Delta }_{10}+\frac{\mathcal{G}{m}_0{m}_2}{|{\Delta }_{20}{|}^3}{\Delta }_{20}+F_0,\hfill \\ m_1\stackrel{\mathrm{..}}{{\text{R}}_1}\hfill & =\hfill & -\frac{\mathcal{G}{m}_0{m}_1}{|{\Delta }_{10}{|}^3}{\Delta }_{10}+\frac{\mathcal{G}{m}_1{m}_2}{|{\Delta }_{21}{|}^3}{\Delta }_{21}+F_1,\hfill \\ m_2\stackrel{\mathrm{..}}{{\text{R}}_2}\hfill & =\hfill & -\frac{\mathcal{G}{m}_0{m}_2}{|{\Delta }_{20}{|}^3}{\Delta }_{20}-\frac{\mathcal{G}{m}_1{m}_2}{|{\Delta }_{21}{|}^3}{\Delta }_{21}+F_2,\hfill \end{array}$$(7)

where for compactness we have denoted the relative position vectors as $${\Delta }_{\text{ij}}\equiv {\text{R}}_{\text{i}}-{\text{R}}_{\text{j}}.$$(8)

In terms of the m₀-centric position vectors, these are simply given by Δ₁₀ = r₁, Δ₂₀ = r₂ and Δ₂₁ = r₂ −r₁. The last terms of the equations of motion are the complete tidal forces acting on each mass. Following Ferraz-Mello et al. (2008), considering the reacting forces, these may be expressed by

Fig. 3 Tidally interacting three-body system. The tidal bulge generated on each body (filled grey ellipsoids) are the sum of the deformations generated by each of its companions. Ri denote the position vectors with respect to a generic inertial reference frame.

$$\begin{array}{ll}{F}_0\hfill & =F_{0,1}+F_{0,2}-F_{1,0}-F_{2,0},\hfill \\ F_1\hfill & =F_{1,0}+F_{1,2}-F_{0,1}-F_{2,1},\hfill \\ F_2\hfill & =F_{2,0}+F_{2,1}-F_{0,2}-F_{1,2},\hfill \end{array}$$(9)

where F_i,j to the tidal force acting on m_i due to the deformation in m_j. We note that the positive contributions in F_i are the direct effect of the deformation of the other bodies while the negative terms corresponds to the reaction of the force due to the deformation of m_i. These have the form $$F_{\text{i,j}}=-\frac{{\mathcal{K}}_{i,j}}{|{\Delta }_{\text{ij}}{|}^{10}}[2({\Delta }_{\text{ij}}\cdot \stackrel{.}{{\Delta }_{\text{ij}}}){\Delta }_{\text{ij}}+{\Delta }_{\text{ij}}{}^2({\Delta }_{\text{ij}}\times {\Omega }_{\text{j}}+\stackrel{.}{{\Delta }_{\text{ij}}})]$$(10)

(Mignard 1979), where ${\mathcal{K}}_{i,j}$ is a measure of the magnitude of the tidal force and is given by $${\mathcal{K}}_{i,j}=3\mathcal{G}{m}_i^2{\mathcal{R}}_j^5{k}_{2,j}\Delta t_j.$$(11)

As before, Ω_j is the spin angularvelocity of m_j and is assumed parallel to the orbital angular momentum. We have neglected the tidal contributions which arise from the zero-order potential since its effect is restricted to a precession of the pericenters and does not introduce any secular changes in the spins, semimajor axes or eccentricities.

2.3 The rotational dynamics

While the orbital dynamics can be obtained solving the equations of motion (Eq. (7)), the time variation of the spins are deduced from the conservation of the total angular momentum L_tot. Since we assumed rotations perpendicular to the common orbital plane, $$L_{\text{tot}}=L_{\text{orb}}+{\sum }_{i=0}^2{C}_i{\Omega }_{\text{i}}=\text{const}.,$$(12)

where C_i is the principal moment of inertia of m_i. In turn, the orbital angular momentum in the inertial reference frame is given by $$L_{\text{orb}}={\sum }_{i=0}^2{m}_i({\text{R}}_{\text{i}}\times \stackrel{.}{{\text{R}}_{\text{i}}}).$$(13)

Differentiating this equation with respect to time and substituting expressions (Eq. (7)) for the accelerations $\stackrel{\mathrm{..}}{{\text{R}}_{\text{i}}}$, we obtain $$\begin{array}{lll}\stackrel{.}{L_{\text{orb}}}\hfill & =\hfill & {\Delta }_{10}\times F_{\mathrm{1,0}}+{\Delta }_{20}\times F_{\mathrm{2,0}}\hfill \\ \hfill & \hfill & +{\Delta }_{01}\times F_{\mathrm{0,1}}+{\Delta }_{21}\times F_{\mathrm{2,1}}\hfill \\ \hfill & \hfill & +{\Delta }_{02}\times F_{\mathrm{0,2}}+{\Delta }_{12}\times F_{\mathrm{1,2}}.\hfill \end{array}$$(14)

Furthermore, assuming that the variation in the spin angular momenta of the body m_j is only due to the terms in ${\stackrel{.}{L}}_{\text{orb}}$ associated with its deformation, we obtain $$C_j\stackrel{.}{{\Omega }_{\text{j}}}=-{\sum }_{i\ne j}{\Delta }_{\text{ij}}\times F_{\text{i,j}}.$$(15)

Note that in the limit where the physical radius of m_j reduces to zero (i.e. ${\mathcal{R}}_j=0$), the tidal terms F_i,j are also zero for all i≠j, and Eq. (15) is automatically satisfied.

Finally, using expression (Eq. (10)) for the tidal forces, the time evolution of the spin vectors are given by $$\frac{\text{d}{\Omega }_{\text{j}}}{\text{d}t}=\frac{1}{C_j}{\sum }_{i\ne j}\frac{K_{i,j}}{|{\Delta }_{\text{ij}}{|}^6}[\frac{{\Delta }_{\text{ij}}\times \stackrel{.}{{\Delta }_{\text{ij}}}}{|{\Delta }_{\text{ij}}{|}^2}-{\Omega }_{\text{j}}].$$(16)

Contrary to the two-body case (e.g. Ferraz-Mello et al. 2008), the time derivative of the spin is given by the sum of two distinct terms. Depending on the magnitudes of each tidal term, it is not immediately obvious what would be the equilibrium rotational frequencies associated with stationary solutions.

3 Numerical simulations

In order to study the dynamical predictions of our model, we analyse the tidal evolution of a 3-body system consisting of a single planet around a binary star. The orbital and rotational evolution will be followed solving the equations of motion (Eq. (7)) for the orbit and Eq. (16) for each of the spins.

As before, we choose the Kepler-38 system as a test case, previously studied in Zoppetti et al. (2018) using a simpler tidal model. Nominal values for system parameters and initial orbital elements are detailed in Table 1. Stellar masses and radii were taken from Orosz et al. (2012), while the value of m₂ was estimated from the semi-empirical mass-radius from Mills & Mazeh (2017). The orbital elements of the secondary star respect to m₀ are those expected during the early stages of the system before tidal interactions had time to act (see Zoppetti et al. 2018), assuming tidal parameters and moments of inertia equal to those given in the table.

The orbital parameters for the planet are similar to those presented by Orosz et al. (2012), while the value of ${Q^{\prime }}_2$ is consistent with rocky bodies (Ferraz-Mello et al. 2008). However, it is important to stress that there is little dynamical constraint on the values of the planetary tidal parameters; the value adopted here is for illustrative purposes only. Finally, the parameters highlighted with an asterisk were varied in our different simulations.

We will focus our attention on two different timescales: (i) an early stage (up to ~1–2 Gyr) characterized by the evolution of the rotation rates towards stationary solutions, and (ii) the subsequent long-term dynamical orbital evolution of the system. In this second part we will concentrate primarily on the orbital migration and eccentricity damping of the planet.

3.1 Early dynamical evolution

Figure 4 shows the early rotational and orbital evolution of the binary stars and the planet. Except for the spin rates, all initial conditions and system parameters were taken equal to the nominal values of summarized in Table 1.

We begin our analysis with the binary components, shown in the left-hand plots of Fig. 4. The blue curves correspond to initial spin rates for both stellar components equal to Ω₀ = Ω₁ = n₁ ∕ 10 (i.e. slow rotators), while the black curves show results where the star were considered initially fast rotators: Ω₀ = Ω₁ = 10n₁. Regardless of the initial spin, both stars reach a pseudo-synchronization state in a few Gyrs, with a final rotational frequency equal to the value predicted by two-body tidal models: ${\Omega }_0/n_1={\Omega }_1/n_1=1+6e_1^2$ (e.g. Ferraz-Mello et al. 2008).

If the stars were initially super-synchronous, the change in spin rates deliver angular momenta to their orbit (Eq. (15)), increasing the semimajor axis a₁ and eccentricity e₁. The opposite effect is observed if the stars were initially sub-synchronous: the orbit delivers angular momenta to the stars to increase their spins, decreasing its semimajor axis and eccentricity. Once the rotational stationary solution is attained, the subsequent dynamical effect of the stellar tides acts to reduce the semimajor axis and damp the eccentricity until the circularization is reached (Correia et al. 2016; Hut 1980). Due to its small mass, the presence of the CB planet has no noticeable influence on the tidal evolution of the binary.

The right-hand panels of Fig. 4 show the dynamical evolution of the planetary spin (top panel) and orbit (middle and bottom plots). As before we considered two different initial spin rates: Ω₂ ∕ n₂ = 10 is shown in black while Ω₂ ∕ n₂ = 0.1 in green. The stellar spins were taken equal to the nominal values. We found no appreciable change in the time evolution of the semimajor axis or eccentricity regardless of the initial spins and, as seen in the middle and lower panels, both curves are practically indistinguishable.

Concerning the evolution of the planetary spin, both initial conditions reach stationary values much faster than the stars (typically ina few Myrs), although the equilibrium value is sub-synchronous and significantly displaced with respect to the two-body expectation (red horizontal line). This behaviour will be discussed in detail in Sect. 4.1 and constitutes a new finding.Instead of the super-synchronous stationary solutions found in classical tidal models for eccentric orbits, the interacting binary system leads to a stable sub-synchronous state which does not change even after the stars themselves evolve towards their rotational stationary spins.

Fig. 4 Early tidal evolution of a circumbinary system. In all the panels, the black curve represents the results of our reference simulation (Table 1). Left: dynamical evolution of the binary, showing the spin rate (top), semimajor axis (middle) and eccentricity (bottom panel). The results depicted in blue consider initially slow-rotating stars with Ω0 ∕ n1 = Ω1 ∕ n1 = 0.1 (at t = 0), while those in black correspond to primordial fast rotators Ω0 ∕ n1 = Ω1 ∕ n1 = 10. Right: evolution of the planetary spin and orbit. Black (respectively green) curves correspond to initial super-synchronous (respectively sub-synchronous) planetary spin rates. Time variation of the semimajor axis a2 and eccentricity e2 are practically equal in both cases (middle and bottom panels).

3.2 Long-term orbital evolution

Figure 5 shows three different long-term simulations, integrated over timescales comparable with the estimated age of Kepler-38 system (Zoppetti et al. 2018). All system parameters and initial conditions were chosen equal to their nominal values (Table 1) except for those described in the left-hand panels of each set. In all cases the planetary spin reached a sub-synchronous stationary solution early in the simulation; thus we concentrate on the orbital elements: semimajor axis a₂ in the left-hand plots, eccentricity e₂ in the centre graphs, and difference between longitudes of pericenter Δϖ = ϖ₂ − ϖ₁ in the right-hand graphs. Results after the application of the low-pass filter are shown in darker curves for a₂ and e₂.

The black curves in the top panels correspond to the results of our reference simulation (see Table 1) while in the cyan curves we consider a more dissipative planet with ${Q^{\prime }}_2=1$. The middle panels show results considering a more eccentric initial orbit e₂ (0) = 0.1, again for the same two values of the tidal parameter. Finally, in the lower panel we analyse the case in which the stars in the binary are not tidally interacting. This scenario correspond to setting F_0,1 = F_1,0 = 0 (see Eq. (10)) in our code. Results with non-tidally interacting stars are shown in magenta, while cyan curves repeat the results of our simulation with tidal effects for the stars.

Independently of the adopted tidal parameter ${Q^{\prime }}_2$, the planet is always observed to migrate outwards, marking a second distinct difference with respect to expectations from classical two-body tidal models. This result was already described in Zoppetti et al. (2018), although in that case we used a simpler and non-consistent tidal model. Lower values of ${Q^{\prime }}_2$ (cyan curvesin the upper and middle panels) lead to more larger excursions in semimajor axis, ultimately leading to scattering in a high-orderMMR and temporary excitation of the eccentricity. The magenta curve in the lower panels show that the outward migration is not a consequence of tidal effects in the stars, but seems to be independent of their tidal evolution.

The planetary eccentricity, on the other hand, always seems to decrease, as long as not mean-motion resonances are encountered. For low initial values of e₂ (upper panels of Fig. 5) the planet and secondary star enter an aligned secular mode (Michtchenko & Malhotra 2004) in which Δϖ librates around zero. The amplitude of oscillation increases for larger initial eccentricities until Δϖ is observed to circulate for e₂(0) = 0.1. However, the libration and circulation are purely kinematic and the difference in behaviour is related to the amplitude of oscillation of the eccentricity. Regardless, these results seem to indicate that an analytical model for the tidal evolution of these type of systems must include terms involving the secular angle Δϖ, even if the tidal evolution timescales are much longer than those associated with the precession of pericenters ϖ₁ and ϖ₂.

4 Analytical secular model

In order to construct an analtical model from the equations of motion (Eqs. (7) and (16)), we first introduce a Jacobi reference frame for the position and velocity vectors of the bodies. In terms of the inertial coordinates R_i, the positions of the masses in Jacobi coordinates are given by: $$\begin{array}{lll}{\rho }_0\hfill & =\hfill & \frac{1}{{\sigma }_2}(m_0{\text{R}}_0+m_1{\text{R}}_1+m_2{\text{R}}_2),\hfill \\ {\rho }_1\hfill & =\hfill & {\text{R}}_1-{\text{R}}_0,\hfill \\ {\rho }_2\hfill & =\hfill & {\text{R}}_2-\frac{1}{{\sigma }_1}(m_0{\text{R}}_0+m_1{\text{R}}_1),\hfill \end{array}$$(17)

where $${\sigma }_i={\sum }_{k=0}^i{m}_k.$$(18)

Analogous expressions relate the velocities vectors in both reference systems.

4.1 Secular evolution of the planetary spin

The rate of change of the rotational frequency of the planet can be obtained by expanding the position and velocity vectors in Eq. (16) in α = a₁ ∕ a₂ and the eccentricities, and averaging with respect to both mean longitudes. Up to second order, we have: $$<\frac{\text{d}{\Omega }_2}{\text{d}t}>=\frac{1}{2C_2{a}_2^6}{\sum }_{i,j,k=0}^2{A}_{i,j,k}^{(s)}{\alpha }^i{e}_1^j{e}_2^k,$$(19)

where the non-zero coefficients are given by $$\begin{array}{lll}{A}_{0,0,0}^{(s)}\hfill & =\hfill & 2(K_{0,2}+K_{1,2})(n_2-{\Omega }_2),\hfill \\ A_{2,0,0}^{(s)}\hfill & =\hfill & 6({\gamma }_0^2{K}_{0,2}+{\gamma }_1^2{K}_{1,2})(4_2-n_1-3{\Omega }_2),\hfill \\ A_{2,2,0}^{(s)}\hfill & =\hfill & 3({\gamma }_0^2{K}_{0,2}+{\gamma }_1^2{K}_{1,2})(12n_2+n_1-9{\Omega }_2),\hfill \\ A_{0,0,2}^{(s)}\hfill & =\hfill & 3(K_{0,2}+K_{1,2})(9n_2-5{\Omega }_2),\hfill \\ A_{2,0,2}^{(s)}\hfill & =\hfill & 12({\gamma }_0^2{K}_{0,2}+{\gamma }_1^2{K}_{1,2})(44n_2-7n_1-21{\Omega }_2),\hfill \\ A_{1,1,1}^{(s)}\hfill & =\hfill & 9(8n_2-5{\Omega }_2)({\gamma }_0{K}_{0,2}+{\gamma }_1{K}_{1,2})\cos (\Delta \varpi ),\hfill \end{array}$$(20)

with $${\gamma }_0=\frac{m_1}{{\sigma }_1};{\gamma }_1=-\frac{m_0}{{\sigma }_1}.$$(21)

The stationary spin rate $<{\Omega }_2{>}_{\text{stat}}$ predicted by this equation can be easily calculated by equating Eq. (19) to zero. The explicit form of the equilibrium rotational frequency was found to be $$<{\Omega }_2{>}_{\text{stat}}\text{ \&}=\text{ \&}(1+6e_2^2)n_2-6\frac{{\gamma }_0^2{\gamma }_1^2}{{\gamma }_0^2+{\gamma }_1^2}\left(n_1-n_2\right){\alpha }^2+3\frac{{\gamma }_0^2{\gamma }_1^2}{{\gamma }_0^2+{\gamma }_1^2}\left(n_1+3n_2\right){\alpha }^2{e}_1^2-3\frac{{\gamma }_0^2{\gamma }_1^2}{{\gamma }_0^2+{\gamma }_1^2}\left(13n_1-41n_2\right){\alpha }^2{e}_2^2+\frac{27}{2}\frac{{\gamma }_0{\gamma }_1({\gamma }_0+{\gamma }_1)}{{\gamma }_0^2+{\gamma }_1^2}{n}_2\alpha e_1{e}_2\cos (\Delta \varpi ).$$(22)

In the limit case in which the mass of one of the stars reduces to zero we recover the classical two-body super-synchronous stationary solution $<{\Omega }_2{>}_{\text{stat}}=(1+6e_2^2)n_2$ (Ferraz-Mello et al. 2008; Correia et al. 2011). On the other hand, we can observe that for low binary and planetary eccentricities (e₁, e₂ → 0), the CB planet stationary solution is sub-synchronous by a factor that decreases proportional to α² as we moveoutwards from the binary, and is maximum for equal-mass stars m₁ = m₀.

Figure 6 shows two colour plots with the value of $<{\Omega }_2{>}_{\text{stat}}$ as a function of different system parameters and eccentricities (assumed constant). The top frame shows the dependence of the equilibrium spin rate of the planet with the distance from the binary system and the mass of the secondary star. Except for initial conditions very close to the binary or m₁ ∕ m₀ ≲ 0.1, the estimated value of $<{\Omega }_2{>}_{\text{stat}}$ is always sub-synchronous with respect to the mean orbital frequency. The dashed black curve corresponds to the equilibrium value of the spin as obtained from the two-body problem, that is $<{\Omega }_2{>}_{\text{stat}}/n_2=1+6e_2^2$. Our model predicts lower values for practically all values of the system parameters, at least for the nominal eccentricities. This seems to imply that even a low-mass secondary, or even a large interior planet may counteract the super-synchronous state deduced from the two-body solution and lead to appreciable differences in the rotational dynamics.

The dependence of $<{\Omega }_2{>}_{\text{stat}}$ with the eccentricities is analysed in the bottom frame of Fig. 6. We note that the sub-synchronous equilibrium state is only observed for low eccentricities of the planet, typically e₂ ≲ 0.1–0.15, while super-synchronous states may be attained form more eccentric planets. However, since we expect tidal effects to damp the eccentricity, it appears that $<{\Omega }_2{>}_{\text{stat}}<n_2$ should probably be the most common configuration in real-life systems. Finally, we observe little sensitivity of the equilibrium spin with respect to the eccentricity of the binary.

In order to test the validity and precision of our analytical model, Fig. 7 shows four sets of different N-body simulations of the evolution of the planetary spin, considering binaries with different mass ratios and planets in orbits with different initial eccentricities. All results were digitally filtered to remove short-period variations.

The top right-hand panel uses initial conditions from Table 1 while the bottom right-hand panel considers a more eccentric CB planet. The left panels explore the case in which the mass of the secondary star is smaller than the nominal value. In every case the black curves correspond to initially super-synchronous planets while the green curves correspond to initially sub-synchronous CB planets. In dashed yellow curve, we show the synchronization spin predicted by our model (Eq. (22)) while in dashed red curve we compare with the stationary 2-body solution.

In accordance with the initial simulations presented in the previous section, the planetary spin reaches a stationary state rapidly, typically in about 10⁵ yr, and our model seems to reproduce the equilibrium behaviour extremely well. In the case of low-massive secondary star (left panels), the synchronization spin is very close to that predicted by the 2-body model. However, when we consider binaries with mass ratios similar to Kepler-38 system, the synchronization spin is very different: sub-synchronous by an amount that can be very large for binaries with stars of comparable mass. Since the gravitational interaction causes long-term (secular) variations in the eccentricity of the planet, the value of Ω₂ also suffers periodic oscillations.

Finally, as can be observed from Eq. (22), the stationary spin solution for the CB planets is not a function of the planetary mass m₂ nor of the physical radii of the bodies. Thus, if we assume that all currently known circumbinary planets have reached their stationary spin, we can predict their current rotational period just from the stellar masses and planetary orbits. As an example, considering its maximum possible eccentricity (Orosz et al. 2012) and that the planet is in an aligned secular mode (Zoppetti et al. 2018), we estimate the rotation period of the planet in the Kepler-38 system in P_K38 ≃ 118 days, about a 12% higher than the one predicted by the 2-body synchronization model.

Fig. 5 Long-term orbital tidal evolution of our planet in our Kepler-38-like system. Except for the parameters inlaid in the left-hand plots, all parameters and initial conditions were taken equal to those in Table 1. Light-tone curves for a2 and e2 show osculating values while darker curves correspond to mean elements obtained from a digital filter. The magenta curves in the lower panels are the result of a simulation disregarding tidal interaction between the stars.

Fig. 6 Stationary planetary spin $<{\Omega }_2{>}_{\text{stat}}/n_2$ as function of semimajor axis ratio α and mass of secondary star (top), and as function of eccentricities (bottom). All other system parameters were chosen equal to those given in Table 1. The nominal parameters for Kepler-38 are highlighted with a filled white circle and marked as “K38”. Dashed curve in the top graph corresponds to $<{\Omega }_2{>}_{\text{stat}}/n_2=1+6e_2^2$.

Fig. 7 N-body simulation of the spin evolution of fictitious CB planets, considering binaries with different mass ratios (different columns) and different initial eccentricity for the planets (different rows). In all the panels, the black curves correspond to the evolution of an initially super-synchronous planet while the green curves represent the initially sub-synchronous case. The dashed yellow curves are the stationary spins predicted by our model (Eq. (22)) while the dashed red curves are the 2-body stationary solution.

4.2 Variational equations for the orbital evolution

Having developed an analytical model for the rotational dynamics, we turn our attention to the time evolution of the semimajor axis a₂ and eccentricity e₂. As before, we will focus on the planetary orbit, although analogous expressions can be found also for the binary.

Following Beutler (2005), the variational equation for the semimajor axis in the Jacobi reference frame may be written as $$\frac{\text{d}{a}_2}{\text{d}t}=\frac{2a_2^2}{\mathcal{G}{\sigma }_2}({\stackrel{.}{\rho }}_2\cdot \delta f_2),$$(23)

where δf₂ is the total tidal force (per unit mass) affecting the 2-body motion of the planet around the center of mass of m₀ and m₁, and has the form: $$\delta f_2=\frac{F_2}{m_2}-\frac{1}{{\sigma }_1}(F_0+F_1).$$(24)

Substituting Eq. (9) in order to express the total force in terms of the individual two-body tidal interactions, we obtain $$\delta f_2=\frac{1}{{\beta }_2}[(F_{\mathrm{2,0}}-F_{\mathrm{0,2}})+(F_{\mathrm{2,1}}-F_{\mathrm{1,2}})],$$(25)

where $${\beta }_i=\frac{m_i{\sigma }_{i-1}}{{\sigma }_i}$$(26)

is the reduced-mass (e.g. Beaugé et al. 2007). An analogous reasoning leads to a similar equation for the binary orbital evolution.

Expression (25) shows that the total tidal force δf₂ may be written in terms of differences of the type (F_2,j −F_j,2), where j = 0, 1. From Eqs (10), each of these differences may be explicitly written as $$\begin{array}{lll}{F}_{\text{2,j}}-F_{\text{j,2}}\hfill & =\hfill & -\frac{{\mathcal{K}}_j^{(+)}}{|{\Delta }_{2j}{|}^{10}}[2({\Delta }_{\text{2j}}\cdot \stackrel{.}{{\Delta }_{\text{2j}}}){\Delta }_{\text{2j}}\hfill \\ \hfill & \hfill & +{\Delta }_{\text{2j}}{}^2({\Delta }_{\text{2j}}\times {\overline{\Omega }}_2^{(j)}+\stackrel{.}{{\Delta }_{\text{2j}}})],\hfill \end{array}$$(27)

where we have defined $$K_j^{(+)}=K_{2,j}+K_{j,2},$$(28)

and a new “averaged” rotational frequency $${\overline{\Omega }}_2^{(\text{j})}=\frac{{\mathcal{K}}_{2,j}{\Omega }_{\text{j}}+{\mathcal{K}}_{j,2}{\Omega }_2}{K_{2,j}+K_{j,2}}.$$(29)

Notice that Eq. (27) has the same functional form as the tidal force in the 2-body problem (Eq. (10)) with a magnitude given by $K_j^{(+)}$ and a rotational frequency defined by ${\overline{\Omega }}_2^{(\text{j})}$. In the limit where m₁ → 0 and $R_1\to 0$, the term in the tidal force associated with $K_1^{(+)}$ becomes negligible and we recover the same expression as found in the classical 2-body case.

Writing the tidal forces in terms of Jacobi coordinates through Δ_2j = ρ₂ + γ_jρ₁, substituting in the Gauss equation (Eq. (23)), expanding in power series of α, e₁ and e₂ and, finally, averaging over the mean longitudes, we obtain: $$<\frac{\text{d}{a}_2}{\text{d}t}>=\frac{n_2}{G{m}_2{\sigma }_2{a}_2^4}{\sum }_{i=0}^4{\sum }_{j,k=0}^2{\sum }_{l=0}^1{A}_{i,j,k,l}^{(a)}{K}_l^{(+)}{\gamma }_l^i{\alpha }^i{e}_1^j{e}_2^k,$$(30)

where the non-zero coefficients are explicitly given by $$\begin{array}{lll}{A}_{0,0,0,l}^{(a)}\hfill & =\hfill & 2[{\overline{\Omega }}_2^{(l)}-n_2],\hfill \\ A_{2,0,0,l}^{(a)}\hfill & =\hfill & 2[12{\overline{\Omega }}_2^{(l)}+5n_1-17n_2],\hfill \\ A_{4,0,0,l}^{(a)}\hfill & =\hfill & 20[6{\overline{\Omega }}_2^{(l)}+5n_1-11n_2],\hfill \\ A_{2,2,0,l}^{(a)}\hfill & =\hfill & [36{\overline{\Omega }}_2^{(l)}-5n_1-51n_2],\hfill \\ A_{4,2,0,l}^{(a)}\hfill & =\hfill & 100[6{\overline{\Omega }}_2^{(l)}+n_1-11n_2],\hfill \\ A_{1,1,1,l}^{(a)}\hfill & =\hfill & 6[12{\overline{\Omega }}_2^{(l)}-19n_2]\cos (\Delta \varpi )\hfill \\ A_{3,1,1,l}^{(a)}\hfill & =\hfill & \frac{25}{2}[96{\overline{\Omega }}_2^{(l)}+32n_1-193n_2],\cos (\Delta \varpi )\hfill \\ A_{0,0,2,l}^{(a)}\hfill & =\hfill & [27{\overline{\Omega }}_2^{(l)}-46n_2],\hfill \\ A_{2,0,2,l}^{(a)}\hfill & =\hfill & [528{\overline{\Omega }}_2^{(l)}+5(44n_1-227n_2)],\hfill \\ A_{4,0,2,l}^{(a)}\hfill & =\hfill & 10[390{\overline{\Omega }}_2^{(l)}+(325n_1-1008n_2)].\hfill \end{array}$$(31)

Figure 8 shows the normalized value of $<\text{d}{a}_2/\text{d}t>$ in the (α, m₁ ∕ m₀) plane for three different values of the binary and planet eccentricities. For each value of m₁ the physical radius of the star was modified following the empirical rule $R_1\simeq 0.9m_1$. The nominal values are shown in the top panel, and the parameters corresponding to Kepler-38 highlighted with a white circle. All initial conditions and physical parameters leading to an inward orbital migration of the planet are coloured in tones of blue, while those leading to a secular increase of a₂ in tones of red. The limit between both is marked with a white curve.

Although the plots show some quantitative differences as function of the eccentricities, in all cases there seems to exist a lower value of m₁ ∕ m₀ above whichthe tidal interaction of the system leads to an outward migration of the planet. The critical value of m₁ appears to be larger for more eccentric binaries and lower for stars in almost circular orbits. As expected, as m₁ → 0 the migration is inward, in accordance with known results for the 2-body case.

It is necessary to point out that our analytical model was obtained through a Legendre expansion of the elliptic functions truncated at fourth-orderof α. Consequently, the results shown here and in Fig. 6 are not expected to be accurate (or even valid) for α → 1. We have nevertheless opted to include the complete range solely for illustrative purposes.

The time variation of the eccentricity e₂ may be found from the orbital angular momentum L₂ in the Jacobi reference frame. In the planar case, we have $$L_2={\beta }_2|({\rho }_2\times {\stackrel{.}{\rho }}_2)|={\beta }_2\sqrt{\mathcal{G}{\sigma }_2{a}_2(1-e_2^2)},$$(32)

whose time derivative due to tidal forces leads to $$\frac{1}{{\beta }_2}{\dot{L}}_2=\frac{\mathcal{G}{\sigma }_2{\beta }_2}{2L_2}((1-e_2^2)\frac{\text{d}{a}_2}{\text{d}t}-a_2\frac{\text{d}{e}_2^2}{\text{d}t})=|({\rho }_2\times \delta f_2)|.$$(33)

Extracting the eccentricity term, we finally obtain: $$\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)|].$$(34)

Introducing elliptic expansions in a similar manner as done for Eq. (30), and averaging over short-period terms, we obtain: $$<\frac{\text{d}{e}_2^2}{\text{d}t}>=\frac{n_2}{4G{m}_2{\sigma }_2{a}_2^5}{\sum }_{i=0}^4{\sum }_{j,k=0}^2{\sum }_{l=0}^1{A}_{i,j,k,l}^{(e)}{K}_l^{(+)}{\gamma }_l^i{\alpha }^i{e}_1^j{e}_2^k,$$(35)

where now the non-zero coefficients are given by $$\begin{array}{ll}{A}_{0,0,2,l}^{(e)}\hfill & =4[11{\overline{\Omega }}_2^{(l)}-18n_2],\hfill \\ A_{2,0,2,l}^{(e)}\hfill & =20[36{\overline{\Omega }}_2^{(l)}+15n_1-74n_2],\hfill \\ A_{4,0,2,l}^{(e)}\hfill & =40[139{\overline{\Omega }}_2^{(l)}+95n_1-282n_2],\hfill \\ A_{1,1,1,l}^{(e)}\hfill & =2[39{\overline{\Omega }}_2^{(l)}-54n_2]\cos (\Delta \varpi ),\hfill \\ A_{3,1,1,l}^{(e)}\hfill & =10[102{\overline{\Omega }}_2^{(l)}+34n_1-185n_2]\cos (\Delta \varpi ).\hfill \end{array}$$(36)

Contrary to da₂ ∕ dt, we found that the eccentricity of the planet is always damped, at least for the initial conditions and system parameters tested here.

Fig. 8 Normalized values of the secular rate of change of the planetary semimajor axis, as function of the binary mass ratio and α. Each panel shows results for different eccentricities, assumed fixed for this calculation. Blue tones denote regions where the planet experiences an inward orbital migration, while red tone identify regions where the migration is outward. The primordial parameters of Kepler-38 are again highlighted in the top panel with a filled white circle and marked as “K38”.

4.3 Comparisons with numerical integrations

To test the accuracy of our analytical model, for given initial conditions we compare the variation in planetary semimajor axis and eccentricity predicted by Eqs. (30) and (35) with the numerical results obtained using the original unaveraged Eqs. (23) and (34). We consider the nominal system parameters detailed in Table 1 but varied the planetary eccentricity and semimajor axis ratio α. For each we computed da₂ ∕ dt and de₂ ∕ dt as a function of the reduced mass $$\tilde{\mu }=\frac{m_1}{m_0+m_1}$$(37)

by varying m₁. Due to the rapid rotational synchronization timescales, we consider stationary spins for the stars and for the planet according to Eq. (22).

Results are shown in Fig. 9. In all the panels the colours represent different planetary eccentricities (e₂ = 0.01 in blue, e₂ = 0.05 in green and e₂ = 0.1 in red) while the type of curve makes reference to the calculation method (full line for numerical and dashed line for analytical). Different rows correspond to different values of α: the reference value in the bottom panels (α = 5 ∕ 16, see Table 1) and half the nominal value in top panels.

From the right panels we note that, as a result of the tidal interaction, the eccentricity of the planet always decreases with a rate that seems weakly dependent on the secondary mass. However, as in the 2-body case, e₂ decays morerapidly for eccentric planets. Thus, the effect of tides on the eccentricity of circumbinary planets is very similar to that in the case of bodies around single stars. In the absence of additional forces we expect the systems to evolve towards quasi-circular orbits. Since our analytical model only included terms up to second order in e_i, the accuracy decreases substantially for larger eccentricities, leading to an relative error of the order of 20% for e₂ ~ 0.1. A more complete model, perhaps including Mignard eccentricity functions Mignard (1980) are necessary for more eccentric orbits.

The rate of change of the semimajor axis (left-hand plots) shows a better agreement between our model and the full unaveraged equations, leading to practically the same magnitude in the derivatives even for moderate eccentricities. In particular, the values of the critical reduced mass ${\tilde{\mu }}_{\text{crit}}$ associated with the limit between inward and outward migration is very well reproduced.

Finally, Fig. 10 shows the dependence of ${\tilde{\mu }}_{\text{crit}}$ as function of α for different eccentricities. As before, calculations performed with the unaveraged equations are plotted in continuous lines, while dashed curves show results with the analytical model including terms up to fourth order in α. To test the necessity of such high orders, the dotted lines show analogous results, this time truncating the expansions at third order in the semimajor-axis ratio. While the precision of the fourth-order analytical model is very good up to α ~ 0.3, the truncated version shows a much smaller region of validity, reduced down to α ~ 0.1. Thus, systems such as Kepler-38 require a high-order model in order to reproduce the dynamics with a fair accuracy.

It is interesting to note that ${\tilde{\mu }}_{\text{crit}}$ increases for smaller values of α. In the limit when α → 0, we expect the system to behave as a planet orbiting a single star of mass m₀ + m₁ and all initial conditions should lead to an inward migration of the semimajor axis.

Fig. 9 Time derivative of the semimajor axis (left panels) and eccentricity variation (right panels) of a circumbinary planet at different distances from the binary: α = 5 ∕ 32 (top panels) and α = 5 ∕ 16 (bottom panels). Different colours are employed for different eccentricities (e2 = 0.01 in blue, e2 = 0.05 in green and e2 = 0.1 in red) and different type of curves make reference to the calculation method: numerical (full line) and analytical (dashed line).

5 Summary and discussion

In this workwe present a model for treating the tides in a circumbinary system with one planet, in which all bodies are assumed to be extended and tidally interacting. To built the model, we consider a weak friction regime where the tidal forces can be approximated by the classical expressions of Mignard (1979) and proceed in two steps:

First, we revisited the Mignard theory and studied which tidal forces have a net effect onto the dynamical evolution of the system. In the classical two-body problem, where we are computing the torques on the same body that exerts the deformation, the zero-order Mignard torques have zero net secular effect. We found that this torques also has a null effect on the third body, as long as there are no MMRs between m₁ and m₂. Thus, in thenon-resonant circumbinary problem, the only forces that should be taken into account are those that are applied on the same body that exerts the deformation. In a resonant case the zero-order torques may have important effects; their consequences will be the focus of a forthcoming work.

Secondly, we incorporate the tidal forces in the gravitational equations of motion in a self-consistent approach. Namely, we consider that each of the bodies is deformed by the other two and there is a reaction force for each tidal force applied. As a result, we obtain the spin evolution equation for the bodies and the orbital evolution equation for the planet.

We have undertaken a series of numerical simulations, considering Kepler-38 system as a working example, in order to compare theresults of this model with our previous work (Zoppetti et al. 2018). We observed that in the short-timescales the dynamics is dominated by the spin synchronization of the bodies: the planet, assumed a rocky body, synchronize very quickly (in ~ Myr) in a stationary spin lower than the orbital mean motion. On the other hand, the stars exhibit super-synchronous spins in values predicted by the 2-body classical problem. The subsequent orbital evolution of the binary is little affected by the planet and proceeds to a decrease in the semimajor axis a₁ and eccentricity e₁.

The long-term orbital evolution of the planet is curiously different: as a result of the tidal interaction the planet migrates outwards and the direction of migration is not dependent on the initial planetary eccentricity or the assumed planetary tidal parameter. Moreover, the outward migration is also not an indirect effect of the migration of the binary, but observed even if the tidal evolution of the stars is neglected.

During the tidal migration, the eccentricity of the planet oscillates around the force eccentricity, which decreases as we move away from the binary (Leung & Lee 2013). For some initial conditions, we found that the difference of pericenter angle Δϖ librates around zero. Thus, when studying the secular tidal evolution of circumbinary planets, the usual procedure of averaging over the longitudes of pericenters may not be accurate.

To better understand the numerical results, we constructed an analytical secular model expanding the full spin and orbital equations of motion and averaging only over the mean longitudes. Regarding the spins, the simplicity of the full equation, allows us to expand only up to second-order in α and the eccentricities e₁ and e₂. The resulting expressions showed a very good agreement with N-body simulations. We furthermore obtained a simple equation estimating the stationary spin of CB planets that is not dependent on the planetary mass. If we assume that their spins have reached their equilibrium state, this would allow us to predict the rotation period of almost all circumbinary systems requiring only knowledge of the stellar masses and the orbital configuration of its members. Our analytical approach was validated comparing the planetary stationary spin of the numerical simulation with those predicted by our analytical equations.

Contrary to the spins, the analytical model for the orbital evolution required an expansion in the semimajor-axis ratio up to fourth-order in α. We maintained the eccentricities up to second order; however, latter simulations showed that higher orders are probably needed in systems with moderate-to-high eccentricities.

Regarding the eccentricity evolution, we found that the tidal forces on the CB planet always seem to act circularizating its orbit. We observed a strong dependence on the eccentricities but only a marginal dependence on the mass ratio of the stellar components. On the other hand, the complex dependence of the planetary semimajor axis evolution with the mass of the stars is reflected in the fact that the direction of migration depends on the binary mass ratio: for binaries in which the secondary star is much less massive, even the case in which the secondary companion is a planet, the tidal migration direction is inward. However, when the mass of both stars are of the same order the planet migrates outwards. The critical value of mass ratio for which the direction of migration changes sign is dependent on the planetary eccentricity and also on the position of the CB planet but can be predicted very accurately with our model.

The magnitude of the semimajor-axis variation is also very sensitive to the planetary eccentricity and proximity to the binary, but mainly dominated by the amount of energy that is dissipated in the planet due to tides. This quantity is very uncertain; however, the unexpected outward tidal migration of CB planet seems to be only dependent on the stellar masses and system configuration. A preliminary application of our model to other observed Kepler systems seems to indicate that many systems could also have suffered an outward tidal migration.

Fig. 10 Critical value of $\tilde{\mu }$ above which tidal effects on the planet lead to outward orbital migration. Different colours represent different eccentricities for the planet (same as in Fig. 9) and different types of curves refer to different calculation method: numerical (continuous curves), analytical up to fourth order in α (dashed) and analytical up to third order in the semimajor-axis ratio (dotted line).

Acknowledgements

The authors are very grateful to M. Efroimsky and S. Ferraz-Mello for their helpful comments and suggestions. We wish to express our gratitude to IATE for an extensive use of their computing facilities, without which this work would not have been possible. This research was funded by CONICET, SECYT/UNC, FONCYT and FAPESP (Grant 2016/20189-9).

References

All Tables

All Figures

	Fig. 1 Tidal-lagged bulge generated on m0 due to m1 and its effecton test body m2, from m0 -centric coordinate frame.
In the text

Fig. 2 Secular normalized torques of zero-order |⟨T0(r2)⟩| (left column) and first-order |⟨T1(r2)⟩| (right column), computed on m2 due to the tidal deformation on m0 induced by m1, plotted as a function of the mean-motion ratio n1 ∕ n2. We considered m0 = 1, a1 = 1 and varied a2 to include orbits both interior and exterior to m1. Upper panels: circular orbits (e1 = e2 = 0). Lower panels: eccentric orbits with e1 = e2 = 0.1. Light brown vertical lines highlight the location of some important MMRs. We note that the first-order torques in the right panels are also normalized respect to the time-lag Δt0.

In the text

	Fig. 3 Tidally interacting three-body system. The tidal bulge generated on each body (filled grey ellipsoids) are the sum of the deformations generated by each of its companions. Ri denote the position vectors with respect to a generic inertial reference frame.
In the text

Fig. 4 Early tidal evolution of a circumbinary system. In all the panels, the black curve represents the results of our reference simulation (Table 1). Left: dynamical evolution of the binary, showing the spin rate (top), semimajor axis (middle) and eccentricity (bottom panel). The results depicted in blue consider initially slow-rotating stars with Ω0 ∕ n1 = Ω1 ∕ n1 = 0.1 (at t = 0), while those in black correspond to primordial fast rotators Ω0 ∕ n1 = Ω1 ∕ n1 = 10. Right: evolution of the planetary spin and orbit. Black (respectively green) curves correspond to initial super-synchronous (respectively sub-synchronous) planetary spin rates. Time variation of the semimajor axis a2 and eccentricity e2 are practically equal in both cases (middle and bottom panels).

In the text

Fig. 5 Long-term orbital tidal evolution of our planet in our Kepler-38-like system. Except for the parameters inlaid in the left-hand plots, all parameters and initial conditions were taken equal to those in Table 1. Light-tone curves for a2 and e2 show osculating values while darker curves correspond to mean elements obtained from a digital filter. The magenta curves in the lower panels are the result of a simulation disregarding tidal interaction between the stars.

In the text

Fig. 6 Stationary planetary spin $<{\Omega }_2{>}_{\text{stat}}/n_2$ as function of semimajor axis ratio α and mass of secondary star (top), and as function of eccentricities (bottom). All other system parameters were chosen equal to those given in Table 1. The nominal parameters for Kepler-38 are highlighted with a filled white circle and marked as “K38”. Dashed curve in the top graph corresponds to $<{\Omega }_2{>}_{\text{stat}}/n_2=1+6e_2^2$.

In the text

Fig. 7 N-body simulation of the spin evolution of fictitious CB planets, considering binaries with different mass ratios (different columns) and different initial eccentricity for the planets (different rows). In all the panels, the black curves correspond to the evolution of an initially super-synchronous planet while the green curves represent the initially sub-synchronous case. The dashed yellow curves are the stationary spins predicted by our model (Eq. (22)) while the dashed red curves are the 2-body stationary solution.

In the text

Fig. 8 Normalized values of the secular rate of change of the planetary semimajor axis, as function of the binary mass ratio and α. Each panel shows results for different eccentricities, assumed fixed for this calculation. Blue tones denote regions where the planet experiences an inward orbital migration, while red tone identify regions where the migration is outward. The primordial parameters of Kepler-38 are again highlighted in the top panel with a filled white circle and marked as “K38”.

In the text

Fig. 9 Time derivative of the semimajor axis (left panels) and eccentricity variation (right panels) of a circumbinary planet at different distances from the binary: α = 5 ∕ 32 (top panels) and α = 5 ∕ 16 (bottom panels). Different colours are employed for different eccentricities (e2 = 0.01 in blue, e2 = 0.05 in green and e2 = 0.1 in red) and different type of curves make reference to the calculation method: numerical (full line) and analytical (dashed line).

In the text

Fig. 10 Critical value of $\tilde{\mu }$ above which tidal effects on the planet lead to outward orbital migration. Different colours represent different eccentricities for the planet (same as in Fig. 9) and different types of curves refer to different calculation method: numerical (continuous curves), analytical up to fourth order in α (dashed) and analytical up to third order in the semimajor-axis ratio (dotted line).

In the text