The dynamical evolution of close-in binary systems formed by a super-Earth and its host star. Case of the Kepler-21 system

TL;DR

This paper develops a formalism to study the long-term orbital and rotational evolution of close-in binary systems, including tidal and triaxial effects, and applies it to the Kepler-21 system to understand its dynamical behavior.

Contribution

It introduces a new set of singularity-free equations of motion for binary systems with tidal and triaxial effects, applicable to systems like Kepler-21.

Findings

Tidal and triaxiality effects can significantly influence orbital evolution.

Time to tidal synchronization depends strongly on rheological parameters.

The formalism accurately models the evolution near spin-orbit resonances.

Abstract

The aim of this work is to develop a formalism for the study of the secular evolution of a binary system which includes interaction due to the tides that each body imparts on the other. We also consider the influence of the $J_2$-related secular terms on the orbital evolution and the torque, caused by the triaxiality, on the rotational evolution, both of which are associated only to one of the bodies. We apply these set of equations to the study of the orbital and rotational evolution of a binary system composed of a rocky planet and its host star in order to characterize the dynamical evolution at work, particularly near spin-orbit resonances. We used the equations of motion that give the time evolution of the orbital elements and the spin rates of each body to study the time evolution of the Kepler-21 system as an example of how the formalism that we have developed can be applied. We obtained a set of equations of motion without singularities for vanishing eccentricities and inclinations. This set gives, on one hand, the time evolution of the orbital elements due to the tidal potentials generated by both members of the system as well as the triaxiality of one of them. On the other hand, it gives the time evolution of the stellar spin rate due to the corresponding tidal torque and of the planet's rotation angle due to both the tidal and triaxiality-induced torques. We found that for the parameters and the initial conditions explored here, the tidally and triaxiality-induced modifications of the tidal modes can be more significative than expected and that the time of tidal synchronization strongly depends on the values of the rheological parameters.

Figures

Tables

Table 1. Table 1 : Inclination functions for l = 2 𝑙 2 l=2 and their first derivatives.

$m$	$p$	$F_{2mp}\left(i\right)$	$\frac{d{F}_{2mp}\left(i\right)}{di}$	$\frac{1}{\sin i}\left(F_{2mp}\left(i\right)\frac{d{F}_{2mp}\left(i\right)}{di}\right)$
0	0	$\frac{3}{8}{\sin }^{2}i$	$\frac{3}{4}\sin i\cos i$	$\frac{9}{32}{\sin }^{2}i\cos i$
0	1	$-\frac{3}{4}{\sin }^{2}i+\frac{1}{2}$	$-\frac{3}{2}\sin i\cos i$	$-\frac{3}{2}\cos i\left(\frac{1}{2}-\frac{3}{4}{\sin }^{2}i\right)$
0	2	$\frac{3}{8}{\sin }^{2}i$	$\frac{3}{4}\sin i\cos i$	$\frac{9}{32}{\sin }^{2}i\cos i$
1	0	$-\frac{3}{4}\left(1+\cos i\right)\sin i$	$-\frac{3}{4}\left[\cos i+\cos \left(2i\right)\right]$	$\frac{9}{16}\left(1+\cos i\right)\left[\cos i+\cos \left(2i\right)\right]$
1	1	$\frac{3}{2}\sin i\cos i$	$\frac{3}{2}\cos \left(2i\right)$	$\frac{9}{4}\cos i\cos \left(2i\right)$
1	2	$\frac{3}{4}\left(1-\cos i\right)\sin i$	$\frac{3}{4}\left[\cos i-\cos \left(2i\right)\right]$	$\frac{9}{16}\left(1-\cos i\right)\left[\cos i-\cos \left(2i\right)\right]$
2	0	$\frac{3}{4}{\left(1+\cos i\right)}^2$	$-\frac{3}{2}\left(1+\cos i\right)\sin i$	$-\frac{9}{8}{\left(1+\cos i\right)}^3$
2	1	$\frac{3}{2}{\sin }^{2}i$	$3\sin i\cos i$	$\frac{9}{2}{\sin }^{2}i\cos i$
2	2	$\frac{3}{4}{\left(1-\cos i\right)}^2$	$\frac{3}{2}\left(1-\cos i\right)\sin i$	$\frac{9}{8}{\left(1-\cos i\right)}^3$

Table 2. Table 2 : RMSD between the left and right-hand sides of Eqs. ( 47 ), ( 48 ), ( 53 ), and ( 60 ). Columns four through seven correspond to the RMSD of the real ( ℜ \Re ) and imaginary ( ℑ \Im ) parts of the left and right hands of Eqs. ( 53 ) and ( 60 ), respectively. This results are obtained with e = 0.1 𝑒 0.1 e=0.1 .

$q_{\max }$	RMSD
	(47)	(48)	$\Re$ (53)	$\Im$ (53)	$\Re$ (60)	$\Im$ (60)
2	$4.756939\times {10}^{-3}$	$1.443418\times {10}^{-1}$	$1.237304\times {10}^{-2}$	$1.237193\times {10}^{-2}$	$3.704831\times {10}^{-1}$	$3.704534\times {10}^{-1}$
4	$9.884915\times {10}^{-5}$	$4.962669\times {10}^{-3}$	$4.149271\times {10}^{-4}$	$4.148743\times {10}^{-4}$	$2.064961\times {10}^{-2}$	$2.064789\times {10}^{-2}$
6	$1.997451\times {10}^{-6}$	$1.399011\times {10}^{-4}$	$1.175513\times {10}^{-5}$	$1.175534\times {10}^{-5}$	$8.183466\times {10}^{-4}$	$8.182644\times {10}^{-4}$
8	$4.064722\times {10}^{-8}$	$3.563261\times {10}^{-6}$	$3.117166\times {10}^{-7}$	$3.051717\times {10}^{-7}$	$2.708817\times {10}^{-5}$	$2.717803\times {10}^{-5}$
10	$1.000686\times {10}^{-8}$	$1.142291\times {10}^{-7}$	$7.273942\times {10}^{-8}$	$1.078445\times {10}^{-7}$	$9.280304\times {10}^{-7}$	$9.066432\times {10}^{-7}$
12	$9.977323\times {10}^{-9}$	$7.538773\times {10}^{-8}$	$7.213027\times {10}^{-8}$	$1.076484\times {10}^{-7}$	$4.663512\times {10}^{-7}$	$4.368208\times {10}^{-7}$

Table 3. Table 3 : Physical and orbital parameters of Kepler-21 system studied in this work (López-Morales et al., 2016 ) .

Parameter	Symbol	Value
Stellar parameters
Mass ($M_{\odot }$)	$m_1$	$1.408$
Radius ($R_{\odot }$)	$R_1$	$1.902$
Rotational period (Days)	$P_{\mathrm{rot},1}$	$12.6$
Planetary parameters
Mass ($M_{\oplus }$)	$m_2$	$5.10$
Radius ($R_{\oplus }$)	$R_2$	$1.639$
Orbital parameters
Orbital period (Days)	$P_{\mathrm{orb}}$	$2.78578$
Major semiaxis (au)	$a$	$0.0427172$
Eccentricity	$e$	$0.02$

Table 4. Table 4 : Values of the rheological and dynamical parameters adopted in the numerical integrations performed in this work, together with the corresponding ID number assigned to ease the identification of the result of each simulation.

	ID	$e$	${\tau }_M$ [yr]	$\alpha$	$\left(B-A\right)/C_2$	$J_2$
	1	$0.02$	$10$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-3}$
	2	$0.02$	$50$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-3}$
Set 1	3	$0.02$	$100$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-3}$
	4	$0.02$	$10$	$0.15$	$7\times {10}^{-4}$	$1\times {10}^{-3}$
	5	$0.02$	$10$	$0.3$	$7\times {10}^{-4}$	$1\times {10}^{-3}$
	6	$0.02$	$10$	$0.2$	$1\times {10}^{-3}$	$1\times {10}^{-3}$
	7	$0.02$	$10$	$0.2$	$1\times {10}^{-5}$	$1\times {10}^{-3}$
	8	$0.02$	$10$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-4}$
Set 2	9	$0.02$	$10$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-5}$
	10	$0.02$	$10$	$0.2$	$0$	$1\times {10}^{-3}$
	11	$0.02$	$10$	$0.2$	$7\times {10}^{-4}$	$0$
	12	$0.02$	$10$	$0.2$	$0$	$0$
Set 3	13	$0.1$	$10$	$0.2$	$7\times {10}^{-4}$	$1\times {10}^{-3}$

Table 5. Table 5 : Orbital decay and eccentricity damping timescales of the Kepler-21 system obtained through Eqs. ( 76 ) and ( 77 ) by evaluating the corresponding equation of motion at θ ˙ 2 / n = 1 subscript ˙ 𝜃 2 𝑛 1 \dot{\theta}_{2}/n=1 for different values of the rheological parameters and the eccentricity.

$e$	${\tau }_M$ [yr]	$\alpha$	${\tau }_a$ [$\times {10}^6$ yr]	${\tau }_e$ [$\times {10}^6$ yr]
$0.01$	$10$	$0.2$	$751.32$	$0.101$
$0.02$	$10$	$0.2$	$629.23$	$0.345$
$0.02$	$50$	$0.2$	$662.78$	$0.361$
$0.02$	$100$	$0.2$	$675.81$	$0.367$
$0.02$	$10$	$0.15$	$617.57$	$0.339$
$0.02$	$10$	$0.3$	$679.20$	$0.368$
$0.1$	$10$	$0.2$	$98.955$	$1.542$

Equations

$$G\left(m_1+m_2\right)=n^2{a}^3,$$

$$\frac{\mathrm{d}a}{\mathrm{d}t}=\frac{2}{na}\frac{\partial R}{\partial M},$$

$$\frac{\mathrm{d}e}{\mathrm{d}t}=\frac{1-e^2}{n{a}^{2}e}\frac{\partial R}{\partial M}-\frac{\sqrt{1-e^2}}{n{a}^{2}e}\left(\frac{\partial R}{\partial {\omega }_1}+\frac{\partial R}{\partial {\omega }_2}\right),$$

$$\frac{\mathrm{d}M}{\mathrm{d}t}=n-\frac{2}{na}\frac{\partial R}{\partial a}-\frac{1-e^2}{n{a}^{2}e}\frac{\partial R}{\partial e},$$

$$\begin{gathered} \frac{\mathrm{d}{i}_1}{\mathrm{d}t}=\frac{\beta }{C_1{\dot{\theta }}_1\sin i_1}\left(\frac{\partial R}{\partial {\omega }_1}-\cos i_1\frac{\partial R}{\partial {\Omega }_1}\right)-\frac{1}{n{a}^2\sqrt{1-e^2}\sin i_1}\left(\frac{\partial R}{\partial {\Omega }_1}-\cos i_1\frac{\partial R}{\partial {\omega }_1}\right) \\ -\frac{\sin ({\omega }_1-{\omega }_2)}{n{a}^2\sqrt{1-e^2}}\frac{\partial R}{\partial i_2}-\frac{\cos ({\omega }_1-{\omega }_2)}{n{a}^2\sqrt{1-e^2}\sin i_2}\left(\frac{\partial R}{\partial {\Omega }_2}-\cos i_2\frac{\partial R}{\partial {\omega }_2}\right), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{i}_2}{\mathrm{d}t}=\frac{\beta }{C_2{\dot{\theta }}_2\sin i_2}\left(\frac{\partial R}{\partial {\omega }_2}-\cos i_2\frac{\partial R}{\partial {\Omega }_2}\right)-\frac{1}{n{a}^2\sqrt{1-e^2}\sin i_2}\left(\frac{\partial R}{\partial {\Omega }_2}-\cos i_2\frac{\partial R}{\partial {\omega }_2}\right) \\ -\frac{\sin ({\omega }_2-{\omega }_1)}{n{a}^2\sqrt{1-e^2}}\frac{\partial R}{\partial i_1}-\frac{\cos ({\omega }_2-{\omega }_1)}{n{a}^2\sqrt{1-e^2}\sin i_1}\left(\frac{\partial R}{\partial {\Omega }_1}-\cos i_1\frac{\partial R}{\partial {\omega }_1}\right), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{\omega }_1}{\mathrm{d}t}=-\frac{\beta }{C_1{\dot{\theta }}_1\sin i_1}\frac{\partial R}{\partial i_1}+\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\partial R}{\partial e}-\frac{\cos i_1}{n{a}^2\sqrt{1-e^2}\sin i_1}\frac{\partial R}{\partial i_1} \\ -\frac{\cos i_1\cos \left({\omega }_1-{\omega }_2\right)}{n{a}^2\sqrt{1-e^2}\sin i_1}\frac{\partial R}{\partial i_2}-\frac{\cos i_1\sin \left({\omega }_1-{\omega }_2\right)}{n{a}^2\sqrt{1-e^2}\sin i_1\sin i_2}\left(\cos i_2\frac{\partial R}{\partial {\omega }_2}-\frac{\partial R}{\partial {\Omega }_2}\right), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{\omega }_2}{\mathrm{d}t}=-\frac{\beta }{C_2{\dot{\theta }}_2\sin i_2}\frac{\partial R}{\partial i_2}+\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\partial R}{\partial e}-\frac{\cos i_2}{n{a}^2\sqrt{1-e^2}\sin i_2}\frac{\partial R}{\partial i_2} \\ -\frac{\cos i_2\cos \left({\omega }_2-{\omega }_1\right)}{n{a}^2\sqrt{1-e^2}\sin i_2}\frac{\partial R}{\partial i_1}-\frac{\cos i_2\sin \left({\omega }_2-{\omega }_1\right)}{n{a}^2\sqrt{1-e^2}\sin i_2\sin i_1}\left(\cos i_1\frac{\partial R}{\partial {\omega }_1}-\frac{\partial R}{\partial {\Omega }_1}\right), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{\Omega }_1}{\mathrm{d}t}=\frac{\beta \cos i_1}{C_1{\dot{\theta }}_1\sin i_1}\frac{\partial R}{\partial i_1}+\frac{\beta \cos {\epsilon }_1}{C_1{\dot{\theta }}_1\sin {\epsilon }_1}\left(\sin {\Omega }_1\cot i_1\frac{\partial R}{\partial {\Omega }_1}-\cos {\Omega }_1\frac{\partial R}{\partial i_1}-\frac{\sin {\Omega }_1}{\sin i_1}\frac{\partial R}{\partial {\omega }_1}\right)+\frac{1}{n{a}^2\sqrt{1-e^2}\sin i_1}\frac{\partial R}{\partial i_1} \\ +\frac{\cos \left({\omega }_2-{\omega }_1\right)}{n{a}^2\sqrt{1-e^2}\sin i_1}\frac{\partial R}{\partial i_2}+\frac{\sin \left({\omega }_1-{\omega }_2\right)}{n{a}^2\sqrt{1-e^2}\sin i_1\sin i_2}\left(\cos i_2\frac{\partial R}{\partial {\omega }_2}-\frac{\partial R}{\partial {\Omega }_2}\right), \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{\Omega }_2}{\mathrm{d}t}=\frac{\beta \cos i_2}{C_2{\dot{\theta }}_2\sin i_2}\frac{\partial R}{\partial i_2}+\frac{\beta \cos {\epsilon }_2}{C_2{\dot{\theta }}_2\sin {\epsilon }_2}\left(\sin {\Omega }_2\cot i_2\frac{\partial R}{\partial {\Omega }_2}-\cos {\Omega }_2\frac{\partial R}{\partial i_2}-\frac{\sin {\Omega }_2}{\sin i_2}\frac{\partial R}{\partial {\omega }_2}\right)+\frac{1}{n{a}^2\sqrt{1-e^2}\sin i_2}\frac{\partial R}{\partial i_2} \\ +\frac{\cos \left({\omega }_1-{\omega }_2\right)}{n{a}^2\sqrt{1-e^2}\sin i_2}\frac{\partial R}{\partial i_1}+\frac{\sin \left({\omega }_2-{\omega }_1\right)}{n{a}^2\sqrt{1-e^2}\sin i_2\sin i_1}\left(\cos i_1\frac{\partial R}{\partial {\omega }_1}-\frac{\partial R}{\partial {\Omega }_1}\right), \end{gathered}$$

$$\text{@vec}g={\left(\frac{\partial \text{@vec}f}{\partial t}\right)}_{c_k=const.},$$

$$\dot{\text{@vec}r}=\frac{\mathrm{d}\text{@vec}r}{\mathrm{d}t}={\left(\frac{\partial \text{@vec}f}{\partial t}\right)}_{c_k=const.}+\sum _{k=1}^6\frac{\partial \text{@vec}f}{\partial c_k}\frac{\mathrm{d}{c}_k}{\mathrm{d}t}=\text{@vec}g+\text{@vec}\Phi ,$$

$$\sum _{k=1}^6\frac{\partial \text{@vec}f}{\partial C_k}\frac{\mathrm{d}{C}_k}{\mathrm{d}t}=\text{@vec}\Phi .$$

$$\frac{\mathrm{d}{\epsilon }_k}{\mathrm{d}t}=-\frac{\beta }{C_k{\dot{\theta }}_k}\left(\cos {\Omega }_k\cot i_k\frac{\partial R}{\partial {\Omega }_k}+\sin {\Omega }_k\frac{\partial R}{\partial i_k}-\frac{\cos {\Omega }_k}{\sin i_k}\frac{\partial R}{\partial {\omega }_k}\right),$$

$$\frac{{\mathrm{d}}^2{\theta }_k}{\mathrm{d}{t}^2}=-\frac{\beta }{C_k}\frac{\partial R}{\partial {\Omega }_k}.$$

$$R=-\frac{1}{\beta }\left[m_1\left[U_2\left(\text{@vec}{r}_{2\to 1},\text{@vec}{r}_{2\to 1}^{\ast }\right)+V_2\left(\text{@vec}{r}_{2\to 1}\right)\right]+m_2{U}_1\left(-\text{@vec}{r}_{1\to 2},-\text{@vec}{r}_{1\to 2}^{\star }\right)\right],$$

$$U_2\left(\text{@vec}{r}_1,\text{@vec}{r}_1^{\ast }\right)=-\frac{G{m}_1}{a^{\ast }}\sum _{l=2}^{\infty }\sum _{m=0}^l\sum _{p=0}^l\sum _{q=-\infty }^{\infty }\sum _{h=0}^l\sum _{s=-\infty }^{\infty }{U}_{lmpqhs}^{(2)},$$

$$\begin{gathered} U_{lmpqhs}^{(2)}={\left(\frac{R_2}{a}\right)}^{l+1}{\left(\frac{R_2}{a^{\ast }}\right)}^l\frac{(l-m)!}{(l+m)!}(2-{\delta }_{0m})F_{lmp}(i_2^{\ast })G_{lpq}(e^{\ast })F_{lmh}(i_2)G_{lhs}(e)[\cos \left[\left(v_{lmpq}^{(2),\ast }-m{\theta }_2^{\ast }\right)-\left(v_{lmhs}^{(2)}-m{\theta }_2\right)\right]K_{\mathrm{R}}^{(2)}\left(l,{\omega }_{lmpq}^{(2)}\right) \\ +\sin \left[\left(v_{lmpq}^{(2),\ast }-m{\theta }_2^{\ast }\right)-\left(v_{lmhs}^{(2)}-m{\theta }_2\right)\right]K_{\mathrm{I}}^{(2)}\left(l,{\omega }_{lmpq}^{(2)}\right)]. \end{gathered}$$

$$U_1\left(-\text{@vec}{r}_{1\to 2},-\text{@vec}{r}_{1\to 2}^{\star }\right)=-\frac{G{m}_2}{a^{\star }}\sum _{l=2}^{\infty }\sum _{m=0}^l\sum _{p=0}^l\sum _{q=-\infty }^{\infty }\sum _{h=0}^l\sum _{s=-\infty }^{\infty }{U}_{lmpqh}^{(1)},$$

$$\begin{gathered} U_{lmpqhs}^{(1)}={\left(\frac{R_1}{a}\right)}^{l+1}{\left(\frac{R_1}{a^{\star }}\right)}^l\frac{(l-m)!}{(l+m)!}(2-{\delta }_{0m})F_{lmp}(i_1^{\star })G_{lpq}(e^{\star })F_{lmh}(i_1)G_{lhs}(e)[\cos \left[\left(v_{lmpq}^{(1),\star }-m{\theta }_1^{\star }\right)-\left(v_{lmhs}^{(1)}-m{\theta }_1\right)\right]K_{\mathrm{R}}^{(1)}\left(l,{\omega }_{lmpq}^{(1)}\right) \\ +\sin \left[\left(v_{lmpq}^{(1),\star }-m{\theta }_1^{\star }\right)-\left(v_{lmhs}^{(1)}-m{\theta }_1\right)\right]K_{\mathrm{I}}^{(1)}\left(l,{\omega }_{lmpq}^{(1)}\right)], \end{gathered}$$

$$v_{lmhs}^{(k)}=(l-2h){\omega }_k+(l-2h+s)M+m{\Omega }_k$$

$$K_{\mathrm{R}}(l,\omega )$$

$$K_{\mathrm{I}}(l,\omega )$$

$${\chi }_{lmpq}^{(k)}=\left|{\omega }_{lmpq}^{(k),\ast }\right|=\left|(l-2p){\dot{\omega }}_k^{\ast }+(l-2p+q){\dot{M}}^{\ast }+m\left({\dot{\Omega }}_k^{\ast }-{\dot{\theta }}_k^{\ast }\right)\right|,$$

$${\omega }_{lmpq}^{(k),\ast }\approx (l-2p+q)n^{\ast }-m{\dot{\theta }}_k^{\ast },$$

$$B_l=\frac{2l^2+4l+3}{l{g}_k{\rho }_k{R}_k},$$

$$\bar{J}(\chi )={\int }_0^{\infty }\dot{J}(t-t^{\prime })\exp \left[-\mathrm{i}\chi (t-t^{\prime })\right]\mathrm{d}{t}^{\prime },$$

$$J(t-t^{\prime })=J(0)\Theta (t-t^{\prime })+\text{mbox}viscousandhereditaryterms,$$

$$2{\bar{u}}_{\gamma \nu }(\chi )=\bar{J}(\chi ){\bar{\sigma }}_{\gamma \nu }(\chi ),$$

$$V_2(\text{@vec}{r}_1)=-\frac{G{m}_1}{a}\sum _{l=2}^{\infty }\sum _{m=0}^l\sum _{p=0}^l\sum _{q=-\infty }^{\infty }{V}_{lmpq}^{(2)},$$

$$V_{lmpq}^{(2)}={\left(\frac{R_2}{a}\right)}^l{F}_{lmp}(i_2)G_{lpq}(e)\left[C_{lm}^{(2)}\cos \left(v_{lmpq}^{(2)}-m{\theta }_2+{\varphi }_{lm}\right)+S_{lm}^{(2)}\sin \left(v_{lmpq}^{(2)}-m{\theta }_2+{\varphi }_{lm}\right)\right],$$

$$C_{20}$$

$$C_{22}$$

$$⟨\frac{\mathrm{d}a}{\mathrm{d}t}⟩=\frac{2}{na}⟨\frac{\partial R}{\partial M}⟩,$$

$$⟨\frac{\mathrm{d}e}{\mathrm{d}t}⟩=\frac{\sqrt{1-e^2}}{n{a}^{2}e}\left[\sqrt{1-e^2}⟨\frac{\partial R}{\partial M}⟩-⟨\frac{\partial R}{\partial {\omega }_1}⟩-⟨\frac{\partial R}{\partial {\omega }_2}⟩\right],$$

$$⟨\frac{\mathrm{d}M}{\mathrm{d}t}⟩=n-\frac{2}{na}⟨\frac{\partial R}{\partial a}⟩-\frac{1-e^2}{n{a}^{2}e}⟨\frac{\partial R}{\partial e}⟩,$$

Full text

*Methods. *We used the equations of motion that give the time evolution of the orbital elements and the spin rates of each body to study the time evolution of the Kepler-21 system as an example of how the formalism that we have developed can be applied.

*Results. *We obtained a set of equations of motion without singularities for vanishing eccentricities and inclinations. This set gives, on one hand, the time evolution of the orbital elements due to the tidal potentials generated by both members of the system as well as the triaxiality of one of them. On the other hand, it gives the time evolution of the stellar spin rate due to the corresponding tidal torque and of the planet’s rotation angle due to both the tidal and triaxiality-induced torques. We found that for the parameters and the initial conditions explored here, the tidally and triaxiality-induced modifications of the tidal modes can be more significative than expected and that the time of tidal synchronization strongly depends on the values of the rheological parameters.

1 Introduction

Recent developments and improvements in the field of exoplanet detection methods have allowed for the discovery of many Earth-sized planets, so-called super-Earths, which have masses between $1$ and $10\,M_{\oplus}$. This has given rise to an interest in other, related topics, such as the dynamical and thermal evolution (Stamenković et al., 2012; Shoji & Kurita, 2014; Makarov et al., 2018; Renaud & Henning, 2018), the presence of atmospheres (Seager, 2010), the internal structure and composition (Valencia et al., 2006, 2007b, 2007a), and the possibility that these types of celestial bodies offer for harboring life.

In this work, we are interested in the first of these problems, that is, the dynamical evolution of a binary system formed by a super-Earth and its host star, modeled as bodies that are capable of experiencing internal stresses and energy dissipation. Our aim is to develop a set of equations of motion to study the orbital and rotational evolution of such a system due to the tides that each celestial body induces on the other and the triaxiality of the planet itself. We are particularly interested in the dynamical evolution near spin orbit-resonances. For this reason, we have derived the equations of motion that give, on one hand, the time evolution of the orbital elements due to the secular terms of both the tidal and the $J_{2}$-related disturbing potentials and, on the other hand, the time evolution of the stellar rotation rate and that of the rotation angle of the planet, which are both due to the secular terms of the corresponding tidal torque and to the triaxiality-induced torque on the latter.

Certainly, we expect that close-in planets ends their rotational history in the 1:1 spin-orbit resonance, that is, when $\dot{\theta}=n,$ where $\dot{\theta}$ is the spin rate of the considered body and $n$ is the the mean orbital frequency, which, in the unperturbed two point-masses problem, is defined by the mathematical expression of the Kepler’s third law:

[TABLE]

where $G$ is the universal gravitation constant. Yet, as we shall see later, a rotating body may get trapped in higher-than-synchronous rotation states (Efroimsky, 2012) and that the eccentricity plays a key role in the possibility of capture in spin-orbit resonances, as shown by Goldreich & Peale (1966); Efroimsky (2012); Makarov (2012); Makarov et al. (2012); Noyelles et al. (2014).

Resonant capture in spin-orbit resonances is a phenomenon that occurs due to the interplay between the tidal and the triaxiality-induced torques. The former torque is due to the tides that each component of the binary system induces on the other, while the latter is characterized by the principal moments of inertia of at least one of them. Naturally, here, we consider super-Earths as triaxial bodies because they probably have a rocky composition, in contraposition to the stars, which are expected to be fluid.

This work is organized as follows: In Sect. 2, we present the equations of motion in its general form as were derived by Boué & Efroimsky (2019). We also explicitly give the tidal and triaxial potentials and the rheological models for each component. In Sect. 3, we derive the equations that have to be solved numerically in order to obtain the time evolution of the particular system under study, which have the particularity of being regular for very small and null eccentricities and inclinations. In Sect. 4, we discuss the choice of the values of the physical and rheological parameters involved and the initial conditions for the numerical integration. In Sect. 5, we study the dynamical evolution of the Kepler-21 system as an example of application of the developed methodology. Also, we present and discuss the results obtained from the numerical integration of the equations of motion mentioned above and evaluate the possibility of detecting tidal effects. Finally, in Sect. 6, we offer a final discussion and draw our conclusions.

2 Generalities

To describe the dynamical evolution, including the orbital and rotational motions, of a binary system constituted of a rocky planet and its host star, we make use of the equations of motion derived by Boué & Efroimsky (2019), which are valid for the case when the tidal dissipation within both bodies is considered, including the disturbing potential due to the permanent-triaxial shape of the planet.

We consider a system of two bodies, the host star and a rocky planet. Each body is thought of as a nearly spherical mass distribution characterized by its principal moments of inertia, total mass, and mean radius. The former is assumed to have total mass, $m_{1}$, and radius, $R_{1}$, and the latter has the total mass, $m_{2}$, and radius, $R_{2}$. The corresponding principal moments of inertia of each body is represented by $A_{k}$, $B_{k}$, and $C_{k}$ (with $k=1,2$ and $A<B<C$). The index, $k$, identifies each body: $k=1$ corresponds to the star-related parameters and variables, while $k=2$ corresponds to the planet-related ones.

The mutual orbit is parameterized by the six classical orbital elements: the major semiaxis ($a$), the eccentricity ($e$), the mean anomaly ($M$), the argument of pericenter ($\omega$), the inclination ($i$), and the longitude of the ascending node ($\Omega$). However, as we allow for the tidal dissipation within both bodies, then the two reference systems corresponding to each one must be considered in order to evaluate the forces each body exerts on the other. In other words, when the force of the planet on the star is to be evaluated, we have to consider the orbit of the former as seen from the latter, and vice versa when evaluating the force of the star on the planet. As both orbits are defined with respect to the corresponding equatorial planes of each body, then the aforementioned reference systems are related by a rotation from one system to an inertial reference system and then to the system attached to the other body. However, this rotations only affect the angles that give the orientation of the orbits in the space, that is, $\omega$, $i$ and $\Omega$, while $a$, $e,$ and $M$ are the same in both systems (Boué & Efroimsky, 2019).

We assume that each body rotates about its maximal inertia axis. The rotational state of the star and the planet are completely determined by the corresponding Euler angles, that is, the precession, $\psi_{k}$, nutation, $\varepsilon_{k}$, and proper rotation angles, $\theta_{k}$ – the latter being reckoned from a fixed direction in an inertial reference system – as well as their time rates $\dot{\psi}_{k}$, $\dot{\varepsilon}_{k}$ and $\dot{\theta}_{k}$, respectively. As we see later in this paper, under the gyroscopic approximation, it is enough to follow the time evolution of $\varepsilon_{k}$ and $\theta_{k}$ (Boué & Efroimsky, 2019).

2.1 Orbital evolution

The orbital dynamics is described by a set of Lagrange-type planetary equations which give, on one hand, the time evolution of the major semiaxis, the eccentricity, and the mean anomaly along with, on the other hand, the time evolution of the inclination, the argument of the pericenter and the longitude of the ascending node of the common orbit as seen from each body. The respective equations were derived by Boué & Efroimsky (2019), which correspond to Eqs. (116)-(121) of the cited work and we reproduce them here:

[TABLE]

where $\beta=m_{1}\,m_{2}/\left(m_{1}+m_{2}\right)$ is the reduced mass, $R$ is the disturbing function, $\omega_{1}$, $i_{1}$, and $\Omega_{1}$ are, respectively, the argument of pericenter, the inclination, and the longitude of the ascending node of the mutual orbit with respect to the star’s equatorial plane, and $\omega_{2}$, $i_{2}$ and $\Omega_{2}$ are the corresponding elements defined with respect to the planet’s equator.

As we are considering a system of two bodies, there is only one orbital plane and the inclinations of the equatorial plane of each body with respect to the orbital plane are, indeed, the obliquities. But, as pointed out by Boué & Efroimsky (2019), these quantities are called inclinations, in terms of the continuity of Kaula’s notation. On the other hand, the names of obliquity and inclination of the orbital plane depends on the problem that is being considered. In our case, what would be called inclination of the planet’s orbit as seen from the star ($i_{1}$) is simply the orbital inclination and, conversely, the inclination of the orbit of the star as seen from the planet ($i_{2}$) can simply be referred to as the obliquity of the planet.

It worthwhile to indicate that the orbital elements parameterizing the positions of each body with respect to the other are not osculating. As it is widely known, in the unperturbed two-body problem, the position, $\@vec{r,}$ and the velocity, $\dot{\@vec{r}}$, vectors are functions of the orbital elements, that is: $\@vec{r}=\@vec{f}(a,e,M,i,\omega,\Omega)$ and $\dot{\@vec{r}}=\@vec{g}(a,e,M,i,\omega,\Omega)$ being, by definition,

[TABLE]

where $c_{k}=a,e,i,M,\omega,\Omega$ for $k=1,\ldots,6$. Under perturbation, the position vector remains equal to $\@vec{f}$ but the velocity vector reads:

[TABLE]

where $\@vec{\Phi}$ is known as the gauge function and it is completely arbitrary (Efroimsky, 2005a, b; Kopeikin et al., 2011). In other words, we can fix a particular functional form of the convective term of Eq. (4):

[TABLE]

The usual practice is to set $\@vec{\Phi}=\@vec{0}$, which is also called the Lagrange’s constraint. In such a case, we have $\dot{\@vec{r}}=\@vec{g}$, meaning that the perturbed velocity is tangential to the conic parameterising the instantaneous position of each celestial body with respect to the other. Of course, the corresponding orbital elements are the osculating ones.

However, within the formalism developed by Boué & Efroimsky (2019), the orbital elements are defined in the equatorial planes of each body, that is, they choose to describe the relative motions of the bodies in the coprecessing reference frames whose origins coincide with the center of mass of each body. In such situations, it is more convenient to choose $\@vec{\Phi}=\@vec{\mu}\times\@vec{r}$, where $\@vec{\mu}$ is the precession velocity vector (Efroimsky & Goldreich, 2004; Efroimsky, 2005a, b). Thus, the velocity vector of each body is given by $\dot{\@vec{r}}=\@vec{g}+\@vec{\mu}\times\@vec{f}$. Later on, we can perform the transformation from the set $(\@vec{r},\dot{\@vec{r}})$ to the set $(a,e,M,i,\omega,\Omega)$ to get the osculating elements, which have a more direct interpretation. The orbital elements appearing in Eqs. (2a) – (2i) are called contact elements.

2.2 Rotational evolution

In order to close the system of equations of motion, at least, in principle, it is necessary to add the following equations describing the rotational evolution, which correspond to Eqs. (122) and (123) of the work by Boué & Efroimsky (2019):

[TABLE]

[TABLE]

where $k=1,2$, $\varepsilon_{k}$ is the angle between each body’s spin axis and the $z$ direction of the inertial reference frame, and, again, $\theta_{k}$ are the sidereal angles. It should be remembered that the last equations and Eqs. (2a) – (2i) are obtained under the so-called gyroscopic approximation. This means that $\left|\dot{\theta}\right|\gg\left|\dot{\psi}\right|\,\mathrm{,}\,\left|\dot{\varepsilon}\right|$ is assumed. The latter assertion is equivalent to stating that the rotational angular momentum vector coincides with the principal axis of inertia the body rotates around.

Another important and necessary assumption is that we also neglect the time variation of the moments of inertia of the considered bodies. This is justified because even though the tidal theory considers the deformation of celestial bodies – as we see later – these deformations (assumed to be linear) are always very small, that is, $\epsilon<10^{-6}$, where $\epsilon$ is a typical value of the strain (Frouard & Efroimsky, 2018).

2.3 Disturbing potentials

The disturbing function, $R$, present in the equations of motion is given by (Boué & Efroimsky, 2019):

[TABLE]

where $U_{2}\left(\@vec{r}_{2\rightarrow 1},\@vec{r}_{2\rightarrow 1}^{*}\right)$ is the tidal disturbing potential generated by the (deformed) planet and $U_{1}\left(-\@vec{r}_{1\rightarrow 2},-\@vec{r}_{1\rightarrow 2}^{\star}\right)$ is the corresponding tidal potential generated by the star. As we are interested in the dynamical evolution near the spin-orbit resonances, we add the permanent-triaxiality-induced potential of the planet $V_{2}\left(\@vec{r}_{2\rightarrow 1}\right)$, as well as in the influence of the dynamical form factor on the secular evolution of the system. The first one is given by (Efroimsky, 2012):

[TABLE]

where

[TABLE]

The tidal potential generated by the star is given by an expression analogous to Ec. (9) and (10):

[TABLE]

with

[TABLE]

where $\delta_{0m}$ is the Kronecker’s delta, $F_{lmp}(i_{k})$ and $G_{lpq}(e)$ are the inclination functions (Gooding & Wagner, 2008) and the eccentricity functions (Giacaglia, 1976), respectively. The $v_{lmpq}^{(k)}$ represent linear combinations of $\omega_{k}$, $\Omega_{k}$ and $M$, expressed via:

[TABLE]

and analog expressions for $v_{lmpq}^{(k),*}$ and $v_{lmpq}^{(k),\star}$.

Within the framework of the Darwin-Kaula formalism of bodily tides, we may consider the presence of three bodies in order to take the tidal dissipation into account: the primary – whose shape and deformation is considered – and two “secondaries,” the tide-rising and the tidally disturbed ones, considered to be point masses. Let $\@vec{r}_{2\rightarrow 1}$ be the position vector of the body 1 (the star) with respect to the body 2 (the planet), and $\@vec{r}_{1\rightarrow 2}$ the corresponding position vector of the planet with respect to the star, both expressed in the corotating frame of one of the bodies. In the case of the planet being the primary, the tide-rising secondary is located at $\@vec{r}_{2\rightarrow 1}^{*}$, while the tide disturbed one is at $\@vec{r}_{2\rightarrow 1}$. Conversely, when the star is considered to be the primary, then the tide-rising secondary is at $-\@vec{r}_{1\rightarrow 2}^{\star}$ while the tide-disturbed one is at $-\@vec{r}_{1\rightarrow 2}$. It worthwhile indicating that as we are considering a two-body system, in both cases the two secondaries coincide with one other, that is, in the first case, both secondaries coincide with the star and, in the second case, both coincide with the planet. By virtue of the Kaula (1961) transformation, the set of coordinates for each secondary is expressed through the corresponding set of orbital elements that parameterizes their orbits (Efroimsky & Williams, 2009).

In the above expressions – Eqs. (10) and (12) – $K_{\mathrm{R}}^{k}(l,\chi^{k}_{lmpq})$ and $K_{\mathrm{I}}^{k}(l,\chi^{k}_{lmpq})$ are known as quality functions (Makarov, 2012) and given by:

[TABLE]

where $\Re(\bar{z})$ and $\Im(\bar{z})$ denotes the real and imaginary parts of the complex number $\bar{z}$ and $\mathrm{sign}(x)$ is the signum function. We omitted the indices for the sake of simplicity. Then, $\chi^{(k)}_{lmpq}$ are the tidal forcing frequencies, that is, they are the physical frequencies at which the deformations are excited. They are defined as the absolute value of the tidal modes, $\omega^{(k),*}_{lmpq}$, (Efroimsky, 2012; Efroimsky & Makarov, 2013), mathematically:

[TABLE]

where $\dot{\theta}^{\,*}_{k}$ is the spin rate of the primary that is being considered. If both the apsidal and the nodal precessions are neglected, the tidal modes can be expressed as:

[TABLE]

where $n^{*}_{k}$ is the mean orbital frequency of the tide-rising secondary. Consequently, the $\omega^{(k),*}_{lmpq}$’s can be either positive or negative depending on whether the primary’s spin rate is slower or faster than the mean rate at which the tide-rising secondary describes its orbit, respectively, while the tidal frequencies are always positive. $B_{l}$ is given by (Efroimsky, 2015):

[TABLE]

where $\rho_{k}$ is the mean density, $R_{k}$ is the mean radius, and $g_{k}=Gm_{k}/R_{k}^{2}$ is the surface gravity of the body, $k$. Finally, $\bar{J}(\chi)$ is the complex compliance and is defined by:

[TABLE]

where the over-dot means differentiation with respect to $t^{\prime}$. The functional form of the kernel $J(t-t^{\prime})$ depends on the rheology of the material that composes the primary or, in other words, on the rheological model that describe the creep response of the body as a whole. Realistic materials are thought to have a mixture of elastic, viscous (or viscoeslastic), and hereditary behaviors (Renaud & Henning, 2018), that is, the deformation is directly proportional to the stress and simultaneously depends on the time derivative of the latter and on its past values, respectively. Thus, $J(t-t^{\prime})$ is given generally by:

[TABLE]

with $J(0)$ as the value of the unrelaxed compliance, corresponding to the elastic response and inverse of the elastic rigidity, $\mu(0)$, and $\Theta(t-t^{\prime})$ is the Heaviside step function (Efroimsky, 2012).

It worthwhile indicating that the complex compliance enters in the constitutive equation that relates the stress and strain tensors together. In a linear regime, and under the assumption that the material the primary is made of is an homogeneous incompressible and isotropic medium, the relationship between the components of the tensors, written in the frequency domain, is:

[TABLE]

where $\bar{u}_{\gamma\nu}$ and $\bar{\sigma}_{\gamma\nu}$ are the complex counterparts of the components of the strain and stress tensors, respectively. We refer the reader to the work by Efroimsky (2012) for a comprehensive explanation of this subject and detailed derivations of the equations given above.

The permanent-triaxiality-induced potential is a special case of the more general expression of the potential generated by the so-called coefficients of the gravitational field (CGF), $C_{lm}$, and $S_{lm}$, which is given by (Frouard & Efroimsky, 2017):

[TABLE]

where

[TABLE]

and where $\phi_{lm}=\left[(-1)^{l-m}-1\right]\frac{\pi}{4}$. The leading terms of the expansion given in Eqs. (22) and (23) are those related with the triaxiality of the planet, namely the terms associated with the $C_{lm}$ coefficients for which $l=2$ and $m=0,2$, specifically, $C_{20}$ and $C_{22}$. The latter are related with the aforementioned principal moments of inertia through:

[TABLE]

where $J_{2}$ is known as the dynamic form factor, which usually is three orders of magnitude larger than the rest of the CGF, at least for the case of the Earth (Torge & Müller, 2012). In this work, we consider the influence of the $J_{2}$-related secular terms on the time evolution of the orbital elements. The $C_{22}$ coefficient is kept because is the dominating factor of the triaxiality-induced torque (Makarov, 2012; Makarov et al., 2012) which, together with the tidal torque, drive the time evolution of the proper rotation angle of the planet. Frouard & Efroimsky (2017) point out that the latter torque is also responsible for the longitudinal libration, which causes extra dissipation. However this effect is out of the scope of this work and we will consider it in the future.

We only consider the triaxiality of the planet for two main reasons: 1) The tidal interaction is expected to affect more strongly its dynamical evolution than that of the star, in particular the rotational evolution where there are probabilities related to the planet getting trapped in a spin-orbit resonance; 2) We can reasonably neglect the triaxiality (i.e. the deviation from a rotational symmetry) of a fluid body, such as a star. In that case, we assume that $A_{1}=B_{1}$ and, thus, the $C_{22}$-related terms of the gravity potential of the star are nullified, as can be confirmed by looking at Eq. (24). As only the triaxiality of the planet is to be considered, we can safely suppress the corresponding subscripts of the principal moments of inertia, $A$ and $B$.

Regarding the $J_{2}$-related terms of the corresponding CGF disturbing potential, even though the oblateness of the star could have a relevant impact on the dynamics of close-in systems, such as the one studied in this work (Correia et al., 2011), we have left the consideration of this effect for a future work.

3 Development of the analytical model

The principal considerations that guides the development of our analytical model consist of the following:

• •

As a first step, we consider the secular evolution of the orbital elements. In this sense, our model takes into account the secular terms of the derivatives of the additional tidal potential and the $J_{2}$-related secular terms of the CGF disturbing potential.

• •

Concerning the time evolution of the star’s spin rate, we consider that its angular acceleration is caused by the secular terms of the tidal torque acting on it.

• •

We then consider that the time evolution of the planet’s spin rate is dominated by both the secular terms of the corresponding tidal and permanent-triaxiality-induced torques.

• •

Finally, we limit the expansions of the disturbing potentials to the quadrupole approximation, i.e. $l=2$. A justification of this assumption is given in Sect. 4.

The equations describing the orbital and rotational dynamics given in Eqs. (2a)–(2i), (6), and (7) are somewhat general and, therefore, we find it necessary at this point to rewrite them in a form that is more convenient to our purposes. As we consider the secular evolution of the orbital elements, we make use of the averaged equations of motion:

[TABLE]

In this context, the angular brackets on the left-hand side of Eqs. (26) represent the averaging over the mean anomaly, the arguments of pericenters, and the longitudes of the ascending nodes, while the angular brackets on the right-hand side of the aforementioned equations indicate that we are considering the secular terms of the disturbing function’s derivatives.

Since the equations of motion depends on the gradient of the disturbing function, $R$, among other factors, the calculation of the derivatives of $R$ with respect to the orbital elements deserves some discussion. As we point out in the previous subsection, $\@vec{r}_{2\rightarrow 1}$ correspond to the position of the tide-disturbed secondary while the tide-rising one is located at $\@vec{r}_{2\rightarrow 1}^{*}$ when the planet is considered to be the primary body. Accordingly, $-\@vec{r}_{1\rightarrow 2}$ and $-\@vec{r}_{1\rightarrow 2}^{\star}$ are the radius-vectors of the tide-disturbed and tide-rising secondaries, respectively, when the roles of the planet and the star are reversed with respect to the former case. It should be stressed that we always consider the motion of the tide-disturbed secondary and, thus, the general rule that is to be followed in order to find the aforementioned derivatives is the following: first, we have to differentiate the disturbing function with respect to the orbital elements of the tide-disturbed secondary and then identify both secondaries with each other, that is, to set $\@vec{r}_{2\rightarrow 1}=\@vec{r}_{2\rightarrow 1}^{*}$ and $\@vec{r}_{1\rightarrow 2}=\@vec{r}_{1\rightarrow 2}^{\star}$ in order to properly evaluate the force per unit mass each body exerts on the other. It may be noted that as the arguments of the sine and cosine functions remain the same, once the identification is performed, we obtain:

[TABLE]

By virtue of Eq. (27), it can be seen that the derivatives of the tidal potentials may be split into two types of terms, namely, the secular and the oscillating ones. The former comprises the terms independent of $M$ and $\omega_{k}$, that is, those for which $h=p$ and $q=s$, with the remainder being part of the latter type of terms. As we can easily see, for the case of the secular terms, the arguments of the sine and cosine given by Eq. (27) are identically null and, as a consequence, the derivatives of the disturbing function become independent of one of the quality functions.

It can be shown by Eq. (23) that the $J_{2}$-related secular term of the CGF disturbing potential corresponds to the one with $(lmpq)=(2010)$ (Kaula, 2000; Efroimsky, 2005b), that is:

[TABLE]

In general terms, analogously to Eqs. (9), (11), and (22), the disturbing potentials can be expressed as the product of a dimensional factor multiplied by a sum of adimensional terms. Taking into account the dependence of each perturbing potential upon the corresponding orbital elements, which can be expressed through:

[TABLE]

the derivatives of the disturbing function can be calculated as follows. Since the three disturbing potentials considered in this work depend on $a$, $e,$ and $M$, we have, on one hand:

[TABLE]

while, on the other hand, the corresponding expression of the derivatives with respect to the inclinations, the arguments of the pericenter and the longitudes of the ascending nodes (as seen from each coordinate system) are:

[TABLE]

Specifically, the secular terms of the derivatives of the disturbing potentials with respect to the major semiaxis can be expressed as:

[TABLE]

where $\left\langle U_{lmpq}^{(1)}\right\rangle$ and $\left\langle U_{lmpq}^{(2)}\right\rangle$ are the secular terms of the tidal potentials of both the star and the planet, respectively. These terms are obtained through Eqs (10) and (12) by setting $h=p$ ans $q=s$, as previously explained. Even though we have pointed out before that we aim to limit the expansions to $l=2$, we still consider all multipoles in the expressions of the derivatives of the disturbing potentials as a matter of convenience in order to ease the exposition. Later, when we give the final expression of the equations of motion, we provide explicit expressions that correspond to the quadrupolar approximation.

The expressions of the derivatives of the disturbing potentials with respect to the eccentricity and the mean anomaly are:

[TABLE]

The rest of the derivatives, those with respect to $i_{k}$, $\omega_{k}$, and $\Omega_{k}$, can be expressed as follows:

[TABLE]

The full expressions of the derivatives of the adimensional terms are given in Appendix A.

The next step is to write down the expressions of the derivatives of the disturbing function with respect to the orbital elements. By inserting Eqs. (31a), (31b), (31c), (32a), (32b), and (32c) in Eq. (29) we have:

[TABLE]

The partial derivatives of the disturbing function with respect to the inclinations, arguments of pericenters and longitudes of the corresponding ascending nodes can be expressed in a similar way. Inserting Eqs. (33a), (33b), and (33c) again in Eqs. (30), we obtain:

[TABLE]

where

[TABLE]

By inserting Eqs. (34) and (35) in Eq. (26) we obtain the equations describing the orbital motion. If we also take into account Eq. (28) we obtain:

[TABLE]

With regard to the rotational evolution, since Eqs. (37) depend neither on $\varepsilon_{1}$ nor $\varepsilon_{2}$, we only need Eqs. (7) to close the set of equations describing the dynamical evolution of binary systems such as the ones we consider in this work. However, if we insert the expression of the derivative of the disturbing function with respect to the corresponding longitudes of the ascending nodes, given by Eq. (35), we obtain just the equations for the angular accelerations due to the secular terms of the tidal potentials, that is, those of the tidal torques. That equation is enough to describe rotational evolution for the star, but not for the planet. Since the latter could cross one or several spin-orbit resonances, we should take into account the permanent-triaxiality-induced torque. This torque arises from the $C_{22}$-related terms of the CGF disturbing potential, which are given by:

[TABLE]

where

[TABLE]

Thus, the equations governing the time evolution of the spin rate of both the star and the planet, respectively, are obtained as follows: By inserting in Eq. (7) the expression of the corresponding derivative, given by Eq. (35a), we arrive at:

[TABLE]

3.1 Simplifications available for low inclinations and eccentricities

Up to this point, we have presented the analytical tools in a general form which allow for a description of the dynamical evolution of a binary system that considers the tidal deformation of both of its components and the triaxiality of one of them. Based on the observation that the disturbing potentials and their corresponding gradients – which build up both the disturbing function and its derivatives, respectively – are represented as an infinite series, it becomes necessary to examine the leading terms in order to perform the mandatory truncation that allows for the implementation of the formulae in a computer code for its practical application. To that end, we take into account the quadrupolar approximation mentioned at the beginning of Sect. 3.

We now focus on the inclination functions. These are defined by (Gooding & Wagner, 2008):

[TABLE]

where the limits of summation are $\max(0,l-m-2p)$ and $\min(l-m,2(l-p))$. In Table 1, where we show the inclination functions for $l=2$, it can be seen that for low inclinations $F_{201}(i)$ and $F_{220}(i)$ tend towards $1/2$ and $3$, respectively, as $i$ tends towards zero while all the others tend towards zero. In consequence, the terms corresponding to $(lmpq)=(201q)$ and $(lmpq)=(220q)$ of the disturbing potentials dominate the dynamical evolution in the case of low inclinations. For this reason, the latter are given in Appendix A.

At this point, it is necessary to briefly discuss on the truncation of the sums over $q$ in order to make practical computation possible. In general terms, the eccentricity functions are defined as:

[TABLE]

where

[TABLE]

where $\nu$ is the true anomaly. Thus, the eccentricity functions are defined as the coefficients of the expansion in Fourier series of the left-hand side of Eq. (42). Of course, the derivatives of the eccentricity functions are also needed. These can be defined in an analogous way as the coefficients of the Fourier series expansion of the derivative of the left-hand side of Eq. (42) with respect to the eccentricity (Giacaglia, 1976):

[TABLE]

where

[TABLE]

and

[TABLE]

The expressions corresponding to $(lmpq)=(201q)$ are:

[TABLE]

while those corresponding to $(lmpq)=(220q)$ are:

[TABLE]

The question of how many terms are to be kept in order to obtain a reliable representation of the functions on the left-hand side of Eqs. (47), (48), (53), and (60) can be answered by comparing them with their corresponding Fourier series representation. The criterion adopted in this work is the computation of the root mean squared difference (RMSD) between the left and right-hand sides of the latter equations corresponding to different values of $q_{\mathrm{max}}$, which is the maximum number of coefficients. Mathematically expressed, let $g(e,\nu)$ be a function that represents each of those on the left-hand side of Eqs. (47) and (48) as well as the real and imaginary parts of Eqs (53) and (60). Also, let $\tilde{g}(e,M)$ be the function representing the corresponding right-hand side of the same equations. The RMSD is given by:

[TABLE]

where $N$ is the number of values of $M$ and $\nu$ within the interval $[-\pi,\pi]$, each of them labeled with $M_{i}$ and $\nu_{i}$ being the respective value of the true anomaly.

Of course, the number of terms of each of the Fourier series expansion strongly depend on the eccentricity. As can be observed in Table 2, for $e=0.1$, there is no significative reduction of the RMSD for $q_{\mathrm{max}}>10$. We have considered the maximum value of the eccentricity from the estimated uncertainty for the Kepler-21 system. Thus, we consider that it is safe to take $q_{\mathrm{max}}=10$. It is worth mentioning that the eccentricity functions and its first derivatives with respect to the eccentricity were computed by using the Newcomb operators, as explained in Appendix B.

3.2 Final form of the equations of motion

The equations governing the orbital evolution, Eqs. (37) can be rewritten in an even more convenient form that could be suitable for its implementation in a computational code specially for low (or even zero) eccentricities and inclinations. As pointed out by Brouwer & Clemence (1961), the divergences in the aforementioned equations “are more apparent than real” and that observation remains valid in this context. First, it can be noted that in Eqs. (37f), (37g), (37h), and (37i) each derivative of the tidal potentials, generated by the star and by the planet, with respect to the inclination, appear divided by the sine of the corresponding inclination. By observing the expressions of the latter derivatives, given by Eqs. (85b), then we see that the inclinations functions appear multiplied by its derivative. In consequence, observing the third column of Table 1, we also note that there are no singularities in the respective terms of the aforementioned equations, provided that the dividing sine of the inclination is canceled.

Something analogous occurs with the eccentricity dividing the derivatives of the tidal potentials with respect to the eccentricity. We can observe in Eqs. (83b) that each term of the aforementioned derivatives contains products of the corresponding eccentricity function and its derivative. It can be shown that the these factors, $G_{2pq}(e)\mathrm{d}G_{2pq}(e)/\mathrm{d}e$, do not contain terms of order $e^{0}$.

We can introduce further modifications to Eqs. (37b), (37d) and (37e). In order to avoid divergences that are also present in the equations mentioned previously, we define the following variables: $\xi=\sqrt{1-e^{2}}$ and $x=\cos i$. Therefore, we have:

[TABLE]

[TABLE]

Thus, the final form of the equations describing the orbital evolution, upon considering the simplifications discussed before, are:

[TABLE]

Concerning the time evolution of the spin rate of both the star and the planet, the respective equations of motion take the form:

[TABLE]

Evidently, the lack of angular brackets in the last equation is due to the presence of oscillatory terms originated by the triaxiality-induced torque, which does not average out in the neighborhood of spin-orbit resonances.

As can be observed, the equations of motion thus obtained contain no singularities in the limit of vanishing eccentricities and inclinations. However, we must stress the fact that this particular form of the equations of motion is intended to be translated to a computer code. In this sense, we were interested in the elimination of the singularities present in the equations for the time rate of the eccentricity and of the inclinations – by defining the variables $\xi$ and $x$ – in order to avoid possible instabilities during the numerical integrations. By no means would we pretend to assert that through the aforementioned transformation, the arguments of pericenters and the longitudes of the ascending nodes become non-singular. These parameters still become meaningless when the eccentricities or the inclinations or both are zero. The only non-singular set of variables are the Cartesian-like ones (Correia et al., 2014; Boué et al., 2016).

Regarding the secular evolution of the system, only the values of the major semiaxis, the eccentricity, the inclinations, the argument of the pericenter, and the longitude of the ascending node defined on the planet’s equator, as well as the rotation angle of the latter, are needed. The time derivatives of $M$, $\omega_{1}$, $\omega_{2}$, $\Omega_{1}$ and $\Omega_{2}$ are necessary to compute the tidal forcing frequencies by virtue of Eq. (16). It is a normal practice to directly use Eq. (17), but as we show later in this work, this approximation is not always valid.

We can summarize the principal characteristic of the dynamical model proposed, whose mathematical expression is given by Eqs. (64) and (65), as follows:

• •

The major semiaxis, the eccentricity and the inclinations evolve in time strictly due to tidal effects as well as the longitude of the ascending node as seen from the star.

• •

The time evolution of the mean anomaly, the arguments of pericenters, and the longitude of the ascending node of the orbit seen from the planet is dominated by both the tidal interaction between both bodies and the $J_{2}$-related secular terms of the gravitational field of the planet.

• •

The stellar spin rate is dominated only by the secular terms of the corresponding tidal torque.

• •

The spin rate of the planet is driven by the combined action of the secular terms of the tidal and triaxiality-induced torques.

4 Preliminary discussion

In this section, we briefly present the exoplanetary system studied in this work and we discuss the rheology of each component and the values of physical and rheological parameters we assumed for it, as well as the initial conditions for the numerical simulation we perform in this study, especially for those concerning the rotational evolution of the planet.

4.1 The Kepler-21 system

The Kepler-21 star, also known as HD 179070, is a F6IV sub-giant located at a distance of 352 light-years in the Lyra constellation (Howell et al., 2012). In Table 3, we show the orbital, planetary, and stellar parameters relevant to our work. From the mass and radius of the planet, it can be inferred that it has rocky composition, with a mean density of approximately $6.4\pm 2.1\,\mbox{g}\,\mbox{cm}^{-3}$ (López-Morales et al., 2016). In Fig. 1, we show the mass-radius relationships according to several ice-mass fraction values (which is a measure of the water content of the planet) and considering a planetary composition similar to that of the Earth (Valencia et al., 2007a). It also shows the corresponding position of some super-Earths on the graph. The Earth would be placed at the origin. As may be noted, Kepler-21b is expected to have an Earth-like composition with about 10 % made up of water.

We selected this system based on two main reasons: the proximity of the components suggest that tidal effects are strong at least on the planet and all the relevant physical and orbital parameters are well-determined.

4.2 The quadrupolar approximation

We mentioned in Sect. 3 that we want to limit the expansion of the tidal potentials of both the star and the planet up to $l=2$. Despite the fact that the quadrupolar approximation is widely used in the literature of tidal dynamics, such as in studies of the Earth-Moon system (Mignard, 1979, 1980), the Mars-Phobos system (Efroimsky & Lainey, 2007), and of close-in binary systems (Correia et al., 2011, 2014), we still need to verify if it is secure to use that approximation for the case of the Kepler-21 system as well.

Taking into account Table 3, we have for the Kepler-21 system: $R_{1}/a\approx 0.21$ and $R_{2}/a\approx 0.0016$. The first of these values seems to be somewhat high. However, we have to keep in mind that in order to neglect the $l=3$ term we have to compare it with the $l=2$ one. Recalling that the $l$-th term of the tidal potential scales as $(R/a)^{2l+1}$, starting with $l=2$ as is shown by Eqs. (10) and (12), for the case of the star we have $(R_{1}/a)^{5}\approx 3.8\times 10^{-4}$ and $(R_{1}/a)^{7}\approx 1.63\times 10^{-5}$; while for the case of the planet: $(R_{2}/a)^{5}\approx 1.71\times 10^{-14}$ and $(R_{2}/a)^{7}\approx 3.13\times 10^{-20}$. Thus, it seems that it could be necessary to consider terms up to $l=3$ in the expansion of the tidal potential of the star. However, in order to avoid mathematical and numerical complications in our model, we still limit ourselves to the quadrupolar approximation and consider that the order of magnitude of difference between the $l=2$ and the $l=3$ terms of the tidal potential of the star is enough to neglect the latter.

4.3 Rheology of the components

Because the tidal theory we adopt in this work is applicable to binary systems, we find it useful to present a discussion of the rheology of each component. In particular, we consider those that are conformed by a rocky planet and its host star. The reason we limit our study to Earth-like planets is twofold: On one hand, most recent developments in the field of planetary detection methods and their accuracy have led to the discovery of such worlds and, consequently, there is a growing interest in researching whether these celestial bodies offer the conditions that makes life possible. On the other hand, if we can infer that the compositions of the studied planets are similar to that of the Earth, then we can assume that their rheology is also similar and then propose reliable models to describe its internal dissipation.

In Table 3, we show the values of the physical parameters involved in the simulations we performed, that is, the mass and the radius of both components together with those of the major semiaxis and the eccentricity of the system we have selected for applying the model described in Sect. 3, as well as the rotational period of the star and the orbital period of the planet. In what follows, we first discuss the rheology of the host star and then the rheology of the planet.

4.3.1 Rheology of Sun-like stars

As discussed by Veras et al. (2019) stars can be considered to be viscous bodies, despite the fact that they can present some magnetic rigidity. Therefore, the dominant rheological parameter is the shear viscosity $\eta_{1}$. From the work by Efroimsky (2012) we know that the complex compliance corresponding to viscous deformation is given by:

[TABLE]

Then, we have that its real and imaginary parts evidently are:

[TABLE]

As a consequence, the quality functions of the star $K_{\mathrm{I}}^{(1)}\left(l,\omega^{(1)}_{lmpq}\right)$ and $K_{\mathrm{R}}^{(1)}\left(l,\omega^{(1)}_{lmpq}\right)$, given in general form by Eq. (14) and (15), respectively, become (Efroimsky, 2015):

[TABLE]

4.3.2 Rheology of Earth-like planets

The rheological response of rocky and icy celestial bodies like planets and satellites – including the possibility that they can have a subsurface ocean – is best modeled by the Maxwell-Andrade model, or even better, perhaps, by the Sundberg & Cooper (2010) model, as asserted by Renaud & Henning (2018). Due to the complexity of the latter model and the fact that, in general, the rheological parameters are poorly constrained even for the Earth, we chose the former rheological model to describe the dissipation within the planet. The complex compliance for this rheology is (Efroimsky, 2012, 2015):

[TABLE]

where $\alpha$ is known as the Andrade parameter, $\Gamma$ is the gamma function, $\tau_{\mathrm{M}}$ and $\tau_{\mathrm{A}}$ are the Maxwell and Andrade characteristic time scales, respectively. It worthwhile indicating that $\tau_{\mathrm{M}}$ is defined as the ratio of the shear viscosity over the mean rigidity $\mu$, that is:

[TABLE]

It should be noted that the first two terms in the first line of Eq. (69) correspond to a viscoeslastic response which is described by the Maxwell rheology, while the last term corresponds to an hereditary reaction caused by dislocation process within the material under stress (Efroimsky, 2012), as prescribed in Eq. (20). Andrade’s rheology can be applied to a wide range of materials, from rocky planets or satellites partially molten up to icy bodies, even if these bodies have a subsurface ocean. The Andrade’s parameter ranges between $[0.15,0.2]$ for the former kind of celestial bodies and between $[0.2,0.4]$ for the latter.

The real and imaginary parts of the complex compliance given by Eq. (69) are:

[TABLE]

and

[TABLE]

Then, again by the virtue of Eqs. (14) and (15), the quality functions of the planet $K_{\mathrm{R}}^{(2)}\left(l,\omega^{(2)}_{lmpq}\right)$ and $K_{\mathrm{I}}^{(2)}\left(l,\omega^{(2)}_{lmpq}\right)$, respectively, are given by:

[TABLE]

4.4 Initial conditions and physical parameters

Apart from the values of the parameters given in Table 3 we must, on one hand, set some other physical and rheological parameters and, on the other hand, choose the values of the initial conditions of the orbital and rotational parameters in order to solve the differential equations of motion.

With regard to the rheological parameters of the host star and the planet, we only have to set the value of the shear viscosity $\eta_{1}$ of the former and the shear rigidity $\mu$, the Maxwell time $\tau_{\mathrm{M}}$, and the Andrade’s parameter $\alpha$ of the latter. Because we lack a better model, we set $\tau_{\mathrm{A}}=\tau_{\mathrm{M}}$. In the first case, we take the value of the shear viscosity of the host star to be that estimated by Veras et al. (2019) for the Sun and Sun-like stars, which is $\eta_{1}\approx 10^{12}\,\mathrm{Pa}\,\mathrm{s}$. In the second, based on Fig. 1, we can expect that the considered planet has a composition similar to that of the Earth and, therefore, we can take the corresponding values of its rheological parameters. For the Earth, $\mu=0.8\times 10^{11}$ Pa, which is the most reliable value among the rheological parameters involved. Keeping in mind both Eq. (70) and the fact that the viscosity depends on the temperature, such that the former decreases when the latter increases, we then have to expect that the Maxwell time decreases as the planet’s internal temperature increases. At the same time, planets that are internally differentiated, such as the Earth, with partially molten cores, are expected to have the corresponding values of Andrade’s parameter, that is, $\alpha\in[0.15,0.2]$. Thus, due to this uncertainty, here we perform numerical simulations with different values of the less constrained parameters as follows: $\tau_{\mathrm{M}}=10,50,100$ years and $\alpha=0.15,0.2,0.3$.

Concerning the rest of the physical parameters, the values of the mean densities and surface gravity of both the planet and the host star are computed using the data given in Table 3. The maximum moments of inertia $C_{1}$ and $C_{2}$ of the star and the planet, respectively, are computed with:

[TABLE]

where $\xi_{k}$ is called the moment of inertia factor. For a homogeneous solid sphere, its value is $2/5$. But for internally differentiated rocky planets, its value is a little lower. In this sense, Zeng & Jacobsen (2017) showed that for both the Earth and Earth-like rocky planets, its value is approximately equal to $1/3$. In consequence, we can take $\xi_{1}=2/5$ and $\xi_{2}=1/3$. The last physical parameter to be set is the degree of triaxiality, indicating that we need a value for $(B-A)/C$. Such a parameter can be estimated following the work by Rodríguez et al. (2012), namely:

[TABLE]

For the Kepler-21 system, we obtain $7\times 10^{-4}$. However, in order to evaluate the effect of this parameter on the rotational evolution, we have also set $(B-A)/C_{2}=10^{-3}$ and zero.

With respect to the values of $J_{2}$, we have adopted three values, namely: $10^{-3}$, which is on the same order of magnitude of that of the Earth, along with $10^{-4}$, $10^{-5}$ and zero.

Lastly, regarding the orbital parameters we have set the initial values of $a$ and $e$ as those given by Table 3, correspondingly, as we are interested in the study of the fate of the Kepler-21 system and not in its past values. In this sense, we set $i_{1}(0)=0.5\degr$ and $i_{2}(0)=5\degr$ in order to estimate planarization and obliquity damping timescales.

Of course, the rotational state of the planet is completely unknown. However, it is believed that rocky planets rotate very fast after their formation (Miguel & Brunini, 2010). For exploratory purposes, we set $\theta_{2}(0)=0$ and $\dot{\theta}_{2}/n(0)=4.1$ for the planet. The reason behind the selection of this particular value for the planet’s initial spin rate is soon to become clear. As we shall see, it is very likely that the planet is already captured in synchronous rotation. The initial spin rate of the star is computed from its current estimated rotation period, which is given in Table 3.

5 Semianalytical and numerical characterization of the dynamical evolution

As we have already pointed out, the methodology developed in this work was applied to the study of the dynamical evolution of the Kepler-21 exoplanetary system, which is known to consist of a planet – presumably a rocky one – and its host star. Its orbital parameters and certain physical parameters are shown in Table 3.

We performed 13 numerical experiments varying the rheological and dynamical parameters of the planet because, as we have pointed out, these parameters are the less constrained ones. In these experiments the equations of motion – Eqs. (64a), (64b), (64c), (64d), (64e), (64f), (64g), (64h), (64i), (65a), and (65b) – were numerically integrated over a time span of $10^{5}$ years. The adopted values of the rheological and dynamical parameters in each integration (labeled with an ID number) are shown in Table 4. The aforementioned integrations were divided into three sets: in the first one, only the rheological parameters of the planet ($\alpha$ and $\tau_{M}$) were varied; in the second one, only the dynamical parameters of the latter were varied ($(B-A)/C_{2}$ and $J_{2}$); and in the third one, we only performed a numerical simulation setting $e=0.1$, which is within the uncertainty of the eccentricity value (see Table 3) in order to explore the possibility of entrapment in higher-than-synchronous spin-orbit resonances. As a general check, we verified that the total angular momentum of the systems is conserved in all the simulations we performed.

5.1 Results from the numerical simulations

We focus first on the rotational evolution of the system. In Fig. 2 we show the graphs of the normalized spin rate of the planet obtained from the first (Fig. 2(a)) and second (Fig. 2(b)) sets of integrations. It can be noted from Fig. 2(a) that the time necessary to reach the synchronous rotation is greatly affected by the variation of the rheological parameters, while the modification of the dynamical parameter have no consequences in this sense, as can seen in Fig. 2(b). However, we found that the synchronism is not always exact. In this sense, we are observing the influence of apsidal and nodal precession rates, as seen from the planet, on the corresponding tidal forcing frequencies, which were computed using the exact expression given by Eq. (16) in contraposition to the widely used approximation presented in Eq. (17).

The third set of numerical simulations consist only of an integration setting of $e=0.1$. It can be noted in Fig. 3 that the planet is trapped in the 3:2 spin-orbit resonance as can be expected from Fig. 4, where we plot the angular acceleration of the planet due to the secular terms of the tidal torque vs. its normalized spin rate (Makarov & Efroimsky, 2013). In that figure, the solid blue line (which correspond to $e=0.1$) crosses the horizontal reference axis near $\dot{\theta}_{2}/n=1.5$. This is the reason why we pointed out in Sect. 1 that eccentricity plays a key role in the possibility of capture in higher-than-synchronous resonances. This observation is in agreement with the results obtained by Noyelles et al. (2014).

In Fig. 4 we also may find a justification of this choice in the initial value of the planet’s spin rate (i.e., $\dot{\theta}_{2}=4.1\,n$). As can be noted in the aforementioned figure, the angular acceleration of the planet – due to the secular terms of the tidal torque – is practically constant for $\dot{\theta}_{2}/n\geq 4$ and thus, if we consider higher initial values of the planet’s spin rate, the only difference would be that the latter ends up reaching the 1:1 spin-orbit resonance some time later.

Now, we shift the focus to the orbital evolution of the Kepler-21 system. As a general observation, it can be noted that the major semiaxis and the eccentricity evolved in time in a similar way in each simulation when compared to the others, only differing in the final value of the corresponding time rates, as noted in Figs 5 and 6. This similarity can be understood by examining Fig. 7, where we can verify that the system behaves as expected from the same figure, provided that the planet’s rotational evolutionary timescale is much shorter than the orbital one. In this sense, we can know in advance how the system is expected to evolve over time before performing any numerical integration. In the case of Kepler-21, we see that the major semiaxis, the eccentricity, and the obliquity grow in time up to the 1:1 spin-orbit resonance and then start decreasing. This result is also consistent with Noyelles et al. (2014), who found that Mercury’s spin evolved so swiftly that the planet got trapped in its end-state only some 10-20 million of years after its accretion.

Perhaps it would be interesting to see the individual contributions of the tidal potentials of the star and of the planet onto the dynamics of the system that they are both a part of. In Figs. 8(a) and 8(b), we show the time derivatives of the major semiaxis and of the eccentricity, respectively, due to the secular terms of the star’s tidal potential, labeled with superscript (1), and to the planet’s tidal potential, labeled with a superscript (2). It can be observed that in both cases, the star’s contribution is approximately four orders of magnitude weaker than that of the planet’s.

Concerning the time evolution of the inclination ($i_{1}$), we obtained a seemingly counter-intuitive result, provided that its time derivative has a positive value of about $1.35\times 10^{-10}\,\degr\,\mbox{yr}^{-1}$, which is practically independent of the planet’s normalized spin rate. This implies that

[TABLE]

for the case of the Kepler-21 system, as can be noted from Eqs. (64d) and (63). As we can also note in Eq. (64d), the time evolution of $i_{1}$ depends only on the stellar parameters. Thus, the time derivative of the orbital inclination can be plotted against the normalized spin rate of the star, as shown in Fig. 9, where we consider two different values of the stellar viscosity. From Table 3 we can verify that $\dot{\theta}_{1}/n<1$ and, as a consequence, the time derivative of $i_{1}$ is inevitably positive by virtue of the dynamical model developed in this work.

With respect to the obliquity ($i_{2}$), it begins to decrease exponentially after the planet gets trapped in the synchronous rotation, as can be observed in Figs. 10(a) and 10(b). By virtue of Fig. 7, we can also know in advance the general behavior of the time evolution of the obliquity, at least qualitatively as in the case of the major semiaxis and the eccentricity discussed before.

It may seem that in Fig. 7 the time derivatives of the major semiaxis, the eccentricity, and the obliquity are equal to zero when the planet is trapped in the 1:1 spin-orbit resonance but, as seen in Figs. 5, 6, and 10, these orbital elements decrease in time after the planet has reached the synchronous rotation. There are two reasons behind this fact. On one hand, as pointed out before, the synchronism is not exact. On the other hand, if we take a look at a very narrow interval around the 1:1 spin orbit resonance in Fig. 7, as shown in Fig. 11, we note that the time derivatives of $a$, $e$ and $i_{2}$ are not exactly zero for $\dot{\theta}_{2}=n$.

As a final comment on the numerical simulations, it should be noted that this kind of problem involves different timescales of the dynamical evolution. As noted in Figs. 12 and 13(a), the time rates of the mean anomaly at the epoch, the argument of pericenter, and the longitude of the ascending node, both defined on the stellar equator, are on the order of $10^{-4}$; while the apsidal and nodal precessions as seen from the planet are on the order of $10^{1}$–$10^{2}$ (Fig. 13(b)). Then it is necessary for the numerical integration scheme to treat these timescales separately.

Taking into account that $n\approx 4.8\times 10^{4}{}\degr\mbox{yr}^{-1}$, as well, we then have that $\dot{M}\approx n$ and, in consequence, we consider a time step of $10^{-3}P_{\mathrm{orb}}$ ($\Delta t_{\mathrm{fast}}$) to integrate $\theta_{2}$, $M$, $\omega_{2}$, and $\Omega_{2}$, and a time step of one year ($\Delta t_{\mathrm{slow}}$) for the rest of the dynamical parameters, that is, $a$, $e$, $i_{1}$ and $i_{2}$.

The dynamical variables were integrated together but in different stages within the numerical integration scheme in order to implement the separation of the fast-varying and the slow-varying variables described before. More specifically, the former were integrated using fourth-order Runge-Kutta method ($\Delta t_{\mathrm{slow}}$), while the latter were integrated using the Bulirsch-Stoer method (Press et al., 1992) with the corresponding time step $\Delta t_{\mathrm{fast}}$.

The aforementioned integration scheme starts with an evaluation of the time derivatives of the slow-varying variables at the beginning of the integration interval, that is, the first stage of the Runge-Kutta method, for example, at $t_{i}$. Then the fast-varying variables are integrated using the Bulirsch-Stoer method up to $t_{i}+\Delta t_{\mathrm{slow}}/2.$ At that point, the second stage of the Runge-Kutta method is performed with the updated value of the planet’s spin rate in order to evaluate properly the tidal frequencies and the time derivatives of the slow-varying variables. Immediately, the fast-varying variables are integrated again up to the end of the integration interval, that is, $t_{i}+\Delta t_{\mathrm{slow}}$ where the new value of $\dot{\theta}_{2}$ is taken to compute the new values of both the tidal frequencies and the slow-varying variables’ time derivatives. Finally, the last stage of the Runge-Kutta method is carried out in order to compute the new values of the slow-varying variables; these are used to start the integration on the next interval together with the last computed values of the planet’s spin rate and the rest of the fast-varying variables.

5.2 Timescales of orbital decay and circularization

In general terms, we define the characteristic timescale of the evolution of an orbital element as that which is required for a certain change in that element to take place. Since the orbital inclination and the obliquity of the planet evolve in the way described earlier in this work, we only consider the timescales that correspond to the major semiaxis and the eccentricity. Mathematically expressed, the characteristic timescale of the former, $\tau_{a}$, and the latter, $\tau_{e}$, are defined as:

[TABLE]

[TABLE]

where $r_{\mathrm{Roche}}$ is the Roche limit and the values of $a$ and $e$ are the current ones. It can be noted that $\tau_{a}$ correspond to the time interval that elapses from the current value of $a$ up to the planet reaches the Roche radius, beyond which it would be completely destroyed by tidal forces. The latter is defined as (Murray & Dermott, 1999):

[TABLE]

Of course, this definition does not take into account the rigidity of the material that compose the body but only considers self-gravity. Even so, it can be considered as an upper limit. For the Kepler-21 system, $r_{\mathrm{Roche}}\approx 4.54\times 10^{-3}$ au.

Replacing in Eqs. (76) and (77) the corresponding values given in Table 3 and that of $\left\langle\frac{\mathrm{d}a}{\mathrm{d}t}\right\rangle$ and $\left\langle\frac{\mathrm{d}e}{\mathrm{d}t}\right\rangle$ from Eqs. (64a) and (64b) evaluated at $\dot{\theta}_{2}=n$ for different values of the rheological parameters, we obtain the values of orbital decay and circularization timescales. The results are shown in Table 5 together with the corresponding values of the rheological parameters. As can be observed in that table, since the orbital decay timescale is much larger than the others, when $a\approx r_{\mathrm{Roche}}$, the orbit will be completely circularized.

5.3 Possibility of detecting tidal effects

It may be an interesting exercise to consider the possibility of contrasting the results obtained from numerical simulation with those from observations. Within the limits of this work, we can answer the question of whether tidal effects can be detected only in part. In recent years, great advances have been made in observational techniques that allow for highly precise measurements of the physical and orbital parameters of planets orbiting stars other than the Sun. One of these quantities is the orbital period, $P_{\mathrm{orb}}$, of planets, which can be measured by observing several successive transits and is defined by:

[TABLE]

If we assume that the orbital parameters of a two body system evolve in time only due to the tides each body rises on the other – as was supposed in the last subsection – then the orbital period is also expected to evolve as a consequence of the time variation of the major semiaxis. This can be expressed mathematically as:

[TABLE]

From Eq. (1), we have:

[TABLE]

Then, Eq. (80) can be rewritten as:

[TABLE]

If we take the values of the parameters for the Kepler 21 system given in Table 3, compute $n$ with Eq. (1), and also take the value of $\left\langle\frac{\mathrm{d}a}{\mathrm{d}t}\right\rangle$, obtained from direct evaluation of Eq. (64a), assuming $\dot{\theta}_{2}=n$ and evaluating the tidal frequencies with Eq. (17), and then we replace them in Eq. (81), we get the following result:

[TABLE]

which is equivalent to a decrement of nearly $3.88\,\mbox{ms}$ per year. As a consequence, we would need more than 1000 years to measure a difference of about 4 seconds. The measured orbital period of planet Kepler-21b is $2.78578$ days with an uncertainty of $0.00003$ days (López-Morales et al., 2016), which is equivalent to approximately $2.6\,\mbox{s}$. This result clearly shows the impossibility of detecting tidal effects, if simply by measuring the orbital period of the planet, at least in this case.

In this regard, the question of whether it is possible to detect tidal effects by measuring the orbital period of certain exoplanets has motivated the publication of several recent works. Depending on each particular case, some authors conclude that it is not possible to detect tidal effect, as in Petrucci et al. (2019), and others not only assert that they measure a decrease of the orbital period, but also that it is a consequence of the orbital decay due to tidal interaction (Yee et al., 2019).

6 Final discussion and conclusions

In this work, we develop an accurate method for studying the dynamical evolution of a binary system formed by a rocky planet and its host star. We follow the work by Boué & Efroimsky (2019) to obtain a set of equations of motion which, on one hand, describe the secular time evolution of the orbital elements attributed to the secular terms of the tidal potentials of both partners and the secular terms related to the $J_{2}$ coefficient. On the other hand, it gives the time rate of the stellar spin rate due to the secular terms of the corresponding tidal torque and the time evolution of the proper rotation angle of the planet due to the combined action of the secular terms of the tidal torque and the periodic terms of the triaxiality-induced torque. Thus, the main results of this investigation are: Eqs. (64a), (64b), (64c), (64d), (64e), (64f), (64g), (64h), (64i), (65a), and (65b) which contain no singularities in the limit of vanishing eccentricities, inclination, and obliquity. Bearing in mind that the transformation introduced to obtain the equations only moves the latter singularity to $e=1$ and $i=90\degr$, we would like to stress the fact that it is a convenient substitution because it allows for the possibility that the eccentricity and the inclinations reach very low values (even zero) without affecting the stability of the numerical integration.

As a case study, we chose the Kepler-21 system. We find that the commonly accepted approximate expression of the tidal modes, given by Eq. (17), could no longer be valid in studies such as that carried out in this paper. In other words, at least in the case of the Kepler-21 system, the apsidal and nodal precessions induced by tides and the $J_{2}$-related terms is appreciable. Also, we find that the time of tidal synchronization strongly depends on the values of the rheological parameters, as is shown in Fig. 2(a).

Regarding the rotational state of planet Kepler-21b, even though it could be or have previously been trapped in the 3:2 spin-orbit resonance, we expect that it has already reached the synchronous rotation, not only because this is its most likely state, but as it orbits very close to a massive star with a radius nearly twice grater than that of our Sun, we can reasonably suppose that the planet is undergoing a strong internal dissipation which can cause an increase of its internal temperature with the consequent decrease of the viscosity and of the Maxwell time. As explained by Makarov et al. (2018), this decrease of the Maxwell time can cause the escape from higher-than-synchronous rotation states in a time interval shorter that the corresponding decrease of the corresponding capture probabilities due to the decrease of the eccentricity in time (Makarov, 2015).

It worth mentioning that from the point of view of tidal theory, the case of Kepler-21 – after the planet has reached the synchronous rotation – is analogous to the case of Phobos, which is gradually spiraling inwards towards Mars (Efroimsky & Lainey, 2007). In this sense, the timescales of orbital decay and circularization given in Table 5 are upper bounds, that is, the rates of orbital decay and circularization become greater as the planet approaches its host star. Even though the results shown in Table 5 suggest that the planet has a short time to exist, we can note that the more rapidly the eccentricity gets damped, the more the time it takes the planet to reach the Roche limit.

The time evolution of the orbital inclination ($i_{1}$) seems to be counter-intuitive. However, this is a particularity of the Kepler-21 system, but there could be other systems in which a more “expectable” conclusion is reached. As shown in Subsect. 5.1, tidal interaction can lead the major semiaxis, the eccentricity, and the inclinations to increase or decrease over time, depending on the particular combinations of the values of the rheological and dynamical parameters and variables involved (see Veras et al., 2019).

We would like to emphasize that even though we have addressed the study of a particular system, the formalism developed in this work is equally applicable to other analogous systems, formed by a star and a terrestrial planet. Moreover, this methodology could be adapted to model the dynamical evolution of binary asteroids. We will carry out the former and the latter types of studies in future works.

The analytical model that we propose in this work takes into account the tidal dissipation of both components of a binary system and the influence of both the $J_{2}$-related secular terms in the orbital evolution and the permanent-triaxiality-induced torque in the rotational evolution of one of the bodies, while the rotational evolution of the other is considered to be dominated only by the tidal interaction. We are aware that there are other effects that could be equally important to the dynamical evolution and could later be incorporated into the formalism, such as the oblateness of the star, General Relativity (Correia et al., 2011), dynamical tides in the star (Bolmont & Mathis, 2016; Hodžić et al., 2018), atmospheric tides in the planet (Correia & Laskar, 2003; Leconte et al., 2015), magnetic couplings between the star and the planet (Ahuir et al., 2019), and stellar evolution (Benbakoura et al., 2019).

Acknowledgements.

We acknowledge the funding from ANPCyT through PICT 1144-13 and CONICET through PIP 699-15. This investigation was supported by the institutional cluster HOPE at IAFE (http://www.iafe.uba.ar/HOPE/), where all the numerical simulations were performed, and by the SAO/NASA Astrophysics Data System. Santiago Luna wants to acknowledge the fruitful discussions with Michael Efroimsky, Gwenaël Boué, Dimitri Veras and Valeri Makarov, which helped to reach a better understanding of the dynamical evolution near spin-orbit resonances and to improve the quality of this work. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Expression of the sums of adimensional terms of the derivatives of the disturbing potentials

In the following, we give the complete expressions of the derivatives of the disturbing potentials considered in this work. Regarding the gradient of the tidal potentials of both bodies, we also give the secular terms of the corresponding derivatives. The latter is identified by the superscript ”(sec)”. Among the aforementioned expressions, the terms corresponding to $(lmpq)=(201q)$ and $(lmpq)=(220q)$ are shown, in particular. As earlier in this work, such terms are identified with the respective subscripts.

A.1 Major semiaxis

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.2 Eccentricity

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.3 Mean anomaly

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.4 Inclination with respect to the star’s and planet’s equators

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.5 Argument of pericenter with respect to the star’s and planet’s equators

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.6 Longitude of ascending node with respect to the star’s and planet’s equators

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Eccentricity functions

With regard to the eccentricity functions, they are identical to the well-known Hansen’s coefficients, $X_{\kappa}^{\nu,\mu}(e)$, which, in general, are defined as:

[TABLE]

where $\nu$, $\mu$ and $\kappa$ are integers, $r$ is the distance, $f$ is the true anomaly and $\textrm{i}=\sqrt{-1}$. For low to moderate eccentricities they can be computed as a power series of the eccentricity (Hughes 1981):

[TABLE]

where $\rho=j+\max(0,\kappa-\mu)$, $\sigma=j+\max(0,\mu-\kappa)$ and $N_{\rho,\sigma}^{\nu,\mu}$ are known as the Newcomb’s operators which can be computed recursively by (Proulx & McClain 1988):

[TABLE]

If $\sigma>\rho$ then $N^{\nu,\mu}_{\rho,\sigma}=N^{\nu,-\mu}_{\sigma,\rho}$. The recursion is initialized by considering that $N^{\nu,\mu}_{\rho,\sigma}=0$ if $\rho<0$ or $\sigma<0$ and by $N^{\nu,\mu}_{0,0}=1$.

The relationship between Kaula’s eccentricity functions, $G_{lpq}(e)$, and the Hansen coefficients is given by:

[TABLE]

Using Eq. (91) and Ec. (89), the $G_{lpq}(e)$ can be computed by means of:

[TABLE]

where $\rho$ and $\sigma$ are given now by: $\rho=j+\max(0,q)$ and $\sigma=j+\max(0,-q)$, respectively. Of course, the summation over $j$ must be also truncated at some maximum order in eccentricity, that is, $2j+|q|\leq j_{\mathrm{max}}$.

The calculation of the first derivatives of the eccentricity functions are straightforward by virtue of Eq. (91) is:

[TABLE]