Revisiting the dynamics of planets in binaries: evolutionary timescales and the effect of early stellar evolution

TL;DR

This paper develops advanced models for the long-term dynamics of planets in binary systems, incorporating stellar evolution effects, and demonstrates their importance through specific case studies and an open-source computational tool.

Contribution

It introduces explicit formulas for the secularized octupolar Hamiltonian with relativistic and non-conservative effects, and applies them to study stellar evolution impacts on planetary dynamics.

Findings

Octupolar expansion is essential for accurate evolution predictions of HD 80606b.

Stellar radius changes significantly affect planetary dynamical timescales.

The open-source package SecDev3B enables reproducible and improved dynamical modeling.

Abstract

The discovery of planets in binaries is one the most interesting outcomes of planetary research. With the growing number of discoveries has also grown the interest on describing their formation, long-term evolution and potential habitability. In this work we revisit the dynamics of planets in S-type binary systems. For that purpose we develop explicit formulas for the secularized octupolar Hamiltonian, coupled with general relativistic corrections and non-conservative interactions. We implemented those formulas in an open-source package \texttt{SecDev3B}, that can be used to reproduce our results or test improved versions of the models. In order to test it, we study the long-term dynamical evolution of S-type binary planets during the pre-main-sequence phase of stellar evolution. During that phase, stellar radius significantly changes in timescales similar to secular timescales. We hypothesize that when close-encounters between the planet and its host star happens (e.g. via Lidov-Kozai effect), particularities in the secular formalism plus changes in stellar radius may alter significantly the dynamical evolution. We study the well-known binary planet HD 80606b and found that an octupolar expansion of the conservative Hamiltonian is required to properly predict its dynamical evolution. We also apply the dynamical model, enriched with results coming from stellar evolutionary models, to demonstrate that in S-type systems around low-mass stars, with relative high inclinations ($\itot\ge 60^\circ$), moderate eccentricities ($0.2\le e\le 0.4$) and planets located around 1 AU, the evolution of stellar radius during the first few hundreds of Myr, alters significantly the timescales of dynamical evolution.

Tables

Table 1. Table 1: Initial conditions of each experiment described in this work. Quantities marked with (-) takes a series of values depending on the grid properties described in Section 5.3 .

Property	Experiment 1	Experiment 2	Experiment 3	Experiment 4
Type	S-type	S-type	S-type	S-type
Protoype	HD 80606b	Fictitious	Fictitious	Fictitious
$a_1;a_2$ (AU)	5; 1000	$2.5;\mathrm{\hspace{0.17em}300}$	-;  50	1;  50
$e_1;e_2$	0.1;  0.5	$0.6;\mathrm{\hspace{0.17em}0.3}$	-;  0.3	0.35;  0.3
$i_{\mathrm{tot}}$ (deg)	85.6	$80$	-	80
$g_1;g_2$ (deg)	45;  0	$0;\mathrm{\hspace{0.17em}0}$	0;  0	0;  0
$m_1;m_2;m_3$ ($M_{\odot }$)	1.1;  0.00744;  1.1	$1;\mathrm{\hspace{0.17em}0.001};\mathrm{\hspace{0.17em}0.667}$	0.6;  0.01;  0.6	0.6;  0.01;  0.6
$R_1;R_2$ ($R_{\odot }$)	$1;\mathrm{\hspace{0.17em}0.1}$	$1;\mathrm{\hspace{0.17em}0.1}$	0.59;  0.1	0.59;  0.1
$P_{\mathrm{rot},A};P_{\mathrm{rot},B}$ (day)	$20;\mathrm{\hspace{0.17em}0.416}$	$20;\mathrm{\hspace{0.17em}1}$	20;  1	20;  1
$t_{v,A};t_{v,B}$ (years)	$50;\mathrm{\hspace{0.17em}0.001}$	$50;\mathrm{\hspace{0.17em}0.001}$	50;  0.001	50;  0.001
$k_s;k_p$	$0.014;\mathrm{\hspace{0.17em}0.255}$	$0.014;\mathrm{\hspace{0.17em}0.255}$	0.014;  0.255	0.014;  0.255
${\kappa }_A;{\kappa }_B$	$0.08;\mathrm{\hspace{0.17em}0.25}$	$0.08;\mathrm{\hspace{0.17em}0.25}$	0.08;  0.25	0.08;  0.25

Equations

$$\begin{split}\displaystyle\mathcal{H}\equiv&\mathcal{H}_{\mathrm{kep},1}+\mathcal{H}_{\mathrm{kep},2}+\mathcal{H}_{\mathrm{int}}\\ =&\frac{\mathcal{G}m_{1}m_{2}}{2a_{1}}+\frac{\mathcal{G}(m_{1}+m_{2})m_{3}}{2a_{2}}+\\ &+\frac{\mathcal{G}}{a_{2}}\sum_{n=2}^{\infty}\alpha^{n}M_{n}\bigg{(}\frac{r_{1}}{a_{1}}\bigg{)}^{n}\bigg{(}\frac{a_{2}}{r_{2}}\bigg{)}^{n+1}P_{n}(\cos\Phi)\end{split}$$

$$l_j$$

$$h_j$$

$$g_j$$

$$\begin{split}\displaystyle\mathcal{H}=&\frac{\beta_{1}}{2L_{1}^{2}}+\frac{\beta_{2}}{2L_{2}^{2}}+8\beta_{3}\frac{L_{1}^{4}}{L_{2}^{6}}\bigg{(}\frac{r_{1}}{a_{1}}\bigg{)}^{2}\bigg{(}\frac{a_{2}}{r_{2}}\bigg{)}^{3}(3\cos^{2}\Phi-1)\\ &+2\beta_{4}\frac{L_{1}^{6}}{L_{2}^{8}}\bigg{(}\frac{r_{1}}{a_{1}}\bigg{)}^{3}\bigg{(}\frac{a_{2}}{r_{2}}\bigg{)}^{4}(5\cos^{3}\Phi-3\cos\Phi)\end{split}$$

$$r_{j,\mathrm{orb}}=\frac{a_j(1-e_j^2)}{1+e_j\cos f_j}$$

$${\mathbf{r}}_{j,\mathrm{in}}=\mathbf{M}({\Omega }_j,I_j,{\omega }_j){\mathbf{r}}_{j,\mathrm{orb}}$$

$$\begin{array}{rl}\cos \Phi & ={\hat{\mathbf{r}}}_{2,\mathrm{in}}\cdot {\hat{\mathbf{r}}}_{1,\mathrm{in}}\\ & ={\hat{\mathbf{r}}}_{1,\mathrm{orb}}^{⊺}{\mathbf{M}}^{⊺}(h_1,I_1,g_1)\mathbf{M}(h_2,I_2,g_2){\hat{\mathbf{r}}}_{2,\mathrm{orb}}\end{array}$$

$$\begin{split}\displaystyle\cos\Phi=&\sin(f_{2}+g_{2})\big{[}A_{1}\sin(f_{1}+g_{1})+A_{3}\cos(f_{1}+g_{1})\big{]}-\\ &\cos(f_{2}+g_{2})\big{[}A_{2}\sin(f_{1}+g_{1})-A_{4}\cos(f_{1}+g_{1})\big{]}\end{split}$$

$$⟨⟨\mathcal{H}⟩⟩=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }{\int }_0^{2\pi }d{l}_{1}d{l}_2{\mathcal{H}}_{\mathrm{int}}$$

$$\begin{split}\displaystyle\langle\langle\mathcal{H}\rangle\rangle\approx&\frac{2\beta_{3}}{\big{(}1-e_{2}^{2}\big{)}^{3/2}}\frac{L_{1}^{4}}{L_{2}^{6}}\bigg{[}-2(2+3e_{1}^{2})+3\big{(}1-e_{1}^{2}\big{)}\big{(}K_{1}^{2}+K_{4}^{2}\big{)}+\\ &+3\big{(}1+4e_{1}^{2}\big{)}\big{(}K_{2}^{2}+K_{3}^{2}\big{)}\bigg{]}+\frac{15\beta_{4}}{16\big{(}1-e_{2}^{2}\big{)}^{5/2}}\frac{L_{1}^{6}}{L_{2}^{8}}\times\\ &\times\Bigg{\{}e_{1}e_{2}\bigg{[}4K_{2}\big{(}4+3e_{1}^{2}\big{)}-5\big{(}3+4e_{1}^{2}\big{)}K_{2}\big{(}K_{2}^{2}+K_{3}^{2}\big{)}-\\ &5\big{(}1-e_{1}^{2}\big{)}\big{(}2K_{1}K_{3}K_{4}+3K_{1}^{2}K_{2}+K_{4}^{2}K_{2}\big{)}\bigg{]}\Bigg{\}}\end{split}$$

$$\begin{split}\displaystyle&\langle\langle\mathcal{H}\rangle\rangle=\frac{\beta_{3}}{2\big{(}1-e_{2}^{2}\big{)}^{3/2}}\frac{L_{1}^{4}}{L_{2}^{6}}\Bigg{\{}\big{(}2+3e_{1}^{2}\big{)}\big{(}1+3\cos 2i_{\mathrm{tot}}\big{)}+\\ &+30e_{1}^{2}\cos 2g_{1}\sin^{2}i_{\mathrm{tot}}\Bigg{\}}+\frac{15\beta_{4}e_{1}e_{2}}{64\big{(}1-e_{2}^{2}\big{)}^{5/2}}\frac{L_{1}^{6}}{L_{2}^{8}}\times\\ &\times\Bigg{\{}70e_{1}^{2}\cos 2g_{1}\sin^{2}i_{\mathrm{tot}}\big{(}\cos g1\cos g2+\cos i_{\mathrm{tot}}\sin g_{1}\sin g_{2}\big{)}+\\ &+\bigg{(}6-13e_{1}^{2}+5\big{(}2+5e_{1}^{2}\big{)}\cos 2i_{\mathrm{tot}}\bigg{)}\cos g_{1}\cos g_{2}+\\ &+\bigg{(}5\big{(}6+e_{1}^{2}\big{)}\cos 2i_{\mathrm{tot}}-7\big{(}2-e_{1}^{2}\big{)}\bigg{)}\cos i_{\mathrm{tot}}\sin g_{1}\sin g_{2}\Bigg{\}}\end{split}$$

$$\frac{d{L}_j}{dt}$$

$$\frac{d{G}_j}{dt}$$

$$\frac{d{H}_j}{dt}$$

$$\ddot{\vec{r}}=-\frac{\mathcal{G}(m_1+m_2)}{r_1^3}{\vec{r}}_1+{\vec{a}}_{\mathrm{TF}}+{\vec{a}}_{\mathrm{QD}}$$

$${\vec{a}}_{\mathrm{TF}}$$

$${\vec{a}}_{\mathrm{QD}}$$

$${\vec{J}}_1=\sqrt{\mathcal{G}(m_1+m_2)a_1(1-e_1^2)}{\hat{z}}_1,$$

$$\displaystyle C=\frac{9}{2}\bigg{(}\frac{1+2k_{s}}{2k_{s}m_{1}R_{s}}\bigg{)}^{2}m_{2}(m_{1}+m_{2})$$

$$\displaystyle t_{F,s}=\frac{t_{V,s}}{9}\bigg{(}\frac{a_{1}}{R_{s}}\bigg{)}^{8}\frac{m_{1}^{2}}{(m_{1}+m_{2})m_{2}}\frac{1}{(1+2k_{s})^{2}}$$

$${\mathcal{H}}_{\mathrm{GR}}={\gamma }_1{\vec{P}}_1^4+{\gamma }_2\frac{{\vec{P}}_1^2}{r_1}+{\gamma }_3\frac{({\vec{r}}_1\cdot {\vec{P}}_1{)}^2}{r_1^3}+{\gamma }_4\frac{1}{r_1^2}$$

$${\vec{P}}_1={\vec{v}}_1+\frac{1}{c^2}\left(4{\gamma }_1({\vec{v}}_1\cdot {\vec{v}}_1){\vec{v}}_1+\frac{2{\gamma }_2}{r_1}{\vec{v}}_1+\frac{2{\gamma }_4}{r_1^3}({\vec{r}}_1\cdot {\vec{v}}_1){\vec{r}}_1\right)$$

$$\begin{split}\displaystyle\langle\langle\mathcal{H}_{\mathrm{GR}}\rangle\rangle=&\gamma_{1}\langle\langle\vec{v}_{1}^{4}\rangle\rangle+\gamma_{2}\bigg{\langle}\bigg{\langle}\frac{\vec{v}_{1}^{2}}{r_{1}}\bigg{\rangle}\bigg{\rangle}+\gamma_{3}\bigg{\langle}\bigg{\langle}\frac{(\vec{r}_{1}\cdot\vec{v}_{1})^{2}}{r_{1}^{3}}\bigg{\rangle}\bigg{\rangle}\\ &+\gamma_{4}\bigg{\langle}\bigg{\langle}\frac{1}{r_{1}^{2}}\bigg{\rangle}\bigg{\rangle}\end{split}$$

$$⟨⟨{\mathcal{H}}_{\mathrm{GR}}⟩⟩=-\frac{3{\mu }^4{\beta }^5}{c^2{L}_1^3{G}_1}$$

$$\displaystyle\bigg{(}\frac{dg_{1}}{dt}\bigg{)}_{\mathrm{GR}}=-\frac{\partial\langle\langle\mathcal{H}_{GR}\rangle\rangle}{\partial G_{1}}=\frac{3\mu^{4}\beta^{5}}{c^{2}L_{1}^{3}G_{1}^{2}}$$

$$\displaystyle\bigg{(}\frac{dg_{1}}{dt}\bigg{)}_{\mathrm{GR}}=\frac{3\mu^{3/2}}{c^{2}a_{1}^{5/2}(1-e_{1}^{2})}$$

$${\tau }_{\mathrm{th}}=\frac{\mathcal{G}{m}_s^2}{R_s{L}_s}$$

$$a_1(t=t_{a,1/2})=\frac{1}{2}{a}_1(t=0)$$

$$R_{L1}\approx 0.44a_{\mathrm{bin}}\frac{q^{0.33}}{(1+q{)}^{0.2}}$$

$$\displaystyle\bigg{(}\frac{da_{\mathrm{bin}}}{dt}\bigg{)}_{\mathrm{no-acc}}=-\frac{\dot{m}_{s1}}{m_{s1}+m_{s2}}a_{\mathrm{bin}}$$

$$L_1$$

$$L_2$$

Full text

††pubyear: 2018††pagerange: Revisiting the dynamics of planets in binaries: evolutionary timescales and the effect of early stellar evolution–A.3

1 Introduction

Planets in binaries are among the most interesting and unexpected discoveries in the era of exoplanets. These triple (or multiple) gravitational systems are divided, according to their hierarchical configuration, into two categories (Dvorak, 1986): 1) S-type systems, where the planet orbits one member of the stellar binary (S stands for satellite) and the other star is an outer perturber; and 2) P-type systems in which both stars are closely orbiting each other with the planet moving on a wider orbit (P stands for planet).

To the date of writing, 132 planets have been discovered around 91 binaries (Schwarz et al., 2016)111See for instance the Catalogue of Exoplanets in Binary Stars Systems https://www.univie.ac.at/adg/schwarz/multiple.html. Around 75% of them correspond to S-type systems.

Planets in binaries are not only interesting because their dynamical properties (which we explore in this work) but also because they challenge our current understanding of planetary formation.

In the case of isolated stars, the formation of planets is the result of two possibles mechanism: core accretion or disk instabilities (Armitage, 2009). The same mechanisms are in action in S-type systems, provided the separation of the outer companion is large enough to make its gravitational effect on the circumstellar disk, negligible. However, since we have discovered planets around compact S-type systems, namely, binaries with separations as small as $20-25$ AU (Hatzes et al. 2003; Zucker et al. 2004; Mugrauer & Neuhäuser 2005), the outer companion is expected to have had a considerable effect on the coalescence of planetesimals. Still, the existence of planets around this type of systems, evidences the fact that binarity does not preclude the subsequent formation of planets (Haghighipour, 2010).

Independently of the mechanism that lead to their formation, planets in binaries appears to be more common than expected. The apparent lack of them in what seems to be planet-less binary systems, could be explained due to the biases of the discovering techniques (Pelupessy & Portegies Zwart, 2013).

Numerous works have studied the dynamics of circumbinary systems, both using direct $N$-body calculations (see e.g. Holman & Wiegert 1999 and Kostov et al. 2016) and secular theories based on quadrupolar and octupolar order expansions of the perturbative Hamiltonian (see e.g. Wu & Murray 2003; Mardling 2013; Naoz et al. 2013; Correia et al. 2016). The secular approach has proved to be significantly more efficient than direct $N$-body methods for long integration times (Naoz, 2016). One of the major advantages of secular theories relies in the fact that relatively large time-steps (even greater than the orbital periods of the components) can be used in the computations. This is possible thanks to the fact that rapidly varying quantities (e.g. mean anomalies) are properly removed from the Hamiltonian and hence from the equations of motion, by averaging-out over relevant time-scales. However, a major disadvantage of this approach is the need for the inclusions of octupolar order terms when both members of the inner binary are not equal in mass.

Mathematical expressions in this formalism are complex. This fact motivated us to develop a software package SecDev3B that implements all those formulas, while adding new functionalities (see below). This package will contribute to the community making our models and results (1) reproducible and (2) improvable. Its modular design and internal documentation will allow the future inclusion of new and better approximations.

But this paper is not only about the software, but also on what we discovered while developing and testing it. We report, for instance, a major difference in the evolution timescale of the well-known HD 80606b planet when the integration is made retaining the octupolar term of the secular Hamiltonian, in contraposition to those previously studied cases in which a quadrupolar expansion was used.

When integrations of hundreds of million years (Myr) are performed on this type of systems, there could be noticeable changes in the radius and internal properties of the stars and the planet. This is particularly notorious during the pre and post main sequence evolution of the stellar components (Carroll & Ostlie, 2007) and during the first hundreds of Myr following planetary accretion, especially in the case of giant planets (Fortney et al., 2007).

In the past, several works have explored the interplay between dynamics, stellar evolution and other aspects of stellar structure, especially for the case of triple stars (see e.g. Shappee & Thompson 2013; Toonen et al. 2016). Similar effects have been studied in the case of planets in binaries using $N$-body simulations (Kratter & Perets, 2012). These works, however, have focused on the impact that wind-driven mass loss may have on the dynamic of these systems.

Before achieving hydrostatic equilibrium (i.e. zero main sequence or ZAMS for short), stars experience a relatively slow contraction. Depending on the stellar mass, the radius of a pre-main-sequence star may vary by a factor of up to 3 while “traveling” along the Hayashi track (Harpaz, 1994). The time of these changes are of the order of the thermal timescale, which for the case of a $0.6\,M_{\odot}$ star is about $220$ Myr. Other changes in the stellar structure, relevant to their interaction with potential companions, also occur during these phases. Thus for instance, the internal mass distribution and its response to tidal effects, which affect the interchange of angular momentum, also change while the star is contracting (see e.g. Zuluaga et al. 2016). On the dynamical side and depending on the relative orbital inclination of the binary and the planet (and to a less extent, on the initial values of inner eccentricty and periastron argument, see e.g. Bataille et al. 2018), high-amplitude excitations of the planetary eccentricity will happen (Valtonen & Karttunen, 2006). During these eccentricity excursions, the planet may reach temporarily small periapse distances and strong tidal interaction with the host star are triggered. Since tidal torques are very sensitive to stellar and planetary radius (Murray & Dermott, 1999), the outcome of the dynamical evolution may be significantly different than in the case when we assume that those radius are constant or vary in time.

We hypothesize that including changes in stellar and planetary radius into the dynamics will significantly impact the outcome of the early dynamical evolution of planets in binaries, especially in S-type systems. Since the orbital “architecture” of these systems, their survival and stability, depend on what happen at the beginning of their life, including stellar and planetary evolution will have an impact on our interpretation of the observations. This paper is also aimed to (partially) test this hypothesis. For the sake of conciseness, we will focus here on testing the effect of changes in stellar radius and left the effects of the evolution of stellar interior structure and planetary radius, for a future paper.

As we will show, we can identify a limited region in the space of inner semi-major axis and eccentricity, where the inclusion of changes in radius alters significantly the timescales of evolution.

This paper is structured as follows. In Section 3 we present the conservative, non-conservative and relativistic formalisms that constitute the building blocks of the complete secular theory. Section 4 presents some relevant aspects of the pre-main-sequence stellar evolution with an special focus on the timescales involved. In Section 5 we show and discuss the results of four numerical experiments we performed to test our hypothesis. Section 6 is intended to discuss the impact that other effects could have in the secular dynamics of these systems, the limitations and prospect of our models.

2 Secular theory of hierarchical three-body systems

The secular theory of hierarchical three-body systems have been widely developed and presented in literature (see eg. Harrington 1968; Ford et al. 2000; Naoz et al. 2013; Carvalho et al. 2016; Correia et al. 2016 and references there in). For the sake of completeness, we reproduce in this section the most relevant results of the theory, paying special attention on presenting explicit expressions of the relevant formulas. Many of these formulas are absent in literature, due (partially) to the complexity of the mathematical involved expressions. In an effort to make all our results fully reproducible, we also refer the reader to Section A for further details of the mathematical expressions.

3 Secular theory of hierarchical three-body systems

The secular theory of hierarchical three-body systems have been widely developed and presented in literature (see e.g. Harrington 1968; Ford et al. 2000; Naoz et al. 2013; Carvalho et al. 2016; Correia et al. 2016 and references there in). For the sake of completeness, we reproduce in this section the most relevant results of the theory, paying special attention on presenting explicit expressions of the relevant formulas. Many of these formulas are absent in literature, due (partially) to the complexity of the mathematical involved expressions. In an effort to make all our results fully reproducible, we also refer the reader to Appendix A for further details of the mathematical expressions.

3.1 Conservative motion

In the absence of dissipative forces and relativistic effects, the Hamiltonian for the orbital motion of a general three body system is (Harrington, 1968)222We have adopted the convention of defining the Hamiltonian as the negative of the mechanical energy to agree with the notation of previous works (e.g. Ford et al. 2000; Naoz et al. 2013; Carvalho et al. 2016).:

[TABLE]

Here $M_{n}=m_{1}m_{2}m_{3}[m_{1}^{n-1}-(-m_{2})^{n-1}]/(m_{1}+m_{2})^{n}$, $P_{n}$ is the Legendre polynomial of degree $n$ and $\Phi$ is the angle between the two Jacobian coordinates $\vec{r}_{1}\equiv\mathbf{\vec{r}_{1}}$, $\vec{r}_{2}\equiv\mathbf{\vec{r}_{2}}$ (see Figure 1) and $\alpha\equiv a_{1}/a_{2}$ is the ratio of the semi major axis of the inner (star 1 and planet) and the outer system (star 1-planet and star 2). Hereafter, quantities with subindex $1$ and $2$ will be referred to the inner and outer orbits respectively. To avoid confusion with other dynamical quantities, the gravitational constant will be written as $\mathcal{G}$.

The first two terms of the Hamiltonian ($\mathcal{H}_{\mathrm{kep},1}$, $\mathcal{H}_{\mathrm{kep},2}$) correspond to the contribution of the relative Keplerian motions of $m_{2}$ around $m_{1}$, and that of $m_{3}$ around the center of mass of $m_{1}$ and $m_{2}$, respectively. The third term, $\mathcal{H}_{\mathrm{int}}$ correspond to the interaction Hamiltonian, which accounts for the perturbative effect of the outer orbit on the inner one. The strength of the latter term is proportional to the relative size $\alpha$ of the inner and outer orbits. If $\alpha$ is small enough, namely $\alpha<1/2$, secular perturbation techniques can be used (Correia et al., 2016). In this hierarchical regime, high order terms in the summation can be neglected. If the only term retained is that corresponding to $n=2$, the approximation is made to quadrupole level. When summation is expanded up to $n=3$, we have the octupole approximation.

Since the Hamiltonian does not depends explicitly on time, we can use Hamilton-Jacobi theory to find a better set of canonical variables (Goldstein et al., 2002). We will denote the new generalized coordinates as $l_{j},h_{j},g_{j}$ (where $j=1,2$ are indexes for the inner and outer system) and they are a function only of the initial conditions and time. The corresponding conjugated momenta will be $L_{j},H_{j},G_{j}$, and they are all constants of motion.

A particularly useful set of canonical coordinates for this system, fulfilling the previous conditions, are the Dalaunay elements (Valtonen & Karttunen, 2006):

[TABLE]

where $\mu_{j}={\cal G}{\cal M}_{j}$, ${\cal M}_{j}$ is the total mass of each of the systems (${\cal M}_{1}=m_{1}+m_{2}$, ${\cal M}_{2}={\cal M}_{1}+m_{3}$), $m_{r,j}$ are the reduced masses: $m_{r,1}=m_{1}m_{2}/\mathcal{M}_{1}$, $m_{r,2}=m_{3}\mathcal{M}_{1}/\mathcal{M}_{2}$ and $M_{j},\Omega_{j},\omega_{j}$ are the mean anomalies, the longitudes of the ascending nodes and the arguments of the periapse with respect to an inertial whose fundamental plane is the Laplace plane of the three body system.

In terms of these elements, the octupolar Hamiltonian is now written:

[TABLE]

with $\beta_{i}$ being auxiliar functions of the particle masses, whose explicit expressions are presented in Appendix A.

All the varying quantities in the Hamiltonian, namely $\Phi$, $r_{1}$ and $r_{2}$, implicitly depend the generalized-mean anomalies $l_{j}$ (see equation 2) which are related to the true anomalies $f_{j}$.

In the “natural” (non-inertial) reference frame of each orbit ($x_{j}$ axis pointing to the periapse and $z_{j}$ axis being perpendicular to the orbital plane), the relative position of each system is simply (Murray & Dermott, 1999): $\mathbf{r}_{j,\mathrm{orb}}=r_{j,\mathrm{orb}}\mathbf{\hat{r}}_{j,\mathrm{orb}}$, with $\mathbf{\hat{r}}_{j,\mathrm{orb}}=(\cos f_{j},\sin f_{j},0)$ and:

[TABLE]

The position vector of each body in the inertial system will be:

[TABLE]

with $\mathbf{M}(\Omega_{j},I_{j},w_{j})=\mathbf{M}(h_{j},I_{j},g_{j})$ an Eulerian $3\times 3$ rotation matrix ($I_{j}$ is the inclination of each orbit relative to the Laplace plane). Here, subscripts “orb" and “in" emphasizes that the quantity is measured in the orbital or inertial reference system, respectively. In the following derivation $\mathbf{M}^{\intercal}$ will be the tranpose matrix.

Using this conventions, the relative angle $\Phi$ will be:

[TABLE]

or explicitly, in terms of the Delaunay elements:

[TABLE]

where we have introduced auxiliary functions $A_{i}(h_{j},I_{j},g_{j})$ ($i=1,2,3,4$), whose explicit expressions are in Appendix A.

When equations (6) and (9) are replaced in (5), the octupolar Hamiltonian becomes an explicit function of the true anomalies of both orbits.

To construct a secular theory, we need to get rid from the Hamiltonian those rapidly oscillating variables. This is done with aid of the Von Zeipel transformation (Marchal, 2012), namely, if the system does not exhibit mean motion resonances, the dynamics could be described with the doubly-averaged interaction Hamiltonian:

[TABLE]

The resulting doubly-averaged Hamiltonian expanded to octupolar order reads333It is interesting to notice that the integration on the inner orbit is easier to be done if we instead of the true anomaly $l_{1}$, use the eccentric anomaly as the integration variable.:

[TABLE]

where explicit expressions for the auxiliary functions $K_{i}(h_{j},I_{j},g_{j})$ ($i=1,2,3,4$) are also provided in Appendix A. This is the first time that an analytical formula for this Hamiltonian without any simplification is provided explicitly in literature.

This expression can be greatly simplified if we have in mind that when working in the Laplace plane $\Delta h=h_{1}-h_{2}=180^{\circ}$ (Valtonen & Karttunen, 2006) a condition called “elimination of the nodes”. We need to be very cautious with this step. As discussed in Naoz et al. (2013), this procedure is not correct because by applying it, we are reducing the total volume of the phase space over which are we will construct the variational trajectories to derive the equations of motion through Hamilton’s principle. The elimination of the nodes directly in the Hamiltonian leads to erroneous conclusions like the independent conservation of the $z-$component of the angular momentum of both orbits: $H_{j}=G_{j}\cos I_{j}$ (Kozai 1962). To derive the equations of motion for inner and outer inclinations, we rely in the geometric properties arising from the conservation of the total angular momentum as explained in Appendix A. Fortunately, for the orbital elements whose defining equation of motion does not depends on the partial differential operator with respect to $h_{i}$, we can use the simplified Hamiltonian. Therefore, once this caution has been made, we can express the “node-eliminated" form of the secular Hamiltonian as:

[TABLE]

Once the secular Hamiltonian has been derived, the secular equations of motion are given by the canonical Hamilton’s equations:

[TABLE]

All this procedure produces a set of 8 (non-trivial) first order, nonlinear and coupled differential equations. The explicit expression of the equation of motion can be found in Appendix A.1 where we have followed the same notation of Naoz et al. (2013).

3.2 Non-conservative motion

As orbital evolution develops, the Lidov-Kozai effect generates high-amplitude oscillations of the inner eccentricity and mutual inclination (Lidov, 1962; Kozai, 1962). Depending on initial conditions, the eccentricity values reached can be high enough to produce close encounters at periastron passages, which in turn may trigger strong tidal interactions between the inner bodies.

Tides have a direct impact on the orbital elements of the inner system and on the rotational properties of $m_{1}$ and $m_{2}$. On the other hand, since we are interested only on those systems which are “hierarchically” stable, i.e. systems where the hierarchy of the mean distances are satisfied along the evolution, we can safely neglect the tidal effects experimented by the outer body.

Recently, many works have addressed the dynamics of tidal interaction with the aid of sophisticated rheological models (e.g Efroimsky & Williams 2009; Ferraz-Mello 2013; Correia et al. 2014). Others have applied those models to the context of exoplanets dynamics (e.g Makarov et al. 2012; Cuartas-Restrepo et al. 2016). In this work, and for the sake of simplicity, we implemented the classical tidal formalism of Hut (1981) and Eggleton et al. (1998).

In the kind of S-type systems we are studying here, the planet ($m_{2}$) raises a tidal bulge on its host star ($m_{1}$). Moreover, since the star is rotating (very fast at early ages), it also has a rotational distortion (that for simplicity we will express only to quadrupolar order). The total two-body acceleration will be then given by:

[TABLE]

where the induced non-gravitational tidal $\vec{a}_{\mathrm{TF}}$ and quadrupolar-deformation $\vec{a}_{\mathrm{QD}}$ accelerations are given by (Eggleton et al., 1998):

[TABLE]

where $m_{r}$ is the reduced mass of the inner system and $\vec{\Omega}_{s}$ is the rotational velocity vector of the star ($m_{1}$) whose orientation is completely arbitrary in this formalism.

Hereafter, quantities labeled with subscripts $s$ and $p$ correspond to properties of the star and the planet, respectively.

The tidal component of the acceleration in Equation (16) depends on the angular momentum of the inner system:

[TABLE]

on the viscous timescale $t_{\mathrm{V1}}$ (which is assumed constant) and a rheological constant $C$, which is defined as:

[TABLE]

where $k_{s}$ is the apsidal motion constant (a measure of the quadrupolar deformability). This constant is related to the Love number of the star by $k_{L,s}=2k_{s}$ (Naoz, 2016).

From Equation (17), it is clear that tidal acceleration has a strong dependence on the instantaneous separation of the inner bodies $r_{1}$, $a_{\rm TF}\propto r_{1}^{-10}$, which make its effect very relevant at short distances.

On the other hand the quadrupolar deformation contribution444It is worth noticing that the quadrupolar dissipative force arises in the regime when the mechanical energy is dissipated at a rate proportional to the square of the rate of change of the quadrupolar moment tensor(Eggleton et al., 1998). is proportional to the factor $A_{s}=2k_{s}R_{s}^{5}$ which also depends strongly on the stellar radius.

After those forces are included in the equations of motion for $m_{1}$ and $m_{2}$, the next step is to derive the averaged (secularized) effect of them on the dynamical evolution of the inner system. For further details of the secularization processes we refer the reader to the work of Eggleton et al. (1998).

Here we just summarized the contribution of these dissipative forces to the inner secular dynamics. The complete set of secular equations that rules the orbital and rotational evolution due to the dissipative effects are given by explicit expressions in Appendix A.

The magnitude of the tidal interaction on the orbital dynamics of the inner system, is better illustrated using the so-called tidal friction timescale for the star (Eggleton et al. 1998; Naoz 2016):

[TABLE]

This quantity is a measure of the efficiency of energy dissipation due to tides. In general, the smallest $t_{Fs}$, the more efficient is the dissipation and therefore tidal evolution takes place more rapidly.

Please notice the strong dependence of this quantity on stellar radius. This dependence implies, for instance, that an increment in the star’s size by a factor of two, increases the rate of tidal dissipation by a factor of $256$. This sensitivity on the stellar radius is precisely what motivate us to study the effect of stellar evolution in the dynamics of these systems.

3.3 General relativity correction

Due to the excursions of the inner eccentricity to high values, relativistic effects are expected to have a major impact, especially on the evolution of the argument of periastron of the inner orbit (periastron precession).

In order to include this effect, we implemented the following Post-Newtonian correction to the Hamiltonian (Richardson & Kelly, 1988):

[TABLE]

with $\vec{v}_{1}$ being the relative velocity of the inner system and $\vec{P}_{1}$ the specific (relativistic) momentum, defined as:

[TABLE]

with $c$ the speed of light. In Appendix A we provide explicit expressions for the auxiliary functions $\gamma_{i}$.

If we approximate the momentum to first order in $v_{1}/c$, i.e. $\vec{P}_{1}\approx\vec{v}_{1}$, the relativistic doubly-averaged Hamiltonian correction will have error terms of the order $\mathcal{O}(v_{1}^{4}/c^{4})$ and can be written as:

[TABLE]

For a comprehensive explanation of the secularization process of this Hamiltonian, we refer the reader to the work of Migaszewski & Goździewski (2009). the final secularized relativistic Hamiltonian will be:

[TABLE]

and, from Equation (14), the corresponding general relativity contribution to equation of motion for the argument of periastron will be:

[TABLE]

that can be expressed explicitly in terms of the eccentricity and semi-major axes, as:

[TABLE]

In summary, the dynamical evolution of a hierarchical triple system obeys equations of motion with contributions coming from a “conservative” Hamiltonian, a non conservative contribution only associated to the interaction of $m_{1}$ and $m_{2}$, and last but not least, a relativistic correction on the argument of the pericenter of the inner orbit $g_{1}$.

Among all these effects, the non-conservative one strongly depends on the stellar radius. This dependency, and as we will show later, can not be neglected when studying the evolution of the system during those phases when the radius of the star changes in timescales similar to that of the orbital dynamics. We will focus here on the early phases (pre-main-sequence) when the radius decreases by a factor of 2-3 in low-mass stars.

4 Pre-main-sequence stellar evolution

Stars are formed in overdensed regions of a primogenitor molecular cloud. The progressive release of gravitational potential energy leads to the formation of a protostar which is still embedded inside the cloud making it visible only as an infrared source. When the cloud is dispersed and the protostar has reached its final mass, a pre-main-sequence star was born (Carroll & Ostlie, 2007).

For masses lower than $3\,M_{\odot}$, the radius and observed properties of the star (effective temperature, luminosity, surface gravity, etc.) evolves following what is called a Hayashi track (Hayashi, 1966). The particular track followed by a star will depends on its mass and composition.

Along its Hayashi track the luminosity of the pre-main-sequence stars drops while maintaining a roughly constant effective temperature. These changes are a by product of the sustained decrease in radius of the star (see Figure 2).

While in its Hayashi track, the star is characterized by having a fully convective interior with a thin radiative photosphere. However, due to the continuous increase in central temperature, the opacity of the core decreases leading to the formation of a radiative nucleus. At this point the central region rapidly heats reaching temperatures high enough to trigger nuclear fusion of hydrogen-1 into helium-4 (Harpaz, 1994). Hydrogen fusion provides enough thermal energy to guarantee hydrostatic equilibrium and the star reaches his final radius. At this point the star is said to be a zero age main sequence star (ZAMS).

The duration of the pre-main-sequence phase (from the arrival to the Hayashi track to ZAMS) is of the order of the thermal (or Kelvin-Helmholtz) timescale (Toonen et al., 2016):

[TABLE]

where here $L_{S}$ is the stellar instantaneous luminosity. Assuming the ZAMS values of $R_{s}$ and $L_{s}$ of a $0.6\,M_{\odot}$ mass star, $\tau_{\mathrm{th}}\approx{\cal O}(100$ Myr) which is consistent with the results of precise simulations. Since $R_{s}\propto L_{s}^{1/2}$ (for almost constant effective temperature) and $L_{s}\propto m_{s}^{3.5}$ (Carroll & Ostlie 2007; Demircan & Kahraman 1991), the timescale of pre-main-sequence collapse is $\tau_{\rm th}\propto m_{s}^{-3.25}$. This implies that more massive stars will have a considerably shorter pre-main-sequence, but also that low mass stars will spend several hundred of Myrs in this evolutionary phase.

In Figure 2 we see the evolution of stellar radius for three solar-metallicity stars with masses $1,\,0.8$ and $0.6$ solar masses. The results were obtained from the PARSEC v1.2555https://people.sissa.it/~sbressan/ evolutionary tracks (Bressan et al., 2012).

Vertical green line represents the onset of the zero age main sequence (ZAMS) phase in the $0.6\,M_{\odot}$ case. At the left of the vertical line, we can see how the radius evolves in the pre-main-sequence phase, changing by a factor of about $3$ during a time interval similar to that predicted by equation (25).

All what has been said here correspond to the evolution of an isolated star. In a binary system, the pre-main sequence phase could be altered by the presence of a companion, specially in the case of compact S-type systems (Toonen et al., 2016). We further discuss the potential implications of these effects in Section 6.

5 Numerical experiments

In order to test the effect that the change in stellar radius has in the dynamical evolution of S-type binaries, we have devised 4 numerical experiments:

• •

Experiment 1: Secular dynamical evolution of a real binary planet, namely HD 80606 b. This system has been already studied in literature and we can use it for verification purposes.

• •

Experiment 2: A generic S-type system assuming constant stellar radius. This experiment will allow us to illustrate the typical behavior of some dynamical quantities and define the “proxies” we will use to evaluate the effect of stellar radius evolution.

• •

Experiment 3: A systematic exploration of the parameter space (initial orbital parameters). We are interested on to identify the initial values $a_{1}$, $e_{1}$ and $i_{\mathrm{tot}}$ for which the role of a variable stellar radius is more important.

• •

Experiment 4: A new case of study. With the initial conditions identified in Experiment 3 we fully simulate a hypothetical system.

To have a configurable set of numerical scripts to perform these experiments, but also to allow other researchers to reproduce our results, we developed and published an open C package, SecDev3B666https://github.com/bayronportilla/SecDev3B. Provided a set of initial conditions, SecDev3B solves the complete set of differential equations that rules the secular evolution of a hierarchical system of three bodies. A detailed description of this package is included in the README of the GitHub repository.

5.1 Experiment 1: revisiting the case of HD 80606b

Here, we revisited the long-term dynamical evolution of one of the most studied binary planets, HD 80606b. For this purpose, we solve the secular equation of motion coming from the conservative (expanded to octupolar degree) and the relativistic terms of the Hamiltonian as well as the non-conservative equations of motion. For this particular experiment and in order to compare our results to previous ones, we assume constant stellar radius.

HD 80606b is an eccentric ($e=0.93$) hot Jupiter orbiting the G5-type star HD 80606 which is located at $190$ light-years away in the Ursa Major constellation. The host star is member of a binary system with the HD 80607. The average separation between the components of the binary system is $\sim 1200$ AU (Moutou et al., 2009).

Radial velocity measurements and the analysis of the Rossiter-McLaughlin anomaly has suggested that the planet is in a prograde, highly inclined orbit around its host star.

Wu & Murray (2003) has suggested that the observed configuration is the result of the perturbative effect of the outer star. They performed simulations of the system using a secular theory expanded to quadrupolar order including the same non-conservative effects described in Section 3.2 and a general relativity correction to the inner argument of pericenter.

With the set of initial conditions in Table 1 of Wu & Murray 2003, a solution was obtained. Their results suggest that the system underwent Lidov-Kozai oscillations during the first $2.7$ Gyr. Using the observed values of eccentricity and semi-major axis, they also conclude that the planet is not experiencing this effect anymore. Independently other authors (Fabrycky & Tremaine, 2007; Correia et al., 2011) studied the system with analogues results.

We use our package SecDev3B to reproduce these results. In Table 1 we present the initial conditions we use in our simulations. We prepare two different simulations in this case: 1) using the expansion of the conservative Hamiltonian up to quadrupolar order; and 2) using the expansion up to the octupolar term. The latter numerical experiment was never attempted before.

In order to see why the octupolar terms may become relevant in this particular case, we should notice that the octupolar term in equation (12) is proportional to the factor $\beta_{4}\sim(m_{1}-m_{2})$ (see Appendix A). This term can only be neglected when both members of the inner binary have similar masses. In this case, however, since $m_{1}/m_{2}\sim 150$, the octupolar term may become relevant. Our numerical results confirm this expectation.

In Figure 3 we present the result of our experiment with HD 80606b. The results to quadrupolar order are identical to those reported in previous works (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Correia et al., 2011). This confirm that our package is well-suited for this problem.

As a by product of this test, we find that the inclusion of the octupolar term leads to a significant modification of the dynamical evolution timescales. In particular, the time required to suppress the Lidov-Kozai oscillations, or in other terms, the time required for the oscillations in eccentricity to reach an amplitude less than $1\%$ of its initial value (see second row in Fig. 3) at quadrupolar-order is about $2.7$ Gyr whereas when the octupolar term is included this time becomes about $5.3$ Gyr.

This significant increase of the dynamical timescale is triggered by the emergence of longer “plateaus” in the inner semi-major axis and eccentricity oscillations (steps in the yellow shaded lines and in the semi-major axis curve in the first row of Fig. 3). When the system is subject to strong Lidov-Kozai oscillations, the inner semi-major axis is temporarily “captured” into quasi-stable states (where the value of the semi-major axis becomes constant for tens of Myrs). This “stable” periods are observed both when the quadrupolar and the octupolar terms are included. However they are much longer in the latter case, delaying the orbit decay. It is also worth noticing that the reached final values for all dynamical quantities are the same under both approximations.

Octupolar terms in the Hamiltonian implies the inclusion of extra harmonics into a finite series development. In the octupolar integration, we can notice the appearance of additional oscillation frequencies for each quantity plotted. This additional frequencies shift the solution and again produces larger timescales.

The actual impact this “octupolar delay” may have in our understanding of planets in binaries is diverse. Thus for instance, it could affect our expectation of the time spent for a planet in the habitable zone or the time required to reach or leave the unstable region (Kane & Hinkel 2013; Mason et al. 2015).

Should the octupolar order be included in every simulation? Not necessarily. In the case of S-type binary planets, since the inner system will be always composed by two bodies of masses considerably different, the octupolar term should always be included (unless the outer eccentricity is zero). In the case of P-type systems, where the components of the inner system have similar masses, this condition is relaxed reducing the contribution of the octupolar term to the total energy (see Equation 12).

5.2 Experiment 2: a generic S-type system

Once we verify that our code runs properly, we want to identify the characteristic features of the secular evolution of our hierarchical systems. For this purpose we simulate a hypothetical S-type binary starting with the initial conditions given in Table 1. In Figure 4 we present the result of our simulation in this case.

The first noticeable feature is the “ladder-shaped” evolution of the semi-major axis. This behavior was already noticed in Experiment 1, however here it is easier to identify the possible reasons of this behavior. The almost sudden drop in the semi-major axis are correlated with excursions to high-eccentricity values in the inner system. Thus, for instance, the first maximum of $e_{1}$ is $0.985$ reached after $1.2$ Myr, which is the same time when the first drop in $a_{1}$ happens. We must keep in mind that in this secular formalism the only phenomena able to modify the value of the inner semi-major axis are the non-conservative effects. The correlation between the drops in $a_{1}$ and the excursion of $e_{1}$ becomes evident if we have in mind that the magnitude of the tidal acceleration goes as $r^{-10}$ (see equation 17), and therefore, it is much bigger in each periapse passage. When eccentricity gets its highest value, the periapse distance gets very small (for a given value of the semi-major axis). For instance, at $t=1.2$ Myr, the periastron distance is just $\sim 0.04$ AU. This distance is small enough to trigger a strong tidal interaction which is the responsible of the sudden drop in the value of $a_{1}$.

Our results also evidence the strong anticorrelation between the inner eccentricity and the mutual inclination. This anticorrelation arise from the interchange of angular momentum between both the inner and outer orbits. This is precisely what we call the Lidov-Kozai mechanism (Lidov 1962; Kozai 1962).

In our experiments, we have identified that the time required for the semi-major axis to reduce to the half of its initial value ($t_{a,1/2}$):

[TABLE]

is a good proxy of the inner orbit evolution.

We could also define similar time scales for other quantities such as eccentricity or inclination. However, since the evolution in eccentricity is mainly controlled by the conservative terms of the Hamiltonian, the introduction of the evolution in the stellar radius will not modify significantly the value of those other proxies of dynamical evolution.

5.3 Experiment 3: the effect of stellar evolution

In order to test how the evolution of the stellar radius impact the timescales of orbital evolution in S-type binaries, we should explore the space of all possible initial conditions. If we consider only the orbital evolution the space to explore has 11 dimensions (eight possible initial conditions and three masses). The size of the parameter space is significantly increased if we include the rotational parameters (moment of inertial, radius, rotational rates, etc.)

Since this paper is intended only to test the hypothesis that the evolution of stellar radius may have a significant effect on the dynamical evolution of hierarchical systems, we restrict our selves to “subspaces” of the whole parameter space. In particular we will restrict to two sections $a_{1}\,\mathrm{vs.}\,e_{1}$ and $a_{1}\,\mathrm{vs.}\,i_{\mathrm{tot}}$. For this purpose we fixed all the remaining relevant physical and orbital parameters.

For the $a_{1}$ vs. $e_{1}$ subspace we set initial mutual inclination in $80^{\circ}$ (which is as large as that observed in HD 180606b). For the $a_{1}$ vs. $i_{\mathrm{tot}}$ subspace the initial value of the inner eccentricity was set in $0.6$.

5.3.1 The $a_{1}$ vs. $e_{1}$ space

The exploration of this subspace begins with the construction of a $10\times 10$ grid of $(a_{1},e_{1})$ values. The semi-major axes cover a range between $0.85$ AU to $4.0$ AU with a constant separation of $0.35$ AU. On the other hand, inner eccentricity goes from $0.18$ to $0.9$ in steps of $0.08$.

Each point in the grid represents an individual simulation with initial conditions given by the coordinates of the point. The remaining initial parameters are summarized in Table 1.

A first set of simulations was made maintaining a constant value of the stellar radius. The decay time of the semi-major axes ($t_{a,1/2}$) was recorded for each simulation. Next, a second set of $100$ simulations was ran with the same initial conditions, but this time letting the radius to change according to the predicted evolution of a $0.6\,M_{\odot}$ star (see Fig. 2). Once again, the decay time for $a_{1}$ was also recorded.

If after $5000$ Myr, the inner semi-major does not match the decay condition, namely $a_{1}(t=5000)\leq a_{1}(t=0)/2$, the simulation was stopped and the value assigned to $t_{a,1/2}$ was $5000$ Myr. Once all the simulations were performed, each point in the grid gets a value of $t_{a,1/2}$. In other words, we obtained a set of discrete sample points of the function $t_{a,1/2}=t_{a,1/2}(a_{1},e_{1})$. Finally, we interpolate the discrete values and obtain contour plots of $t_{a,1/2}(a_{1},e_{1})$. The result of our simulations are shown in the contours of Figure 5.

There are noticeable differences in decay times for low initial semi-major axes. This is not surprising since the effect of variable radius is relevant for the non-conservative interaction, which has a strong dependence on the average inner system separation (Equation 19).

To analyze more closely this region, we made a “zoom” in the area defined by $0.8\text{ AU }\leq a_{1}\leq 1.5\text{ AU }$ and $0.18\leq e_{1}\leq 0.6$. The results are shown in the panels of the lower row of Figure 5. Differences in decay times are again quite evident. Thus, for instance, for the selected initial conditions $(a_{1},e_{1})=(1,0.35)$ (white dots), the decay time with constant stellar radius is of the order of $700$ Myr. Including the effect of a variable stellar radius, the decay time drops to $150$ Myr.

To explain this behavior, we should notice that at the beginning of the orbital evolution in the case of variable stellar radius, the initial value of this quantity was $\sim 2R_{\odot}$. This is almost 3 times bigger than the ZAMS radius which was precisely the value used in the constant stellar radius case. Naturally, the higher the stellar radius, higher is the energy dissipation due to tides (see Equation 19). Details of the evolution of a hypothetical system having this initial conditions are presented and discussed in Section 5.4.

5.3.2 The $a_{1}$ vs. $i_{\mathrm{tot}}$ space

We repeat the previous procedure but now in the $a_{1}$ vs. $i_{\mathrm{tot}}$ subspace. The grid has the same range and periodicity in $a_{1}$ than before and the inclination was varied from $50^{\circ}$ to $85^{\circ}$ with a constant step of $4^{\circ}$. A zoomed area was selected, spanning $0.6\text{ AU }\leq a_{1}\leq 1.4\text{ AU}$ and $62.5^{\circ}\leq i_{\mathrm{tot}}\leq 85^{\circ}$. The initial value of the inner eccentricity was fixed in $0.6$ while the rest of binary parameters are in Table 1. The resulting contour plots of $t_{a,1/2}(a_{1},i_{\mathrm{tot}})$ are shown in Figure 6.

Our proxy is now less sensitive to the selected set of parameters. At low inclinations, the inner semi-major axis stays stable for a long-time (it does not decay). Lidov-Kozai oscillations have very small-amplitude and are unable to bring the planet close enough to the star for the tidal interaction to have a real impact on the timescale of orbital evolution.

More importantly, at this specific eccentricity and range of inclinations, the orbital decay times are much less sensitive to the evolution in stellar radius.

5.4 Experiment 4: a hypothetical evolving S-type system

For our last experiment, we simulate the orbital evolution of a system having the properties of a selected spot (white point in Figure 5) in the $a_{1}$ vs. $e_{1}$ subspace. Initial conditions for this system are enumerated in Table 1. Results of this experiment are summarized in Figure 7.

Our planet belongs to a “twin” binary of 0.6 $M_{\odot}$ stars which orbit each other in an eccentric orbit ($e_{2}=0.3$) with a semi-major axis of 50 AU. The planet is initially on a mildly eccentric orbit ($e_{1}=0.35$) at $1$ AU from its host star.

According to stellar evolution models the pre-main-sequence of the host star takes around 200 Myr (see upper panel in Fig. 7) during which the stellar radius varied from about $2\,R_{\odot}$ to $0.6\,R_{\odot}$ (Bressan et al., 2012). Interestingly, for the particular orbital properties, the timescales of stellar evolution are of the same order than the timescales of orbital evolution.

It is significant the effect that an evolving stellar radius has in the timescale of dynamical evolution. While in the case of a constant radius star (left panel in Fig. 7) it takes almost 200 Myr for the orbit to decay at about half of its initial value, when the stellar evolution is turned on the time-scale is reduced to less than 100 Myr. Even more dramatic is reduction in 362 Myr in the case of an evolving star, of the total time required for the end of the high-amplitude Lidov-Kozai oscillations.

In both cases however, the orbital properties at the end of the simulation was almost identical, namely, $a_{1,\mathrm{end}}=0.15$ AU, while the final eccentricity were similar but not identical ($e_{1,\mathrm{end}}=0.83$ in the constant-radius case and $e_{1,\mathrm{end}}=0.86$ in the evolving one).

6 Discussion

Stellar and planetary evolution in binaries does not happen in isolation. In the presence of a binary companion, different physical processes may affect the formation and evolution not only of the planet but the stars itself.

In the case of P-type circumbinary systems, the distance between both stars are in general relatively small. When the radius of one of the stellar member is larger than its Roche lobe, mass could be transferred to the companion. This phenomena is known as Roche lobe overflow (RLOF) and may occur during the protostellar evolution and even the pre-main-sequence phase (Ivanova, 2015). When matter is transferred from one star to another, the hydrostatic and thermal equilibrium of the donor star are disturbed and the radius and timescales shown in Figure 2 are modified.

To evaluate under which circumstances an early RLOF may affect pre-main-sequence evolution we need to estimate the Roche radius $R_{L1}$. For a binary system with a mass ratio $q=m_{s1}/m_{s2}$ it is given by (Eggleton, 1983):

[TABLE]

For the system we simulate in experiment 4 (see Section 5.4), $q=1$ and $a_{\mathrm{bin}}=50$ AU, the Roche radius is about $19$ AU. Therefore, even when the highest expected possible value of the stellar radius is assumed, it is less than $0.04\%$ the Roche radius. RLOF effects on stellar evolution in this systems can be safely neglected and we can assume that both stars evolve following the models of isolated stars.

Stars transfer angular momentum not only via tidal interactions with stellar and planetary companions. Mass-loss and magnetic breaking are other mechanisms in action, and they could be very strong, especially during the early phases of stellar evolution (Huang, 1963; Yakut et al., 2008; Zuluaga et al., 2016). These effects may also impact the total angular momentum budget in the system affecting orbital evolution.

Mass-losses induced by strong stellar winds at the beginning of the life of a low mass-star, may have an impact on the evolution of the binary orbit itself. Thus, for instance, For a $0.6\,M_{\odot}$ T-tauri star the mass loss rates are as large as $3.5\times 10^{-8}\,M_{\odot}/yr$ (Ezer & Cameron, 1971). If the ejection process is adiabatic and none of the wind-matter is accreted by the companion star, the average separation between both bodies changes according to (Toonen et al., 2016):

[TABLE]

when we use the notation adopted in Section A.1, $a_{\rm bin}$ will be $a_{1}$ or $a_{2}$ depending on the nature of the system (S or P type). In this adiabatic case, the eccentricity of the binary orbit ($e_{1}$ or $e_{2}$ depending on the type) is not modified by the mass-loss.

On the other hand, if the expelled material is additionally accreted by the companion star, there will be a transfer of angular momentum between them and binary eccentricity will not be constant any more (Huang 1963; Toonen et al. 2016).

None of these effects were taken into account in our simplified models. This simple estimations however, contribute to reveal the complex role that stellar evolution may play at determining the orbital fate of planet in binaries.

Planetary formation is also an important factor in the dynamical evolution of planets in binaries. It will not only determine the initial conditions at which our simulations start, but several key early processes may affect orbital evolution. Circumstellar (or circumbinary) disks may last between $0.1$ Myr to $10$ Myr (Haisch et al., 2001; Pelupessy & Portegies Zwart, 2013; Haghighipour, 2010). Although these times are relatively small with respect to the few hundreds of Myr involved in the orbital evolution of planet in binaries, the presence of a (perturbed) disk during the first phases of orbital evolution should be taken into account.

Haghighipour (2010) performed numerical experiments to study the effect on orbital evolution of an S-type planet when it was embedded in a circumstellar disk. They performed simulations of giant planets with and without accretion. Their results showed that binary planets in this type of systems may migrate inwards and get their eccentricity excited both by the disk and the outer companion (see Fig. 6.5 of Haghighipour 2010). If the planet accretes gas from the disk the timescales of semi-major axis and eccentricity evolution, as well as their asymptotic values were significantly modified with respect to the non-accretion case.

7 Summary and conclusions

We revisit in this paper the well-known problem of hierarchical three body systems applied to the case of planets in binaries. For this purpose we deduce and provide the full general secularized Hamiltonian up octupolar order in Delaunay variables. In the equations of motion, we also included non-conservative (tidal interaction and quadrupolar deformation) and general relativity terms of the interaction between the inner bodies. To solve the equations under general initial conditions we developed an open source package that is properly documented and can be used to reproduce our results. During a given integration, SecDev3B is well suited to include or not the contributions from non-conservative and relativistics effects as desired by the user. It is also possible to explore how the change of the inner bodies’ bulk properties affects the dynamics as those quantities can be treated as variables whose value will be interpolated from a provided data set. We apply this model to study the dynamical evolution of S-type binary planets.

Using this formalism we investigated the secular dynamics of the well-known binary planet HD 80606b. We find that using a conservative Hamiltionian expanded to quadrupolar order our software reproduces properly the results found in literature. However, and as a novel contribution, when we used the expansion to octupolar order, the evolutionary timescales were significantly modified. This result highlight the importance to use high-order developments when dealing with hierarchical S-type binaries which does not match the requirement of having inner bodies of equal mass.

We then apply the model to study the role that the evolution of stellar radius during the pre-main-sequence phase could have in the secular dynamics of the system. For this purpose we explore a small part of a $>$12-dimensional parameter space, in particular, subspaces of semi-major axes, eccentricities and mutual inclinations. We found that in S-type systems around low-mass stars, with relative high inclinations ($i_{\mathrm{tot}}\geq 60\deg$), moderate eccentricities ($0.2\leq e\leq 0.4$) and for planets located around 1 AU (typical distances of the Habitable zones of the host stars), the changes experienced by the star during the first few hundreds of Myr alter significantly the dynamical evolution timescales. In a particular case of study we found that a planet decay from an initial low eccentric orbit of 1 AU to a final eccentric orbit of 0.1 AU around its host star in a time almost 400 Myr shorter than that calculated assuming the star have a constant radius.

This first approach is still incomplete. Other regions of the parameter space should be explored and non-trivial effects such as early stellar mass-loss, binary angular momentum interchange, planetary radius evolution, the role of a circumstellar disk during the first Myrs, among other should be considered in future developments of this exciting topic.

The more important conclusion derived from this work is that for assessing the fate and history of the orbital properties of planets in binaries, not only dynamical issues should be at play. It is worth realizing that the evolution of the orbits of planets and stars in binaries, especially during the first hundred of Myrs, actually happens to evolving non-static objects.

Acknowledgements

We are grateful with all the colleagues that through the discussion in seminars and presentations of this work in conferences directly and inderectly contribute to create its final form. We especially thanks to Josué Cardoso dos Santos for his insightful comments. B.P. was funded by the program Jóvenes Investigadores e Innovadores 2016 of COLCIENCIAS through grant JIC 018-2017. He wants to thanks to the people in the Vicerrectoría de Investigación for his support during the program.

Appendix A Detailed mathematical expressions

In this appendix we will give explicit expressions for the auxiliary constants, functions and equation of motion mentioned through the paper.

A.1 Conservative components

Quantities given in expressions 29 appears explicitly in the expansion of the three body (non-secularized) Hamiltonian. In order to avoid confussion, please notice we have adopted the caligraphyc letter $\mathcal{G}$ to represent the Cavendish gravitational constant.)

[TABLE]

The literal expression of the doubly averaged Hamiltonian (before the node elimination) is cumbersome. To simplify it, we provide in the paper the complete expansion in terms of the auxiliary quantities $K_{i}\,(i=1,2,3,4)$ that are a function of the orbital elements of the inner and outer orbits.

[TABLE]

The direct application of the Hamilton canonical equations to the leads to the following

The equation of motion derived from the doubly averaged Hamiltonian expanded to octupolar order are provided below. Subindex “orb" emphasizes the fact that the expression is valid only for the conservative (or orbital) regime. To simplify the expression we have adopted the notation of Naoz et al. (2013).

[TABLE]

As was explained in Section A.1, the proper way to derive the equations of motion for the orbital inclinations is through the conservation law of the total angular momentum $\vec{G}_{\mathrm{tot}}$. For the sake of completeness, we must to strees the fact that this quantity is the vector summation of the orbital plus the spin angular momentum. However, the former component is by far bigger than the latter and we can work with the approximation: $\vec{G}_{\mathrm{tot}}\approx\vec{G}_{1}+\vec{G}_{2}$. The The magnitudes of the inner and outer orbital angular momentum are given by

[TABLE]

and the $z-$components are just:

[TABLE]

The magnitude of the total angular momentum can be written in terms of those $z-$ components as $G_{\mathrm{tot}}=H_{1}+H_{2}$. The rate of change of those angular momenta are given by:

[TABLE]

In terms of those quantities, the equations of motion for the inclinations are:

[TABLE]

In previous equations, the following auxiliary quantities have been used:

[TABLE]

A.2 Relativistic correction

Many of the definitions given in this section are based in the notation of Richardson & Kelly (1988) and Migaszewski & Goździewski (2009). If we define $\sigma\equiv m_{1}m_{2}/(m_{1}+m_{2})^{2}$ and $\mu\equiv G(m_{1}+m_{2})$. The auxiliary factors $\gamma_{i}$ appearing in equation (22) are functions of the masses as follows:

[TABLE]

A.3 Non-conservative motion

We define the reduced mass of the inner system as $m_{r}=m_{1}m_{2}/(m_{1}+m_{2})$ and the mean motion as $n=\sqrt{\mu/a_{1}^{3}}$. The equations of motion of the non-conservative part are written in terms of the following auxiliary functions of the inner eccentricity.

[TABLE]

Please notice that the definition of $f_{i}(e_{1})$ differs from that adopted in Correia et al. (2011).

The following auxiliary functions of the rheological parameters and rotational velocity of each body should be defined to write down the non-conservative components of the equation of motion.

[TABLE]

In terms of the previous definitions, the secularized non-conservative (nc) effects on the inner orbit and on the spin of the inner bodies are given by:

[TABLE]

Where the moment of inertia of the inner components are:

[TABLE]

with $\kappa_{i}$ ($i=s,p$) are the so-called gyration radius of the involved bodies.