GRIT: A Package for Structure-Preserving Simulations of Gravitationally Interacting Rigid Bodies

Under a specific fixed reference frame of Euclidean space ${{\mathbb{R}}}^{3}$ with basis $\left({{\boldsymbol{E}}}_{1},{{\boldsymbol{E}}}_{2},{{\boldsymbol{E}}}_{3}\right)$, the spatial position and the translational speed of a body can be expressed by vectors in ${{\mathbb{R}}}^{3}$. Meanwhile, the body frame (Figure 1) attached to the body gives fixed coordinates of each small particle of the body. As can be seen from Figure 1, orthogonal bases $\left({{\boldsymbol{e}}}_{{\bf{1}}},{{\boldsymbol{e}}}_{{\bf{2}}},{{\boldsymbol{e}}}_{{\bf{3}}}\right)$ form a body frame of this rigid body. As the body moves along the dashed trajectory following the arrows, as well as self rotating from time t₀ to time t₁, the coordinates of points P, Q under the body frame stay the same, without being subject to the motion of the rigid body.

Zoom In Zoom Out Reset image size

The configuration of a rigid body is described by both the position of its center of mass and its rotational orientation. The orientation in the reference frame can be expressed as a rotation by an orthogonal matrix R (t) ∈ SO(3) from the body frame (e.g., z-axis of the body frame at time t will be ${\boldsymbol{R}}(t)\cdot {\left[0\ \ 0\ \ 1\right]}^{T}$ in the reference frame). To switch between the inertia frame and the body frame, one can simply left multiply the rotation matrix R or R ⁻¹. Note that R ∈ SO(3) and if a numerical method can keep R exactly in this Lie group, then its inverse will be equal to its transpose; that is, R ⁻¹ = R ^T .

By denoting ${\boldsymbol{\Omega }}={\left[{{\rm{\Omega }}}_{1}\ \ {{\rm{\Omega }}}_{2}\ \ {{\rm{\Omega }}}_{3}\right]}^{T}\in {{\mathbb{R}}}^{3}$ the angular velocity of the rigid body under the body frame, the direction of Ω matches the rotational axis and ${\parallel {\boldsymbol{\Omega }}\parallel }_{2}$ represents the rotational speed. Consider a mass point ${\boldsymbol{x}}={\left[{x}_{1}\ \ {x}_{2}\ \ {x}_{3}\right]}^{T}$ in one rigid body under the body frame, its speed under the body frame can be expressed as

$$\begin{eqnarray}{\boldsymbol{\Omega }}\times {\boldsymbol{x}}=\left[\begin{array}{ccc}0 & -{{\rm{\Omega }}}_{3} & {{\rm{\Omega }}}_{2}\\ {{\rm{\Omega }}}_{3} & 0 & -{{\rm{\Omega }}}_{1}\\ -{{\rm{\Omega }}}_{2} & {{\rm{\Omega }}}_{1} & 0\end{array}\right]\left[\begin{array}{c}{x}_{1}\\ {x}_{2}\\ {x}_{3}\end{array}\right]=\hat{{\boldsymbol{\Omega }}}{\boldsymbol{x}},\end{eqnarray} \tag{ 1 }$$

where the hat-map $\hat{\cdot }$ is an isomorphism from the Lie algebra ${\mathfrak{so}}(3)$ to 3-by-3 skew-symmetric matrices, defined by

$$\begin{eqnarray}\hat{{\boldsymbol{\Omega }}}:= \left[\begin{array}{ccc}0 & -{{\rm{\Omega }}}_{3} & {{\rm{\Omega }}}_{2}\\ {{\rm{\Omega }}}_{3} & 0 & -{{\rm{\Omega }}}_{1}\\ -{{\rm{\Omega }}}_{2} & {{\rm{\Omega }}}_{1} & 0\end{array}\right].\end{eqnarray} \tag{ 2 }$$

In addition, the inverse map of $\hat{\cdot }$ is denoted by $\check{\cdot }$.

With the angular velocity, the rotational kinetic energy of this rigid body can be expressed as

$$\begin{eqnarray}&&{{\boldsymbol{T}}}^{{\boldsymbol{rot}}}\left({\boldsymbol{\Omega }}\right)=\displaystyle \frac{1}{2}{{\boldsymbol{\Omega }}}^{T}{\boldsymbol{J}}{\boldsymbol{\Omega }}\end{eqnarray} \tag{ 3 }$$

with J the (standard) moment of inertia tensor. Specifically, for an ellipsoid with semiaxes a, b, c and mass M, and by choosing the principal axes as the body frame such that x, y, z-axes matches semiaxes and taking the integral, we have the following moment of inertia tensor for a uniform density object.

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{J}} & = & {\displaystyle \int }_{{ \mathcal B }}\rho ({\boldsymbol{x}}){\hat{{\boldsymbol{x}}}}^{T}\hat{{\boldsymbol{x}}}\,{d}^{3}{\boldsymbol{x}}\\ & = & \left[\begin{array}{ccc}\displaystyle \frac{1}{5}M({b}^{2}+{c}^{2}) & 0 & 0\\ 0 & \displaystyle \frac{1}{5}M({a}^{2}+{c}^{2}) & 0\\ 0 & 0 & \displaystyle \frac{1}{5}M({a}^{2}+{b}^{2})\end{array}\right].\end{array}\end{eqnarray} \tag{ 4 }$$

Note that one may substitute this with the principal moment of inertia directly.

Alternatively, the rotational kinetic energy can also be expressed as

$$\begin{eqnarray}&&{{\boldsymbol{T}}}^{{\boldsymbol{rot}}}\left(\hat{{\boldsymbol{\Omega }}}\right)=\displaystyle \frac{1}{2}{Tr}\left[\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}{\hat{{\boldsymbol{\Omega }}}}^{T}\right].\end{eqnarray} \tag{ 5 }$$

with ${{\boldsymbol{J}}}^{(d)}={\int }_{{ \mathcal B }}\rho ({\boldsymbol{x}}){\boldsymbol{x}}{{\boldsymbol{x}}}^{T}\,d{\boldsymbol{x}}$ (nonstandard) moment of inertia. We also have ${{\boldsymbol{J}}}^{(d)}=\tfrac{1}{2}\,{Tr}\left[{\boldsymbol{J}}\right]{{\boldsymbol{I}}}_{3\times 3}-{\boldsymbol{J}}$ (${\boldsymbol{J}}={Tr}\left[{{\boldsymbol{J}}}^{(d)}\right]{{\boldsymbol{I}}}_{3\times 3}\,-{{\boldsymbol{J}}}^{(d)}$ ).

By definition, the angular momentum in the body frame is Π = J Ω. By left multiplying the rotation matrix R , the angular velocity and the angular momentum in the inertia frame are ω = R Ω and π = R Π respectively.

Express ${\boldsymbol{R}}(t)=\left[{{\boldsymbol{c}}}_{1}(t)| {{\boldsymbol{c}}}_{2}(t)| {{\boldsymbol{c}}}_{3}(t)\right]$ where c ₁(t), c ₂(t), c ₃(t) are columns of R (t). We have c ₁(t), c ₂(t), c ₃(t) representing directions of three axes of the body in the reference frame, respectively. By the definition of angular velocity, we have ${\dot{{\boldsymbol{c}}}}_{i}(t)=\hat{{\boldsymbol{\Omega }}}{c}_{i}(t)$ for i = 1,2,3, thus

$$\begin{eqnarray}&&\dot{{\boldsymbol{R}}}(t)=\hat{{\boldsymbol{\omega }}}{\boldsymbol{R}}(t).\end{eqnarray} \tag{ 6 }$$

By multiplying both sides of Equation (6) with R (t)^T , we have $\dot{{\boldsymbol{R}}}{(t){\boldsymbol{R}}(t)}^{T}=\hat{{\boldsymbol{\omega }}}$, which is a skew-symmetric matrix. Considering the speed of an arbitrary mass point x , ${{\boldsymbol{v}}}_{{\boldsymbol{x}}}={\boldsymbol{R}}\hat{{\boldsymbol{\Omega }}}{\boldsymbol{x}}=\hat{{\boldsymbol{\omega }}}\left[{\boldsymbol{R}}{\boldsymbol{x}}\right]$. Thus ${\boldsymbol{R}}\hat{{\boldsymbol{\Omega }}}=\hat{{\boldsymbol{\omega }}}{\boldsymbol{R}}$, which implies $\hat{{\boldsymbol{\Omega }}}={{\boldsymbol{R}}}^{T}\hat{{\boldsymbol{\omega }}}{\boldsymbol{R}}={{\boldsymbol{R}}}^{T}\dot{{\boldsymbol{R}}}$.

To summarize, the angular velocity and angular momentum of a rigid body in different frames are denoted as follows,

with Π = J Ω and π = J ω . Specifically, we have $\hat{{\boldsymbol{\Omega }}}={{\boldsymbol{R}}}^{T}\dot{{\boldsymbol{R}}}$ and $\hat{{\boldsymbol{\Omega }}}=\dot{{\boldsymbol{R}}}{{\boldsymbol{R}}}^{T}$.

The rotation matrix R and the angular momentum Π will be utilized to describe a rigid body when we design an N-rigid-body integrator later (further details can be found in Section 3).

The mass of the ith body is denoted by m_i ; ${{\boldsymbol{q}}}_{i}\in {{\mathbb{R}}}^{3}$ denotes the position of the ith body; ${{\boldsymbol{p}}}_{i}\in {{\mathbb{R}}}^{3}$ denotes the linear momentum of the ith body; R _i ∈ SO(3) denotes the rotation matrix of the ith body; ${{\boldsymbol{\Pi }}}_{i}\in {{\mathbb{R}}}^{3}$ denotes the angular momentum of the ith body; and ${{\boldsymbol{J}}}_{i}\in {{\mathbb{R}}}^{3\times 3}$ denotes the (standard) moment of inertia tensor for the ith body.

The Hamiltonian of this system consists of the linear kinetic energy ${T}^{\mathrm{linear}}={\sum }_{i}\tfrac{1}{2}{{\boldsymbol{p}}}_{i}^{T}{{\boldsymbol{p}}}_{i}/{m}_{i}$, the rotational kinetic energy ${T}^{\mathrm{rot}}={\sum }_{i}\tfrac{1}{2}{{\boldsymbol{\Pi }}}_{i}^{T}{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}$ and the potential energy

$$\begin{eqnarray}&&V({\mathsf{q}},{\mathsf{R}})=\sum _{i\lt j}{V}_{{ij}}\left({{\boldsymbol{q}}}_{i},{{\boldsymbol{q}}}_{j},{{\boldsymbol{R}}}_{{\boldsymbol{i}}},{{\boldsymbol{R}}}_{{\boldsymbol{j}}}\right).\end{eqnarray} \tag{ 7 }$$

Denote ${\mathsf{q}}\,=\,\left\{{{\boldsymbol{q}}}_{1},\,{{\boldsymbol{q}}}_{2},\ldots ,{{\boldsymbol{q}}}_{N}\right\}$, ${\mathsf{p}}\,=\,\left\{{{\boldsymbol{p}}}_{1},\,{{\boldsymbol{p}}}_{2},\ldots ,{{\boldsymbol{p}}}_{N}\right\}$, ${\rm{\Pi }}\,=\left\{{{\boldsymbol{\Pi }}}_{1},{{\boldsymbol{\Pi }}}_{2},\ldots ,{{\boldsymbol{\Pi }}}_{N}\right\}$, ${\mathsf{R}}=\left\{{{\boldsymbol{R}}}_{1},{{\boldsymbol{R}}}_{2},\ldots ,{{\boldsymbol{R}}}_{N}\right\}$. The Hamiltonian can be expressed as

$$\begin{eqnarray}&&\begin{array}{l}H({\mathsf{q}},{\mathsf{p}},{\rm{\Pi }},{\mathsf{R}})=\sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{p}}}_{i}^{T}{{\boldsymbol{p}}}_{i}/{m}_{i}\\ \ \ +\ \sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{\Pi }}}_{i}^{T}{{\boldsymbol{I}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}+V({\mathsf{q}},{\mathsf{R}})\ \mathrm{with}\ {{\boldsymbol{R}}}_{{\boldsymbol{i}}}\in {\mathsf{SO}}(3).\end{array}\end{eqnarray} \tag{ 8 }$$

The true potential energy between ith body and jth body is

$$\begin{eqnarray}&&{\int }_{{{ \mathcal B }}_{i}}{\int }_{{{ \mathcal B }}_{j}}-\displaystyle \frac{{ \mathcal G }\rho ({{\boldsymbol{x}}}_{{\boldsymbol{i}}})\rho ({{\boldsymbol{x}}}_{{\boldsymbol{j}}})}{\parallel \left({{\boldsymbol{q}}}_{i}+{{\boldsymbol{R}}}_{{\boldsymbol{i}}}{{\boldsymbol{x}}}_{{\boldsymbol{i}}}\right)-\left({{\boldsymbol{q}}}_{j}+{{\boldsymbol{R}}}_{{\boldsymbol{j}}}{{\boldsymbol{x}}}_{{\boldsymbol{j}}}\right)\parallel }\,d{{\boldsymbol{x}}}_{{\boldsymbol{j}}}\,d{{\boldsymbol{x}}}_{{\boldsymbol{i}}}.\end{eqnarray} \tag{ 9 }$$

We may approximate this as V_ij (in Equation (7)) by Taylor expanding the denominator. By expanding to the 2nd order with respect to the radius of the planet over the distance between two bodies (see Appendix A), the approximated potential is,

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{{{ \mathcal B }}_{i}}{\displaystyle \int }_{{{ \mathcal B }}_{j}}-\displaystyle \frac{{ \mathcal G }\rho ({{\boldsymbol{x}}}_{i})\rho ({{\boldsymbol{x}}}_{j})}{\parallel \left({{\boldsymbol{q}}}_{i}+{{\boldsymbol{R}}}_{{\boldsymbol{i}}}{{\boldsymbol{x}}}_{{\boldsymbol{i}}}\right)-\left({{\boldsymbol{q}}}_{j}+{{\boldsymbol{R}}}_{{\boldsymbol{j}}}{{\boldsymbol{x}}}_{{\boldsymbol{j}}}\right)\parallel }\,d{{\boldsymbol{x}}}_{{\boldsymbol{j}}}\,d{{\boldsymbol{x}}}_{{\boldsymbol{i}}}\\ \approx \ -\displaystyle \frac{{ \mathcal G }{m}_{i}{m}_{j}}{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\parallel }-\displaystyle \frac{{ \mathcal G }\left({m}_{i}{Tr}[{{\boldsymbol{J}}}_{i}]+{m}_{j}{Tr}[{{\boldsymbol{J}}}_{{\boldsymbol{j}}}]\right)}{2{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\parallel }^{3}}\\ \ +\ \displaystyle \frac{3{ \mathcal G }{\left({{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\right)}^{T}\left({m}_{j}{{\boldsymbol{R}}}_{i}{{\boldsymbol{J}}}_{i}{{\boldsymbol{R}}}_{i}^{T}+{m}_{i}{{\boldsymbol{R}}}_{j}{{\boldsymbol{J}}}_{{\boldsymbol{j}}}{{\boldsymbol{R}}}_{j}^{T}\right)\left({{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\right)}{2{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\parallel }^{5}}\end{array}\end{eqnarray} \tag{ 10 }$$

with $-\tfrac{{ \mathcal G }{m}_{i}{m}_{j}}{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\parallel }$ being the potential of purely point-mass interactions and the rest being the corrections of the potential due to the body i and j not being point masses. If we further expand the potential to the 4th order (see Appendix A), then rigid-body rigid-body interactions will also be included as higher order corrections. For example, fourth order potential has recently been considered for binary asteroids with large non-spherical terms, and leads to interesting effects (Hou et al. 2017).

The Lagrangian for a system consisting of one rigid body is a function of R (t) and $\dot{{\boldsymbol{R}}}(t)$ by plugging in $\hat{{\boldsymbol{\Omega }}}={{\boldsymbol{R}}}^{T}\dot{{\boldsymbol{R}}}$ in Equation (5),

$$\begin{eqnarray}&&L\left({\boldsymbol{R}},\dot{{\boldsymbol{R}}}\right)=\displaystyle \frac{1}{2}{Tr}\left[\dot{{\boldsymbol{R}}}{{\boldsymbol{J}}}_{d}{\dot{{\boldsymbol{R}}}}^{T}\right]-V({\boldsymbol{R}}).\end{eqnarray} \tag{ 11 }$$

By utilizing the constraint R ^T R − I = 0 (Appendix B.1) or using the variational principle of Hamilton's for Lie group (Appendix B.2), one can derive the equations of motion

$$\begin{eqnarray}&&\left\{\begin{array}{l}\dot{{\boldsymbol{R}}}={\boldsymbol{R}}\widehat{{{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}},\\ \dot{{\boldsymbol{\Pi }}}={\boldsymbol{\Pi }}\times {{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

Similarly, the equations of motion of the N-rigid-body system for the Hamiltonian (Equation (8)) are,

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{{\boldsymbol{q}}}}_{i}=\displaystyle \frac{{{\boldsymbol{p}}}_{i}}{{m}_{i}},\\ {\dot{{\boldsymbol{p}}}}_{i}=-\displaystyle \frac{\partial V}{\partial {{\boldsymbol{q}}}_{i}},\\ {\dot{{\boldsymbol{R}}}}_{i}={{\boldsymbol{R}}}_{{\boldsymbol{i}}}\widehat{{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}},\\ {\dot{{\boldsymbol{\Pi }}}}_{i}={{\boldsymbol{\Pi }}}_{i}\times {{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}-{\left({{\boldsymbol{R}}}_{{\boldsymbol{i}}}^{T}\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}-{\left(\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}\right)}^{T}{{\boldsymbol{R}}}_{{\boldsymbol{i}}}\right)}^{\vee }.\end{array}\right.\end{eqnarray} \tag{ 13 }$$

In this section, we utilize the splitting method to construct symplectic integrators. A diverse range of symplectic integrators with different accuracies and time complexities can be designed because the splitting method is quite flexible in terms of splitting and composition. Based on our Hamiltonian of the N-rigid-body system, we will explore three different types of integrators. One integrator splits the Hamiltonian into two parts with comparable size, the other two integrators split the Hamiltonian into, respectively, three and four parts corresponding to various magnitudes and hence different timescales.

In terms of the shape of rigid bodies, we make the axis-symmetric assumption in this section for simplicity; that is, without loss of generality, ${{\boldsymbol{J}}}_{i}=\left[\begin{array}{ccc}{J}_{i}^{(1)} & 0 & 0\\ 0 & {J}_{i}^{(1)} & 0\\ 0 & 0 & {J}_{i}^{(3)}\end{array}\right]$. Different splitting mechanisms can be applied for general rigid bodies that are not axis-symmetric (see Section 3.5).

3.3.1. Classical Splitting for Rigid Body: H = H₁ + H₂ with $\tfrac{{H}_{1}}{{H}_{2}}={ \mathcal O }(1)$

One way of splitting is H = H₁ + H₂ following Dullweber et al. (1997), with

$$\begin{eqnarray}&&\left\{\begin{array}{l}{H}_{1}({\mathsf{q}},{\mathsf{p}},{\rm{\Pi }},{\mathsf{R}})=\displaystyle \sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{p}}}_{i}^{T}{{\boldsymbol{p}}}_{i}/{m}_{i}+\displaystyle \sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{\Pi }}}_{i}^{T}{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i},\\ {H}_{2}({\mathsf{q}},{\mathsf{p}},{\rm{\Pi }},{\mathsf{R}})=V({\mathsf{q}},{\mathsf{R}})=\displaystyle \sum _{i\ne j}{V}_{{ij}}\left({{\boldsymbol{q}}}_{i},{{\boldsymbol{q}}}_{j},{{\boldsymbol{R}}}_{{\boldsymbol{i}}},{{\boldsymbol{R}}}_{{\boldsymbol{j}}}\right).\end{array}\right.\end{eqnarray} \tag{ 14 }$$

For H₁, the equations of motion are

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{{\boldsymbol{q}}}}_{i}=\displaystyle \frac{{{\boldsymbol{p}}}_{i}}{{m}_{i}},\\ {\dot{{\boldsymbol{p}}}}_{i}=0,\\ {\dot{{\boldsymbol{R}}}}_{i}={{\boldsymbol{R}}}_{i}\widehat{{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}},\\ {\dot{{\boldsymbol{\Pi }}}}_{i}={{\boldsymbol{\Pi }}}_{i}\times {{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

In Equation (15), the 4th equation is the Euler equation for a free rigid body. This is exactly solvable, and the solution expression is particularly simple for axial-symmetric bodies:

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{i}(t)=\exp \left(-\theta t\widehat{\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]}\right){{\boldsymbol{\Pi }}}_{i}(0):= {R}_{z}^{T}(\theta t){{\boldsymbol{\Pi }}}_{i}(0)\end{eqnarray} \tag{ 16 }$$

with $\theta =\left(\tfrac{1}{{J}_{i}^{(3)}}-\tfrac{1}{{J}_{i}^{(1)}}\right){{\boldsymbol{\Pi }}}_{i}^{T}(0)\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$ and R_z being the rotation matrix. Take Π_i (t) back to the 3rd equation of Equation (15), we can also obtain R _i (t).

Therefore, the flow ${\phi }_{t}^{[1]}$ of H₁ is,

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{q}}}_{i}(t)={{\boldsymbol{q}}}_{i}(0)+\displaystyle \frac{{{\boldsymbol{p}}}_{i}}{{m}_{i}}t,\\ {{\boldsymbol{p}}}_{i}(t)={{\boldsymbol{p}}}_{i}(0),\\ {{\boldsymbol{R}}}_{{\boldsymbol{i}}}(t)={{\boldsymbol{R}}}_{i}(0)\,{{\rm{R}}}_{{\boldsymbol{\Pi }}(0)}\left(\displaystyle \frac{\parallel {\boldsymbol{\Pi }}(0)\parallel }{{J}_{i}^{(1)}}t\right)\,{{\rm{R}}}_{z}\left(\theta t\right),\\ {{\boldsymbol{\Pi }}}_{i}(t)={{\rm{R}}}_{z}^{T}\left(\theta t\right)\,{\boldsymbol{\Pi }}(0),\end{array}\right.\end{eqnarray} \tag{ 17 }$$

with R_z and R_Π(0) being rotation matrices representing the rotations around the z-axis and Π(0), respectively.

For H₂, the equations of motion are

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{{\boldsymbol{q}}}}_{i}=0,\\ {\dot{{\boldsymbol{p}}}}_{i}=-\displaystyle \frac{\partial V}{\partial {{\boldsymbol{q}}}_{i}},\\ {\dot{{\boldsymbol{R}}}}_{i}=0,\\ {\dot{{\boldsymbol{\Pi }}}}_{i}=-{\left({{\boldsymbol{R}}}_{{\boldsymbol{i}}}^{T}\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}-{\left(\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}\right)}^{T}{{\boldsymbol{R}}}_{{\boldsymbol{i}}}\right)}^{\vee }.\end{array}\right.\end{eqnarray} \tag{ 18 }$$

Because q _i and R _i stay constant, we have p _i and Π_i changing at constant rates. Therefore, the flow ${\phi }_{t}^{[2]}$ for H₂ is given by

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{q}}}_{i}(t)={{\boldsymbol{q}}}_{i}(0),\\ {{\boldsymbol{p}}}_{i}(t)={{\boldsymbol{p}}}_{i}(0)-\displaystyle \frac{\partial V}{\partial {{\boldsymbol{q}}}_{i}}t,\\ {{\boldsymbol{R}}}_{{\boldsymbol{i}}}(t)={{\boldsymbol{R}}}_{{\boldsymbol{i}}}(0),\\ {{\boldsymbol{\Pi }}}_{i}(t)={{\boldsymbol{\Pi }}}_{i}(0)-{\left({{\boldsymbol{R}}}_{{\boldsymbol{i}}}^{T}\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}-{\left(\displaystyle \frac{\partial V}{\partial {{\boldsymbol{R}}}_{{\boldsymbol{i}}}}\right)}^{T}{{\boldsymbol{R}}}_{{\boldsymbol{i}}}\right)}^{\vee }t.\end{array}\right.\end{eqnarray} \tag{ 19 }$$

We may compose ${\phi }_{t}^{[1]}$ and ${\phi }_{t}^{[2]}$ via ${ \mathcal C }$ to construct different symplectic integrators (see, e.g., McLachlan & Quispel 2002), see Figure 2. To name a few, set ${ \mathcal C }$ as ${{ \mathcal C }}_{\mathrm{Euler}}$, ${\phi }_{h}={\phi }_{h}^{[1]}\circ {\phi }_{h}^{[2]}$ is a 1st order scheme with h being the step size (see Appendix D for ${{ \mathcal C }}_{\mathrm{Euler}}$ and the following composition methods ${{ \mathcal C }}_{\mathrm{Verlet}}$ and ${{ \mathcal C }}_{S6}$).

Zoom In Zoom Out Reset image size

By applying symmetric composition ${{ \mathcal C }}_{\mathrm{Verlet}}$, a 2nd order integrator ${{ \mathcal T }}_{2}$ is in the form of

$$\begin{eqnarray}&&{\phi }_{h}^{\mathrm{Verlet}}:= {{ \mathcal C }}_{\mathrm{Verlet}}\left({\phi }_{h}^{[1]},{\phi }_{h}^{[2]}\right)={\phi }_{\displaystyle \frac{h}{2}}^{[1]}\circ {\phi }_{h}^{[2]}\circ {\phi }_{\displaystyle \frac{h}{2}}^{[1]}\end{eqnarray} \tag{ 20 }$$

By applying ${\phi }_{h}^{\mathrm{TriJump}}$ (Suzuki 1990), we have the following 4th-order scheme ${{ \mathcal T }}_{4}$,

$$\begin{eqnarray}&&{\phi }_{{\gamma }_{1}h}^{\mathrm{Verlet}}\circ {\phi }_{{\gamma }_{2}h}^{\mathrm{Verlet}}\circ {\phi }_{{\gamma }_{3}h}^{\mathrm{Verlet}},\end{eqnarray} \tag{ 21 }$$

with ${\gamma }_{1}={\gamma }_{3}=\tfrac{1}{2-{2}^{1/3}}$, γ₂ = 1 − 2γ₁. Similarly, a 6th-order scheme ${{ \mathcal T }}_{6}$ can be constructed by composing ${\phi }_{h}^{[1]}$, ${\phi }_{h}^{[2]}$ with ${{ \mathcal C }}_{S6}$

In the package, ${{ \mathcal T }}_{2}$, ${{ \mathcal T }}_{4}$ and ${{ \mathcal T }}_{6}$ are implemented.

3.3.2. Tailored Splitting I: H = H₁ + H₂ + H₃ + H₄ with $\tfrac{{H}_{1}}{{H}_{2}}={ \mathcal O }(1)$ and $\tfrac{{H}_{3}}{{H}_{1}},\tfrac{{H}_{4}}{{H}_{2}}={ \mathcal O }(\varepsilon )$

In contrast from point-mass systems, which can already exhibit dynamics over multiple timescales, the N-rigid-body system can have additional timescales created by the rotational dynamics.

Thus, we further split the Hamiltonian into more terms of different magnitudes, which produce flows at different timescales, and we then carefully compose them. ⁴ More specifically, consider H = H₁ + H₂ + H₃ + H₄ with

$$\begin{eqnarray}&&\left\{\begin{array}{l}{H}_{1}({\mathsf{q}},{\mathsf{p}})=\displaystyle \sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{p}}}_{i}^{T}{{\boldsymbol{p}}}_{i}/{m}_{i}\\ {H}_{2}({\mathsf{q}},{\mathsf{p}})=-\displaystyle \sum _{i\lt j}\displaystyle \frac{{ \mathcal G }{m}_{i}{m}_{j}}{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{j}\parallel },\\ {H}_{3}({\rm{\Pi }})=\displaystyle \sum _{i}\displaystyle \frac{1}{2}{{\boldsymbol{\Pi }}}_{i}^{T}{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i},\\ {H}_{4}({\mathsf{q}},{\mathsf{R}})=V({\mathsf{q}},{\mathsf{R}})-{H}_{2}.\end{array}\right.\end{eqnarray} \tag{ 22 }$$

Here, H₁, H₂ have comparable size and $\tfrac{{H}_{3}}{{H}_{1}},\tfrac{{H}_{4}}{{H}_{2}}={ \mathcal O }\left(\varepsilon \right)$, with ε being a small scaling parameter determined by the properties of the system. Based on scales of the dynamics, we denote H^fast = H₁ + H₂ and H^slow = H₃ + H₄. For example, consider the solar system, setting all the bodies to be point masses except the Earth, ε ≈ 10⁻⁶.

The flows $\left\{{\varphi }_{t}^{[1]},{\varphi }_{t}^{[2]},{\varphi }_{t}^{[3]},{\varphi }_{t}^{[4]}\right\}$ of $\left\{{H}_{1},{H}_{2},{H}_{3},{H}_{4}\right\}$ can be derived similarly to Section 3.3.1 and the schemes are built by hierarchically composing ${\left\{{\varphi }_{t}^{[i]}\right\}}_{i=1}^{4}$ together. Specifically, as shown in Figure 3, we first group the flows of the fast dynamics (${\varphi }_{t}^{[1]}$ and ${\varphi }_{t}^{[2]}$) as a sub-scheme ${\varphi }_{h}^{\mathrm{fast}}$ via ${{ \mathcal C }}_{\mathrm{fast}}$ and the flows of the slow dynamics (${\varphi }_{t}^{[3]}$ and ${\varphi }_{t}^{[4]}$) as a sub-scheme ${\varphi }_{h}^{{\rm{s}}{low}}$ via ${{ \mathcal C }}_{\mathrm{slow}}$, respectively. Then, composing ${\varphi }_{h}^{{\rm{f}}{ast}}$ and ${\varphi }_{h}^{{\rm{s}}{low}}$ together as the final scheme ${\varphi }_{h}^{\mathrm{multi}}$ via ${{ \mathcal C }}_{\mathrm{multi}}$. ${{ \mathcal C }}_{\mathrm{fast}}$ and ${{ \mathcal C }}_{\mathrm{slow}}$ are composition methods of composing two Hamiltonian flows with comparable scales (McLachlan & Quispel 2002). ${{ \mathcal C }}_{\mathrm{multi}}$ is a composition method that is specialized in perturbative Hamiltonian systems of the form H = A + ε B (McLachlan 1995; Laskar & Robutel 2001; Blanes et al. 2013). Note that the flows ${\varphi }_{h}^{\mathrm{fast}}$, ${\varphi }_{h}^{\mathrm{slow}}$ are not exact, so the order of ${\varphi }_{h}^{\mathrm{multi}}$ is not the same as the order of ${{ \mathcal C }}_{\mathrm{multi}}$ applied for exact flows. In fact, the global error of ${\varphi }_{h}^{\mathrm{multi}}$ is the summation of the global errors of all three methods ${{ \mathcal C }}_{\mathrm{fast}}$, ${{ \mathcal C }}_{\mathrm{slow}}$ and ${{ \mathcal C }}_{\mathrm{multi}}$ (see Appendix C for the proof).

Zoom In Zoom Out Reset image size

For example, if we set ${{ \mathcal C }}_{\mathrm{fast}}$, ${{ \mathcal C }}_{\mathrm{slow}}$ and ${{ \mathcal C }}_{\mathrm{multi}}$ as ${{ \mathcal C }}_{S6}$, ${{ \mathcal C }}_{\mathrm{Verlet}}$ and ${{ \mathcal C }}_{{ABA}42}$ (see Appendix D), respectively. The global error of the above method is ${ \mathcal O }({h}^{6})+{ \mathcal O }({\varepsilon }^{2}{h}^{2})+{ \mathcal O }(\varepsilon {h}^{4}+{\varepsilon }^{2}{h}^{2})$, i.e., ${ \mathcal O }\left({h}^{6}+\varepsilon {h}^{4}+{\varepsilon }^{2}{h}^{2}\right)$. We name this the ${{ \mathcal M }}_{642}$ scheme with 6, 4, 2 representing the power of h of each term in the order and ${ \mathcal M }$ representing multiscale splitting.

Similarly, we design the ${{ \mathcal M }}_{42}$ scheme by choosing ${{ \mathcal C }}_{\mathrm{fast}}$, ${{ \mathcal C }}_{\mathrm{slow}}$ and ${{ \mathcal C }}_{\mathrm{multi}}$ as ${{ \mathcal C }}_{\mathrm{TriJump}}$, ${{ \mathcal C }}_{\mathrm{Verlet}}$ and ${{ \mathcal C }}_{{ABA}22}$, respectively, with the global error being ${ \mathcal O }({h}^{4}+\varepsilon {h}^{2})$.

Compared with schemes in Section 3.3.1, tailored splitting is able to mix the fast and slow flows flexibly, and thus is able to control the time complexity. In fact, $T({\phi }_{h}^{[1]})=T({\varphi }_{h}^{[1]})+T({\varphi }_{h}^{[3]})$, $T({\phi }_{h}^{[2]})=T({\varphi }_{h}^{[2]})+T({\varphi }_{h}^{[4]})$ with T( · ) being the number of operations of the one-step forward flow and evolving ${\varphi }_{h}^{[3]}$, ${\varphi }_{h}^{[4]}$ are much more expensive than evolving ${\varphi }_{h}^{[1]}$, ${\varphi }_{h}^{[2]}$. Meanwhile, ${\varphi }_{h}^{[3]}$ and ${\varphi }_{h}^{[4]}$ are (expensive) slow dynamics that can be evolved with less effort (e.g., larger step size, less stages) than fast dynamics when evolving together and tailored splitting makes it possible to control the number of expensive stages. To compare, the number of expensive stages and the global errors of all schemes mentioned (in Section 3.3.1 and Section 3.3.2) are listed in Table 1. In Table 1, the order index (o₀, o₁, ...) represents the power of h in front of ε⁰, ε¹, ... (e.g., a scheme of order (o₀, o₁, o₂) has a global error of ${ \mathcal O }({h}^{{o}_{1}}+\varepsilon {h}^{{o}_{2}}+{\varepsilon }^{2}{h}^{{o}_{3}})$).

Table 1. Comparisons of Different Schemes with Respect to the Number of Dominating Expensive Stages and the Global Error Order

Scheme	Expensive Stages	Order
${{ \mathcal T }}_{2}$	3	(2)
${{ \mathcal T }}_{4}$	7	(4)
${{ \mathcal T }}_{6}$	15	(6)
${{ \mathcal M }}_{42}$	3	(4, 2)
${{ \mathcal M }}_{642}$	6	(6, 4, 2)

Note. The number of expensive stages are counted in an isolated step without considering the concatenation of the last stage with the first stage of the next step. The notation in the "order" column is explained in the main text.

Download table as:  ASCIITypeset image

Moreover, since the hierarchical composition is a general framework, one can easily extend the family of numerical schemes, such as to construct higher order schemes, by applying a variety of existing splitting and composition methods.

3.3.3. Tailored Splitting II: H = K₁ + K₂ + K₃ with $\tfrac{{K}_{3}}{{K}_{1}},\tfrac{{K}_{2}}{{K}_{1}}={ \mathcal O }({\varepsilon }_{K})$

We also provide an option to use the popular Wisdom-Holman (Wisdom & Holman 1991) scheme for the orbital part, which works well for the specific but common setup of near-Keplerian orbits. These systems usually correspond to N − 1 well-separated bodies orbiting around a massive central body (indexed by 1 in our following description). This method is similar to the approach by Touma & Wisdom (1994), except that their coordinates are set using the body frame and we provided a higher-order implementation.

By isolating the Keplerian dynamics as K₁, combining the rotational kinetic energy with the rest translational kinetic energy as K₂, and putting the rest potential energy to K₃, H = K₁ + K₂ + K₃ with

$$\begin{eqnarray}&&\left\{\begin{array}{l}{K}_{1}({\mathsf{q}},{\mathsf{p}})={H}_{{Kepler}}({\mathsf{Q}},{\mathsf{P}})=\sum _{i=2}^{N}\displaystyle \frac{1}{2}{{\boldsymbol{P}}}_{i}^{T}{{\boldsymbol{P}}}_{i}/{m}_{i}-\displaystyle \frac{{ \mathcal G }{m}_{1}{m}_{i}}{\parallel {{\boldsymbol{Q}}}_{i}\parallel },\\ {K}_{2}({\mathsf{p}},{\rm{\Pi }})=\sum _{i=1}^{N}\displaystyle \frac{1}{2}{{\boldsymbol{\Pi }}}_{i}^{T}{{\boldsymbol{J}}}_{i}^{-1}{{\boldsymbol{\Pi }}}_{i}+\displaystyle \frac{{\parallel {{\boldsymbol{p}}}_{1}-\tfrac{{m}_{1}}{{m}_{\mathrm{tot}}}\sum _{i=1}^{N}{{\boldsymbol{p}}}_{i}\parallel }^{2}}{2{m}_{1}},\\ {K}_{3}({\mathsf{q}},{\mathsf{R}})=V({\mathsf{q}},{\mathsf{R}})+\sum _{i=2}^{N}\displaystyle \frac{{ \mathcal G }{m}_{1}{m}_{i}}{\parallel {{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{1}\parallel },\end{array}\right.\end{eqnarray} \tag{ 23 }$$

and K₂, K₃ ≪ K₁. Here, V(q, R) is defined in Equation (7). Note that K₁ represents Keplerian orbits in Q, P variables, which are canonical democratic heliocentric variables (Duncan et al. 1998) with

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{i}=\left\{\begin{array}{cl}{{\boldsymbol{q}}}_{i}-{{\boldsymbol{q}}}_{1} & i\ne 1,\\ \displaystyle \frac{\sum _{j=1}^{N}{m}_{j}{{\boldsymbol{q}}}_{j}}{{m}_{\mathrm{tot}}} & i=1,\end{array}\right.\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}{{\boldsymbol{P}}}_{i}=\left\{\begin{array}{cl}{{\boldsymbol{p}}}_{i}-\displaystyle \frac{{m}_{i}}{{m}_{\mathrm{tot}}}\displaystyle \sum _{j=1}^{N}{{\boldsymbol{p}}}_{i} & i\ne 1,\\ \displaystyle \sum _{j=1}^{N}{{\boldsymbol{p}}}_{i} & i=1.\end{array}\right.\end{eqnarray} \tag{ 25 }$$

So when evolving K₁ dynamics, additional steps of switching back and forth between (q, p) and (Q, P) coordinates are necessary. In terms of compositions, similarly, we first compose the flows of K₂ and K₃ together as ${\varphi }_{h}^{K,{slow}}$ via ${{ \mathcal C }}_{\mathrm{slow}}^{K}$, we then compose the flow of K₁ (${\varphi }_{h}^{K,\mathrm{fast}}$) with ${\varphi }_{h}^{K,\mathrm{slow}}$ via a multiscale compositing method ${{ \mathcal C }}_{\mathrm{multi}}^{K}$. The error of this composition is the summation of the global errors of two methods ${{ \mathcal C }}_{\mathrm{slow}}^{K}$, ${{ \mathcal C }}_{\mathrm{multi}}^{K}$ (and the numerical error of evolving Keplerian orbits).

For instance, ${{ \mathcal K }}_{\cdot 2}$ method in our package is based on choosing ${{ \mathcal C }}_{\mathrm{slow}}^{K}$, ${{ \mathcal C }}_{\mathrm{multi}}^{K}$ as ${{ \mathcal C }}_{\mathrm{Verlet}}$ and ${{ \mathcal C }}_{{ABA}22}$, and its global error is ${ \mathcal O }({\varepsilon }_{K}{h}^{2})$.

3.3.4. Which One to Use, the ${ \mathcal T }$-series, the ${ \mathcal M }$-series, or the ${ \mathcal K }$-series Methods?

Non-conservative forces such as tidal forces and post-Newtonian corrections are incorporated in the package. Because the implemented schemes are based on symmetric splitting and composition, the corresponding non-conservative momentum update is inserted in the middle of the composition. This is similar to how dissipative forces were added in REBOUNDx (Tamayo et al. 2020).

3.4.1. Tidal Forces

We model the tidal dissipation between each pair of bodies using the constant time lag equilibrium tide model, following Hut (1981), Eggleton et al. (1998). Note that we only adopted the dissipative component in the tidal force here. The expression of the acceleration of the tidal force is

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{a}}}_{\mathrm{host},\mathrm{guest}}^{\mathrm{tidal}} & = & -\displaystyle \frac{9\sigma {m}_{\mathrm{guest}}^{2}{A}^{2}}{2{\mu }_{\mathrm{host},\mathrm{guest}}{d}^{10}}\left[3{\boldsymbol{d}}\left({\boldsymbol{d}}\cdot \dot{{\boldsymbol{d}}}\right)\right.\\ & & \left.+\left(\left({\boldsymbol{d}}\times \dot{{\boldsymbol{d}}}\right)-{\boldsymbol{\omega }}{d}^{2}\right)\times {\boldsymbol{d}}\right].\end{array}\end{eqnarray} \tag{ 26 }$$

Here, m_host, m_guest denote the masses of the host and the guest body, respectively; d = q _guest − q _host denotes the relative position of the guest body; $d=\parallel {\boldsymbol{d}}\parallel$ denotes the distance between two bodies; ${\mu }_{\mathrm{host},\mathrm{guest}}=\tfrac{{m}_{\mathrm{host}}\cdot {m}_{\mathrm{guest}}}{{m}_{\mathrm{host}}+{m}_{\mathrm{guest}}}$ denotes the reduced mass; ω denotes the angular velocity of the host body under the reference frame (the inertia frame); the constant σ denotes the dissipation rate; A is defined as

$$\begin{eqnarray}&&A=\displaystyle \frac{{d}^{5}Q}{1-Q},\end{eqnarray} \tag{ 27 }$$

with Q the constant that measures quadrupolar deformability of the objects.

The dissipation rate σ is related to the time lag τ by the following formula,

$$\begin{eqnarray}&&\tau =\displaystyle \frac{3\sigma {d}^{5}}{4{ \mathcal G }}\cdot \displaystyle \frac{Q}{1-Q}.\end{eqnarray} \tag{ 28 }$$

We may integrate the tidal acceleration ${{\boldsymbol{a}}}_{i,j}^{\mathrm{tidal}}$ to our integrator after each time step by considering all pairs of bodies under tidal interactions. Note that each ${{\boldsymbol{a}}}_{i,j}^{\mathrm{tidal}}$ only calculates the force of each (host, guest) pair, where each pair (i, j) treat i as the extended object and j as the point-mass object. Thus, the equations of motion due to tidal dissipation are listed below:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{p}}}_{i}={{\boldsymbol{p}}}_{i}+h\displaystyle \sum _{j\ne i}\left(-{\mu }_{i,j}\,{{\boldsymbol{a}}}_{i,j}^{\mathrm{tidal}}+{\mu }_{j,i}\,{{\boldsymbol{a}}}_{j,i}^{\mathrm{tidal}}\right),\\ {{\boldsymbol{\Pi }}}_{i}={{\boldsymbol{\Pi }}}_{i}-h\sum _{j\ne i}{\mu }_{i,j}\,{{\boldsymbol{R}}}_{i}^{T}\left(\left({{\boldsymbol{q}}}_{j}-{{\boldsymbol{q}}}_{i}\right)\times {{\boldsymbol{a}}}_{i,j}^{\mathrm{tidal}}\right).\end{array}\right.\end{eqnarray} \tag{ 29 }$$

3.4.2. General Relativistic Effects

We added the first-order post-Newtonian correction for general relativistic effects following (for example) Blanchet (2006). For planetary systems, we assumed that the central object (the host star) is much more massive when compared to the surrounding objects (the planets). Thus, we only included the correction due to the star. The acceleration can be expressed as follows (e.g., Anderson et al. 1975; Benitez & Gallardo 2008):

$$\begin{eqnarray}&&{\boldsymbol{a}}=\displaystyle \frac{{{GM}}_{\mathrm{star}}}{{r}^{3}{c}^{2}}\left[\left(\displaystyle \frac{4{{GM}}_{\mathrm{star}}}{r}-{{\boldsymbol{v}}}^{2}\right){\boldsymbol{r}}+4({\boldsymbol{v}}\cdot {\boldsymbol{r}}){\boldsymbol{v}}\right]\end{eqnarray} \tag{ 30 }$$

For planets with close-in orbits, both rotational flattening and tidal force distort the shape of the planets, and lead to non-axial symmetric distortions. Thus, we include the option to study non-axial symmetric planets here, where one could specify the principal moment of inertia or the semiaxes of the planets directly. In this case, ${J}_{i}^{(1)}\ne {J}_{i}^{(2)}\ne {J}_{i}^{(3)}$ in J _i , and our splitting of the Hamiltonian is modified as the previous Hamiltonian plus H_asymmetric, where

$$\begin{eqnarray}&&{H}_{\mathrm{asymmetric}}({\mathsf{R}},{\rm{\Pi }})=\sum _{i}\left(\displaystyle \frac{1}{{J}_{i}^{(2)}}-\displaystyle \frac{1}{{J}_{i}^{(1)}}\right)\cdot \displaystyle \frac{{{{\rm{\Pi }}}_{i}^{(2)}}^{2}}{2},\end{eqnarray} \tag{ 31 }$$

and H_asymmetric ≪ H₃ in Equation (22).

The dynamics of Equation (31) are

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{R}}}_{i}(t)={R}_{y}\left(\delta {{\rm{\Pi }}}_{i}^{(2)}t\right)\,{{\boldsymbol{R}}}_{i}(0),\\ {{\boldsymbol{\Pi }}}_{i}(t)={R}_{y}(-\delta {{\rm{\Pi }}}_{i}^{(2)}t)\,{{\boldsymbol{\Pi }}}_{i}(0),\end{array}\right.\end{eqnarray} \tag{ 32 }$$

with $\delta =\tfrac{1}{{J}_{i}^{(2)}}-\tfrac{1}{{J}_{i}^{(1)}}$. Based on the symmetric schemes in Section 3.3.2, we simply evolve H_asymmetric half step at the beginning and the end of each step.

4.1.1. Conservation Properties

The conservation properties of the integrators are tested for ${{ \mathcal T }}_{4}$ and ${{ \mathcal M }}_{42}$ schemes in the Sun–Earth–Moon system with all three bodies being rigid. As shown in Figure 4, both schemes conserve linear momentum and angular momentum (except that there are arithmetic inaccuracies due to machine precision), and the energies exhibit no drift but only fluctuate at magnitudes ${ \mathcal O }({h}^{4})$ and ${ \mathcal O }({h}^{4}+\varepsilon {h}^{2})$ for ${{ \mathcal T }}_{4}$ and ${{ \mathcal M }}_{42}$, respectively. In the simulations, tides are not included (otherwise the system is no longer conservative) and initial conditions are set to be the data of epoch J2000 from JPL HORIZONS System.

Zoom In Zoom Out Reset image size

Here, the floating-point format is set to double precision, although our package can also use long-double or single-precision.

Our integrators also (exactly) preserve symplecticity when tidal dissipation is excluded because they are Hamiltonian splitting schemes. The definition of symplecticity in a non-Euclidean setup is not completely trivial, but the symplecticity of splitting approaches considered here has been established in (for example) Tao & Ohsawa (2020) (with r(t) = 0; otherwise one gets a more general result, namely conformal symplecticity).

4.1.2. Convergence Tests and Accuracy Comparisons

We now numerically illustrate how the integration error depends on h for different numerical schemes, which include both the methods that we implemented in GRIT and SMERCURY-T. SMERCURY-T is a concurrent simulation package that can evolve an object's spin axis under obliquity tide (Kreyche et al. Submitted), which is based on the Mercury simulation package (Chambers 1999). Specifically, it includes a subroutine to evolve the spin-axis dynamics following the procedure outlined in Lissauer et al. (2012), which is based on the Lie–Poisson integrator of rigid-body dynamics developed by Touma & Wisdom (1994). In addition, it includes a subroutine for obliquity tide following the algorithms outlined in Bolmont et al. (2015). The model for tidal interaction of SMERCURY-T is different from what we included in GRIT, which naturally contains both obliquity tide and tidal effects due to non-tidally synchronized orbits. Thus, we focus here on the rigid-body dynamics, where we do not include tidal interactions in our convergence test. In the comparisons presented here, we also turned off our rigid-body rigid-body interaction option (which is mainly for accurate simulations of rigid bodies' close encounters) because these interactions are supported only in GRIT.

We first test on the Sun–Earth–Moon system (Figure 5). One observation in this case is that if the step size is too large so that splitting into H₁ + H₂ + H₃ + H₄ (GRIT's ${{ \mathcal M }}_{42}$, ${{ \mathcal M }}_{642}$) does not work, then SMERCURY-T does not work either (unlike expected by some). More precisely, with h = 2 · 10⁻² yr, ${{ \mathcal M }}_{42}$ and SMERCURY-T cannot resolve the the motion of the Moon orbiting around the Earth, whose period is a month, and even the performance of the 6th-order method ${{ \mathcal M }}_{642}$ is not ideal and significant errors are observed in all methods. Accuracy is improved for step sizes below this stability limit, and the rate of improvement is (as expected) dependent on the order of the numerical scheme. Consequently, higher order methods, such as ${{ \mathcal M }}_{42}$ and ${{ \mathcal M }}_{642}$, show substantially smaller errors when smaller step sizes are applied (readers interested in understanding this together with computational costs are referred to Section 4.1.3).

Zoom In Zoom Out Reset image size

We then test on a non-Keplerian system (note SMERCURY-T performs well for near-Keplerian problems as designed) using an Earth-like planet orbiting around two stars alternatively in a stellar binary system (Figure 6). Given that there is no single body that has the dominant mass of the system and the planet is alternatively captured by the two stars, the planetary orbit is not nearly Keplerian, and splitting into H₁ + H₂ + H₃ + H₄ is more accurate than SMERCURY-T for all choices of step size here. Specifically, as shown in Figure 6, the orbital position of SMERCURY-T saturates to ${ \mathcal O }(1)$ relative error after a relatively short period of time, no matter if h = 10⁻³, 10⁻⁴, or 10⁻⁵ yr. The orbital inaccuracy naturally also affects the spin angle. Meanwhile, ${{ \mathcal M }}_{42}$ and ${{ \mathcal M }}_{642}$ do not have this issue.

Zoom In Zoom Out Reset image size

For reproducibility, the initial condition used is

$$\begin{eqnarray*}\begin{array}{ll}{{\boldsymbol{q}}}_{{\mathrm{star}}_{1}}={\left[-0.5\ \ 0\ \ 0\right]}^{T}, & {{\boldsymbol{v}}}_{{\mathrm{star}}_{1}}={\left[0\ \ -0.0086012119\ \ 0\right]}^{T},\\ {{\boldsymbol{q}}}_{{\mathrm{star}}_{2}}={\left[0.5\ \ 0\ \ 0\right]}^{T}, & {{\boldsymbol{v}}}_{{\mathrm{star}}_{2}}={\left[0\ \ 0.0086012119\ \ 0\right]}^{T},\\ {{\boldsymbol{q}}}_{\mathrm{planet}}={\left[1.16\ \ 0\ \ 0\right]}^{T}, & {{\boldsymbol{v}}}_{\mathrm{planet}}={\left[0\ \ 0.0164271047\ \ 0\right]}^{T}\end{array}\end{eqnarray*}$$

in units of au and au/day, and ${m}_{{\mathrm{star}}_{1}}={m}_{{\mathrm{star}}_{2}}=0.5{m}_{\odot }$, m_planet = m_⊕.

4.1.3. Investigation of Efficiency

We now demonstrate the improved computational efficiency of the tailored splitting schemes. A comparison of the time efficiency among the traditional splitting method ${{ \mathcal T }}_{4},{{ \mathcal T }}_{6}$ and the tailored splitting scheme ${{ \mathcal M }}_{42},{{ \mathcal M }}_{642},{{ \mathcal K }}_{\cdot 2}$ in the 10 rigid body (Sun with 8 planets and the Moon) is shown in Table 2. ${{ \mathcal M }}_{42}$ (${{ \mathcal M }}_{642}$) is about twice the speed of ${{ \mathcal T }}_{4}$ (${{ \mathcal T }}_{6}$) with comparable integrating accuracy. Note that SMERCURY-T cannot be compared against here because its currently available version ⁶ can only set one of the objects as a rigid-body.

Table 2. Efficiency Comparison Among Scheme ${{ \mathcal T }}_{4}$, ${{ \mathcal T }}_{6}$, ${{ \mathcal M }}_{42}$, ${{ \mathcal M }}_{642}$ and ${{ \mathcal K }}_{\cdot 2}$

h = 10−3 yr	Wall Time (s)	MAE of Earth's Obliquity
${{ \mathcal T }}_{4}$	30.573	1.996646e-05
${{ \mathcal M }}_{42}$	15.488	1.997454e-05
${{ \mathcal T }}_{6}$	72.55	1.728156e-08
${{ \mathcal M }}_{642}$	40.626	4.365093e-10
${{ \mathcal K }}_{\cdot 2}$	14.673	2.186379e-05
SMERCURY-T	N/A	N/A
h = 10−4 yr
${{ \mathcal T }}_{4}$	273.71	2.680897e-09
${{ \mathcal M }}_{42}$	140.58	3.817091e-09
${{ \mathcal T }}_{6}$	708.26	8.689218e-11
${{ \mathcal M }}_{642}$	395.52	2.039980e-10
${{ \mathcal K }}_{\cdot 2}$	131.56	1.292609e-05
SMERCURY-T	N/A	N/A
h = 10−5 yr
${{ \mathcal K }}_{\cdot 2}$	1299.9	1.378428e-07

Note. The Solar system with eight planets and the Moon (10 rigid bodies in total) is simulated until 1000 yr with h = 10⁻³ yr and h = 10⁻⁴ yr for all schemes using a single thread. The benchmark is simulated using the ${{ \mathcal T }}_{6}$ scheme with h = 10⁻⁵ yr and long-double precision. Mean absolute errors (MAE) of the Earth's obliquity (rad) are measured. Data is output every 0.1 yr.

Download table as:  ASCIITypeset image

To gain additional understanding of the performance of GRIT, complementary results that include comparisons to SMERCURY-T are also provided. For a fair comparison, we continue using the Solar system example, which is a near-Keplerian problem that SMERCURY-T specializes in, but we had to alter it by setting only the Earth to be a rigid body and all others as point masses. The results are given in Table 3, where ${{ \mathcal M }}_{42}$ shows improved accuracy over SMERCURY-T, while ${{ \mathcal M }}_{642}$ is even more accurate but with traded-off time complexity.

Table 3. Efficiency Comparison Among Scheme ${{ \mathcal M }}_{42}$, ${{ \mathcal M }}_{642}$ and SMERCURY-T

h = 10−3 yr	Wall Time (s)	MAE of Earth's Obliquity
${{ \mathcal M }}_{42}$	6.408	1.997119e-05
${{ \mathcal M }}_{642}$	23.122	3.841649e-10
SMERCURY-T	8.638	2.157662e-05
h = 10−4 yr
${{ \mathcal M }}_{42}$	53.782	3.833661e-09
${{ \mathcal M }}_{642}$	216.09	1.990336e-10
SMERCURY-T	39.079	1.903458e-05

Note. The Solar system with eight planets and the Moon (nine point masses and one rigid body (the Earth) in total) is simulated until 1000 yr with h = 10⁻³ yr and h = 10⁻⁴ yr for all schemes using a single thread. The benchmark is simulated using the ${{ \mathcal T }}_{6}$ scheme with h = 10⁻⁵ yr and long-double precision. Mean absolute errors (MAE) of the Earth's obliquity (rad) are measured. Data is output every 0.1 yr.

Download table as:  ASCIITypeset image

In addition, for the sake of fairness, note that wall-clock counts are platform dependent and therefore should only be used as a qualitative (not quantitative) indicator. The experiments reported here are conducted on a machine with an AMD Ryzen 7 3700X 8-Core Processor, 16 GB memory and the Linux distribution of openSUSE Leap 15.2. GRIT was compiled using GNU C++ compiler and SMERCURY-T using GNU Fortran compiler, both with the default compilation options. A single thread is used for experiments in both Tables 2 and 3 for fairness (note that a parallelization option is available in GRIT and we recommend turning it on when the simulated system has large numbers of rigid objects). We also noted that SMERCURY-T slows down more significantly than GRIT when its integration is outputted more frequently, and thus chose a large output step size to reduce SMERCURY-T's I/O overhead so that the focus can be on the integration time itself.

4.1.4. Summary of the Numerical Tests in Section 4.1

To further verify the accuracy of our integration package, we compare our simulation results to secular theory here. We include two examples: the first example integrate the obliquity variation of a moonless Earth without the influence of tidal interactions, and the second example considers tidal interactions between a hypothetical Earth–Moon system. We find good agreement between our simulation package with the results of the secular theory.

4.2.1. Obliquity Variations of a Moonless Earth

Spin-orbit resonances lead to large obliquity variations for a moonless Earth (Laskar et al. 1993). This classical example can serve as a test case for our simulation package. Specifically, planetary companions of the Earth (from Mercury to Neptune) all perturb Earth's orbit and lead to forced oscillations in the orbital plane of the Earth. At the same time, torquing from the Sun leads to precession of the Earth's spin axis. The natural precession frequency coincides with the forcing frequencies and drives resonant obliquity variations of the Earth. Tidal interactions are weak in this case, so we neglected tidal effects in our code, and consider the dynamical coupling between the planetary spin axes and its orbit.

We include the eight Solar System planets in this system, and we adopt the position and velocity of the Solar System planets from JPL database (Giorgini et al. 1996). We only treat the Earth as a rigid object with oblateness of 0.00335, and we set the other planets and the Sun as point particles.

Figure 7 shows the comparison of the obliquity variations of the moonless Earth with that from the secular theory shown in Laskar et al. (1993) and Li & Batygin (2014). We included three examples, starting with different initial obliquities, and all of them show good agreement with the secular results. In particular, below ∼40°, large obliquity variations can be seen, due to the spin–orbit resonances. We chose a time step of 10⁻⁴ yrs to resolve the spin of the Earth. The fractional change in energy is of the order of 10⁻¹⁴ and the fractional change in angular momentum is of the order of 10⁻¹² for all the three runs with different initial obliquities.

Zoom In Zoom Out Reset image size

4.2.2. Tidal Interactions of a Hypothetical Earth–Moon System

To illustrate the accuracy of our simulation package including tidal interactions, we use a simple hypothetical Earth–Moon two-body system here. We set the initial semimajor axis and eccentricity to be 0.0018au and 0.4. For the Earth, we set the spin period to 1 day, oblateness to 0.00335, love number to 0.305 and tidal time lag to 698 sec. For the Moon, we set the spin period to 14 days, oblateness to 0.0012, love number to 0.02416 and tidal time lag to 8639 s.

Figure 8 shows the agreement between our simulation package (solid lines) with the secular results (dashed lines). The secular results are obtained following Eggleton et al. (1998). The upper panel plots orbital eccentricity versus time and the lower panel plots the spin rate of the Moon versus time. This shows that the spin rate of the Moon increases to the pseudo-synchronized state within a few hundred years, and then slowly decreases as orbital eccentricity decays due to tide. We chose a time step of 10⁻⁴ yr to resolve the spin of the Earth, and the total fractional change in angular momentum is 7 × 10⁻¹².

Zoom In Zoom Out Reset image size

We use the same orbital initial condition and physical properties for the planets in Trappist-I following Bolmont et al. (2020) (Table A.2 in Bolmont et al. 2020) to compare the magnitude of TTVs, and we use the same reference tidal parameters for the star and the planet (e.g., k_2f,* = 0.307, k_2f,p = 0.9532, Δτ_p = 712.37 sec). The Q coefficient in the tidal model can then be calculated (k₂ = Q/(1 − Q)) (Eggleton et al. 1998; Eggleton & Kiseleva-Eggleton 2001; Fabrycky & Tremaine 2007).

To calculate the moment of inertia along the three principal axes (A, B, C), we follow the derivation by Van Hoolst et al. (2008), while assuming a homogeneous model for simplicity and assuming that the rotation velocity of the planet is close to the orbital velocity. Specifically, the moment of inertia can be expressed as follows:

$$\begin{eqnarray*}\begin{array}{rcl}A & = & {I}_{0}\left(1-\displaystyle \frac{1}{3}\alpha -\displaystyle \frac{1}{2}\beta \right)\\ B & = & {I}_{0}\left(1-\displaystyle \frac{1}{3}\alpha +\displaystyle \frac{1}{2}\beta \right)\\ C & = & {I}_{0}\left(1+\displaystyle \frac{2}{3}\alpha \right)\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}\alpha & = & \displaystyle \frac{5}{4}q(1+{k}_{f})\\ \beta & = & \displaystyle \frac{3}{2}q(1+{k}_{f})\end{array}\end{eqnarray*}$$

and k_f is the love number, and q is the ratio of the centrifugal acceleration to the gravitational acceleration. We assume all the planets have the same radius of gyration squared ${{rg}}_{p}^{2}=0.3308$ following Bolmont et al. (2020), and we include in Table 4 the moment of inertia of the planets.

Moreover, because the planets are very close to their host star, general relativistic precession plays a non-negligible role in the transit time. Thus, we also included the first order post-Newtonian correction in our simulation code (see Section 3.4.2).

The measurement of TTVs is a powerful method to derive the physical properties of planets, in particular masses and eccentricity of planets (Agol & Fabrycky 2018). Although most studies consider only point-mass dynamics, full-body dynamics including tidal effects and distortion of the planets could also play an important role (Miralda-Escudé 2002; Heyl & Gladman 2007; Ragozzine & Wolf 2009; Maciejewski et al. 2018). It was recently shown that new measurements of the TTV of the Trappist-I system lead to significant increase in the mass estimate for planets b and c, which may be due to unaccounted physical processes including tidal effects and rotational distortion of the planets (Grimm et al. 2018; Agol et al. 2020; Bolmont et al. 2020). Thus, we use our simulation package to estimate the TTV of the inner planets in Trappist-I as an example, in comparison with the study by Bolmont et al. (2020).

We include the result of the TTVs for Trappist-I b,c and d in 1500 days in Figure 9, to compare our results with those in Bolmont et al. (2020). Similar to Figure 1 in Bolmont et al. (2020), the upper panels show the TTVs assuming that the planets are all point-mass particles, and the lower panels show the differences in the TTVs due to different effects. The differences due to GR and rotational flattening of the planets computed using our simulation package are agreeable with those in Bolmont et al. (2020). In contrast to Bolmont et al. (2020), we assume the objects are rigid bodies when considering tidal interactions with the central star using our rigid-body simulator. This leads to slightly larger TTV differences. We note that the magnitude of the differences in the TTVs depend on the misalignment between the elongated principal axis and the location direction of the planet from the central star. For the illustrative example, we assume that the planets all start with their long axes perfectly aligned to the direction of the central star.

Zoom In Zoom Out Reset image size

As the system evolves further, the misalignment could be excited to larger values (as discussed further in Section 5.3). The differences in TTVs could reach ∼5 s for 1500 days, and a few minutes in 10-year measurements, as shown in Figure 10. A detailed study of how the TTVs depend on the physical properties of the planets (e.g., the love number, tidal time lag, etc.) is out of the scope of this paper, and will be discussed in a follow up project.

Zoom In Zoom Out Reset image size

The long-term dynamics of the spin axes of planets, in particular their synchronized states, play an important role in the atmosphere circulation of the planets. When the planets are tidally locked, the extreme temperature differences on one side of the planet facing the star from the other side may lead to the collapse of the planetary atmosphere (Kasting et al. 1993; Joshi et al. 1997; Wordsworth 2015). For Trappist-I, Vinson et al. (2019) developed a framework studying the spin-axis variations of the planets and found that the mean-motion resonant chain could drive the habitable zone planets out of the synchronized state.

Specifically, Vinson et al. (2019) evolves the longitude of the substellar point separately based on results of orbital evolution of Trappist-I using the Rebound simulation package (Tamayo et al. 2017). This does not include effects of the variation of the spin axis on the orbits and the developed framework neglected the 3-D variations of the planetary spin axis (i.e., assuming zero planetary obliquities) for simplicity. We use our simulation package to evaluate the spin-axis dynamics more accurately, which allows backreactions of the spin-axis dynamics on the orbit, as well as the full 3-D dynamics of the planetary spin axis.

We use the same initial condition as those in Section 5.2 for the long-term dynamical simulation over 100,000 yr. We start the planets in synchronized configurations, and we calculate the misalignment between the long axes of the planets and their radial direction from the host star. ψ is illustrated in Figure 11.

Zoom In Zoom Out Reset image size

Figure 12 shows this misalignment (ψ) of the planets. Planets b, c, d and e are closer to the host star, and allow stronger tidal interactions. This leads to low variations in the long axes of the planets. However, planets f, g and h are further away, where planetary interactions could compete with tidal re-alignment and drive larger spin-axis variations. We note that the obliquities of these planets still remain low (within a few degrees). The detailed dependence of the spin-axis variations on the parameters of the planets are beyond the scope of this article, and will be invested in a follow up paper.

Zoom In Zoom Out Reset image size

With ${\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}\in {{\mathbb{R}}}^{3}$, ${\boldsymbol{D}}=\left[\begin{array}{ccc}{d}_{1} & 0 & 0\\ 0 & {d}_{2} & 0\\ 0 & 0 & {d}_{3}\end{array}\right]$, we have

• 1.  
  $\hat{{\boldsymbol{u}}}{\boldsymbol{v}}={\boldsymbol{u}}\times {\boldsymbol{v}}$.
• 2.  
  $\widehat{{\boldsymbol{u}}\times {\boldsymbol{v}}}=\hat{{\boldsymbol{u}}}\hat{{\boldsymbol{v}}}-\hat{{\boldsymbol{v}}}\hat{{\boldsymbol{u}}}$.
• 3.  
  $\hat{{\boldsymbol{u}}}{\boldsymbol{D}}-{\boldsymbol{D}}{\hat{{\boldsymbol{u}}}}^{T}={Tr}\left[{\boldsymbol{D}}\right]\hat{{\boldsymbol{u}}}-\widehat{{\boldsymbol{D}}{\boldsymbol{u}}}$.
• 4.  
  ${\hat{{\boldsymbol{u}}}}^{T}\hat{{\boldsymbol{u}}}{\boldsymbol{D}}-{\boldsymbol{D}}{\hat{{\boldsymbol{u}}}}^{T}\hat{{\boldsymbol{u}}}=\widehat{{\boldsymbol{u}}\times {\boldsymbol{D}}{\boldsymbol{u}}}$

We can view R to be in the embedded Euclidean space ${{\mathbb{R}}}^{3\times 3}\hookleftarrow {\mathsf{SO}}(3)$ and use R ∈ SO(3) as a holonomic constraint. The Lagrangian L (Equation (11)) has 9-DOF before applying the constraint R ∈ SO(3). The conjugate variable of R (t) will be denoted by P (t).

The constraint of a system forces the evolution of the system in a specific manifold, and the manifold can be directly calculated from the constraint (one may refer Chapter VII of Hairer et al. (2006) for details). For rigid-body dynamics represented by a rotation matrix R (t), the constraint is R (t)^T R (t) − I _3×3 = 0_3×3. Reich & Zentrum (1996) and Hairer et al. (2006) have shown the procedure of finding equations of motion by utilizing the constraint for a rigid-body system with a R dependent potential. By using Lagrange multipliers (Hairer et al. 2006) for the constraint R ^T R − I _3×3 = 0, we have the following Lagrangian,

$$\begin{eqnarray}\begin{array}{rcl}L\left({\boldsymbol{R}},\dot{{\boldsymbol{R}}}\right) & = & \displaystyle \frac{1}{2}{Tr}\left[\dot{{\boldsymbol{R}}}{{\boldsymbol{J}}}_{d}{\dot{{\boldsymbol{R}}}}^{T}\right]-V({\boldsymbol{R}})\\ & & -\displaystyle \frac{1}{2}{Tr}\left({{\boldsymbol{\Lambda }}}^{T}\left({{\boldsymbol{R}}}^{T}{\boldsymbol{R}}-{{\boldsymbol{I}}}_{3\times 3}\right)\right),\end{array}\end{eqnarray} \tag{ B1 }$$

with 6-dim Lagrange multipliers ${\boldsymbol{\Lambda }}=\left[\begin{array}{ccc}{\lambda }_{1} & {\lambda }_{4} & {\lambda }_{6}\\ {\lambda }_{4} & {\lambda }_{2} & {\lambda }_{5}\\ {\lambda }_{6} & {\lambda }_{5} & {\lambda }_{3}\end{array}\right]\in {{\mathbb{R}}}^{3\times 3}$ a symmetric matrix.

Doing Legendre transform for Equation (B1), we have

$$\begin{eqnarray}&&{\boldsymbol{P}}=\displaystyle \frac{\partial L\left({\boldsymbol{R}},\dot{{\boldsymbol{R}}}\right)}{\partial \dot{{\boldsymbol{R}}}}=\dot{{\boldsymbol{R}}}{{\boldsymbol{J}}}_{d},\end{eqnarray} \tag{ B2 }$$

and the corresponding Hamiltonian,

$$\begin{eqnarray}\begin{array}{rcl}H({\boldsymbol{R}},{\boldsymbol{P}}) & = & \displaystyle \frac{1}{2}{Tr}\left[{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}{{\boldsymbol{P}}}^{T}\right]+V({\boldsymbol{R}})\\ & & +\displaystyle \frac{1}{2}{Tr}\left({{\boldsymbol{\Lambda }}}^{T}\left({{\boldsymbol{R}}}^{T}{\boldsymbol{R}}-{{\boldsymbol{I}}}_{3\times 3}\right)\right).\end{array}\end{eqnarray} \tag{ B3 }$$

As the constraint for R is R ^T R = I _3×3, according to Hairer et al. (2006), the constraint for P can be obtained by taking time derivative for R ^T R − I _3×3 = 0_3×3, i.e., ${{\boldsymbol{J}}}_{d}^{-1}{{\boldsymbol{P}}}^{T}{\boldsymbol{R}}+{{\boldsymbol{R}}}^{T}{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}\,={{\bf{0}}}_{3\times 3}$.

So,

$$\begin{eqnarray}&&\left\{\begin{array}{l}\dot{{\boldsymbol{R}}}=\displaystyle \frac{\partial H}{\partial {\boldsymbol{P}}}={\boldsymbol{P}}{{{\boldsymbol{J}}}_{d}}^{-1},\\ \dot{{\boldsymbol{P}}}=-\displaystyle \frac{\partial H}{\partial {\boldsymbol{R}}}=-\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\boldsymbol{R}}{\boldsymbol{\Lambda }},\end{array}\right.\end{eqnarray} \tag{ B4 }$$

on the manifold

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal M }=\left\{\left({\boldsymbol{R}},{\boldsymbol{P}}\right)| {{\boldsymbol{R}}}^{T}{\boldsymbol{R}}={{\boldsymbol{I}}}_{3\times 3},{{\boldsymbol{J}}}_{d}^{-1}{{\boldsymbol{P}}}^{T}{\boldsymbol{R}}\right.\\ \left.\ \ \,+\ {{\boldsymbol{R}}}^{T}{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}={{\bf{0}}}_{3\times 3}\right\}.\end{array}\end{eqnarray} \tag{ B5 }$$

Note that $\hat{{\boldsymbol{\Omega }}}={{\boldsymbol{R}}}^{T}{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}$ with Ω being the body's angular velocity. Taking time derivative for $\hat{{\boldsymbol{\Omega }}}$, we have

$$\begin{eqnarray}&&\hat{\dot{{\boldsymbol{\Omega }}}}={{\boldsymbol{J}}}_{d}^{-1}{{\boldsymbol{P}}}^{T}{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}+{{\boldsymbol{R}}}^{T}\left(-\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\boldsymbol{R}}{\boldsymbol{\Lambda }}\right){{\boldsymbol{J}}}_{d}^{-1}.\end{eqnarray} \tag{ B6 }$$

Physically, we want to find dynamics of R and the body's angular momentum Π. Since $\hat{{\boldsymbol{\Pi }}}=\widehat{{\boldsymbol{J}}{\boldsymbol{\Omega }}}={Tr}\left[{{\boldsymbol{J}}}_{d}\right]\hat{{\boldsymbol{\Omega }}}-\widehat{{{\boldsymbol{J}}}_{d}{\boldsymbol{\Omega }}}\,=\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}-{{\boldsymbol{J}}}_{d}{\hat{{\boldsymbol{\Omega }}}}^{T}$ (see Appendix A.1), we may find dynamics of Π,

$$\begin{eqnarray}\begin{array}{rcl}\widehat{\dot{{\boldsymbol{\Pi }}}} & = & \left({{\boldsymbol{J}}}_{d}^{-1}{{\boldsymbol{P}}}^{T}{\boldsymbol{P}}-{{\boldsymbol{P}}}^{T}{\boldsymbol{P}}{{\boldsymbol{J}}}_{d}^{-1}\right)\\ & & +\left({\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}-{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right),\end{array}\end{eqnarray} \tag{ B7 }$$

with the symmetric Λ vanished. ⁷

As ${\boldsymbol{P}}={\boldsymbol{R}}\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}$, properties of hat-map (see Appendix A.1) lead to

$$\begin{eqnarray}\begin{array}{rcl}\widehat{\dot{{\boldsymbol{\Pi }}}} & = & \left({\hat{{\boldsymbol{\Omega }}}}^{T}\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}-{{\boldsymbol{J}}}_{d}{\hat{{\boldsymbol{\Omega }}}}^{T}\hat{{\boldsymbol{\Omega }}}\right)+\left({\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}-{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)\\ & = & \widehat{{\boldsymbol{\Omega }}\times {{\boldsymbol{J}}}_{d}{\boldsymbol{\Omega }}}+\left({\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}-{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right).\end{array}\end{eqnarray} \tag{ B8 }$$

Thus

$$\begin{eqnarray}\begin{array}{rcl}\dot{{\boldsymbol{\Pi }}} & = & {\boldsymbol{\Omega }}\times {{\boldsymbol{J}}}_{d}{\boldsymbol{\Omega }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }\\ & = & {\boldsymbol{\Omega }}\times \left({Tr}[{{\boldsymbol{J}}}_{d}]-{\boldsymbol{J}}\right){\boldsymbol{\Omega }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }\\ & = & -{\boldsymbol{\Omega }}\times {\boldsymbol{J}}{\boldsymbol{\Omega }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }\\ & = & {\boldsymbol{\Pi }}\times {{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }\end{array}\end{eqnarray} \tag{ B9 }$$

So, the equations of motion with respect to R and Π for a one rigid-body system are

$$\begin{eqnarray}&&\left\{\begin{array}{l}\dot{{\boldsymbol{R}}}={\boldsymbol{R}}\widehat{{{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}},\\ \dot{{\boldsymbol{\Pi }}}={\boldsymbol{\Pi }}\times {{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}-{\left({{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}-{\left(\displaystyle \frac{\partial V\left({\boldsymbol{R}}\right)}{\partial {\boldsymbol{R}}}\right)}^{T}{\boldsymbol{R}}\right)}^{\vee }.\end{array}\right.\end{eqnarray} \tag{ B10 }$$

Obtaining a Euler–Lagrange equation for the Hamilton's variational principle on a Lie group has been well studied (e.g., Marsden & Ratiu 1994; Holm et al. 2009). Here, we summarize the results for the special case of rigid bodies from the expository part of Lee et al. (2005).

Denote the infinitesimally varied rotation by ${{\boldsymbol{R}}}_{\epsilon }={\boldsymbol{R}}\exp (\epsilon \hat{{\boldsymbol{\eta }}})$ with $\epsilon \in {\mathbb{R}}$ and ${\boldsymbol{\eta }}\in {{\mathbb{R}}}^{3}$, where $\exp (\cdot )$ is a mapping from ${\mathfrak{so}}(3)$ to SO(3). The varied angular velocity is

$$\begin{eqnarray}\begin{array}{rcl}{\hat{{\boldsymbol{\Omega }}}}_{\epsilon } & = & {{\boldsymbol{R}}}_{\epsilon }^{T}{\dot{{\boldsymbol{R}}}}_{\epsilon }={e}^{-\epsilon \hat{{\boldsymbol{\eta }}}}{{\boldsymbol{R}}}^{T}\left(\dot{{\boldsymbol{R}}}{e}^{\epsilon \hat{{\boldsymbol{\eta }}}}+{\boldsymbol{R}}\cdot {e}^{\epsilon \hat{{\boldsymbol{\eta }}}}\epsilon \hat{\dot{{\boldsymbol{\eta }}}}\right)\\ & = & {e}^{-\epsilon \hat{{\boldsymbol{\eta }}}}\hat{{\boldsymbol{\Omega }}}{e}^{\epsilon \hat{{\boldsymbol{\eta }}}}+\epsilon \hat{\dot{{\boldsymbol{\eta }}}}\\ & = & \hat{{\boldsymbol{\Omega }}}+\epsilon \left\{\hat{\dot{{\boldsymbol{\eta }}}}+\hat{{\boldsymbol{\Omega }}}\hat{{\boldsymbol{\eta }}}-\hat{{\boldsymbol{\eta }}}\hat{{\boldsymbol{\Omega }}}\right\}+{ \mathcal O }({\epsilon }^{2}).\end{array}\end{eqnarray} \tag{ B11 }$$

Consider the action

$$\begin{eqnarray}\begin{array}{rcl}S\left({\boldsymbol{\Omega }},{\boldsymbol{R}}\right) & = & {\displaystyle \int }_{{t}_{0}}^{{t}_{1}}L\left({\boldsymbol{\Omega }},{\boldsymbol{R}}\right)\,{dt}\\ & = & {\displaystyle \int }_{{t}_{0}}^{{t}_{1}}\displaystyle \frac{1}{2}{Tr}\left[\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}{\hat{{\boldsymbol{\Omega }}}}^{T}\right]-V({\boldsymbol{R}})\,{dt}.\end{array}\end{eqnarray} \tag{ B12 }$$

Taking the variation of the action S, we have

$$\begin{eqnarray}&&\begin{array}{l}{S}_{\epsilon }\left({\boldsymbol{\Omega }},{\boldsymbol{R}}\right)=S\left({{\boldsymbol{\Omega }}}_{\epsilon },{{\boldsymbol{R}}}_{\epsilon }\right)\\ \ \ =\ S\left({\boldsymbol{\Omega }},{\boldsymbol{R}}\right)+\epsilon {\displaystyle \int }_{{t}_{0}}^{{t}_{1}}\displaystyle \frac{1}{2}{Tr}\left[-\hat{\dot{{\boldsymbol{\eta }}}}\left({{\boldsymbol{J}}}_{d}\hat{{\boldsymbol{\Omega }}}+\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}\right)\right.\\ \ \ \,\left.+\ \hat{{\boldsymbol{\eta }}}\hat{{\boldsymbol{\Omega }}}\left({{\boldsymbol{J}}}_{d}\hat{{\boldsymbol{\Omega }}}+\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}\right)-\hat{{\boldsymbol{\eta }}}\left({{\boldsymbol{J}}}_{d}\hat{{\boldsymbol{\Omega }}}+\hat{{\boldsymbol{\Omega }}}{{\boldsymbol{J}}}_{d}\hat{{\boldsymbol{\Omega }}}\right)\right]\\ \ \ \,+\ {Tr}\left[\hat{{\boldsymbol{\eta }}}{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial {\boldsymbol{R}}}{\partial {\boldsymbol{R}}}\right]\,{dt}+{ \mathcal O }({\epsilon }^{2}).\end{array}\end{eqnarray} \tag{ B13 }$$

Using Hamilton's Principle, we have ${\left.\tfrac{d}{d\epsilon }\right|}_{\epsilon =0}{S}_{\epsilon }=0$, i.e.

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{\int }_{{t}_{0}}^{{t}_{1}}{Tr}\left[\hat{{\boldsymbol{\eta }}}\left\{\widehat{{\boldsymbol{J}}\dot{{\boldsymbol{\Omega }}}}+\widehat{{\boldsymbol{\Omega }}\times {\boldsymbol{J}}{\rm{\Omega }}}+2{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V}{\partial {\boldsymbol{R}}}\right\}\right]=0\end{eqnarray} \tag{ B14 }$$

for any ${\boldsymbol{\eta }}\in {{\mathbb{R}}}^{3}$. Therefore, $\left\{\widehat{{\boldsymbol{J}}\dot{{\boldsymbol{\Omega }}}}+\widehat{{\boldsymbol{\Omega }}\times {\boldsymbol{J}}{\rm{\Omega }}}+2{{\boldsymbol{R}}}^{T}\tfrac{\partial V}{\partial {\boldsymbol{R}}}\right\}$ must be skew-symmetric, which gives us

$$\begin{eqnarray}&&\widehat{{\boldsymbol{J}}\dot{{\boldsymbol{\Omega }}}}=-\widehat{{\boldsymbol{\Omega }}\times {\boldsymbol{J}}{\boldsymbol{\Omega }}}+\left({\displaystyle \frac{\partial V}{\partial {\boldsymbol{R}}}}^{T}{\boldsymbol{R}}-{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V}{\partial {\boldsymbol{R}}}\right).\end{eqnarray} \tag{ B15 }$$

Thus

$$\begin{eqnarray}&&\widehat{\dot{{\boldsymbol{\Pi }}}}=\widehat{{\boldsymbol{\Pi }}\times {{\boldsymbol{J}}}^{-1}{\boldsymbol{\Pi }}}+\left({\displaystyle \frac{\partial V}{\partial {\boldsymbol{R}}}}^{T}{\boldsymbol{R}}-{{\boldsymbol{R}}}^{T}\displaystyle \frac{\partial V}{\partial {\boldsymbol{R}}}\right).\end{eqnarray} \tag{ B16 }$$