Numerical Simulation of Tidal Evolution of a Viscoelastic Body Modelled with a Mass-Spring Network

0.45 MB·

by  Julien Frouard

MNRAS000,000–000(0000) Preprint4March2016 CompiledusingMNRASLATEXstylefilev3.0
Numerical Simulation of Tidal Evolution of a Viscoelastic
Body Modelled with a Mass-Spring Network
Julien Frouard1⋆, Alice C. Quillen2, Michael Efroimsky1, David Giannella2
1USNaval Observatory,3450 Massachusetts Ave NW, Washington DC 20392 USA
2Department of Physics and Astronomy, Universityof Rochester, Rochester, NY 14627 USA
6
1
0
2 4March2016

r
a
M ABSTRACT
We use a damped mass-spring model within an N-body code to simulate the tidal
evolutionofthe spinandorbitofaself-gravitatingviscoelasticsphericalbodymoving
3
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
P number k2 dividedbythetidalqualityfactor Q)fromthenumericallycomputedtidal
E drift of the semimajor axis of the binary. The shape of k2/Q, as a function of the
. principal tidal frequency, reproduces the kink shape predicted by Efroimsky (2012a;
h
CeMDA 112:283) for the tidal response of near-spherical homogeneous viscoelastic
p
- rotators.Wedemonstratethatwecandirectlysimulatethetidalevolutionofspinning
o viscoelastic objects. In future, the mass-spring N-body model can be generalised to
r inhomogeneous and/or non-spherical bodies.
t
s
a Key words: planets and satellites: dynamical evolution and stability, interiors –
[
methods: numerical

2
v
2
2
1 MOTIVATION AND PLAN Ostoja-Starzewski (2002) and Kot et al. (2014) haveshown
2
that mass-spring systems can accurately model elastic ma-
8 The analytical theory of tidal interaction of solid bodies
terials. As we shall show, viscoelastic response and grav-
0 has a long and rich history — from the early mathemat-
itational forces can be incorporated into the mass-spring
1. ical development by Darwin (1879),1 to its generalisation particle-based simulation technique. We demonstrate nu-
0 by Kaula (1964), to an avalanche of more recent results. mericalsimulationsofthetidalresponseofaspinninghomo-
6 Verification of tidal theories through direct measurements
geneoussphericalbodyexhibitingspin-downorupandasso-
1 is not easy because the tidal evolution is slow and requires
ciated drift in semi-major axis, without using an analytical
: eitherhighobservationalprecision(Williams & Boggs2015)
v tidalevolutionmodel.Thenwecomparetheoutcomeofour
oranextendedobservationaltimespan(Lainey et al.2012).
i numerics with analytical calculations. The results are simi-
X Purely analytical theories of tidal evolution describe homo-
lar,justifyingournumericalapproach.Ournumericalmodel
r geneous or, at best, two-layered (as in Remuset al. 2015)
may in future be used to take into account more complex
a near-sphericalbodiesof linear rheology.This makesnumer-
effects that cannot be easily computed by analytical means
ical simulations attractive (e.g., Henning& Hurtford 2014)
(like inhomogeneity, compressibility, or a complex shape of
as they may make it possible to explore the tidal evolution
thetidally perturbed body).
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. Be- 2 TIDES IN A KELVIN-VOIGT
causeoftheirsimplicityandspeed(comparedtomorecom- VISCOELASTIC SOLID
putationally intensive grid-based or finite element meth-
2.1 How tides work
ods), mass-spring computations are a popular method for
simulatingsoft deformablebodies(e.g., Nealen et al.2006). Consideranextendedsphericalbodyofmass M andradius
∗
R,tidallydeformedbyaperturberofmass M residingata
position r,where |r|>R.Thebinary’sorbithasthesemi-
⋆ Contact e-mail:[email protected] majoraxis a andthemeanmotion n= G(M +M∗)/a3,
1 Darwin’s work is presented, in the modern notation, in where G is the gravitational constant.pFor a distant per-
Ferraz-Melloetal.(2008). turber (a ≫ R), the quadrupole part in the Fourier ex-
(cid:13)c 0000TheAuthors
2 Frouard et al.
pansion (A5) for the perturbing potential W is dominant. where
Thispartcomprisesseveralterms.Ofthese,thetermcalled 3(2l2+4l+3) 1 4π(2l2+4l+3)
semidiurnal is usually leading.2 This term is a function of Bl = 4lπGρ2R2 = eg 3l , (5)
theprincipaltidalFouriermode ω whichwedenotesim-
2200
ply as ω: while
ω ≡ ω2200 = 2(n − θ˙) , (1) eg ≡ GRM42 (6)
where θ and θ˙ are thebody’s rotation angle and spin rate
is (to order of magnitude) the gravitational energy density
in the equatorial plane. (See the equation (A14) in the Ap-
of the body.Here Jr =1/µr is thestatic (relaxed) compli-
pendix.)
anceofthematerial, whichisinversetothestatic(relaxed)
The secular part of the semidiurnal term of the polar
rigidity µr. Switching from a static to an evolving config-
tidal torqueacting on thebody is
uration, we invokethe equivalence principle for viscoelastic
hT(z) i = 3 GM∗2 R5 k (ω) sinǫ (ω) , (2) materials, in order to obtain the complex Love number in
2200 2 a6 2 2 thefrequencydomain:
where k2(ω) and ǫ2(ω) arethequadrupoledynamicalLove k¯(χ)= 3 1 =|k¯(χ)|e−iǫl(χ) , (7)
number and the quadrupole phase lag, both taken at the l 2(l−1)1+Bl/J¯ l
semidiurnal frequency given by the above expression (1).
χ≡|ω| beingthetidalfrequency,and J¯(χ) beingthecom-
The product k (ω) sinǫ (ω) is often called the quality
2 2 plexcomplianceofthematerial.Oncethecompliance J¯(χ)
function (Makarov 2012, 2013; Efroimsky 2015) or, some-
is prescribed by a rheological model, we can compute the
times, kvalitet (Makarov 2015; Makarov et al. 2016).
quantityfunction
Whentheinclinationandeccentricityaresmall,conser-
vationofangularmomentumgivesanestimateofthesecular k(χ) sinǫ(χ) = −Im[k¯(χ)] , (8)
l l l
drift rate of the semi-major axis (see the equation (A17) in
where k(χ) ≡ |k¯(χ)| and ǫ(χ) are,correspondingly,the
theAppendix): l l l
degree-l dynamical Love number and phase lag (the latter
a˙ 2T(z) a beinglinkedtothedegree-l tidalqualityfactorthroughthe
= − 2200
na GM∗M equation A13).
M∗ R 5 However, an lmpq term in the expansion (A10 - A11)
= −3(cid:18) M (cid:19)(cid:18)a(cid:19) k2(ω) sinǫ2(ω) . (3) forthetidaltorqueincludesthefactor kl sinǫl written not
asafunctionofthetidalfrequency χ≡|ω| butasafunction
Below we compute the drift rate of the semi-major axis a˙ of the tidal mode ω — see, e.g., the semidiurnal term
throughadirectnumericalsimulation,andthencomparethe given by the expression (2). It can be demonstrated that
result with that obtained analytically from thetidal theory thisbrings an extra factor equal tothe sign of themode:
usingahomogenousKelvin-Voigtviscoelasticrheology.This
comparison will demonstrate that we can directly simulate [klsinǫl](ω) =
thetidalevolutionofviscoelasticobjectswithamass-spring (9)
model. 3 B Im(J¯(χ))
− l ×Signω .
2(l−1) Re(J¯(χ))+B 2+ Im(J¯(χ)) 2
l
2.2 The shape of the quality function Whatever re(cid:2)alistic rheologica(cid:3)l com(cid:2)pliance J¯((cid:3)χ) is inserted
As is demonstrated in the Appendix, the Fourier decom- into the above formula, the shape of the quality function is
positions of the disturbing potential W, the tidal re- similar for all viscoelastic bodies. It exhibits a sharp kink
sponsepotential U,andthetidaltorque T compriseterms with two peaks having opposite signs (e.g., Noyelles et al.
that are numbered with the four indices lmpq. An lmpq 2014, Efroimsky 2015).3
term of the torque is proportional to the quality function Thegenericshapeofthequalityfunctioncanbeunder-
k(ω ) sinǫ(ω ). The quality function depends on stood by comparing the frequency χ to the inverse of the
l lmpq l lmpq
the degree l, the composition of the body, the rheology of viscoelastic relaxation time (e.g., Ferraz-Mello 2015b.). At
its layers, and the size and mass of the body. The size and afixedpoint inthebody,thetidal stressingin thematerial
mass are important, because the tidal response is defined is oscillating at the frequency χ ≡ |ω|. When χ is small
notonlybytheinternalstructureandrheology,butalsoby compared to the inverse timescale of viscoelastic relaxation
self-gravitation. in the material, the body deformation stays almost exactly
Theprocessofderivingthequalityfunctionforanylin- in phase with the tidal perturbation. The reaction and ac-
ear viscoelastic rheology is described in Efroimsky (2015). tionbeingvirtuallyinphase,noworkiscarriedoutandthe
Here we provide a short account of that derivation. We be- tidal effects are minimal. On the other hand, at very high
gin with the expression for the static Love number for an frequencies,thebody’sviscosity preventsit from deforming
homogeneous,incompressible,self-gravitatingelasticsphere, during the short forcing period 2π/χ. The reaction cannot
3 1
k(static) = , (4)
l 2(l−1) 1+Bl/Jr 3 A somewhat different approach to tides, explored by
Ferraz-Mello(2013,2015a),doesnotemployaconstitutiveequa-
tion explicitly. That model also predicts a similar shape for the
2 SeeSectionA5oftheAppendix. qualityfunction.
MNRAS000,000–000(0000)
3
catchupwiththeaction,andstaysclosetozero—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
Thegoalofourpaperisnottofavouraparticularrheological
model,buttotestoursimulationmethod.TheKelvin-Voigt
rheologicalmodelissimplistic,4butitistheeasiesttomodel
with a mass-spring model. Subsequent work will be aimed
at extendingour simulation approach to more realistic rhe-
ologies.
The Kelvin-Voigt model of a viscoelastic solid can be
Figure 1. Green circles are mass elements. In the Kelvin-Voigt
representedbyapurelyviscousdamperandapurelyelastic
model (the top drawing), spring and damping elements are set
springconnectedtotwomasselementsinparallel.Ifwecon-
inparallel.IntheMaxwellmodel(thebottomdrawing),thetwo
nect thesetwoelements in series ratherthanin parallel, we
elements are arranged in series. Since it is easier to represent
describetheMaxwellmodel;Figure1.AKelvin-Voigtbody aKelvin-Voigtviscoelasticmaterialwithamass-springnetwork,
is easier to model with a mass-spring system, because both wecomputethequalityfunctionfortheKelvin-Voigtmodel,thus
the spring and damping forces are directly applied to each allowing a direct comparison between the quality functions cal-
node particle. This can be seen from the expression for the culatedfromtheoryandthatcomputedthroughsimulation.
stresstensor σpq asafunctionofthestrainratetensor ε˙pq
in thetime domain (see Efroimsky 2012a):
theKelvin-Voigtbehaviorofthedampedsprings,described
t above,andthebulkandshearpropertiesofthematerial(i.e.
′ ′ ′
σpq(t)=2Z−∞µ(t−t)ε˙pq(t)dt , (10) theconstitutiveequations)isnotobvious,ourcomputations
demonstratethatthesimulatedresolvedbodybehavessimi-
′
where µ(t−t) isthestress-relaxationfunction.Inthecon- lartothatexpectedforaKelvin-Voigtsolidwiththeactual
textofamass-springmodel, wheretheforce betweenparti-
elastic modulus µ ≈ µI, viscosity η ≈ ηI, and relaxation
cles is likened to a uniaxial stress, the normal force applied
time
to theparticle i dueto a spring and dashpot connectingit
η
to particle j is given by τ = . (14)
µ
t
′ ′ ′ There may be a difference between our a priori estimated
Fi(t)=2 µ(t−t)ε˙(t)dt . (11)
Z−∞ values µI, ηI oftheshearrigidityandviscosity(computed
from the network), and the values µ, η of the simulated
For theKelvin-Voigtmodel (e.g., Mase et al. 2010),
material (see thediscussion byKot et al. 2014).
′ ′
µ(t−t)=µ+ηδ(t−t) (12) We now derive the quality function (9) for the case of
a Kelvin-Voigt sphere. We insert into the equation (9) the
where µ and η aretheunrelaxedshearrigidityandviscosity
complex compliance of a Kelvin-Voigtbody:
of the link between i and j. Inserting that rheology into
theequation (11),we find J¯(χ, µ, η) ≡ µ−iχη = µ−1 1−iχτ . (15)
µ2+χ2η2 1+χ2τ2
Fi(t)=2µε(t)−2µε(−∞)+2ηε˙(t) . (13)
The complex compliance J¯ can be presented either as a
Inneglectofthestrainat t=−∞,thisbecomesequivalent functionoftheunrelaxedshearmodulus µ andviscosity η,
toequations(26-28) presentedbelow andemployedin our oras a function of theshear modulus µ and therelaxation
code. time τ. Introducingthedimensionless Fourier tidal mode
Inmass-springmodelsimulations,massiveparticlesare
interlinked with a network of massless springs. To each of ω˜ ≡ ωτ , (16)
the two particles linked by a spring, a damping force is ap-
and following χ=|ω|, thedimensionless physicalfrequency
pliedthatdependsonthespringstrainrate(seeSection2.2
in Quillen et al. 2015). The shear elastic modulus, µI, can χ˜ ≡ χτ , (17)
be computed for an isotropic, initially random mass-spring
we write down the complex compliance as
model from the strength, lengths and distribution of the
springs(Kot et al.2014).Weproposeinthenextsectionan J¯(χ˜, µ) ≡ µ−1 1−iχ˜ . (18)
estimateforthesimulatedmaterialshearviscosity ηI based 1+χ˜2
onthespringdampingforces.Althoughtherelationbetween
Wedefineadimensionless function
¯j(l)(χ˜,µ) ≡ B−1J¯(χ˜,µ) = (B µ)−11−iχ˜ , (19)
4 While it is still unknown what rheological models should de- l l 1+χ˜2
scribecomets andrubble-pileasteroids, bothseismicand geode-
with B givenbytheexpression(5).Usingthatexpression,
l
ticdataindicatethattheEarth’smantlebehavesviscoplastically
we write down the function of thedegree l=2 as
as an Andrade body, see Efroimsky&Lainey (2007); Efroimsky
(2012a)andEfroimsky(2012b).Atverylowfrequencies,itsbe- ¯j(l=2)(χ˜, µ) = 3 eg 1 − iχ˜ . (20)
haviourbecomesMaxwell(Karato&Spetzler 1990). 38 π µ 1 + χ˜2
MNRAS000,000–000(0000)
4 Frouard et al.
This is a function of the shear modulus µ given in the tidalviscoelastic response weusethemass-spring modelby
units of eg, and of the frequency χ given in the units of Quillen et al. (2015), that is based on the modular N-body
theinverserelaxationtime τ.Intermsofthedimensionless coderebound (Rein & Liu2012).Arandomspringnetwork
Fourier mode ω˜ and the dimensionless frequency χ˜ ≡|ω˜|, model, rather than a lattice network model, was chosen so
theequation (9) becomes: that the modeled body is approximately isotropic and ho-
mogenousand lacks planesassociated with crystalline sym-
[k sinǫ ](ω˜, µ) =
l l metry. We work in units of radius and mass of the tidally
(21) perturbed body: R = 1, M = 1. Time is given in units
3 Im(¯j(l)(χ˜)) of tgrav = R3/GM referred to as the gravitational time
− ×Signω˜ .
2(l−1)[Re(¯j(l)(χ˜))+1]2+[Im(¯j(l)(χ˜))]2 scale.Presspure,energydensityandelasticmoduliarespeci-
fied in units of eg ≡ GM2/R4. In these units,thevelocity
This function attains its extrema at ofamassless particleon acircularorbit,whichisjust graz-
µB + 1 1 1 ing the surface of the body, is 1, the period of the grazing
ω(l) = ± l = ± 1 + .
peak Bl η τ (cid:18) µBl(cid:19) orbit being 2π.
Modelingthetidallyperturbedbody,werandomlygen-
For l=2, thedimensionless peak mode and frequencyare
erate an initial spherical distribution of particles of equal
ω˜(l=2) ≡ ω(l=2)τ = ± 1 + 3 eg (22a) mass (as described in Section 2.2 in Quillen et al. 2015);
peak peak (cid:18) 38 π µ(cid:19) see our Figure 2 for an illustration. Particle positions are
randomly generated using a uniform distribution in three
and
dimensionsbutaparticleisadded(asanode)tothespring
χ˜(l=2) ≡ χ(l=2)τ = |ω(l=2)|τ = 1 + 3 eg . (22b) network only if it is sufficiently separated from other parti-
peak peak peak (cid:18) 38 π µ(cid:19) cles (at a distance greater than minimum distance dI) and
Elasticbodiesshouldobey µ>eg,lesttheycollapsedueto withinaradiusof R=1 from thebodycentre.Springsare
self-gravity. Hencewe expect that addedbetween twonodesifthedistancebetweenthenodes
ωp(le=a2k) τ ≈ ±1 . (23) is lesTshtehapnardtiisctlaesncaeredssu.bjectedtothreetypesofforces: the
gravitational forces acting on every pair of particles in the
Usingthe shorthandnotation ω˜ = ω˜(l=2) and
body and with the massive companion, and the elastic and
y(µ,ω˜) ≡ 3 eg (1 + ω˜2)−1 , (24) dampingspringforces actingonly betweensufficientlyclose
38 π µ particle pairs. Springs with a spring constant kI intercon-
theequation (21) for l=2 can bewritten as necteachpairofparticlescloserthansomedistance ds.The
rest length of each spring is initially set equal to its initial
3 yω˜
[k sinǫ ](ω˜, µ) = . (25) length.Thespringsexperiencecompressionorextensionrel-
2 2 2 y2(1+ω˜2)+2y+1
ative to their rest length. The number of springs and their
We will use this expression for the quality function to ana- rest lengths stay fixed during our simulations. The springs
lytical predict tidal response. havedifferent rest lengths,and we denotetheiraverage ini-
Classical tidal theory is valid only for incompressible tial (rest) length with LI. The total number of particles is
materials(thosehavingthePoissonratio ν =0.5)—which NI, and the total number of springs is NSI. The mass of
is why the standard expression for the static Love number each particle is mI = 1.0/NI, and the initial mass density
k of an incompressible sphere depends only on the shear is approximately uniform.
2
modulusof rigidity,not on thebulkmodulus. However,the The number and distribution of particles and links in
mass-spring models approximate a material with Poisson’s a mass-spring model is mainly set by the minimum dis-
ratio ν =0.25 (Kot et al. 2014). However, Love (1911) de- tances dI,ds, however because the initial particle positions
rived a general formula for the static Love number k also arerandomlygenerated,therearelocalvariationsindensity
2
for a compressible homogeneous sphere.5 We have numer- of nodes and the spring network. The ratio ds/dI sets the
ically checked that his formula for k in the compressible mean numberof springs permass node.
2
case is only very weakly dependent on the Poisson’s ratio, Consider a spring linking the particle i residing at xi
withinbroadrangesofthebodysizeanddensity.Thismakes withaparticle j locatedin xj.Theforceduetothespring
usconfidentthatthestandardtidalformulae(derivedforthe oneachmassnodeiscomputedasfollows.Thevectorpoint-
Poisson ratio ν =0.5) can accurately describe a compress- ingfromoneoftheseparticlestoanother, xi−xj,givesthe
ible body with Poisson ratio ν =0.25 spring length Lij = |xi −xj| that we compare with the
spring rest length Lij,0. The elastic force exerted upon on
particle i by theparticle j is
3 DAMPED MASS-SPRING MODEL Feilastic =−kI(Lij−Lij,0)nˆij (26)
SIMULATIONS
where kI is the spring constant, while the unit vector is
We will compare the tidal spin down rate of a simulated nˆij = (xi−xj)/Lij. The appropriate force acting on the
viscoelasticbodytothatpredictedanalytically.Tosimulate particle j from the particle i is equal in magnitude and
opposite in direction.
5 Keepinmindthat theassumptionsofhomogeneity andcom- The strain rate of a spring with length Lij is
pthriessstihbeiolirtyycinanLboveeu’ssetdheoonrlyyaarseamnuatpuparlolxyicmoanttiroandi(cMtievlec,hwiohre1r9e7fo2r)e. ǫ˙ij = LL˙ijij,0 = LijL1ij,0(xi−xj)·(vi−vj) (27)
MNRAS000,000–000(0000)
5
modulus)andashearviscosity(inanalogytotheshearmod-
ulus). For a mass-spring model comprised of equal masses
mI andparameterisedwiththedampingcoefficients γI and
spring constants kI, 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
alsoassumethat ηI scaleswiththenetworkinthesameway
µI does in equations (30) and (29). The resulting estimate
for theviscosity is:
ηI ≈µI γImI . (31)
(cid:18) kI (cid:19)
Therelaxation timescale of a Kelvin-Voigt solid is
τ = ηI ≈ γI mI . (32)
relax µI kI
While the fidelity of the computed elastic modulus of
amass-springmodelhasbeencheckedbynumericalsimula-
tionsusingstaticappliedforces(Kot et al.2014),theviscos-
ityofadampedmass-springmodelhasnotbeentested.We
Figure 2. A snapshot of one of our simulations as viewed in consider the equations (31) and (32) as approximate. They
the open-GL viewer of rebound. We only show the tidally per- may need to be amended with factors of order unity. We
turbed body, as the perturbing one is distant from it. The ren- shall discuss this possibility later, when we compare mea-
dered spheres areshown to illustratethe random distributionof surementsfrom oursimulations topredictionsfrom anana-
node point masses, not imply that the body behaves as a rub- lytical model of tidal evolution.
ble pile (e.g., Richardsonetal. 2009; S´anchez &Scheeres 2011).
In our simulations, the body was given an initial spin
Nodesareconnected withanetworkofdampedelasticsprings. θ˙=σ perpendiculartotheorbitalplane.Thiswasdoneby
0
setting theinitial velocities of particle equal to
whereviandvj aretheparticlevelocitiesandL˙ij istherate vi = xi × σ0ˆz , (33)
ofchangeofthespringlength.Totheelasticforceactingon
the particle i, we add a damping (viscous) force propor- vi and xi being the velocity and position vectors of the
tional to thestrain rate: i-thparticle,withrespecttothebody’scentreofmass,and
the unit vector zˆ being orthogonal to the orbit. In most
Fdiamping =−γIǫ˙ijLij,0mInˆij , (28) runs,wechosetheinitialrotationtoberetrograde(σ0 <0),
to allow for a larger range of values of the tidal frequency
withadampingcoefficient γI equaltotheinversedamping χ to be simulated. This choice always led to a decrease in
time scale. The parameter γI is independent of kI. As all the semimajor axis (see Table 3) and to acceleration of the
ourparticleshavethesamemass,wedonotusethereduced
body’s spin rate. Variations in the initial particle distribu-
mass in equation (28), as did Quillen et al. (2015).
tioncreatesonlynegligiblenon-diagonaltermsintheinertia
Howdoes thespringconstant kI andthedampingpa- tensorin thiscoordinate system.
rameter γI relate totheglobal rigidity and viscosity of the Fromthereboundcode,version2asofNovember2015,
body? According to Kot et al. (2014), the static Young’s
we used the open-GL display with open boundary condi-
modulusis given bya sum over thesprings Lij,0: tions, the direct all-pairs gravitational force computation,
1 and the leap-frog integrator needed to advance particle po-
EI = 6V kIL2ij,0 , (29) sitions.Totheparticles’accelerationscausedbythegravity,
X
we added the additional spring and spring-damping forces,
where V is the total volume. For an initially random
aswasexplainedabove.Tomaintainnumericalstability,the
isotropic mass-spring system, the Poisson ratio is ν =0.25
time step was chosen to be smaller than the elastic oscilla-
(Kot et al. 2014). Our mass-spring model allows us to di-
tionfrequencyofasinglenodeparticleinthespringnetwork,
rectly estimate the Young’s modulus EI, from which we
or equivalentlythetime it takesvibrational waves totravel
can compute the shear elastic rigidity µI commonly used
between two neighbouringnodes.
for tidal evolution calculations. The relation between the
two relaxed (static) moduli is given by
µI = EI = EI , (30) 4 TIDAL EVOLUTION OF THE SEMIMAJOR
2(1+ν) 2.5
AXIS
where we set thePoisson ratio to be ν =0.25.
4.1 The semimajor axis’ tidal drift rate computed
Withnon-zerodampingcoefficients,thestressisasum
from the mass-spring model
of an elastic term proportional to the strain and a viscous
termproportionaltothestrainrate,sothemodelshouldlo- Commonparametersforourfirstsetofsimulationsarelisted
callyapproximatealinearKelvin-Voigtrheology.Themass- in Table 1, along with their chosen or computed values. A
spring model is compressible. So, when damped springs are listofvariedandmeasuredparametersispresentedinTable
used, it exhibits a bulk viscosity (in analogy to the bulk 2. The values used in different individual simulations with
MNRAS000,000–000(0000)
6 Frouard et al.
commonvaluesinTable1arelistedinTable3.Wechosethe Foreachsimulationwegenerateanewmass-springnet-
∗
valuesforthemassratio M /M andtheinitialsemi-major work. Hence each simulation has slightly different numbers
axis a to be the same in the two sets of simulations per- ofmassnodesandspringsandvariationsinthenodedistri-
0
formed. However, in each set, the spring damping rate γI butionand associated springnetwork. Two simulations run
and the initial body spin rate σ were chosen to sample a with identical input parameters will differ slightly in their
0
rangeofvaluesoftheFouriertidalmodefrequency ω andof measured semi-major axis drift rates. We set the minimum
the viscoelastic relaxation time τrelax. Using the equation distance between nodes dI and maximum spring length ds
(29),wecomputedtheYoung’smodulusbyonlyconsidering (setting the mean number of mass nodes and numbers of
thespringswithamidpointradiilessthan0.9R,andweused springs)sothatthedifferencesinsemi-majoraxisdriftrate
the volume within the same radius 0.9R. The region near for simulations drawn from thesame parameterset differed
thesurface was discarded so that theYoung’smodulus was bylessthan10%.Wewilldiscussthesensitivityofthesimu-
computedin aregion where thespringnetwork isisotropic. lationstothenumbersofnodesandspringspernodefurther
Afairly soft bodyunderstrongtidalforcing(corresponding below.
toalargemassratioandweaksprings)waschosen,toreduce The simulations were run on a MacBook Pro (early
theintegrationtimerequiredtoobserveasignificanttidally 2015) with 3.1 GHzIntelCore i7 microprocessor. Thecom-
inducedchangeinthesemimajoraxis.Atthesametime,we putation time for an individual simulation listed in Tables
made sure that the Young’s modulus was sufficiently large, 1 and 3 was about 7 minutes. If the number of mass nodes
so that the body was strong enough to maintain a nearly is doubled then the number of direct all-pairs gravity force
constant radial density profile. In the absence of exterior computations increases by a factor of 4. However the to-
pressure, the body is held up against self-gravity by spring tal computation time increases by a slightly larger factor
forces only — so the springs in the interior are under com- than 4 because the number of springs also increases (scales
pression.Inallourruns,theparticlesweredisplaced,inthe with the number of nodes) and the timestep decreases. For
bodyframe,byatmostafewpercentsoftheunitlength R. simulationswithsimilarvibrationwavespeeds,thetimestep
Wechosetheinitialsemi-majoraxis a largeenough,to scaleswiththeinterparticlespacingandsothetotalnumber
0
ensurethatthequadrupoletidalpotentialtermwoulddom- of particles to the-1/3 power.
inate. We also set the initial relative velocity of the bodies
such that the orbit would be circular. In the frame coro-
4.2 Comparison of numerically measured and
tatingwiththetidallyperturbedbody,theperturberorbits
predicted quality functions
withaperiodof Po = 2π/χ.Wecarriedouteachrunover
the time span of t = 11Po and recorded the semimajor Inverting the equation (3), we obtain the following expres-
axis’ values 10 times, with even intervals of Po. Since the sion for thequality function:
springs’lengthshadinitially beensettohavetheirrestval-
a˙ M a 5
ues, gravity caused the system to bounce at the beginning [k sinǫ ](χ˜) = − . (34)
ofeachsimulation.Dampedoscillations areexpectedin our 2 2 3na(cid:18)M∗(cid:19)(cid:16)R(cid:17)
model as each of the individual links between pairs of par- From the semimajor axis’ drift rates computed numerically
ticles is acting as a damped harmonicoscillator. The initial in our simulations, we computed the values of the qual-
oscillations are just the consequences of abruptly“turning ity function over a range of values of frequency. This was
on”self-gravity at the beginning of the simulations. Thus, performed using the equation (34) and thequantities listed
duringthefirsttimeintervalof Po,weintegrated thebody in Table 3. The numerically measured values of the quality
with a higher damping parameter, to dissipate the initial function are plotted as a function of the dimensionless fre-
vibrational oscillations. We did not take into account the quency χ˜ in Figures 3 and 4. The dimensionless frequency
semimajor axis’ value computed during that time interval. χ¯wascomputedbyusingtherelaxationtimethatwaseval-
Subsequently,werecordedthesemi-majoraxis’valuesatan uated for each simulation individually.
intervalof Po,sothattheirregularitiesoftheparticledistri- ThenumericallygeneratedpointsinFigure3showthat
butiondidnotaffectthemeasurementoftheslowlydrifting the simulated quality function is proportional to the fre-
semi-major axis. quency, at small frequencies, but decays at large frequen-
The semimajor axis was computed using the distance cies.Predictedanalyticallyforvariousrheologies(Efroimsky
between thecentre of mass of thebodies, and theirrelative 2012a;Noyelles et al.2014;Efroimsky2015),thisbehaviour
velocities.Tomeasurethesemimajoraxis’driftrate a˙,wefit has never been reproduced by direct numerical computa-
alinetothetenmeasurementsofthesemimajor axisatthe tions.Thepointsusedtobuildthisplotweremeasuredfrom
ten time intervals Po. Each fit was individually inspected, simulations with different spin rates and relaxation times.
to ensure that the ten points lay on a line. The standard Nevertheless, they all lie near the same curve, suggesting
errorofthefittedslopevaluethatprovides a˙ was .1%.For that the function is primarily dependent on the relaxation
eachsimulation,wealsocomparedtheinitialandfinalvalues time.Wethushavenumerically confirmedtheexpectedbe-
of each component of the total angular-momentum vector haviour and sensitivity of the quality function to the fre-
(relativetothecentreofmassofthebinary),andfoundthe quencyand tothe viscoelastic timescale.
absolutedifferencewasbelow 10−12 forallthecomponents. With the points obtained through simulations, we plot
Similarly, we checked that the evolution of the measured the quality function predicted analytically using equation
spinrateandsemimajoraxisweretightlyanti-correlated,as (25). We set the shear modulus µ = µI, with the value of
expected when angular momentum is conserved. The total µI given by the expression (30). This value is proportional
energy was not conserved because of the spring damping totheYoung’smoduluslistedinTable1computedfromthe
forces. springnetwork.Themassratio andmean motion aretaken
MNRAS000,000–000(0000)
7
from Table 3. The result is the grey line (the lowest one)
in Figure 3. This is the analytically predicted quality func-
tion, 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 thepredicted one
Table 1.Common simulation parameters
as qfratio and from this plot we findqfratio ∼1.3.
Figure 3 demonstrates that the numerically obtained
tidal drift rate of the semimajor axis is maximal at a fre- NI 800 Numberofparticlesinresolvedbody
quency χ˜ ≈ 1.0, which is consistent with the analytical NSI 9254 Numberofinterconnecting springs
model.Hadtherelaxationtimescale(equation32)beenmis- LI 0.2623 Meanrestspringlength
computedinoursimulations,thenumericallyobtainedpeak kI 0.08 Springconstant
EI 2.3 MeanYoung’smodulus
wouldhavebeendisplacedawayfromthatpredictedanalyt-
dI 0.15 Minimuminitialinterparticledistance
ically. The good match between predicted and numerically ds 0.345 Springformationdistance
measuredpeakfrequencysupportsourestimatefortheshear dt 0.001 timestep
viscosity (equation 31) and the associated relaxation time
(equation 32) for themass-spring model. NI, NSI, EI and LI vary slightlybetween simulations as parti-
Both the peak magnitude and the peak location of the cle distributions are randomly generated. EI is computed using
equation (29) and is the average value for all the simulations.
qualityfunctionmayshiftifhigher-orderterms(withhigher
TheseparametersarecommontosimulationslistedinTable3.
valuesof l) are included into theanalytical calculation. As
thismightexplainthedifferenceinheightofournumerically
computed quality function, compared to that predicted, we
tested this possibility by running simulations with a larger Table2.Description of varied and measuredsimulation
value of the semimajor axis. Simulations with a larger ini- parameters
tial semimajor axis are also listed in Table 3, and thequal-
ity function for these runs is plotted in Figure 4. We find a0 initialsemi-majoraxis
that the amplitude correction factor required to match the M∗/M massratio
numerical resultsis thesame as for theprevious set of sim- χ˜ Unitlessfrequency
ulations — and, again, the frequency does not need to be a˙ Rateoforbitaldecayinsemi-majoraxis
rescaled. From this,we concludethat higher-order terms in γI Springrelaxationtime
thetidalpotentialdonotexplaintheamplitudediscrepancy σ0 Initialspin
betweenournumericalsimulationsmodelandouranalytical τrelax Estimatedrelaxationtimeofviscoelasticsolid
χ Initialtidalfrequency(semi-diurnal)
predictions.
Po Timebetweenintegrationoutputs
Since we are using a random mass-spring model, both EI ComputedYoung’smodulus
the particle distribution and the spring network differ be-
tween simulations. We computed the Young’s modulus and The viscoelastic relaxation time τrelax is computed using the
equations (29) and (32). The frequency χ is defined by the ex-
relaxation time for each run, and these are listed in Table
pression (A14), while the tidal forcing frequency χ˜ is given by
3.Wefoundsmallvariationsintheelasticmodulusbetween
the equation (17). The period Po = 2π/χ is also a part of the
different simulations. In Figures 3 and 4, we also show the
simulationoutput.
analytically predicted quality function factored by 1.3 and
offset by ±10% (green curves). Points with higher values
of EI, corresponding to harder bodies, systematically lie mass nodes, the second 10 (columns C, D) about 800 mass
lower than the appropriate points with lower values of EI nodesandthelast10(columnsE,F)about1600massnodes.
corresponding to softer bodies experiencing stronger tidal The first 5 in each group of 10 have about 10 springs per
deformation.Thescatterinourpointsaboveandbelowthe node (columns A, C and E) whereas the second 5 of each
blue line can be attributed to variations in the particle dis- groupof10(columnsB,DandF)haveabout20springsper
tribution and in theassociated springnetwork. node. The spring constants, kI, and damping parameters,
γI,foreachsetwerechosensothatχ¯≈0.225,EI ≈2.3and
τ ≈0.155. Foreach simulation separatelywecomputed
4.3 Discrepancy in amplitude of quality function relax
and sensitivity to the number of masses and ratio qfratio of the numerically measured value of the qual-
ity function k sinǫ to the predicted one computed using
springs simulated 2 2
thenormalizedfrequencyχ¯andYoung’smodulusmeasured
To see how the discrepancy in amplitude is related to the fromthenodesandspringsinthesimulations.Wethencom-
numberof mass nodesand springs, we carried out asecond puted the mean and standard deviation of qfratio for each
series of simulations each with the same initial spin, esti- setof5simulationswithidenticalrunparametersandthese
matedelasticmodulusandviscoelasticrelaxationtimescale, are listed in the bottom rows of Table 4. The standard de-
buthavingdifferentnumbersofnodemassesandnumbersof viation in semi-major axis drift rate divided by the mean
springspernode.Parametersforthesesimulationsarelisted value of the drift rate is also computed for each set of 5
in Table4. Weran 6sets, each of5 simulations, all approx- simulations and listed in Table 4.
imately matching the third row in Table 3 with a = 10 In the simulations listed in Table 1 and 3, the ratio
0
, M∗/M = 100, σ0 = −0.4, and Po = 4.377. The 5 simu- of the number of springs per node is about 11, which is
lations in each set have identical run parameters. The first slightly lower than the number recommended by Kot et al.
10 simulations (columns A, B in Table 4) have about 400 (2014,Figure5).Forthisratio,Kot et al.(2014)foundthat
MNRAS000,000–000(0000)
8 Frouard et al.
One possible reason for the discrepancy between pre-
dicted and measured quality function is that near thebody
surface thenumberof springs perparticle is lower than the
interior,andthusthespringnetworkisanisotropicnearthe
Table 3. Quantities either set or computed in different
surface. The simulated body is weaker (and floppier) than
mass-spring N-body simulations
we estimated by integrating the spring properties over the
volume. Because the overall rigidity of the body is weaker,
χ˜ a˙ γI σ0 τrelax χ Po EI
its tidal response would be larger than we predicted ana-
witha0=10,M∗=100,n=0.318 lytically. This would be consistent with our greater than
0.069 -3.897e-05 20 0.1 0.158 0.436 14.424 2.30 unity qfratio. If this were the primary cause of our am-
plitude discrepancy then the ratio of numerical computed
0.162 -8.505e-05 20 -0.2 0.156 1.036 6.067 2.35
toanalyticalpredictedqualityfunctionshouldinverselyde-
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 pend on the number of simulated masses. Table 4 shows
0.645 -2.309e-04 50 -0.5 0.394 1.636 3.841 2.33 that this ratio qfratio does decrease as the particle number
0.896 -2.517e-04 80 -0.4 0.624 1.436 4.377 2.36 is increased from 400 to 1600, however it only decreases by
1.123 -2.288e-04 100 -0.4 0.782 1.436 4.377 2.36 about a percent between 800 and 1600 particles. Conver-
1.467 -2.325e-04 100 -0.6 0.799 1.836 3.423 2.26 gencehasnotbeenachieved,butthestandarddeviationsin
witha0=20,M∗=200,n=0.159 quantities measured from groups of simulations imply that
thesimulation results are reproducibleand that thescatter
0.080 -2.239e-06 20 -0.1 0.156 0.517 12.153 2.40 is smaller than the amplitude discrepancy. If doubling the
0.160 -4.540e-06 40 -0.1 0.309 0.517 12.153 2.41 particle number reduces the quality function ratio by 2%
0.227 -7.124e-06 40 -0.2 0.317 0.717 8.763 2.29
thentoreach a qualityfunction ratio of 1 wewould require
0.350 -8.601e-06 40 -0.4 0.314 1.117 5.625 2.36
a million particles in the simulation and we are not yet set
0.534 -1.157e-05 60 -0.4 0.478 1.117 5.625 2.25
upto run simulations this large.
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
Thesurfaceofourbodylieswithinaradiusof1,conse-
1.188 -1.372e-05 100 -0.6 0.783 1.517 4.142 2.31
quently our simulated body effectively has an outer radius
1.771 -1.197e-05 150 -0.6 1.167 1.517 4.142 2.28
that is smaller than 1. As the semi-major axis drift rate is
Therightmostcolumncontains thevalues oftheYoung’s modu- fasterforlargerbodies,wewouldhaveexpectedoursimula-
lus,computedforeachsimulationbymeansoftheequation(29). tionstoexhibitslowerratherthanfasterorbitaldecayrates
The simulations listed here have common parameters listed in (andevidentfromthefactorofR−5 inequation34).Were-
Table1.
callthattheanalyticallypredictedqualityfunctiondepends
on the ratio µ/eg (see equations 24 and equation 25). As
the spring network behaved 10% weaker 6 than what was theenergydensityeg ∝R4,foran effectiveradiusless than
estimated by equation (29).In Table 4 we can compare the 1,theenergydensityishigherthanwepreviouslyestimated
ratioofnumericallymeasuredtopredictedqualityfunctions andtheratioµ/eg wouldbelowerthanweused.Thismeans
for simulations with similar numbers of nodes but differ by we have underestimated thesize of the tidal response. This
the number of springs per node. Comparing columns A to isintherightdirectiontoaccountfortheamplitudediscrep-
BandCtoDandEtoFweseethattheratioofnumerical ancy.HoweverifwetakeintoaccountboththefactorofR5
measured to predicted quality function is only slightly less inequation 34when computingournumericalqualityfunc-
(about 2% lower) when the number of springs per node is tionandthefactorofR4inequation25whencomputingour
about 20 rather than 10. We find that for greater than 10 analyticalqualityfunction,thetwocorrections moreorless
springspernodes,theratioofmeasuredtopredictedquality cancel each other and we do not resolve the source of dis-
functionisrelativelyinsensitivetothenumberofspringsper crepancy in quality function amplitude. So we find that we
node. cannot resolve our discrepancy in amplitude by correcting
Aseachsimulationgeneratesanewparticledistribution thebody radiusto an effectively smaller radius.
andspringnetwork,wecanmeasureeffectsduetovariations
Equations(29-30)characterisestaticpropertiesofthe
inthesepropertiesbycomparingsimulationsgeneratedfrom
mass-springmodel.However,theanalyticaltheoryofevolv-
identicalinputparameters.Table4showsthatthestandard
ing tides (including that based on the viscoelastic rheol-
deviation of the quality function ratio and semi-major axis
ogy 15) employs the unrelaxed shear rigidity µ. We have
driftratedecreaseswithincreasedparticlenumberandthat
assumed that the two values are equivalent, although this
thescatterinthesemeasuredquantitiesdiffersbyonlyafew
is not perfectly true for real materials. As any rheological
percent. We conclude that variations in the spring network
parameter, the shear elasticity modulus depends upon fre-
are unlikely to explain the discrepancy between predicted
quency. The unrelaxed or frequency dependent shear mod-
and numerically measured quality function.
ulus in real materials can be lower than the relaxed or
staticcounterpart(e.g.,Faul& Jackson2005).Perhapsthis
is also true in our simulated material. The numerical tests
6 Because of a misprint, Figure 5 in Kotetal. (2014) actually
byKot et al.(2014)usedstaticforcesandsowouldnothave
displays (E0−E)/E0 insteadof (E−E0)/E0,with E0 and E
been sensitive to frequency dependence in the shear modu-
beingtheestimatedandmeasuredYoung’smodulus,respectively
(Kot et al., private communication). A smaller average number lus.Hadweemployedasmallervaluefortheshearmodulus,
hSi of springs per node gives a smaller Young’s modulus and a thiswouldhaveincreasedthevaluesofthetheoreticallyob-
weakerbody. tained k2 sinǫ2 , and thus would have reduced the offset
MNRAS000,000–000(0000)
9
betweenournumericalandtheoreticalvaluesforthequality
function.
Each numerical run starts out with a sphere of parti-
cles, with springs at their rest lengths. However, after the
simulation begins, the body is compressed by self-gravity.
Theresultingdensityprofileisnotperfectlyflat,asthecen-
trebecomesmorecompressed thantheregionsclosertothe
surface. Furthermore,theinitial spin of thebody causes its Table4.Comparisonof simulations with different num-
shapetodeviatefrom asphere.Tosampleabroad rangeof bers of node masses and springs per node
frequencies (in units of the relaxation time), and to have a
strongertidalresponse(i.e.,toreducethesimulation time), witha0=10,M∗=100,n=0.318,σ0=−0.4,Po=4.377
weuseasoftbodythatisparticularlypronetodeformation andχ¯∼0.225,EI ∼2.3,τrelax∼0.155
whenspunandcompressedbyself-gravity.Neverthelessour
Set A B C D E F
simulated bodies are not strongly deformed and we do not
hNIi 403.8 405.8 795.2 804.2 1597.8 1596.2
expectbodydeformationandcompressiontoaccountforour
hNSI/NIi 10.92 18.7 11.5 20.34 12.02 21.34
quality function amplitudediscrepancy. dI 0.19 0.19 0.15 0.15 0.118 0.118
Since our simulated body is comprised of randomly- ds/dI 2.3 2.8 2.2 2.8 2.3 2.8
distributedpointparticles,thebodyisneitherexactlyspher- kI 0.102 0.04 0.08 0.0301 0.062 0.023
icalnoruniform.Initially,itsprincipalaxesofinertia(com- γI 12.858 5.05 20.0 7.48 31.3 11.5
puted from the moment of inertia tensor) are oriented ran- hqfratioi 1.48 1.39 1.41 1.36 1.40 1.38
domly, and the three moments of inertia are not exactly σ[qfratio] 0.069 0.052 0.053 0.026 0.030 0.012
σ[a˙]/ha˙i 0.061 0.064 0.055 0.024 0.026 0.0094
equal. So, initially, the body’s spin angular momentum is
not exactly perpendicular to the orbit. In the course of the Each column represents a group of 5 simulations. From each
simulations, we measured the values of the x and y com- group of 5 the mean number of nodes and springs per node are
ponentsof thespin angular momentum,findingthem to be listed in the second and third rows. The rows labelled hqfratioi
a few hundredthsof the initial z component. This is small and σ[qfratio] give the mean and standard deviations, respec-
enough thatthespin aboutnon-polaraxes isunlikelyto be tively, of the ratio of the numerically measured to analytically
predicted quality function. The bottom row gives the ratio of
thecause of our quality function amplitudediscrepancy.
the standard deviation in the semi-major axis drift rate divided
Tosummarize, wehavecompared simulations with dif-
by the means. Each of the means and standard deviations are
ferent numbers of mass nodes and springs per node and
computed from 5 simulations. The simulations listed here differ
foundthatthescatterandsensitivityoftheratioofnumeri-
fromthosedescribedbyTables1and3.
callymeasuredtoanalyticallypredictedqualityfunctionare
onlyweaklydependentonthesequantities.Convergencehas
notbeenreachedat1600simulatedparticlesbutthescatter axisandthespinratewereadirectoutcomeofthedamped
in the simulations is low enough that the variations in the mass-springmodel.Wecomparedtheresultsobtainedbythe
spring network are small enough that we have reproducible twomethods,andconcludedthatamass-springnetworkcan
numerical measurements and scatter within these measure- serve as a faithful model of a tidally deformed viscoelastic
mentaresmallerthanourmeasuredamplitudediscrepancy. celestialbody.Specifically,wecomputedthequalityfunction
We have not identified the source of the 30% discrepancy (theratioofLovenumberk andthequalityfactorQ)from
2
in theamplitude of our numerically measured quality func- thenumerically simulated tidal drift of the semimajor axis.
tion, however we suspect variations in the spring network, Thequalityfunctionshowedastrongfrequency-dependence
floppiness in the outer surface and a possible frequency de- that is close to the dependence derived analytically for a
pendenceinthebehaviorofthenumericallysimulatedshear Kelvin-Voigt sphere by a method suggested by Efroimsky
modulus. (2015).
While direct estimates for the global shear and bulk
rigiditycanbederivedfromthespringconstantsandspring
lengths, the shear viscosity has not been estimated in the
5 CONCLUSION
literature.Consequently,wewereuncertainofourestimates
In this article, we have used a self-gravitating damped for the shear viscosity and the associated relaxation time
mass-springmodel,withinanN-bodysimulation,todirectly (equation 32). However, a comparison between our com-
model the tidal orbital evolution and rotational spin-up of putedqualityfunction(specificallythepeakfrequency)and
a viscoelastic body. thatpredictedanalytically foraKelvin-Voigtsolid suggests
Weconsideredabinarycomprisedofanextendedbody that our estimates for the numerically simulated shear vis-
(assumed spherical and homogeneous) and a point-mass cosity and viscoelastic relaxation time were correct.
companion.Withinthissetting,wehavetestedournumeri- The magnitude of our numerically predicted quality
calapproachagainstananalyticalcalculationforthissimple function is about 30% larger than the one predicted ana-
case. The semimajor axis’ tidal evolution rate was calcu- lytically. By comparing simulations performed for two dif-
lated by two methods. One was a direct simulation based ferent initial values of the semi-major axis, we concluded
on simulating the first body with a mass-spring network. thatthecauseisnottheneglectofhigherordertermsinthe
Another, analytical, method was based on a preconceived quadrupoleexpansionforthepotential.Wefindtheratioof
viscoelasticmodeloftidalfriction(theKelvin-Voigtmodel). numerically measured to predicted quality function is only
Thenumericallycomputedtidalevolutionofthesemi-major weakly dependent on numbers of mass nodes and springs
MNRAS000,000–000(0000)
10 Frouard et al.
(cid:1)(cid:1)(cid:2)(cid:1)(cid:6)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:5)(cid:3)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:5)
(cid:1)
(cid:5)
(cid:1)(cid:4) (cid:1)(cid:1)(cid:2)(cid:1)(cid:4)(cid:3)
(cid:2)(cid:3)
(cid:1)(cid:1)
(cid:1)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:4)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:1)(cid:3)
(cid:10)(cid:11)(cid:4)(cid:2)(cid:6)
(cid:1)(cid:1)
(cid:1)(cid:1) (cid:1)(cid:1)(cid:2)(cid:5) (cid:1)(cid:1)(cid:2)(cid:7) (cid:1)(cid:1)(cid:2)(cid:8) (cid:1)(cid:1)(cid:2)(cid:9) (cid:1)(cid:4) (cid:1)(cid:4)(cid:2)(cid:5) (cid:1)(cid:4)(cid:2)(cid:7) (cid:1)(cid:4)(cid:2)(cid:8)
(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:3)(cid:6)(cid:7)(cid:8)(cid:1)(cid:9)(cid:10)
Figure3.Acomparisonofthenumericallycomputedqualityfunction,shownaspoints,tothatcalculatedanalyticallyforahomogeneous
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
shearmodulusby10%intheoffsetanalyticalcalculation.Boththenumericalandanalyticalcalculationswerecarriedoutforaperturber
ofmass M∗=100 andinitialvalueofthesemi-majoraxis a0=10.Wefindthatthenumericallycomputedqualityfunctionhasashape
andpeakfrequencyconsistentwiththatobtainedanalytically,buttheamplitudeistoohighbyabout30%.Weattributethescatterof
the points offthe linetovariations inthe value ofthe shear modulus of therandom mass-springnetwork due tonon-uniformityinthe
particledistributionandspringnetwork(seesection4.2).
(cid:1)(cid:1)(cid:2)(cid:1)(cid:6)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:5)(cid:3)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:5)
(cid:1)
(cid:5)
(cid:1)(cid:4) (cid:1)(cid:1)(cid:2)(cid:1)(cid:4)(cid:3)
(cid:2)(cid:3)
(cid:1)(cid:1)
(cid:1)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:4)
(cid:1)(cid:1)(cid:2)(cid:1)(cid:1)(cid:3)
(cid:10)(cid:11)(cid:4)(cid:2)(cid:6)
(cid:1)(cid:1)
(cid:1)(cid:1) (cid:1)(cid:1)(cid:2)(cid:5) (cid:1)(cid:1)(cid:2)(cid:7) (cid:1)(cid:1)(cid:2)(cid:8) (cid:1)(cid:1)(cid:2)(cid:9) (cid:1)(cid:4) (cid:1)(cid:4)(cid:2)(cid:5) (cid:1)(cid:4)(cid:2)(cid:7) (cid:1)(cid:4)(cid:2)(cid:8) (cid:1)(cid:4)(cid:2)(cid:9)
(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:3)(cid:6)(cid:7)(cid:8)(cid:1)(cid:9)(cid:10)
Figure 4. The same as in Figure 3, except that the perturber’s mass is M∗ = 200, and the initial value of the semi-major axis is
a0=100.
per node simulated but does slowly decrease with an in- than the orbital time scale. Nonetheless, mass-spring mod-
crease in the numbers of particles simulated. We have not els can be used to explore the tidal evolution of inhomoge-
yet identifiedthecause of thediscrepancy.Wesuspect that neousandanisotropicbodies,andtostudyhowtheirquality
thenon-uniformityofthespringnetworkatthebodysurface functions depend on the rheological properties and internal
couldhavecausedalargerthanexpectedtidalresponse.Al- structure of the body. This fully numerical approach also
ternativelythesimulatedshearmoduluscouldbefrequency permits study of tidal phenomenathat are not easy to pre-
dependentand overestimated by its computed static value. dict analytically — such as capture into (and crossing of)
spin-orbit and spin-spin resonances, tidally induced orbital
Our study demonstrates that we can directly simu-
evolution, the distribution of tidal heating, and the effects
late the tidal evolution of viscoelastic bodies. We currently
of non-linearrheological behaviour.
achieveanaccuracy of30% butwith ongoing effort wemay
improve upon this. It would be difficult to model this pro-
cess over long (Myr or Byr) timescales, because the time
step is determined by the number of particles within the
body and by the spring constants for the links connecting
theparticles.Thisrequiresthetimesteptobemuchshorter
MNRAS000,000–000(0000)