Complete spin and orbital evolution of close-in bodies using a Maxwell viscoelastic rheology

Abstract

In this paper, we present a formalism designed to model tidal interaction with a viscoelastic body made of Maxwell material. Our approach remains regular for any spin rate and orientation, and for any orbital configuration including high eccentricities and close encounters. The method is to integrate simultaneously the rotation and the position of the planet as well as its deformation. We provide the equations of motion both in the body frame and in the inertial frame. With this study, we generalize preexisting models to the spatial case and to arbitrary multipole orders using a formalism taken from quantum theory. We also provide the vectorial expression of the secular tidal torque expanded in Fourier series. Applying this model to close-in exoplanets, we observe that if the relaxation time is longer than the revolution period, the phase space of the system is characterized by the presence of several spin-orbit resonances, even in the circular case. As the system evolves, the planet spin can visit different spin-orbit configurations. The obliquity is decreasing along most of these resonances, but we observe a case where the planet tilt is instead growing. These conclusions derived from the secular torque are successfully tested with numerical integrations of the instantaneous equations of motion on HD 80606 b. Our formalism is also well adapted to close-in super-Earths in multiplanet systems which are known to have non-zero mutual inclinations.

Appendices

Appendix 1: Spherical harmonic

By convention, Legendre associated polynomials are defined as

$$\begin{aligned} P_{l,m}(x) = \frac{1}{2^l l!} (1-x^2)^{m/2} \frac{d^{l+m}}{dx^{l+m}}(x^2-1)^l, \end{aligned}$$

(28)

with the symmetry

$$\begin{aligned} P_{l,-m} (x) = (-1)^m \frac{(l-m)!}{(l+m)!}P_{l,m}(x). \end{aligned}$$

(29)

The Schmidt semi-normalized spherical harmonics are defined as

$$\begin{aligned} Y_{l,m}(\theta ,\phi ) = (-1)^m \sqrt{\frac{(l-m)!}{(l+m)!}} P_{l,m}(\cos \theta ) \mathrm {e}^{\mathrm {i}m\phi } \end{aligned}$$

(30)

with the symmetry

$$\begin{aligned} Y_{l,-m}(\theta ,\phi ) = (-1)^m \bar{Y}_{l,m}(\theta ,\phi ) . \end{aligned}$$

(31)

Using the complex Cartesian coordinate system as defined in Varshalovich et al. (1988), for any unit vector $\hat{x}$, we have

$$\begin{aligned}&\displaystyle Y_{0,0}(\hat{\mathbf{x}}) = 1 , \end{aligned}$$

(32a)

$$\begin{aligned}&\displaystyle Y_{1,0}(\hat{\mathbf{x}}) = \hat{x}_0 , \end{aligned}$$

(32b)

$$\begin{aligned}&\displaystyle Y_{1,1}(\hat{\mathbf{x}}) = \hat{x}_+ , \end{aligned}$$

(32c)

$$\begin{aligned}&\displaystyle l\, Y_{l,0}(\hat{\mathbf{x}}) = (2l-1)\hat{x}_0 Y_{l-1,0}(\hat{\mathbf{x}}) - (l-1)Y_{l-2,0}(\hat{\mathbf{x}}), \end{aligned}$$

(32d)

$$\begin{aligned}&\displaystyle \sqrt{l+m} Y_{l,m}(\hat{\mathbf{x}}) = \sqrt{l-m}\, \hat{x}_0 Y_{l-1,m}(\hat{\mathbf{x}}) + \sqrt{2(l+m-1)} \hat{x}_+ Y_{l-1,m-1}. \end{aligned}$$

(32e)

The last two Eqs. (32d) and (32e) allow to recursively compute all spherical harmonics of order $m\ge 0$. Those with $m<0$ are deduced from the symmetry relation (31). Up to the degree 3 included, we have

$$\begin{aligned} \left\{ \begin{array}{ll} Y_{2,0} &{}= \displaystyle \frac{1}{2}(3 \hat{x}_0^2 -1) \\ Y_{2,1} &{}= \displaystyle \sqrt{3} \hat{x}_0 \hat{x}_+ \\ Y_{2,2} &{}= \displaystyle \frac{\sqrt{6}}{2} \hat{x}_+^2 \end{array} \right. , \qquad \left\{ \begin{array}{ll} Y_{3,0} &{}= \displaystyle \frac{5}{2}\hat{x}_0^3 - \frac{3}{2}\hat{x}_0 \\ Y_{3,1} &{}= \displaystyle \frac{\sqrt{6}}{4} \left( 5 \hat{x}_0^2 \hat{x}_+ - \hat{x}_+\right) \\ Y_{3,2} &{}= \displaystyle \frac{\sqrt{30}}{2} \hat{x}_0 \hat{x}_+^2 \\ Y_{3,3} &{}= \displaystyle \frac{\sqrt{10}}{2} \hat{x}_+^3 \end{array} \right. . \end{aligned}$$

(33)

Appendix 2: Ladder operators

Regular solid harmonics $x^l Y_{l,m}(\hat{\mathbf{x}})$ and irregular ones $Y_{l,m}(\hat{\mathbf{x}})/x^{l+1}$ are eigenvectors of each component of the gradient operator $\varvec{\nabla } = (\nabla _+, \nabla _0, \nabla _-)$ and of the angular momentum operator $\mathbf{J} = (J_+, J_0, J_-)$. The respective eigenvalues can be found in (e.g., Varshalovich et al. 1988). We have

$$\begin{aligned} \nabla _+ \left( x^{l}Y_{l,m}(\hat{\mathbf{x}})\right)&= -\,\sqrt{\frac{(l-m-1)(l-m)}{2}} \,x^{l-1}Y_{l-1,m+1}(\hat{\mathbf{x}}) \nonumber \\ \nabla _0 \left( x^{l}Y_{l,m}(\hat{\mathbf{x}})\right)&= +\,\sqrt{(l+m)(l-m)} \,x^{l-1} Y_{l-1,m}(\hat{\mathbf{x}}) \nonumber \\ \nabla _- \left( x^{l}Y_{l,m}(\hat{\mathbf{x}})\right)&= -\,\sqrt{\frac{(l+m-1)(l+m)}{2}} \,x^{l-1}Y_{l-1,m-1}(\hat{\mathbf{x}}), \end{aligned}$$

(34)

$$\begin{aligned} \nabla _+ \left( \frac{1}{x^{l+1}}Y_{l,m}(\hat{\mathbf{x}})\right)&= -\,\sqrt{\frac{(l+m+1)(l+m+2)}{2}} \frac{1}{x^{l+2}}Y_{l+1,m+1}(\hat{\mathbf{x}}) \nonumber \\ \nabla _0 \left( \frac{1}{x^{l+1}}Y_{l,m}(\hat{\mathbf{x}})\right)&= -\,\sqrt{(l+m+1)(l-m+1)} \frac{1}{x^{l+2}} Y_{l+1,m}(\hat{\mathbf{x}}) \nonumber \\ \nabla _- \left( \frac{1}{x^{l+1}}Y_{l,m}(\hat{\mathbf{x}})\right)&= -\,\sqrt{\frac{(l-m+1)(l-m+2)}{2}} \frac{1}{x^{l+2}}Y_{l+1,m-1}(\hat{\mathbf{x}}), \end{aligned}$$

(35)

and

$$\begin{aligned} J_+ \bigg ( f(x) Y_{l,m}(\hat{\mathbf{x}}) \bigg )&= -\,\sqrt{\frac{l(l+1)-m(m+1)}{2}} f(x) Y_{l,m+1}(\hat{\mathbf{x}}) \nonumber \\ J_0 \bigg ( f(x) Y_{l,m}(\hat{\mathbf{x}}) \bigg )&= m f(x) Y_{l,m}(\hat{\mathbf{x}}) \nonumber \\ J_- \bigg ( f(x) Y_{l,m}(\hat{\mathbf{x}}) \bigg )&= +\,\sqrt{\frac{l(l+1)-m(m-1)}{2}} f(x) Y_{l,m-1}(\hat{\mathbf{x}}), \end{aligned}$$

(36)

where f(x) is any function of the modulus $x=\Vert \vec {x}\Vert$.

Appendix 3: Rotation and Wigner matrices

Let a vector $\vec {x}$ and two coordinate systems $\mathcal{B}$ and $\mathcal{B}'$ such that $\mathbf{x}$ and $\mathbf{x}'$ are the coordinates of $\vec {x}$ in $\mathcal{B}$ and $\mathcal{B}'$, respectively. Let us further assume that $\mathbf{x}$ and $\mathbf{x}'$ are related to each other by a rotation of the form

$$\begin{aligned} \mathbf{x} = {{\mathsf {R}}}_3(\alpha ) {{\mathsf {R}}}_2(\beta ) {{\mathsf {R}}}_3(\gamma ) \mathbf{x}', \end{aligned}$$

where ${{\mathsf {R}}}_3$ and ${{\mathsf {R}}}_2$ are the matrices of rotation around the third and the second axes, respectively. Wigner D matrix ${\mathsf {D}}^l_{m,m'}(\alpha ,\beta ,\gamma )$ is defined such that (e.g., Varshalovich et al. 1988)

$$\begin{aligned} Y_{l,m}(\hat{\mathbf{x}}') = \sum _{m'=-l}^l {\mathsf {D}}^l_{m',m}(\alpha ,\beta ,\gamma ) Y_{l,m'}(\hat{\mathbf{x}}). \end{aligned}$$

(37)

Each element ${\mathsf {D}}^l_{m,m'}(\alpha ,\beta ,\gamma )$ can be written as (e.g., Varshalovich et al. 1988)

$$\begin{aligned} {\mathsf {D}}^l_{m,m'}(\alpha ,\beta ,\gamma ) = \mathrm {e}^{-\mathrm {i}m \alpha } {\mathsf {d}}^l_{m,m'}(\beta ) \mathrm {e}^{-\mathrm {i}m'\gamma }, \end{aligned}$$

(38)

where ${\mathsf {d}}^l_{m,m'}(\beta )$ is the Wigner d matrix. The inverse ${\mathsf {D}}^l_{m,m'}(-\gamma ,-\beta ,-\alpha )$ is given by the adjoint $\bar{{\mathsf {D}}}^l_{m',m}(\alpha ,\beta ,\gamma )$ of ${\mathsf {D}}^l_{m,m'}(\alpha ,\beta ,\gamma )$:

$$\begin{aligned} {\mathsf {D}}^l_{m,m'}(-\gamma ,-\beta ,-\alpha ) = \mathrm {e}^{\mathrm {i}m'\alpha } {\mathsf {d}}^l_{m',m}(\beta )\mathrm {e}^{\mathrm {i}m\gamma }. \end{aligned}$$

The convention 3-2-3 of the rotation (Eq. 37) is such that ${\mathsf {d}}^l_{m,m'}(\beta )$ is a real function. Wigner d matrix possesses many symmetries, among which (e.g., Varshalovich et al. 1988)

$$\begin{aligned} {\mathsf {d}}^l_{m,m'}(\beta ) = (-1)^{m-m'}{\mathsf {d}}^l_{-m,-m'}(\beta ) = (-1)^{m-m'}{\mathsf {d}}^l_{m',m}(\beta ) = {\mathsf {d}}^l_{-m',-m}(\beta ). \end{aligned}$$

Wigner d matrix can be constructed recursively using the hereinabove symmetries, the following initialization (e.g., Varshalovich et al. 1988)

$$\begin{aligned} {\mathsf {d}}^0_{0,0}(\beta )&= 1\,,\ {\mathsf {d}}^1_{0,0}(\beta ) = \cos \beta \,,\ {\mathsf {d}}^1_{1,-1}(\beta ) = \frac{1-\cos \beta }{2} \,,\ {\mathsf {d}}^1_{1,0}(\beta ) = -\frac{\sin \beta }{\sqrt{2}} \,,\ {\mathsf {d}}^1_{1,1}(\beta )\nonumber \\&= \frac{1+\cos \beta }{2} \end{aligned}$$

(39)

and the recurrence relation (Gimbutas and Greengard 2009)

$$\begin{aligned} \begin{aligned} {\mathsf {d}}^l_{m,m'}(\beta ) =&+\sqrt{\frac{(l+m')(l+m'-1)}{(l+m)(l+m-1)}} {\mathsf {d}}^1_{1,1}(\beta ) {\mathsf {d}}^{l-1}_{m-1,m'-1}(\beta ) \\&- \sqrt{\frac{(l+m')(l-m' )}{(l+m)(l+m-1)}} \sin (\beta ) {\mathsf {d}}^{l-1}_{m-1,m'}(\beta ) \\&+ \sqrt{\frac{(l-m')(l-m'-1)}{(l+m)(l+m-1)}} {\mathsf {d}}^1_{1,-1}(\beta ) {\mathsf {d}}^{l-1}_{m-1,m'+1}(\beta ) \end{aligned} \end{aligned}$$

(40)

which also implies

$$\begin{aligned} {\mathsf {d}}^l_{l,l}(\beta ) = {\mathsf {d}}^1_{1,1}(\beta ) {\mathsf {d}}^{l-1}_{l-1,l-1}(\beta ) \qquad \text {and} \qquad {\mathsf {d}}^l_{l,-l}(\beta ) = {\mathsf {d}}^1_{1,-1}(\beta ) {\mathsf {d}}^{l-1}_{l-1,1-l}(\beta ). \end{aligned}$$

(41)

The algorithm is the following:

For completeness, we also provide the explicit terms at order $l=2$ in Table 2.

Appendix 4: Time derivatives

Let a function $f(\vec {x}, t)$ developed in spherical harmonics as

$$\begin{aligned} f(\vec {x}, t) = \sum _{l,m} \bar{z}_{l,m}(t) Y_{l,m}(\hat{\mathbf{x}}) \end{aligned}$$

(42)

Table 2 Explicit Wigner d matrix ${\mathsf {d}}^2_{m,m'}(\beta )$

in the inertial frame $\mathcal{F}_0$, and as

$$\begin{aligned} f(\vec {x}, t) = \sum _{l,m} \bar{Z}_{l,m}(t) Y_{l,m}(\hat{\mathbf{X}}) \end{aligned}$$

(43)

in the body frame $\mathcal{F}_p$. For any constant vector $\vec {x}$ in $\mathcal{F}_p$, we have

$$\begin{aligned} \dot{\mathbf{x}} = {\varvec{\omega }} \times \mathbf{x} \qquad \mathrm {and} \qquad \dot{\mathbf{X}} = \mathbf{0} , \end{aligned}$$

(44)

with respect to the frame $\mathcal{F}_p$. By consequence, in $\mathcal{F}_p$, on the one hand,

$$\begin{aligned} \dot{f}(\vec {x}, t) = \sum _{l,m} \left( \dot{\bar{z}}_{l,m}(t) Y_{l,m}(\hat{\mathbf{x}}) + \bar{z}_{l,m}(t) \dot{\mathbf{x}}\cdot \varvec{\nabla } Y_{l,m}(\hat{\mathbf{x}}) \right) , \end{aligned}$$

(45)

and on the other hand,

$$\begin{aligned} \dot{f}(\vec {x}, t) = \sum _{l,m} \dot{\bar{Z}}_{l,m}(t) Y_{l,m}(\hat{\mathbf{X}}). \end{aligned}$$

(46)

But given that the time derivative of $\mathbf{x}$ is $\dot{\mathbf{x}} = \varvec{\omega }\times \mathbf{x}$, we get

$$\begin{aligned} \dot{\mathbf{x}}\cdot \varvec{\nabla } = ({\varvec{\omega }} \times \mathbf{x}) \cdot \varvec{\nabla } = \mathrm {i}({\varvec{\omega }} \cdot \mathbf{J}) \end{aligned}$$

(47)

where $\mathbf{J} = -\mathrm {i}\mathbf{x} \times \varvec{\nabla }$ is the angular momentum operator and where, by construction of the scalar product (Varshalovich et al. 1988),

$$\begin{aligned} {\varvec{\omega }} \cdot \mathbf{J} = - \omega _{+} J_{-} + \omega _0 J_0 - \omega _{-} J_{+}. \end{aligned}$$

(48)

We then define a matrix ${\mathsf {J}}({\varvec{\omega }})$ of size $(2l+1)\times (2l+1)$ such that

$$\begin{aligned} ({\varvec{\omega }} \cdot \mathbf{J}) Y_{l,m}(\hat{\mathbf{x}}) = \sum _{m'=-l}^l [{\mathsf {J}}({\varvec{\omega }})]^l_{m',m} Y_{l,m'}(\hat{\mathbf{x}}), \end{aligned}$$

(49)

where all non-zero coefficients are

$$\begin{aligned} \begin{array}{lll} \left[ {\mathsf {J}}({\varvec{\omega }})\right] ^l_{m-1,m} &{}=&{} - \sqrt{\frac{l(l+1)-m(m-1)}{2}} \omega _+, \\ \left[ {\mathsf {J}}({\varvec{\omega }})\right] ^l_{m,m} &{}=&{} m\, \omega _0 , \\ \left[ {\mathsf {J}}({\varvec{\omega }})\right] ^l_{m+1,m} &{}=&{} + \sqrt{\frac{l(l+1)-m(m+1)}{2}} \omega _-. \end{array} \end{aligned}$$

(50)

Combining Eqs. (10), (45–47), and (49), we obtain

$$\begin{aligned} \sum _{m'} {\mathsf {D}}^l_{m,m'}\dot{\bar{Z}}_{l,m'} = \dot{\bar{z}}_{l,m} + \mathrm {i}\sum _{m'}\,[{\mathsf {J}}({\varvec{\omega }})]^l_{m,m'} {\bar{z}}_{l,m'} . \end{aligned}$$

(51)

Appendix 5: Fourier transform

Let two functions $f(\vec {x}, t)$ and $g(\vec {x}, t)$ expanded in spherical harmonics as $f = \sum _l f_l$ and $g=\sum _l g_l$ with

$$\begin{aligned} f_l(\vec {x}, t) = \sum _{m=-l}^l \bar{Z}_{l,m}(t) Y_{l,m}(\hat{\mathbf{X}}) \qquad \text {and} \qquad g_l(\vec {x}, t) = \sum _{m=-l}^l \bar{Z}'_{l,m}(t) Y_{l,m}(\hat{\mathbf{X}}) \end{aligned}$$

in the frame $\mathcal{F}_p$ and

$$\begin{aligned} f_l(\vec {x}, t) = \sum _{m=-l}^l \bar{z}_{l,m}(t) Y_{l,m}(\hat{\mathbf{x}}) \qquad \text {and} \qquad g_l(\vec {x}, t) = \sum _{m=-l}^l \bar{z}'_{l,m}(t) Y_{l,m}(\hat{\mathbf{x}}) \end{aligned}$$

in $\mathcal{F}_0$. Let $\alpha$, $\beta$, and $\gamma =\omega t$ be the three angles such that

$$\begin{aligned} \mathbf{x} = {\mathsf {R}}_3(\alpha ) {\mathsf {R}}_2(\beta ) {\mathsf {R}}_3(\gamma ) \mathbf{X}. \end{aligned}$$

We have then

$$\begin{aligned} z_{l,m}(t) = \sum _{m'=-l}^l \bar{{\mathsf {D}}}^l_{m,m'}(t) Z_{l,m'}(t) \qquad \text {and} \qquad z'_{l,m}(t) = \sum _{m'=-l}^l \bar{{\mathsf {D}}}^l_{m,m'}(t) Z'_{l,m'}(t) \end{aligned}$$

(52)

with

$$\begin{aligned} {\mathsf {D}}^l_{m,m'} (t) = {\mathsf {D}}^l_{m,m'}(0) \mathrm {e}^{-\mathrm {i}m' \omega t}. \end{aligned}$$

Let us further assume that the two functions are related to each other in $\mathcal{F}_p$ by

$$\begin{aligned} f_l(\vec {x}, t) = h_l(t) * g_l(\vec {x}, t) \qquad \text {for all } l , \end{aligned}$$

where $h_l(t) \in \mathbb {R}$ is a real distribution. The symbol $*$ denotes the convolution product. As the convolution is done with respect to time, the orthogonality of the spherical harmonics implies that for all l and m,

$$\begin{aligned} Z_{l,m}(t) = h_l(t) * Z'_{l,m}(t). \end{aligned}$$

(53)

Combining Eqs. (52) and (53), we get

$$\begin{aligned} z_{l,m}(t)&= \sum _{m'=-l}^l \sum _{m''=-l}^l \int _{-\infty }^ \infty \bar{{\mathsf {D}}}^l_{m,m'}(t) h_l(t-t') {{\mathsf {D}}}^l_{m'',m'}(t') z'_{l,m''}(t')\,dt' \nonumber \\&= \sum _{m''=-l}^l {\mathsf {h}}^l_{m,m''}(t) * z'_{l,m''}(t), \end{aligned}$$

(54)

where

$$\begin{aligned} {\mathsf {h}}^l_{m,m''}(t) = \sum _{m'=-l}^l \bar{{\mathsf {D}}}^l_{m,m'}(0) h_l(t) \mathrm {e}^{\mathrm {i}m' \omega t} {\mathsf {D}}^l_{m'',m'}(0). \end{aligned}$$

In particular, if the rotation axis $\vec {\omega }$ is aligned with the third axis of $\mathcal{F}_0$ and $\mathcal{F}_p$, i.e., if $\alpha =\beta =0$, ${\mathsf {h}}^l_{m,m''}(t)$ is diagonal and we obtain

$$\begin{aligned} z_{l,m}(t) = \left( h_l(t)\mathrm {e}^{\mathrm {i}m \omega t}\right) * z'_{l,m}(t) \qquad \text {if } \quad \alpha =\beta =0. \end{aligned}$$

(55)

Taking the Fourier transform of Eqs. (54) and (55), we get

$$\begin{aligned} \underline{z}_{l,m}(\nu )= & {} \sum _{m''=-l}^l \underline{{\mathsf {h}}}^l_{m,m''}(\nu ) \underline{z}'_{l,m''}(\nu ) \quad \text {} \\ \end{aligned}$$

with

$$\begin{aligned} \underline{{\mathsf {h}}}^l_{m,m''}(\nu )= & {} \sum _{m'=-l}^l \bar{{\mathsf {D}}}^l_{m,m'}(0) \underline{h}_l(\nu -m'\omega ) {\mathsf {D}}^l_{m'',m'}(0), \end{aligned}$$

on the one hand, and

$$\begin{aligned} \underline{z}_{l,m}(\nu ) = \underline{h}_l(\nu - m\omega ) \underline{z}'_{l,m}(\nu ) \qquad \text {if}\quad \alpha =\beta =0, \end{aligned}$$

on the other.