Numerical simulation of tidal evolution of a viscoelastic body modelled with a mass-spring network

Frouard, Julien; Quillen, Alice C.; Efroimsky, Michael; Giannella, David

doi:10.1093/mnras/stw491

Abstract

We use a damped mass-spring model within an N-body code to simulate the tidal evolution of the spin and orbit of a self-gravitating viscoelastic spherical body moving around a point-mass perturber. The damped mass-spring model represents a Kelvin–Voigt viscoelastic solid. We measure the tidal quality function (the dynamical Love number k₂ divided by the tidal quality factor Q) from the numerically computed tidal drift of the semimajor axis of the binary. The shape of k₂/Q, as a function of the principal tidal frequency, reproduces the kink shape predicted by Efroimsky for the tidal response of near-spherical homogeneous viscoelastic rotators. We demonstrate that we can directly simulate the tidal evolution of spinning viscoelastic objects. In future, the mass-spring N-body model can be generalized to inhomogeneous and/or non-spherical bodies.

methods: numerical, planets and satellites: dynamical evolution and stability, planets and satellites: interiors

1 MOTIVATION AND PLAN

The analytical theory of tidal interaction of solid bodies has a long and rich history – from the early mathematical development by Darwin (1879),¹ to its generalization by Kaula (1964), to an avalanche of more recent results. Verification of tidal theories through direct measurements is not easy because the tidal evolution is slow and requires either high observational precision (Williams & Boggs 2015) or an extended observational time span (Lainey et al. 2012). Purely analytical theories of tidal evolution describe homogeneous or, at best, two-layered (as in Remus et al. 2015) near-spherical bodies of linear rheology. This makes numerical simulations attractive (e.g. Henning & Hurtford 2014) as they may make it possible to explore the tidal evolution of more complex objects.

Here we explore computations that treat the solid medium as a set of mutually gravitating massive particles connected with a network of damped massless springs. Because of their simplicity and speed (compared to more computationally intensive grid-based or finite element methods), mass-spring computations are a popular method for simulating soft-deformable bodies (e.g. Nealen et al. 2006). Ostoja-Starzewski (2002) and Kot, Nagahashi & Szymczak (2015) have shown that mass-spring systems can accurately model elastic materials. As we shall show, viscoelastic response and gravitational forces can be incorporated into the mass-spring particle-based simulation technique. We demonstrate numerical simulations of the tidal response of a spinning homogeneous spherical body exhibiting spin-down or spin-up and associated drift in semimajor axis, without using an analytical tidal evolution model. Then we compare the outcome of our numerics with analytical calculations. The results are similar, justifying our numerical approach. Our numerical model may in future be used to take into account more complex effects that cannot be easily computed by analytical means (like inhomogeneity, compressibility, or a complex shape of the tidally perturbed body).

2 TIDES IN A KELVIN–VOIGT VISCOELASTIC SOLID

2.1 How tides work

Consider an extended spherical body of mass M and radius R, tidally deformed by a perturber of mass M* residing at a position |${{\boldsymbol {r}}}$|⁠, where |$|{\boldsymbol {r}}|\ge R$|⁠. The binary's orbit has the semimajor axis a and the mean motion |$n = \sqrt{G(M + M^{\ast })/a^3}$|⁠, where G is the gravitational constant. For a distant perturber (a ≫ R), the quadrupole part in the Fourier expansion (A5) for the perturbing potential W is dominant. This part comprises several terms. Of these, the term called semidiurnal is usually leading.² This term is a function of the principal tidal Fourier mode ω₂₂₀₀ which we denote simply as ω:

\begin{eqnarray} \omega \ \equiv \ \omega _{2200}\,=\ 2(n\ -\ \dot{\theta }), \end{eqnarray}

(1)

where θ and |$\dot{\theta }$| are the body's rotation angle and spin rate in the equatorial plane. (See the equation A14 in Appendix A.)

The secular part of the semidiurnal term of the polar tidal torque acting on the body is

\begin{eqnarray} \left\langle {\cal {T}}^{(z)}_{2200}\right\rangle \;=\;\frac{3}{2}\,G\,M^{{\ast }2}\,\frac{R^5}{a^6}\;k_2(\omega )\;\sin \epsilon _2(\omega ), \end{eqnarray}

(2)

where k₂(ω) and ϵ₂(ω) are the quadrupole dynamical Love number and the quadrupole phase lag, both taken at the semidiurnal frequency given by the above expression (1). The product k₂(ω) sin ϵ₂(ω) is often called the quality function (Makarov 2012, 2013; Efroimsky 2015) or, sometimes, kvalitet (Makarov 2015; Makarov, Frouard & Dorland 2016).

When the inclination and eccentricity are small, conservation of angular momentum gives an estimate of the secular drift rate of the semimajor axis (see the equation A17 in Appendix A):

\begin{eqnarray} \frac{\dot{a}}{n\,a} &=&\ -\frac{2\,{\cal {T}}^{(z)}_{2200}\,a}{G M^{\ast } M} \nonumber \\ &=& - 3 \left( \frac{M^{\ast }}{M} \right)\left( \frac{R}{a} \right)^5\;k_2(\omega )\;\sin \epsilon _2(\omega ) .\quad \end{eqnarray}

(3)

Below we compute the drift rate of the semimajor axis |$\dot{a}$| through a direct numerical simulation, and then compare the result with that obtained analytically from the tidal theory using a homogenous Kelvin–Voigt viscoelastic rheology. This comparison will demonstrate that we can directly simulate the tidal evolution of viscoelastic objects with a mass-spring model.

2.2 The shape of the quality function

As is demonstrated in Appendix A, the Fourier decompositions of the disturbing potential W, the tidal response potential U, and the tidal torque |${\cal {T}}$| comprise terms that are numbered with the four indices lmpq. An lmpq term of the torque is proportional to the quality function k_l(ω_lmpq) sin ϵ_l(ω_lmpq). The quality function depends on the degree l, the composition of the body, the rheology of its layers, and the size and mass of the body. The size and mass are important, because the tidal response is defined not only by the internal structure and rheology, but also by self-gravitation.

The process of deriving the quality function for any linear viscoelastic rheology is described in Efroimsky (2015). Here we provide a short account of that derivation. We begin with the expression for the static Love number for a homogeneous, incompressible, self-gravitating elastic sphere:

\begin{equation} k_l^{({\rm static})}\,=\;\frac{3}{2(l-1)}\;\frac{1}{1+B_l/J_{\rm r}}, \end{equation}

(4)

where

\begin{equation} B_l\,=\;\frac{3(2 l^2 + 4l + 3)}{4 l \pi G \rho ^2 R^2}\;=\;\frac{1}{e_{\rm g}}\;\frac{4\pi (2l^2+ 4l + 3)}{3l}, \end{equation}

(5)

while

\begin{equation} e_{\rm g}\,\equiv \;\frac{GM^2}{R^4} \end{equation}

(6)

is (to order of magnitude) the gravitational energy density of the body. Here J_r = 1/μ_r is the static (relaxed) compliance of the material, which is inverse to the static (relaxed) rigidity μ_r. Switching from a static to an evolving configuration, we invoke the equivalence principle for viscoelastic materials, in order to obtain the complex Love number in the frequency domain:

\begin{equation} \bar{k}_l(\chi ) = \frac{3}{2(l-1)} \frac{1}{1+B_l/\bar{J}} = |\bar{k}_l(\chi )| {\rm e}^{-{\rm i} \epsilon _l(\chi )}, \end{equation}

(7)

χ ≡ |ω| being the tidal frequency, and |$\bar{J}(\chi )$| being the complex compliance of the material. Once the compliance |$\bar{J}(\chi )$| is prescribed by a rheological model, we can compute the quantity function

\begin{equation} {k}_l(\chi )\;\sin \epsilon _l(\chi )\;=\;-\;{\rm Im} [\bar{k}_l(\chi )], \end{equation}

(8)

where |${k}_l(\chi )\,\equiv \,|\,\bar{k}_l(\chi )\,|$| and ϵ_l(χ) are, correspondingly, the degree l dynamical Love number and phase lag (the latter being linked to the degree l tidal quality factor through the equation A13).

However, an lmpq term in the expansion (A10)–(A11) for the tidal torque includes the factor k_l sin ϵ_l written not as a function of the tidal frequency χ ≡ |ω| but as a function of the tidal mode ω – see e.g. the semidiurnal term given by the expression (2). It can be demonstrated that this brings an extra factor equal to the sign of the mode:

\begin{eqnarray} &&{[k_l \sin \epsilon _l] (\omega )}\nonumber \\ &&{\quad=-\;\frac{3}{2(l-1)}\; \frac{B_l \; {\rm Im}({\bar{J}}(\chi ))}{[\, {\rm Re} (\,{\bar{J}} (\chi )\,) + B_l \,]^2 + [\,{\rm Im} (\,{\bar{J}}(\chi )\,)\,]^2 } {\rm Sign}\,\omega .} \end{eqnarray}

(9)

Whatever realistic rheological compliance |$\bar{J}(\chi )$| is inserted into the above formula, the shape of the quality function is similar for all viscoelastic bodies. It exhibits a sharp kink with two peaks having opposite signs (e.g. Noyelles et al. 2014; Efroimsky 2015).³

The generic shape of the quality function can be understood by comparing the frequency χ to the inverse of the viscoelastic relaxation time (e.g. Ferraz-Mello 2015b.). At a fixed point in the body, the tidal stressing in the material is oscillating at the frequency χ ≡ |ω|. When χ is small compared to the inverse time-scale of viscoelastic relaxation in the material, the body deformation stays almost exactly in phase with the tidal perturbation. The reaction and action being virtually in phase, no work is carried out and the tidal effects are minimal. On the other hand, at very high frequencies, the body's viscosity prevents it from deforming during the short forcing period 2π/χ. The reaction cannot catch up with the action, and stays close to zero – so, once again, little work is being done, and the tidal effects are again minimal.

2.3 The quality function for a Kelvin–Voigt sphere

The goal of our paper is not to favour a particular rheological model, but to test our simulation method. The Kelvin–Voigt rheological model is simplistic,⁴ but it is the easiest to model with a mass-spring model. Subsequent work will be aimed at extending our simulation approach to more realistic rheologies.

The Kelvin–Voigt model of a viscoelastic solid can be represented by a purely viscous damper and a purely elastic spring connected to two mass elements in parallel. If we connect these two elements in series rather than in parallel, we describe the Maxwell model; Fig. 1. A Kelvin–Voigt body is easier to model with a mass-spring system, because both the spring and damping forces are directly applied to each node particle. This can be seen from the expression for the stress tensor σ_pq as a function of the strain rate tensor |$\dot{\varepsilon }_{pq}$| in the time domain (see Efroimsky 2012a):

\begin{equation} \sigma _{pq}(t) = 2 \int _{-\infty }^t \mu (t-t^{\prime }) \dot{\varepsilon }_{pq}(t^{\prime }) {\rm d} t^{\prime }, \end{equation}

(10)

where μ(t − t′) is the stress-relaxation function. In the context of a mass-spring model, where the force between particles is likened to a uniaxial stress, the normal force applied to the particle i due to a spring and dashpot connecting it to particle j is given by

\begin{equation} F_i(t) = 2 \int _{-\infty }^t \mu (t-t^{\prime }) \dot{\varepsilon }(t^{\prime }) {\rm d} t^{\prime }. \end{equation}

(11)

For the Kelvin–Voigt model (e.g. Mase, Smelser & Mase 2010),

\begin{equation} \mu (t-t^{\prime }) = \mu + \eta \delta (t-t^{\prime }), \end{equation}

(12)

where μ and η are the unrelaxed shear rigidity and viscosity of the link between i and j. Inserting that rheology into the equation (11), we find

\begin{equation} F_i(t) = 2 \mu \varepsilon (t) - 2 \mu \varepsilon (-\infty ) + 2 \eta \dot{\varepsilon }(t). \end{equation}

(13)

In neglect of the strain at t = −∞, this becomes equivalent to equations (26)–(28) presented below and employed in our code.

Figure 1.

Green circles are mass elements. In the Kelvin–Voigt model (the top drawing), spring and damping elements are set in parallel. In the Maxwell model (the bottom drawing), the two elements are arranged in series. Since it is easier to represent a Kelvin–Voigt viscoelastic material with a mass-spring network, we compute the quality function for the Kelvin–Voigt model, thus allowing a direct comparison between the quality functions calculated from theory and that computed through simulation.

Open in new tab Download slide

In mass-spring model simulations, massive particles are interlinked with a network of massless springs. To each of the two particles linked by a spring, a damping force is applied that depends on the spring strain rate (see Escribano et al. (2008) with regards to the spin-orbit problem, and section 2.2 in Quillen et al. 2015). The shear elastic modulus, μ_I, can be computed for an isotropic, initially random mass-spring model from the strength, lengths, and distribution of the springs (Kot et al. 2015). We propose in the next section an estimate for the simulated material shear viscosity η_I based on the spring damping forces. Although the relation between the Kelvin–Voigt behaviour of the damped springs, described above, and the bulk and shear properties of the material (i.e. the constitutive equations) is not obvious, our computations demonstrate that the simulated resolved body behaves similar to that expected for a Kelvin–Voigt solid with the actual elastic modulus μ ≈ μ_I, viscosity η ≈ η_I, and relaxation time

\begin{equation} \tau \;=\;\frac{\eta }{\mu }. \end{equation}

(14)

There may be a difference between our a priori estimated values μ_I, η_I of the shear rigidity and viscosity (computed from the network), and the values μ, η of the simulated material (see the discussion by Kot et al. 2015).

We now derive the quality function (9) for the case of a Kelvin–Voigt sphere. We insert into the equation (9) the complex compliance of a Kelvin–Voigt body:

\begin{eqnarray} {\bar{J}}(\chi ,\,\mu ,\,\eta )\;\equiv \;\frac{\mu - {\rm i} \chi \eta }{\mu ^2 + \chi ^2 \eta ^2 }\;=\;\mu ^{-1}\,\frac{1 - {\rm i} \chi \tau }{1 + \chi ^2 \tau ^2}. \end{eqnarray}

(15)

The complex compliance |$\bar{J}$| can be presented either as a function of the unrelaxed shear modulus μ and viscosity η, or as a function of the shear modulus μ and the relaxation time τ. Introducing the dimensionless Fourier tidal mode

\begin{eqnarray} \tilde{\omega }\;\equiv \;\omega \;\tau , \end{eqnarray}

(16)

and following χ = |ω|, the dimensionless physical frequency

\begin{eqnarray} \tilde{\chi }\;\equiv \;\chi \;\tau , \end{eqnarray}

(17)

we write down the complex compliance as

\begin{eqnarray} {\bar{J}}(\tilde{\chi },\,\mu )\;\equiv \;\mu ^{-1}\;\frac{1 - {\rm i} \tilde{\chi } }{1 + \tilde{\chi }^2 }. \end{eqnarray}

(18)

We define a dimensionless function

\begin{eqnarray} {\bar{j}}^{(l)}({\tilde{\chi }},\mu ) & \equiv & B_l^{-1} {\bar{J}}(\tilde{\chi }, \mu )\;=\;(B_l\,\mu )^{-1} \frac{1 - {\rm i} \tilde{\chi }}{1 + \tilde{\chi }^2 }, \end{eqnarray}

(19)

with B_l given by the expression (5). Using that expression, we write down the function of the degree l = 2 as

\begin{eqnarray} {\bar{j}}^{(l=2)}(\tilde{\chi },\,\mu ) &=& \frac{3}{38\pi } \;\frac{e_{\rm g}}{\mu }\;\frac{1\,-\,{\rm i}\,\tilde{\chi }}{1\,+\,{\tilde{\chi }}^2 }. \end{eqnarray}

(20)

This is a function of the shear modulus μ given in the units of e_g, and of the frequency χ given in the units of the inverse relaxation time τ. In terms of the dimensionless Fourier mode |$\tilde{\omega }$| and the dimensionless frequency |$\tilde{\chi }\equiv |\tilde{\omega }|$|⁠, equation (9) becomes

\begin{eqnarray} &&{[\,k_l \;\sin \epsilon _l\,] (\tilde{\omega },\;\mu )}\nonumber \\ &&{\,=\,-\frac{3}{2(l-1)} \frac{{\rm Im}(\,{\bar{j}}^{(l)}(\tilde{\chi })\,)}{\left[\,{\rm Re} (\,{\bar{j}}^{(l)} (\tilde{\chi })\,) + 1\,\right]^2 + \left[\,{\rm Im} (\,{\bar{j}}^{(l)}(\tilde{\chi })\,)\,\right]^2 } {\rm Sign} \, \tilde{\omega }.} \end{eqnarray}

(21)

This function attains its extrema at

\begin{eqnarray*} \omega _{{\rm peak}}^{(l)}\,=\;\pm \;\frac{\mu \;B_l\,+\;1}{B_l\;\eta }\; =\;\pm \;\frac{1}{\tau }\;\left( 1\;+\;\frac{1}{\mu \;B_l} \right). \end{eqnarray*}

For l = 2, the dimensionless peak mode and frequency are

\begin{eqnarray} \tilde{\omega }_{{\rm peak}}^{(l=2)}\,\equiv \,\omega _{{\rm peak}}^{(l=2)}\,\tau \,=\,\pm \,\left( 1\,+\, \frac{3}{38 \pi } \;\frac{e_{\rm g}}{\mu } \right)\qquad \end{eqnarray}

(22a)

and

\begin{eqnarray} \tilde{\chi }_{{\rm peak}}^{(l=2)}\,\equiv \,\chi _{{\rm peak}}^{(l=2)}\,\tau \,=\,\left|\omega _{{\rm peak}}^{(l=2)}\right|\,\tau \,=\,\left(1\,+\, \frac{3}{38 \pi } \;\frac{e_{\rm g}}{\mu }\right).\ \ \end{eqnarray}

(22b)

Elastic bodies should obey μ > e_g, lest they collapse due to self-gravity. Hence we expect that

\begin{eqnarray} \omega _{{\rm peak}}^{(l=2)}\;\tau \;\approx \;\pm \;1. \end{eqnarray}

(23)

Using the shorthand notation |$\tilde{\omega }\,=\,\tilde{\omega }^{(l=2)}$| and

\begin{eqnarray} y(\mu , \tilde{\omega })\;\equiv \; \frac{3}{38 \pi } \; \frac{e_{\rm g}}{\mu }\;(1\;+\;\tilde{\omega }^2)^{-1}, \end{eqnarray}

(24)

the equation (21) for l = 2 can be written as

\begin{eqnarray} [k_2 \sin \epsilon _2] (\tilde{\omega },\;\mu )\;=\;\frac{3}{2}\;\frac{y\;\tilde{\omega }}{y^2\;(1 + \tilde{\omega }^2) + 2 y + 1}. \end{eqnarray}

(25)

We will use this expression for the quality function to analytical predict tidal response.

Classical tidal theory is valid only for incompressible materials (those having the Poisson ratio ν = 0.5) – which is why the standard expression for the static Love number k₂ of an incompressible sphere depends only on the shear modulus of rigidity, not on the bulk modulus. However, the mass-spring models approximate a material with Poisson's ratio ν = 0.25 (Kot et al. 2015). However, Love (1911) derived a general formula for the static Love number k₂ also for a compressible homogeneous sphere.⁵ We have numerically checked that his formula for k₂ in the compressible case is only very weakly dependent on the Poisson's ratio, within broad ranges of the body size and density. This makes us confident that the standard tidal formulae (derived for the Poisson ratio ν = 0.5) can accurately describe a compressible body with Poisson ratio ν = 0.25.

3 DAMPED MASS-SPRING MODEL SIMULATIONS

We will compare the tidal spin-down rate of a simulated viscoelastic body to that predicted analytically. To simulate tidal viscoelastic response we use the mass-spring model by Quillen et al. (2015), that is based on the modular N-body code rebound (Rein & Liu 2012). A random spring network model, rather than a lattice network model, was chosen so that the modelled body is approximately isotropic and homogenous and lacks planes associated with crystalline symmetry. We work in units of radius and mass of the tidally perturbed body: R = 1, M = 1. Time is given in units of |$t_{{\rm grav}} = \sqrt{R^3/GM\,}$| referred to as the gravitational time-scale. Pressure, energy density, and elastic moduli are specified in units of e_g ≡ GM²/R⁴. In these units, the velocity of a massless particle on a circular orbit, which is just grazing the surface of the body, is 1, the period of the grazing orbit being 2π.

Modelling the tidally perturbed body, we randomly generate an initial spherical distribution of particles of equal mass (as described in section 2.2 in Quillen et al. 2015); see our Fig. 2 for an illustration. Particle positions are randomly generated using a uniform distribution in three dimensions but a particle is added (as a node) to the spring network only if it is sufficiently separated from other particles (at a distance greater than minimum distance d_I) and within a radius of R = 1 from the body centre. Springs are added between two nodes if the distance between the nodes is less than distance d_s.

Figure 2.

A snapshot of one of our simulations as viewed in the open-GL viewer of rebound. We only show the tidally perturbed body, as the perturbing one is distant from it. The rendered spheres are shown to illustrate the random distribution of node point masses, not imply that the body behaves as a rubble pile (e.g. Richardson et al. 2009; Sánchez & Scheeres 2011). Nodes are connected with a network of damped elastic springs.

Open in new tab Download slide

The particles are subjected to three types of forces: the gravitational forces acting on every pair of particles in the body and with the massive companion, and the elastic and damping spring forces acting only between sufficiently close particle pairs. Springs with a spring constant k_I interconnect each pair of particles closer than some distance d_s. The rest length of each spring is initially set equal to its initial length. The springs experience compression or extension relative to their rest length. The number of springs and their rest lengths stay fixed during our simulations. The springs have different rest lengths, and we denote their average initial (rest) length with L_I. The total number of particles is N_I, and the total number of springs is NS_I. The mass of each particle is m_I = 1.0/N_I, and the initial mass density is approximately uniform.

The number and distribution of particles and links in a mass-spring model is mainly set by the minimum distances d_I, d_s, however, because the initial particle positions are randomly generated, there are local variations in density of nodes and the spring network. The ratio d_s/d_I sets the mean number of springs per mass node.

Consider a spring linking the particle i residing at |${\boldsymbol x}_i$| with a particle j located in |${\boldsymbol x}_j$|⁠. The force due to the spring on each mass node is computed as follows. The vector pointing from one of these particles to another, |${\boldsymbol x}_i - {\boldsymbol x}_j$|⁠, gives the spring length |$L_{ij} = |{\boldsymbol x}_i - {\boldsymbol x}_j|$| that we compare with the spring rest length L_{ij, 0}. The elastic force exerted upon particle i by the particle j is

\begin{equation} {\boldsymbol F}_{i}^{{\rm elastic}} = -k_{\rm I} (L_{ij} - L_{ij,0}) \hat{\boldsymbol n}_{ij}, \end{equation}

(26)

where k_I is the spring constant, while the unit vector is |$\hat{\boldsymbol n}_{ij} =({\boldsymbol x}_i - {\boldsymbol x}_j)/ L_{ij}$|⁠. The appropriate force acting on the particle j from the particle i is equal in magnitude and opposite in direction.

The strain rate of a spring with length L_ij is

\begin{equation} \dot{\epsilon }_{ij} = \frac{\dot{L}_{ij}}{L_{ij,0}} = \frac{1}{L_{ij} L_{ij,0}} ({\boldsymbol x}_i - {\boldsymbol x}_j) \cdot ({\boldsymbol v}_i - {\boldsymbol v}_j), \end{equation}

(27)

where |${\boldsymbol v}_i$| and |${\boldsymbol v}_j$| are the particle velocities and |$\dot{L}_{ij}$| is the rate of change of the spring length. To the elastic force acting on the particle i, we add a damping (viscous) force proportional to the strain rate:

\begin{equation} {\boldsymbol F}_{i}^{{\rm damping}} = - \gamma _{\rm I} \dot{\epsilon }_{ij} L_{ij,0} m_{\rm I} \hat{{\boldsymbol n}}_{ij}, \end{equation}

(28)

with a damping coefficient γ_I equal to the inverse damping time-scale. The parameter γ_I is independent of k_I. As all our particles have the same mass, we do not use the reduced mass in equation (28), as did Quillen et al. (2015).

How does the spring constant k_I and the damping parameter γ_I relate to the global rigidity and viscosity of the body? According to Kot et al. (2015), the static Young's modulus is given by a sum over the springs L_{ij, 0}:

\begin{equation} E_{\rm I}\,=\;\frac{1}{6V}\;\sum k_{\rm I} L_{ij,0}^2, \end{equation}

(29)

where V is the total volume. For an initially random isotropic mass-spring system, the Poisson ratio is ν = 0.25 (Kot et al. 2015). Our mass-spring model allows us to directly estimate the Young's modulus E_I, from which we can compute the shear elastic rigidity μ_I commonly used for tidal evolution calculations. The relation between the two relaxed (static) moduli is given by

\begin{equation} \mu _{\rm I} = \frac{E_{\rm I}}{2(1+\nu )} = \frac{E_{\rm I}}{2.5}, \end{equation}

(30)

where we set the Poisson ratio to be ν = 0.25.

With non-zero damping coefficients, the stress is a sum of an elastic term proportional to the strain and a viscous term proportional to the strain rate, so the model should locally approximate a linear Kelvin–Voigt rheology. The mass-spring model is compressible. So, when damped springs are used, it exhibits a bulk viscosity (in analogy to the bulk modulus) and a shear viscosity (in analogy to the shear modulus). For a mass-spring model composed of equal masses m_I and parametrized with the damping coefficients γ_I and spring constants k_I, we can tentatively estimate the shear viscosity η_I as the ratio of the damping and elastic forces given by the equations (28) and (26), correspondingly. We also assume that η_I scales with the network in the same way μ_I does in equations (30) and (29). The resulting estimate for the viscosity is

\begin{equation} \eta _{\rm I} \approx \mu _{\rm I} \left( \frac{\gamma _{\rm I} m_{\rm I}}{k_{\rm I}}\right). \end{equation}

(31)

The relaxation time-scale of a Kelvin–Voigt solid is

\begin{equation} \tau _{{\rm relax}}\ =\;\frac{\eta _{\rm I}}{\mu _{\rm I}}\;\approx \;\frac{\gamma _{\rm I}\;m_{\rm I}}{k_{\rm I}}. \end{equation}

(32)

While the fidelity of the computed elastic modulus of a mass-spring model has been checked by numerical simulations using static applied forces (Kot et al. 2015), the viscosity of a damped mass-spring model has not been tested. We consider the equations (31) and (32) as approximate. They may need to be amended with factors of order unity. We shall discuss this possibility later, when we compare measurements from our simulations to predictions from an analytical model of tidal evolution.

In our simulations, the body was given an initial spin |$\dot{\theta }=\sigma _0$| perpendicular to the orbital plane. This was done by setting the initial velocities of particle equal to

\begin{equation} {\boldsymbol v}_i={\boldsymbol x}_i\times \sigma _0\,\hat{\boldsymbol z}, \end{equation}

(33)

|${\boldsymbol v}_i$| and |${\boldsymbol x}_i$| being the velocity and position vectors of the ith particle, with respect to the body's centre of mass, and the unit vector |$\hat{\boldsymbol z}$| being orthogonal to the orbit. In most runs, we chose the initial rotation to be retrograde (σ₀ < 0), to allow for a larger range of values of the tidal frequency χ to be simulated. This choice always led to a decrease in the semimajor axis (see Table 3) and to acceleration of the body's spin rate. Variations in the initial particle distribution create only negligible non-diagonal terms in the inertia tensor in this coordinate system.

From the rebound code, version 2 as of 2015 November, we used the open-GL display with open boundary conditions, the direct all-pairs gravitational force computation, and the leap-frog integrator needed to advance particle positions. To the particles’ accelerations caused by the gravity, we added the additional spring and spring-damping forces, as was explained above. To maintain numerical stability, the time step was chosen to be smaller than the elastic oscillation frequency of a single node particle in the spring network, or equivalently the time it takes vibrational waves to travel between two neighbouring nodes.

4 TIDAL EVOLUTION OF THE SEMIMAJOR AXIS

4.1 The semimajor axis’ tidal drift rate computed from the mass-spring model

Common parameters for our first set of simulations are listed in Table 1, along with their chosen or computed values. A list of varied and measured parameters is presented in Table 2. The values used in different individual simulations with common values in Table 1 are listed in Table 3. We chose the values for the mass ratio M*/M and the initial semimajor axis a₀ to be the same in the two sets of simulations performed. However, in each set, the spring damping rate γ_I and the initial body spin rate σ₀ were chosen to sample a range of values of the Fourier tidal mode frequency ω and of the viscoelastic relaxation time τ_relax. Using the equation (29), we computed the Young's modulus by only considering the springs with a midpoint radii less than 0.9R, and we used the volume within the same radius 0.9R. The region near the surface was discarded so that the Young's modulus was computed in a region where the spring network is isotropic. A fairly soft body under strong tidal forcing (corresponding to a large mass ratio and weak springs) was chosen, to reduce the integration time required to observe a significant tidally induced change in the semimajor axis. At the same time, we made sure that the Young's modulus was sufficiently large, so that the body was strong enough to maintain a nearly constant radial density profile. In the absence of exterior pressure, the body is held up against self-gravity by spring forces only – so the springs in the interior are under compression. In all our runs, the particles were displaced, in the body frame, by at most a few per cent of the unit length R.

Table 1.

Open in new tab

Common simulation parameters.

N_I

800

Number of particles in resolved body

NS_I

9254

Number of interconnecting springs

L_I

0.2623

Mean rest spring length

k_I

0.08

Spring constant

E_I

2.3

Mean Young's modulus

d_I

0.15

Minimum initial interparticle distance

d_s

0.345

Spring formation distance

dt

0.001

Time step

Notes. N_I, NS_I, E_I, and L_I vary slightly between simulations as particle distributions are randomly generated. E_I is computed using equation (29) and is the average value for all the simulations. These parameters are common to simulations listed in Table 3.

Table 2.

Open in new tab

Description of varied and measured simulation parameters.

a₀

Initial semimajor axis

M*/M

Mass ratio

|$\tilde{\chi }$|

Unitless frequency

|$\dot{a}$|

Rate of orbital decay in semimajor axis

γ_I

Spring relaxation time

σ₀

Initial spin

τ_relax

Estimated relaxation time of viscoelastic solid

χ

Initial tidal frequency (semidiurnal)

P_o

Time between integration outputs

E_I

Computed Young's modulus

Notes. The viscoelastic relaxation time τ_relax is computed using the equations (29) and (32). The frequency χ is defined by the expression (A14), while the tidal forcing frequency |$\tilde{\chi }$| is given by the equation (17). The period P_o = 2π/χ is also a part of the simulation output.

Table 3.

Open in new tab

Quantities either set or computed in different mass-spring N-body simulations.

\|$\tilde{\chi }$\|	\|$\dot{a}$\|	γ_I	σ₀	τ_relax	χ	P_o	E_I
With a₀ = 10, M* = 100, n = 0.318
0.069	−3.897e-05	20	0.1	0.158	0.436	14.424	2.30
0.162	−8.505e-05	20	−0.2	0.156	1.036	6.067	2.35
0.223	−1.055e-04	20	−0.4	0.156	1.436	4.377	2.40
0.447	−1.800e-04	40	−0.4	0.311	1.436	4.377	2.38
0.645	−2.309e-04	50	−0.5	0.394	1.636	3.841	2.33
0.896	−2.517e-04	80	−0.4	0.624	1.436	4.377	2.36
1.123	−2.288e-04	100	−0.4	0.782	1.436	4.377	2.36
1.467	−2.325e-04	100	−0.6	0.799	1.836	3.423	2.26
With a₀ = 20, M* = 200, n = 0.159
0.080	−2.239e-06	20	−0.1	0.156	0.517	12.153	2.40
0.160	−4.540e-06	40	−0.1	0.309	0.517	12.153	2.41
0.227	−7.124e-06	40	−0.2	0.317	0.717	8.763	2.29
0.350	−8.601e-06	40	−0.4	0.314	1.117	5.625	2.36
0.534	−1.157e-05	60	−0.4	0.478	1.117	5.625	2.25
0.701	−1.324e-05	80	−0.4	0.627	1.117	5.625	2.26
0.881	−1.656e-05	100	−0.4	0.789	1.117	5.625	2.28
1.188	−1.372e-05	100	−0.6	0.783	1.517	4.142	2.31
1.771	−1.197e-05	150	−0.6	1.167	1.517	4.142	2.28

\|$\tilde{\chi }$\|	\|$\dot{a}$\|	γ_I	σ₀	τ_relax	χ	P_o	E_I
With a₀ = 10, M* = 100, n = 0.318
0.069	−3.897e-05	20	0.1	0.158	0.436	14.424	2.30
0.162	−8.505e-05	20	−0.2	0.156	1.036	6.067	2.35
0.223	−1.055e-04	20	−0.4	0.156	1.436	4.377	2.40
0.447	−1.800e-04	40	−0.4	0.311	1.436	4.377	2.38
0.645	−2.309e-04	50	−0.5	0.394	1.636	3.841	2.33
0.896	−2.517e-04	80	−0.4	0.624	1.436	4.377	2.36
1.123	−2.288e-04	100	−0.4	0.782	1.436	4.377	2.36
1.467	−2.325e-04	100	−0.6	0.799	1.836	3.423	2.26
With a₀ = 20, M* = 200, n = 0.159
0.080	−2.239e-06	20	−0.1	0.156	0.517	12.153	2.40
0.160	−4.540e-06	40	−0.1	0.309	0.517	12.153	2.41
0.227	−7.124e-06	40	−0.2	0.317	0.717	8.763	2.29
0.350	−8.601e-06	40	−0.4	0.314	1.117	5.625	2.36
0.534	−1.157e-05	60	−0.4	0.478	1.117	5.625	2.25
0.701	−1.324e-05	80	−0.4	0.627	1.117	5.625	2.26
0.881	−1.656e-05	100	−0.4	0.789	1.117	5.625	2.28
1.188	−1.372e-05	100	−0.6	0.783	1.517	4.142	2.31
1.771	−1.197e-05	150	−0.6	1.167	1.517	4.142	2.28

Notes. The rightmost column contains the values of the Young's modulus, computed for each simulation by means of the equation (29). The simulations listed here have common parameters listed in Table 1.

Table 3.

Open in new tab

Quantities either set or computed in different mass-spring N-body simulations.

\|$\tilde{\chi }$\|	\|$\dot{a}$\|	γ_I	σ₀	τ_relax	χ	P_o	E_I
With a₀ = 10, M* = 100, n = 0.318
0.069	−3.897e-05	20	0.1	0.158	0.436	14.424	2.30
0.162	−8.505e-05	20	−0.2	0.156	1.036	6.067	2.35
0.223	−1.055e-04	20	−0.4	0.156	1.436	4.377	2.40
0.447	−1.800e-04	40	−0.4	0.311	1.436	4.377	2.38
0.645	−2.309e-04	50	−0.5	0.394	1.636	3.841	2.33
0.896	−2.517e-04	80	−0.4	0.624	1.436	4.377	2.36
1.123	−2.288e-04	100	−0.4	0.782	1.436	4.377	2.36
1.467	−2.325e-04	100	−0.6	0.799	1.836	3.423	2.26
With a₀ = 20, M* = 200, n = 0.159
0.080	−2.239e-06	20	−0.1	0.156	0.517	12.153	2.40
0.160	−4.540e-06	40	−0.1	0.309	0.517	12.153	2.41
0.227	−7.124e-06	40	−0.2	0.317	0.717	8.763	2.29
0.350	−8.601e-06	40	−0.4	0.314	1.117	5.625	2.36
0.534	−1.157e-05	60	−0.4	0.478	1.117	5.625	2.25
0.701	−1.324e-05	80	−0.4	0.627	1.117	5.625	2.26
0.881	−1.656e-05	100	−0.4	0.789	1.117	5.625	2.28
1.188	−1.372e-05	100	−0.6	0.783	1.517	4.142	2.31
1.771	−1.197e-05	150	−0.6	1.167	1.517	4.142	2.28

\|$\tilde{\chi }$\|	\|$\dot{a}$\|	γ_I	σ₀	τ_relax	χ	P_o	E_I
With a₀ = 10, M* = 100, n = 0.318
0.069	−3.897e-05	20	0.1	0.158	0.436	14.424	2.30
0.162	−8.505e-05	20	−0.2	0.156	1.036	6.067	2.35
0.223	−1.055e-04	20	−0.4	0.156	1.436	4.377	2.40
0.447	−1.800e-04	40	−0.4	0.311	1.436	4.377	2.38
0.645	−2.309e-04	50	−0.5	0.394	1.636	3.841	2.33
0.896	−2.517e-04	80	−0.4	0.624	1.436	4.377	2.36
1.123	−2.288e-04	100	−0.4	0.782	1.436	4.377	2.36
1.467	−2.325e-04	100	−0.6	0.799	1.836	3.423	2.26
With a₀ = 20, M* = 200, n = 0.159
0.080	−2.239e-06	20	−0.1	0.156	0.517	12.153	2.40
0.160	−4.540e-06	40	−0.1	0.309	0.517	12.153	2.41
0.227	−7.124e-06	40	−0.2	0.317	0.717	8.763	2.29
0.350	−8.601e-06	40	−0.4	0.314	1.117	5.625	2.36
0.534	−1.157e-05	60	−0.4	0.478	1.117	5.625	2.25
0.701	−1.324e-05	80	−0.4	0.627	1.117	5.625	2.26
0.881	−1.656e-05	100	−0.4	0.789	1.117	5.625	2.28
1.188	−1.372e-05	100	−0.6	0.783	1.517	4.142	2.31
1.771	−1.197e-05	150	−0.6	1.167	1.517	4.142	2.28

Notes. The rightmost column contains the values of the Young's modulus, computed for each simulation by means of the equation (29). The simulations listed here have common parameters listed in Table 1.

We chose the initial semimajor axis a₀ large enough, to ensure that the quadrupole tidal potential term would dominate. We also set the initial relative velocity of the bodies such that the orbit would be circular. In the frame corotating with the tidally perturbed body, the perturber orbits with a period of P_o = 2 π/χ. We carried out each run over the time span of t = 11 P_o and recorded the semimajor axis’ values 10 times, with even intervals of P_o. Since the springs’ lengths had initially been set to have their rest values, gravity caused the system to bounce at the beginning of each simulation. Damped oscillations are expected in our model as each of the individual links between pairs of particles is acting as a damped harmonic oscillator. The initial oscillations are just the consequences of abruptly ‘turning on’ self-gravity at the beginning of the simulations. Thus, during the first time interval of P_o, we integrated the body with a higher damping parameter, to dissipate the initial vibrational oscillations. We did not take into account the semimajor axis’ value computed during that time interval. Subsequently, we recorded the semimajor axis’ values at an interval of P_o, so that the irregularities of the particle distribution did not affect the measurement of the slowly drifting semimajor axis.

The semimajor axis was computed using the distance between the centre of mass of the bodies, and their relative velocities. To measure the semimajor axis’ drift rate |$\dot{a}$|⁠, we fit a line to the 10 measurements of the semimajor axis at the 10 time intervals P_o. Each fit was individually inspected, to ensure that the 10 points lay on a line. The standard error of the fitted slope value that provides |$\dot{a}$| was ≲ 1 per cent. For each simulation, we also compared the initial and final values of each component of the total angular-momentum vector (relative to the centre of mass of the binary), and found the absolute difference was below 10⁻¹² for all the components. Similarly, we checked that the evolution of the measured spin rate and semimajor axis were tightly anticorrelated, as expected when angular momentum is conserved. The total energy was not conserved because of the spring damping forces.

For each simulation we generate a new mass-spring network. Hence each simulation has slightly different numbers of mass nodes and springs and variations in the node distribution and associated spring network. Two simulations run with identical input parameters will differ slightly in their measured semimajor axis drift rates. We set the minimum distance between nodes d_I and maximum spring length d_s (setting the mean number of mass nodes and numbers of springs) so that the differences in semimajor axis drift rate for simulations drawn from the same parameter set differed by less than 10 per cent. We will discuss the sensitivity of the simulations to the numbers of nodes and springs per node further below.

The simulations were run on a MacBook Pro (early 2015) with 3.1 GHz Intel Core i7 microprocessor. The computation time for an individual simulation listed in Tables 1 and 3 was about 7 min. If the number of mass nodes is doubled, then the number of direct all-pairs gravity force computations increases by a factor of 4. However the total computation time increases by a slightly larger factor than 4 because the number of springs also increases (scales with the number of nodes) and the time step decreases. For simulations with similar vibration wave speeds, the time step scales with the interparticle spacing and so the total number of particles to the −1/3 power.

4.2 Comparison of numerically measured and predicted quality functions

Inverting the equation (3), we obtain the following expression for the quality function:

\begin{equation} [k_2 \sin \epsilon _2](\tilde{\chi })\;= \;-\;\frac{\dot{a}}{3na} \left( \frac{M}{M^{\ast }} \right)\left( \frac{a}{R} \right)^5 . \end{equation}

(34)

From the semimajor axis’ drift rates computed numerically in our simulations, we computed the values of the quality function over a range of values of frequency. This was performed using the equation (34) and the quantities listed in Table 3. The numerically measured values of the quality function are plotted as a function of the dimensionless frequency |$\tilde{\chi }$| in Figs 3 and 4. The dimensionless frequency |$\bar{\chi }$| was computed by using the relaxation time that was evaluated for each simulation individually.

Figure 3.

A comparison of the numerically computed quality function, shown as points, to that calculated analytically for a homogeneous Kelvin–Voigt sphere. The analytically derived frequency dependence is given by the grey curve (obtained with equation 25), and when multiplied by s, in blue, with the scaling factor s listed on bottom right. The green curves show the effect of raising and lowering the shear modulus by 10 per cent in the offset analytical calculation. Both the numerical and analytical calculations were carried out for a perturber of mass M* = 100 and initial value of the semimajor axis a₀ = 10. We find that the numerically computed quality function has a shape and peak frequency consistent with that obtained analytically, but the amplitude is too high by about 30 per cent. We attribute the scatter of the points off the line to variations in the value of the shear modulus of the random mass-spring network due to non-uniformity in the particle distribution and spring network (see Section 4.2).

Open in new tab Download slide

The same as in Fig. 3, except that the perturber's mass is M* = 200, and the initial value of the semimajor axis is a0 = 100.

Figure 4.

The same as in Fig. 3, except that the perturber's mass is M* = 200, and the initial value of the semimajor axis is a₀ = 100.

Open in new tab Download slide

The numerically generated points in Fig. 3 show that the simulated quality function is proportional to the frequency, at small frequencies, but decays at large frequencies. Predicted analytically for various rheologies (Efroimsky 2012a, 2015; Noyelles et al. 2014), this behaviour has never been reproduced by direct numerical computations. The points used to build this plot were measured from simulations with different spin rates and relaxation times. Nevertheless, they all lie near the same curve, suggesting that the function is primarily dependent on the relaxation time. We thus have numerically confirmed the expected behaviour and sensitivity of the quality function to the frequency and to the viscoelastic time-scale.

With the points obtained through simulations, we plot the quality function predicted analytically using equation (25). We set the shear modulus μ = μ_I, with the value of μ_I given by the expression (30). This value is proportional to the Young's modulus listed in Table 1 computed from the spring network. The mass ratio and mean motion are taken from Table 3. The result is the grey line (the lowest one) in Fig. 3. This is the analytically predicted quality function, and we see that its numerically obtained counterpart (given by the red points) is higher. We refer to the ratio of numerically measured quality function to the predicted one as qf_ratio and from this plot we find qf_ratio ∼ 1.3.

Fig. 3 demonstrates that the numerically obtained tidal drift rate of the semimajor axis is maximal at a frequency |$\tilde{\chi } \approx 1.0$|⁠, which is consistent with the analytical model. Had the relaxation time-scale (equation 32) been miscomputed in our simulations, the numerically obtained peak would have been displaced away from that predicted analytically. The good match between predicted and numerically measured peak frequency supports our estimate for the shear viscosity (equation 31) and the associated relaxation time (equation 32) for the mass-spring model.

Both the peak magnitude and the peak location of the quality function may shift if higher order terms (with higher values of l) are included into the analytical calculation. As this might explain the difference in height of our numerically computed quality function, compared to that predicted, we tested this possibility by running simulations with a larger value of the semimajor axis. Simulations with a larger initial semimajor axis are also listed in Table 3, and the quality function for these runs is plotted in Fig. 4. We find that the amplitude correction factor required to match the numerical results is the same as for the previous set of simulations – and, again, the frequency does not need to be rescaled. From this, we conclude that higher order terms in the tidal potential do not explain the amplitude discrepancy between our numerical simulations model and our analytical predictions.

Since we are using a random mass-spring model, both the particle distribution and the spring network differ between simulations. We computed the Young's modulus and relaxation time for each run, and these are listed in Table 3. We found small variations in the elastic modulus between different simulations. In Figs 3 and 4, we also show the analytically predicted quality function factored by 1.3 and offset by ± 10 per cent (green curves). Points with higher values of E_I, corresponding to harder bodies, systematically lie lower than the appropriate points with lower values of E_I corresponding to softer bodies experiencing stronger tidal deformation. The scatter in our points above and below the blue line can be attributed to variations in the particle distribution and in the associated spring network.

4.3 Discrepancy in amplitude of quality function and sensitivity to the number of masses and springs simulated

To see how the discrepancy in amplitude is related to the number of mass nodes and springs, we carried out a second series of simulations each with the same initial spin, estimated elastic modulus, and viscoelastic relaxation time-scale, but having different numbers of node masses and numbers of springs per node. Parameters for these simulations are listed in Table 4. We ran six sets, each of five simulations, all approximately matching the third row in Table 3 with a₀ = 10, M_*/M = 100, σ₀ = −0.4, and P_o = 4.377. The five simulations in each set have identical run parameters. The first 10 simulations (columns A, B in Table 4) have about 400 mass nodes, the second 10 (columns C, D) about 800 mass nodes, and the last 10 (columns E, F) about 1600 mass nodes. The first five in each group of 10 have about 10 springs per node (columns A, C, and E) whereas the second five of each group of 10 (columns B, D, and F) have about 20 springs per node. The spring constants, k_I, and damping parameters, γ_I, for each set were chosen so that |$\bar{\chi }\approx 0.225$|⁠, E_I ≈ 2.3, and τ_relax ≈ 0.155. For each simulation separately we computed ratio qf_ratio of the numerically measured value of the quality function k₂sin ϵ₂ to the predicted one computed using the normalized frequency |$\bar{\chi }$| and Young's modulus measured from the nodes and springs in the simulations. We then computed the mean and standard deviation of qf_ratio for each set of five simulations with identical run parameters and these are listed in the bottom rows of Table 4. The standard deviation in semimajor axis drift rate divided by the mean value of the drift rate is also computed for each set of five simulations and listed in Table 4.

Table 4.

Open in new tab

Comparison of simulations with different numbers of node masses and springs per node.

	Set
	A	B	C	D	E	F
〈N_I〉	403.8	405.8	795.2	804.2	1597.8	1596.2
〈NS_I/N_I〉	10.92	18.7	11.5	20.34	12.02	21.34
d_I	0.19	0.19	0.15	0.15	0.118	0.118
d_s/d_I	2.3	2.8	2.2	2.8	2.3	2.8
k_I	0.102	0.04	0.08	0.0301	0.062	0.023
γ_I	12.858	5.05	20.0	7.48	31.3	11.5
〈qf_ratio〉	1.48	1.39	1.41	1.36	1.40	1.38
σ[qf_ratio]	0.069	0.052	0.053	0.026	0.030	0.012
\|$\sigma [\dot{a}]/ \langle \dot{a} \rangle$\|	0.061	0.064	0.055	0.024	0.026	0.0094

	Set
	A	B	C	D	E	F
〈N_I〉	403.8	405.8	795.2	804.2	1597.8	1596.2
〈NS_I/N_I〉	10.92	18.7	11.5	20.34	12.02	21.34
d_I	0.19	0.19	0.15	0.15	0.118	0.118
d_s/d_I	2.3	2.8	2.2	2.8	2.3	2.8
k_I	0.102	0.04	0.08	0.0301	0.062	0.023
γ_I	12.858	5.05	20.0	7.48	31.3	11.5
〈qf_ratio〉	1.48	1.39	1.41	1.36	1.40	1.38
σ[qf_ratio]	0.069	0.052	0.053	0.026	0.030	0.012
\|$\sigma [\dot{a}]/ \langle \dot{a} \rangle$\|	0.061	0.064	0.055	0.024	0.026	0.0094

Notes. With a₀ = 10, M* = 100, n = 0.318, σ₀ = −0.4, P_o = 4.377, |$\bar{\chi }\sim 0.225$|⁠, E_I ∼ 2.3, and τ_relax ∼ 0.155.

Each column represents a group of five simulations. From each group of five the mean number of nodes and springs per node are listed in the second and third rows. The rows labelled 〈qf_ratio〉 and σ[qf_ratio] give the mean and standard deviations, respectively, of the ratio of the numerically measured to analytically predicted quality function. The bottom row gives the ratio of the standard deviation in the semimajor axis drift rate divided by the means. Each of the means and standard deviations are computed from five simulations. The simulations listed here differ from those described by Tables 1 and 3.

Table 4.

Open in new tab

Comparison of simulations with different numbers of node masses and springs per node.

	Set
	A	B	C	D	E	F
〈N_I〉	403.8	405.8	795.2	804.2	1597.8	1596.2
〈NS_I/N_I〉	10.92	18.7	11.5	20.34	12.02	21.34
d_I	0.19	0.19	0.15	0.15	0.118	0.118
d_s/d_I	2.3	2.8	2.2	2.8	2.3	2.8
k_I	0.102	0.04	0.08	0.0301	0.062	0.023
γ_I	12.858	5.05	20.0	7.48	31.3	11.5
〈qf_ratio〉	1.48	1.39	1.41	1.36	1.40	1.38
σ[qf_ratio]	0.069	0.052	0.053	0.026	0.030	0.012
\|$\sigma [\dot{a}]/ \langle \dot{a} \rangle$\|	0.061	0.064	0.055	0.024	0.026	0.0094

	Set
	A	B	C	D	E	F
〈N_I〉	403.8	405.8	795.2	804.2	1597.8	1596.2
〈NS_I/N_I〉	10.92	18.7	11.5	20.34	12.02	21.34
d_I	0.19	0.19	0.15	0.15	0.118	0.118
d_s/d_I	2.3	2.8	2.2	2.8	2.3	2.8
k_I	0.102	0.04	0.08	0.0301	0.062	0.023
γ_I	12.858	5.05	20.0	7.48	31.3	11.5
〈qf_ratio〉	1.48	1.39	1.41	1.36	1.40	1.38
σ[qf_ratio]	0.069	0.052	0.053	0.026	0.030	0.012
\|$\sigma [\dot{a}]/ \langle \dot{a} \rangle$\|	0.061	0.064	0.055	0.024	0.026	0.0094

Notes. With a₀ = 10, M* = 100, n = 0.318, σ₀ = −0.4, P_o = 4.377, |$\bar{\chi }\sim 0.225$|⁠, E_I ∼ 2.3, and τ_relax ∼ 0.155.

Each column represents a group of five simulations. From each group of five the mean number of nodes and springs per node are listed in the second and third rows. The rows labelled 〈qf_ratio〉 and σ[qf_ratio] give the mean and standard deviations, respectively, of the ratio of the numerically measured to analytically predicted quality function. The bottom row gives the ratio of the standard deviation in the semimajor axis drift rate divided by the means. Each of the means and standard deviations are computed from five simulations. The simulations listed here differ from those described by Tables 1 and 3.

In the simulations listed in Tables 1 and 3, the ratio of the number of springs per node is about 11, which is slightly lower than the number recommended by Kot et al. (2015, fig. 5). For this ratio, Kot et al. (2015) found that the spring network behaved 10 per cent weaker⁶ than what was estimated by equation (29). In Table 4 we can compare the ratio of numerically measured to predicted quality functions for simulations with similar numbers of nodes but differ by the number of springs per node. Comparing columns A to B and C to D and E to F we see that the ratio of numerical measured to predicted quality function is only slightly less (about 2 per cent lower) when the number of springs per node is about 20 rather than 10. We find that for greater than 10 springs per nodes, the ratio of measured to predicted quality function is relatively insensitive to the number of springs per node.

As each simulation generates a new particle distribution and spring network, we can measure effects due to variations in these properties by comparing simulations generated from identical input parameters. Table 4 shows that the standard deviation of the quality function ratio and semimajor axis drift rate decreases with increased particle number and that the scatter in these measured quantities differs by only a few per cent. We conclude that variations in the spring network are unlikely to explain the discrepancy between predicted and numerically measured quality function.

One possible reason for the discrepancy between predicted and measured quality function is that near the body surface the number of springs per particle is lower than the interior, and thus the spring network is anisotropic near the surface. The simulated body is weaker (and floppier) than we estimated by integrating the spring properties over the volume. Because the overall rigidity of the body is weaker, its tidal response would be larger than we predicted analytically. This would be consistent with our greater than unity qf_ratio. If this were the primary cause of our amplitude discrepancy, then the ratio of numerical computed to analytical predicted quality function should inversely depend on the number of simulated masses. Table 4 shows that this ratio qf_ratio does decrease as the particle number is increased from 400 to 1600, however, it only decreases by about a per cent between 800 and 1600 particles. Convergence has not been achieved, but the standard deviations in quantities measured from groups of simulations imply that the simulation results are reproducible and that the scatter is smaller than the amplitude discrepancy. If doubling the particle number reduces the quality function ratio by 2 per cent, then to reach a quality function ratio of 1 we would require a million particles in the simulation and we are not yet set up to run simulations this large.

The surface of our body lies within a radius of 1, consequently our simulated body effectively has an outer radius that is smaller than 1. As the semimajor axis drift rate is faster for larger bodies, we would have expected our simulations to exhibit slower rather than faster orbital decay rates (and evident from the factor of R⁻⁵ in equation 34). We recall that the analytically predicted quality function depends on the ratio μ/e_g (see equations 24 and 25). As the energy density, e_g∝R⁴, for an effective radius less than 1, the energy density is higher than we previously estimated and the ratio μ/e_g would be lower than we used. This means we have underestimated the size of the tidal response. This is in the right direction to account for the amplitude discrepancy. However if we take into account both the factor of R⁵ in equation (34) when computing our numerical quality function and the factor of R⁴ in equation (25) when computing our analytical quality function, the two corrections more or less cancel each other and we do not resolve the source of discrepancy in quality function amplitude. So we find that we cannot resolve our discrepancy in amplitude by correcting the body radius to an effectively smaller radius.

Equations (29) and (30) characterize static properties of the mass-spring model. However, the analytical theory of evolving tides (including that based on the viscoelastic rheology 15) employs the unrelaxed shear rigidity μ. We have assumed that the two values are equivalent, although this is not perfectly true for real materials. As any rheological parameter, the shear elasticity modulus depends upon frequency. The unrelaxed or frequency-dependent shear modulus in real materials can be lower than the relaxed or static counterpart (e.g. Faul & Jackson 2005). Perhaps this is also true in our simulated material. The numerical tests by Kot et al. (2015) used static forces and so would not have been sensitive to frequency dependence in the shear modulus. Had we employed a smaller value for the shear modulus, this would have increased the values of the theoretically obtained k₂ sin ϵ₂, and thus would have reduced the offset between our numerical and theoretical values for the quality function.

Each numerical run starts out with a sphere of particles, with springs at their rest lengths. However, after the simulation begins, the body is compressed by self-gravity. The resulting density profile is not perfectly flat, as the centre becomes more compressed than the regions closer to the surface. Furthermore, the initial spin of the body causes its shape to deviate from a sphere. To sample a broad range of frequencies (in units of the relaxation time), and to have a stronger tidal response (i.e. to reduce the simulation time), we use a soft body that is particularly prone to deformation when spun and compressed by self-gravity. Nevertheless our simulated bodies are not strongly deformed and we do not expect body deformation and compression to account for our quality function amplitude discrepancy.

Since our simulated body is composed of randomly distributed point particles, the body is neither exactly spherical nor uniform. Initially, its principal axes of inertia (computed from the moment of inertia tensor) are oriented randomly, and the three moments of inertia are not exactly equal. So, initially, the body's spin angular momentum is not exactly perpendicular to the orbit. In the course of the simulations, we measured the values of the x and y components of the spin angular momentum, finding them to be a few hundredths of the initial z component. This is small enough that the spin about non-polar axes is unlikely to be the cause of our quality function amplitude discrepancy.

To summarize, we have compared simulations with different numbers of mass nodes and springs per node and found that the scatter and sensitivity of the ratio of numerically measured to analytically predicted quality function are only weakly dependent on these quantities. Convergence has not been reached at 1600 simulated particles but the scatter in the simulations is low enough that the variations in the spring network are small enough that we have reproducible numerical measurements and scatter within these measurement are smaller than our measured amplitude discrepancy. We have not identified the source of the 30 per cent discrepancy in the amplitude of our numerically measured quality function, however, we suspect variations in the spring network, floppiness in the outer surface and a possible frequency dependence in the behaviour of the numerically simulated shear modulus.

5 CONCLUSION

In this paper, we have used a self-gravitating damped mass-spring model, within an N-body simulation, to directly model the tidal orbital evolution and rotational spin-up of a viscoelastic body.

We considered a binary composed of an extended body (assumed spherical and homogeneous) and a point-mass companion. Within this setting, we have tested our numerical approach against an analytical calculation for this simple case. The semimajor axis’ tidal evolution rate was calculated by two methods. One was a direct simulation based on simulating the first body with a mass-spring network. Another, analytical, method was based on a pre-conceived viscoelastic model of tidal friction (the Kelvin–Voigt model). The numerically computed tidal evolution of the semimajor axis and the spin rate were a direct outcome of the damped mass-spring model. We compared the results obtained by the two methods, and concluded that a mass-spring network can serve as a faithful model of a tidally deformed viscoelastic celestial body. Specifically, we computed the quality function (the ratio of Love number k₂ and the quality factor Q) from the numerically simulated tidal drift of the semimajor axis. The quality function showed a strong frequency dependence that is close to the dependence derived analytically for a Kelvin–Voigt sphere by a method suggested by Efroimsky (2015).

While direct estimates for the global shear and bulk rigidity can be derived from the spring constants and spring lengths, the shear viscosity has not been estimated in the literature. Consequently, we were uncertain of our estimates for the shear viscosity and the associated relaxation time (equation 32). However, a comparison between our computed quality function (specifically the peak frequency) and that predicted analytically for a Kelvin–Voigt solid suggests that our estimates for the numerically simulated shear viscosity and viscoelastic relaxation time were correct.

The magnitude of our numerically predicted quality function is about 30 per cent larger than the one predicted analytically. By comparing simulations performed for two different initial values of the semimajor axis, we concluded that the cause is not the neglect of higher order terms in the quadrupole expansion for the potential. We find the ratio of numerically measured to predicted quality function is only weakly dependent on numbers of mass nodes and springs per node simulated but does slowly decrease with an increase in the numbers of particles simulated. We have not yet identified the cause of the discrepancy. We suspect that the non-uniformity of the spring network at the body surface could have caused a larger than expected tidal response. Alternatively the simulated shear modulus could be frequency dependent and overestimated by its computed static value.

Our study demonstrates that we can directly simulate the tidal evolution of viscoelastic bodies. We currently achieve an accuracy of 30 per cent but with ongoing effort we may improve upon this. It would be difficult to model this process over long (Myr or Byr) time-scales, because the time step is determined by the number of particles within the body and by the spring constants for the links connecting the particles. This requires the time step to be much shorter than the orbital time-scale. Nonetheless, mass-spring models can be used to explore the tidal evolution of inhomogeneous and anisotropic bodies, and to study how their quality functions depend on the rheological properties and internal structure of the body. This fully numerical approach also permits study of tidal phenomena that are not easy to predict analytically – such as capture into (and crossing of) spin–orbit and spin–spin resonances, tidally induced orbital evolution, the distribution of tidal heating, and the effects of non-linear rheological behaviour.

We thank Benoît Noyelles for helpful comments that helped improve the quality of the paper. We also wish to thank Harry Braviner and Moumita Das for helpful discussions, as well as Maciej Kot, Piotr Szymczak, Hanno Rein, and Darin Ragozzine. This work was in part supported by the NASA grant NNX13AI27G.

1

Darwin's work is presented, in the modern notation, in Ferraz-Mello, Rodríguez & Hussmann (2008).

2

See Section A5 of Appendix A.

3

A somewhat different approach to tides, explored by Ferraz-Mello (2013, 2015a), does not employ a constitutive equation explicitly. That model also predicts a similar shape for the quality function.

4

While it is still unknown what rheological models should describe comets and rubble-pile asteroids, both seismic and geodetic data indicate that the Earth's mantle behaves viscoplastically as an Andrade body; see Efroimsky & Lainey (2007), and Efroimsky (2012a,b). At very low frequencies, its behaviour becomes Maxwell (Karato & Spetzler 1990).

5

Keep in mind that the assumptions of homogeneity and compressibility in Love's theory are mutually contradictive, wherefore this theory can be used only as an approximation (Melchior 1972).

6

Because of a misprint, fig. 5 in Kot et al. (2015) actually displays (E₀ − E)/E₀ instead of (E − E₀)/E₀, with E₀ and E being the estimated and measured Young's modulus, respectively (Kot et al., private communication). A smaller average number 〈S〉 of springs per node gives a smaller Young's modulus and a weaker body.

7

This happens when the spin of a (tidally despun) host star is faster than the orbital period of a close planet orbiting it.

8

The reason why summation in the equation (A1) goes over l ≥ 2 is explained, e.g. in Efroimsky & Williams (2009, equations 5–11).

9

The caveat ‘in time’ is important. Lagging in time does not necessarily imply geometric lagging of the bulge. The lunar orbit being above synchronous, the main (semidiurnal) tide created by the Moon on the Earth always leads, not lags. This, however, gets along well with causality.

10

With a being the semimajor axis and G the gravity constant, the mean anomaly |${\cal {M}}(t)={\cal {M}}_0(t)+\int ^t_{}\,{\rm d}t\,\sqrt{G(M + M^{\ast })/a^3\,}$| renders the anomalistic mean motion as |$n\equiv {{\dot{\cal {M}}}}={{\dot{\cal {M}}}}_0+\sqrt{G(M + M^{\ast })/a^3\,}$|⁠. In neglect of external perturbations, |${{\dot{\cal {M}}}}_0\approx 0$| and the anomalistic mean motion can be approximated with the Keplerian mean motion: |$n \approx \sqrt{G(M + M^{\ast })/a^3\,}$|⁠.

11

Sometimes m is also referred to as the azimuthal wavenumber (Ogilvie 2014).

12

Generally, the two bodies should be treated on equal footing, so this potential should be amended with a similar potential wherewith the extended body is acted upon due to the tides it is exerting on the perturber.

13

The evolution of the argument of the pericentre ω, the longitude of the node Ω, and the mean motion |${\cal {M}}$| depends overwhelmingly on k_l(ω_lmpq) sin ϵ_l(ω_lmpq), and also contains terms with k_l(ω_lmpq) cos ϵ_l(ω_lmpq). The latter terms, though, are very small. It can also be shown that the secular drift of the semimajor axis a, eccentricity e, and inclination i is defined exclusively by the quality functions k_l(ω_lmpq) sin ϵ_l(ω_lmpq), with no terms containing k_l(ω_lmpq) cos ϵ_l(ω_lmpq).

REFERENCES

Darwin

G. H.

1879

Philos. Trans. R. Soc. Lond.

170

447

Crossref

Search ADS

Efroimsky

M.

2012a

Celest. Mech. Dyn. Astron.

112

283

Crossref

Search ADS

Efroimsky

M.

2012b

ApJ

746

150

Crossref

Search ADS

Efroimsky

M.

2015

AJ

150

98

Crossref

Search ADS

Efroimsky

M.

Lainey

V.

2007

J. Geophys. Res.

112

E12003

Crossref

Search ADS

Efroimsky

M.

Makarov

V. V.

2013

ApJ

764

26

Crossref

Search ADS

Efroimsky

M.

Williams

J. G.

2009

Celest. Mech. Dyn. Astron.

104

257

Crossref

Search ADS

Escribano

B.

Vanyó

J.

Tuval

I.

Cartwright

J. H. E.

González

D. L.

Piro

O.

Tél

T.

2008

Phys. Rev. E

78

036216

Crossref

Search ADS

Faul

U. H.

Jackson

I.

2005

Earth Planet. Sci. Lett.

234

119

Crossref

Search ADS

Ferraz-Mello

S.

2013

Celest. Mech. Dyn. Astron.

116

109

Crossref

Search ADS

Ferraz-Mello

S.

2015a

Celest. Mech. Dyn. Astron.

122

359

Crossref

Search ADS

Ferraz-Mello

S.

2015b

A&A

579

A97

Crossref

Search ADS

Ferraz-Mello

S.

Rodríguez

A.

Hussmann

H.

2008

Celest. Mech. Dyn. Astron.

101

171

Crossref

Search ADS

Henning

W. G.

Hurford

T.

2014

ApJ

789

30

Crossref

Search ADS

Karato

S.

Spetzler

H. A.

1990

Rev. Geophys.

28

399

Crossref

Search ADS

Kaula

M.

1964

Rev. Geophys.

2

661

Crossref

Search ADS

Kot

M.

Nagahashi

H.

Szymczak

P.

2015

The Visual Comput.: Int. J. Comput. Graphics

31

1339

Lainey

V.

et al.

2012

ApJ

752

14

Crossref

Search ADS

Love

A. E. H.

1911

Some Problems of Geodynamics

Cambridge Univ. Press

Cambridge

Also see a later edition by Dover, New York, 1967

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Makarov

V. V.

2012

ApJ

752

73

Crossref

Search ADS

Makarov

V. V.

2013

MNRAS

434

L21

Crossref

Search ADS

Makarov

V. V.

2015

ApJ

810

12

Crossref

Search ADS

Makarov

V. V.

Efroimsky

M.

2013

ApJ

764

27

Crossref

Search ADS

Makarov

V. V.

Berghea

C.

Efroimsky

M.

2012

ApJ

761

83

Crossref

Search ADS

Makarov

V. V.

Frouard

J.

Dorland

B.

2016

MNRAS

456

665

Crossref

Search ADS

Mase

G. T.

Smelser

R. E.

Mase

G. E.

2010

Continuum Mechanics for Engineers

3rd edn

CRC Press

Boca Raton, FL

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Melchior

P.

1972

Physique et Dynamique Planétaires

Vander Éditeur

Bruxelles

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Nealen

A.

Müller

M.

Keiser

R.

Boxerman

E.

Carlson

M.

Ageia

N.

2006

Comput. Graphics Forum

25

809

Crossref

Search ADS

Noyelles

B.

Frouard

J.

Makarov

V. V.

Efroimsky

M.

2014

Icarus

241

26

Crossref

Search ADS

Ogilvie

G. I.

2014

ARA&A

52

171

Crossref

Search ADS

Ostoja-Starzewski

M.

2002

Appl. Mech. Rev.

55

35

Crossref

Search ADS

Quillen

A. C.

Gianella

D.

Shaw

J.

Ebinger

C.

2015

‘Crustal Failure on Icy Satellites and Moons from a Strong Tidal Encounter

Available at http://arxiv.org/abs/1512.02154

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Rein

H.

Liu

S.-F.

2012

A&A

537

A128

Crossref

Search ADS

Remus

F.

Mathis

S.

Zahn

J.-P.

Lainey

V.

2015

A&A

573

A23

Crossref

Search ADS

Richardson

D. C.

Michel

P.

Walsh

K. J.

Flynn

K. W.

2009

Planet. Space Sci.

57

183

Crossref

Search ADS

Sánchez

P.

Scheeres

D. J.

2011

ApJ

727

120

Crossref

Search ADS

Williams

J. G.

Boggs

D. H.

2015

J. Geophys. Res.: Planets

120

689

Crossref

Search ADS

APPENDIX: SEVERAL BASIC FACTS ON BODILY TIDES

Tidal interactions play a key role in evolution of planetary systems and multiple stars. Their signature is observed, e.g. in the synchronized spin of the Moon, the pas de deux of Pluto and Charon, and in the 3:2 spin–orbit resonance of Mercury. Slowly but steadily, tides work to circularize the orbits of planets and moons – or, in some cases, to make orbits eccentric.⁷ Tidal dissipation warms up close-in moons (like the volcanic Io), close-in planets (like bloated Jupiters), and short-period binary stars (which can experience tidal coalescence).

Referring the reader to Efroimsky (2015), Efroimsky & Makarov (2013), and references therein for a more detailed introduction, here we provide a minimal kit of ideas and formulae needed to talk about bodily tides.

A1 Static tides

Consider an extended spherical body of mass M and radius R, tidally distorted by a perturber of mass M* located in an exterior point |${{\boldsymbol {r}}}$|⁠, so |$|{\boldsymbol {r}}|\ge R$|⁠. In a surface point |${\boldsymbol {R}}$| of the body, the potential due to the perturber is⁸

\begin{eqnarray} W({\boldsymbol {R}},\,{\boldsymbol {r}})\ =\ \sum _{{{\it l}}=2}^{\infty }W_{{\it l}}({\boldsymbol {R}},\ {\boldsymbol {r}}) , \end{eqnarray}

(A1)

where the inputs |$W_{{\it l}}({\boldsymbol {R}},\ {\boldsymbol {r}})$| are proportional to the appropriate Legendre polynomials P_l(cos γ), with γ being the angle between the vectors |${{\boldsymbol {r}}}$| and |${\boldsymbol {R}}$| pointing from the body centre. The integers l are termed the degrees.

The l-degree term |$W_{{\it l}}({\boldsymbol {R}},\,{\boldsymbol {r}})$| of the perturber's potential causes a tidal deformation of the perturbed body, assumed to be linear. Then the resulting lth addition U_l to the perturbed body's potential is also linear in W_l:

\begin{equation*} U({\boldsymbol {r}}\,^{\prime })\ =\ \sum _{l=2}^{\infty }\ U_{l}({\boldsymbol {r}}\,^{\prime })\ =\ \sum _{l=2}^{\infty }\ k_{l}\;\left(\,\frac{R}{r\,^{\prime }}\,\right)^{l+1}\;W_{l}({\boldsymbol {R}}\,,\;{\boldsymbol {r}}),\end{equation*}

|${\boldsymbol {r}}\,^{\prime }$| being an exterior point, and k_l being the static Love numbers.

Distorting the extended body, the perturber experiences its response in the form of the incremental potential U taken at the point |${\boldsymbol {r}}\,^{\prime }\,=\,{\boldsymbol {r}}$|⁠:

\begin{eqnarray} U({\boldsymbol {r}})\ =\ \sum _{l=2}^{\infty }\ U_{l}({\boldsymbol {r}})\ =\ \sum _{l=2}^{\infty }\ k_{l}\;\left(\,\frac{R}{r}\,\right)^{l+1}\;W_{l}({\boldsymbol {R}}\,,\;{\boldsymbol {r}}). \end{eqnarray}

(A2)

As the perturber is exterior (⁠|$|{\boldsymbol {r}}|>R$|⁠), the quadrupole part of the expansion for the perturbing potential W is dominant. The same pertains to U.

A2 Evolving tides

The case of evolving tides is more complicated. Owing to the internal friction, the tidal deformation (and the resulting additional potential U) always lags in time⁹ behind the perturbation W. To take into account different lagging at different frequencies, it is necessary to expand both the perturbing potential W and the response U in Fourier series. The linearity of response implies that the same frequencies should emerge in both spectra, when W and U are observed at the same point of space. From the cornerstone work by Kaula (1964), it is easy to derive that the Fourier tidal modes read as

\begin{eqnarray} \omega _{\textstyle {_{lmpq}}}&=&(l-2p)\;\dot{\omega }\,+\,(l-2p+q)\;n\,+\,m\;(\dot{\Omega }\,-\,\dot{\theta })\nonumber \\ &\approx &(l-2p+q)\;n\,-\,m\;\dot{\theta }, \end{eqnarray}

(A3)

where θ and |${\dot{\boldsymbol \theta }}$| are the rotation angle and rotation rate of the extended body, introduced in the equatorial plane. In neglect of the equinoctial precession, θ can be identified with the sidereal angle. The notations ω and Ω stand for the perturber's argument of the pericentre and the longitude of the node, as seen from the extended body. The formula also includes the mean anomaly |${\cal {M}}$| and the anomalistic¹⁰ mean motion |$n\equiv {{\dot{\cal {M}}}}$| (with |${\cal {M}}=0$| at the pericentre). Derivation of the expression (A3) is explained in section 4.3 of Efroimsky & Makarov (2013).

The modes |$\omega _{\textstyle {_{lmpq}}}$| can be of either sign, while their absolute values,

\begin{eqnarray} \chi _{\textstyle {_{lmpq}}}\,=\,|\,\omega _{\textstyle {_{lmpq}}}\,|\ \approx \ |\,(l-2p+q)\;n\,-\,m\;\dot{\theta }\,|, \end{eqnarray}

(A4)

have the meaning of positive definite forcing frequencies of stresses and strains in the distorted body. The Fourier modes are parametrized with the four integers l, m, p, q. The integers l and m are the degree and order of the spherical harmonics employed in the expansion.¹¹

The dynamical analogue to the formula (A1) is

\begin{eqnarray} W({\boldsymbol {R}},\,{\boldsymbol {r}},\,t)=\sum _{l=2}^{\infty } W_{l}({\boldsymbol {R}},\,{\boldsymbol {r}},\,t) =\sum _{lmpq}\ W_{{lmpq}}({\boldsymbol {R}},\,{\boldsymbol {r}},\,t), \end{eqnarray}

(A5)

where a term W_lmpq is proportional to cos (ω_lmpq t + ⋅⋅⋅), with ellipsis denoting some phase:

\begin{eqnarray} W_{{lmpq}}({\boldsymbol {R}},\,{\boldsymbol {r}},\,t)=A_{{lmpq}}({\boldsymbol {R}},\,{\boldsymbol {r}},\,t)\;\cos (\omega _{lmpq}\,t\,+\cdots ). \end{eqnarray}

(A6)

Both the static formula (A1) and its dynamical analogue (A5) render the value of the perturbing potential at a surface point |${\boldsymbol {R}}$|⁠.

Writing down a dynamical analogue to the static expression (A2) turns out to be a highly nontrivial problem. Above we stated that, owing to the linearity of the problem, the spectrum of U should contain the same frequencies as that of W, provided both U and W are observed at the same point of space. Therefore, a Fourier series for U would contain terms proportional to cos (ω_lmpq t + ⋅⋅⋅), had it been written for the (evolving in time) value of U at the same surface point |${\boldsymbol {R}}$|⁠. We however are interested in the values of U in a different point, the point |${\boldsymbol {r}}$| where the moving perturber is located. There, the spectrum of |$U({\boldsymbol {r}},\,t)$| will be richer than that of |$W({\boldsymbol {R}},\,{\boldsymbol {r}},\,t)$|⁠, and will be parametrized with six indices lmpqhj:

\begin{eqnarray} U({\boldsymbol {r}},\,t)\ =\ \sum _{{{\it l}}=2}^{\infty }U_{{l}}({\boldsymbol {r}})\ =\ \sum _{lmpqhj}\ U_{{lmpqhj}}({\boldsymbol {r}},\,t) ; \end{eqnarray}

(A7)

see Efroimsky (2012a, sections 7 and 8). As was pointed out by Kaula (1964), |$U({\boldsymbol {r}},\,t)$| contains a secular part – and that part is parametrized with the four indices lmpq:

\begin{eqnarray} \langle \,U({\boldsymbol {r}},\,t)\,\rangle \ =\ \sum _{{{\it l}}=2}^{\infty }\langle \,U_{l}({\boldsymbol {r}},\,t)\,\rangle \ =\ \sum _{lmpq}\ \langle U_{{lmpq}}({\boldsymbol {r}}) \rangle , \end{eqnarray}

(A8)

where the angular brackets 〈⋅⋅⋅〉 denote time averaging, and the terms on the right-hand side are given by

\begin{eqnarray} \langle U_{{lmpq}}({\boldsymbol {r}})\rangle \ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \nonumber \\ =k_l(\omega _{lmpq})\ \cos \epsilon _l(\omega _{lmpq})\ \left(\,\frac{R}{r}\,\right)^{l+1}\;A_{lmpq}({\boldsymbol {R}}\,,\;{\boldsymbol {r}}), \end{eqnarray}

(A9)

where A_lmpq is the magnitudes from the formula (A6), while k_l(ω_lmpq) and ϵ_l(ω_lmpq) are the degree l dynamical Love numbers and phase lags written as functions of the Fourier modes.

A3 The secular part of the tidal torque acting on the spin of the extended body

The negative gradient of the secular potential (A8) renders the secular part of the orbital torque wherewith the extended body is acting on the perturber. An equal but opposite torque is acting on the extended body and is influencing its spin. The polar component of the secular torque reads as

\begin{eqnarray} \langle \,{\cal {T}}^{(z)}\,\rangle \ =\ \sum _{{{\it l}}=2}^{\infty }\left\langle \,{\cal {T}}_l^{(z)}\,\right\rangle \ =\ \sum _{lmpq}\ \left\langle \,{\cal {T}}_{lmpq}^{(z)}\,\right\rangle , \end{eqnarray}

(A10)

where

\begin{eqnarray} &&{\left\langle \,{\cal {T}}_{lmpq}^{(z)}\,\right\rangle} \nonumber \\ &&{\quad=k_l(\omega _{lmpq})\ \sin \epsilon _l(\omega _{lmpq})\ \left(\,\frac{R}{r}\,\right)^{l+1}m\;A_{lmpq}({\boldsymbol {R}}\,,\;{\boldsymbol {r}}).} \end{eqnarray}

(A11)

We see that an lmpq component of the torque may be either decelerating or accelerating the spin, dependent upon the sign of the phase lag ϵ_l(ω_lmpq) – which always coincides with the sign of the Fourier mode ω_lmpq.

A4 The quality function (‘kvalitet’)

The product k_l(ω_lmpq) sin ϵ_l(ω_lmpq) is sometimes termed as the quality function (Makarov 2013; Efroimsky 2015) or kvalitet (Makarov 2015; Makarov et al. 2016). In the literature, it is conventional to write it as

\begin{eqnarray} k_l(\omega _{lmpq})\;\sin \epsilon _l(\omega _{lmpq})\ =\ \frac{k_l(\omega _{lmpq})}{Q_l(\omega _{lmpq})}\ \,\mbox{Sgn}\,\omega _{lmpq}, \end{eqnarray}

(A12)

where the quality factors are introduced via

\begin{eqnarray} \frac{1}{Q_l(\omega _{lmpq})}\,=\,|\,\sin \epsilon _l(\omega _{lmpq})\,|, \end{eqnarray}

(A13)

and where it is taken into account that the sign of a phase lag ϵ_l(ω_lmpq) always coincides with the sign of the Fourier mode ω_lmpq (e.g. Efroimsky & Makarov 2013).

A5 Which terms are leading, and when

As the perturber is exterior (⁠|$|{\boldsymbol {r}}|>R$|⁠), the quadrupole part of the expansion for the perturbing potential W is dominant. The quadrupole part comprises all the terms with l = 2. For low inclination and eccentricity, the largest terms in the expansions (A5), (A7), and (A10) are those with {lmpq} = {2200}. They correspond to the so-called semidiurnal Fourier mode:

\begin{eqnarray} \omega \, \, \equiv \ \omega _{2200}\,=\ 2\ (n\ -\ \dot{\theta }) . \end{eqnarray}

(A14)

When the semidiurnal, or any other lmpq term is leading in the expansion for W, the corresponding lmpq term is leading also in the expansion for the additional tidal potential U. Up to some reservation, this is true also for the expansions of the tidal torque. A reservation comes from the fact that an lmpq term in the expansion for the torque contains as a multiplier the sine of the phase lag ϵ_l(ω_lmpq). For example, in the case of small inclination i and eccentricity e, the semidiurnal part of the polar torque operating on the spin of the perturbed body reads as (Efroimsky 2012a)

\begin{eqnarray} {\cal {T}}^{(z)}_{2200} &=& \frac{3}{2} G {M^{\ast }}^2 \frac{R^5}{a^6} k_2(\omega _{2200})\;\sin \epsilon _2(\omega _{2200})\nonumber\\ &&+\,O(e^2\,\epsilon )\,+\,O(i^2\,\epsilon ). \end{eqnarray}

(A15)

The quality function k_l(ω_lmpq) sin ϵ_l(ω_lmpq) continuously goes through zero (and changes its sign) when the lmpq spin–orbit resonance is transcended, i.e. when ω_lmpq goes through zero. So, when a rotator is trapped into an lmpq spin–orbit resonance, the quality function stays zero; so the Fourier mode ω_lmpq contributes nothing to the torque. Specifically, in the case of synchronous rotation (known as the 1:1 spin–orbit resonance), the mode ω₂₂₀₀ vanishes – and so does the semidiurnal term of the torque. In the resonance, therefore, it is the higher than semidiurnal terms that are leading.

This ‘acceding of leadership’ in resonances, along with its physical consequences for binaries, is described in detail in Makarov & Efroimsky (2013) and Makarov, Berghea & Efroimsky (2012). Here we shall only mention two simple examples. Since the Moon is synchronized, the semidiurnal input into the torque acting on its spin is zero. It is then the other components (mainly, the term with {lmpq} = {2201}) that define the tidal response of the Moon and influence its libration in longitude (Makarov et al. 2016). As another example, take Mercury in its 3:2 resonance (Noyelles et al. 2014). For this planet, the {lmpq} = {2201} input into its tidal response is zero, and it is the semidiurnal mode that overwhelmingly defines the tidal response and plays a crucial role in longitudinal libration (Makarov, in preparation).

A6 Tidally generated secular orbital evolution

The tidal potential (equation A8) should be inserted into the Lagrange- or Delaunay-type planetary equations to calculate the secular evolution of the orbit.¹² It then turns out after some algebra that the secular orbital evolution is determined mainly by the quality function k_l(ω_lmpq) sin ϵ_l(ω_lmpq), with the cosine of the lag playing a very marginal role.¹³

However, in some situations, approximate secular evolution can be calculated via the tidal torque. For example, consider an orbit with no inclination relative to the equator of the extended body. Using the expression

\begin{eqnarray} \;L^{({\rm orb})}\;=\;\frac{M M^{\ast }}{M + M^{\ast }} \sqrt{G (M + M^{\ast })} \sqrt{a(1-e^2)} \end{eqnarray}

(A16)

for the orbital angular momentum, and setting there e = 0, we can use the conservation of the angular momentum |$ \dot{L}^{({\rm orb})}\,=\,-{\cal {T}}^{(z)} $| to derive the evolution rate of the semimajor axis:

\begin{eqnarray} \frac{\dot{a}}{n\,a} &=&\ -\frac{2\,{\cal {T}}^{(z)} \,a}{G M^{\ast } M} \nonumber \\ &=& -\;3 \left( \frac{M^{\ast }}{M} \right)\left( \frac{R}{a} \right)^5\;k_l(\omega _{lmpq})\;\sin \epsilon _l(\omega _{lmpq}) . \end{eqnarray}

(A17)

The expression (A17) is a semidiurnal approximation. As was mentioned in Section A5, this approximation is valid everywhere except in the 1:1 resonance where the semidiurnal term vanishes.

Download all slides

Month:	Total Views:
January 2017	2
February 2017	3
March 2017	6
April 2017	1
May 2017	3
June 2017	6
July 2017	4
August 2017	3
September 2017	6
October 2017	1
November 2017	8
December 2017	2
January 2018	1
March 2018	4
April 2018	2
May 2018	4
June 2018	4
August 2018	2
October 2018	7
November 2018	4
December 2018	7
January 2019	2
February 2019	3
March 2019	1
April 2019	6
May 2019	6
June 2019	3
July 2019	12
August 2019	16
September 2019	15
October 2019	5
November 2019	17
December 2019	5
January 2020	5
February 2020	8
March 2020	9
April 2020	5
May 2020	7
June 2020	20
July 2020	2
August 2020	7
September 2020	10
October 2020	8
November 2020	3
December 2020	7
January 2021	25
February 2021	14
March 2021	6
April 2021	10
May 2021	26
June 2021	18
July 2021	8
August 2021	14
September 2021	10
October 2021	7
November 2021	27
December 2021	18
January 2022	19
February 2022	11
March 2022	8
April 2022	11
May 2022	8
June 2022	59
July 2022	10
August 2022	8
September 2022	18
October 2022	19
November 2022	8
December 2022	11
January 2023	13
February 2023	11
March 2023	7
April 2023	52
May 2023	9
June 2023	16
July 2023	6
August 2023	54
September 2023	8
October 2023	5
November 2023	12
December 2023	11
January 2024	22
February 2024	16
March 2024	25
April 2024	61
May 2024	16
June 2024	8
July 2024	22
August 2024	13
September 2024	10
October 2024	16
November 2024	9
December 2024	11
January 2025	2
February 2025	21
March 2025	18
April 2025	12
May 2025	2

Article Contents

Numerical simulation of tidal evolution of a viscoelastic body modelled with a mass-spring network

Abstract

1 MOTIVATION AND PLAN

2 TIDES IN A KELVIN–VOIGT VISCOELASTIC SOLID

2.1 How tides work

2.2 The shape of the quality function

2.3 The quality function for a Kelvin–Voigt sphere

3 DAMPED MASS-SPRING MODEL SIMULATIONS

4 TIDAL EVOLUTION OF THE SEMIMAJOR AXIS

4.1 The semimajor axis’ tidal drift rate computed from the mass-spring model

4.2 Comparison of numerically measured and predicted quality functions

4.3 Discrepancy in amplitude of quality function and sensitivity to the number of masses and springs simulated

5 CONCLUSION

REFERENCES

APPENDIX: SEVERAL BASIC FACTS ON BODILY TIDES

A1 Static tides

A2 Evolving tides

A3 The secular part of the tidal torque acting on the spin of the extended body

A4 The quality function (‘kvalitet’)

A5 Which terms are leading, and when

A6 Tidally generated secular orbital evolution

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Numerical simulation of tidal evolution of a viscoelastic body modelled with a mass-spring network

Abstract

1 MOTIVATION AND PLAN

2 TIDES IN A KELVIN–VOIGT VISCOELASTIC SOLID

2.1 How tides work

2.2 The shape of the quality function

2.3 The quality function for a Kelvin–Voigt sphere

3 DAMPED MASS-SPRING MODEL SIMULATIONS

4 TIDAL EVOLUTION OF THE SEMIMAJOR AXIS

4.1 The semimajor axis’ tidal drift rate computed from the mass-spring model

4.2 Comparison of numerically measured and predicted quality functions

4.3 Discrepancy in amplitude of quality function and sensitivity to the number of masses and springs simulated

5 CONCLUSION

REFERENCES

APPENDIX: SEVERAL BASIC FACTS ON BODILY TIDES

A1 Static tides

A2 Evolving tides

A3 The secular part of the tidal torque acting on the spin of the extended body

A4 The quality function (‘kvalitet’)

A5 Which terms are leading, and when

A6 Tidally generated secular orbital evolution

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only