<figure><img src="image_1.png" /></figure>

HAL
open science

# A direct N-body integrator for modelling the chaotic, tidal dynamics of multibody extrasolar systems: TIDYMESS

Tjarda C. N. Boekholt, Alexandre C. M. Correia

► To cite this version:

Tjarda C. N. Boekholt, Alexandre C. M. Correia. A direct N-body integrator for modelling the chaotic, tidal dynamics of multibody extrasolar systems: TIDYMESS. Monthly Notices of the Royal Astronomical Society, 2023, 522, pp.2885-2900. 10.1093/mnras/stad1133. insu-04849159

HAL Id: insu-04849159

https://insu.hal.science/insu-04849159v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A direct *N*-body integrator for modelling the chaotic, tidal dynamics of multibody extrasolar systems: TIDYMESS

Tjarda C. N. Boekholt Ⓔ1★ and Alexandre C. M. Correia Ⓔ2,3★

Accepted 2023 April 11. Received 2023 April 3; in original form 2022 September 7

## ABSTRACT

Tidal dissipation plays an important role in the dynamical evolution of moons, planets, stars, and compact remnants. The interesting complexity originates from the interplay between the internal structure and external tidal forcing. Recent and upcoming observing missions of exoplanets and stars in the galaxy help to provide constraints on the physics of tidal dissipation. It is timely to develop new *N*-body codes, which allow for experimentation with various tidal models and numerical implementations. We present the open-source *N*-body code TIDYMESS, which stands for 'Tidal Dynamics of Multibody ExtraSolar Systems'. This code implements a Creep deformation law for the bodies, parametrized by their fluid Love numbers and fluid relaxation times. Due to tidal and centrifugal deformations, we approximate the general shape of a body to be an ellipsoid. We calculate the associated gravitational field to quadruple order, from which we derive the gravitational accelerations and torques. The equations of motion for the orbits, spins and deformations are integrated directly using a fourth-order integration method based on a symplectic composition. We implement a novel integration method for the deformations, which allows for a time-step solely dependent on the orbits, and not on the spin periods or fluid relaxation times. This feature greatly speeds up the calculations, while also improving the consistency when comparing different tidal regimes. We demonstrate the capabilities and performance of TIDYMESS, particularly in the niche regime of parameter space where orbits are chaotic and tides become non-linear.

**Key words:** methods: numerical – stars: kinematics and dynamics – stars: rotation – planets and satellites: dynamical evolution and stability.

## 1 INTRODUCTION

### 1.1 Pure *N*-body problem

The *N*-body problem consists of some initial configuration of *N* bodies, represented as point-particles, and calculating their orbits in time by solving Newton's equations of motion and law of gravitation (Newton 1687). This challenging task remains an active field of research, and has resulted in numerous different *N*-body codes in the last few decades. This repertoire of codes consists of variations in the balance between accuracy and speed, long term conservation properties, and regimes of applicability. The pure *N*-body problem as described above, provides a first-order description of the dynamical evolution of planetary and stellar systems. Interesting dynamical phenomena in astronomical systems have already been analysed this way, such as exponential sensitivity and chaos (e.g. Poincaré 1905; Miller 1964; Sussman & Wisdom 1992; Goodman, Heggie & Hut 1993; Boekholt et al. 2016; Stone & Leigh 2019; Boekholt, Portegies Zwart & Valtonen 2020; Portegies Zwart et al. 2022), orbital resonances in planetary systems (e.g. Laplace 1784; Peale 1976; Wisdom 1982; Tsiganis et al. 2005; Laskar & Gastineau 2009; Rein et al. 2012), Lidov–Kozai cycles in triple systems (e.g.

Lidov 1961; Kozai 1962; Ford, Kozinsky & Rasio 2000; Perets & Kratter 2012; Naoz 2016; Toonen et al. 2020), core collapse in star clusters (e.g. McMillan, Hut & Makino 1990; Makino 1996; Baumgardt, Hut & Heggie 2002; Tanikawa, Hut & Makino 2012; Fujii & Portegies Zwart 2014), and much more.

Detailed observations of the moons and planets in the Solar System however, as well as recent state of the art observations of exoplanets and gravitational waves, have presented some limitations of the pure *N*-body problem. The precession of Mercury's orbit by 43 arcsec per century has a relativistic origin (Einstein 1915), and can be reproduced by adding Post-Newtonian terms to the accelerations (e.g. Will 2014). The observation that the Moon is receding from Earth by about 4 cm yr⁻<fnref n="1" /> (Dickey et al. 1994) motivated the study and addition of tides to the *N*-body problem. The population of hot Jupiters with orbital periods of a few days is hard to explain using conservative *N*-body simulations, and instead, is thought to result from planet–planet scattering or Lidov–Kozai cycles combined with tidal dissipation (e.g. Wu & Murray 2003; Fabrycky & Tremaine 2007; Correia et al. 2011; Beaugé & Nesvorný 2012; Bataille, Libert & Correia 2018). Alternatively, disc tidal interactions may drive a planet to migrate inwards (e.g. Lin, Bodenheimer & Richardson 1996; Murray et al. 1998). Another example is the observed rate of gravitational wave sources, which is higher than expected from pure *N*-body simulations of stellar systems (e.g. Samsing & D'Orazio 2018; Kocsis 2022). Also here, the inclusion of gravitational tides,

<footnote n="1">¹Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, University of Oxford, Parks Road, Oxford, OX1 3PU, UK</footnote>

<footnote n="2">²CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal</footnote>

<footnote n="3">³IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 Avenue Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote>*E-mail: tjarda.boekholt@physics.ox.ac.uk (TCNB); acor@uc.pt (ACMC)</footnote>
in combination with gravitational wave dissipation, is thought to be a fruitful line of inquiry (e.g. Bildsten & Cutler 1992; Mandel & de Mink 2016).

## 1.2 N-body problem with Creep tides

For increasingly dense multibody systems the point-particle approximation will eventually start to break down. Structural properties of the bodies themselves have to be taken into account. These typically include the radius of the body, the internal density distribution (uniform or centrally concentrated), the spin rate and orientation, and the tidal response, that is, how are the bodies tidally deformed, and how does that affect the orbits and spins? The physics of tidal dissipation is complex due to the interplay between internal structure processes, and external tidal forcing due to other bodies. The N-body environment is dynamic and often chaotic, thus making it challenging to calculate tidal effects and to predict the long-term tidal evolution.

There are numerous tidal models in the literature, motivated by the wide variety of applications related to tides. The tidal response depends on the type of body, for example, small rocky bodies, giant planets, low-mass stars, giant stars, compact remnants, etc. However, models are also distinguished by whether one is interested in the effects of tides on the detailed internal structure, or rather on the orbital evolution, which affects the structure of the N-body system as a whole. Numerically, the most accurate method is to resolve each body by its own ensemble of bodies. These ‘sub-bodies’ could represent hydrodynamical gas particles, or pebbles part of a larger asteroid (e.g. Quillen et al. 2016; Wang, Mora & Liang 2022). The net deformation is then calculated directly from the sum of their mutual forces, together with the external force from the other bodies. Alternatively, the internal structure can be approximated by a core and an envelope, with an associated spectrum of natural oscillation modes (e.g. Ogilvie 2009). Particularly interesting is the excitation of discrete oscillation frequencies by a second orbiting body. This might leave an imprint in the gravitational wave signal of double white dwarfs, as well as in the flux of heartbeat stars.

From the point of view of modelling N-body systems with $N \gg 3$ and/or long-term integrations, it is preferable to adopt fast and parametrized tidal models, which on average produce the correct dynamical evolution. A commonly adopted model is the linear tide model (e.g. Mignard 1979; Hut 1981; Mardling & Lin 2002; Correia 2018). Here, the bodies are assumed to take on their tidally deformed equilibrium shapes, but with a very small tidal time lag. The current-day satellite systems on near-circular orbits are an example in which linear tides are sufficient to describe their short-term evolution (e.g. Ćuk, El Moutamid & Tiscareno 2020). On the other hand, for highly eccentric orbits, the linear approximation might not hold during pericentre passages, that is, the tidal time lag might exceed the orbital and spin time-scales. In this inverse regime, a non-linear or fast tide model is required (e.g. Correia et al. 2014). Transitions from one tidal regime to another are expected to occur frequently in compact and chaotic multibody systems. This motivates a tidal model which self-consistently handles the complete spectrum from linear to non-linear tides.

A more realistic description of the tidal response is to assume that the body always tends to creep towards the equilibrium by the action of the gravitational forces acting on it. The Creep tide theory, or simply the Creep model, describes the response to tidal stress by supposing a linear Creep equation for the shape evolution of the body (Ferraz-Mello 2013). For low-frequency perturbations, the Creep model is consistent with the linear tide model. However, the Creep model exhibits a peak dissipation when the tidal forcing period

resonates with the bodies' fluid relaxation time. Further into the non-linear regime, tidal dissipation reduces again. Although there is a characteristic tidal frequency set by the fluid relaxation time, we note that this is different than the discrete oscillation frequencies considered by dynamic tide models.

The Creep model also exhibits the same tidal dissipation as the visco-elastic Maxwell model, which is particularly well suited for reproducing the long-term deformation of the planets (e.g. Correia et al. 2014). Therefore, it can be used to model the tidal evolution of these bodies more accurately than the linear model. For instance, it can handle capture in asynchronous spin-orbit resonances naturally (e.g. Valente & Correia 2022). For stars and gaseous planets, the fluid relaxation time is expected to be several orders of magnitude smaller than for rocky planets (see Ferraz-Mello 2013, Table 1), and so the Creep model usually remains in the linear regime. However, during the early stages of the formation of hot-Jupiters and close binary stars, the orbits are expected to be very eccentric. As a result, the short-period passage at the pericentre, where tides are stronger, may require the use of the Creep model. Furthermore, despite dynamic tides being a more accurate description of the dissipation within this class of bodies (see Ogilvie 2009, Fig. 2), the envelope of the multiple frequency perturbations, that is, the averaged trend in the dissipation spectra, is consistent with a Creep model.

TIDYMESS is designed for long-term studies, where the orbits can present chaotic evolution and change from very eccentric to circular. Therefore, the Creep model is better suited than the linear model, as it can handle different dissipation regimes. Although the Creep model is still a simple tidal model characterized by a single parameter, the relaxation time $\tau$, for long-term studies we do not focus on subtle details such as a given dissipation peak at a specific forcing frequency. On the contrary, we mostly care about the overall picture. For short and long-term perturbations, we know that the dissipation has to vanish, but at some point in the middle there is a maximum dissipation, which can be calibrated by means of $\tau$.

## 1.3 Outline and novelty

We present the new N-body code TIDYMESS,<fnref n="1" /> which stands for the 'Tidal Dynamics of Multibody ExtraSolar Systems'. This code covers a niche regime of parameter space where orbits are chaotic and tides can become non-linear. It handles the tidal deformation using a methodology similar to the ones proposed by Correia et al. (2014); Boué, Correia & Laskar (2016); Ragazzo & Ruiz (2017).

In TIDYMESS we implement the Creep tide model, with a fourth-order direct integration method. The interplay between orbital, spin and shape evolution is calculated self-consistently. Arbitrary number of bodies and configurations can be simulated within a Cartesian inertial frame, including full information on the spin vectors and inertia tensors of the bodies. The main numerical novelty of our method is the ability to integrate the deformation over an orbital time-step, irrespective of its fluid relaxation time or spin period (see Section 3.2 for the theoretical derivation and Section 4.1 for a numerical demonstration). Our code complements the repertoire of tidal N-body codes in the literature (e.g. Bolmont et al. 2015; Tamayo et al. 2020; Chen, Li & Tao 2021; Trani & Spera 2023; Lu et al. 2023) with a more realistic tidal model. It facilitates comparison and scrutiny of tidal models, which can be tested against current and

<footnote n="1">¹TIDYMESS is open-source and available through GitHub (https://github.com/tidymess-code) and through the Astrophysical Multipurpose Software Environment (https://github.com/amusecode).</footnote>
upcoming observations of compact, multibody systems. TIDYMESS is also incorporated in the AMUSE framework, where it can be coupled to other codes for stellar evolution and hydrodynamics (Portegies Zwart & McMillan 2018).

We present the theoretical framework of the tidal model in Section 2. The numerical implementation is described in Section 3. A gallery of demonstrations of TIDYMESS is presented in Section 4, including planetary and stellar applications. The main take away messages are gathered in Section 5.

## 2 MODEL

In this section, we derive a very general model that is suited to study a system of *N* bodies. Our model is valid in 3D for the orbital planes and individual spins. We use Cartesian inertial coordinates.

### 2.1 Tidal model

#### 2.1.1 Gravitational potential of an extended body

We consider an extended body of mass *m*, and choose as reference the Cartesian inertial frame (*i*, *j*, *k*). In this frame, the angular velocity and rotational angular momentum vectors of the body are given by $\mathbf{\Omega} = (\Omega_1, \Omega_2, \Omega_3)$ and $\mathbf{L} = (L_1, L_2, L_3)$, respectively, which are related through the inertia tensor $\mathcal{I}$ as

$$ \mathbf{L} = \mathcal{I} \cdot \mathbf{\Omega} \Leftrightarrow \mathbf{\Omega} = \mathcal{I}^{-1} \cdot \mathbf{L}, \quad (1) $$

where

$$ \mathcal{I} = \begin{bmatrix} I_{11} & I_{12} & I_{13} \\ I_{12} & I_{22} & I_{23} \\ I_{13} & I_{23} & I_{33} \end{bmatrix}, \quad (2) $$

$$ \mathcal{I}^{-1} = \frac{1}{\Delta\mathcal{I}} \begin{bmatrix} I_{22}I_{33} - I_{23}^2 & I_{13}I_{23} - I_{12}I_{33} & I_{12}I_{23} - I_{22}I_{13} \\ I_{13}I_{23} - I_{12}I_{33} & I_{11}I_{33} - I_{13}^2 & I_{12}I_{13} - I_{11}I_{23} \\ I_{12}I_{23} - I_{22}I_{13} & I_{12}I_{13} - I_{11}I_{23} & I_{11}I_{22} - I_{12}^2 \end{bmatrix}, \quad (3) $$

with

$$ \Delta\mathcal{I} = I_{11}I_{22}I_{33} + 2I_{12}I_{13}I_{23} - I_{11}I_{23}^2 - I_{22}I_{13}^2 - I_{33}I_{12}^2. \quad (4) $$

The gravitational potential of the body at a generic position *r* from its centre-of-mass is given by (e.g. Goldstein 1950)

$$ V(\mathbf{r}) = -\frac{Gm}{r} + \frac{3G}{2r^3} \left[ \hat{\mathbf{r}} \cdot \mathcal{I} \cdot \hat{\mathbf{r}} - \frac{1}{3} \mathrm{tr}(\mathcal{I}) \right], \quad (5) $$

where *G* is the gravitational constant, $\hat{\mathbf{r}} = \mathbf{r}/r = (\hat{x}, \hat{y}, \hat{z})$ is the unit vector, and tr($\mathcal{I}$) = $I_{11} + I_{22} + I_{33}$. We neglect terms in $(\mathbf{R}/r)^3$, where *R* is the average radius of the body (quadrupolar approximation). Adopting the Lagrange polynomial $P_2(x) = (3x^2 - 1)/2$, we can rewrite the previous potential as

$$ V(\mathbf{r}) = -\frac{Gm}{r} + \frac{G}{r^3} \left[ (I_{22} - I_{11}) P_2(\hat{y}) + (I_{33} - I_{11}) P_2(\hat{z}) + 3(I_{12}\hat{x}\hat{y} + I_{13}\hat{x}\hat{z} + I_{23}\hat{y}\hat{z}) \right]. \quad (6) $$

#### 2.1.2 Deformation of the extended body

The mass distribution inside an extended body, characterized by its inertia tensor, is a result of the forces acting on it, that is, self-gravity and the body's response to any perturbing potential. Here, we consider that the extended body deforms under the action of the centrifugal and tidal potentials.

For a body with rotation rate $\mathbf{\Omega}$, a mass element $dm$ at a location $\mathbf{r}'$ is subject to a centrifugal potential (e.g. Goldstein 1950)

$$ V_c(\mathbf{r}') = -\frac{1}{2}(\mathbf{\Omega} \times \mathbf{r}')^2 = -\frac{1}{2}\mathbf{\Omega}^2 r'^2 + \frac{1}{2}(\mathbf{\Omega} \cdot \mathbf{r}')^2. \quad (7) $$

Assuming also that the body moves in the field of a point-mass $m_0$ located at $\mathbf{r}$, a mass element $dm$ at a location $\mathbf{r}'$ is subject to a tidal potential (e.g. Lambeck 1980)

$$ V_t(\mathbf{r}') = -\frac{Gm_0}{r} \left(\frac{\mathbf{r}'}{r}\right)^2 P_2(\hat{\mathbf{r}}' \cdot \hat{\mathbf{r}}), \quad (8) $$

where we have neglected again terms higher than $(\mathbf{r}'/r)^3$.

The resulting perturbing potential is then given by

$$ V_p(\mathbf{r}') = V_c(\mathbf{r}') + V_t(\mathbf{r}'), \quad (9) $$

that can also be rearranged as

$$ V_p(\mathbf{r}') = -\frac{r'^2\mathbf{\Omega}^2}{3} + \frac{3G}{2r'^3} \left[ \hat{\mathbf{r}}' \cdot \mathcal{I}_p' \cdot \hat{\mathbf{r}}' - \frac{1}{3} \mathrm{tr}(\mathcal{I}_p') \right], \quad (10) $$

with

$$ \frac{\mathcal{I}_p'}{mr'^2} = \frac{r'^3}{3Gm} \mathbf{\Omega} \cdot \mathbf{\Omega}^T - \frac{m_0}{m} \left(\frac{\mathbf{r}'}{r}\right)^3 \hat{\mathbf{r}}' \hat{\mathbf{r}}^T, \quad (11) $$

where $^T$ denotes the transpose. Alternatively, we can write

$$ V_p(\mathbf{r}') = -\frac{r'^2\mathbf{\Omega}^2}{3} + \frac{3G}{2r'^3} \left[ \hat{\mathbf{r}}' \cdot \mathcal{I}_p \cdot \hat{\mathbf{r}}' \right], \quad (12) $$

with

$$ \mathcal{I}_p = \mathcal{I}_p' - \frac{1}{3} \mathrm{tr}(\mathcal{I}_p') \mathbb{I}, \quad (13) $$

where $\mathbb{I}$ is the identity tensor. Thus, on the body's surface, $\mathbf{R}$, the non-spherical contribution of the perturbing potential is given by

$$ V_p(\mathbf{R}) = \frac{3G}{2R^3} \left[ \hat{\mathbf{R}} \cdot \mathcal{I}_p \cdot \hat{\mathbf{R}} \right], \quad (14) $$

where $\mathbf{R} = |\mathbf{R}|$. For simplicity, when the surface of the body is nearly spherical, we can assume $\mathbf{R}$ to be constant and equal to the average radius of the body.

A most convenient way of handling the deformation is through the Love number approach (e.g. Love 1927; Peltier 1974). As long as the distortions are small, we can simplify the problem by ignoring the small interaction terms between the centrifugal and tidal potentials (Zharkov & Trubitsyn 1978). The equilibrium potential $V_e$ is then obtained for static perturbations or instantaneous response simply as

$$ V_e(\mathbf{R}) = k_f V_p(\mathbf{R}) = \frac{3G}{2R^3} (\hat{\mathbf{R}} \cdot \mathcal{I}_e \cdot \hat{\mathbf{R}}), \quad (15) $$

with (equation 13)

$$ \frac{\mathcal{I}_e}{mR^2} = \frac{k_f R^3}{3Gm} \left( \mathbf{\Omega} \cdot \mathbf{\Omega}^T - \frac{\mathbf{\Omega}^2}{3} \right) - k_f \frac{m_0}{m} \left( \frac{\mathbf{R}}{r} \right)^3 \left( \hat{\mathbf{r}} \cdot \hat{\mathbf{r}}^T - \frac{1}{3} \mathbb{I} \right), \quad (16) $$

where $k_f$ is the fluid second Love number for potential. By comparing expression (15) with the gravitational potential in expression (5) at the body's surface, we see that $\mathcal{I}_e$ corresponds to the non-spherical contribution to the inertia tensor. The values for each coefficient of $\mathcal{I}_e$ are

$$ \frac{I_{11}^e}{mR^2} = \frac{k_f R^3}{3Gm} \left( \Omega_1^2 - \frac{\Omega_2^2}{3} \right) - k_f \frac{m_0}{m} \left( \frac{R}{r} \right)^3 \left( \hat{s}^2 - \frac{1}{3} \right), \quad (17) $$

$$ \frac{I_{22}^e}{mR^2} = \frac{k_f R^3}{3Gm} \left( \Omega_2^2 - \frac{\Omega_3^2}{3} \right) - k_f \frac{m_0}{m} \left( \frac{R}{r} \right)^3 \left( \hat{s}^2 - \frac{1}{3} \right), \quad (18) $$
$$ \frac{I_{33}^{e}}{m R^{2}} = \frac{k_{f} R^{3}}{3Gm} \left( \Omega_{3}^{2} - \frac{\Omega^{2}}{3} \right) - k_{f} \frac{m_{0}}{m} \left( \frac{R}{r} \right)^{3} \left( \hat{z}^{2} - \frac{1}{3} \right), \quad (19) $$

$$ \frac{I_{12}^{e}}{m R^{2}} = \frac{k_{f} R^{3}}{3Gm} \Omega_{1} \Omega_{2} - k_{f} \frac{m_{0}}{m} \left(\frac{R}{r}\right)^{3} \hat{x} \hat{y}, \quad (20) $$

$$ \frac{I_{13}^{e}}{m R^{2}} = \frac{k_{f} R^{3}}{3Gm} \Omega_{1} \Omega_{3} - k_{f} \frac{m_{0}}{m} \left(\frac{R}{r}\right)^{3} \hat{x} \hat{z}, \quad (21) $$

$$ \frac{I_{23}^{e}}{m R^{2}} = \frac{k_{f} R^{3}}{3Gm} \Omega_{2} \Omega_{3} - k_{f} \frac{m_{0}}{m} \left(\frac{R}{r}\right)^{3} \hat{y} \hat{z}. \quad (22) $$

In general, the perturbing potential is not constant and the deformation of the body is not instantaneous. However, for the Creep model, in the body frame B, the inertia tensor of the body $\mathcal{I}^B$, and the instantaneous equilibrium inertia tensor $\mathcal{I}_e^B$, can be related through a simple differential equation as (Correia et al. 2014):

$$ \mathcal{I}^B + \tau \dot{\mathcal{I}}^B = \mathcal{I}_0 + \mathcal{I}_e^B, \quad (23) $$

where $\tau = 19\eta/(2g\rho R)$ is the fluid relaxation time, $\eta$ is the viscosity, $g$ is the surface gravity and $\rho$ is the mean density. $\mathcal{I}_0 = \xi m R^2$ corresponds to the spherical contribution to the inertia tensor of the body, where $\xi$ is the inner structure factor (we have $\xi = 2/5$ for an homogenous body).

We let $\mathcal{S}$ be the rotation matrix that allows us to convert any vector $\mathbf{u}_B$ in a frame attached to the body into the Cartesian inertial frame $\mathbf{u}$, such that $\mathbf{u} = \mathcal{S} \mathbf{u}_B$. The time evolution of $\mathcal{S}$ is given by

$$ \dot{\mathcal{S}} = \tilde{\mathcal{S}} \mathcal{S}, \quad \text{and} \quad \dot{\mathcal{S}}^T = -\mathcal{S}^T \tilde{\mathcal{S}}, \quad (24) $$

with

$$ \tilde{\mathcal{S}} = \begin{bmatrix} 0 & -\Omega_3 & \Omega_2 \\ \Omega_3 & 0 & -\Omega_1 \\ -\Omega_2 & \Omega_1 & 0 \end{bmatrix}. \quad (25) $$

In the inertial frame,

$$ \mathcal{I} = S \mathcal{I}^B S^T, \quad (26) $$

and so we can rewrite expression (23) as

$$ \mathcal{I} + \tau S \dot{\mathcal{I}}^B S^T = \mathcal{I}_0 + \mathcal{I}_e. \quad (27) $$

Now, differentiating equation (26), and using the derivatives of $\mathcal{S}$ given by equation (24), we have

$$ \dot{\mathcal{I}} = \tilde{\mathcal{S}} \mathcal{I} - \mathcal{I} \tilde{\mathcal{S}} + S \dot{\mathcal{I}}^B S^T, \quad (28) $$

or, using equation (27)

$$ \dot{\mathcal{I}} = \tilde{\mathcal{S}} \mathcal{I} - \mathcal{I} \tilde{\mathcal{S}} + (\mathcal{I}_0 + \mathcal{I}_e - \mathcal{I})/\tau. \quad (29) $$

For incompressible bodies, $\text{tr}(\mathcal{I}) = \text{tr}(\mathcal{I}_0) = 3\xi m R^2$. We thus have, $I_{33} = 3\xi m R^2 - (I_{11} + I_{22})$, that is, the inertia tensor has only five independent coefficients:

$$ \left\{ \begin{array}{l} \dot{I}_{11} = 2I_{13}\Omega_2 - 2I_{12}\Omega_3 + (\xi m R^2 + I_{11}^e - I_{11})/\tau \\ \dot{I}_{22} = -2I_{23}\Omega_1 + 2I_{12}\Omega_3 + (\xi m R^2 + I_{22}^e - I_{22})/\tau \\ \dot{I}_{12} = -I_{13}\Omega_1 + I_{23}\Omega_2 + (I_{11} - I_{22})\Omega_3 + (I_{12}^e - I_{12})/\tau \\ \dot{I}_{13} = I_{12}\Omega_1 - (2I_{11} + I_{22})\Omega_2 - I_{23}\Omega_3 + (I_{13}^e - I_{13})/\tau \\ \dot{I}_{23} = (I_{11} + 2I_{22})\Omega_1 - I_{12}\Omega_2 + I_{13}\Omega_3 + (I_{23}^e - I_{23})/\tau \end{array} \right. . \quad (30) $$

### 2.1.3 Linear approximation

When the relaxation time, $\tau$, is much smaller than the orbital and spin evolution time-scales, the integration time of the system is dominated by the deformation equations (29). However, for small $\tau$ values, we

can assume that the deformations are linear and directly obtain an expression for $\mathcal{I}$. From expression (23) we compute

$$ \mathcal{I}^B(t) = \mathcal{I}_0 + \frac{e^{-t/\tau}}{\tau} \int_0^t \mathcal{I}_e^B(x) e^{x/\tau} dx + C_1 e^{-t/\tau}, \quad (31) $$

where $C_k$ is an integration constant. Integrating by parts we get

$$ \int_0^t \mathcal{I}_e^B(x) e^{x/\tau} dx = C_2 + \tau \mathcal{I}_e^B(t) e^{t/\tau} - \tau^2 \dot{\mathcal{I}}_e^B(t) e^{t/\tau} + \mathcal{O}(\tau^3), \quad (32) $$

and thus,

$$ \mathcal{I}^B(t) = \mathcal{I}_0 + \mathcal{I}_e^B(t) - \tau \dot{\mathcal{I}}_e^B(t) + C_0 e^{-t/\tau} + \mathcal{O}(\tau^2). \quad (33) $$

In the following, we neglect the term $C_0 e^{-t/\tau}$ because it vanishes after some time. Replacing in expression (26) we get

$$ \mathcal{I} = \mathcal{I}_0 + \mathcal{I}_e - \tau S \dot{\mathcal{I}}_e^B S^T + \mathcal{O}(\tau^2), \quad (34) $$

or, similarly (equation 28)

$$ \mathcal{I} = \mathcal{I}_0 + \mathcal{I}_e + \tau (\tilde{\mathcal{S}} \mathcal{I}_e - \mathcal{I}_e \tilde{\mathcal{S}} - \dot{\mathcal{I}}_e) + \mathcal{O}(\tau^2). \quad (35) $$

## 2.2 N-body problems

### 2.2.1 Point-mass problem

We now consider that the extended body orbits a point-mass $m_0$ located at $\mathbf{r}$. The force between the two bodies is easily obtained from the potential energy of the system $U(\mathbf{r}) = m_0 V(\mathbf{r})$ (equation 6) as

$$ F = -\nabla U(\mathbf{r}) = f(m_0, m, \mathbf{r}) + g(m_0, \mathcal{I}, \mathbf{r}) + h(m_0, \mathcal{I}, \mathbf{r}), \quad (36) $$

with

$$ f(m_0, m, \mathbf{r}) = -\frac{Gm_0m}{r^3} \mathbf{r}, \quad (37) $$

$$ g(m_0, \mathcal{I}, \mathbf{r}) = \frac{15Gm_0}{r^5} \left[ \frac{I_{22} - I_{11}}{2} (\hat{y}^2 - \frac{1}{5}) + \frac{I_{33} - I_{11}}{2} (\hat{z}^2 - \frac{1}{5}) + I_{12}\hat{x}\hat{y} + I_{13}\hat{x}\hat{z} + I_{23}\hat{y}\hat{z} \right] \mathbf{r}, \quad (38) $$

$$ h(m_0, \mathcal{I}, \mathbf{r}) = -\frac{3Gm_0}{r^4} \left[ (I_{22} - I_{11})\hat{y}\mathbf{j} + (I_{33} - I_{11})\hat{z}\mathbf{k} + I_{12}(\hat{x}\mathbf{j} + \hat{y}\mathbf{i}) + I_{13}(\hat{x}\mathbf{k} + \hat{z}\mathbf{i}) + I_{23}(\hat{y}\mathbf{k} + \hat{z}\mathbf{j}) \right]. \quad (39) $$

We thus obtain for the orbital evolution of the system

$$ \ddot{\mathbf{r}} = \mathbf{F}/\beta, \quad (40) $$

where $\beta = m_0 m/(m_0 + m)$ is the reduced mass. The spin evolution of the extended body can also be obtained from the force, by computing the gravitational torque. In the inertial frame we have

$$ \dot{\mathbf{L}} = \mathbf{T}(m_0, \mathcal{I}, \mathbf{r}) = -\mathbf{r} \times \mathbf{F} = -\mathbf{r} \times \mathbf{h}, \quad (41) $$

that is,

$$ \mathbf{T}(m_0, \mathcal{I}, \mathbf{r}) = \frac{3Gm_0}{r^3} \hat{\mathbf{r}} \times \left[ (I_{22} - I_{11})\hat{y}\mathbf{j} + (I_{33} - I_{11})\hat{z}\mathbf{k} + I_{12}(\hat{x}\mathbf{j} + \hat{y}\mathbf{i}) + I_{13}(\hat{x}\mathbf{k} + \hat{z}\mathbf{i}) + I_{23}(\hat{y}\mathbf{k} + \hat{z}\mathbf{j}) \right], \quad (42) $$
or

$$ T = \frac{3Gm_0}{r^3} \left[ \begin{pmatrix} I_{33} - I_{22} & \hat{y}\hat{z} - I_{12}\hat{x}\hat{z} + I_{13}\hat{x}\hat{y} + I_{23}(\hat{y}^2 - \hat{z}^2) \\ I_{11} - I_{33} & \hat{x}\hat{z} + I_{12}\hat{y}\hat{z} + I_{13}(\hat{z}^2 - \hat{x}^2) - I_{23}\hat{x}\hat{y} \\ I_{22} - I_{11} & \hat{x}\hat{y} + I_{12}(\hat{x}^2 - \hat{y}^2) - I_{13}\hat{y}\hat{z} + I_{23}\hat{x}\hat{z} \end{pmatrix} \right] . \quad (43) $$

To solve the spin-orbit motion, we need to integrate equations (29), (40), and (41).

### 2.2.2 Two-body problem

Consider now that two extended bodies with masses $m_0$ and $m_1$, and inertia tensors $\mathcal{I}_0$ and $\mathcal{I}_1$, respectively, orbit around each other at a distance $\mathbf{r}$ from their centres-of-mass. The total potential energy can be written from expression (5) as

$$ U = -\frac{Gm_0m_1}{r} + \frac{3G}{2r^3} \left[ \hat{\mathbf{r}} \cdot \mathcal{J} \cdot \hat{\mathbf{r}} - \frac{1}{3} \mathrm{tr}(\mathcal{J}) \right], \quad (44) $$

with $\mathcal{J} = m_0\mathcal{I}_0 + m_1\mathcal{I}_0$. This potential is very similar to the previous point-mass problem and the equations of motion are simply

$$ \ddot{\mathbf{r}} = \mathbf{F}_{01}(\mathbf{r})/\beta_{01}, \quad (45) $$

$$ \dot{\mathbf{L}}_0 = \mathbf{T}_{01}(\mathbf{r}), \quad \dot{\mathbf{L}}_1 = \mathbf{T}_{10}(\mathbf{r}), \quad (46) $$

where

$$ \mathbf{F}_{kl}(\mathbf{r}) = f(m_k, m_l, \mathbf{r}) + g(m_k, \mathcal{I}_l, \mathbf{r}) + g(m_l, \mathcal{I}_k, \mathbf{r}) + \mathbf{h}(m_k, \mathcal{I}_l, \mathbf{r}) + \mathbf{h}(m_l, \mathcal{I}_k, \mathbf{r}), \quad (47) $$

$$ \mathbf{T}_{kl}(\mathbf{r}) = \mathbf{T}(m_l, \mathcal{I}_k, \mathbf{r}), \quad (48) $$

$\beta_{kl} = m_k m_l (m_k + m_l)$, $\mathbf{L}_k = \mathcal{I}_k \mathbf{\Omega}_k$ is the rotational angular momentum vector of the body with mass $m_k$, and $\mathbf{\Omega}_k$ its angular velocity.

### 2.2.3 N-body problem

The previous equations can be easily generalized to the motion of several bodies. For a system of $N+1$ extended bodies, with masses $m_k$ and inertia tensors $\mathcal{I}_k$ ($k=0, 1, ..., N$), the total potential energy can be written from expression (44) as

$$ U = \sum_{k=0}^{N} \sum_{l>k}^{N} \left( -\frac{Gm_k m_l}{r_{kl}} + \frac{3G}{2r_{kl}^3} \left[ \hat{\mathbf{r}}_{kl} \cdot \mathcal{J}_{kl} \cdot \hat{\mathbf{r}}_{kl} - \frac{1}{3} \mathrm{tr}(\mathcal{J}_{kl}) \right] \right), \quad (49) $$

with $\mathcal{J}_{kl} = m_k \mathcal{I}_l + m_l \mathcal{I}_k$, and $\mathbf{r}_{kl} = \mathbf{R}_k - \mathbf{R}_l$, where $\mathbf{R}_k$ is the position of the centre-of-mass of the body $k$ in the inertial frame. The equations of motion for each body in the inertial frame are thus

$$ \ddot{\mathbf{R}}_k = \mathbf{a}_k = \frac{1}{m_k} \sum_{l \neq k} \mathbf{F}_{kl}(\mathbf{r}_{kl}), \quad (50) $$

$$ \dot{\mathbf{L}}_k = \sum_{l \neq k} \mathbf{T}_{kl}(\mathbf{r}_{kl}), \quad (51) $$

where $\mathbf{F}_{kl}(\mathbf{r}_{kl})$ and $\mathbf{T}_{kl}(\mathbf{r}_{kl})$ are given by expressions (47) and (48), respectively. As for the two-body problem, we can also express the motion in the relative frame, for instance, with respect to the body with mass $m_0$. We let $\mathbf{r}_k = \mathbf{r}_{k0} = \mathbf{R}_k - \mathbf{R}_0$ ($k=1, ..., N$). Thus,

$$ \ddot{\mathbf{r}}_k = \frac{1}{\beta_{0k}} \mathbf{F}_{0k}(\mathbf{r}_k) + \sum_{l \neq k} \left[ \frac{1}{m_k} \mathbf{F}_{kl}(\mathbf{r}_k - \mathbf{r}_l) + \frac{1}{m_0} \mathbf{F}_{0l}(\mathbf{r}_l) \right]. \quad (52) $$

Note also that for the equilibrium inertia tensor (equation 16), we now need to take into account the tidal perturbations from all bodies

in the system:

$$ \mathcal{I}_{e, k} = \frac{k_f R_k^5}{3G} \left( \mathbf{\Omega}_k \mathbf{\Omega}_k^T - \frac{\mathbf{\Omega}_k^2}{3} \right) - \sum_{l \neq k} k_f m_l R_k^2 \left( \frac{R_k}{r_{kl}} \right)^3 \left( \hat{\mathbf{r}}_{kl} \hat{\mathbf{r}}_{kl}^T - \frac{1}{3} \mathbf{\mathbf{I}} \right), \quad (53) $$

where $k_f$ also pertains to the body with index $k$.

## 2.3 Other physical phenomena

### 2.3.1 General relativity corrections

Relativistic effects become important in close binary systems, through apsidal precession and gravitational wave emission leading to orbital shrinkage. Relativistic perturbations also play a role in large stellar systems by affecting the global chaotic properties, for example, the Lyapunov time-scale (Portegies Zwart et al. 2022). Furthermore, interesting phenomena are expected to result from the interplay between relativity and tides in compact multiple systems.

A common method for implementing relativistic effects in N-body codes is by adding Post-Newtonian (PN) acceleration terms. The force is expanded in a series based on the small ratio of velocity divided by the speed of light. We will include the pairwise, conservative 1PN and 2PN terms, as well as the first dissipative term, 2.5PN. The Post-Newtonian force experienced by body $k$ needed to correct expression (50) is given by

$$ \mathbf{F}_{\mathrm{PN},k} = m_k (\mathbf{a}_{1\mathrm{PN},k} + \mathbf{a}_{2\mathrm{PN},k} + \mathbf{a}_{2.5\mathrm{PN},k}). \quad (54) $$

The 1PN acceleration term is given by

$$ \mathbf{a}_{1\mathrm{PN},k} = \frac{1}{c^2} \sum_{l \neq k} \frac{Gm_l}{r_{kl}^2} (B_{1,kl} \hat{\mathbf{r}}_{kl} + C_{1,kl} \mathbf{v}_{kl}), \quad (55) $$

with $\mathbf{v}_{kl} = \dot{\mathbf{r}}_{kl}$, and $c$ the speed of light. The coefficients are

$$ B_{1,kl} = 5 \frac{Gm_k}{r_{kl}} + 4 \frac{Gm_l}{r_{kl}} - v_k^2 - 2v_l^2 + 4\mathbf{v}_k \cdot \mathbf{v}_l + \frac{3}{2} (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^2, \quad (56) $$

$$ C_{1,kl} = 4\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k - 3\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l. \quad (57) $$

The 2PN acceleration term is given by

$$ \mathbf{a}_{2\mathrm{PN},k} = \frac{1}{c^4} \sum_{l \neq k} \frac{Gm_l}{r_{kl}^2} (B_{2,kl} \hat{\mathbf{r}}_{kl} + C_{2,kl} \mathbf{v}_{kl}) + D_{2,kl} \hat{\mathbf{r}}_{kl}, \quad (58) $$

with

$$ \begin{aligned} B_{2,kl} = & -2v_l^4 + 4v_l^2 \mathbf{v}_k \cdot \mathbf{v}_l - 2(\mathbf{v}_k \cdot \mathbf{v}_l)^2 - \frac{15}{8} (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^4 \\ & + \left( \frac{3}{2} v_k^2 + \frac{9}{2} v_l^2 - 6 \mathbf{v}_k \cdot \mathbf{v}_l \right) (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^2 \\ & + \frac{Gm_l}{r_{kl}} \left( 4v_l^2 - 8 \mathbf{v}_k \cdot \mathbf{v}_l + 2 (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k)^2 \right. \\ & \left. - 4 (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k) (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l) - 6 (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^2 \right) \\ & + \frac{Gm_k}{r_{kl}} \left( -\frac{15}{4} v_k^2 + \frac{5}{4} v_l^2 - \frac{5}{2} \mathbf{v}_k \cdot \mathbf{v}_l + \frac{39}{2} (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k)^2 \right. \\ & \left. - 39 (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k) (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l) + \frac{17}{2} (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^2 \right) \end{aligned} $$

$$ \begin{aligned} C_{2,kl} = & v_k^2 \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l + 4v_l^2 \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k - 5v_l^2 \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l - 4(\mathbf{v}_k \cdot \mathbf{v}_l)(\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k) \\ & +4(\mathbf{v}_k \cdot \mathbf{v}_l)(\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l) - 6(\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k)(\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^2 + \frac{9}{2}(\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l)^3 \\ & + \frac{Gm_k}{r_{kl}} \left( -\frac{63}{4} \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k + \frac{55}{4} \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l \right) \\ & + \frac{Gm_l}{r_{kl}} \left( -2\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k - 2\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l \right) \end{aligned} \quad (59) $$
<table><caption>**Table 1.** Tidal models implemented in TIDYMESS. Each model inherits the particle attributes from the models above: *m* = mass, *r* = position vector, *v* = velocity vector, *R* = mean radius, *ξ* = moment of inertia factor, *I* = inertia tensor, *Ω* = spin frequency vector, *L* = angular momentum vector, *k*<sub>f</sub> = fluid Love number, and *τ* = fluid relaxation time.</caption><thead><tr><th>Tidal model</th><th>Attributes</th><th>Notes</th></tr></thead><tbody><tr><td>0. No tides</td><td><i>m</i>, <i>r</i>, <i>v</i></td><td>Pure <i>N</i>-body code</td></tr><tr><td>1. Conservative tides</td><td><i>R</i>, <i>ξ</i>, <i>I</i>, <i>Ω</i>, <i>L</i>, <i>k</i><sub>f</sub></td><td>Equilibrium shape</td></tr><tr><td>2. Linear tides</td><td><i>τ</i></td><td>Small tidal time lag</td></tr><tr><td>3. Creep tides direct</td><td></td><td>Arbitrary time lag</td></tr><tr><td>4. Creep tides tidymess</td><td></td><td>Orbital time-steps only</td></tr></tbody></table>

$$D_{2,kl} = \frac{G^3 m_l}{r_{kl}^4} \left( -\frac{57}{4} m_k^2 - 9 m_l^2 - \frac{69}{2} m_k m_l \right). \qquad (60)$$

The 2.5PN acceleration term is given by

$$\mathbf{a}_{2.5PN,k} = \frac{1}{c^5} \sum_{l \neq k} \frac{4 G^2 m_k m_l}{5 r_{kl}^3} (B_{25,kl} \hat{\mathbf{r}}_{kl} + C_{25,kl} \mathbf{v}_{kl}), \qquad (61)$$

with

$$B_{25,kl} = \left( 3v_{kl}^2 - 6\frac{Gm_k}{r_{kl}} + \frac{52}{3}\frac{Gm_l}{r_{kl}} \right) (\hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_k - \hat{\mathbf{r}}_{kl} \cdot \mathbf{v}_l), \qquad (62)$$

$$C_{25,kl} = -v_{kl}^2 + 2\frac{Gm_k}{r_{kl}} - 8\frac{Gm_l}{r_{kl}}. \qquad (63)$$

### 2.3.2 Magnetic braking

Low-mass stars lose angular momentum due to spin-down magnetic braking through the action of magnetized stellar winds (Skumanich 1972; Kawaler 1988). Cool stars are thought to be more efficient in generating magnetic winds and thus lose angular momentum more rapidly than hot stars (e.g. Weber & Davis 1967). An empirical relation can be derived from the observations, where $\tilde{\Omega} \propto -\Omega^3$. For Sun-like stars we can use this simple law to compute the contribution of magnetic braking to the angular momentum variation of each star in the model:

$$\dot{\mathbf{L}}_k = -\alpha_k \Omega_k^2 \mathbf{L}_k, \qquad (64)$$

where $\alpha_k$ is a proportionality constant determined from the observations. For G and K-dwarf stars, we estimate $\alpha_k \approx 1.5 \times 10^{-14}$ yr (Barker & Ogilvie 2009). For F-dwarfs the magnetic braking is less efficient and may be only 10 per cent of the previous value. Therefore, in TIDYMESS it is possible to redefine the value of $\alpha_k$ for each star.

### 2.3.3 Mergers and collisions

Mergers and collisions are natural outcomes in compact, multibody systems (e.g. Vergara et al. 2021; Toonen, Boekholt & Portegies Zwart 2022). For example, in resonant multiplanet systems, planet-planet interactions can lead to one planet being scattered towards the host star. There is a thin line between either a collision or tidal capture. In the latter case, the subsequent tidal evolution might lead to the formation of a hot Jupiter (see Section 4.5). We will therefore include stopping conditions for collisions, based on the criterion of overlapping radii.

Optionally, the simulations can be continued by replacing the collisional components by an approximate collision product. However, this requires further details on the internal structure properties and corresponding tidal response parameters as a function of mass and impact properties. TIDYMESS implements a generalization of the MNRAS 522, 2885–2900 (2023)

sticky sphere model, which assumes the conservation of mass, centre of mass, linear momentum, and total angular momentum (orbit and spin). Internal properties are adopted from the most massive collision component, for example, the radius is determined by preserving the average density. Alternatively, the code can be halted when a collision is detected, and a collision product can be generated externally, for example, using the AMUSE framework (Portegies Zwart & McMillan 2018) and specialized codes for collisions (e.g. Gaburov et al. 2018).

## 3 NUMERICAL IMPLEMENTATION

In this section we provide details on the numerical implementation of TIDYMESS, based on the physical model given in the previous section. First, we discuss the tidal models implemented in TIDYMESS, followed by our general and optimized approximation of the Creep model integration method. Then, we describe the integration maps, time-step functions and the user interface.<fnref n="2" />

### 3.1 Tidal models

TIDYMESS includes five tidal models/implementations, which are listed in Table 1. The most general and accurate model is the tidal Creep model with direct integration of the deformation equations (equation 29). This model is best used in the regime where the orbital time-scales are of the same order or shorter than the spin and fluid relaxation time-scales. For fast spinning bodies or very short fluid relaxation times, this method becomes excessively slow (see example in Section 4.1). This is because the integration step is required to be the minimum of all orbital periods, spin periods and fluid relaxation times. Therefore, we also implement a simplified Creep tidal model based on the algorithm described in Section 3.2, which we name tidymess. Using this model, the performance in the linear and non-linear tidal regimes is the same. A (small) approximation is introduced however, by the splitting of the tidal and rotational deformation parts (see example in Section 4.1). The third tidal model is the linear tide model (equation 35). If the fluid relaxation time is always much shorter than the orbital and spin time-scales in the system, then the Creep model reduces to the linear tide model. In this regime however, the linear tide model has a better performance. In the asymptotic limit of $\tau \to 0$, we obtain the conservative tide model, which assumes the instantaneous equilibrium shape for the deformations (equation 53). Both angular momentum and energy are conserved in this model, and therefore also allows for backward integration. For completeness, we also implement a model for perfectly rigid spheres, which is the regime of $k_f \to 0$, that is, the pure N-body point-mass problem.

### 3.2 Creep tides tidymess

Analogous to the linear approximation (Section 2.1.3), we start with the general solution to equation (23):

$$\mathcal{I}^B(\Delta t) = \mathcal{I}_0 + (\mathcal{I}^B - \mathcal{I}_0) e^{-\Delta t/\tau} + \frac{e^{-\Delta t/\tau}}{\tau} \int_0^{\Delta t} \mathcal{I}_e^B(\alpha) e^{\alpha/\tau} dx. \qquad (65)$$

Here, $\Delta t$ is a time-step into the future, and for the integral we set the current time to $t=0$. Since $\mathcal{I}_0$ is constant, we define for simplicity $\mathcal{I}_v$ as the non-spherical part of the inertia tensor

$$\mathcal{I}_v = \mathcal{I} - \mathcal{I}_0 \quad \Rightarrow \quad \mathcal{I} = \mathcal{I}_0 + \mathcal{I}_v. \qquad (66)$$

<footnote n="2">²Further code tutorials and examples are provided with the source code.</footnote>
Using integration by parts, we can write

$$ \mathcal{I}_v^B (\Delta t) = \mathcal{I}_e^B (\Delta t) - (\mathcal{I}_e^B - \mathcal{I}_v^B) e^{-\Delta t/\tau} - e^{-\Delta t/\tau} \int_0^{\Delta t} \dot{\mathcal{I}}_e^B (x) e^{x/\tau} dx. \quad (67) $$

We perform another integration by parts, so that we can write

$$ \mathcal{I}_v^B (\Delta t) = \mathcal{I}_e^B (\Delta t) - \tau \dot{\mathcal{I}}_e^B (\Delta t) - (\mathcal{I}_e^B - \tau \dot{\mathcal{I}}_e^B - \mathcal{I}_v^B) e^{-\Delta t/\tau} \\ + \tau e^{-\Delta t/\tau} \int_0^{\Delta t} \ddot{\mathcal{I}}_e^B (x) e^{x/\tau} dx. \quad (68) $$

We could continue integrating by parts and add terms of increasing order in $\tau$. However, when $\tau \ll \Delta t$, the higher order terms will be very small, and we can therefore ignore terms of order $\tau^2$ and higher (linear approximation). If $\tau \gg \Delta t$ then we cannot ignore the higher order terms. But in that case, $\Delta t$ can be considered very small, which allows us to estimate the final integral in equation (68) as

$$ \tau e^{-\Delta t/\tau} \int_0^{\Delta t} \dot{\mathcal{I}}_e^B (x) e^{x/\tau} dx \approx \tau e^{-\Delta t/\tau} \dot{\mathcal{I}}_e^B \Delta t \approx \tau e^{-\Delta t/\tau} (\dot{\mathcal{I}}_e^B (\Delta t) - \dot{\mathcal{I}}_e^B), \quad (69) $$

where we have estimated that $\dot{\mathcal{I}}_e^B \approx (\dot{\mathcal{I}}_e^B (\Delta t) - \dot{\mathcal{I}}_e^B) / \Delta t$. The full expression for the inertia tensor in the body frame then becomes

$$ \mathcal{I}_v^B (\Delta t) = \mathcal{I}_e^B (\Delta t) + (\mathcal{I}_v^B - \mathcal{I}_e^B) e^{-\Delta t/\tau} - \tau \dot{\mathcal{I}}_e^B (\Delta t) (1 - e^{-\Delta t/\tau}). \quad (70) $$

The derivative of the inertia tensor is further approximated as $\dot{\mathcal{I}}_e^B (\Delta t) \approx (\mathcal{I}_e^B (\Delta t) - \mathcal{I}_e^B) / \Delta t$, resulting in

$$ \mathcal{I}_v^B (\Delta t) = \mathcal{I}_e^B (\Delta t) + (\mathcal{I}_v^B - \mathcal{I}_e^B) e^{-\Delta t/\tau} \\ - \frac{\tau}{\Delta t} (\mathcal{I}_e^B (\Delta t) - \mathcal{I}_e^B) (1 - e^{-\Delta t/\tau}). \quad (71) $$

This equation is general and a good approximation for every value of the fraction $\Delta t/\tau$. To demonstrate this, we give the simplifying expressions in the different asymptotic regimes:

(i) $\mathcal{I}_v^B (\Delta t) = \mathcal{I}_v^B + (\mathcal{I}_e^B - \mathcal{I}_v^B) \Delta t/\tau$, if $\Delta t/\tau \ll 1$,

(ii) $\mathcal{I}_v^B (\Delta t) = \mathcal{I}_e^B (\Delta t) - \tau \dot{\mathcal{I}}_e^B (\Delta t)$, if $\Delta t/\tau \gg 1$.

The expression above is for the inertia tensor in the body frame, that is, it is rotating with the body. In order to get the inertia tensor in the inertial frame, we need to perform a rotation from the body to the inertial frame (see equation 26). Furthermore, the change of the inertia tensor occurs on two separate time-scales: the spin time-scale which is proportional to $\Omega^{-1}$, and the fluid relaxation time-scale, $\tau$ (see equation 29). Since these time-scales can be very different, it is natural to split up the differential equation of the inertia tensor into two parts, in a manner that is often done in Hamiltonian splittings:

$$ \dot{\mathcal{I}} = (\mathcal{I}_e - \mathcal{I}_v)/\tau \quad (\text{Deformation part}). \quad (72) $$

$$ \dot{\mathcal{I}} = \Omega \mathcal{I} - \Omega \dot{\mathcal{I}} \quad (\text{Rotational part}), \quad (73) $$

Each of these parts separately can be solved accurately and without limits on the time-step size. The solution to the deformation part is given in equation (71) without the body frame superscripts. The rotational part can be problematic if the body is spinning very rapidly, for example, when $\Omega \Delta t > 1$ or $\Omega \tau > 1$, then we cannot simply do $\mathcal{I}(\Delta t) = \mathcal{I} + \dot{\mathcal{I}} \Delta t$. This can be solved accurately using a (quaternion) rotation operator. However, the deformation of a fluid body does not rotate in the same way as a rigid, aspherical deformation. In the case of a rigid body we would need to rotate for the full time-step $\Delta t$ around the spin vector axis given by $\Omega$. In the fluid body case, deformations damp and fade away on the time-scale $\tau$. Therefore, when $\Delta t/\tau \ll 1$, the body rotation has to be over a time-step $\Delta t$, but when $\Delta t/\tau \gg 1$, then the rotation has to be over a

time-scale given by $\tau$. This constraint is secured by introducing the following time-step for the rotational part:

$$ \Delta t_{\text{rot}} = \tau (1 - e^{-\Delta t/\tau}). \quad (74) $$

Given the spin vector, $\Omega$, and the rotational time-step, $\Delta t_{\text{rot}}$, we calculate the associated quaternion as

$$ \mathbf{q} = (q_0, q_1, q_2, q_3) = \left[ \cos\left(\frac{\Omega \Delta t_{\text{rot}}}{2}\right), \sin\left(\frac{\Omega \Delta t_{\text{rot}}}{2}\right) \frac{\Omega}{\Omega} \right]. \quad (75) $$

This quaternion is subsequently used to perform a rotation of the inertia tensor, by first calculating the rotation matrix

$$ S(\mathbf{q}) = \begin{bmatrix} 2(q_0^2 + q_1^2) - 1 & 2(q_1q_2 - q_0q_3) & 2(q_1q_3 + q_0q_2) \\ 2(q_1q_2 + q_0q_3) & 2(q_0^2 + q_2^2) - 1 & 2(q_2q_3 - q_0q_1) \\ 2(q_1q_3 - q_0q_2) & 2(q_2q_3 + q_0q_1) & 2(q_0^2 + q_3^2) - 1 \end{bmatrix}, \quad (76) $$

and then performing the actual rotation according to equation (26). Note that we are not integrating quaternions, as for the deformation of fluid bodies it is unnecessary to keep track of the precise orientation of the body.

### 3.3 Integration maps

The conservation of angular momentum plays a key role in simulations of multibody systems with tides. Tidal effects are usually small and therefore the angular momentum exchange between orbits and spins has to be calculated very accurately. The numerical error in the total angular momentum has to be kept to a minimum for the duration of the simulation. The time-symmetric, second-order Verlet-Leapfrog integrator exhibits excellent angular momentum conservation by design. The relative error in the angular momentum is kept close to machine-precision, and due to the symplectic nature of the algorithm, there is no drift in time. The base integrator of TIDYMESS is therefore the time-symmetric, second-order Verlet-Leapfrog (Drift-Kick-Drift) integrator:

$$ \begin{align} \mathbf{r}_{\frac{1}{2}} &= \mathbf{r}_0 + \mathbf{v}_0 \frac{h}{2} \quad (\text{Drift}) \notag \\ \mathbf{a}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) &= \mathbf{a}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) \quad (\text{Acceleration}) \notag \\ \mathbf{v}_{1} &= \mathbf{v}_{0} + \mathbf{a}_{\frac{1}{2}} \frac{h}{2} \quad (\text{Kick orbit}) \notag \\ \mathbf{r}_{1} &= \mathbf{r}_{\frac{1}{2}} + \mathbf{v}_{1} \frac{h}{2} \quad (\text{Drift}) \tag{77} \end{align} $$

A fourth-order integrator is constructed by a composition of Leapfrog steps (e.g. Suzuki 1990; Yoshida 1993; McLachlan 1995). We adopt a method by McLachlan (1995), where a single 'McLachlan step' of size $h$, is composed of 5 Leapfrog steps with sizes: $c_1h, c_2h, c_3h, c_2h$, and $c_1h$, with $c_1 = 0.28, c_2 = 0.62546642846767004501$ and $c_3 = 1 - 2c_1 - 2c_2$. In combination with constant time-steps, the pure N-body model in TIDYMESS is thus fourth-order and symplectic.

For the conservative tidal model, we need to extend the Leapfrog step to include an update of the spin angular momenta using the torque. Furthermore, the accelerations and torques are not only functions of positions, but also of the inertia tensors:

$$ \begin{align} \mathbf{r}_{\frac{1}{2}} &= \mathbf{r}_0 + \mathbf{v}_0 \frac{h}{2} \quad (\text{Drift}) \notag \\ \mathbf{a}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) &= \mathbf{a}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) \quad (\text{Acceleration}) \notag \\ \mathbf{T}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) &= \mathbf{T}_{\frac{1}{2}} \left( \mathbf{r}_{\frac{1}{2}} \right) \quad (\text{Torque}) \notag \\ \mathbf{v}_{1} &= \mathbf{v}_{0} + \mathbf{a}_{\frac{1}{2}} \frac{h}{2} \quad (\text{Kick orbit}) \notag \\ \mathbf{L}_{1} &= \mathbf{L}_{0} + \dot{\mathbf{T}}_{\frac{1}{2}} \frac{h}{2} \quad (\text{Kick spin}) \notag \\ \mathbf{r}_{1} &= \mathbf{r}_{\frac{1}{2}} + \mathbf{v}_{1} \frac{h}{2} \quad (\text{Drift}) \tag{78} \end{align} $$

The problem here is that in the calculation of the acceleration and torque, the inertia tensor, $\mathcal{I}_{\frac{1}{2}}$, is not known. This is because the
inertia tensor is not constant, contrary to models which consider a constant principal moment of inertia. In the conservative model, the inertia tensor can be split into a tidal part and a centrifugal part: $\mathcal{I}_{\frac{1}{2}} = \mathcal{I}_e(\mathbf{r}_{\frac{1}{2}}) + \mathcal{I}_e(\mathbf{\Omega}_{\frac{1}{2}})$. The positions, $\mathbf{r}_{\frac{1}{2}}$, are known, but not the spin frequencies, $\mathbf{\Omega}_{\frac{1}{2}}$. But in order to estimate $\mathbf{\Omega}_{\frac{1}{2}}$, through the relation $\mathbf{\Omega}_{\frac{1}{2}} = \mathcal{I}_{\frac{1}{2}}^{-1}\mathbf{L}_{\frac{1}{2}}$, we require $\mathcal{I}_{\frac{1}{2}}$ again, thus producing an endless loop. Furthermore, $\mathbf{L}_{\frac{1}{2}}$ is not known either. Inspired by the Auxiliary-Vector-Algorithm (ÅVA; Hellström & Mikkola 2010), we resolve these issues by defining an auxiliary angular momentum, $\mathbf{K}$. A reversible map is constructed as follows:

$$ 
\begin{aligned}
\mathbf{r}_{\frac{1}{2}} &= \mathbf{r}_0 + \mathbf{v}_0 \frac{h}{2} && (\text{Drift}) \\
\mathbf{T}_{\frac{1}{2},0} &= (\mathbf{r}_{\frac{1}{2}}, \mathbf{L}_0) && (\text{Torque}) \\
\mathbf{K}_{\frac{1}{2}} &= \mathbf{K}_0 + \mathbf{T}_{\frac{1}{2},0} \frac{h}{2} && (\text{Kick spin}) \\
\mathbf{a}_{\frac{1}{2}} &= (\mathbf{r}_{\frac{1}{2}}, \mathbf{K}_{\frac{1}{2}}) && (\text{Acceleration}) \\
\mathbf{T}_{\frac{1}{2}} &= (\mathbf{r}_{\frac{1}{2}}, \mathbf{K}_{\frac{1}{2}}) && (\text{Torque}) \\
\mathbf{v}_1 &= \mathbf{v}_0 + \mathbf{a}_{\frac{1}{2}} \frac{h}{2} && (\text{Kick orbit}) \\
\mathbf{L}_1 &= \mathbf{L}_0 + \dot{\mathbf{T}}_{\frac{1}{2}} \frac{h}{2} && (\text{Kick spin}) \\
\mathbf{T}_{\frac{1}{2},1} &= (\mathbf{r}_{\frac{1}{2}}, \mathbf{L}_1) && (\text{Torque}) \\
\mathbf{K}_1 &= \mathbf{K}_{\frac{1}{2}} + \mathbf{T}_{\frac{1}{2},1} \frac{h}{2} && (\text{Kick spin}) \\
\mathbf{r}_1 &= \mathbf{r}_{\frac{1}{2}} + \mathbf{v}_1 \frac{h}{2} && (\text{Drift})
\end{aligned}
\quad (79) $$

Note that for the acceleration and torque steps, we used the angular momentum as the argument instead of the inertia tensor. The inertia tensor is solved for iteratively. For example, the inertia tensor corresponding to $\mathbf{T}_{\frac{1}{2},0}(\mathbf{r}_{\frac{1}{2}}, \mathbf{L}_0)$ is calculated as follows. We initially set $\mathbf{\Omega} = \mathbf{\Omega}_0$, and then we iterate:

$$ 
\left\{ 
\begin{aligned}
\mathcal{I} &= \mathcal{I}_e(\mathbf{r}_{\frac{1}{2}}) + \mathcal{I}_e(\mathbf{\Omega}) \\
\mathbf{\Omega} &= \mathcal{I}^{-1}\mathbf{L}_0
\end{aligned}
\right. 
\quad (80) $$

In TIDYMESS, the number of iterations can be set by the user, and is a balance between the required level of reversibility, and computational expense. In principle, with a sufficient number of iterations, the conservative tide model is symplectic. After performing several tests however, we noticed that the angular momenta, $\mathbf{L}$ and $\mathbf{K}$, can sometimes rapidly diverge. We resolve this by resetting $\mathbf{K} = \mathbf{L}$ after each ‘McLachlan’ step. Although this affects the level of reversibility, in practice we still obtain excellent results. The final inertia tensor, $\mathcal{I}_1$, is only calculated at the snapshot output interval, using the same iteration scheme.

The integration scheme for the linear tide model is very similar to that of the conservative model. The update of the inertia tensor now also includes a rotational component due to the tidal time lag, as well as the derivative of the equilibrium tensor, $\dot{\mathcal{I}}^e$. We approximate the derivatives by

$$ 
\begin{aligned}
\dot{\mathcal{I}}_{\frac{1}{2},0}^e &= (\mathcal{I}_e(\mathbf{r}_{\frac{1}{2}}, \mathbf{L}_0) - \mathcal{I}_e(\mathbf{r}_0, \mathbf{L}_0)) \frac{2}{h} \\
\dot{\mathcal{I}}_{\frac{1}{2}}^e &= (\mathcal{I}_e(\mathbf{r}_{\frac{1}{2}}, \mathbf{K}_{\frac{1}{2}}) - \mathcal{I}_e(\mathbf{r}_0, \mathbf{L}_0)) \frac{2}{h} \\
\mathcal{I}_1^e &= (\mathcal{I}_e(\mathbf{r}_1, \mathbf{L}_1) - \mathcal{I}_e(\mathbf{r}_{\frac{1}{2}}, \mathbf{K}_{\frac{1}{2}})) \frac{2}{h}
\end{aligned}
\quad (81) $$

The Creep tidal model with direct integration has a different integration map. Before the first integration step, we calculate the accelerations, $\mathbf{a}$, torques, $\mathbf{T}$, and deformation tensor, $\dot{\mathbf{K}}$. We also introduce the auxiliary inertia tensor, $\mathcal{K}$. We subsequently integrate

the system as follows:

$$ 
\begin{aligned}
\mathbf{r}_{\frac{1}{2}} &= \mathbf{r}_0 + \mathbf{v}_0 \frac{h}{2} && (\text{Drift}) \\
\mathbf{K}_{\frac{1}{2}} &= \mathbf{K}_0 + \mathbf{T}_{0,0} \frac{h}{2} && (\text{Kick spin}) \\
\mathcal{K}_{\frac{1}{2}} &= \mathcal{K}_0 + \dot{\mathcal{I}}_0 \frac{h}{2} && (\text{Kick shape}) \\
\mathbf{\Omega}_{\frac{1}{2}} &= \mathcal{K}_{\frac{1}{2}}^{-1} \mathbf{K}_{\frac{1}{2}} && (\text{Spin frequency}) \\
\frac{\mathbf{a}_{\frac{1}{2}}}{\dot{\mathbf{I}}_{\frac{1}{2}}} &= (\mathcal{K}_{\frac{1}{2}}, \mathbf{r}_{\frac{1}{2}}) && (\text{Acceleration}) \\
\mathbf{T}_{\frac{1}{2}} &= (\mathcal{K}_{\frac{1}{2}}, \mathbf{r}_{\frac{1}{2}}) && (\text{Torque}) \\
\dot{\mathbf{I}}_{\frac{1}{2}} &= (\mathcal{K}_{\frac{1}{2}}, \mathbf{r}_{\frac{1}{2}}, \mathbf{\Omega}_{\frac{1}{2}}) && (\text{Deformation}) \\
\mathbf{v}_1 &= \mathbf{v}_0 + \mathbf{a}_{\frac{1}{2}} \frac{h}{2} && (\text{Kick orbit}) \\
\mathbf{L}_1 &= \mathbf{L}_0 + \dot{\mathbf{T}}_{0,0} \frac{h}{2} && (\text{Kick spin}) \\
\mathcal{I}_1 &= \mathcal{I}_0 + \dot{\mathcal{I}}_{0,0} \frac{h}{2} && (\text{Kick shape}) \\
\mathbf{r}_1 &= \mathbf{r}_{\frac{1}{2}} + \mathbf{v}_1 \frac{h}{2} && (\text{Drift}) \\
\mathbf{\Omega}_1 &= \mathcal{I}_1^{-1} \mathbf{L}_1 && (\text{Spin frequency}) \\
\frac{\mathbf{a}_1}{\dot{\mathcal{I}}_1}, \mathbf{r}_1) && (\text{Acceleration}) \\
\mathbf{T}_1(\mathcal{I}_1, \mathbf{r}_1) && (\text{Torque}) \\
\dot{\mathcal{I}}_1(\mathcal{I}_1, \mathbf{r}_1, \mathbf{\Omega}_1) && (\text{Deformation}) \\
\mathbf{K}_1 &= \mathbf{K}_{\frac{1}{2}} + \mathbf{T}_{1,0} \frac{h}{2} && (\text{Kick spin}) \\
\mathcal{K}_1 &= \mathcal{K}_{\frac{1}{2}} + \dot{\mathcal{I}}_{1,0} \frac{h}{2} && (\text{Kick shape})
\end{aligned}
\quad (82) $$

A consequence of directly using the derivative of the inertia tensor, is that the time-step cannot be made too large, since it has to resolve the orbital, spin and deformation time-scales. This becomes problematic for rapidly spinning bodies, or tides in the linear regime. As a resolution to this limitation, we implement a tidal model based on the general and optimized approximation of the Creep model (see Section 3.2). The integration map is the same as for the linear tide model, both requiring the equilibrium tensor at the beginning and end of a time-step. The net deformation is split into a tidal deformation and rotational part. In the positive time direction, we first perform the tidal deformation followed by the rotation. In the negative time direction, the inverse sequence is applied.

Although the tidal deformation equation (71) is time reversible, in a numerical sense it is not. For values of $\Delta t/\tau \gg 1$, the exponents become small enough, such that any memory of initial conditions gets lost in the numerical noise. This affects the performance of the fourth-order composition of the integrator, which includes a negative time-step. In order to resolve this problem, we check the ratio of $\Delta t/\tau$ for each individual body. If this ratio is larger than some critical value, $f_{\text{crit}}$, than we use the equation from the linear tide model instead. We find that a value of $f_{\text{crit}} = 7$ to be optimal. In a future version of TIDYMESS, we plan to implement a symplectic integration map applicable to arbitrary N-body configurations, where all substeps are in the same direction as the integration (e.g. Chambers & Murison 2000; Chin & Chen 2005).

### 3.4 Time-step functions

TIDYMESS implements three time-step criteria:

(i) Constant

(ii) Adaptive,

(iii) Adaptive, weighted, and symmetrized.

In the adaptive case, we adopt the minimum of the pairwise free fall and interparticle flyby time-scales (e.g. Pelupessy, Jänes & Portegies
Zwart 2012). The magnitude of the time-steps can be controlled by the constant time-step multiplicative factor, $\eta$.

Adapting the time-step size affects the advantageous conservation properties of symplectic integrators. One method to improve this is by symmetrizing the time-steps (e.g. Dehnen 2017). Symmetrization methods are compatible with time-step functions which are smooth and continuous. The latter can be accomplished by using adaptive time-steps which are weighted over all pairs. We adopt the geometric symmetrization scheme in combination with the direct weighting method, that is, methods A2 and W2 as defined by Boekholt, Vaillant & Correia (2023). For the direct Creep model, we extend the pairwise, orbital time-steps to include spin periods and fluid relaxation time-scales. The weighting method however, is only applied to the orbital time-step functions. If the symmetrized time-step, $h$, diverges from the time-step function, $T$, by a factor 4 or more, then we reset the symmetrized time-step to $h = T$ (Hands et al. 2019).

## 3.5 Input/output

TIDYMESS reads in two files: the parameter file and the initial condition file. Their default names are *tidymess.par* and *tidymess.ic*, but these can be user defined. The parameter file specifies various input parameters necessary to run a simulation. Here, the user defines which tidal model to use (Table 1), the Post-Newtonian order (Section 2.3.1), whether to turn on magnetic braking effects (Section 2.3.2), and how to handle collisions (Section 2.3.3). The simulation time is specified here, as well as the time-step function and time-step parameter. The user can specify the output directory, the snapshot time interval, and the unit system to be used in the output snapshots. Various output formats are available, including an output file per body, output file per snapshot, or all snapshots inside a single file. The snapshots are always in Cartesian coordinates in the inertial frame specified by the initial condition. Other functionalities include a flag for restarting a simulation, setting the initial shapes of the bodies (sphere or equilibrium shape), varying the number of bodies to include from the initial condition file, specifying the type of coordinates used in the initial condition file, for positions and velocities (Cartesian, elliptical astrocentric, or elliptical Jacobian), as well as for spins (which can be given in the inertial frame or with respect to the bodies' orbital plane around the primary body with mass $m_0$). The number of iterations to improve the reversibility of the integrator can also be defined (see Section 3.3). All input parameters can also be defined on the terminal command line. In case parameters are defined in both the parameter file and the terminal, then the latter value is assigned.

The initial condition file consists of a table with the properties of the bodies in the system. The internal properties (mass $m$, mean radius $R$, moment of inertia factor $\tilde{\epsilon}$, fluid Love number $k_f$, fluid relaxation time $\tau$, and magnetic braking constant ($\alpha$) can all be set in various physical units. The positions and velocities, as well as the initial spin vectors, can be defined here, in the coordinate frame specified in the parameter file (Cartesian coordinates or orbital elements). Each body can also be assigned a name tag. Hence, in order to start a new simulation, the user must first gather the necessary data from an external source, and convert that to a table of initial conditions in the *TIDYMESS* format. Then, the user updates the parameter file in order to define the simulation that they wish to run. The outputted snapshots are text files, which are easily read by other programs, such as *PYTHON*. Binary files are also written out for backup and restart purposes. A summary log file is also kept, with information such as the duration of the simulation, and the level of conservation of various quantities.

# 4 DEMONSTRATIONS OF TIDYMESS

In this section we provide six demonstrations of TIDYMESS, each highlighting a different tidal feature or application. In Section 4.1, we consider the current-day Earth-Moon system, and model the orbital recession speed of the Moon due to tides. In Section 4.2, we focus on the deformations of the heartbeat star KOI-54. In Section 4.3, we model the eccentric exoplanet HD80606b, and compare to previous results from the literature. In Section 4.4, we consider the three-body problem consisting of the early Sun-Earth-Moon system, and we reproduce the ejection and eviction events. In Section 4.5, we model a resonant chain of three planets, and the formation of a hot Jupiter. In Section 4.6, we consider the tidal evolution of a stellar association, until the formation of the first, synchronized binary system. We note that the aim of these examples is to provide a gallery of demonstrations of TIDYMESS, rather than providing in-depth analyses of the results, which can mostly be found in the given references.

## 4.1 Earth-Moon recession speed

In our first demonstration of TIDYMESS, we consider the current-day Earth-Moon system. We simplify the system to a one-body system, that is, we treat the Moon as a point-mass body and only consider tides on Earth. We simplify the problem somewhat more by assuming a zero obliquity, and a circular orbit. The initial conditions are taken from NASA's fact sheets (Williams 2016) and given in Appendix Table B1. We evolve the system for 1000 yr, and measure the semimajor axis as a function of time. The slope of a linear fit gives the average recession speed of the Moon from the Earth. TIDYMESS also outputs a log file, including the simulation wall-clock time.

In Fig. 1, we plot the recession speed as a function of the tidal parameter, $\delta = \omega\tau$. Here, $\omega$ is the tidal forcing frequency given by $\omega = 2(\Omega - n)$, with $n$ the orbital frequency (mean motion). Starting with the linear tide model, we observe a linear increase of the recession speed of the Moon as a function of tidal parameter. Intuitively, a larger tidal parameter corresponds to a larger lag of the tidal bulge, resulting in a larger torque. In the linear model, the increase in the recession speed is maintained, even up to very large tidal parameters. However, in the limit of infinite fluid relaxation times, the shapes of the bodies do not change, thus resulting in a conservative system without tidal dissipation. Hence, there should be a non-linear regime in which the tidal dissipation (and therefore the recession speed) decreases again. The Creep models exhibit such a feature. The peak recession speed occurs when the forcing frequency resonates with the fluid relaxation time, that is, $\omega\tau = 1$. The curve is symmetric around this peak. We confirm that in the linear regime, the Creep model reduces to the linear tide model. Furthermore, the two Creep model implementations are mutually consistent, except for a small deviation near the peak of the curve. This follows from the approximations made in the Creep TIDYMESS algorithm. We also ran the simulations with the conservative tidal model, and we confirm a negligible residual recession speed $\le 10^{-5}$ [cm yr$^{-1}$], independent of the tidal parameter.

A large difference is observed in the performance of the two Creep models. In Fig. 1, we also plot the wall-clock time of the simulation, normalized by the wall-clock time for the pure N-body model. The linear model is about a factor 3 slower than a normal N-body code. The Creep TIDYMESS model is 4–5 times slower,<fnref n="3" />

<footnote n="3">³This benchmark was obtained using 1 iteration in the integration map. Without iterations, the TIDYMESS Creep model is 3–4 times slower.</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1.** Tides in the current Earth-Moon system. As a function of the tidal parameter, $\delta = \omega\tau$, we plot the slow down factor of the wall-clock time, relative to simulations without tides (top panel), and the recession speed of the Moon in terms of the derivative of the semimajor axis, $\dot{a}$ (bottom panel). For the Creep models, we observe a maximum recession speed at $\omega\tau = 1$. We confirm the linear model is consistent for $\omega\tau \ll 1$, but blows up in the opposite regime. This example also illustrates the performance of the TIDYMESS Creep model, which is independent of the fluid relaxation time. This is in contrast to the direct Creep model, which becomes excessively expensive in the linear regime.

but the performance is independent of the tidal parameter. The direct Creep model is about 30 times slower in the non-linear regime, which is explained by the rapid spin of Earth compared to its orbit (1 month/1 d $\approx$30). However, into the linear regime, the wall-clock time increases linearly with decreasing fluid relaxation time-scale. This example demonstrates the advantage of the new TIDYMESS implementation of the Creep model.

## 4.2 KOI-54: eccentric binary with a heartbeat

Apart from the usual phase space coordinates, TIDYMESS also integrates and outputs the spin and inertia tensor coordinates in an inertial reference frame. This makes it possible to project the extended bodies onto a plane, and measure their projected surface areas as a function of time.

Heartbeat stars are known for their periodic flux increases. They consist of eccentric binary stars, and during pericentre passage, tidal deformations cause a brief spike in the observed flux. One observed example is the system KOI-54 (e.g. Welsh et al. 2011). We adopt their derived system parameters, which are repeated in Appendix Table B2. We perform brief simulations for a single orbital period, using the direct Creep tidal model. The projected surface area is calculated according to the method described in Appendix A. Observed heartbeat flux profiles can be explained by contributions

<figure><img src="image_4.png" /></figure>

**Figure 2.** Heartbeat of the eccentric binary system KOI-54. We plot the normalized, projected surface area of the stars during pericentre passage. The normalization is taken to be the initial value at apocentre starting from their equilibrium shapes, which is not necessarily the same as the sum of the two unperturbed stars. The profile depends sensitively on the eccentricity (top panel, $\tau = 10^{-4}$ [yr]), as well as the fluid relaxation time (bottom panel, $e = 0.8335$). The latter introduces an asymmetry in the profile.

from both tidal distortion and irradiation/reflection effects (Welsh et al. 2011). Here, we will only focus on the contribution from the tidal deformation of the bodies as predicted by the Creep model.

In Fig. 2, we plot the normalized increase of the projected surface area as a function of mean anomaly. We first fix the fluid relaxation time of the stars to $\tau = 10^{-4}$ [yr], and vary the eccentricity (top panel). We confirm the finding by Welsh et al. (2011), that the flux (which is proportional to surface area), depends sensitively on the eccentricity. This is explained by the steep dependence of tides on separation. Secondly, we fix the eccentricity to $e = 0.8335$, and systematically vary the fluid relaxation time (bottom panel). We observe that the fluid relaxation time of the stars affects the symmetry of the profile. The larger the value of $\tau$, the longer it takes for the star to restore to its non-perturbed shape after pericentre passage. Detailed observations of stellar shapes thus hold stringent information on their tidal response.

## 4.3 HD80606b: a very eccentric Jupiter

The exoplanet HD80606b is highly eccentric and its future tidal evolution will most likely evolve it into a hot Jupiter. Here, we reproduce simulations by Correia et al. (2014) and Boué et al. (2016),
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Tidal evolution of exoplanet HD80606b. We plot the time evolution of various initial obliquities (top panel), semimajor axis and eccentricity (2nd row), ratio of spin to orbital frequency (3rd row), and two parameters describing oblateness ($J_2$) and prolateness ($\epsilon$) of the planet (bottom panel). We vary the fluid relaxation time per column, such that we sample different tidal regimes: linear (left), peak dissipation (middle), and non-linear/inverse (right). In each case, the planet attains the synchronized, equilibrium state.</figcaption></figure>

who adopted the same tidal Creep model as TIDYMESS. The initial conditions are repeated in Appendix Table B3. We evolve the planet HD80606b with three different values for the fluid relaxation time: $\tau = 10^{-4}$, $10^{-2}$, and 1 yr. These values represent the linear regime, the non-linear regime, and the transitional regime in between.

Starting with a zero obliquity, we evolve the systems until tidal synchronization is approximately reached. We confirm in Fig. 3 (2nd row), that the semimajor axes shrink until a fixed value determined by angular momentum conservation. The eccentricities also decrease

to zero. At the same time, we observe that the ratio between the spin and orbital frequencies evolves to unity (Fig. 3, 3rd row). We successfully reproduce the ‘staircase’ pattern found by Correia et al. (2014), which is a consequence of passing through a spectrum of spin-orbit resonances in the non-linear regime. The deformation of the planet evolves over a similar time-scale as the orbits and spins. We measure the planet’s Stokes coefficients $J_2$ and $\epsilon = (C_{22}^2 + S_{22}^2)^{1/2}$, which represent the oblateness and prolateness, respectively (Fig. 3, bottom row). Large variations are observed in the prolateness
due to the eccentric orbit and the steep dependence of tides on separation. The variations are physical, although the visualization exhibits some small irregularities due to a finite sampling of the output data. Towards the end, an equilibrium figure is reached, which is approximately independent of fluid relaxation time. A non-zero initial obliquity rapidly decreases to zero on a time-scale much shorter than the orbital time-scale (Fig. 3, top row). Hence, large initial obliquities do not change the subsequent evolution much.

## 4.4 Early evolution of the Sun-Earth-Moon system

The Moon is currently receding from Earth with a rate of about $4 \text{ cm yr}^{-1}$ (Dickey et al. 1994). This implies that in the past the Moon was closer to Earth. Origin scenarios, such as a giant impact (e.g Canup 2004), indeed consider an ‘initial’ orbit after formation, just beyond the Roche limit. Supersynchronous rotation of the Earth and Moon subsequently resulted in an orbital expansion. During this expansion, multiple interesting resonances are crossed, which excite the eccentricity and the inclination.

We reproduce the results from Touma & Wisdom (1998), who analyse the resonant phenomena of ejection and eviction. In the case of ejection, the eccentricity of the Earth-Moon orbit is excited due to a resonance between Earth’s orbital frequency around the Sun and the Moon’s pericentre precession rate. This is followed by eviction, where the inclination of the Moon is excited due to an eccentricity–inclination resonance (Touma & Wisdom 1998). These two resonance passages provide a stringent validation test for TIDYMESS.

We adopt the initial conditions from Touma & Wisdom (1998), which are repeated here in Appendix Table B4. The Moon is initially in the equatorial plane of the Earth in a near circular orbit. The Earth–Moon separation starts off at $3.5 \text{ Earth radii}$, while Earth’s obliquity with respect to its orbit around the Sun is set to $10^\circ$. The initial rotation period of Earth is chosen to be $5.0 \text{ h}$. In our experiment, we will only consider the tidal deformation of the Earth, whose fluid relaxation time-scale is set to one second, that is, linear tidal regime. Although this time-scale seems rather short, it reproduces the average recession rate of the Moon of $1 \text{ km yr}^{-1}$, as adopted by Touma & Wisdom (1998).

In Fig. 4, we plot the evolution of the eccentricity and inclination as a function of semimajor axis for the Earth-Moon orbit. We reproduce the increase of the eccentricity to a value of about 0.5, near a separation of ~4.5 Earth radii. We also confirm that the eviction phase occurs at a separation between 5.5 and 6 Earth radii, resulting in an increase in the inclination of the Moon to a value between 2° and 3°. Hence, TIDYMESS can be used to model resonant phenomena in the Solar System.

## 4.5 Chaotic planet–planet scattering

TIDYMESS covers a niche regime of parameter space where orbits are chaotic and tides can become non-linear. As a demonstration, we consider a host star and a three-planet system, which are all mutually tidally interacting. If the planets start out in a resonant chain, that is, each neighbouring pair of planets is in orbital resonance, instabilities can grow and excite eccentricities. When two planets are nearly orbit-crossing, strong gravitational deflections will lead to a variety of possible outcomes. These include planet ejection, rearrangement of orbits, giant impacts, and the formation of a hot Jupiter (e.g. Beaugé & Nesvorný 2012). In the latter case, one planet is scattered towards the host star. If the perihelion distance is of

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Evection resonance leads to eccentricity excitation (top), and eviction resonance leads to inclination excitation (bottom), during the early evolution of the Sun-Earth-Moon system. These results are consistent with those from Touma & Wisdom (1998).</figcaption></figure>

order a few times the stellar radius, then the planet can be tidally captured by the star. The subsequent evolution might resemble that of exoplanet HD80606b (see Section 4.3). We note that TIDYMESS can model both the initial multiplanet interaction phase, as well as the subsequent tidal synchronization of the hot Jupiter, without any need for a transition between numerical methods. Furthermore, it is possible to model bodies in the same system with different tidal models. In our experiment, we will assign the host star the conservative tide model, while the planets respond according to the TIDYMESS Creep model.

We define our own initial condition for this experiment, which might be common among young exoplanetary systems. The host star is taken to be a Sun-like star. The planets have masses 1, 2, and 4 times the mass of Jupiter, with the lightest planet nearest the star, and the heaviest planet on the outside. The resonant chain consists of a 2:1 and a 3:2 orbital resonance between the inner and outer pair, respectively. The inner planet has a semimajor axis of 2 au. The spin periods of the bodies are taken to be the current-day values of the Sun and Jupiter. Due to the presence of chaos, small changes in the initial realization have a big influence on the final outcome. We randomly varied the eccentricities between 0 and 0.04, and inclinations between $\pm 0.4^\circ$. The initial orbital phases were also randomly sampled. For our demonstration, we chose a system which produced both a hot Jupiter and a planet ejection.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** We plot the chaotic orbital separations of three planets from their host star (top left). One planet is ejected after about 0.18 Myr (red curve), while another planet is scattered inwards, after which its orbit is tidally synchronized, that is, it has become a hot Jupiter (black curve). This planet's final spin-orbit angle is about 16° (top right), its spin frequency is synchronized to its orbital frequency (bottom left), and its obliquity is zero (bottom right). The planet with the relatively unaffected orbital separation (blue curve), obtains a final obliquity of about 24° demonstrates that planet-planet scattering can change the obliquities of planets. This planet also exhibits an oscillation of the spin-orbit angle (top right, blue curve), which is explained by the spin precession and nutation of the host star, due to its tidal coupling with the hot Jupiter.</figcaption></figure>

In Fig. 5, we plot the time evolution of various spin-orbit properties of the three planets. We observe in the top left panel that the distance of the planets to the host star varies chaotically until about 0.18 Myr, when one planet is ejected from the system (red curve, initially the middle planet). Some time later just before 0.20 Myr, a planet is scattered inwards, followed by a rapid tidal synchronization with the host star (black curve, initially the inner and lightest planet). The final semimajor axis of the hot Jupiter is about 4–5 stellar radii. The initial outer planet, which is also the highest mass, remains relatively unaffected in its orbit.

In parallel to the orbital evolution, the hot Jupiter also synchronizes its spin (bottom left panel), while the spin-orbit angle becomes fixed at approximately 16° (top right panel). The hot Jupiter's obliquity was greatly excited during the multiplanet scattering phase, up to about 100°. Tidal evolution after the scatter however, synchronizes the obliquity rapidly to 0° (bottom right panel). Interestingly, the most massive planet obtains a final obliquity of about 24° (blue curve, bottom right panel). This result shows that planet-planet scattering is a channel for producing planets with large obliquities. Furthermore, the same planet shows large variations in the final spin-orbit angle by about 20°. This is explained by the spin precession and nutation of the star, which in turn is driven by the tidal coupling between the star and the hot Jupiter.

## 4.6 Close binary formation in a stellar association

TIDYMESS can model the tidal evolution of an arbitrary configuration of stars and planets. These include moons around planets, a host star with multiple planets, binary stars with planets, hierarchical triple stars, or even stellar associations up to a few hundred stars.<fnref n="4" /> Tidal effects play an important role in stellar systems, and can result in the formation of binary stars by tidal capture (e.g. Press & Teukolsky 1977), and potentially rearrange stellar spins.

Here, we will consider a stellar association consisting of 32 single, solar mass stars in a Plummer configuration (Plummer 1911). This is the minimal number of stars for a system with two-body relaxation effects (Portegies Zwart et al. 2021). The virial radius of the association is set to $10^5$ solar radii (for denser systems, we observed the common occurrence of collisions, that is, formation of blue stragglers). The tidal model is the TIDYMESS Creep model, with a fluid relaxation time arbitrarily set to 0.01 yr (3.65 d), which is of the same order as the orbital period of very close binary stars. We evolve the association for 10 kyr (approximately 12 dynamical times).

<footnote n="4">⁴A future version of TIDYMESS is planned to have a parallel implementation allowing larger values for N.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Formation of a close binary star from a stellar association. A three-body encounter produces an ejected binary star after about 4 kyr. This eccentric binary subsequently tidally evolves into a circular and synchronized binary star with an orbital period of about 6 d.</figcaption></figure>

In Fig. 6, we visualize the formation of a close binary star through a two-part process. Up to about 4 kyr, the stars are moving chaotically through the association, ocassionally encountering other stars. The separation between the two stars which eventually form the binary, varies chaotically around a value given by the size of the association (10⁵ solar radii, blue curve in top panel). Then a three-body encounter occurs, resulting in the ejection of a star and an eccentric binary. In Fig. 6, we observe that the distance of the binary from the center of mass of the association starts to linearly increase after 4 kyr (top panel, orange curve). At the same time, the tidal evolution of the binary results in orbital shrinkage and circularization. The final orbit has an orbital period of about 6 d. For the two binary components, we also measure the spin periods and orbital periods in time (bottom panel). The spins of both stars are indistinguishable. We estimate the orbital period in the initial unbound phase as the separation divided by the relative speed times 2π. In the initial phase up to 4 kyr, the stars are supersynchronous. Spins are greatly affected however in the binary system, through tidal synchronization. The spins of the remaining stars in the association remained barely affected on this short time-scale. It remains an open question whether tides in large-N stellar systems can redistribute the spin vector distribution.

## 5 CONCLUSION

We present the new N-body code TIDYMESS, which stands for the 'Tidal Dynamics of Multibody ExtraSolar Systems'. Due to the combination of a direct integration method and the tidal Creep model, this code covers a niche regime of parameter space where orbits become chaotic and where tides can be non-linear. A new implementation of the Creep model, which we call the TIDYMESS Creep model, consists of an integration method for the tidal deformations, which is independent of the spin and fluid relaxation time-scales. This greatly benefits the performance of the Creep model in the linear regime, or in the presence of rapidly spinning bodies.

We provide six different demonstrations of TIDYMESS, with applications to the tidal evolution of moons, planets, and stars. The code is open-source and provides the community with a new method for studying the chaotic, tidal evolution of Solar System bodies, exoplanets, stellar multiples, and associations. Comparison to other tidal N-body codes, combined with stringent constraints from current and upcoming observations, will allow us to scrutinise the physics of tidal dissipation and its manifestations in various astronomical bodies.

## ACKNOWLEDGEMENTS

We kindly thank the referee, John Chambers, for the suggested improvements to the paper. This work was supported by funds from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement number 638 435 (GalNUC). Other support includes CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal. The numerical simulations were performed at the Laboratory for Advanced Computing at University of Coimbra (https://www.uc.pt/lca). The diagrams were produced using PYTHON (Van Rossum & Drake Jr 1995).

## DATA AVAILABILITY

Our new N-body code TIDYMESS is open-source and available through GitHub (https://github.com/tidymess-code) and through the Astrophysical Multipurpose Software Environment (https://github.com/amusecode).

## REFERENCES

* Barker A. J., Ogilvie G. I., 2009, MNRAS, 395, 2268

* Bataille M., Libert A. S., Correia A. C. M., 2018, MNRAS, 479, 4749

* Baumgardt H., Hut P., Heggie D. C., 2002, MNRAS, 336, 1069

* Beaugé C., Nesvorný D., 2012, ApJ, 751, L119

* Bildsten L., Cutler C., 1992, ApJ, 400, L175

* Boekholt T. C. N., Pelupessy F. I., Heggie D. C., Portegies Zwart S. F., 2016, MNRAS, 461, 3576

* Boekholt T. C. N., Portegies Zwart S. F., Valtonen M., 2020, MNRAS, 493, 3932

* Boekholt T. C. N., Vaillant T., Correia A. C. M., 2023, MNRAS, 519, 3281

* Bolmont E., Raymond S. N., Leconte J., Hersant F., Correia A. C. M., 2015, A&A, 583, 116

* Boué G., Correia A. C. M., Laskar J., 2016, Cel. Mech. Dyn. Astron., 126, 31

* Canup R. M., 2004, Icarus, 168, 433

* Chambers J. E., Murison M. A., 2000, AJ, 119, 425

* Chen R., Li G., Tao M., 2021, ApJ, 919, L50

* Chin S. A., Chen C. R., 2005, Cel. Mech. Dyn. Astron., 91, 301
Correia A. C. M., 2018, Icarus, 305, 250

Correia A. C. M., Laskar J., Farago F., Boué G., 2011, Cel. Mech. Dyn. Astron., 111, 105

Correia A. C. M., Boué G., Laskar J., Rodríguez A., 2014, A&A, 571, 50

Ćuk M., El Moutamid M., Tiscareno M. S., 2020, Plan. Sci. J., 1, 22

Dehnen W., 2017, MNRAS, 472, 1226

Dickey J. O., et al., 1994, Science, 265, 482

Einstein A., 1915, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften, 831

Fabrycky D., Tremaine S., 2007, ApJ, 669, L1298

Ferraz-Mello S., 2013, Cel. Mech. Dyn. Astron., 116, 109

Ford E. B., Kozinsky B., Rasio F. A., 2000, ApJ, 535, L385

Fujii M. S., Portegies Zwart S., 2014, MNRAS, 439, 1003

Gaburov E., Lombardi James C., J., Portegies Zwart S., Rasio F. A., 2018, StarSmasher: Smoothed Particle Hydrodynamics code for smashing stars and planets, Astrophysics Source Code Library, preprint (ascl:1805.0100)

Goldstein H., 1950, Classical mechanics. Addison-Wesley, Reading

Goodman J., Heggie D. C., Hut P., 1993, ApJ, 415, L715

Hands T. O., Dehnen W., Gration A., Stadel J., Moore B., 2019, MNRAS, 490, 21

Hellström C., Mikkola S., 2010, Cel. Mech. Dyn. Astr., 106, 143

Hut P., 1981, A&A, 99, 126

Kawaler S. D., 1988, ApJ, 333, L236

Kocsis B., 2022, in Bambi C., ed., Handbook of Gravitational Wave Astronomy. Springer, Singapore, p. 15

Kozai Y., 1962, AJ, 67, 591

Lambeck K., 1980, The Earth's Variable Rotation: Geophysical Causes and Consequences. Cambridge Univ. Press, Cambridge

Laplace P. S., 1784, in Œuvres complètes, XI. Gauthier-Villars, Paris (1895), p. 49

Laskar J., Gastineau M., 2009, Nature, 459, 817

Lidov M. L., 1961, Iskus. sputniky Zemly (in Russian), 8, 5

Lin D. N. C., Bodenheimer P., Richardson D. C., 1996, Nature, 380, 606

Love A. E. H., 1911, Some Problems of Geodynamics. Cambridge Univ. Press, Cambridge

Love A. E. H., 1927, A Treatise on the Mathematical Theory of Elasticity. Dover Publications, New York

Lu T., Rein H., Tamayo D., Hadden S., Mardling R., Millholland S. C., Laughlin G., 2023, preprint (arXiv:2303.00006)

Makino J., 1996, ApJ, 471, L796

Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634

Mardling R. A., Lin D. N. C., 2002, ApJ, 573, L829

McLachlan R., 1995, SIAM J. Sci. Comput., 16, 151

McMillan S., Hut P., Makino J., 1990, ApJ, 362, L522

Mignard F., 1979, Moon and Planets, 20, 301

Miller R. H., 1964, ApJ, 140, L250

Murray N., Hansen B., Holman M., Tremaine S., 1998, Science, 279, 69

Naoz S., 2016, ARA&A, 54, 441

Newton I., 1687, Philosophiae Naturalis Principia Mathematica. Auctore Js. Newton. Londini, Jussu Societatis Regiæ ac Typis Josephi Streater. Prostat apud plures Bibliopolas. Anno

Ogilvie G. I., 2009, MNRAS, 396, 794

Peale S. J., 1976, ARA&A, 14, 215

Peltier W. R., 1974, Reviews of Geophysics and Space Physics, 12, 649

Pelupessy F. I., Jänes J., Portegies Zwart S., 2012, New A, 17, 711

Perets H. B., Kratter K. M., 2012, ApJ, 760, L99

Plummer H. C., 1911, MNRAS, 71, 460

Poincaré H., 1905, Leçons de Mécanique Céleste, Tome I. Gauthier-Villars, Paris

Portegies Zwart S., McMillan S., 2018, Astrophysical Recipes: The art of AMUSE, IOP Publishing Ltd, Bristol

Portegies Zwart S. F., Boekholt T. C. N., Por E., Hamers A. S., McMillan S. L. W., 2022, A&A, 659, A86

Portegies Zwart S. F., Boekholt T. C. N., Por E. H., Hamers A. S., McMillan S. L. W., 2022, A&A, 659, 86

Press W. H., Teukolsky S. A., 1977, ApJ, 213, L183

Quillen A. C., Kueter-Young A., Frouard J., Ragozzine D., 2016, MNRAS, 463, 1543

Ragazzo C., Ruiz L. S., 2017, Cel. Mech. Dyn. Astron., 128, 19

Rein H., Payne M. J., Veras D., Ford E. B., 2012, MNRAS, 426, 187

Samsing J., D'Orazio D. J., 2018, MNRAS, 481, 5445

Skumanich A., 1972, ApJ, 171, L565

Stone N. C., Leigh N. W. C., 2019, Nature, 576, 406

Sussman G. J., Wisdom J., 1992, Science, 257, 56

Suzuki M., 1990, Physics Letters A, 146, 319

Tamayo D., Rein H., Shi P., Hernandez D. M., 2020, MNRAS, 491, 2885

Tanikawa A., Hut P., Makino J., 2012, New A, 17, 272

Toonen S., Portegies Zwart S., Hamers A. S., Bandopadhyay D., 2020, A&A, 640, 16

Toonen S., Boekholt T. C. N., Portegies Zwart S., 2022, A&A, 661, 61

Touma J., Wisdom J., 1998, Astron. J., 115, 1653

Trani A. A., Spera M., 2023, IAU Symp., 362, 404

Tsiganis K., Gomes R., Morbidelli A., Levison H. F., 2005, Nature, 435, 459

Valente E. F. S., Correia A. C. M., 2022, A&A, 665, A130

Van Rossum G., Drake Jr F. L., 1995, Python reference manual. Centrum voor Wiskunde en Informatica, Amsterdam, Netherlands

Vergara M. Z. C., Schleicher D. R. G., Boekholt T. C. N., Reinoso B., Fellhauer M., Klessen R. S., Leigh N. W. C., 2021, A&A, 649, 160

Wang Y., Mora P., Liang Y., 2022, Cel. Mech. Dyn. Astron., 134, 42

Weber E. J., Davis L., Jr, 1967, ApJ, 148, L217

Welsh W. F., Orosz J. A., Aerts C., et al., 2011, Astrophys. J. (Supp.), 197, 4

Will C. M., 2014, Phys. Rev. D, 89, 044043

Williams D. R., 2016, NASA NSSDCA Planetary Fact Sheets. https://nssdc.gsfc.nasa.gov/planetary/planetfact.html

Wisdom J., 1982, AJ, 87, 577

Wu Y., Murray N., 2003, ApJ, 589, 605

Yoshida H., 1993, Cel. Mech. Dyn. Astron., 56, 27

Zharkov V. N., Trubitsyn V. P., 1978, Physics of planetary interiors, Astronomy and Astrophysics Series. Pachart Publishing House, Tucson

## APPENDIX A: SHAPE OF THE DEFORMED STAR

The shape (or the figure) of a large celestial body is usually well described by a reference ellipsoid (quadrupolar approximation), that is, a mathematically-defined closed quadric surface that is a 3D analogue of an ellipse. The standard equation of a triaxial ellipsoid centred on the origin of a Cartesian coordinate system (X, Y, Z) and aligned with the axes is

$$ \frac{X^2}{a^2} + \frac{Y^2}{b^2} + \frac{Z^2}{c^2} = 1, \qquad (A1) $$

where *a*, *b*, and *c* are called the semiprincipal axes. More generally, for any reference frame, equation (A1) can be rewritten as

$$ R^T A R = 1, \quad \text{with} \quad A = S \begin{bmatrix} a^{-2} & 0 & 0 \\ 0 & b^{-2} & 0 \\ 0 & 0 & c^{-2} \end{bmatrix} S^T, \qquad (A2) $$

where $R = (X, Y, Z)$ is a generic point at the surface of the ellipsoid, and $S$ is the rotation matrix that connects the general frame with the aligned axes frame (body frame).

Assuming that the deformation of the ellipsoid is small with respect to a sphere of radius $R_0$, we can express $a = R_0(1 + \delta_a)$, $b = R_0(1 + \delta_b)$ and $c = R_0(1 + \delta_c)$, with $\delta \ll 1$. Then, we get

$$ A \approx \frac{1}{R_0^2} \mathbf{I} - \frac{2}{R_0^2} S \begin{bmatrix} \delta_a & 0 & 0 \\ 0 & \delta_b & 0 \\ 0 & 0 & \delta_c \end{bmatrix} S^T, \qquad (A3) $$

and the general equation for the ellipsoid becomes

$$ (1 + \delta_x)X^2 + (1 + \delta_y)Y^2 + (1 + \delta_z)Z^2 + \delta_1XY + \delta_2XZ + \delta_3YZ = R_0^2. \qquad (A4) $$
Using previous expression, the radial distance of the ellipsoid surface, $R = |\mathbf{R}|$, can be expressed as

$$ 
\begin{aligned}
R^2 &= \mathbf{R} \cdot \mathbf{R} = X^2 + Y^2 + Z^2 \\
&= R_0^2 - (\delta_x X^2 + \delta_y Y^2 + \delta_z Z^2 + \delta_1 XY + \delta_2 XZ + \delta_3 YZ) ,
\end{aligned} 
\quad (A5) $$

and since all $\delta \ll 1$, we have $R \approx R_0 + \Delta R$, with

$$ \Delta R \approx -\frac{2}{R_0} (\delta_x X^2 + \delta_y Y^2 + \delta_z Z^2 + \delta_1 XY + \delta_2 XZ + \delta_3 YZ) . \quad (A6) $$

The deformation of the body is a response to the perturbing potential $V_p$ (equation 9). A very convenient way to define this deformation is through the Love number approach (e.g. Love 1911) in which the equilibrium radial displacement $\Delta R_e$ is proportional to the equipotential perturbing surface

$$ \Delta R_e = -h_f V_p(\mathbf{R})/g, \quad (A7) $$

where $g = Gm/R_0^2$ is the surface gravity, $G$ is the gravitational constant, $m$ is the mass of the star, and $h_f = 1 + k_f$ is the fluid second Love number for radial displacement. In general, for stars we have $k_f \ll 1$, thus $h_f \approx 1$. Using expression (15) we get

$$ \Delta R_e \approx -\frac{3}{2k_f} \left[ \hat{\mathbf{R}} \cdot \frac{\mathcal{I}_e}{m R_0} \cdot \hat{\mathbf{R}} \right] . \quad (A8) $$

The equilibrium displacements of the star are computed directly from expressions (17) – (22). However, in a non-equilibrium more general situation, the instantaneous deformation of the star will obey to a similar rheologic law as the one provided by equation (23):

$$ \Delta R + \tau \frac{d}{dt} \Delta R = \Delta R_e . \quad (A9) $$

Thus, the instantaneous deformation of the star is simply given by

$$ \frac{\Delta R}{R_0} = -\frac{3}{2k_f} \left[ \hat{\mathbf{R}} \cdot \frac{\delta \mathcal{I}}{m R_0^2} \cdot \hat{\mathbf{R}} \right] , \quad (A10) $$

which allows us to straightforwardly compute the instantaneous shape of the ellipsoid (equation A4) from the inertia tensor:

$$ 
\begin{aligned}
\delta_x &= \frac{2I_{11} - I_{22} - I_{33}}{4k_f m R_0^2}, \quad \delta_1 = \frac{3I_{12}}{4k_f m R_0^2}, \\
\delta_y &= \frac{2I_{22} - I_{11} - I_{33}}{4k_f m R_0^2}, \quad \delta_2 = \frac{3I_{13}}{4k_f m R_0^2}, \\
\delta_z &= \frac{2I_{33} - I_{11} - I_{22}}{4k_f m R_0^2}, \quad \delta_3 = \frac{3I_{23}}{4k_f m R_0^2}.
\end{aligned} 
\quad (A11) $$

## APPENDIX B: INITIAL CONDITIONS OF THE EXPERIMENTS

For the demonstrations in Section 4 based on observed systems, we provide tables with the initial conditions. These initial conditions and code parameters are also provided in the examples folder in the source code.

<table><caption>**Table B1.** Initial conditions for the Earth-Moon system (Williams 2016) in Section 4.1. We assume that (1) the bodies are initially spherical, (2) Earth's obliquity is zero, (3) the Moon is a rigid body, and (4) the orbit is circular.</caption><thead><tr><th>Quantity</th><th>Unit</th><th>Earth</th><th>Moon</th></tr></thead><tbody><tr><td>m</td><td>[10<sup>24</sup> kg]</td><td>5.9724</td><td>0.07346</td></tr><tr><td>R</td><td>[km]</td><td>6371.0</td><td>1737.4</td></tr><tr><td>ξ</td><td></td><td>0.3308</td><td>0.394</td></tr><tr><td>k<sub>f</sub></td><td></td><td>0.933</td><td>0.0</td></tr><tr><td>τ</td><td>[s]</td><td>[2<sup>0</sup> - 2<sup>26</sup>]</td><td></td></tr><tr><td>P<sub>spin</sub></td><td>[day]</td><td>1</td><td></td></tr><tr><td>a</td><td>[km]</td><td></td><td>0.3844e6</td></tr><tr><td>e</td><td></td><td></td><td>0.0</td></tr></tbody></table>

<table><caption>**Table B2.** Initial conditions for the heartbeat star KOI-54 (Section 4.2), adopted from Welsh et al. (2011). The tidal response parameters are chosen to be in the linear tide regime, while the obliquites are set to zero.</caption><thead><tr><th>Quantity</th><th>Unit</th><th>Primary</th><th>Secondary</th></tr></thead><tbody><tr><td>m</td><td>[M<sub>☉</sub>]</td><td>2.39</td><td>2.33</td></tr><tr><td>R</td><td>[R<sub>☉</sub>]</td><td>2.33</td><td>2.20</td></tr><tr><td>ξ</td><td></td><td>0.052</td><td>0.052</td></tr><tr><td>k<sub>f</sub></td><td></td><td>0.012</td><td>0.012</td></tr><tr><td>τ</td><td>[yr]</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td></tr><tr><td>P<sub>spin</sub></td><td>[day]</td><td>2.55101914403</td><td>0.728862612579</td></tr><tr><td>a</td><td>[AU]</td><td></td><td>0.3956</td></tr><tr><td>e</td><td></td><td></td><td>0.8335</td></tr><tr><td>I</td><td>[deg]</td><td></td><td>5.50</td></tr><tr><td>Ω</td><td>[deg]</td><td></td><td>0</td></tr><tr><td>ω</td><td>[deg]</td><td></td><td>36.70</td></tr><tr><td>M</td><td>[deg]</td><td></td><td>180</td></tr></tbody></table>

<table><caption>**Table B3.** Initial conditions for exoplanet HD80606b (Section 4.3), which are reproduced from Correia et al. (2014).</caption><thead><tr><th>Quantity</th><th>Unit</th><th>Star</th><th>Planet</th></tr></thead><tbody><tr><td>m</td><td>[10<sup>30</sup> kg]</td><td>2.0088092</td><td>0.0077459434</td></tr><tr><td>R</td><td>[km]</td><td>702455.0</td><td>68488.3446</td></tr><tr><td>ξ</td><td></td><td>0.070</td><td>0.25</td></tr><tr><td>k<sub>f</sub></td><td></td><td>0</td><td>0.5</td></tr><tr><td>τ</td><td>[yr]</td><td>0</td><td>[10<sup>-4</sup> - 1]</td></tr><tr><td>P<sub>spin</sub></td><td>[day]</td><td>24.47</td><td>0.5</td></tr><tr><td>θ</td><td>[deg]</td><td>0</td><td>[0-180]</td></tr><tr><td>a</td><td>[AU]</td><td></td><td>0.455</td></tr><tr><td>e</td><td></td><td></td><td>0.9330</td></tr></tbody></table>

<table><caption>**Table B4.** Initial conditions for the Sun-Earth-Moon experiment in Section 4.4 adopted from Touma & Wisdom (1998) and Williams (2016). Earth's obliquity is zero with respect to the Earth-Moon orbit.</caption><thead><tr><th>Quantity</th><th>Unit</th><th>Earth</th><th>Moon</th><th>Sun</th></tr></thead><tbody><tr><td>m</td><td>[10<sup>24</sup> kg]</td><td>5.9724</td><td>0.07346</td><td>1989000.</td></tr><tr><td>R</td><td>[km]</td><td>6371.0</td><td>1737.4</td><td>695700</td></tr><tr><td>ξ</td><td></td><td>0.3308</td><td>0.394</td><td>0.070</td></tr><tr><td>k<sub>f</sub></td><td></td><td>0.933</td><td>0</td><td>0</td></tr><tr><td>τ</td><td>[s]</td><td>1</td><td></td><td></td></tr><tr><td>P<sub>spin</sub></td><td>[hour]</td><td>5</td><td></td><td></td></tr><tr><td>a</td><td>[km]</td><td></td><td>22298.5</td><td>149598073</td></tr><tr><td>e</td><td></td><td></td><td>0.01</td><td>0</td></tr><tr><td>I</td><td>[deg]</td><td></td><td>0</td><td>10</td></tr></tbody></table>

This paper has been typeset from a TeX/LaTeX file prepared by the author.

© 2023 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.