A self-consistent weak friction model for the tidal evolution of circumbinary planets

TL;DR

This paper introduces a comprehensive weak-friction tidal model for circumbinary planets, analyzing their orbital and spin evolution, and providing analytical expressions for their long-term behavior.

Contribution

It develops a self-consistent tidal evolution model considering all bodies, deriving analytical formulas for spin rates and orbital migration in circumbinary systems.

Findings

Planet reaches a sub-synchronous equilibrium spin faster than stars.

Eccentricities are damped over time for all bodies.

Planet's semimajor axis can migrate inward or outward depending on parameters.

Abstract

We present a self-consistent model for the tidal evolution of circumbinary planets. Based on the weak-friction model, we derive expressions of the resulting forces and torques considering complete tidal interactions between all the bodies of the system. Although the tidal deformation suffered by each extended mass must take into account the combined gravitational effects of the other two bodies, the only tidal forces that have a net effect on the dynamic are those that are applied on the same body that exerts the deformation, as long as no mean-motion resonance exists between the masses. We apply the model to the Kepler-38 binary system. The evolution of the spin equations shows that the planet reaches a stationary solution much faster than the stars, and the equilibrium spin frequency is sub-synchronous. The binary components evolve on a longer timescale, reaching a super-synchronous solution very close to that derived for the 2-body problem. After reaching spin stationarity, the eccentricity is damped in all bodies and for all the parameters analyzed here. A similar effect is noted for the binary separation. The semimajor axis of the planet, on the other hand, may migrate inwards or outwards, depending on the masses and orbital parameters. In some cases the secular evolution of the system may also exhibit an alignment of the pericenters, requiring to include additional terms in the tidal model. Finally, we derived analytical expressions for the variational equations of the orbital evolution and spin rates based on low-order elliptical expansions in the semimajor axis ratio and the eccentricities. These are found to reduce to the 2-body case when one of the masses is taken equal to zero. This model allow us to find a close and simple analytical expression for the stationary spin rates of all the bodies, as well as predicting the direction and magnitude of the orbital migration.

Figures

Tables

Table 1. Table 1: Initial conditions for our reference numerical simulation, representing the primordial Kepler 38 system (Orosz et al., 2012 ; Zoppetti et al., 2018 ) . Orbital elements are given in a Jacobi reference frame. The parameters highlighted with an asterisk were varied in different simulations as indicated in the text.

body	$m_0$	$m_1$	$m_2$
mass	$0.949M_{\odot }$	$0.249M_{\odot }$	$10M_{\oplus }$
radius	$0.84R_{\odot }$	$0.272R_{\odot }$	$4.35R_{\oplus }$
$C_i/\left(m_i{ℛ}_i^2\right)$	$0.07$	$0.25$	$0.25$
$Q_i^{\prime }$	$1\times {10}^6$	$1\times {10}^6$	$1\times {10}^1$(*)
${\Omega }_i$	$10n_1$(*)	$10n_1$(*)	$10n_2$(*)
$a_i$ [AU]		$0.15$	$0.48$
$e_i$		$0.15$	$0.05$(*)

Equations

$$U(\mathbf{r},{\mathbf{r}}_1)=U^{(0)}(\mathbf{r},{\mathbf{r}}_1)+U^{(1)}(\mathbf{r},{\mathbf{r}}_1)+\mathcal{O}(\Delta t_0^2)$$

$$\displaystyle\begin{split}U^{(0)}({\bf{r}_{2}},{\bf{r}_{1}})&=\frac{\mathcal{G}m_{1}\mathcal{R}_{0}^{5}}{2r_{1}^{5}r_{2}^{5}}k_{2,0}\bigg{[}3{({\bf{r}_{2}}\cdot{\bf{r}_{1}})}^{2}-r_{2}^{2}r_{1}^{2}\bigg{]}\\ U^{(1)}({\bf{r}_{2}},{\bf{r}_{1}})&=\frac{3\mathcal{G}m_{1}\mathcal{R}_{0}^{5}}{r_{1}^{5}r_{2}^{5}}k_{2,0}\Delta t_{0}\bigg{[}\frac{({\bf{r_{1}}}\cdot\dot{{\bf{r}_{1}}})}{2r_{1}^{2}}[5{({\bf{r}_{2}}\cdot{\bf{r}_{1}})}^{2}-r_{2}^{2}r_{1}^{2}]\\ &\hskip 51.21504pt-({\bf{r}_{2}}\cdot{\bf{r}_{1}})[{\bf{r}_{1}}\cdot({\bf{\Omega}}_{0}\times{\bf{r}_{2}})+{\bf{r}_{2}}\cdot\dot{{\bf{r}_{1}}}]\bigg{]}\\ \end{split}$$

$$\mathbf{f}={\nabla }_{\mathbf{r}}(U^{(0)}+U^{(1)})={\mathbf{f}}^{(0)}+{\mathbf{f}}^{(1)}$$

$$\displaystyle\begin{split}{\bf f}^{(0)}&=\frac{3\mathcal{G}m_{1}\mathcal{R}_{0}^{5}}{2r_{2}^{5}r_{1}^{5}}k_{2,0}\bigg{[}2({\bf{r}_{2}}\cdot{\bf{r}_{1}}){\bf{r}_{1}}+\bigg{(}{r_{1}}^{2}-\frac{5}{r_{2}^{2}}{({\bf{r}_{2}}\cdot{\bf{r}_{1}})}^{2}\bigg{)}{\bf{r_{2}}}\bigg{]}\\ {\bf f}^{(1)}&=\frac{3\mathcal{G}m_{1}\mathcal{R}_{0}^{5}}{r_{2}^{5}r_{1}^{5}}k_{2,0}\Delta t_{0}\bigg{[}\frac{({\bf{r}_{1}}\cdot\dot{{\bf{r}_{1}}})}{r_{1}^{2}}[5{\bf{r}_{1}}({\bf{r}_{2}}\cdot{\bf{r}_{1}})-{\bf{r}_{2}}{r_{1}}^{2}]\\ &-[{\bf{r}}_{1}\cdot({\bf{\Omega}}_{0}\times{\bf{r}_{2}})+{\bf{r_{2}}}\cdot\dot{{\bf{r}_{1}}}]{\bf{r}_{1}}-({\bf{r}_{1}}\times{\bf{\Omega}}_{0}+\dot{{\bf{r}_{1}}})({\bf{r}_{2}}\cdot{\bf{r}_{1}})\\ &+\frac{5{\bf{r}_{2}}}{r_{2}^{2}}\bigg{[}({\bf{r}_{2}}\cdot{\bf{r}_{1}})[{\bf{r}_{1}}\cdot({\bf{\Omega}}_{0}\times{\bf{r_{2}}})+{\bf{r}_{2}}\cdot\dot{{\bf{r}_{1}}}]\\ &-\frac{({\bf{r}_{1}}\cdot\dot{{\bf{r}_{1}}})}{2r_{1}^{2}}[5{({\bf{r}_{2}}\cdot{\bf{r}_{1}})}^{2}-{r}^{2}_{2}{r_{1}}^{2}]\bigg{]}\bigg{]}.\\ \end{split}$$

$$\displaystyle\begin{split}{\bf{T}}^{(0)}({\bf{r}_{2}},{\bf{r}_{1}})&=\frac{3\mathcal{G}m_{1}k_{2,0}\mathcal{R}_{0}^{5}}{r_{2}^{5}r_{1}^{5}}({\bf{r}_{2}}\cdot{\bf{r}_{1}})({\bf{r}_{2}}\times{\bf{r}_{1}})\\ {\bf{T}}^{(1)}({\bf{r}_{2}},{\bf{r}_{1}})&=\frac{3\mathcal{G}m_{1}k_{2,0}\mathcal{R}_{0}^{5}}{r_{2}^{5}r_{1}^{5}}\Delta t_{0}\bigg{[}5\frac{({\bf{r}_{1}}\cdot\dot{{\bf{r}}_{1}})}{r_{1}^{2}}({\bf{r}_{2}}\cdot{\bf{r}_{1}})({\bf{r}_{2}}\times{\bf{r}_{1}})\\ &-[{\bf{r}_{1}}\cdot({\bf{\Omega}_{0}}\times{\bf{r}_{2}})+{\bf{r}_{2}}\cdot\dot{{\bf{r}}_{1}}]({\bf{r}_{2}}\times{\bf{r}_{1}})\\ &-({\bf{r}_{2}}\cdot{\bf{r}_{1}})[({\bf{r}_{2}}\cdot{\bf{\Omega}}_{0}){\bf{r}_{1}}-({\bf{r}_{2}}\cdot{\bf{r}_{1}}){\bf{\Omega}_{0}}+{\bf{r}_{2}}\times\dot{{\bf{r}}_{1}}]\bigg{]}.\\ \end{split}$$

$$⟨{\mathbf{T}}^{(i)}({\mathbf{r}}_2,{\mathbf{r}}_1)⟩=\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }{\mathbf{T}}^{(i)}({\mathbf{r}}_2(t),{\mathbf{r}}_1(t))dt,$$

$$\begin{array}{rl}{m}_0{\ddot{\mathbf{R}}}_0& =\frac{\mathcal{G}{m}_0{m}_1}{|{\mathbf{\Delta }}_{10}{|}^3}{\mathbf{\Delta }}_{10}+\frac{\mathcal{G}{m}_0{m}_2}{|{\mathbf{\Delta }}_{20}{|}^3}{\mathbf{\Delta }}_{20}+{\mathbf{F}}_0\\ m_1{\ddot{\mathbf{R}}}_1& =-\frac{\mathcal{G}{m}_0{m}_1}{|{\mathbf{\Delta }}_{10}{|}^3}{\mathbf{\Delta }}_{10}+\frac{\mathcal{G}{m}_1{m}_2}{|{\mathbf{\Delta }}_{21}{|}^3}{\mathbf{\Delta }}_{21}+{\mathbf{F}}_1\\ m_2{\ddot{\mathbf{R}}}_2& =-\frac{\mathcal{G}{m}_0{m}_2}{|{\mathbf{\Delta }}_{20}{|}^3}{\mathbf{\Delta }}_{20}-\frac{\mathcal{G}{m}_1{m}_2}{|{\mathbf{\Delta }}_{21}{|}^3}{\mathbf{\Delta }}_{21}+{\mathbf{F}}_2\end{array}$$

$${\mathbf{\Delta }}_{\mathbf{ij}}\equiv {\mathbf{R}}_{\mathbf{i}}-{\mathbf{R}}_{\mathbf{j}}.$$

$$\begin{array}{rl}{\mathbf{F}}_0& ={\mathbf{F}}_{0,1}+{\mathbf{F}}_{0,2}-{\mathbf{F}}_{1,0}-{\mathbf{F}}_{2,0}\\ {\mathbf{F}}_1& ={\mathbf{F}}_{1,0}+{\mathbf{F}}_{1,2}-{\mathbf{F}}_{0,1}-{\mathbf{F}}_{2,1}\\ {\mathbf{F}}_2& ={\mathbf{F}}_{2,0}+{\mathbf{F}}_{2,1}-{\mathbf{F}}_{0,2}-{\mathbf{F}}_{1,2},\end{array}$$

$$\displaystyle{\bf{F}_{i,j}}=-\frac{\mathcal{K}_{i,j}}{{|{\bf\Delta_{ij}}|}^{10}}\bigg{[}2({\bf\Delta_{ij}}\cdot{\bf\dot{\Delta}_{ij}}){\bf\Delta_{ij}}+{\bf\Delta_{ij}}^{2}({\bf\Delta_{ij}}\times{\bf{\Omega}_{j}}+{\bf\dot{\Delta}_{ij}})\bigg{]}$$

$${\mathcal{K}}_{i,j}=3\mathcal{G}{m}_i^2{\mathcal{R}}_j^5{k}_{2,j}\Delta t_j.$$

$${\mathbf{L}}_{\mathrm{tot}}={\mathbf{L}}_{\mathrm{orb}}+\sum _{i=0}^2{C}_i{\mathbf{\Omega }}_{\mathbf{i}}=const.,$$

$${\mathbf{L}}_{\mathrm{orb}}=\sum _{i=0}^2{m}_i({\mathbf{R}}_{\mathbf{i}}\times {\dot{\mathbf{R}}}_{\mathbf{i}}).$$

$$\begin{array}{rl}{\dot{\mathbf{L}}}_{\mathrm{orb}}& ={\mathbf{\Delta }}_{10}\times {\mathbf{F}}_{1,0}+{\mathbf{\Delta }}_{20}\times {\mathbf{F}}_{2,0}\\ & +{\mathbf{\Delta }}_{01}\times {\mathbf{F}}_{0,1}+{\mathbf{\Delta }}_{21}\times {\mathbf{F}}_{2,1}\\ & +{\mathbf{\Delta }}_{02}\times {\mathbf{F}}_{0,2}+{\mathbf{\Delta }}_{12}\times {\mathbf{F}}_{1,2}.\end{array}$$

$$C_j{\dot{\mathbf{\Omega }}}_{\mathbf{j}}=-\sum _{i\ne j}{\mathbf{\Delta }}_{\mathbf{ij}}\times {\mathbf{F}}_{\mathbf{i},\mathbf{j}}.$$

$$\frac{d{\bf\Omega_{j}}}{dt}=\frac{1}{C_{j}}\sum_{i\neq j}\frac{{\cal K}_{i,j}}{|{\bf\Delta_{ij}}|^{6}}\bigg{[}\frac{{\bf\Delta_{ij}}\times{\bf{\dot{\Delta}}_{ij}}}{|{\bf\Delta_{ij}}|^{2}}-{\bf\Omega_{j}}\bigg{]}.$$

$$\begin{array}{rl}{\mathbf{\rho }}_0& =\frac{1}{{\sigma }_2}(m_0{\mathbf{R}}_0+m_1{\mathbf{R}}_1+m_2{\mathbf{R}}_2)\\ {\mathbf{\rho }}_1& ={\mathbf{R}}_1-{\mathbf{R}}_0\\ {\mathbf{\rho }}_2& ={\mathbf{R}}_2-\frac{1}{{\sigma }_1}(m_0{\mathbf{R}}_0+m_1{\mathbf{R}}_1),\end{array}$$

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

$$⟨\frac{d{\Omega }_2}{dt}⟩=\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,$$

$$\begin{array}{rl}{A}_{0,0,0}^{(s)}& =2({\mathcal{K}}_{0,2}+{\mathcal{K}}_{1,2})(n_2-{\Omega }_2)\\ A_{2,0,0}^{(s)}& =6({\gamma }_0^2{\mathcal{K}}_{0,2}+{\gamma }_1^2{\mathcal{K}}_{1,2})(4n_2-n_1-3{\Omega }_2)\\ A_{2,2,0}^{(s)}& =3({\gamma }_0^2{\mathcal{K}}_{0,2}+{\gamma }_1^2{\mathcal{K}}_{1,2})(12n_2+n_1-9{\Omega }_2)\\ A_{0,0,2}^{(s)}& =3({\mathcal{K}}_{0,2}+{\mathcal{K}}_{1,2})(9n_2-5{\Omega }_2)\\ A_{2,0,2}^{(s)}& =12({\gamma }_0^2{\mathcal{K}}_{0,2}+{\gamma }_1^2{\mathcal{K}}_{1,2})(44n_2-7n_1-21{\Omega }_2)\\ A_{1,1,1}^{(s)}& =9(8n_2-5{\Omega }_2)({\gamma }_0{\mathcal{K}}_{0,2}+{\gamma }_1{\mathcal{K}}_{1,2})\cos (\Delta \varpi ),\end{array}$$

$${\gamma }_0=\frac{m_1}{{\sigma }_1};{\gamma }_1=-\frac{m_0}{{\sigma }_1}.$$

$$\displaystyle\begin{split}\big{<}\Omega_{2}\big{>}_{\rm stat}&=(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).\\ \end{split}$$

$$\frac{d{a}_2}{dt}=\frac{2a_2^2}{\mathcal{G}{\sigma }_2}({\dot{\mathbf{\rho }}}_2\cdot \delta {\mathbf{f}}_2)$$

$$\delta {\mathbf{f}}_2=\frac{{\mathbf{F}}_2}{m_2}-\frac{1}{{\sigma }_1}({\mathbf{F}}_0+{\mathbf{F}}_1).$$

$$\delta{\bf{f}_{2}}=\frac{1}{\beta_{2}}\bigg{[}({\bf{F}_{2,0}}-{\bf{F}_{0,2}})+({\bf{F}_{2,1}}-{\bf{F}_{1,2}})\bigg{]},$$

$${\beta }_i=\frac{m_i{\sigma }_{i-1}}{{\sigma }_i}$$

$$\displaystyle\begin{split}{\bf{F}_{2,j}}-{\bf{F}_{j,2}}&=-\frac{\mathcal{K}^{(+)}_{j}}{|{\bf{\Delta_{2j}}}|^{10}}\bigg{[}2({\bf\Delta_{2j}}\cdot{\bf\dot{\Delta}_{2j}}){\bf\Delta_{2j}}\\ &\hskip 56.9055pt+{\bf\Delta_{2j}}^{2}({\bf\Delta_{2j}}\times{\bf{\bar{\Omega}}^{(j)}_{2}}+{\bf\dot{\Delta}_{2j}})\bigg{]}\\ \end{split}$$

$${\mathcal{K}}_j^{(+)}={\mathcal{K}}_{2,j}+{\mathcal{K}}_{j,2}$$

$${\bar{\mathbf{\Omega }}}_2^{(\mathbf{j})}=\frac{{\mathcal{K}}_{2,j}{\mathbf{\Omega }}_{\mathbf{j}}+{\mathcal{K}}_{j,2}{\mathbf{\Omega }}_2}{{\mathcal{K}}_{2,j}+{\mathcal{K}}_{j,2}}.$$

$$\bigg{<}\frac{da_{2}}{dt}\bigg{>}=\frac{n_{2}}{{\cal G}m_{2}\sigma_{2}a_{2}^{4}}\sum_{i=0}^{4}\sum_{j,k=0}^{2}\sum_{l=0}^{1}A^{(a)}_{i,j,k,l}\,{\cal K}_{l}^{(+)}\gamma_{l}^{i}\alpha^{i}e_{1}^{j}e_{2}^{k},$$

Full text

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 $\sim 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 inward 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 modeled 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 Section 2 we present the model in two steps: in Section 2.1 we first discuss which tidal forces have a net effect onto the dynamical evolution of an 3-extended-body system while in Section 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 concentrate 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 Section 4, we construct 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 compare these predictions with full N-body simulations. Finally, Section 5 summarizes our main results and discusses their implications.

2 The model

Let us consider a binary system in which $m_{0}$ and $m_{1}$ are the masses of the stellar components and $m_{2}$ is a circumbinary planet. We suppose that all the bodies lie in the same orbital plane and their spins are perpendicular to it. We also assume that all the bodies are extended masses with physical radii ${\cal R}_{i}$ and are deformable due to tidal effects between them.

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

2.1 The Mignard forces revisited

We begin considering our three-body system with two simplifications. First, we will neglect the gravitational perturbations generated by $m_{2}$ on the other two bodies, as well as the effects of $m_{1}$ on $m_{2}$. Second, only $m_{0}$ will be assumed to be an extended mass while $m_{1}$ and $m_{2}$ will be taken as point masses. As a consequence of these approximations, the dynamics of both $m_{1}$ and $m_{2}$ around $m_{0}$ will be defined by the point-mass approximation plus the tidal deformation of $m_{0}$ generated solely by $m_{1}$. The role of $m_{2}$ is 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 Figure 1, where ${\bf r_{i}}$ are the $m_{0}$-centric position vectors of the other masses. Following Mignard (1979), the tidal bulge of $m_{0}$ considered is displaced with respect to the instantaneous position of $m_{1}$ by a constant time-lag $\Delta t_{0}$. We assume the lag is sufficiently small to expand the gravitational potential $U$ generated by $m_{0}$ in anywhere in the space up to first-order in $\Delta t_{0}$, such that

[TABLE]

where $U^{(0)}$ and $U^{(1)}$ are the expanded tidal potentials of order $\mathcal{O}(0)$ and $\mathcal{O}(\Delta t_{0})$, respectively. In particular, if we evaluate (1) on the position of $m_{2}$ (i.e. ${\bf{r}}={\bf{r}_{2}}$), we obtain

[TABLE]

where $\mathcal{G}$ is the gravitational constant, ${\bf{\Omega}}_{0}$ is the spin vector of $m_{0}$ and $k_{2,0}$ its the second degree Love number.

The tidal force per unit mass ${\cal{\bf f}}$ generated by $m_{0}$ at a generic position vector ${\bf r}$ can be obtained as

[TABLE]

where explicit expressions evaluated on $m_{2}$ are given by

[TABLE]

Finally, the torques per unit mass can be calculated as ${\bf{T}}({\bf{r}},{\bf{r}_{1}})\simeq{\bf{r}}\times({\bf f}^{(0)}+{\bf f}^{(1)})={\bf{T}}^{(0)}({\bf{r}},{\bf{r}_{1}})+{\bf{T}}^{(1)}({\bf{r}},{\bf{r}_{1}})$. As before, evaluating on the position of $m_{2}$ yields

[TABLE]

In the classical two-body tidal problem, the acceleration ${\bf f}$ and the torque ${\bf{T}}$ are computed on the position of the deforming body $m_{1}$. It is easy to see that in such a case, the zero-order torque reduces to zero (i.e. ${\bf{T}}^{(0)}({\bf{r}}={\bf{r}_{1}},{\bf{r}_{1}})={\bf{0}}$) and the only net contribution to the orbital and spin evolution (notwithstanding a precession term) stems from the the first-order expressions ${\bf f}^{(1)}$ and ${\bf{T}}^{(1)}$ (see equations (5) and (6) of Mignard (1979)). However, it is not immediately clear what occurs if ${\bf r}\neq{\bf r_{1}}$. In other words, we wish to analyze what are the (long-term) dynamical effects of a tidal bulge generated on $m_{0}$, due to the perturbing potential of $m_{1}$, on the orbit of another body $m_{2}$.

To address this question, let the $m_{0}$-centric orbits of $m_{i}$ be characterized by semimajor axis $a_{i}$, eccentricity $e_{i}$, mean longitude $\lambda_{i}$ and longitude of pericenter $\varpi_{i}$. Furthermore, let $n_{i}$ denote the mean-motion (orbital frequency) of each body. We then compute the net secular torques ($\langle\bf{T}^{(0)}({\bf{r}}_{2},{\bf{r}}_{1})\rangle$ and $\langle\bf{T}^{(1)}({\bf{r}}_{2},{\bf{r}}_{1})\rangle$) for different values of $a_{2}$, assuming fixed values for $(a_{1},e_{1},e_{2},\varpi_{1},\varpi_{2})$. The secular torques are calculated averaging over the short-period terms associated to $\lambda_{1}$ and $\lambda_{2}$. Since we will not restrict our analysis to non-resonant configurations between $m_{1}$ and $m_{2}$, we cannot assume that both mean longitudes are necessarily independent. We thus substitute the classical double averaging over $\lambda_{i}$ with a time average over time, such that

[TABLE]

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_{1}$ and $m_{2}$.

Results are shown in Figure 2 for $m_{0}=1$, $a_{1}=1$, and $\varpi_{1}=\varpi_{2}$ and $e_{1}=e_{2}$. The values of the averaged torques are normalized with respect to $m_{1}$ and $\Delta t_{0}$, and plotted as function of the mean-motion ratio $n_{1}/n_{2}$. The left-hand panels correspond to the zero-oder torque $\langle\bf{T}^{(0)}\rangle$, while the first-order contributions $\langle\bf{T}^{(1)}\rangle$ are depicted in the right-hand graphs.

In the upper panels we analyze the circular case ($e_{1}=e_{2}=0$) and in the lower panel eccentric orbits ($e_{1}=e_{2}=0.1$). Eccentricities are assume 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_{1}/n_{2}\rightarrow\infty$. This background trend is a consequence of the numerical approximation employed to evaluate the time integral (6), which basically consisted in 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_{2}$ for which $m_{2}$ receives a non-zero net torque corresponds to a $1/1$ MMR, that is in a coorbital position with the deforming body $m_{1}$. In particular, the case in which the position of both masses coincide (i.e. $\lambda_{1}=\lambda_{2}$) yields results analogous to those obtained from the 2-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_{2}$ located in a Trojan-like orbit with the other masses. Finally, although $|\langle{\bf{T}}^{0}({\bf{r}}_{2})\rangle|$ 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 mean-motion resonances and not only in the coorbital region. This seems to imply that the tidal deformation generated by $m_{1}$ on $m_{0}$ should affect the dynamical evolution of $m_{2}$ 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 Figure 2 were confirmed introducing elliptic expansions for the position and velocity vectors in the equations (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_{1}$ and $m_{2}$, the only tidal forces that need to be considered on $m_{i}$ are those stemming from the deformation that $m_{i}$ generates on $m_{j}$ and $m_{k}$ (with $i\neq j\neq 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 equations (5) and (6) of Mignard (1979). The effect of the tidal torques on resonant orbits will be investigated in a forthcoming work.

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 Figure 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 Section 2.1, we only need to consider the direct distortion between pairs.

Let us denote by ${\bf 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:

[TABLE]

where for compactness we have denoted the relative position vectors as

[TABLE]

In terms of the $m_{0}$-centric position vectors, these are simply given by ${\bf\Delta_{10}}={\bf r_{1}}$, ${\bf\Delta_{20}}={\bf r_{2}}$ and ${\bf\Delta_{21}}={\bf r_{2}}-{\bf r_{1}}$. 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

[TABLE]

where ${\bf{F}_{i,j}}$ to the tidal force acting on $m_{i}$ due to the deformation in $m_{j}$. Note that the positive contributions in $\bf{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

[TABLE]

(Mignard 1979), where $\mathcal{K}_{i,j}$ is a measure of the magnitude of the tidal force and is given by

[TABLE]

As before, ${\bf\Omega_{j}}$ is the spin angular velocity 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 (7), the time variation of the spins are deduced from the conservation of the total angular momentum ${\bf L_{\rm tot}}$. Since we assumed rotations perpendicular to the common orbital plane,

[TABLE]

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

[TABLE]

Differentiating this equation with respect to time and substituting expressions (7) for the accelerations ${\bf\ddot{R}_{i}}$, we obtain

[TABLE]

Furthermore, assuming that the variation in the spin angular momenta of the body $m_{j}$ is only due to the terms in ${\bf{\dot{L}}_{\rm orb}}$ associated to its deformation, we obtain

[TABLE]

Note than in the limit where the physical radius of $m_{j}$ reduces to zero (i.e. $\mathcal{R}_{j}=0$), the tidal terms ${\bf{F}_{i,j}}$ are also zero for all $i\neq j$, and equation (15) is automatically satisfied.

Finally, using expression (10) for the tidal forces, the time evolution of the spin vectors are given by

[TABLE]

Contrary to the 2-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 to stationary solutions.

3 Numerical simulations

In order study the dynamical predictions of our model, we analyze 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 (7) for the orbit and equation (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_{2}$ was estimated from the semi-empirical mass-radius from (Mills & Mazeh (2017)). The orbital elements of the secondary star respect to $m_{0}$ 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 asterisks were varied in our different simulations.

We will focus our attention on two different timescales: (i) an early stage (up to $\sim 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 Figure 4. The blue curves correspond to initial spin rates for both stellar components equal to $\Omega_{0}=\Omega_{1}=n_{1}/10$ (i.e. slow rotators), while the black curves show results where the star were considered initially fast rotators: $\Omega_{0}=\Omega_{1}=10n_{1}$. 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 2-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_{1}$ and eccentricity $e_{1}$. 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 Figure 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: $\Omega_{2}/n_{2}=10$ is shown in black while $\Omega_{2}/n_{2}=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 in a few Myrs), although the equilibrium value is sub-synchronous and significantly displaced with respect to the 2-body expectation (red horizontal line). This behavior will be discussed in detail in section 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.

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_{2}$ in the left-hand plots, eccentricity $e_{2}$ in the center graphs, and difference between longitudes of pericenter $\Delta\varpi=\varpi_{2}-\varpi_{1}$ in the right-hand graphs. Results after the application of the low-pass filter are shown in darker curves for $a_{2}$ and $e_{2}$.

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_{2}(0)=0.1$, again for the same two values of the tidal parameter. Finally, in the lower panel we analyze the case in which the stars in the binary are not tidally interacting. This scenario correspond to setting ${\bf{F}_{0,1}}={\bf{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 2-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 curves in the upper and middle panels) lead to more larger excursions in semimajor axis, ultimately leading to scattering in a high-order MMR and temporary excitation of the eccentricity. The magenta curve in the lower panels show that the outwards 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_{2}$ (upper panels of Figure 5) the planet and secondary star enter an aligned secular mode (Michtchenko & Malhotra (2004)) in which $\Delta\varpi$ librates around zero. The amplitude of oscillation increases for larger initial eccentricities until $\Delta\varpi$ is observed to circulate for $e_{2}(0)=0.1$. However, the libration/circulation is purely kinematic and the difference in behavior 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 $\Delta\varpi$, even if the tidal evolution timescales are much longer than those associated to the precession of pericenters $\varpi_{1}$ and $\varpi_{2}$.

4 Analytical secular model

In order to construct an analtical model from the equations of motion (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 ${\bf R_{i}}$, the positions of the masses in Jacobi coordinates are given by:

[TABLE]

where

[TABLE]

Analogous expressions relate the velocities vectors in both reference systems.

4.1 Secular evolution of the planetary spin

Expanding the position and velocity vectors in equation (16) up to second order in $\alpha=a_{1}/a_{2}$ and the eccentricities, and averaging with respect to both mean longitudes, we finally obtain the rate of change of the rotational frequency of the planet as:

[TABLE]

where the non-zero coefficients different are given by

[TABLE]

with

[TABLE]

The stationary spin rate \big{<}\Omega_{2}\big{>}_{\rm stat} predicted by this equation can be easily calculated by equating expression (19) to zero. The explicit form of the equilibrium rotational frequency was found to be

[TABLE]

In the limit case in which the mass of one of the stars reduces to zero we recover the classical 2-body super-synchronous stationary solution \big{<}\Omega_{2}\big{>}_{\rm 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_{1},e_{2}\to 0$), the CB planet stationary solution is sub-synchronous by a factor that decreases proportional to $\alpha^{2}$ as we move outward from the binary, and is maximum for equal-mass stars $m_{1}=m_{0}$.

Figure 6 shows two color plots with the value of \big{<}\Omega_{2}\big{>}_{\rm 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_{1}/m_{0}\lesssim 0.1$, the estimated value of \big{<}\Omega_{2}\big{>}_{\rm 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 2-body problem, i.e. \big{<}\Omega_{2}\big{>}_{\rm 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 2-body solution and lead to appreciable differences in the rotational dynamics.

The dependence of \big{<}\Omega_{2}\big{>}_{\rm stat} with the eccentricities is analyzed in the bottom frame of Figure 6. We note that the sub-synchronous equilibrium state is only observed for low eccentricities of the planet, typically $e_{2}\lesssim 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 \big{<}\Omega_{2}\big{>}_{\rm stat}<n_{2} should probably 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, Figure 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^{5}$ years, and our model seems to reproduce the equilibrium behavior 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 $\Omega_{2}$ also suffers periodic oscillations.

Finally, as can be observed from equation (22), the stationary spin solution for the CB planets is not a function of the planetary mass $m_{2}$ 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}\simeq 118$ days, about a $12\%$ higher than the one predicted by the 2-body synchronization model.

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_{2}$ and eccentricity $e_{2}$. 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

[TABLE]

where $\delta{\bf{f}_{2}}$ is the total tidal force (per unit mass) affecting the 2-body motion of the planet around the center of mass of $m_{0}$ and $m_{1}$, and has the form:

[TABLE]

Substituting equation (9) in order to express the total force in terms of the individual two-body tidal interactions, we obtain

[TABLE]

where

[TABLE]

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 $\delta{\bf{f}}_{2}$ may be written in terms of differences of the type $({\bf{F}_{2,j}}-{\bf{F}_{j,2}})$, where $j=0,1$. From equations (10), each of these differences may be explicitly written as

[TABLE]

where we have defined

[TABLE]

and a new “averaged” rotational frequency

[TABLE]

Notice that expression (27) has the same functional form as the tidal force in the 2-body problem (eq. 10) with a magnitude given by ${\cal K}_{j}^{(+)}$ and a rotational frequency defined by ${\bf{\bar{\Omega}}^{(j)}_{2}}$. In the limit where $m_{1}\rightarrow 0$ and ${\cal R}_{1}\rightarrow 0$, the term in the tidal force associated to ${\cal 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 ${\bf\Delta_{2j}}=\boldsymbol{\rho}_{2}+\gamma_{j}\boldsymbol{\rho}_{1}$, substituting in the Gauss equation (23), expanding in power series of $\alpha$, $e_{1}$ and $e_{2}$ and, finally, averaging over the mean longitudes, we obtain:

[TABLE]

where the non-zero coefficients are explicitly given by

[TABLE]

Figure 8 shows the normalized value of \big{<}da_{2}/dt\big{>} in the $(\alpha,m_{1}/m_{0})$ plane for three different values of the binary and planet eccentricities. For each value of $m_{1}$ the physical radius of the star was modified following the empirical rule ${\cal 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 colored in tones of blue, while those leading to a secular increase of $a_{2}$ 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_{1}/m_{0}$ above which the tidal interaction of the system leads to an outward migration of the planet. The critical value of $m_{1}$ appears to be larger for more eccentric binaries and lower for stars in almost circular orbits. As expected, as $m_{1}\rightarrow 0$ the migration is inwards, 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-order of $\alpha$. Consequently, the results shown here and in Figure 6 are not expected to be accurate (or even valid) for $\alpha\rightarrow 1$. We have nevertheless opted to include the complete range solely for illustrative purposes.

The time variation of the eccentricity $e_{2}$ may be found from the orbital angular momentum ${\bf{L}}_{2}$ in the Jacobi reference frame. In the planar case, we have

[TABLE]

whose time derivative due to tidal forces leads to

[TABLE]

Extracting the eccentricity term, we finally obtain:

[TABLE]

Introducing elliptic expansions in a similar manner as done for (30), and averaging over short-period terms, we obtain:

[TABLE]

where now the non-zero coefficients are given by

[TABLE]

Contrary to $da_{2}/dt$, we found that the eccentricity of the planet is always damped, at least for the initial conditions and system parameters tested here.

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 equations (30) and (35) with the numerical results obtained using the original unaveraged equations (23) and (34). We consider the nominal system parameters detailed in Table 1 but varied the planetary eccentricity and semimajor axis ratio $\alpha$. For each we computed $da_{2}/dt$ and $de_{2}/dt$ as a function of the reduced mass

[TABLE]

by varying $m_{1}$. Due to the rapid rotational synchronization timescales, we consider stationary spins for the stars and for the planet according to equation (22).

Results are shown in Figure (9). In all the panels the colors represent different planetary eccentricities ($e_{2}=0.01$ in blue, $e_{2}=0.05$ in green and $e_{2}=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 $\alpha$: the reference value in the bottom panels ($\alpha=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_{2}$ decays more rapidly 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_{2}\sim 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}_{crit}$ associated to the limit between inward and outward migration is very well reproduced.

Finally, Figure (10) shows the dependence of $\tilde{\mu}_{crit}$ as function of $\alpha$ 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 $\alpha$. 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 $\alpha\sim 0.3$, the truncated version shows a much smaller region of validity, reduced down to $\alpha\sim 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}_{crit}$ increases for smaller values of $\alpha$. In the limit when $\alpha\rightarrow 0$, we expect the system to behave as a planet orbiting a single star of mass $m_{0}+m_{1}$ and all initial conditions should lead to an inward migration of the semimajor axis.

5 Summary and discussion

In this work we 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 2-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 mean-motion resonances between $m_{1}$ and $m_{2}$. Thus, in the non-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. 2. 2.

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 the results 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 $\sim$ 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_{1}$ and eccentricity $e_{1}$.

The long-term orbital evolution of the planet is curiously different: as a result of the tidal interaction the planet migrates outward 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 $\Delta\varpi$ 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 $\alpha$ and the eccentricities $e_{1}$ and $e_{2}$. 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 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 $\alpha$. 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 outward. 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.

Acknowledgements.

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).