On the elastodynamics of rotating planets

Equations of motion are derived for (visco)elastic, self-gravitating, and variably-rotating planets. The equations are written using a decomposition of the elastic motion that separates the body's elastic deformation from its net translational and rotational motion as far as possible. This separation is achieved by introducing degrees of freedom that represent the body's rigid motions; it is made precise by imposing constraints that are physically motivated and should be practically useful. In essence, a Tisserand frame is introduced exactly into the equations of solid mechanics. The necessary concepts are first introduced in the context of a solid body, motivated by symmetries and conservation laws, and the corresponding equations of motion are derived. Next, it is shown how those ideas and equations of motion can readily be extended to describe a layered fluid--solid body. A possibly new conservation law concerning inviscid fluids is then stated. Thereafter the equilibria and linearisation of the fluid--solid equations of motion are discussed, along with new equations for use within normal-mode coupling calculations and other Galerkin methods. Finally, the extension of these ideas to the description of multiple, interacting fluid--solid planets is qualitatively discussed.


INTRODUCTION
The continuum mechanics of rotating, self-gravitating, viscoelastic, layered fluid-solid planets is relevant in studies of: seismology, co-and post-seismic deformation, Earth rotation, tides and glacial isostatic adjustment (GIA). Since the underlying equations can be formulated in different ways, it is natural to ask how best to do so. The standard approach as summarised in Chapter 3 of Dahlen & Tromp (1998, 'D&T' hereafter) is to write the exact equations of finite elasticity (e.g. Coleman & Noll 1961;Coleman 1964;Biot 1965;Noll 1974;Truesdell & Noll 2004) with respect to a uniformly rotating reference frame whose origin is fixed in inertial space, to take a referential ('Lagrangian') view of particle motions and a spatial ('Eulerian') view of the gravitational potential, and then to linearise about equilibrium. This approach can be seen as a culmination of geophysical work that began in the late 1960s (Backus 1967;Dahlen 1972a) and was continued through the 1970s by Dahlen (1972b), Smith (1974, Dahlen & Smith (1975) and Woodhouse & Dahlen (1978). In this paper we present a different approach: we formulate the equations of continuum mechanics using an entirely referential description and with respect to a dynamically evolving reference frame. This idea has been directly inspired by work on Euler-Poincaré reduction within geometric mechanics (e.g. Holm et al. 2009;Marsden & Ratiu 2013). Our approach clarifies the physical unity of Earth's rotation and elastodynamics, providing a general framework that bridges the respective theoretical approaches within geodesy and seismology and that generalises both (c.f. Smith 1977, p. 104). We also emphasise the exact nonlinear formulation of the problem, in contrast to earlier studies that focused on linearised theories.
Before discussing the novel aspects of this work, it is worth emphasising that the final equations of motion we obtain take a broadly similar form to, and are not appreciably more complex than, those of conventional formulations (e.g. D&T, Chs. 2&3). For example, we find a momentum equation almost identical to D&T's, just with a dynamical angular velocity and associated Euler force, and subject to some constraints. This momentum equation is coupled to simple ODEs governing the evolution of the reference frame. There is therefore no real penalty incurred in using our reformulation, with near-identical numerical methods being applicable in both cases. On the other hand, as will be explained below, our equations have a greater range of applicability while offering certain conceptual advantages.
Our treatment of rotation is based upon an exact implementation of the Tisserand frame On the elastodynamics of rotating planets 3 (e.g. Munk & McDonald 1960) within continuum mechanics. By definition, the deformation of a planet viewed from such a frame possesses no net linear or angular momentum. Instead, the evolution of the frame captures any translational and rotational aspects of the planet's motion.
Although discussion and use of such a frame is not uncommon within observational studies, we are not aware of its exact implementation within the equations of solid mechanics (though see Rose & Buffett (2017) for comparable ideas in the context of mantle convection). In contrast to conventional formulations based on a steadily rotating reference-frame, the Tisserand frame simplifies the study of rotating bodies' elastodynamics. For example, consider how the melting of an ice-sheet will affect a planet's motion. Not only will the change in surface loading cause (visco)elastic deformation but, through changes in the planet's moment-of-inertia, there will also be an associated rotational perturbation. Viewed from a steadily rotating reference frame, such a rotational perturbation manifests as a secular growth in the displacement. From the Tisserand frame, however, that perturbation is absorbed by the frame's motion and there is no secular component to the planet's internal deformation. More generally, as discussed by Canavin & Likins (1977), the Tisserand frame is preferable to a steadily rotating frame for bodies undergoing large rotational variations because the latter frame can make small internal deformations look larger than they are. Within conventional geophysical applications there is perhaps no great need to consider large rotational perturbations, but we reiterate that the use of the Tisserand frame does not entail major complications. Furthermore, such large perturbations are relevant within certain planetary science problems (e.g. Ćuk et al. 2016;Lock et al. 2018), and there the use of a steadily rotating frame would be inconvenient within an exact theory and could lead to potentially significant errors upon linearisation.
As noted above, the standard geophysical approach to continuum mechanics mixes spatial and referential variables. The reason for this mixed formulation seems to be largely historical, perhaps motivated by a desire to obtain familiar-looking equations. Here, in contrast, we use a consistently referential description of the problem. While either approach is valid, we feel that a fully referential treatment is clearer even if initially less familiar. In the context of numerical work, a spatial formulation is often just as convenient as a referential one, but for certain methods (e.g. Maitra & Al-Attar 2019; Leng et al. 2019) a referential formulation is essential.
A notable feature of this paper is its emphasis on the exact nonlinear equations of continuum mechanics. We develop such equations in a form suitable for numerical implementation, taking the view that nonlinear rather than traditional linearised equations will be most useful in future numerical work. To give an example, recent work on modelling GIA has focused on the effect of lateral variations in Earth structure, including that of realistic surface topography (e.g. Latychev et al. 2005). The amplitude of such topography is, however, comparable to the surface deformation produced close to continental ice sheets, and so one could ask whether the use of linearised theory is appropriate. A similar point could be made in the context of coand post-seismic deformation modelling (e.g. Gharti et al. 2019) if the cumulative effects of many earthquake cycles are considered. Importantly, approaching such problems using the full nonlinear equations should not entail unnecessary complications because with modern numerical methods the cost of solving a weakly nonlinear problem is similar to that of the corresponding linearised problem. Thus, if nonlinear effects do not matter then one might as well solve the nonlinear problem, while if such effects do matter then one is obliged to solve the nonlinear problem.
This paper places particular emphasis on symmetries and associated conservation laws.
For this reason we start from Hamilton's principle of Least Action and use Noether's theorem to derive conservation laws. The Euclidean symmetries that underlie conservation of linearand angular-momentum directly motivate our introduction of the Tisserand frame. We also discuss an infinite-dimensional symmetry-group associated with inviscid fluid regions, and use this to derive a novel conservation law. It is the latter symmetry that allows us to propose a new, fully referential approach to modelling quasi-static deformation in planets with a fluid core (c.f. Dahlen 1974). As a special case, we show that the standard potential formulation for the displacement within neutrally-stratified fluid regions (e.g. Komatitsch & Tromp 2002) can be derived from the conservation law. We anticipate that further study of this conservation law could help address complications related to core undertone modes (e.g. Crossley 1975;Valette 1989;Rogister & Valette 2009).
An apparent drawback of using Hamilton's Principle is that it is only valid for nondissipative materials. Nevertheless, Hamilton's Principle leads to equations of motion of precisely the same form as do derivations starting from, say, balance-laws or a weak formulation.
The only practical difference between the two approaches is that the latter does not require the first Piola-Kirchhoff (FPK) stress to be derived from a strain-energy function. Instead, the existence of the FPK stress follows from Cauchy's theorem (e.g. Marsden & Hughes 1994, p.127) and constitutive behaviour is imposed by relating the FPK stress to the deformation gradient history (e.g. Coleman & Noll 1961). On a pragmatic level, all that this means is that we can use this paper's equations of motion to model dissipative bodies, but we must just remember that the stress is not related to a strain-energy function.
The paper is organised as follows. We begin in Section 2 by considering solid bodies. Then in Section 3 we show how Section 2's methods can readily be extended to fluid-solid bodies On the elastodynamics of rotating planets 5 through the work of Al-Attar et al. (2018, 'A18' henceforward). This presentation allows us to introduce the Tisserand frame in a simple context before considering complications due to fluid-solid boundaries. In Section 4 we discuss qualitatively how the methods developed for a single fluid-solid body can be generalised to problems of multiple interacting bodies. We give two examples: first we discuss elasticity's exact effect on the Keplerian two-body problem; secondly, we suggest a new, exact approach to the Earth's quasi-rigid motions. Finally, in Section 5 we sum up the preceding work and speculate on paths forward.

ROTATING SOLID ELASTIC PLANET
Here we show how Hamilton's Principle is used to derive the equations of motion and study the symmetries of a solid elastic planet. In Section 2.1 we state this paper's initial version of Hamilton's principle, then discuss symmetries and conservation laws in Sections 2.2 and 2.3. Section 2.4 introduces this paper's main idea of separating elastic from rigid motion, then in Sections 2.5 and 2.6 we restate the action accordingly and derive new equations of motion.
Unless otherwise stated, Marsden & Hughes (1994) is our source for standard results from continuum mechanics.

Kinematics
To describe the deformation of an elastic body we require the concepts of a reference body and a configuration. A reference body is a compact subset B ⊆ R 3 with open interior and smooth boundary, ∂B, whose points correspond bijectively to particles of the body. The form of the reference body encodes topological information (e.g. the body's connectedness) but is otherwise arbitrary. The possible states of the body are described by mappings from B into R 3 , with each such mapping known as a configuration. We require that these configurations are smooth embeddings of B into R 3 , this meaning that such a mapping, φ, is smooth, injective, and that there exists a smooth inverse mapping φ −1 : φ(B) → B defined on its image. Physically, these conditions mean that the body cannot be torn nor interpenetrate during its motion. We will write the set of all such embeddings as Emb B, R 3 , with this set forming the configuration space within our problem. Without dwelling on any technical details, we note that Emb B, R 3 forms an infinite-dimensional manifold (e.g. Abraham et al. 2012); this geometric point of view is sometimes useful. In Section 2.1.3 we will slightly modify this picture so as to facilitate a referential approach to self-gravity. Further changes are also necessary in Section 3.2 when we consider internal boundaries within the body across which the material parameters fail to be smooth. In both cases, the result is that the configuration space needs to be suitably generalised, but the following definitions remain valid.
As time progresses, the body adopts different configurations. Its motion can then be viewed as a time-parametrised curve within the configuration space, Emb B, R 3 , that we denote by The configuration φ(t) carries the particle labelled by x ∈ B to the point y = φ(t)(x) in physical space. We can equivalently consider the motion as a mapping which returns the position of the particle x at time t. By a harmless abuse of notation, we will also denote this latter mapping, by φ, and hence write φ(x, t) = φ(t)(x). Here we see that the motion can be viewed as either (i) a function of time that takes values within a certain function space, or (ii) a function of space and time that takes values in R 3 . Both perspectives are useful in different contexts, and the same holds for some other quantities introduced below.
The velocity, v, is the time-derivative of the motion, with v(t) taking its values in the tangent space to Emb B, R 3 based at φ(t). Note that because R 3 is flat, this latter space can be identified with the space of vector fields on B. The velocity of a particle, x, at time, t, can be written as v(t)(x) or equivalently, and more conventionally, as v(x, t).
For a configuration, φ ∈ Emb B, R 3 , we can form its derivative, Dφ, which is defined through the first-order Taylor expansion Here we see that this derivative takes values in the space of linear operators on R 3 . More accurately, Dφ(x) is linear mapping between the tangent spaces T x B and T φ(x) R 3 , but this distinction can be ignored because R 3 's being flat allows tangent vectors based at different points to be trivially identified. Due to its importance within continuum mechanics, the derivative of a configuration is denoted by the symbol F and is known as the deformation gradient. This definition extends trivially to a motion, φ : I → Emb B, R 3 , with the time-t deformation gradient, F(t), being the derivative of the instantaneous configuration, φ(t). The value of the deformation gradient at a particle x and time t can be written either as F(t)(x) or, as will generally be done, F(x, t). Finally, the Jacobian of a configuration or motion is then On the elastodynamics of rotating planets 7 defined as J(x, t) = det F(x, t).
(2.4) Due to our assumption that φ(t) has a smooth inverse for each t ∈ I, it follows from the inverse function theorem (Marsden & Hughes 1994, p.31) that F(x, t) takes values in the general linear group GL(3) (Appendix A1). We assume without loss of generality that J is everywhere positive, meaning that the motion is orientation preserving.

Variational principle
To write down a variational principle for the body, we require expressions for its kinetic and potential energy. Recall that due to conservation of mass (Marsden & Hughes 1994, p.87, Theorem 5.7), the referential density, ρ, can be defined through where ϱ(y, t) denotes the density within the instantaneous body, B t . The kinetic energy at time t is given by Meanwhile, we subdivide the body's potential energy into two parts: self-gravitational (V G ) and elastic (V E ). The former is just where the gravitational potential, ϕ, satisfies the usual Poisson equation (e.g Landau & Lifshitz 1971;Dahlen & Tromp 1998). As for the elastic potential energy, we have where W denotes the body's strain-energy function. The form of W is constrained by the principle of material-frame indifference (e.g. Marsden & Hughes 1994;Truesdell & Noll 2004), which requires that for all rotation matrices Q ∈ SO(3) and F ∈ GL(3) (see Appendix A1). This implies that the strain-energy depends on F through the right Cauchy-Green tensor so that the strain-energy function takes the form (2.11) Following Al-Attar & Crawford (2016, 'AC16' hereafter), we will include a stress-glut representation of a seismic source by allowing the strain-energy function to have the explicit where the stress glut S is a symmetric second-order tensor field.
All in all, the action describing the elastic body is the integrated difference of the body's kinetic and potential energies: with ϕ understood to satisfy Poisson's equation (defined below, eq. 2.14) and where we have pulled V G back to the reference body. Hamilton's principle states that the motion of the body is a stationary point of the action subject to fixed endpoint conditions (e.g. Marsden & Hughes 1994, Ch. 7).

Accounting for gravity variationally
There are different ways of incorporating into the variational principle the constraint that ϕ satisfy the Poisson equation. One approach is to solve the Poisson equation explicitly, using Green's functions to express ϕ in terms of the motion and thus allowing the action to be seen as a functional of φ alone within a purely referential picture (e.g. AC16). Alternatively, one can follow Woodhouse & Deuss (2015) and incorporate Poisson's equation into the problem via Lagrange multipliers, thus treating both the motion and the gravitational potential on an equal footing. Within this work we find the latter approach preferable, showing, in contrast to Woodhouse & Deuss (2015), how it can be implemented within a fully referential framework.
with ψ now viewed as a time-dependent Lagrange multiplier field. We then vary the action with respect to ϕ to give which is a weak-form PDE for ψ whose unique solution is by inspection Substituting this back into the action gives As noted by Woodhouse, variation of this action yields the correct equations of motion. We have thus eliminated the Lagrange multiplier to give a 'reduced' variational principle that has φ and ϕ as its independent variables. Action (2.18) still mixes the spatial with the referential through its dependence on ϕ.
Within the reference body, we can define a referential potential, ζ = ϕ • φ (e.g. AC16), but this definition does not apply within the exterior domain. Our solution is simply to extend the domain of φ. Instead of taking φ(t) to be an embedding of B into R 3 , we take it to be a diffeomorphism that maps R 3 onto itself. For later reference, we note that a diffeomorphism is a smooth mapping with a smooth inverse. The collection of such mappings defined on a smooth manifold, M , will be written Diff (M ) and forms a group under composition. We now see that the instantaneous configurations of the body lie in the diffeomorphism group of Euclidean space, Diff R 3 , while the motion is a time-parametrised curve taking values within this group.
In Section 2.2.1 we will show that such an extension of the motion is well-defined.
With the domain of the motion so-extended, we can now define the referential potential globally by along with the tensor (see Maitra & Al-Attar 2019) and rewrite action (2.18) as ⟨∇ζ, a · ∇ζ⟩ dV dt. (2.21) This form of the action is central to the rest of the paper.

Some symmetries of the action
We just extended φ arbitrarily from B to the whole of space, a manoeuvre that might raise eyebrows. To see that this is a legitimate step, we need to establish the invariance of the action with respect to a certain transformation group. We will also take this opportunity to discuss other symmetries of the action, which will lead to a discussion of conservation laws and ultimately motivate our introduction of the Tisserand frame.

Particle-relabelling transformations (or, 'Right-actions of the diffeomorphism group')
Consider a body with motion φ and referential potential ζ. Given a diffeomorphism, ξ ∈ Diff R 3 we can define new fields by This mapping defines a right-action of the diffeomorphism group, and we will denote it by (2.23) In the terminology of AC16 and A18, this right-action is known as a particle-relabelling transformation. From the kinematic identities discussed above, ξ's right-action on φ induces transformations of the problem's other variables, for example Starting from the action in eq. (2.21), we can use ξ to change variables to obtain ∇ζ,ã · ∇ζ dV dt, (2.25) where for convenience we have set (φ,ζ) = (φ, ζ) · ξ,ρ = ρ · ξ,W = W · ξ andB = ξ −1 (B), and where we have introduced the right-actions for arbitrary (x, t) ∈ R 3 × I and F ∈ GL(3). This argument has not, in general, demonstrated an invariance of the action because it has been necessary to transform the material parameters (ρ, W ); we will return later to the question of whether these parameters can be left unchanged by the action of a non-trivial diffeomorphism. Nevertheless, if we restrict ξ such that it equals the identity on an open subset containing B then the material parameters are unaltered and we have identified a symmetry of the action. Denoting by Diff B R 3 the subgroup of such diffeomorphisms, we conclude that if (φ, ζ) is a solution of the equations of motion, then so is (φ, ζ) · ξ for any ξ ∈ Diff B R 3 . These two solutions agree within the reference body, while in the exterior domain they are such that the spatial potential ϕ = ζ • φ −1 is uniquely defined.
Here we see an example a gauge invariance, with the solution of the problem only being defined up to the action of a certain transformation group, but with all physically meaningful quantities being unaffected by such transformations. We can, in solving the problem, fix these additional degrees of freedom in any manner that is convenient. In more formal terms, we have extended the configuration space of the body to Diff R 3 , but have shown that the dynamics is defined only up to an element of Diff B R 3 . This means that a unique solution can be found within the quotient space Diff R 3 /Diff B R 3 which precisely equals the initial configuration space Emb B, R 3 .

Left-actions of the diffeomorphism group
One can also carry out transformations where the diffeomorphism group acts from the left.
Similarly to right-actions, such left-actions of Diff R 3 are not generally a symmetry of the action. But we do have a symmetry when acting with the Euclidean group E(3), a subgroup of Diff R 3 . That symmetry underlies the rest of this paper.
Consider some φ ∈ Diff R 3 and let χ ∈ Diff R 3 . By the diffeomorphism group's definition χ • φ is also an element of Diff R 3 , so we can define a left-action of Diff R 3 on itself by The corresponding left-action of χ on the spatial potential ϕ is (2.29) To be able to transform action (2.21) using such left-actions we must also consider how derivatives of φ are transformed. Therefore consider some time-dependent embedding (i.e. a motion) φ : I → Diff R 3 . The associated tangent-vector at t = 0 is Given χ ∈ Diff R 3 we can then consider with associated tangent vector (2.32) Thus we see that χ induces a linear mapping from T φ(0) Diff R 3 to T χ·φ(0) Diff R 3 . This is of course true for any t, so for some φ ∈ Diff R 3 and v ∈ T φ Diff R 3 we set which defines the left-action of Diff R 3 on TDiff R 3 . One may also show readily that the deformation gradient transforms to from which one may easily find Diff R 3 's action on the referential Poisson tensor a → χ · a.
The action on φ of Diff R 3 's Lie algebra diff R 3 is also important in what follows. We need only focus on the motion φ because ζ is left-invariant. To proceed, fix φ ∈ Diff R 3 and where χ(0) = id (and with id denoting the identity mapping on a set). We can definė and since χ(0) = id the RHS is by definition Diff R 3 's Lie algebra, diff R 3 . We can now form the curve in Diff R 3 s → χ(s) · φ = χ(s) • φ, (2.37) whose tangent vector at the origin is d ds Given that χ(0) = id, the LHS is an element of T φ Diff R 3 while the RHS represents the action of diff R 3 on Diff R 3 . We have obtained a mapping of Diff R 3 into T φ Diff R 3 by diff R 3 , which we write aṡ This is the desired left-action of diff R 3 .
All in all, action (2.21) transforms under Diff R 3 's left-action into ⟨∇ζ, (χ · a) · ∇ζ⟩ dV dt. (2.40) For general χ ∈ Diff R 3 it is clear that this does not equal the original action; hence the action is not left-invariant under the whole of Diff R 3 . It would be invariant if we could somehow 'remove' the terms in Dχ -but to do that we need to restrict attention to a certain subgroup of Diff R 3 : the Euclidean group E(3).

The Euclidean group and its left-actions
The elements of E(3) ⊆ Diff R 3 are the isometries of 3D space: translations and rotations.

E(3)
is the semi-direct product of T(3) and SO(3), the groups of translations and rotations respectively, which are discussed in Appendix A1. We write an arbitrary element Λ ∈ E(3) as for some Φ ∈ R 3 and R ∈ SO(3); Λ's action on a motion φ ∈ Diff R 3 is (c.f. eq. 2.27) simultaneously translating φ by Φ and rotating it by R. The inverse transformation is Correspondingly, the left-action of the Euclidean group's Lie algebra on the motion is given by for some curve s → Λ(s) with Λ(0) the identity. We show in Appendix A1.2 that the translation group's Lie algebra consists of all velocities V ∈ R 3 , and that τ ∈ t(3) acts like On the other hand, SO(3)'s Lie algebra so(3) contains all possible antisymmetric matrices Ω (Appendix A1.1). In the context of rigid body mechanics these are known as angular velocities, acting on φ like Ω · φ. It follows readily that e(3)'s action is φ → V + Ω · φ. We therefore define an element λ ∈ e(3) to be with its action on a motion described by (2.47) Recall from Appendix A1 (eq. A.8) that the rotation group's adjoint representation Ad allows SO(3) to rotate elements of so(3) into other elements of so (3). Analogously, we may define the adjoint representation of the Euclidean group, Ad : The action of Ad Λ on λ = { V, Ω } produces a new element of e(3) that acts on an arbitrary motion as With these preliminaries we can return to E(3)'s left-actions on S. We see immediately from eq. (2.42) that DΛ • φ = R ∈ SO(3), so by restricting eqs. (2.27), (2.33) and (2.34) from Crucially, looking back to expression (2.40) and recalling the Principle of Material Frame Indifference (eq. 2.9), it follows that the action is invariant under the left-action of the whole Euclidean group.

Noether's theorem
Now that we have established various symmetries of the action, we can use Noether's theorem to study their associated conservation laws. We start by deriving the general result in the appropriate context, then we apply it to obtain some familiar results. Later (Section 3.5) we show that the same ideas can be applied within fluid regions to obtain an apparently novel conservation law.
Consider an action taking the general form where on the LHS we make explicit the action's functional dependence on the motion. Its variation is and, given that (2.54) If the variation δφ vanishes at the start-and endpoints, then we may appeal to Hamilton's principle and demand that δS = 0, leading to the weak-form Euler-Lagrange equations for time-independent test functions δφ and with δF = D (δφ).

Now consider the one-parameter family of fields
where φ 0 is a particular solution to eq. (2.55); we describe φ 0 as 'on-shell'. If the condition is satisfied for all s, then the family φ s is a one-parameter continuous symmetry of the action.
Differentiating this with respect to s about s = 0 gives which is sometimes more practically useful. This derivative is given explicitly by eq. (2.54), but with (2.59) This δφ need not satisfy the vanishing-endpoint condition, but φ 0 itself is on-shell, so eq. (2.58) Given that the start-and end-times are arbitrary it follows that the quantity is conserved in time. This is Noether's theorem. The result shows concretely how the conserved quantity corresponds to the symmetry generating it, and gives a means of computing it.

(Angular-) momentum conservation
For the specific action (2.21) we now consider a definite solution φ 0 and ζ 0 . Then we define the curve s → Λ(s) ∈ E(3), and thus define the one-parameter family with Λ(0) equal to the identity. Given that the action is invariant under the left-action of E(3), φ s and ζ s satisfy and are therefore a one-parameter symmetry of the action (c.f. eq. 2.57). Since Noether's theorem (eq. 2.61) shows that λ is arbitrary, so we have shown the existence of a conserved quantity π ∈ e(3), which we define implicitly through for arbitrary λ ∈ e(3) and where the e(3)-inner-product on the LHS is defined in the obvious manner. We refer to π as the generalised momentum. From eq. (2.47) and we may use the wedge product defined by eq. (A.13) to rearrange this as Since V and Ω are arbitrary, B ρv dV and B ρv ∧ φ dV are both conserved. Conservation of the generalised momentum π is thus seen to encompass conservation of both the linear and angular momenta of the elastic body. For later reference we define We remark that this exercise has also demonstrated that the wedge product is analogous to the standard cross product.

Energy conservation
If we temporarily neglect the strain-energy's stress-glut term, then action (2.21) is also invariant when its fields are shifted by a time τ and a corresponding shift is made in where we are temporarily viewing the action as a function of the time-interval as well. The one-parameter family is therefore another one-parameter continuous symmetry of the action. It follows that with the term in L arising from the variation of I. There is no term here in ζ τ because the action does not depend on ζ's time-derivatives. Concretely, Noether's theorem thus implies that the total energy is conserved.

Decomposition of the motion
Motivated by the action's invariance under a constant Euclidean motion, we now introduce the Tisserand frame by writing φ as an elastic motion upon which a time-dependent Euclidean motion is imposed: is the Euclidean motion, which captures the motion of the frame, and where φ : I → Diff R 3 is what we term the internal motion. Now, Λ(t)'s left-action will alter the action via T , so this might seem to complicate the problem. However, introducing the new function Λ(t) frees us -indeed, obliges us -to impose constraints upon φ, and we will choose those constraints so that φ has no secular behaviour and couples minimally to the Euclidean motion. The remainder of this section lays all this out. In what follows we will generally suppress the respective motions' time-arguments.

The decomposition's kinematics
Given eq. (2.77), the corresponding velocity is To write the velocity more conveniently we define both of which take values in e(3). λ and λ are related by Λ's adjoint representation (eq. 2.48): It now follows that v may be written in two distinct ways: These expressions show that one may choose to express the velocity such that it features either the internal motion φ or the original motion φ.
Similarly to the velocity, we may write the conserved generalised momentum π (eq. 2.66) in two different ways. Using eqs. (2.77, 2.81a), where κ denotes an arbitrary element of e(3). We now define the spatial and internal effective mass tensors J, J to be with κ ′ ∈ e(3) also arbitrary. They map e(3) onto itself; they are clearly symmetric and positive-definite, and therefore invertible. If we then set the second term on eq. (2.82)'s RHS may be written either as (J · λ, κ) or J · λ, κ . Next, we After defining the internal momentum π by the momentum π may finally be written as: Eq. (2.87a) is the more useful for Section 2.4.2 while (2.87b) will be used almost exclusively from Section 2.5 onwards.

The decomposition's uniqueness
As written, decomposition (2.77) is non-unique. That is, given a definite φ one cannot uniquely construct Λ and φ. This is for two similar but distinct reasons. Firstly, we have introduced six extra degrees of freedom at each time t through the new function Λ. Secondly, the decomposition is only defined up to an arbitrary Euclidean motion: if eq. (2.77) holds for some Λ and φ, then it also holds for ΛΛ −1 andΛ · φ withΛ ∈ E(3) arbitrary. However, requiring the internal motion's velocity not to contribute to the internal momentum is sufficient to make decomposition (2.77) unique up to initial conditions on Λ, and we can in fact choose those to be whatever is convenient.
To begin, consider eq. (2.87a)'s statement that the conserved momentum π is equal to the internal momentum π plus a term related to λ and φ. π represents the linear and angular momentum associated with the internal elastic motion φ, so in order to satisfy the Tisserand condition it seems desirable to demand that π vanish. We therefore set to be a constraint associated with decomposition (2.77).
To prove that this constraint is sound, we start by noting that it reduces eq. (2.87a) to the statement π = J · λ. (2.89) Using J's invertibility and λ's definition (2.79a), this implies thaṫ Given a definite motion φ we can directly construct both π and J, then solve this ODE to give Λ. The final step is to write inverting eq. (2.77) to find φ. The solution to eq. (2.90) is unique only up to an arbitrary initial condition Λ 0 , and therefore so is expression (2.91) for φ. However, an arbitrary change of initial condition Λ 0 → Λ 0Λ merely transforms Λ and φ by Such a transformation reflects the decomposition's inherent non-uniqueness, mentioned just above. Crucially, it has no effect on the physically observable motion φ. The arbitrary initial condition Λ 0 is thus a further example of a gauge freedom. Therefore we may reasonably describe eqs. (2.77) and (2.88) as together constituting a unique decomposition of the motion.

A pragmatic notation and choice of initial conditions
It will now be convenient to reintroduce an explicit separation between the translational and rotational components of the Euclidean motion, writing with Φ : I → R 3 and R : I → SO(3). The motion is then written as To write the velocity we need the actions of λ and λ. This leads us to define the spatial angular velocity and body angular velocity by both take values in so (3). Then λ's action is (eq. 2.79a) Substituting these explicit actions of e(3) into eqs.
(2.81) shows that the velocity is written as Next, to consider the generalised momentum we should start by introducing the spatial and internal moment of inertia (MoI) tensors: for arbitrary Ω ′ ∈ so (3). We also note that the constraint π = 0 is equivalent to conservation of internal linear and angular momentum: It then follows upon substitution of eqs. (2.94, 2.98) into eqs. (2.69) that the conserved spatial linear and angular momenta take the forms These expressions, which are equivalent to eq. (2.89), clarify that we still have not achieved a 'minimal coupling' of the internal and Euclidean motions inasmuch as P and J still depend on φ. Fortunately though, they depend on φ only through B ρφ dV which, courtesy of constraint (2.100a), is constant in time and therefore only represents three degrees of freedom.
We therefore declare that B ρφ dV = 0, (2.102) and will now show that this is allowed by our freedom to choose initial conditions on Φ and R.
Given definite initial conditions on the spatial motion φ(x, 0) = φ 0 (the initial condition on the velocity is irrelevant at present) we trivially have where Φ 0 and R 0 are (so far) arbitrary initial conditions on Φ and R. Integrating over B gives It follows that a necessary and sufficient condition for eq. (2.102) to hold is i.e. that Φ 0 gives the initial position of the body's CoM; R 0 remains arbitrary. As long as we use these initial conditions, then we may replace constraint (2.100a) with the stronger constraint (2.102) henceforward.
This new constraint means that Φ is the CoM of the body at all times and carries all the body's linear momentum: eq. (2.101a) becomes which we may integrate trivially to give The total angular momentum also factors neatly into the orbital angular momentum of the CoM's linear motion, and spin angular momentum associated with rotation about the CoM: (2.108) With the moment of inertia and angular velocity as defined, the form of the spin angular momentum is even the same as for a rigid body, but with an 'effective' moment of inertia that accounts for the body's deformation.

Decomposed action
We now substitute eq. (2.77) into the action in order to obtain a variational principle that depends on both the Euclidean and internal motions. From now on we will use 'internal' expressions almost exclusively. Using eqs. (2.79b, 2.81b) to write The original constraint π = 0 knocks out the second term, while the choice of Φ 0 eliminates the cross-terms between v and Ω in the kinetic energy: (2.112) Thus we may write the action as ⟨∇ζ, a · ∇ζ⟩ dV dt, and with the motion given by eq. (2.94). Intuitively, we can see eqs.
(2.94) and (2.114) as describing elastic motion that takes place in a rotating frame whose orientation is given by R relative to an origin at Φ, with the elastic motion within the rotating frame carrying neither linear nor angular momentum and with its CoM always at the origin.
We could of course have chosen the constraints differently, but the choice we have made seems to be the most convenient. Indeed, the action has factored into a 'fully internal' part (the third and fourth terms) cosmetically identical to eq. (2.21)'s action, but now augmented by two terms which feature the Euclidean motion. Those first two terms are identical in form to a rigid body's action (e.g. Holm et al. 2009), but with the one subtle difference that I here is not constant due to its dependence on φ. Most importantly, the constraints ensure that φ can have no secular behaviour.

The reconstruction equation
The decomposed action has no explicit dependence on the Euclidean motion. Λ only enters implicitly via λ ≡ Λ −1Λ . The variational principle will provide differential equations for λ, φ and ζ, but not for Λ. In order to write down the full, spatial motion φ = Λ · φ, we need to be able to relate Λ to λ, φ and ζ. That relationship is the reconstruction equatioṅ an ODE on E(3) obtained by rearranging eq. (2.79b). To couch eq. (2.115) in terms of Φ, R etc., we use the identities for arbitrary a ∈ R 3 . It follows that eq. (2.115) implies the trivial ODEΦ = V, as well aṡ Henceforward we refer to eq. (2.117) without ambiguity as "the reconstruction equation"and we emphasise that R 0 may still be chosen according to convenience.
The reconstruction equation is a kinematic identity that may be stated without reference to the internal dynamics. And as alluded to just above, we see concretely that R only enters via an ODE once the rest of the equations have been solved and Ω is known fully. We do not have to integrate the reconstruction equation simultaneously with the others. Nevertheless, eq. (2.117) acts as a further constraint on the action that must be accounted for when deriving equations of motion.

Constraining the action
We will impose the reconstruction equation strongly, that is without the use of Lagrange multipliers. On the other hand we choose to impose the Euclidean constraints weakly, introducing the Lagrange multipliers The action now takes the final form ⟨∇ζ, a · ∇ζ⟩ dV dt. (2.119) The action is to be viewed as a functional of V, Ω, φ and ζ, as well as A and Ω.
We have chosen to impose the constraints in this mixed weak-strong way for the following reasons.
Regarding Ω, it is easy to impose its constraint (2.117) directly upon the variation δΩ at the point of applying Hamilton's principle, so there is no need for a Lagrange multiplier (see below, Section 2.6.1). The analogous procedure does not work (not to our knowledge, at least) for the angular-momentum constraint, so that must be applied through a Lagrange multiplier. If only for aesthetic reasons we impose the translational constraint the same way.

Equations of motion
We are now ready to derive the differential equations governing the Euclidean and internal motions. As a final piece of pragmatism we will drop all the overlines on internal variables henceforward, reinstating them only on the rare occasions that they are needed to avoid ambiguity.

Variation of the action
In order to obtain the EoMs we vary with respect to all variables and demand that the first variation of the action vanish. Starting with the Lagrange multipliers, variation in A and Ω trivially returns the Euclidean constraints, Varying with respect to Φ is almost as easy: we find that As expected for a system with no external forcing, the CoM moves at constant velocity.
Variation with respect to Ω is more involved because the variation δΩ ∈ so(3) cannot be chosen freely due to the constraint Ω = R TṘ . Instead we must consider arbitrary variations in R. We may write a small change in R as R → Re δθ (2.125) for arbitrary small δθ ∈ so(3); this way the group structure is respected and the right hand side is guaranteed to be an element of SO(3). The exponential is then Taylor-expanded to give For the variation of Ω we therefore have where the little-adjoint ad is defined in eq. (A.10). We can schematically write any term in the action involving Ω as I f (Ω) dt for some scalar-valued function f : so(3) → I; recalling that all fields are varied subject to fixed endpoint conditions, it then follows that any term in the action involving Ω satisfies so that the corresponding variation is For the present case f (Ω) = 1 2 ⟨Ω, I · Ω⟩, so this implies We have derived the equation which is the celebrated Euler equation of rigid body dynamics, but with a body MoI that depends on φ and therefore on time.
Next, variation with respect to the referential potential is accomplished readily: ⟨a · ∇ζ, ∇δζ⟩ dV dt = 0. (2.132) Noting that this expression contains no time-derivatives of the variation δζ, it is equivalent to B ρ δζ dV + 1 4πG R 3 ⟨a · ∇ζ, δζ⟩ dV = 0, (2.133) where δζ is a time-independent test-function. This is of course the referential Poisson equation in weak-form.
Now for variation in φ. From definition (2.99b), the MoI's variation satisfies (2.138) The analogous first Piola-Kirchoff gravitational stress is ⟨∇ζ, a · ∇ζ⟩ . (2.139) If we define the referential gravitational acceleration 140) then N is given explicitly by with 1 the identity tensor. It is convenient to split the integral over R 3 into respective integrals over B and R 3 /B: ⟨N, δF⟩ dV dt = 0. (2.142) The exterior term may be integrated by parts to give where δφ now denotes a time-independent test-function.

Identification of the Lagrange multipliers
To complete the derivation we eliminate the Lagrange multipliers A and Ω. For Ω take δφ = θ · φ so that δF = θF for arbitrary θ ∈ so (3) and T is such that TF T is symmetric; an analogous result holds for the gravitational stress.
The kinetic term vanishes because of the rotational constraint (and θ's antisymmetry): and, after a little algebra, one may show that We conclude that the Lagrange multiplier satisfies Given eq. (2.131), one may obtain a solution to the present equation by setting Ω = Ω. (2.151)

Summary
We can summarise the results of the preceding derivation as follows: We have renamed δφ to w and δζ to χ in order to help clarify that we now regard them as arbitrary test-functions. The first three equations describe the evolution of the Tisserand frame. the corresponding strong-form equations are ever needed they could be derived readily, but weak-forms are generally more useful for numerical work and are therefore our focus here and elsewhere in the paper.

EXTENSION TO FLUID-SOLID PLANETS
If we wish to model Earth-like planets, then we must adapt the preceding approach to account for (visco)elastic bodies composed of alternating fluid and solid layers. That is the subject of Sections 3.2-3.4. Section 3.5 discusses a symmetry group associated with fluid regions along with the associated conservation law.

The form of the fluid strain-energy function
Within an elastic fluid, the strain-energy function takes the form (e.g. Noll 1974) for some function V and where we recall that J = det F. It follows that the FPK stress is where we have defined the pressure and made use of Jacobi's formula (e.g. Holzapfel 2000) for the derivative of a determinant. Figure 1. A three-layer fluid-solid planet (after D&T, Fig. 3.1). An arbitrary point x in the reference body B is mapped to φ(x, t) in the instantaneous configuration B t . B S and B F respectively denote the solid and fluid regions of the reference body. Σ F S refers to the union of all internal fluid-solid boundaries, while ∂B is the free-surface. A unit vectorn on a boundary is mapped by φ ton t ∝ JF −T ·n. In this illustration the uppermost layer is solid, but this Section's theoretical developments do not require that.

Tangential slip
Consider a body composed of N distinct layers, alternating between fluid and solid; Fig. 1 shows a special case consisting of three layers. It is not relevant at present which of the layers are fluid, only that the layers possess mutual fluid-solid boundaries where the fluid is able to slide freely past the solid. Such tangential slip is the crucial kinematical difference between solid and fluid-solid bodies. Otherwise all of our earlier kinematical considerations carry over to a layered fluid-solid planet.

Illustrative special case: two layers
We assume that the reference body B contains a single referential fluid-solid boundary Σ F S .
Consider two particles on either side of Σ F S that are instantaneously adjacent at some time.
During the subsequent motion of the body these particles can separate from one another as the fluid slides tangentially past the moving solid interface. The motion of such a body can therefore be discontinuous across its fluid-solid boundaries, but as long as the layers neither separate nor interpenetrate, the motion is subject to the constraint that where φ − (resp. φ + ) gives the motion evaluated immediately below (resp. above) Σ F S . This relation does not require the point-wise equality φ − (x, t) = φ + (x, t), but only the weaker condition that for each x ∈ Σ F S there is at time t a unique point for some χ ∈ Diff (Σ F S ). Due to this constraint, the configuration space for the body is no longer the diffeomorphism group, Diff R 3 , but must be extended to allow for suitable piecewise-diffeomorphisms.

Extension to N layers
For an N -layer planet the reference body B is decomposed in the form with referential fluid-solid boundaries (3.7) Following D&T, we now use the symbol Σ F S to denote the union of those boundaries: The same kinematic considerations as above hold for each of the fluid-solid boundaries individually, so the motion's restriction to Σ F S must satisfy where once again χ ∈ Diff (Σ F S ) and we use a subscript plus and minus to denote evaluation of φ just above and below Σ F S . This notation allows us to deal with any number of fluidsolid boundaries without any cosmetic difference from the case of a single such boundary. In particular, it avoids the need to explicitly decompose B into multiple sub-bodies. Nevertheless, in Sections 3.5 and 3.7 it will be helpful to consider B's fluid and solid regions separately. To that end we import from D&T the notation B F (resp. B S ) to denote the union of all of B's fluid (resp. solid) regions, with B = B F ∪ B S .

Variational principle (or, 'Everything else is basically the same')
There is no potential energy associated specifically with such tangential slip at a fluid-solid boundary. Although the tangential slip modifies the problem's configuration space, the action of a fluid-solid body is therefore unchanged from eq. (2.21). Indeed, the tangential slip constraint has no impact on the derivations of Sections 2.2, 2.3 and 2.4. The action is of course subject to constraint (3.9), but it still exhibits the same symmetries under left-and right-actions of Diff R 3 and E(3); the body's total momentum and energy are conserved; and the decomposition may clearly still be made without violating eq. (3.9). In fact, the tangential slip constraint on φ takes the same form when applied to φ: for some χ ∈ Diff (Σ F S ). We now revert to notating the internal motion by φ and the body angular velocity by Ω.
It is crucial to recognise that both the motion and the test-functions are constrained, and moreover that the test-functions satisfy a different constraint from the motion. Not only do we have φ − • χ = φ + (3.12) on Σ F S , but also, taking this equation's first-variation, Since χ ∈ Diff (Σ F S ), δχ is tangent to the fluid-solid boundary at χ(x) (see A18, eq. 55), so this is equivalent to withn − (x) the outward-pointing normal to Σ F S at x (that is, the outward-pointing normal with respect to the sub-body just beneath the boundary). This translates straight into the constraint satisfied by the test-functions: All in all, we may take Section 2.5's action (2.119) as the starting point for our analysis of the fluid-solid body, always remembering that it is subject to the strongly-imposed reconstruction equation (2.117) and tangential slip constraints (3.12) and (3.15), as well as to the weaklyimposed Euclidean constraints.

Variation of the action
To derive the fluid-solid Euler-Lagrange equations we may at first vary the action exactly as we did in Section 2.6. The necessity to include the tangential slip constraint on the motion poses no particular problem; given that we will impose it directly on the equations of motion, the algebra required to vary the action is the same as before. Moreover, the forms that we chose for the test-functions in Section 2.6.2 trivially respect the tangential slip constraint, so it remains the case that the Lagrange multipliers are Eqs. (2.152) are therefore still valid in the present case. We must now impose the tangential slip constraint suitably on the test functions.

Imposing tangential slip
On a practical level, we will impose eq. (3.15) weakly even as we continue to impose eq. (3.12) directly on the motion. To do that we define new test-functions ϖ : Σ F S × I → R and add to the equations of motion the term Henceforward, an unsubscripted variable in the vicinity of Σ F S that could be evaluated on either side of the boundary should be understood to take the value on the lower side; for example, F| Σ F S = F − | Σ F S . Also, we will loosely refer to ϖ as a 'Lagrange multiplier' because, just like A and Ω, it is there to enforce a constraint weakly.
We now augment the equations of motion to incorporate the tangential slip constraints on both the motion and test-functions. Nothing changes except in the momentum equation: where we recall that χ is defined through φ − • χ = φ + .

Identification of the Lagrange multiplier
To find the multiplier enforcing tangential slip on the test-functions we integrate eq. (3.19)'s term in Dw by parts. By requiring that the strong-form of the momentum equation hold within each sub-region, and by using the traction-free boundary conditions on the surface, we are then left with (3.20) Now,n + = −n − becausen + points outward with respect to the body on the upper side of the boundary. Therefore if we 'pair up' the various instances of w − and w + we find Because w − and w + are independent, we obtain from the first term where we have defined the traction vectors t − = T − ·n and n − = N − ·n for convenience.
Similarly, from the second part of the equality we find where we have used χ to change variables in part of the surface integral, with J χ being the associated Jacobian. Comparing these equalities, we firstly have This states the equality of mechanical and gravitational tractions at points that are instantaneously adjactent on the fluid-solid boundary (note that the surface Jacobian term is due to the tractions being measured per unit area on different parts of the reference boundary).

From Poisson's equation, it can be show that continuity of the gravitational traction follows
automatically (e.g. Dahlen & Tromp 1998, Section 2.8.3), and hence this equality reduces to

Summary
The full set of equations of motion for a single, variably-rotating, fluid-solid planet are noẇ where we recall that χ : After having solved for the Lagrange multipliers these equations are invariant under any transformation of the form for arbitrary a and θ (c.f. Section 2.6.2). We are therefore free to constrain the test-functions to satisfy Such constraints will be useful within certain applications. In particular, they will lead in Section 3.8 to the test-functions satisfying identical constraints to the linearised motion; this being necessary within Galerkin methods.
To summarise, we have obtained equations of motion that describe the deformation of a rotating, self-gravitating, fluid-solid elastic body in terms of: (i) the motion of its centre of mass, Φ; (ii) a rotation about that centre of mass, R; (iii) and elastic internal motion φ with respect to the frame defined by Φ and R.
Comparing these results to AC18, who discussed variational principles for fluid-solid bodies relative to a steadily rotating reference frame, we see that no additional complications arise from the use of the Tisserand frame.

The fluid relabelling group
Recall from Section 2.2.1 that the action of a solid elastic body is form-invariant under particle-relabelling transformations. Such transformations do not, in general, constitute a full symmetry of the action because they change the material parameters (ρ, W ). It is interesting to ask whether there exist non-trivial transformations, ξ ∈ Diff R 3 , such that these parameters are invariant, this meaning where we recall eq. (2.26a) and (2.26b) for the definition of these group actions.
To examine this question, we start with invariance of density. For an arbitrary subset where we have used ξ to change variables in obtaining the right hand side of the equality.
Supposing that ξ leaves ρ invariant, we then obtain the relation While such an equality cannot hold in general, there are situations where it does. For example, when ρ is a constant and ξ is volume preserving. More generally, it is not difficult to see that the stated equality holds if there exists a reference configuration of the body with respect to which the density it constant. In physical terms, this situation describes a body that has a uniform composition, meaning that spatial variations in density within the chosen reference configuration are due solely to deformation from this hypothetical undeformed state. In fact, the existence of such a reference configuration follows from a theorem of Moser (1965). These arguments extend to piecewise continuous densities so long as the support of ξ is suitably limited. We conclude that there do exist non-trivial particle relabelling transformations that leave the density invariant.
Turning to the invariance of the strain energy function we can apply a similar method.
Again for arbitrary U ⊆ B, we can use ξ to change variables to arrive at where F is here an arbitrary field taking values in GL(3). If W is invariant under ξ, then this equality reduces to (3.33) and because F is arbitrary we are free to choose F = 1. In this manner we have arrived at the (3.34) that must hold if the action of ξ leaves the strain energy function invariant. This identity has a simple physical interpretation. On the left hand side we have the elastic potential energy of the particles contained within U within the reference configuration, while on the right is the energy of these particles following their deformation under ξ. It follows that a particle relabelling transformation leaves W invariant only if this transformation, regarded as a configuration of the body, would be associated with no change in the elastic potential energy. For an elastic solid, this is only possible for elements of the Euclidean group, but these can be discounted because we also require that B be mapped onto itself (a rigid rotation of a spherically symmetric reference body about its centre would be a counter example, but this is not interesting in practice). It follows that for solid elastic bodies there are no non-trivial particle relabelling transformations that leave the strain energy invariant.
The preceding argument can be specialised to a body containing a fluid sub-region, B F , with more interesting results. Supposing that U ⊆ B F , and recalling eq. (3.1), we find that the above necessary identity reduces to The interpretation of this equation is as before: that the elastic potential energy of the particles within U is unchanged under a deformation by ξ. Here, however, this could hold for non-trivial mappings. Indeed, the characteristic property of an elastic fluid is that its potential energy is determined by changes to volume and not to shape (e.g. Noll 1974). To go further, we assume without essential loss of generality, that the reference configuration is stress free. This means that J = 1 is a minimum of the V , and we are free to adjust V such that this minimum value is zero. If, further, we suppose that V is a non-negative and strictly convex function of its second argument then the equality requires J ξ = 1. Note that the latter assumptions are met for conventional forms of V such as a fluid Saint-Venant-Kirchhoff material (e.g. Holzapfel

2000) for which
V (x, J) = 1 2 κ(x)(ln J) 2 , (3.36) with κ a positive definite constant (equal to the bulk modulus within the corresponding linearised theory). We conclude that non-trivial particle relabelling transformations that leave V invariant exist so long as they are volume preserving and satisfy Such conditions can, in general, be met.
Summarising the above arguments, we have shown that the action is invariant under a particle relabelling transformation if its support is restricted to B F , and if the following conditions are met where in the final equation x ∈ B F and J is arbitrary. These conditions define a sub-group of the particle relabelling transformations that we term the fluid relabelling group FRL(B).
Such transformations can be said to map level surfaces of both ρ and V into themselves, and hence following such a transformation it would appear to an observer shown only the initial and final states that nothing had happened. This group is equivalent to the geostrophic modes discussed by D&T (p. 118) within the context of linearised elasticity.
To proceed, we characterise the Lie algebra of FRL(B). Varying ξ → id + δξ, we find quickly that Thus, the Lie algebra comprises divergence-free vector fields supported in B F that are everywhere tangent to the level surfaces of density and the strain energy. In fact, it would be usual for these level surface to be coincident (e.g. Stevenson 1987), in which case we could drop one of the final two conditions. As a special case, there are neutrally stratified fluids for which FRL(B) is equal to all volume preserving diffeomorphisms with support contained in B.
Noether's theorem yields the conservation law arising from invariance of the action under where s → ξ s is a curve in FRL(B) with ξ 0 = id. We know from earlier that B ⟨D v L, δφ⟩ dV is conserved, with (3.43) For the chosen form of φ s , where we define (3.45) Finally, since D v L = ρ (v + Ω · φ), the conserved quantity is for arbitrary δξ within the Lie algebra of FRL(B). The derivation of such a conservation law from a particle relabelling symmetry resembles work within the context of fluid mechanics (e.g. Müller 1995;Albert 1997;Névir 2004;Bridges et al. 2005), but a precise correspondence has not been established.
The full implications of this conservation law will be discussed elsewhere. But we conclude this discussion with the case of a neutrally stratified fluid. The conservation law now states that the restriction of ρF T · (v + Ω · φ) to B F is orthogonal to the space of divergence-free vector fields on this domain with vanishing normal components. From the Helmholtz decomposition theorem (e.g. Abraham et al. 2012), it follows that for some scalar field, χ say, we have

Accounting for the oceans and other surface loads
We now consider how the foregoing equations can be modified to incorporate a thin surface layer such as an ocean and/or an icesheet. Let us regard the layer as an 'extra' body B imposed on top of the original reference body B. In effect, we rename B → B ∪ B. We then integrate the momentum-equation in B by parts, and impose that the strong-form of the momentum-equation holds within B. This gives the following modified momentum-equation where t is the traction applied by the surface layer. This is an exact result that could be used in modelling, say, the coupled evolution of the solid Earth and an ice-sheet. However, in many situations it is sufficient to approximate t so as to avoid explicit consideration of the surface layer's dynamics. For example, Komatitsch & Tromp (2002) use the water column approximation to derive a form of t that is useful for describing the ocean's effect on seismic waves within the solid Earth at sufficiently long periods. Another example is discussed below in the context of quasi-static problems. Note that we must also take account of B in the Poisson equation, MoI and total mass, but that simple thin-layer approximations are sufficient for that purpose.

Equations of motion
Relative to the Tisserand frame, a relative equilibrium for a planet is characterised by V, Ω, φ, and ζ being independent of time. Such a planet can still be moving at a constant velocity through inertial space and rotating steadily, but there is no internal deformation. In this situation, the weak form of the momentum equation reduces to while the Euler equation becomes ad Ω (I · Ω) = 0. (3.50) As ever, it is assumed that the tangential slip constraint on the motion has been imposed strongly. We also have the translational constraint in eq. (3.26f) on the internal motion, but note that the angular momentum constraint in eq. (3.26g) holds trivially at equilibrium.
Within this section, we consider quasi-static deformation associated with prescribed variations in a thin surface layer, this situation being, for example, relevant to the study of GIA. As discussed above, the presence of such a surface layer can be incorporated into the momentum equation by prescribing a suitable surface traction, t, while corresponding modifications are made to other equations (detailed below). Within the surface layer, the strong form of the momentum equation reads or equivalently as where we have written γ = −F −T ∇ζ for the referential gravitational acceleration. To first-order accuracy in the layer thickness, we can integrate this equation to obtain the expression for the applied traction, where σ denotes a mass per unit area that we will call the load. Using this result, we can modify the quasi-static momentum equation to read The load also enters into the weak form of the Poisson equation as into the definition of the MoI tensor (defined relative to the Tisserand frame) In this manner, we have obtained a complete set of equations for modelling quasi-static deformation due to prescribed variations within a surface layer subject to a thin-sheet approximation.
Within these equations, deformation is driven by the load, σ, which is a given function on ∂B × I.

Invariance under the Euclidean group
Suppose that (Ω, φ, ζ) solve the quasi-static equations of motion subject to a given load, σ.
For arbitrary Q ∈ SO(3), it is readily seen that provides another solution. Indeed, the only non-trivial step is checking that the MoI tensor transforms as which follows directly from its definition. Here we have arrived at the familiar fact that solutions of quasi-static problems are defined only up to an arbitrary rotation (e.g. Marsden & Hughes 1994). The fact that Φ drops out of the equilibrium equations shows that the solution is also defined only up to an arbitrary translation.
There are two ways to address this non-uniqueness. First, one can restrict attention to observations that are invariant under E(3) such as relative distances. The non-uniqueness is then simply a gauge-freedom within the problem, and can be fixed in whatever manner is convenient. The alternative is to remember that quasi-static problems arise as an approximation to the full dynamics in cases where the forcing is slow relative to the time-scales of internal deformation. We can, therefore, appeal to conservation laws from the full problem to fix the undetermined Euclidean transformation. To see how this is done suppose for definiteness that initially the planet was at rest relative to inertial space (an inertial frame can always be chosen for which this is true because V is constant). It follows that the net linear momentum is always zero, and hence there can be no translational part to the quasi-static solution. Regarding rotation, it follows from the quasi-static Euler equation that Ω is always parallel to the spatial angular momentum. If, in solving the quasi-static problem, this condition has not been met, then we simply rotate Ω along with the other variables so that it holds. The necessary rotation is determined uniquely up to a rotation parallel to the angular momentum whose value is determined, if necessary, by appeal to the initial conditions on the rotation matrix, R, that described the Tisserand frame.

Invariance under the fluid relabelling group
Quasi-static calculations in planets possessing fluid sub-regions present additional challenges (e.g. Longman 1963;Dahlen 1974;Wunsch 1974;Crossley & Gubbins 1975). A satisfactory resolution was provided by Dahlen (1974), who observed that linearisation of the displacement was not valid within fluid regions. Instead, he presented a mixed formulation of the problem, with solid regions described referentially, and fluid regions spatially. In particular, it was shown that all relevant dynamic variables within fluid regions could be described in terms of the perturbed gravitational potential. Such an approach provides a practical method for quasistatic calculations that has been applied widely (e.g. Tromp & Mitrovica 1999;Al-Attar & Tromp 2014;Crawford et al. 2018). It is not clear whether Dahlen's approach can be extended to non-linear quasi-static calculations. Here, instead, we outline an alternative method that should work for both linearised and non-linear calculations.
Recall from Section 3.5 that the action of an elastic body is invariant under the action of the fluid relabelling group, FRL(B). Within the context of quasi-static calculations, this implies that solutions can only be defined up to an element of this infinite-dimensional group.
Indeed, this is, precisely the reason that linearisation of the displacement in fluid regions cannot be justified, as Dahlen described in different language. Crucially however, we know from our earlier discussion that the action of FRL(B) leaves φ and ζ unchanged outside B F , and has no effect on Ω. From an observational perspective, therefore, this is a gauge symmetry, with any solution of the problem being as good as another. There remains, however, the practical question as to how any solution can be obtained numerically in a well-posed manner. To proceed, we can seek the solution of the quasi-static equations that minimises the following where µ is a positive-definite scalar field. Note that this functional is, by design, invariant under SO (3), and so our earlier discussion of that group remains valid. In essence, this choice amounts to picking from amongst all possible solutions the one that is closest to a Euclidean transformation within fluid regions. Other functionals could be chosen, with all the matters being that they lead to a well-posed problem. To implement this idea practically, a Lagrange multiplier method could be applied (c.f. Seliger & Whitham 1968;Al-Attar & Woodhouse 2010). A full discussion of this idea will be presented elsewhere, with its practical viability to be established.

Linearisation
If a planet at equilibrium is disturbed by some small force, then we can approximate its motion by linearising about that equilibrium state. For definiteness, we imagine that a small stress-glut sets the planet in motion, so we introduce a dimensionless perturbation parameter ϵ and write To linearise about the equilibrium described in Section 3.7.1 we substitute the perturbative into the EoMs (3.26). In principle we could perturb the test-functions too. After all, the constraints cause w and ϖ to intertwine closely with φ. However, it is easy to show that the O (ϵ) part of the test-functions would simply repeat the zeroth-order equations of motion.
because it decouples fully from the other variables.

Derivation of the linearised equations
First we linearise the Euclidean constraints. Given that we take the zeroth-order motion to satisfy the equilibrium equations we find readily that We also recall that with R (0) having been set at equilibrium, the initial condition θ (1) (0) constitutes a gauge-freedom of the first-order problem; we set B ρu (1) ∧ φ (0) dV = 0 (3.65) by choosing θ (1) (0) suitably.
Next, the reconstruction equation reduces to while the Euler equation becomes upon using the equilibrium solution I (0) · Ω (0) = λΩ (0) . Linearising the momentum and Poisson equations involves tedious algebra that we have consigned to Appendix B. We find the Poisson and the momentum equation Dw dV,(3.69) with the elastic tensor defined as and where Q is defined in eq. (B.5), G in eq. (B.11) and Γ in eq. (B.14). The linearised motion u (1) and test-functions w are elements of the function-space F φ (0) defined by eq. (B.7).
Although the boundary terms have been formulated differently, these equations are equivalent to A18's eq. (120) aside from our use of the Tisserand frame. The latter equations are equivalent in turn to those of Woodhouse & Dahlen (1978) when the equilibrium configuration is the identity mapping.

A streamlined notation
As in A18, it will be helpful to establish a more streamlined notation before proceeding. For a start, we drop superscript ones and zeroes, and rename φ (0) → Ψ. Then we define the bilinear (3.73) For the linearised rotational parts of the momentum equation we may also usefully define For the Euler equation, on the other hand, we need (3.76) Next we define the bilinear form for the Poisson term, while the terms that couple the motion to the gravity require With this we can rewrite the linearised equations more conveniently as 3.9 An application: generalised normal-mode coupling for a rotating N -layer fluid-solid planet To derive the mode-coupling equations we now expand the first-order fields and test-functions on a suitable set of basis functions. For NMC calculations those basis functions are essentially the eigenfunctions of a reference, spherically-symmetric Earth model, but altered so that they satisfy the strong-form constraints imposed on the problem. In this paper we will not dwell on the details of imposing the constraints numerically, deferring such considerations to future work on applications. For now we just expand on some set of basis-functions {û n } and {ζ p } , with the resulting equations assumed to incorporate the necessary constraints. The discussion is applicable to any Galerkin method.
We express the linearised motion and test-functions as for the gravity we set then we substitute these expansions into eqs. (3.79). Next we take each of the w n and c p to be unity in turn -with all of the others vanishing -thus projecting out the infinite-dimensional linear-algebraic ODEs to give the following frequency-domain equations of GNMC: These equations, together with their time-domain counterparts (3.84), are one of this paper's main results.

INTERACTIONS BETWEEN MULTIPLE BODIES
In some problems one will need to consider multiple interacting bodies. Such problems can be approached using the formalism we have just presented, as we now show by discussing two problems: one of orbital dynamics (Section 4.1), the other related to Earth's rotational variations (Section 4.2). This section is by no means intended to be exhaustive. The examples we give are for methodological illustration -and even speculation -so we mainly present actions and comment on physics with minimal mathematical formality. We plan to revisit both of the following problems in future application-based work, at which time we will develop these ideas more fully.

Application: generalised orbital dynamics
Let two fluid-solid planets interact gravitationally. The action describing such a system is the sum of the planets' respective 'single-body actions' (eq. 2.119) but with an extra term added to account for the planets' gravitational binding energy. That binding energy is just the spatial gravitational potential of one body integrated against the other's density, so if we define and recall definition (2.120) of the internal elastic action R, then the total action is We have added subscripts in an obvious way. Each φ i of course satisfies the tangential-slip constraint (3.11) while Ω i is related to R i by the reconstruction equation (2.117).
Now let us make a change of variables learned from the Keplerian two-body problem. We repackage m 1 and m 2 into a total mass M and a reduced mass µ, then we define the bodies' relative position r and the barycentre position d by (4.4b) These transformations hit the top line of eq. (4.2) to give We see immediately that d is a cyclic variable, i.e. the two-body system's overall CoM moves at constant velocity. r on the other hand satisfies as follows quickly by varying action (4.5). Since r will generally have much greater magnitude than φ 1,2 this can clearly be approximated as (4.7) The planets' relative translational motion thus satisfies a modified Keplerian two-body problem.
Since r couples to φ 1,2 and R 1,2 , equation (4.6) exhibits in principle both spin-orbit coupling and coupling between translation and elastic deformation. This is applicable to studies of long-timescale orbital evolution accounting for tidal dissipation (Ćuk et al. 2016;Lock & Stewart 2017).
As for the elastic motions φ 1 and φ 2 , not only are they coupled to r, but also, in principle, to each other. Still, body 2's elastic deformation should generally have a negligible effect on body 1's. This approximation would be relevant to present-day tidal deformation on Earth. In other words, in order to compute φ 1 one could reasonably approximate the binding energy as Noting that the integration in x ′ trivially yields body 2's mass, and rearranging Γ's argument, it follows that φ 1 's dynamics would be described by the action The third term, through which φ 1 feels body 2, represents a tidal potential. That term will also provide a torque for Ω 1 to respond to. The extension to N -body systems is clear. Smith 1976). It motivates a more general decomposition of the motion whereby each layer has its own Euclidean motion.

Application
We start with eqs. (3.26) describing an N -layer fluid-solid planet. Then we regard the restrictions φ| B i ≡ φ i as distinct motions and decompose each of them in the same way as we decomposed the overall motion (c.f. eqs. 3.10): Recall that the φ of eqs. (3.26) really represents internal motion with respect to the Tisserand frame, so in expression (4.10) Φ i and R i represent relative Euclidean motions, and φ i relative internal motions. The frame defined by each of the Φ i , R i is defined by imposing the constraints once again in precise analogy with eqs. (3.10). Within the original Tisserand frame, we are in effect defining 'relative Tisserand frames', one for each layer. Recalling from eq. (3.28) that we constrained the test-functions analogously to the motion, we make a corresponding decomposition of the test-functions: where each w i is an arbitrary test-function with support on B i while Ψ i and Θ i are arbitrary elements of R 3 and so (3) respectively. With decompositions (4.10) and their constraints (4.11), decomposition (3.10)'s original Euclidean constraints -which of course define the Tisserand Ad R i I i · Ω i = 0, (4.14) where Ω i = R T iṘ i and each layer's relative MoI is defined by (c.f. eq. 2.99b) Similarly for the test functions: The tangential-slip constraints mix together the relative Euclidean and internal motions: Now we substitute all these expressions into eqs. (3.26). We continue to regard the referential potential as a single field, although we will sometimes use ζ i as a shorthand for the restriction Substituting these expressions into eqs. (3.26), and defining the useful quantities showing that it is indeed unaffected by the planet's internal interactions. Meanhwile, the N relative Tisserand frames each obeẏ for arbitrary Ψ i and Θ i . Here we see the boundary tractions contributing net forces and torques that influence the respective layers' CoMs' motions and orientations. As for the relative internal motions, they obey

DISCUSSION
We have derived equations of motion for variably rotating, self-gravitating, (visco)elastic, fluid-solid planets relative to the Tisserand frame. The main results are the exact equations presented in Section 3.4.4, their quasi-static form in Section 3.7, and the resulting linearised equations of motion in Section 3.8. The final form of the equations of motion obtained differs only rather slightly from standard approaches based on a steadily rotating reference frame.
There would, therefore, be no substantial cost or difficult in implementing this new formulation, but by doing so the potential range of applicability is extended. Throughout this paper we have emphasised the role of symmetry, and within Section 3.5 used this approach to identify a seemingly new conservation law related to fluid regions. Section 4 extended the 'Tisserand approach' to problems concerning the interactions of multiple bodies. This paper's new theoretical framework will form the basis for a range of future studies.

ACKNOWLEDGMENTS
MM's work on this paper has been supported by an EPSRC studentship and a CASE award from BP, and more recently by the European Research Council (agreement 833848-UEMHP) under the Horizon 2020 programme. DA has been supported through the Natural Environment Research Council grant numbers NE/V010433/1 and NE/X013804/1.

Data availability statement
There are no data or code relating to this work. where × represents the standard cross-product between 3D vectors, and where Ω is a vector defined through the final equality. While the use of vectors and cross products is perhaps more common that anti-symmetric matrices, the latter approach simplifies various arguments and so is preferred within this work.
Finally, let us ask how torques map from R 3 into so(3). Torques commonly take the form Thus t(3) too is isomorphic to R 3 and has the same action on all motions. We also remark that the adjoint representation of T(3), Ad T · τ ≡ Tτ T −1 , (A.21) acts on t(3) as the identity: for an arbitrary motion, so We have used the fact that τ 's action on any motion is just V, on which T acts as the identity.

A2 Inner-products
This paper uses several different inner-products. We have the standard inner-product between two Cartesian vectors, which may be defined by ⟨a, b⟩ = 1 4 ∥a + b∥ 2 − ∥a − b∥ 2 (A.24) for arbitrary a, b ∈ R 3 . We define the inner product between arbitrary matrices A, B ∈ GL (3) as ⟨A, B⟩ GL(3) = tr AB T , (A.25) where tr (·) denotes the trace operation. This leads us to define the inner-product between elements of so (3) as Ω, Ω ′ so(3) = −tr ΩΩ ′ (A.26) for arbitrary Ω, Ω ′ ∈ so(3). Our earlier definition (A.13) of the wedge-product means that ⟨b, ω · a⟩ ≡ ⟨ω, b ∧ a⟩ for arbitrary a, b ∈ R 3 , ω ∈ so (3), (A.27) which demonstrates a manipulation that we will employ often. While, within this section, we have used subscripts to indicate the space over which aninner product is defined, within the main text we leave this implicit to avoid notational clutter.

A3 Miscellaneous
We use the notation • to denote composition of two functions. For example, If functions depend on both space and time, then composition will always be on the spatial variable unless stated otherwise: while at first order Now, χ (1) ∈ T id Diff (Σ F S ), so by analogy with the reasoning that led to eq. (3.14), the linearised tangential slip constraint on the motion is Now consider the tangential slip constraint as applied to the test-functions. We will need the third-rank tensor Q that is defined as for arbitrary a ∈ R 3 (see A18). It is then a matter of simple algebra to show that where we have used eq. (B.3) to produce the terms in Q.
Finally, let F φ (0) be the space of motions u satisfying (B.6) so too is w at leading order. From now on we therefore write: